Sage Journals: Discover world-class research

Abstract

This paper aims to investigate the effect of viscoelastic behavior of polymer matrix of unidirectional fiber-reinforced laminated composite on stress distribution around the pin-loaded hole under tensile loading. The Laplace transform is used to prevent the integral form of matrix governing stress-strain relation. Applying a micromechanical model, all equilibrium equations for the fibers are written analytically in the Laplace domain. The numerical algorithm of Gaver–Stehfest is implemented, and the governing equations were solved at any given time to extract the concerned results in the time domain. The obtained results are validated against the Finite Element Method results obtained through ANSYS software. Moreover, a comparison of the results of this study at the time equal zero with elastic solutions of other references showed a good agreement. The results revealed that in the long term, the maximum tensile load in the intact fiber around the pinhole was enlarged and the tensile load in fibers far from the pinhole slightly was decreased. Moreover, the location of the maximum axial load that had occurred on pinhole edges was moved slightly toward the center over time.

Keywords

Viscoelasticity PMC pin-loaded joints micromechanical model long-term stress distribution

Introduction

Mechanical joints are introduced as one of the most vulnerable points of a structure, among which is the pin joint that is a frequently used joint in composite materials. Pin joint is utilized if there is no need for preload. If a joint in a composite structure has not been properly designed, it then acts as the damage initiation factor leading to the piece failure at that point. Many studies have been carried out experimentally, analytically, and numerically in this field. Echavarría et al.¹ developed an analytical solution to determine stress distribution around pin-loaded holes in mechanically fastened composites with elastic behavior and showed that the presented solution reduces the computation. Also, a closed-form solution method was extended for stress analysis of pin-loaded joints of laminated composite materials by Grüber et al.² In a review article, the typical failure modes of pin-loaded polymeric matrix composites (PMC) were investigated by Thoppul et al.³

Stress analysis around the pinhole has been the subject of many researches. Aluko and Whitworth⁴ applied the complex stress function method to investigate the contact stress distribution around the pinholes and circular/elliptical pinhole shape dependency to friction. In this respect, Aluko et al.⁵ presented an analytical approach by assuming coulomb friction on the contact areas to obtain the peak of radial, hoop, and shear stress in the pin-loaded joints of laminated composite materials and assumed the pin as a rigid body and neglected from any clearance between pin and hole. Wu et al.⁶ presented a closed-form stress analysis to determine hoop and radial stresses around interference fit pinholes of composite plate joints. Zhou et al.⁷ studied the influence of pin profiles on the net-tension and bearing stresses under tensile loading using FEM and investigated the effect of material orthotropy and load tolerance capacity. Prakash and Hithendra⁸ studied the effect of an additional hole around the pin-loaded hole on the stress-relieving through a 2-D FEA using ANSYS.

Micromechanical modeling has been one of the approaches of pin-loaded joint analysis. Shishesaz et al.⁹ studied stress distribution around pinhole of unidirectional laminated composite using a micromechanical model and shear-lag theory. Also, the influence of fiber arrangement on displacement and stress distribution was investigated. Using micromechanical analysis based on shear-lag theory and experimental tests, failure modes of two serials pin-loaded joint of unidirectional laminated PMC was investigated by Attar et al.¹⁰ Another study about micromechanical modeling of unidirectional laminated composite based on shear-lag theory can be referred to the work of Barati et al.,¹¹ in which the load transfer mechanism in PMC fiber-reinforced with triangular cross-section was studied.

Polymer-based composite materials mainly show creep and relaxation behavior or other forms of viscoelastic behavior, which are more noticeably magnified once a composite material is subjected to the effect of humidity or increased temperature. Viscoelastic behavior can be significantly observed in the matrix of polymer composite materials. Since polymer materials have distinctive viscoelastic properties at different temperatures such as environment temperature, viscoelasticity has an important influence on their behavior and strength and it is needed to be incorporated into the design process. Abadi¹² micromechanically analyzed the stress relaxation behavior of PMC by considering viscoelastic properties of a polymeric matrix to study the time-dependent responses. Andrianov et al.¹³ presented an analytical method to study load transfer in a composite material with elastic steel fiber and viscoelastic poly-methyl-methacrylate matrix using the correspondence principle. Using shear-lag theory, Mondali et al.¹⁴ developed an analytical model in which they predict steady-state creep and stress distribution of short fiber composites. Gusev and Kern¹⁵ developed a micromechanics design of unidirectional carbon/epoxy and glass/epoxy composites using finite element homogenization method in a frequency domain to evaluate the viscoelastic moduli. Hofer et al.¹⁶ analytically extended the so-called Chamis-equations to determine the effective viscoelastic properties of unidirectional fiber-reinforced composite materials using fractional viscoelastic models for polymeric matrix. Kotelnikova-Weiler et al.¹⁷ used micromechanical model based on shear-lag theory to investigate the progressive damage of unidirectional fiber-reinforced composite with viscoelastic polymeric matrix.

The literature review indicates that matrix behavior has been considered elastic in most studies on pin-loaded joints of polymer matrix composite materials. Furthermore, the majority of applied methods in the analysis of pin joints have been experimental in which analytical methods were impractical. Also, it is worth noting that all numerical methods have been performed in the finite element framework using software simulations. In the present study, the effect of the viscoelastic behavior of polymeric matrix of laminated composite materials on the stress distribution around the pinhole has been studied. Using shear-lag theory and a micromechanical approach, the governing equilibrium equations of fibers have been obtained in terms of axial displacement in the Laplace domain. Governing differential equations solved analytically and simultaneously inversed by a numeric algorithm in the time domain.

Micromechanical modeling

Development of laminate displacement relations in Laplace domain

In this investigation, the micromechanical model of the laminated composites was used. Regarding the assumption of viscoelastic behavior for polymer matrix and to avoid the integro-differential form of relations, all of the relations were expressed in the Laplace domain from the beginning. According to the low modulus of polymer matrix against the fiber glass of laminated composite in this research, displacement relations in the entire model were derived based on the shear-lag theory. By applying the boundary conditions corresponding to the pin effect, the desired pin hole is modeled. The following assumptions were made to derive the governing equations.

Fiber and matrix are homogenous.

Unidirectional fibers are considered along the loading direction.

Normal Stress is applied only in the x-direction.

Polymer matrix has a linear viscoelastic behavior and the fibers are linear elastic.

A thorough bond is established between each fiber and the surrounding matrix.

For more simplified modeling, the cross-sections of fibers are approximated as squares.

Considering only axial displacement in the fibers, the tensile force exerted on each fiber (n-th fiber from m-th lamina) is defined in the Laplace domain as equation (1)

{\hat{p}}_{n, m} = E_{f} A_{f} \frac{\partial {\hat{u}}_{n, m}}{\partial x}

(1)

where

E_{f}

and

A_{f}

is the elastic modulus and cross-section area of fibers, respectively. Also,

u_{n, m}

and

p_{n, m}

stand for displacement and load in the n-th fiber from the m-th lamina, respectively. The symbol

"^"

indicates that the variable has been transformed into the Laplace domain. Furthermore, shear stresses in the viscoelastic matrix are expressed as given in equation (2)¹⁸

τ_{y x} = \int_{0}^{t} G_{m} (t - ξ) \frac{\partial γ_{y x}}{\partial ξ} d ξ, τ_{z x} = \int_{0}^{t} G_{m} (t - ξ) \frac{\partial γ_{z x}}{\partial ξ} d ξ

(2)

Here, $G_{m}$ denotes the matrix shear relaxation modulus function. In the present study, the Prony-series model was used to model the viscoelastic behavior of the phenolic-epoxy matrix. The shear relaxation modulus function for this model is defined as equation (3)¹⁸

G_{m} (t) = G_{\infty} + \sum_{i = 1}^{N} G_{i} e^{- t / τ_{i}}

(3)

Given the following shear strain relations

γ_{y x} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}, γ_{z x} = \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x}

(4)

and assuming that the displacement in the y- and z-directions of the matrix is not a function of x, the shear stresses in the viscoelastic matrix will be in the form of equation (5)

τ_{y x} = \int_{0}^{t} G_{m} (t - ξ) \frac{\partial}{\partial ξ} (\frac{\partial u}{\partial y}) d ξ, τ_{z x} = \int_{0}^{t} G_{m} (t - ξ) \frac{\partial}{\partial ξ} (\frac{\partial u}{\partial z}) d ξ

(5)

Applying the Laplace transform to the convolution integral of equations (5) and (8) is derived

{\hat{τ}}_{y x} = s {\hat{G}}_{m} (s) \frac{\partial \hat{u}}{\partial y}, {\hat{τ}}_{z x} = s {\hat{G}}_{m} (s) \frac{\partial \hat{u}}{\partial z}

(6)

Here, ${\hat{G}}_{m} (s)$ is the transformed shear relaxation modulus of the viscoelastic matrix in the Laplace domain, which is as follows for the Prony-series model¹⁸

{\hat{G}}_{m} (s) = \frac{G_{\infty}}{s} + \sum_{i = 1}^{N} \frac{G_{i} τ_{i}}{1 + τ_{i} s}

(7)

Now, writing the shear stresses in the matrix bays in terms of the displacement of adjacent fibers, one can extract the shear stress around the n-th fiber of the m-th lamina according to equation (8)

\begin{array}{l} {({\hat{τ}}_{y x})}_{n + 1, m} = s {\hat{G}}_{m} (s) \frac{({\hat{u}}_{n + 1, m} - {\hat{u}}_{n, m})}{t_{m}}, {({\hat{τ}}_{y x})}_{n - 1, m} = s {\hat{G}}_{m} (s) \frac{({\hat{u}}_{n, m} - {\hat{u}}_{n - 1, m})}{t_{m}} \\ {({\hat{τ}}_{z x})}_{n, m + 1} = s {\hat{G}}_{m} (s) \frac{({\hat{u}}_{n, m} - {\hat{u}}_{n, m + 1})}{t_{m}}, {({\hat{τ}}_{z x})}_{n, m - 1} = s {\hat{G}}_{m} (s) \frac{({\hat{u}}_{n, m - 1} - {\hat{u}}_{n, m})}{t_{m}} \end{array}

(8)

where,

t_{m}

is the width of matrix bays and

t_{f}

is the width of fiber cross-section. The equilibrium equations of fibers neighboring four other fibers and embedded in the laminate (Figure 1), not on the edges or corners of the laminate, can be estimated as equation (9).

Figure 1.

The free-body diagram of the fibers in the middle of laminate and its surrounding matrix bays.

Replacing shear stresses using equation (8) in the equilibrium equations for these fibers yields the following equations in terms of displacements

\frac{\partial {\hat{p}}_{n, m}}{\partial x} d x + {{({\hat{τ}}_{y x})}_{n + 1, m} - {({\hat{τ}}_{y x})}_{n - 1, m} + {({\hat{τ}}_{z x})}_{n, m - 1} - {({\hat{τ}}_{z x})}_{n, m + 1}} t_{f} d x = 0

(9)

\frac{\partial^{2} {\hat{u}}_{n, m}}{\partial x^{2}} + s {\hat{G}}_{m} (s) \frac{t_{f}}{E_{f} A_{f} t_{m}} [{\hat{u}}_{n + 1, m} + {\hat{u}}_{n - 1, m} - 4 {\hat{u}}_{n, m} + {\hat{u}}_{n, m + 1} + {\hat{u}}_{n, m - 1}] = 0; f o r 2 \leq n \leq N - 1, 2 \leq m \leq M - 1

(10)

where N is the number of total fibers in each layer and M is the number of layers. Considering the following vectors and using a similar procedure for the fibers located at the edges and corners of the laminate^19,20

\begin{array}{l} \hat{u} = {[{\hat{u}}_{1,1}, {\hat{u}}_{2,1}, \dots, {\hat{u}}_{N, 1}, {\hat{u}}_{1,2}, \dots, {\hat{u}}_{N - 1, M}, {\hat{u}}_{N, M}]}^{T} \\ \hat{u}'' = {[\hat{u}''_{1,1}, \hat{u}''_{2,1}, \dots, \hat{u}''_{N, 1}, \hat{u}''_{1,2}, \dots, \hat{u}''_{N - 1, M}, \hat{u}''_{N, M}]}^{T} \end{array}

(11)

The governing equations of laminate can be written in the form of equation (12)

\hat{u}'' (x, s) - L (s) \hat{u} (x, s) = 0

(12)

where the matrix L is defined as follows

{\underline{L}}_{(N \times M) \times (N \times M)} (s) = s {\hat{G}}_{m} (s) δ \times [\begin{array}{l} {\underline{L}}_{N \times N}^{2} & {\underline{I}}_{N \times N} & {\underline{0}}_{N \times N} & \dots & {\underline{0}}_{N \times N} \\ {\underline{I}}_{N \times N} & {\underline{L}}_{N \times N}^{3} & {\underline{I}}_{N \times N} & ⋱ & ⋮ \\ {\underline{0}}_{N \times N} & ⋱ & ⋱ & ⋱ & {\underline{0}}_{N \times N} \\ ⋮ & ⋱ & {\underline{I}}_{N \times N} & {\underline{L}}_{N \times N}^{3} & {\underline{I}}_{N \times N} \\ {\underline{0}}_{N \times N} & \dots & {\underline{0}}_{N \times N} & {\underline{I}}_{N \times N} & {\underline{L}}_{N \times N}^{2} \end{array}]

(13)

Here, $δ = \frac{t_{f}}{E_{f} A_{f} t_{m}}$ , ${\underline{\underline{I}}}_{N \times N}$ is the identity matrix, and ${\underline{\underline{L}}}_{N \times N}^{2}$ and ${\underline{\underline{L}}}_{N \times N}^{3}$ are defined below

{\underline{L}}_{N \times N}^{2} = [\begin{array}{l} - 2 & 1 & 0 & \dots & 0 \\ 1 & - 3 & 1 & ⋱ & ⋮ \\ 0 & ⋱ & ⋱ & ⋱ & 0 \\ ⋮ & ⋱ & 1 & - 3 & 1 \\ 0 & \dots & 0 & 1 & - 2 \end{array}]

(14)

{\underline{L}}_{N \times N}^{3} = [\begin{array}{l} - 3 & 1 & 0 & \dots & 0 \\ 1 & - 4 & 1 & ⋱ & ⋮ \\ 0 & ⋱ & ⋱ & ⋱ & 0 \\ ⋮ & ⋱ & 1 & - 4 & 1 \\ 0 & \dots & 0 & 1 & - 3 \end{array}]

(15)

The solution $\underline{\hat{u}} = \underline{\hat{R}} (s) e^{{\hat{λ}}_{i} (s) x}$ is used to solve the system of differential equations of equation (12) in the Laplace domain. Substituting this solution in the governing equations leads to equation (16)

\underline{\hat{R}} (s) {\hat{λ}}_{i}^{2} (s) e^{{\hat{λ}}_{i} (s) x} - \underline{L} (s) \underline{\hat{R}} (s) e^{{\hat{λ}}_{i} (s) x} = \underline{0}

(16)

Assuming, one can obtain eigenvalues and eigenvector $\underline{\hat{R}} (s)$ of the matrix $\underline{\underline{L}} (s)$ , which are a function of the Laplace variable $s$ . Therefore, given the eigenvalues, eigenvectors, and boundary conditions in the Laplace domain, the solution of the system can be obtained.

Unidirectional laminated composite with single pin-loaded holes

The laminate under tensile loading is connected to another part by a pin joint. The geometric dimensions of the laminate can be calculated via the following relations

\begin{array}{l} t_{c} = M t_{f} + (M - 1) t_{m} \\ b = N t_{f} + (N - 1) t_{m} \\ d^{p} = r t_{f} + (r - 1) t_{m} \end{array}

(17)

where

t_{c}

b

, and

d^{p}

stand for laminate thickness, laminate width, and pinhole diameter, respectively. Also, f is the first cut fiber and r is the total number of cut fibers by pinhole. A single pinned laminate is divided into parts A and C. Part A located between the laminate edge and back of the pin is called the short-fiber region, and part C is located between the pin and the under-tension edge of laminate is named the long-fiber region. This way, the fibers displacement equations in these two parts are obtained along with the axial load relation in the Laplace domain.

For the short-fiber region, equations (18) and (19) pertain to the displacement and axial load, respectively

{\hat{u}}_{n, m}^{A} = \sum_{i = 1}^{N \times M} A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x}

(18)

{\hat{p}}_{n, m}^{A} = E_{f} A_{f} \sum_{i = 1}^{N \times M} A_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{{\hat{λ}}_{i} x} - B_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x}

(19)

For the long-fiber region, the displacement and axial load are as equations (20) and (21), respectively

{\hat{u}}_{n, m}^{C} = \sum_{i = 1}^{N \times M} C_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x} + \frac{P_{0} x}{s (N \times M) E_{f} A_{f}}

(20)

{\hat{p}}_{n, m}^{C} = - E_{f} A_{f} \sum_{i = 1}^{N \times M} C_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x} + \frac{P_{0}}{s (N \times M)}

(21)

Boundary conditions

For a single pin-loaded hole case, there are four groups of boundary conditions, as illustrated in Figure 2.

I. On the right-side free surface of region A (α boundary), where x_n,m^α=0, the axial load has a zero value, that is

{{\hat{p}}_{n, m}^{A} |}_{x = x_{n, m}^{α} = 0} = 0 f o r 1 \leq n \leq N, 1 \leq m \leq M

(22)

Figure 2.

Boundaries classification of laminated composite with single pin-loaded holes.

Thus

\sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{{\hat{λ}}_{i} x_{n, m}^{α}} - B_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n, m}^{α}}) = 0, 1 \leq n \leq N, 1 \leq m \leq M

(23)

II. At boundary β and the intersection of intact short and long fibers, the continuity condition of axial load and displacement must be satisfied. Thus, the displacement equality condition for these fibers results in the following equation

\begin{array}{l} {{\hat{u}}_{n, m}^{A} |}_{x = x_{n, m}^{β}} = {{\hat{u}}_{n, m}^{C} |}_{x = x_{n, m}^{β}} f o r {\begin{array}{l} 1 \leq n < f, 1 \leq m \leq M \\ f + r \leq n \leq N, 1 \leq m \leq M \end{array} \\ \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n, m}^{β}} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{β}} - C_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{β}}) = \frac{P_{0} x_{n, m}^{β}}{s (N \times M) E_{f} A_{f}} \end{array}

(24)

Also, equation (25) is established for the axial load continuity condition

{{\hat{p}}_{n, m}^{A} |}_{x = x_{n, m}^{β}} = {{\hat{p}}_{n, m}^{C} |}_{x = x_{n, m}^{β}} f o r {\begin{array}{l} 1 \leq n < f, 1 \leq m \leq M \\ f + r \leq n \leq N, 1 \leq m \leq M \end{array} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{{\hat{λ}}_{i} x_{n, m}^{β}} - B_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n, m}^{β}} + C_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n, m}^{β}}) = \frac{P_{0}}{s (N \times M) E_{f} A_{f}}

(25)

III. On the semi-cylinder surface at the back of the pin, that is, boundary γ, fibers are in contact with the pin. At this boundary, two cases of the rigid and elastic pin can be considered. Considering the rigid pin, radial displacement of fibers in x_n,m^γ can be assumed zero. According to Figure 3(a), the following relation can be written

x_{n, m}^{γ} = x_{p} - d_{n, m}

(26)

Figure 3.

The cross-section view of the joint and elements of the fibers cut by the pinhole. (a) Location of the fibers cut by the pinhole, (b) hole-pin of side left at element fiber, (c) hole-pin of side right at element fiber.

Here, $d_{n, m}$ denotes the edge distance of the n-th fiber of the m-th lamina along the x-axis, and $x_{p}$ is the coordinate of the pinhole center along the x-axis. However, if the pin is not rigid and its stiffness is not large compared to the fibers of laminate, the equilibrium of forces can be written in the radial direction according to Figure 3(b) and (c).

This equilibrium equation in the radial direction is as follows

cos θ_{n, m} \times [{\hat{p}}_{n, m}^{A} - t_{f}^{2} tan θ_{n, m} {({\hat{τ}}_{y x})}_{n + 1, m} + \frac{1}{2} t_{f}^{2} tan θ_{n, m} {({\hat{τ}}_{z x})}_{n, m + 1} - \frac{1}{2} t_{f}^{2} tan θ_{n, m} {({\hat{τ}}_{z x})}_{n, m - 1}] + {\hat{F}}_{n, m}^{r} = 0

(27)

It should be noted that the obtained relations at this boundary are related to the fibers on the left side of the pinhole. For the right-side fibers, it is sufficient to replace $n + 1$ with $n - 1$ in the former relations. The pin stiffness is assumed as linear springs radially positioned.⁹ As shown in Figure 3(a), the following relation exists for $tan θ_{n, m}$

tan θ_{n, m} = \frac{Y_{n, m}}{d_{n, m}}

(28)

where

d_{n, m} = \sqrt{{(\frac{d^{p}}{2})}^{2} - Y_{n, m}^{2}}, Y_{n, m} = | (n - \frac{1}{2}) t_{f} + (n - 1) t_{m} - \frac{b}{2} |

(29)

{\hat{F}}_{n, m}^{r}

is also obtained by equation (30)²¹

{\hat{F}}_{n, m}^{r} = K {\hat{u}}_{n, m} cos θ_{n, m}

(30)

Here, $K$ is the stiffness between the pin and hole edge. Now, replacing equation (30) into equation (27), the boundary condition of $γ$ can be written as follows

{\hat{p}}_{n, m}^{A} - t_{f}^{2} tan θ_{n, m} {({\hat{τ}}_{y x})}_{n + 1, m} + \frac{1}{2} t_{f}^{2} tan θ_{n, m} {({\hat{τ}}_{z x})}_{n, m + 1} - \frac{1}{2} t_{f}^{2} tan θ_{n, m} {({\hat{τ}}_{z x})}_{n, m - 1} + K {\hat{u}}_{n, m} = 0, [f + \frac{r}{2}] \leq n \leq f + r - 1, 2 \leq m \leq M - 1

(31)

Replacing displacement-based shear stress expression in equation (31) yields the following equation

{\hat{p}}_{n, m}^{A} - (tan θ_{n, m} \frac{s {\hat{G}}_{m} (s) t_{f}^{2}}{t_{m}}) [{\hat{u}}_{n + 1, m} - 2 {\hat{u}}_{n, m} + \frac{1}{2} {\hat{u}}_{n, m + 1} + \frac{1}{2} {\hat{u}}_{n, m - 1}] + K {\hat{u}}_{n, m} = 0

(32)

Finally, by replacing axial load and displacement of equations (18)–(21) in the former relation, the boundary condition of the elastic pin for $2 \leq m \leq M - 1$ is as follows

\begin{array}{l} E_{f} A_{f} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} - B_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}}) \\ \begin{array}{l} - (tan θ_{n, m} \frac{s {\hat{G}}_{m} (s) t_{f}^{2}}{t_{m}}) {\sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + 1 + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n + 1, m}^{γ}} + B_{i} {\hat{R}}_{n + 1 + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n + 1, m}^{γ}}) \\ - 2 \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}}) + \frac{1}{2} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N m}^{i} e^{{\hat{λ}}_{i} x_{n, m + 1}^{γ}} + B_{i} {\hat{R}}_{n + N m}^{i} e^{- {\hat{λ}}_{i} x_{n, m + 1}^{γ}}) \\ + \frac{1}{2} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 2)}^{i} e^{{\hat{λ}}_{i} x_{n, m - 1}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 2)}^{i} e^{- {\hat{λ}}_{i} x_{n, m - 1}^{γ}}) \\ + K \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}})} = 0, 2 \leq m \leq M - 1 \end{array} \end{array}

(33)

Similarly, this condition for $m = 1$ is as the following

{\hat{p}}_{n, m}^{A} - (tan θ_{n, m} \frac{s {\hat{G}}_{m} (s) t_{f}^{2}}{t_{m}}) [{\hat{u}}_{n + 1, m} - \frac{3}{2} {\hat{u}}_{n, m} + \frac{1}{2} {\hat{u}}_{n, m + 1}] + K {\hat{u}}_{n, m} = 0

(34)

Replacing axial load and displacement relations in equation (34), the boundary condition of the elastic pin for $m = 1$ is as follows

\begin{array}{l} \begin{array}{l} \begin{array}{l} E_{f} A_{f} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} - B_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}}) \end{array} \\ - (tan θ_{n, m} \frac{s {\hat{G}}_{m} (s) t_{f}^{2}}{t_{m}}) {\sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + 1 + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n + 1, m}^{γ}} + B_{i} {\hat{R}}_{n + 1 + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n + 1, m}^{γ}}) \\ - \frac{3}{2} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}}) \\ + \frac{1}{2} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N m}^{i} e^{{\hat{λ}}_{i} x_{n, m + 1}^{γ}} + B_{i} {\hat{R}}_{n + N m}^{i} e^{- {\hat{λ}}_{i} x_{n, m + 1}^{γ}}) \\ + K \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}})} \end{array} \\ = 0, m = 1 \end{array}

(35)

Similarly, for $m = M$

{\hat{p}}_{n, m}^{A} - (tan θ_{n, m} \frac{s {\hat{G}}_{m} (s) t_{f}^{2}}{t_{m}}) [{\hat{u}}_{n + 1, m} - \frac{3}{2} {\hat{u}}_{n, m} + \frac{1}{2} {\hat{u}}_{n, m - 1}] + K {\hat{u}}_{n, m} = 0

(36)

Finally, the boundary condition for $m = M$ is expressed as below

\begin{array}{l} E_{f} A_{f} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} - B_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}}) \\ \begin{array}{l} - (tan θ_{n, m} \frac{s {\hat{G}}_{m} (s) t_{f}^{2}}{t_{m}}) \\ {\sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + 1 + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n + 1, m}^{γ}} + B_{i} {\hat{R}}_{n + 1 + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n + 1, m}^{γ}}) \\ - \frac{3}{2} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}}) \\ + \frac{1}{2} \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 2)}^{i} e^{{\hat{λ}}_{i} x_{n, m - 1}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 2)}^{i} e^{- {\hat{λ}}_{i} x_{n, m - 1}^{γ}}) \\ + K \sum_{i = 1}^{N \times M} (A_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{{\hat{λ}}_{i} x_{n, m}^{γ}} + B_{i} {\hat{R}}_{n + N (m - 1)}^{i} e^{- {\hat{λ}}_{i} x_{n, m}^{γ}})} \end{array} \\ = 0, m = M \end{array}

(37)

IV. The fourth boundary condition pertains to the long fibers cut by the pinhole, that is, boundary μ. These fibers are in no contact with the pin’s outer surface and are considered as stress-free surface. Therefore

{{\hat{p}}_{n, m}^{C} |}_{x = x_{n, m}^{μ}} = 0 f o r f \leq n \leq f + r - 1, 1 \leq m \leq M

(38)

such that

x_{n, m}^{μ} = x_{p} + d_{n, m}

(39)

equation (38) can be expressed as follows

\sum_{i = 1}^{N \times M} C_{i} {\hat{R}}_{n + N (m - 1)}^{i} {\hat{λ}}_{i} e^{- {\hat{λ}}_{i} x_{n, m}^{μ}} = \frac{P_{0}}{s (N \times M) E_{f} A_{f}}

(40)

Inverse Laplace transform

Due to the viscoelastic behavior of the polymeric matrix of laminated composite, the displacement and stresses in the fibers and matrix would be a function of time and location. As demonstrated in the previous sections, axial displacement in the fibers can be obtained in terms of location and Laplace variable $s$ . Regarding the complexity of the relations derived for displacements and stresses in the Laplace domain, it is impossible to perform the inverse Laplace transform procedure by conventional analytical methods. This is one of the typical difficulties in solving viscoelastic problems for complex models. Thus, to evaluate displacements and stresses at a specific time, special inverse Laplace transform algorithms must be implemented.^13,22,23 One of the most useful and beneficial methods for the numerical inverse Laplace transform is to use equation (41)^24–26

f (t) \approx f_{n} (t) = \frac{1}{t} \sum_{k = 0}^{n} Re [ω_{k} \hat{f} (\frac{α_{k}}{t})], t > 0

(41)

If $\hat{f} (s)$ is the Laplace transform of $f (t)$ for positive real values ( $t > 0$ ), $f (t)$ can be obtained through multiplying the inverse of the time by a finite linear combination of the transformed function and replacing Laplace variable with $\frac{α_{k}}{t}$ . Accordingly, to calculate $f (t)$ , it is sufficient to obtain the value of the concerned function at the value of $\frac{α_{k}}{t}$ for $k = 0 \dots n$ . In this study, $\hat{f} (s)$ is the fibers axial displacement function $\hat{u} (x, s)$ , which is not directly and explicitly accessible. Hence, the governing differential equation needs to be solved through the respective boundary conditions as well as the simultaneous implementation of the inverse Laplace transform algorithm for a limited number of terms. Coefficients $α_{k}$ and $ω_{k}$ in the numerical inverse Laplace transform method are complex numbers depending on the value of n. The degree of accuracy to estimate $f_{n} (t)$ depends on the number of employed terms.

Different algorithms can be put into equation (41) within the mentioned general framework, namely, Gaver-Stehfest^13,24, Fixed Talbot,²⁶ and Durbin^27,28. In the present study, Gaver–Stehfest methods have been applied for the inverse transform of the function. Replacing $n = 2 M$ in equation (41), where M stands for the number of Gaver functions, and considering $α_{k} = ln (2)$ and weighting factors $ω_{k} = ξ_{k} ln (2)$ , one can estimate the concerned function at any given time using a limited number of summation terms. Thus, to calculate the inverse of a transformed function in the Laplace domain by the Gaver–Stehfest method, the following relation is used

f (t) \approx f_{M} (t) = \frac{ln (2)}{t} \sum_{k = 1}^{2 M} ξ_{k} \hat{f} (\frac{k ln (2)}{t}), t > 0

(42)

where the corresponding coefficients are as follows

ξ_{k} = {(- 1)}^{k + M} \sum_{j = [(k + 1) / 2]}^{m i n {k, M}} \frac{j^{M + 1}}{M!} (\begin{array}{l} M \\ j \end{array}) (\begin{array}{l} 2 j \\ j \end{array}) (\begin{array}{l} j \\ k - j \end{array}), (\begin{array}{l} n \\ m \end{array}) = \frac{n!}{m! (n - m)!}

(43)

Finite element method

Finite element simulation via ANSYS software has been used to verify the obtained results from the semi-analytical method. A finite element model was prepared and solved for the stress distribution around the pin-loaded hole of laminated composite using ANSYS v11. The element SOLID186 has been employed to mesh this model. This is a higher-order 3D 20-node solid element that exhibits quadratic displacement behavior and is well suited to modeling irregular meshes. Each node of this element has 3 degrees of freedom. In the FEM simulation by ANSYS, the displacement of the semi-cylinder surface in the back of the pinhole is constrained in the radial direction. Furthermore, the displacement of all nodes in the model is constrained in the z-direction. In the ANSYS software, the generalized Maxwell model and Prony series are used to model the relaxation modulus of the viscoelastic material. The Prony series must be rewritten as equation (44) to incorporate the required constants into the ANSYS software

G (t) = G_{0} [ξ_{\infty} + \sum_{i = 1}^{N} ξ_{i} e^{- t / τ_{i}}], ξ_{i} = \frac{G_{i}}{G_{\infty}}

(44)

where

G_{0}

is inserted as the initial shear elastic modulus and

ξ_{i}

is the relative shear modulus.

Results and discussions

The proposed semi-analytical approach in this study has been programmed and calculated by MATLAB software. The effect of some parameters, including volume fraction of laminated composite and some geometric characteristics of laminate, is evaluated in this study. In all of the obtained results, the applied tensile load is equal to 100 N. The laminated composite is also included of five laminae of glass fiber and phenolic-epoxy matrix. Table 1 lists the constants of the Prony series employed to model the viscoelastic behavior of the phenolic-epoxy polymer matrix.

Table 1.

Constants of Prony series for phenolic epoxy (G₀ = 9.8 GPa).

i	$G_{i}$ (MPa)	$η_{i}$ (Pa s)	I	$G_{i}$ (MPa)	$η_{i}$ (Pa s)
1	123	1.1 × 10⁻⁶	6	589	0.2
2	386	1.1 × 10⁻⁵	7	593	2
3	120	9.8 × 10⁻⁵	8	262	21.6
4	260	1.2 × 10⁻³	9	11.5	4.4 × 10³
5	666	1.0 × 10⁻²	10	66.9	193.4

First, the results of the proposed semi-analytical method in this study are validated against the results obtained by the FEM using the ANSYS software. Tables 2 and 3 compare the results obtained by ANSYS and the results of the proposed semi-analytical method, respectively, for the maximum tensile stress (on side edges of the pinhole) and maximum compressive stress (in the back of the pinhole). The presented results imply a significantly good agreement between both methods, particularly for tensile stress on the pinhole edge.

Table 2.

Comparison of the maximum tensile stress around the pinhole between the FEM and the presented semi-analytical method.

Time (sec.)		0	10	100	1000
Tensile stress (MPa)	Present study	1284	1176	1136	1127
Tensile stress (MPa)	ANSYS	1260	1140	1111	1099
Error (%)		2	3	2	3

Table 3.

Comparison of the maximum compressive stress around the pinhole between the FEM and the presented semi-analytical method.

Time (sec.)		0	10	100	1000
Compressive stress (MPa)	Present study	481	608	650	658
Compressive stress (MPa)	ANSYS	547	690	731	733
Error (%)		12	11	11	10

Figure 4 shows the axial load distribution in the fibers of a laminated composite with a pin-loaded hole at t = 0, 100, and 1000 s. According to this figure, as time passes, the intensity of the axial load increases around the pinhole.

Figure 4.

Distribution of axial load in the fibers for time zero, 100, and 1000 s. (a) t = 0, (b) t = 100 s, (c) t = 1000 s.

Figure 5 depicts the axial load of the fiber passing through the middle of the pinhole at t = 0, 10, 100, and 1000s, followed by a comparison with the elastic solution of reference 10⁹ at t = 0. It can be observed in the compressive region in the back of the pinhole that the compressive load increases at first until about 10 s and then decreases over time, while the axial load of the under-tension region undergoes a descending trend. Besides, the elastic solution by the corresponding elastic model of reference⁹ agrees well with the viscoelastic solution at t = 0. Similarly, this evaluation is performed for the axial load of the first intact fiber neighboring the pinhole in Figure 6, showing that the tensile load available in the intact fiber around the pinhole is elevated over time, unlike the fiber passing through the middle of the pinhole.

Figure 5.

The axial load of the fiber passing through the middle of the pinhole at different times and compared with the elastic solution of reference 9.

Figure 6.

The axial load of the first intact fiber adjacent to the pinhole for different times and compared with the elastic solution of the reference 9.

Figure 7 displays the matrix shear stress distribution in the shear-out plane for different times along with a comparison with the elastic response of the presented model. As can be observed, the shear stress concentration decays gradually over time. Furthermore, the elastic response is fully consistent with the viscoelastic response at t = 0.

Figure 7.

Shear stress distribution of the matrix in the shear-out plane for different times.

Figure 8 depicts the compressive load distribution in the back of the pin in terms of the angle $θ$ . The maximum compressive load is enlarged as time passes and its location slightly tends from the edge toward the inside of the hole. Figure 9 indicates the axial load in the fibers neighboring the pinhole at t = 0. According to this figure, the effect of axial load concentration is only observed in three fibers adjacent to the pinhole edge. Figure 10 presents a similar investigation on shear stress concentration. Matrix shear stress concentration is observed at a further distance from the pinhole edge compared to the axial load concentration of fibers.

Figure 8.

Distribution of compressive load in the back half of the pinhole in terms of angle θ for different times.

Figure 9.

Axial load distribution in the fibers adjacent to the pinhole at t = 0.

Figure 10.

Shear stress distribution in the fibers adjacent to the pinhole at t = 0.

Conclusion

The purpose of this study was to investigate the effect of the viscoelastic behavior of unidirectional fiber-reinforced polymer matrix composites in the pin-loaded joint by a micromechanical approach. The results of this study indicate relatively large variations in the shear stress distribution in the matrix and axial load on the fibers over time, so that, in the long term, the tensile load in the fibers adjacent to the pinhole increases and the shear stress concentration in the matrix around the pinhole decreases.

Footnotes

Acknowledgments

The authors wish to thank the Ahvaz Branch, Islamic Azad University for the financial support.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Arash Reza

References

Echavarría

Haller

Salenikovich

. Analytical study of a pin-loaded hole in elastic orthotropic plates. Compos Struct 2007; 79: 107–112. 10.1016/j.compstruct.2005.11.038.

Grüber

Hufenbach

Kroll

, et al Stress concentration analysis of fibre-reinforced multilayered composites with pin-loaded holes. Compos Sci Technol 2007; 67: 1439–1450. 10.1016/j.compscitech.2006.08.018 .

Thoppul

Finegan

Gibson

. Mechanics of mechanically fastened joints in polymer-matrix composite structures—A review. Compos Sci Technol 2009; 69: 301–329. 10.1016/j.compscitech.2008.09.037 .

Aluko

Whitworth

. Contact stress analysis around elliptic pin-loaded holes in orthotropic plates. Volume 9: Mechanics of Solids, Structures and Fluids. ASME 2010 International Mechanical Engineering Congress and Exposition, Vancouver, British Columbia, Canada, November 12–18, 2010, pp. 277–283. 10.1115/IMECE2010-37662.

Aluko

Whitworth

Owolabi

. Effect of friction on contact stress distribution in pin-loaded orthotropic plates. Volume 8: Mechanics of Solids, Structures and Fluids. ASME 2012 International Mechanical Engineering Congress and Exposition, Houston, Texas, USA, November 9–15, 2012, pp. 79–83. 10.1115/IMECE2012-86913.

Zhang

Cheng

, et al Analytical modelling for stress distribution around interference fit holes on pinned composite plates under tensile load. Compos B: Eng 2016; 100: 176–185. 10.1016/j.compositesb.2016.06.011 .

Zhou

Fei

Tao

. Profile design of loaded pins in composite single lap joints: from circular to non-circular. Results Phys 2016; 6: 471–480. 10.1016/j.rinp.2016.07.010 .

Prakash

Hithendra

. Combined influence of pin-hole interference and additional hole on the stress intensity factor for a pin loaded finite plate. Proced Struct Integr 2020; 28: 1125–1133. 10.1016/j.prostr.2020.11.127 .

Shishesaz

Attar

Robati

. The effect of fiber arrangement on stress concentration around a pin in a laminated composite joint. ASME 2010 10th Biennial Conference of Engineering Systems Design Analysis, Volume 2. ASME 2010 10th Biennial Conference on Engineering Systems Design and Analysis, Istanbul, Turkey, July 12–14, 2010, 2010, pp. 153–161. 10.1115/ESDA2010-24587.

, et al. Failure analysis of unidirectional polymeric matrix composites with two serial pin loaded-holes. J Mech Sci Technol 2016; 30: 2583–2591. 10.1007/s12206-016-0519-5.

, et al. Micromechanical analysis of unidirectional polymeric composites material with triangular fibers. J Mech Sci Technol 2017; 31: 4755–4762. 10.1007/s12206-017-0922-6.

12.

Abadi

. Micromechanical analysis of stress relaxation response of fiber-reinforced polymers. Compos Sci Technol 2009; 69: 1286–1292. 10.1016/j.compscitech.2009.02.036 .

. Load transfer in fibre-reinforced composites with viscoelastic matrix: an analytical study. Arch Appl Mech 2009; 79: 999–1007. 10.1007/s00419-008-0265-y.

. A new analytical shear-lag based model for prediction of the steady state creep deformations of some short fiber composites. Mater Des 2009; 30: 1075–1084. 10.1016/j.matdes.2008.06.039 .

15.

Gusev

Kern

. Frequency domain finite element estimates of viscoelastic stiffness of unidirectional composites. Compos Structures 2018; 194: 445–453. 10.1016/j.compstruct.2018.02.027 .

, et al. Closed-form expressions for effective viscoelastic properties of fiber-reinforced composites considering fractional matrix behavior. Mech Mater 2018; 127: 14–25. 10.1016/j.mechmat.2018.08.005 .

, et al. Progressive damage of a unidirectional composite with a viscoelastic matrix, observations and modelling. Compos Structures 2018; 188: 297–312. 10.1016/j.compstruct.2017.12.067 .

18.

Brinson

. Polymer Engineering Science and Viscoelasticity: An Introduction. 2nd ed.. NY: Springer, 2015, p. 482.

19.

Reza

Shishesaz

. Transient load concentration factor due to a sudden break of fibers in the viscoelastic PMC under tensile loading. Int J Sol Struct 2016; 88-89: 1–10. 10.1016/j.ijsolstr.2016.04.005 .

20.

Reza

Shishesaz

. The effect of viscoelasticity on the stress distribution of adhesively single-lap joint with an internal break in the composite adherends. Mech Time-Depend Mater 2018; 22: 373–399. 10.1007/s11043-017-9362-z.

. Application of ANSYS to contact between machine elements. Engineering Analysis with ANSYS Software. Oxford: Butterworth-Heinemann, 2006, pp. 331–451.

. A new class of Laplace inverses and their applications. Appl Math Lett 2000; 13: 97–104. 10.1016/S0893-9659(00)00062-8 .

23.

Temel

Şahan

. Transient analysis of orthotropic, viscoelastic thick plates in the Laplace domain. Eur J Mech - A/Solids 2013; 37: 96–105. 10.1016/j.euromechsol.2012.05.008 .

24.

Montella

. LSV modelling of electrochemical systems through numerical inversion of Laplace transforms. I-The GS-LSV algorithm. J Electroanal Chem 2008; 614: 121–130. 10.1016/j.jelechem.2007.11.010 .

25.

Montella

. Re-examination of the potential step chronoamperometry method through numerical inversion of Laplace transforms. II. Application examples. J Electroanal Chem 2009; 633: 45–56. 10.1016/j.jelechem.2009.04.025 .

. New approach of electrochemical systems dynamics in the time-domain under small-signal conditions. I. A family of algorithms based on numerical inversion of Laplace transforms. J Electroanal Chem 2008; 623: 29–40. 10.1016/j.jelechem.2008.06.015 .

. Dual reciprocity boundary element method in Laplace domain applied to anisotropic dynamic crack problems. Comput Struct 2003; 81: 1703–1713. 10.1016/S0045-7949(03)00184-6 .

28.

Wen

Adetoro

. The fundamental solution of Mindlin plates with damping in the Laplace domain and its applications. Eng Anal Boundary Elem 2008; 32: 870–882. 10.1016/j.enganabound.2007.12.005 .

Micromechanical approach to viscoelastic stress analysis of a pin-loaded hole in unidirectional laminated PMC

Abstract

Keywords

Introduction

Micromechanical modeling

Development of laminate displacement relations in Laplace domain

Unidirectional laminated composite with single pin-loaded holes

Boundary conditions

Inverse Laplace transform

Finite element method

Results and discussions

Conclusion

Footnotes

Acknowledgments

Declaration of conflicting interests

Funding

ORCID iD

References