Evaluation of Interlaminar Stresses in Composite Laminates with a Bolt-Filled Hole Using a Linear Elastic Traction-Separation Description

Determination of the local interlaminar stress distribution in a laminate with a bolt-filled hole is helpful for optimal bolted joint design, due to the three-dimensional (3D) nature of the stress field near the bolt hole. A new interlaminar stress distribution phenomenon induced by the bolt-head and clamp-up load, which occurs in a filled-hole composite laminate, is investigated. In order to efficiently evaluate interlaminar stresses under the complex boundary condition, a calculation strategy that using zero-thickness cohesive interface element is presented and validated. The interface element is based on a linear elastic traction-separation description. It is found that the interlaminar stress concentrations occur at the hole edge, as well as the interior of the laminate near the periphery of the bolt head. In addition, the interlaminar stresses near the periphery of the bolt head increased with an increase in the clamp-up load, and the interlaminar normal and shear stresses are not at the same circular position. Therefore, the clamp-up load cannot improve the interlaminar stress distribution in the laminate near the periphery of the bolt head, although it can reduce the magnitude of the interlaminar shear stress at the hole edge. Thus, the interlaminar stress distribution phenomena may lead to delamination initiation in the laminate near the periphery of the bolt head, and should be considered in composite bolted joint design.


Introduction
Interlaminar stresses can lead to damage in the form of delamination and matrix cracking of laminate, especially at the free-edges of a laminate [1,2].Accurate and proper design of bolted composite joints requires the determination of the local stress field within the laminate [3].Furthermore, the stress field in the vicinity of a bolt hole is generally considered to be three-dimensional (3D) [4][5][6].Thus, it is important to take interlaminar stresses into account in composite bolted joint design.
Due to various kinds of bolted composite joints, a filled-hole laminate was used as a means to characterize the stress distribution of bolted joints in this paper.Filled-hole strength has been widely used to generate design data to establish allowables for composite bolted joint analysis [7].Numerous investigators have studied the stress fields in laminates with a bolt.Camanho and Matthews [8] presented a review of the previous investigation into the stress and strength analysis of mechanically fastened joints in fiber-reinforced plastics.They concluded that significant issues such as the effects of friction, clearance or interference and the contact area on the stress distribution around a pin-loaded hole.They also suggested that a 3D model was required to accurately predict the strength of mechanically fastened joints.Following that study, Ireman [9] developed a 3D finite element model (FEM) of composite bolted joints to determine non-uniform stress distribution through the thickness of laminate in the vicinity of a bolt hole.McCarthy et al. [10] analyzed the effects of bolt hole clearance on the mechanical behavior of composite bolted joints using a 3D FEM, and studied radial stress through thickness variations at the hole.Caccese et al. [11] investigated the influence of stress relaxation on clamp-up load in hybrid composite-to-metal bolted connections by experiment, and the pressure distribution at the free face of the laminate under different clamp-up load was obtained.Egan et al. [12] provided an analysis of the stress distribution at a countersunk hole boundary using the nonlinear finite element code Abaqus.Feo et al. [13] examined the stress distribution among the different bolts by varying the number of rows of bolts as well as the number of bolts per row.The results of their study indicated that in the presence of washers, the stress distributions in the fiber direction and varying fiber inclinations decreased for each value of maximum bearing stresses.Zhao et al. [14] proposed a determination method of stress concentration relief factors for failure prediction of composite joints.Moreover, the in-plane stress distribution around the hole edge for open-, filled-and loaded-hole laminates was discussed in their study.
These studies were mainly aimed at the investigation of stress distributions at bolt hole edge or in-plane stress distributions.However, the experiments outlined in References [15,16] have shown that damage might initiate near the edge of the bolt head or washer in composite laminates with bolts.Whitney et al. [3] investigated the singular stress fields near the contact surfaces of the bolt head and laminate in a composite bolted joint, and gave all six stress distributions in the position; however, they did not focus on interlaminar stresses.It is also noteworthy that interlaminar stresses are a stress component along the thickness direction of the laminate, therefore, the interlaminar stress distribution along the thickness direction of laminate with a bolt-filled hole should also be investigated.An understanding of geometry-based effects and interlaminar stress distributions are helpful in explaining the failure mode of laminates with a bolt.
Several levels of modeling can be used for the evaluation of interlaminar stresses.Pagano [17] provided an exact solution for a simply supported plate subject to sinusoidal pressure load, which had been considered a benchmark problem.The approaches based on the equivalent single layer (ESL) theories, layer-wise theories and zigzag theories are introduced in these References [1,18,19].The sets of ESL refer to the Classical Laminated Plate Theory (CLPT), the First Order Shear Deformation Plate Theory (FSDT) and the Higher Order Deformation Plate Theory (HSDT).The limitations of these theories, especially the ESL model, are introduced in References [20][21][22].In addition, some other reviews and extensive assessments were investigated in References [23,24].In a filled-hole laminate, the stress state near the hole is a complex 3D problem.Generally, numerical methods based on 3D elasticity theory are usually used to accurately evaluate interlaminar stresses under complex boundary.The 3D elasticity theory includes traditional 3D elasticity formulations and layerwise theories [1,[25][26][27].In this method, the laminate is modeled as 3D elements [20,28].However, the element type and mesh density are significant factors that can affect the calculation accuracy of FEM.If the laminate is discretized with the 3D solid element, a further challenge is that the transverse shear stresses are still discontinuous at the layer interface without enough elements in thickness direction.To overcome these limitations, some multiple model methods, such as the interface model, are used.Dakshina Moorthy and Reddy [29] developed an interface model based on the penalty function method and obtained the interlaminar stresses using this approach.At present, the cohesive zone model (CZM) is widely used in the interface modeling [30][31][32].Some of the CZMs are described by the so-called traction-separation law (TSL).The TSL is derived from trial-and-error finite element computations [33] and originally designed to be used to model complex fracture mechanisms at the crack tip.
Against this background, there are two main contributions in this paper.Firstly, a new interlaminar stress distribution phenomenon induced by the bolt-head and clamp-up load was investigated, in accordance with experimental failure evidence of composite laminate with a bolt-filled hole [16].This analysis can yield a better understanding of stress distributions in bolt joints and improve reliability and part life.Secondly, the interface model based on TSL was employed to evaluate interlaminar stresses.This modeling approach has been widely applied in composites failure analyses; however, less attention has been dedicated to its role in stress calculation.With this approach, the interlaminar stress distribution can be evaluated efficiently.

Linear Elastic Traction-Separation Constitutive Behavior of the Interface Model
In this approach, a zero-thickness cohesive interface element was introduced into the finite element computations for composite materials.The zero-thickness means that the geometric thickness of cohesive elements is equal to zero, and not the constitutive thickness.Since the initial geometric thickness of the interface element is zero, the deformation state of the zero-thickness element cannot be described by the classical definition of strain.Instead, the constitutive response of the zero-thickness interface element is established in terms of a linear elastic traction-separation description.In general, the available traction-separation model consists of three parts: linear elastic traction-separation; damage initiation criteria; and damage evolution laws [34,35].However, only initially linear elastic behavior is considered due to stress distribution being the only concern in this paper.
The constitutive response of this interface was obtained using a procedure introduced by Camanho and Davila [36].Consider a typical interface that uses a traction-separation description between two plies as shown in Figure 1.The two coincident nodes of the interface element are a pair of nodes, and each interface element has four pairs of nodes, and each node has three degrees of freedom (DOFs); thus, there are 24 DOFs in total for each interface element.The nodes global displacement vector u N is defined as follows: where iz are the displacements of the i nodes of the interface element in the x, y, and z directions.employed to evaluate interlaminar stresses.This modeling approach has been widely applied in composites failure analyses; however, less attention has been dedicated to its role in stress calculation.With this approach, the interlaminar stress distribution can be evaluated efficiently.

Linear Elastic Traction-Separation Constitutive Behavior of the Interface Model
In this approach, a zero-thickness cohesive interface element was introduced into the finite element computations for composite materials.The zero-thickness means that the geometric thickness of cohesive elements is equal to zero, and not the constitutive thickness.Since the initial geometric thickness of the interface element is zero, the deformation state of the zero-thickness element cannot be described by the classical definition of strain.Instead, the constitutive response of the zero-thickness interface element is established in terms of a linear elastic traction-separation description.In general, the available traction-separation model consists of three parts: linear elastic traction-separation; damage initiation criteria; and damage evolution laws [34,35].However, only initially linear elastic behavior is considered due to stress distribution being the only concern in this paper.
The constitutive response of this interface was obtained using a procedure introduced by Camanho and Davila [36].Consider a typical interface that uses a traction-separation description between two plies as shown in Figure 1.The two coincident nodes of the interface element are a pair of nodes, and each interface element has four pairs of nodes, and each node has three degrees of freedom (DOFs); thus, there are 24 DOFs in total for each interface element.The nodes global displacement vector uN is defined as follows: where , , , , , u u u u u u are the displacements of the i nodes of the interface element in the x, y, and z directions.

Element of lower layer
Interface element z,3 x,1 y,2 The relative displacement u for each pair of nodes can be described as where I  is the identity diagonal matrix.
Converting u to the relative displacement at the interface local coordinate, defining the relative displacements tensor in local coordinates.Subscript n, s and t represent the local 3-direction, the 1-and 2-direction, respectively, in which δ k can be written as where B is the geometric matrix.Unlike continuum elements, the interface element stiffness matrix before softening onset requires the penalty stiffness K of the interface material [36].The penalty stiffness K relates to the element traction stresses to the element relative displacements tensor δ k : The relative displacement ∆u for each pair of nodes can be described as where I 12×12 is the identity diagonal matrix.
Converting ∆u to the relative displacement at the interface local coordinate, defining the relative displacements tensor δ k = {δ n , δ s , δ t } T in local coordinates.Subscript n, s and t represent the local 3-direction, the 1-and 2-direction, respectively, in which δ k can be written as where B is the geometric matrix.Unlike continuum elements, the interface element stiffness matrix before softening onset requires the penalty stiffness K of the interface material [36].The penalty stiffness K relates to the element traction stresses σ k = {σ n , σ s , σ t } T to the element relative displacements tensor δ k : From Equation (4), it can be seen that traction-separation constitutive behavior is not stress-strain constitutive behavior.The interface element can be visualized as a spring between the initially coincident nodes of the element.After softening onset, the initially coincident nodes will open or slide relative to each other [37].Due to the zero-thickness interface element describing interlaminar behavior, the three traction stresses can be considered as corresponding interlaminar stresses.In addition, the penalty stiffness K is the interface element constitutive parameter, which has been introduced by some researchers [38].The value of K is usually expected to be higher than the elastic modulus of the material, but it should also be noted that a larger value for K in FEM might lead to a bad convergence [31].

Verification of the Interface Model with a Benchmark Problem
The interlaminar stress calculation strategy based on linear elastic traction-separation behavior was validated by the consideration of a benchmark problem.This benchmark problem was analyzed by Pagano [17], where the 3D exact elasticity solution of a square, simply supported, symmetric cross-ply [0/90/0] laminated plate under a bisinusoidally distributed transverse load q, was provided.As shown in Figure 2, the span and thickness of the laminate is denoted by "a" and "h", and q can be described as q = q 0 sin(πx/a) sin(πy/a) (5) where q 0 is a constant.A laminate of span-to-thickness ratio S = a/h = 4 are considered, and can be considered as a thick laminate.The thickness of the ply is 0.25 mm when a = 3 mm.The following lamina material properties are used: E 1 = 25 × 10 6 psi, E 2 = E 3 = 10 6 psi, G 12 = G 23 = 0.2 × 10 6 psi, ν 12 = ν 13 = ν 23 = 0.25.The resulting stresses are normalized according to the following formulae: Appl.Sci.2017, 7, 93 4 of 16 From Equation (4), it can be seen that traction-separation constitutive behavior is not stress-strain constitutive behavior.The interface element can be visualized as a spring between the initially coincident nodes of the element.After softening onset, the initially coincident nodes will open or slide relative to each other [37].Due to the zero-thickness interface element describing interlaminar behavior, the three traction stresses can be considered as corresponding interlaminar stresses.In addition, the penalty stiffness K is the interface element constitutive parameter, which has been introduced by some researchers [38].The value of K is usually expected to be higher than the elastic modulus of the material, but it should also be noted that a larger value for K in FEM might lead to a bad convergence [31].

Verification of the Interface Model with a Benchmark Problem
The interlaminar stress calculation strategy based on linear elastic traction-separation behavior was validated by the consideration of a benchmark problem.This benchmark problem was analyzed by Pagano [17], where the 3D exact elasticity solution of a square, simply supported, symmetric cross-ply [0/90/0] laminated plate under a bisinusoidally distributed transverse load q, was provided.As shown in Figure 2, the span and thickness of the laminate is denoted by "a" and "h", and q can be described as sin( / )sin( / ) where q0 is a constant.A laminate of span-to-thickness ratio S = a/h = 4 are considered, and can be considered as a thick laminate.The thickness of the ply is 0.25 mm when a = 3 mm.The following lamina material properties are used: E1 = 25 × 10 6 psi, E2 = E3 = 10 6 psi, G12 = G23 = 0.2 × 10 6 psi, ν12 = ν13 = ν23 = 0.25.The resulting stresses are normalized according to the following formulae: The interlaminar stresses analyses were performed using the commercial FE code, Abaqus (Version 6.11, Dassault Systemes Simulia Corp., Providence, RI, USA, 2011).A 3D 8-node linear brick, with an incompatible modes element, was used to model the ply.The interface element based on the linear traction-separation description was introduced into two adjacent plies.This constitutive response is readily available in Abaqus and therefore does not require user-defined subroutines.This 3D FE model is generated using uniform meshes, where the in-plane mesh consisted of 30 × 30 elements.For comparison purposes, two kinds of mesh density in the thickness directions are considered, namely coarse mesh and fine mesh, respectively.For the case, a = 3 mm, S = 4, the former of which is just one element for each ply in the thickness direction, and the latter of which indicates four elements for each ply.According to the recommended rage of the interface stiffness by Reference [38], K is taken as 5 × 10 6 N/mm 3 .The distribution of various non-dimensionalized stresses τ xz , τ yz respectively, through the thickness of the plate at the The interlaminar stresses analyses were performed using the commercial FE code, Abaqus (Version 6.11, Dassault Systemes Simulia Corp., Providence, RI, USA, 2011).A 3D 8-node linear brick, with an incompatible modes element, was used to model the ply.The interface element based on the linear traction-separation description was introduced into two adjacent plies.This constitutive response is readily available in Abaqus and therefore does not require user-defined subroutines.This 3D FE model is generated using uniform meshes, where the in-plane mesh consisted of 30 × 30 elements.For comparison purposes, two kinds of mesh density in the thickness directions are considered, namely coarse mesh and fine mesh, respectively.For the case, a = 3 mm, S = 4, the former of which is just one element for each ply in the thickness direction, and the latter of which indicates four elements for each ply.According to the recommended rage of the interface stiffness by Reference [38], K is taken as 5 × 10 6 N/mm 3 .The distribution of various non-dimensionalized stresses τ xz , τ yz respectively, through the thickness of the plate at the position where each stress attained a maximum is shown in Figure 3.In this procedure, the interlaminar stresses were calculated by the interface element, which used the traction-separation description.The transverse stresses in layers were analyzed based on 3D elasticity and extrapolated from the element integration point.position where each stress attained a maximum is shown in Figure 3.In this procedure, the interlaminar stresses were calculated by the interface element, which used the traction-separation description.The transverse stresses in layers were analyzed based on 3D elasticity and extrapolated from the element integration point.As shown in Figure 3, for the coarse mesh, the percentage error of τ xz at the point of interface, compared with the exact solution is 4.3%, and τ yz is 2.9%.However, the percentage error at the points of maximum value of τ xz and τ yz compared with the exact solution is 21%, 20.9%, respectively.Refined mesh in thickness direction can improve the precision of the transverse stress calculation precision within the ply, for the maximum value of τ xz and τ yz , the percentage error compared with the exact solution is 3.1% and 5.3%, respectively, and the stresses obtained at interface are very close to the exact solution.In summary, for the transverse stresses at interface (namely interlaminar stresses), both of these meshing strategies showed a reasonable accuracy (within 5%).It was noted that the stress distributions within the ply predicted by the coarse mesh model showed considerable error for thick plate.As the transverse stresses at the interface are only considered in this paper, the modeling strategy with coarse mesh in thickness direction is still used in this paper.Furthermore, the interlaminar shear stress τ xz and τ yz distribution at the 0/90 interface without additional normalization based on the coarse mesh strategy is shown in Figure 4.This stress distribution states are consistent with the results in Reference [20].As shown in Figure 3, for the coarse mesh, the percentage error of τ xz at the point of interface, compared with the exact solution is 4.3%, and τ yz is 2.9%.However, the percentage error at the points of maximum value of τ xz and τ yz compared with the exact solution is 21%, 20.9%, respectively.Refined mesh in thickness direction can improve the precision of the transverse stress calculation precision within the ply, for the maximum value of τ xz and τ yz , the percentage error compared with the exact solution is 3.1% and 5.3%, respectively, and the stresses obtained at interface are very close to the exact solution.In summary, for the transverse stresses at interface (namely interlaminar stresses), both of these meshing strategies showed a reasonable accuracy (within 5%).It was noted that the stress distributions within the ply predicted by the coarse mesh model showed considerable error for thick plate.As the transverse stresses at the interface are only considered in this paper, the modeling strategy with coarse mesh in thickness direction is still used in this paper.Furthermore, the interlaminar shear stress τ xz and τ yz distribution at the 0/90 interface without additional normalization based on the coarse mesh strategy is shown in Figure 4.This stress distribution states are consistent with the results in Reference [20].

FEM of Composite Laminates with a Bolt-Filled Hole
A filled-hole composite laminate was selected to analyze the interlaminar stress distributions within the laminate near the bolt.The geometric model is designed according to the geometry of the filled-hole specimens in References [16,39].Henceforth in this paper, "filled-hole laminate" will be used to refer to the laminate with tightened bolts, and the specimens with untightened bolts will be special notes.
A 3D FE model was developed in Abaqus to predict the interlaminar stress distributions near the hole edge of the filled-hole composite laminate.This 3D FE model consists of the bolt, laminate, and the interface based on the linear elastic traction-separation behavior.The bolt and nut simplified as a whole in the FE model is shown in Figure 5.The stacking sequence of the laminate is [02/45/90/−45/0/45/0/−45/0]s.One linear solid element per ply in the thickness direction was used to model the laminate, and the interface element was introduced between plies.The bolt was defined as isotropy material, and linear solid element was also employed.As high interlaminar stress concentrations were expected near the hole, a high mesh refinement was used in this area, as seen in Figure 5.The Coulomb friction model was adopted for the contact between the laminate and the bolt.Based on previous investigations [40], the friction coefficients between bolt shank and laminate is μ = 0.1, and between bolt head and laminate is μ = 0.43.The boundary and loading conditions are

FEM of Composite Laminates with a Bolt-Filled Hole
A filled-hole composite laminate was selected to analyze the interlaminar stress distributions within the laminate near the bolt.The geometric model is designed according to the geometry of the filled-hole specimens in References [16,39].Henceforth in this paper, "filled-hole laminate" will be used to refer to the laminate with tightened bolts, and the specimens with untightened bolts will be special notes.
A 3D FE model was developed in Abaqus to predict the interlaminar stress distributions near the hole edge of the filled-hole composite laminate.This 3D FE model consists of the bolt, laminate, and the interface based on the linear elastic traction-separation behavior.The bolt and nut simplified as a whole in the FE model is shown in Figure 5.The stacking sequence of the laminate is [0 2 /45/90/−45/0/45/0/−45/0]s. One linear solid element per ply in the thickness direction was used to model the laminate, and the interface element was introduced between plies.The bolt was defined as isotropy material, and linear solid element was also employed.As high interlaminar stress concentrations were expected near the hole, a high mesh refinement was used in this area, as seen in Figure 5.

FEM of Composite Laminates with a Bolt-Filled Hole
A filled-hole composite laminate was selected to analyze the interlaminar stress distributions within the laminate near the bolt.The geometric model is designed according to the geometry of the filled-hole specimens in References [16,39].Henceforth in this paper, "filled-hole laminate" will be used to refer to the laminate with tightened bolts, and the specimens with untightened bolts will be special notes.
A 3D FE model was developed in Abaqus to predict the interlaminar stress distributions near the hole edge of the filled-hole composite laminate.This 3D FE model consists of the bolt, laminate, and the interface based on the linear elastic traction-separation behavior.The bolt and nut simplified as a whole in the FE model is shown in Figure 5.The stacking sequence of the laminate is [02/45/90/−45/0/45/0/−45/0]s.One linear solid element per ply in the thickness direction was used to model the laminate, and the interface element was introduced between plies.The bolt was defined as isotropy material, and linear solid element was also employed.As high interlaminar stress concentrations were expected near the hole, a high mesh refinement was used in this area, as seen in Figure 5.The Coulomb friction model was adopted for the contact between the laminate and the bolt.Based on previous investigations [40], the friction coefficients between bolt shank and laminate is μ = 0.1, and between bolt head and laminate is μ = 0.43.The boundary and loading conditions are The Coulomb friction model was adopted for the contact between the laminate and the bolt.Based on previous investigations [40], the friction coefficients between bolt shank and laminate is µ = 0.1, and between bolt head and laminate is µ = 0.43.The boundary and loading conditions are also shown in Figure 5.The in-plane surface traction load P = 100 N/mm 2 are applied to the both ends of the model.The bolt load is defined in terms of a concentrated force, and different pre-load F 0 are applied on the bolt.

Interlaminar Stress Distributions in the Filled-Hole Laminate near the Hole
A numerical solution of interlaminar stress distributions at the interface in the filled-hole laminate near the hole is shown in Figure 6.Analyses were performed in the condition of two levels of clamp-up load, namely non-tightened bolt case and tightened bolt case.In the former case, the clamping force is assumed as F 0 = 10 N, the latter as F 0 = 2000 N. The interlaminar shear stress τ yz is not present in this paper because of the similar distribution feature of τ xz .The following remarks can be highlighted: 1.
For the non-tightened bolt case (clamping force F 0 = 10 N), interlaminar stress concentrations mainly occur around the hole edge (Figure 6a,b, Position I), and tapers out at the interior of the laminate.This prediction is consistent with the previous phenomenon that laminated composites often exhibit transverse stress concentrations near material and geometric discontinuities (the so-called free-edge effect) [1].In spite of the bolt, Position I can still be considered as geometric discontinuities, and the free-edge effects are observed here.

2.
For the tightened bolt case (clamping force F 0 = 2000 N), the results show that the interlaminar stress concentrations in the filled-hole laminate occur at Position I, as well as the interior of the laminate near the periphery of the bolt head (Figure 6c,d, Position II).Note that Position II of the laminate is not a free edge.However, traditionally, the free edge effect is only a concern near the free edge of the laminate, and decay to zero as the distance from the free edge increases.σ z is a negative value in the clamp-up region, in other words, σ z is interlaminar compressive stress.According to Reference [41], the compression interlaminar stress is able to delay the delamination initiation.

3.
It should also be noted that the stress concentrations of σ z and τ xz near Position II are not in the same circumference locations.The former at Position II is shown in Figure 6c,d.The latter is located near Position II, which is away from the hole in the in-plane direction, this position in the laminate is defined as Position III, as presented in Figure 6d.
Appl.Sci.2017, 7, 93 7 of 16 also shown in Figure 5.The in-plane surface traction load P = 100 N/mm 2 are applied to the both ends of the model.The bolt load is defined in terms of a concentrated force, and different pre-load F0 are applied on the bolt.

Interlaminar Stress Distributions in the Filled-Hole Laminate near the Hole
A numerical solution of interlaminar stress distributions at the interface in the filled-hole laminate near the hole is shown in Figure 6.Analyses were performed in the condition of two levels of clamp-up load, namely non-tightened bolt case and tightened bolt case.In the former case, the clamping force is assumed as F0 = 10 N, the latter as F0 = 2000 N. The interlaminar shear stress τ yz is not present in this paper because of the similar distribution feature of τ xz .The following remarks can be highlighted: 1.For the non-tightened bolt case (clamping force F0 = 10 N), interlaminar stress concentrations mainly occur around the hole edge (Figure 6a,b, Position I), and tapers out at the interior of the laminate.This prediction is consistent with the previous phenomenon that laminated composites often exhibit transverse stress concentrations near material and geometric discontinuities (the so-called free-edge effect) [1].In spite of the bolt, Position I can still be considered as geometric discontinuities, and the free-edge effects are observed here.2. For the tightened bolt case (clamping force F0 = 2000 N), the results show that the interlaminar stress concentrations in the filled-hole laminate occur at Position I, as well as the interior of the laminate near the periphery of the bolt head (Figure 6c,d, Position II).Note that Position II of the laminate is not a free edge.However, traditionally, the free edge effect is only a concern near the free edge of the laminate, and decay to zero as the distance from the free edge increases.σ z is a negative value in the clamp-up region, in other words, σ z is interlaminar compressive stress.According to Reference [41], the compression interlaminar stress is able to delay the delamination initiation.3. It should also be noted that the stress concentrations of σ z and τ xz near Position II are not in the same circumference locations.The former at Position II is shown in Figures 6c,d.The latter is located near Position II, which is away from the hole in the in-plane direction, this position in the laminate is defined as Position III, as presented in Figure 6d.

Stress Distributions at Each Interface around the Hole Edge and the Periphery of Bolt Head
According to our numerical calculations, interlaminar stress concentrations occurred at Position I, as well as the interior of the laminate near Position II.The interlaminar shear stress distributions at each interface around Position I with the tightened bolt (F 0 = 2000 N) are shown in Figure 7.The interlaminar normal stress is not presented in this subsection, due to the assumption of that interlaminar compressive stress is able to delay delamination.This assumption is in agreement with the view of Reference [41].In addition, shear stress distributions play a significant role in determining the mechanical behavior of multi-direction laminates.Therefore, we only present the shear stress distributions here.Furthermore, the data points of the stress distribution curve are the absolute values of the magnitude of the interlaminar shear stresses.

Stress Distributions at Each Interface around the Hole Edge and the Periphery of Bolt Head
According to our numerical calculations, interlaminar stress concentrations occurred at Position I, as well as the interior of the laminate near Position II.The interlaminar shear stress distributions at each interface around Position I with the tightened bolt (F0 = 2000 N) are shown in Figure 7.The interlaminar normal stress is not presented in this subsection, due to the assumption of that interlaminar compressive stress is able to delay delamination.This assumption is in agreement with the view of Reference [41].In addition, shear distributions play a significant role in determining the mechanical behavior of multi-direction laminates.Therefore, we only present the shear stress distributions here.Furthermore, the data points of the stress distribution curve are the absolute values of the magnitude of the interlaminar shear stresses.As seen in Figure 7, the magnitude of interlaminar shear stresses at the 0/0 interface are lower than those of other interfaces in most of the area around Position I.These 0/0 interfaces not only refer to the middle interface, but also the outer interface.The value of the interlaminar shear  As seen in Figure 7, the magnitude of interlaminar shear stresses at the 0/0 interface are lower than those of other interfaces in most of the area around Position I.These 0/0 interfaces not only refer to the middle interface, but also the outer interface.The value of the interlaminar shear stresses of these 0/0 interface do not exceed 2 MPa in most of the area around the hole edge.In contrast, the outer 0/45 interface or the middle −45/0 interface exhibits a higher magnitude of interlaminar shear stresses.In general, a higher shear stresses occurs at all non 0/0 interface at Position I.This phenomenon is similar to the classical problems of interlaminar stress that exist at the free edge of laminate.The two classical problems of interlaminar stress are included in Reference [42].[±θ] laminates exhibit only shear-extension coupling, and [0/90] laminates exhibit only a Poisson mismatch between layers.Both of these two problems exist at the 0/45 interface where adjacent layers of the laminate are subjected to an axial load, so the 0/45 interface exhibits obvious interlaminar shear stress concentrations.In addition, due to a match in the elastic properties between the 0/0 interface adjacent layers, the interlaminar shear stress at 0/0 interface is lower than at other interfaces.
The interlaminar shear stress distributions curve at Position III is shown in Figure 8.As seen in Figure 8a, the magnitude of τ xz in the outer interface of the laminate is higher than that of the middle interfaces at Position III, except near θ = 90 • and −90 • , and there is a decrease in thickness direction from outside to inside.τ yz shows a similar tendency excepted near θ = 0 • .A larger magnitude of the interlaminar shear stress τ yz is also exhibited at the 0/45 interface.Therefore, the maximum value of τ xz or τ yz at the outer 0/0 and 0/45 interfaces are much higher than those of the middle −45/0 and 0/0 interfaces.It is interesting to note that, at Position I (Figure 7), regardless of whether it is the outer interface or the middle, large interlaminar shear stresses at the 0/45 interface and small interlaminar shear stresses at 0/0 interface are found.However, at Position III, the magnitude of interlaminar stresses appears to have little relation to the elastic properties of adjacent layers and exhibit a higher value at the outer interfaces.stresses of these 0/0 interface do not exceed 2 MPa in most of the area around the hole edge.In contrast, the outer 0/45 interface or the middle −45/0 interface exhibits a higher magnitude of interlaminar shear stresses.In general, a higher shear stresses occurs at all non 0/0 interface at Position I.This phenomenon is similar to the classical problems of interlaminar stress that exist at the free edge of laminate.The two classical problems of interlaminar stress are included in Reference [42].
[±θ] laminates exhibit only shear-extension coupling, and [0/90] laminates exhibit only a Poisson mismatch between layers.Both of these two problems exist at the 0/45 interface where adjacent layers of the laminate are subjected to an axial load, so the 0/45 interface exhibits obvious interlaminar shear stress concentrations.In addition, due to a match in the elastic properties between the 0/0 interface adjacent layers, the interlaminar shear stress at 0/0 interface is lower than at other interfaces.
The interlaminar shear stress distributions curve at Position III is shown in Figure 8.As seen in Figure 8a, the magnitude of τ xz in the outer interface of the laminate is higher than that of the middle interfaces at Position III, except near θ = 90° and −90°, and there is a decrease in thickness direction from outside to inside.τ yz shows a similar tendency excepted near θ = 0°.A larger magnitude of the interlaminar shear stress τ yz is also exhibited at the 0/45 interface.Therefore, the maximum value of z τ x or τ yz at the outer 0/0 and 0/45 interfaces are much higher than those of the middle −45/0 and 0/0 interfaces.It is interesting to note that, at Position I (Figure 7), regardless of whether it is the outer interface or the middle, large interlaminar shear stresses at the 0/45 interface and small interlaminar shear stresses at 0/0 interface are found.However, at Position III, the magnitude of interlaminar stresses appears to have little relation to the elastic properties of adjacent layers and exhibit a higher value at the outer interfaces.Generally, the plies elastic properties will have a momentous effect on the interlaminar stresses.However, for the laminate with tightened bolts, bolt clamp-up load seems to play a more important role for interlaminar stress distributions in the laminates near the bolt head, than in the Generally, the plies elastic properties will have a momentous effect on the interlaminar stresses.However, for the laminate with tightened bolts, bolt clamp-up load seems to play a more important role for interlaminar stress distributions in the laminates near the bolt head, than in the ply properties.To further understand the clamp-up load F 0 effect on the stress distribution around Position I and Position III, numerical calculation was performed at different clamping forces.The interlaminar stress distributions in the 0/45 interface at Position I and Position III are only given in this subsection, due to the 0/45 interface exhibiting a high magnitude of the interlaminar shear stress in this position.The variation of the interlaminar stresses around Position I at the 0/45 interface are presented in Figure 9, and the interlaminar stresses around Position III at the 0/45 interface are shown in Figure 10.

Stress Distributions at the 0/45 Interface near the Hole under Different Clamp-Up Load
Generally, the plies elastic properties will have a momentous effect on the interlaminar stresses.However, for the laminate with tightened bolts, bolt clamp-up load seems to play a more important role for interlaminar stress distributions in the laminates near the bolt head, than in the ply properties.To further understand the clamp-up load F0 effect on the stress distribution around Position I and Position III, numerical calculation was performed at different clamping forces.The interlaminar stress distributions in the 0/45 interface at Position I and Position III are only given in this subsection, due to the 0/45 interface exhibiting a high magnitude of the interlaminar shear stress in this position.The variation of the interlaminar stresses around Position I at the 0/45 interface are presented in Figure 9, and the interlaminar stresses around Position III at the 0/45 interface are shown in Figure 10.As seen in Figure 9,  The interlaminar normal stress σ z rises with increased clamp-up load in all regions of Position I at the 0/45 interface (Figure 9a).In addition, σ z also shows interlaminar compressive stress at Position I for the case of F0 = 6000 N.  A higher clamping force results in lower interlaminar shear stresses for most of the region around Position I (Figure 9b,c).This area is in the region around −80° < θ <−10° and 10° < θ < 55°.For the other areas around the hole, such as the region near θ = 0°, the interlaminar shear stresses almost become zero, regardless of the kind of bolt clamp-up load mentioned in this paper.
It should be noted that the above two stress distribution phenomena are conductive to the improvement of the hole edge strength.
Correspondingly, for the interlaminar stresses around Position III at the 0/45 interface as depicted in Figure 10:  For the clamping force F0 = 10 N, both the interlaminar normal and shear stresses are almost zero in most areas at Position III, except τ xz in region −50° < θ < 45°.This is consistent with the results obtained in Section 5.1. With the increase in clamp-up load, both interlaminar normal stress (Figure 10a) and interlaminar shear stresses (Figure 10b,c) increase in most of the circular region.This is As seen in Figure 9,

•
The interlaminar normal stress σ z rises with increased clamp-up load in all regions of Position I at the 0/45 interface (Figure 9a).In addition, σ z also shows interlaminar compressive stress at Position I for the case of F 0 = 6000 N.

•
A higher clamping force results in lower interlaminar shear stresses for most of the region around Position I (Figure 9b,c).This area is in the region around −80 • < θ <−10 • and 10 For the other areas around the hole, such as the region near θ = 0 • , the interlaminar shear stresses almost become zero, regardless of the kind of bolt clamp-up load mentioned in this paper.
It should be noted that the above two stress distribution phenomena are conductive to the improvement of the hole edge strength.
Correspondingly, for the interlaminar stresses around Position III at the 0/45 interface as depicted in Figure 10:

•
For the clamping force F 0 = 10 N, both the interlaminar normal and shear stresses are almost zero in most areas at Position III, except τ xz in region −50 • < θ < 45 • .This is consistent with the results obtained in Section 5.1.

•
With the increase in clamp-up load, both interlaminar normal stress (Figure 10a) and interlaminar shear stresses (Figure 10b,c) increase in most of the circular region.This is different from the interlaminar stress distributions in the laminate at Position I.At Position III, the phenomenon of the interlaminar shear stress reduction does not appear with the increase in clamp-up load.
In contrast, shear stress increases with the increase of the clamping force.
As a result, the bolt clamp-up load cannot ameliorate the interlaminar stress distributions at the interface near Position II.The resulting phenomenon in this subsection can be summarized as follows: First, increasing the clamping force can reduce the magnitude of the interlaminar shear stresses in most regions around the hole.Thus, improvements can be made in the ability of the hole edge of the laminate to resist delamination.Second, the bolt clamp-up force cannot improve the interlaminar stress distribution near Position II.The former is consistent with the understanding of Reference [43].For the latter, this paper will further analyze its causation.

Discussion
Interlaminar stresses exist within the laminate that cannot be directly measured by experiment.As an alternative, the interlaminar stress distribution characteristics of filled-hole compression tests can be identified by the failure mode.Castanie et al. [16] analyzed failure modes occurring during the filled-hole compression (FHC) test, and observed that the crack initiated at the periphery of the bolt with a little off-set from the net section for high 0° oriented lay-ups, and the final failure mode was a static failure of the outer ply oriented at 0° (Figure 11).They also supposed that the reason for this phenomenon was that the tightening modified the local triaxial stress concentration and load transfer between the hole and the bolt.The interlaminar stress distributions trend in the filled-hole laminate near the hole that was obtained in this paper is in accordance with Castanie's inference.The numerical analysis in this paper shows that the tightened bolt does change the interlaminar stress distribution state in the filled-hole laminate, and both Position I and Position II exhibit interlaminar stress concentration.Furthermore, the numerical analysis also shows that initial clamping forces can improve the interlaminar stress distribution state at Position I.However, the bolt clamping force did not make the interlaminar stress distribution better in the laminate near Position II, according to the following two inferences made in this paper.First, the interlaminar shear stresses near Position II increased with an increase in the clamp-up load.Second, the The resulting phenomenon in this subsection can be summarized as follows: First, increasing the clamping force can reduce the magnitude of the interlaminar shear stresses in most regions around the hole.Thus, improvements can be made in the ability of the hole edge of the laminate to resist delamination.Second, the bolt clamp-up force cannot improve the interlaminar stress distribution near Position II.The former is consistent with the understanding of Reference [43].For the latter, this paper will further analyze its causation.

Discussion
Interlaminar stresses exist within the laminate that cannot be directly measured by experiment.As an alternative, the interlaminar stress distribution characteristics of filled-hole compression tests can be identified by the failure mode.Castanie et al. [16] analyzed failure modes occurring during the filled-hole compression (FHC) test, and observed that the crack initiated at the periphery of the bolt with a little off-set from the net section for high 0 • oriented lay-ups, and the final failure mode was a static failure of the outer ply oriented at 0 • (Figure 11).They also supposed that the reason for this phenomenon was that the tightening modified the local triaxial stress concentration and load transfer between the hole and the bolt.The interlaminar stress distributions trend in the filled-hole laminate near the hole that was obtained in this paper is in accordance with Castanie's inference.The numerical analysis in this paper shows that the tightened bolt does change the interlaminar stress distribution state in the filled-hole laminate, and both Position I and Position II exhibit interlaminar stress concentration.Furthermore, the numerical analysis also shows that initial clamping forces can improve the interlaminar stress distribution state at Position I.However, the bolt clamping force did not make the interlaminar stress distribution better in the laminate near Position II, according to the following two inferences made in this paper.First, the interlaminar shear stresses near Position II increased with an increase in the clamp-up load.Second, the compressive interlaminar normal stress and the interlaminar shear stresses are not at the same circular position.These interlaminar stress distribution characteristics may induce delamination in composite joints, and the delamination position may be initiated near the bolt head prior to the hole edge, which leads to ultimate delamination failure modes offset from center of hole.According to the mechanism of interlaminar stresses, there are two most important factors for determining the presence and magnitude of interlaminar stresses.First, when load is applied to a laminate, different layers in the laminate tend to slide over each other due to the differences in their elastic constants or load boundary.Interlaminar stresses are produced along with the process of deformation, for example, the mismatches between the Poisson's ratios of adjacent plies at the free edges.These load conditions may be uniaxial load, in-plane shear, out-of-plane shear/bending, or in-plane bending.Second, a transverse load applied to a laminate along the thickness direction will create transverse stresses, especially σ z for filled-hole composite laminates, that were subjected to uniaxial compression and clamp-up load.To understand the distribution of transverse stresses within a filled-hole composite laminates under compression, Figure 12 depicts the deformation that takes place in that laminate.As seen in Figure 12, the laminate is divided into three regions that are denoted as Regions A, B, and C. Region A is the clamp-up region, Region C is the non clamp-up region, and Region B is the transition region of deformation.There are some differences between the normal deformations of the laminate in these three regions.In Region A, normal deformation near the laminate hole was induced by the bolt clamp-up load.In Region C, the symmetrical laminate was only subjected to in-plane load, so normal deformation was mainly caused by the Poisson effect of the composite laminates.In Region B, normal deformation was the result of multiple causes, including the clamp-up load, the Poisson effect, and the couplings between extension and bending/twisting deformations of composite laminates.Considering that laminates are mainly subjected to in-plane compressive load, the clamp-up load plays a more important role in normal deformation than the Poisson and coupling effect, which will lead to normal deformation in Region A that is more obvious than that in Region C. Thus, normal deformation of the laminate has relatively higher gradients in Region B. In Region B, normal deformation decreases gradually both in the thickness direction from outside to inside and in the in-plane direction away from the bolt-head edge.The difference in normal deformation results in the interaction between different lamina at Region B, followed by the rise in magnitude of interlaminar stresses in this area, especially the interlaminar shear stress.It should be noted that Region B is not right under the bolt head, and is near Position II, away from the hole.According to the mechanism of interlaminar stresses, there are two most important factors for determining the presence and magnitude of interlaminar stresses.First, when load is applied to a laminate, different layers in the laminate tend to slide over each other due to the differences in their elastic constants or load boundary.Interlaminar stresses are produced along with the process of deformation, for example, the mismatches between the Poisson's ratios of adjacent plies at the free edges.These load conditions may be uniaxial load, in-plane shear, out-of-plane shear/bending, or in-plane bending.Second, a transverse load applied to a laminate along the thickness direction will create transverse stresses, especially σ z for filled-hole composite laminates, that were subjected to uniaxial compression and clamp-up load.To understand the distribution of transverse stresses within a filled-hole composite laminates under compression, Figure 12 depicts the deformation that takes place in that laminate.As seen in Figure 12, the laminate is divided into three regions that are denoted as Regions A, B, and C. Region A is the clamp-up region, Region C is the non clamp-up region, and Region B is the transition region of deformation.There are some differences between the normal deformations of the laminate in these three regions.In Region A, normal deformation near the laminate hole was induced by the bolt clamp-up load.In Region C, the symmetrical laminate was only subjected to in-plane load, so normal deformation was mainly caused by the Poisson effect of the composite laminates.In Region B, normal deformation was the result of multiple causes, including the clamp-up load, the Poisson effect, and the couplings between extension and bending/twisting deformations of composite laminates.Considering that laminates are mainly subjected to in-plane compressive load, the clamp-up load plays a more important role in normal deformation than the Poisson and coupling effect, which will lead to normal deformation in Region A that is more obvious than that in Region C. Thus, normal deformation of the laminate has relatively higher gradients in Region B. In Region B, normal deformation decreases gradually both in the thickness direction from outside to inside and in the in-plane direction away from the bolt-head edge.The difference in normal deformation results in the interaction between different lamina at Region B, followed by the rise in magnitude of interlaminar stresses in this area, especially the interlaminar shear stress.It should be noted that Region B is not right under the bolt head, and is near Position II, away from the hole.
in Region A that is more obvious than that in Region C. Thus, normal deformation of the laminate has relatively higher gradients in Region B. In Region B, normal deformation decreases gradually both in the thickness direction from outside to inside and in the in-plane direction away from the bolt-head edge.The difference in normal deformation results in the interaction between different lamina at Region B, followed by the rise in magnitude of interlaminar stresses in this area, especially the interlaminar shear stress.It should be noted that Region B is not right under the bolt head, and is near Position II, away from the hole.Based on the above analysis, for interlaminar stress distributions in the laminate around Position II, the dominant factors that cause interlaminar normal stress and interlaminar shear stress are different.Thus, the stress concentration of interlaminar normal stress and interlaminar shear stress near Position II are not in the same circumference location.The former is mainly caused by the normal deformation in Region A, hence the appearance of interlaminar normal stress in Region A. The latter originates from the difference of normal deformation in Region B; therefore, the shear stress concentration mainly occurs in the outer interfaces in Region B. The larger the clamp-up load, the more obvious the normal deformation in Region B, and the interlaminar shear stresses increase with the increase of the clamp-up load accordingly.

Conclusions
The interlaminar stress distributions in a filled-hole laminate with clamp-up load near the hole, especially near the bolt head, were investigated.To evaluate interlaminar stresses under the complex boundary condition efficiently, a calculation strategy using zero-thickness interface element was presented and validated.The interface element was based on a linear elastic traction-separation description.The resulting analysis is conducive to an understanding of the geometry effects and the failure mechanisms of filled-hole compression.Furthermore, it can provide advice for the appropriate design of mechanically fastened joints in a composite structure.According to the numerical simulation for filled-hole compression (FHC) test, the following conclusions were reached.

1.
For the filled-hole laminate with tightened clamp-up load, the tightened bolt will change the interlaminar stress distribution state in the filled-hole laminate.The interlaminar stress concentration occurs both at the hole edge and the laminate interior near the outer edge of the bolt head.This phenomenon is different from the traditional free edge effect of an open-hole laminate.Correspondingly, for the case of non-tightened bolts, the interlaminar stress concentration only occurs near the hole edge.

2.
The interlaminar normal and shear stress concentrations in the laminate near the periphery of the bolt head are not distributed in the same circumference location.The former is located under the periphery of the bolt head, whereas the latter occurs near the bolt-head edge away from the hole in the in-plane direction, which is at the outer interface of the filled-hole laminate in the thickness direction.

3.
The clamp-up load plays a more important role in the interlaminar stress distribution in the laminate near the periphery of the bolt head.With an increase in the clamp-up load, both the interlaminar normal stress and interlaminar shear stresses rise around the bolt head in a large region.However, the magnitude of interlaminar shear stresses is reduced for most of the region around the hole.
According to the above conductions, the clamp-up load cannot improve the interlaminar stress distribution of the filled-hole laminate near the bolt-head edge, which may lead to the initiation of delamination in bolt-filled laminates near the periphery of the bolt head.However, the specific location of this failure mode is still to be investigated by numerical and experimental analysis.

FFigure 5 .
Figure 5. Finite element mesh and modeling strategy for a filled hole compression (FHC) specimen, and showing boundary conditions.

FFigure 5 .
Figure 5. Finite element mesh and modeling strategy for a filled hole compression (FHC) specimen, and showing boundary conditions.As shown in Figure 5, the laminate dimensions are L = 32 mm by W = 32 mm, the hole of the laminate has a diameter of D = 6.35 mm, and the thickness of the lamina is t = 0.25 mm.The bolt has a shank diameter Dbolt = 6.35 mm, and head diameter d = 11.3 mm.Elastic, isotropic material properties are used for the bolt section, with E = 109 GPa and ν = 0.3.The laminate is made of material plies that are idealized as homogeneous, elastic and orthothropic.The following material properties are used: E11 = 181 GPa, E22 = E33 = 10.3GPa, G12 = G13 = 7.17 GPa, G23 = 4.13 GPa, and ν12 = ν13 = ν23 = 0.25.Subscripts 1, 2, and 3 denote the fiber, transverse, and thickness directions, respectively.The Coulomb friction model was adopted for the contact between the laminate and the bolt.Based on previous investigations[40], the friction coefficients between bolt shank and laminate is μ = 0.1, and between bolt head and laminate is μ = 0.43.The boundary and loading conditions are

Figure 5 .
Figure 5. Finite element mesh and modeling strategy for a filled hole compression (FHC) specimen, and showing boundary conditions.

Figure 6 .
Figure 6.Stress distributions at the interface in the filled-hole laminate near the hole: (a) Interlaminar normal stress σ z , with non-tightened bolt; (b) Interlaminar shear stress τ xz , with non-tightened bolt; (c) σ z , with tightened bolt; (d) τ xz , with tightened bolt.1-Bolt; 2-Interface; 3-Quarter clamp-up region; I-The hole edge; II-The projection boundary of the periphery of the bolt head on the laminate.

3 -
Quarter clamp-up region; I-The hole edge; II-The projection boundary of the periphery of the bolt head on the laminate.

Figure 7 .
Figure 7. Interlaminar shear stress distributions around Position I: (a) τ xz ; (b) τ yz .Outer interface-The interfaces located outside the filled-hole laminate, such as the 0/0 and 0/45 interface.Middle interface-The interfaces near the symmetry plane of laminate, such as the −45/0 and 0/0 interface.
l e  , d e g

Figure 7 .
Figure 7. Interlaminar shear stress distributions around Position I: (a) τ xz ; (b) τ yz .Outer interface-The interfaces located outside the filled-hole laminate, such as the 0/0 and 0/45 interface.Middle interface-The interfaces near the symmetry plane of laminate, such as the −45/0 and 0/0 interface.