A simplified mesoscale model for predicting the mechanical behavior of stitched CFRP laminates

This paper presents a finite element (FE) analysis of the mechanical responses for stitched CFRP laminates under different mechanical loads. Firstly, the through-thickness stitch was simplified to z-pin like reinforcement with a uniform displacement constraint on the upper and lower surfaces of the laminate. Then, a mesoscale 3D representative volume element (RVE) of the stitched composite was proposed and modeled in the FE code ABAQUS, where the reinforcing stitch, composite layers and interfaces were built. A 3D Hashin damage model and built-in cohesive elements were respectively used to predict the mechanical failure of the stitch and the damageable behavior of cohesive interfaces. Simulation results reveal the progressive damage and rupture processes of the RVE under tensile and shear mechanical loads, and macroscopic nonlinear load-displacement responses of the mesoscale model are also captured.


Introduction
Carbon fiber reinforced plastic (CFRP) laminates receive wide use in the aerospace industry. However, their usage is affected by their low delamination resistant ability, and consequently the impact damage and the manufacture defect between laminate layers becomes the main safety concern. Adding throughthe-thickness stitches or z-pins [1,2], although with different strengthening mechanisms, is the effective way to improve the delamination fracture toughness of the composite laminates. In the stitched laminates, the stitches are continuous and constantly stretched, and when delamination happens bridging load provided by the stitch will increase with the delamination displacement until the tension stress in the stitch reaches its ultimate strength, then the ruptured stitch can generates little bridging load during the frictional pullout process. In comparison, the frictional bridging load during z-pin pullout process is relatively larger, which is the main mechanism increasing the mode I delamination toughness of the zpinned laminates.
As the stitching and z-pinning technique become increasingly used in the composite structures, the strength analysis and prediction of these reinforced laminates attracts much attentions along with their application. Many models have been developed in the open literature, which can be roughly categorized into two groups. The first one is to propose theoretical models by analyzing the mechanical process during the tension, bending or shear load of the reinforced structure to obtain closed form solutions [ and the other is to build FE models, simplifying the stitches or z-pins into spring elements [4], connecting elements [5], beam elements [2,6], three-nodes rod elements [7,8] or cohesive zone elements [9,10,11] with material parameters and corresponding load (or traction) versus displacement relationships to predict the mechanical behaviors and load-bearing capacity of the reinforced laminate structures. However, one can't find generally recognized theoretical or FE methods in this research field.
This study aims at proposing a simplified mesoscale FE model for analyzing the stitched laminates, which could be used to predict the strengthening effect of the stitch and its failure modes under different loads, and also capture the macroscopic mechanical responses of the stitched laminate.

Model description
In the paper, we take the typical stitched T-shaped sandwich joints as shown in Fig. 1 for analysis. A representative volume element (RVE) was selected from the stitched area which shows periodical geometrical characteristics. Because the stretched stitch restricts the deformation on the upper and lower surfaces of the RVE, we simplify its geometry into a z-pin like stitch, plus a uniform displacement or strain constraint on these surfaces. Meanwhile, the simplified z-pin like stitch keeps the same cross-section as the original stitch, then the simplified RVE should yield the same through-thickness deformation ΔL under the tensile load F as the original one. Take the cross-section of the RVE for analysis (see Fig. 2), the stitch length of two segments on the upper surface is assumed L 1 and L 2 to respectively, but with only half of the section area of stitch segment in the through-thickness area, which is assumed as A/2 and with tensile modulus E. In the simplified RVE, the total length of the through-thickness stitch is comprised by L 3 and L 4, and its modulus E′ should satisfy the following equation based on the equivalent displacement assumption under the transverse tension load F (see Fig. 3): Then, the tensile modulus of the stitch E′ in the simplified RVE model can be deduced as: In fact, the reinforcing stitch is composed by the Kevlar stitch and resin, and the latter component usually take 60% in the volume, so the longitudinal modulus of the composed stitch should multiply a factor of 0.4 by the modulus provided by the manufactures. In this paper, according to the geometric configuration the stitch length is L 1 = L 2 =2.5mm, and the thickness of the web and flange is L 3 =1.6mm, L 4 = 8.8mm, respectively. As a result, one can obtain E' ≈ 0.51E in the simplified RVE.
The elastic and strength properties of the reinforcing stitch in the simplified RVE are presented in Table 1and Table 2, where the units of modulus and strength parameters are in MPa.    Table 2. Strength properties of the stitch.  Fig. 4. The model is composed by the flange and web along the thickness direction, and an interface between them was introduced for delamination analysis. In the center part of the model, a cylinder represents the reinforcement stitch with a diameter of 0.72mm, and it was surrounded by a 0.1mm thickness resin rich cylinder area. 5

Material model
Under mechanical loads, the stitch firstly sustains elastic stretch and would finally rupture. Hence, the damage and failure behaviors of the stitch should be considered. In this paper, we introduce the wellknown Hashin damage criterion to capture these non-linear mechanical behaviors. In the 3D stress state, they have the following expressions considering the characteristics of the stitch: 1) Stitch tensile failure mode (when 11 0   ) 2) Matrix tensile failure mode ( 22  As the FE analysis code (ABAQUS/Standard) can only provide a 2D Hashin damage material model, a user-defined field subroutine (USDFLD) was programmed to implement the 3D Hashin damage model mentioned above. For simplicity, two state of the material integration point was defined: the elastic state (which is flagged as 0) and failure state (flagged as 1), and three state dependent variables (SDV1, SDV2, and SDV3) were introduced to record the numerical values on the left-hand side of the above three equations, respectively. Besides, two field variables (FV1 and FV2) were introduced with a zero value at the elastic state. When the stitch tensile failure mode is triggered, which means SDV1 reaches and/or exceeds 1.0, FV1 is then set as 1.0. Similarly, when matrix tensile or compressive failure happens (SDV2 or SDV3 ≥ 1.0), FV2 is then set as 1.0.
The related stiffness reduction strategy is listed in Table 3. For numerical reasons, the stiffness at the failure material point is set to 1% of the elastic one.  The ABAQUS/Standard provides cohesive elements capable of simulating the mechanical behaviors of the interface layers between the web and flange, as well as those between the stitch and laminates. Elastic constants of the cohesive elements used in the present analysis are: E nn = 3.0GPa, G ss = G tt = 1.1GPa, and the damage initiation criterion is: where the symbol <> denotes the Heviside operator, t n represents normal tensile stress, t t and t s represents the shear stresses along two perpendicular direction of the cohesive element, and the tension strength Besides, linear damage evolution equation is used to capture the initiation of stiffness change under the mechanical loads. Meanwhile, BK criterion [12] is selected to determine the mixed mode fracture toughness of the cohesive element by setting its power law index η = 1.45, and fracture toughness parameters: C n G = 0.25N/mm and C C s t G G  =0.753N/mm. Finally, the resin area surrounded the stitch can be treated as a homogonous and isotropic material. It's modulus is 3.5GPa, and Poisson's ratio is 0.35.

Finite element discretization
Based on the above geometric and material model, the mesoscale RVE model was built in the FE software. In which an 8-node linear brick element with reduced integration (C3D8R) was selected to discrete the flange and web, as well as the resin area, and an 8-node 3D cohesive element (COH3D8) was used for delamination analysis. The geometrical model was appropriately meshed with a relative finer area around the stitch, see Fig. 5. The total number of the solid elements is 245140. Meanwhile, the viscosity factor is set as 5.0×10 -3 in order to achieve better convergence rate after some simulation tests. As the stitched elements are typically subjected to out-of-plane tension and/or transverse shear loads, two models with different loading conditions are modeled to predict their mechanical responses, respectively.  The transverse shear model adopted the same geometrical model and material parameters as the above tension model, and also the same type of multi-points constraints on both the upper and lower surfaces during the shear loading process, see Fig. 7a. In which one more multi-point constraint surrounded the flange laminate was added to conveniently output the reaction force. The transverse shear load is realized by introducing a pair of transverse displacement along and contrary to the Y-axis, see Fig. 7b.
Similarly, the 3D Hashin damage model and cohesive elements were used to simulate the failure behavior of the stitch and delamination of the structure under the shear load.

Out-of-plane tension
From the simulation results, we obtained the relationship between the reaction tension force and the displacement on the loading surface, as shown in Fig. 8. One can observe that the reaction force increases linearly as the out-of-plane tensile displacement increases at the first stage. However, the loading carry capacity of the RVE cell begins to decrease after the reaction force reaches its maximum value. Fig. 9a reveals that at the maximum load the damage factor of almost all the cohesive elements between the web and flange, i. e. the net section reaches 1.0, which means they are closed to failure and delamination begins to occur. The tensile load in this section gradually transfers to the reinforcing stitch, although it has already experienced some minor damages, so the RVE cell still maintains some residual loadcarrying capacity. As the tensile displacement grows, the stitch damages gradually. The state dependent variable SDV1 contour in the section view of the RVE cell (see Fig. 9b) shows that the stitch section near the delaminated surface mainly suffers tensile damage. When the net section of the stitch finally ruptured, the stiffness of RVE cell reduced to 0.
Besides, it can be deduced from Fig. 8 that the nominal tensile strength of the stitched composite is close to 171.04MPa, which is much higher than the tensile strength property of the cohesive element 65.6MPa. This comparison result indicates the introduced stitch can remarkably promote the delamination resistance of the traditional laminates.

Transverse shear
The Mises stress contour and deformation results of the transversely shear loaded RVE cell at different loading stages are presented in Fig. 10 (given in a symmetric section view for observation purpose). Fig.  10a shows that, at the initial loading stage, the upper and lower plates dislocate under the transverse shear loads, and the stress distribution is not uniform along the loading direction. The cohesive elements near two side-loaded surfaces are subjected to higher stress, and would firstly damage and even failure, then the high stress area gradually transfers to the stitch area (see Fig. 10b). When the shear load reaches its maximum (corresponding to Fig. 10c), the load sheared by the stitch already induces damage and failure in some stitch elements. Finally, the RVE cell ruptures as the effective load area shrinks rapidly.
(a) Cohesive damage state at the maximum load (b) Stitch tensile damage at failure load   11 shows the macroscopic shear load .vs. displacement relationship of the stitched RVE. It should be noted that the computation process is aborted due to convergence difficulties, however, one can still clearly see the load-carry tendency of the RVE cell. The relationship between the transverse shear load and the dislocation displacement can be roughly treated as bilinear, which is similar to the tension simulation results. Besides, the maximum averaged shear stress on the loading surface reaches 57.32MPa, it is much higher than the shear strength of the cohesive element (15.50MPa). This result also reveals that the stitch greatly reinforces the shear stress bear capacity of the laminates.

Conclusion
This paper proposed a simplified z-pin like RVE cell for the through-thickness stitch reinforced laminate structure at the mesoscopic level. The model uses a 3D Hashin damage criteria for damage analysis of the stitch and introduces cohesive elements for delamination prediction. The FE simulation reveals the damage and failure process of different components under different mechanical loads.
Non-linear macroscopic load-displacement relationships are captured for the RVE subjected to outof-plane tension and transverse shear loads. The predicted strengths indicate that the stitch remarkably promote the tension and shear load-carrying capacity compared with the traditional laminates. The proposed mesoscale model is appropriate for microscopic material parameter design, and the predicted strength parameters can be used for structural strength analysis.