Progressive Analysis of Bearing Failure in Pin-Loaded Composite Laminates Using an Elasto-Plastic Damage Model

Bearing failure of composite laminate is very complicated due to the complexity of different failure mechanisms and their interactions. In this paper, an elasto-plastic damage model is built up to describe the process of failure in composite laminates subjected to bearing load. Non-linear behavior of composite before failure is taken into consideration by using a modified Sun-Chen one parameter plasticity model. LaRC05 failure criteria are employed to predict the initiation of failure and the evolution of failure is described by a CDM based stiffness degradation model. Both theory and some application issues like parameter determination are discussed according to phenomenon of experiments. The model is firstly validated by several experiment results of unidirectional laminate and then applicated into the progressive analysis of bearing failure in pin-loaded multidirectional laminates, both intralaminar and interlaminar damage are taken into consideration. The result of finite element analysis is compared with experiment results; it shows good agreements in both mechanical response and progress of failure, so the model can be evaluated to be effective and practical in bearing failure analysis of composite laminates.


Introduction
At present, the application of composite materials in aircrafts has been transited gradually from secondary parts to main parts.In the main parts, the connections between different components are mostly in the form of mechanically fastened was incorporated into this analysis which means that the effect of intralaminar damage was ignored.Shen et al. [8] discussed about two issues which played important roles in the simulation of delamination in composite laminates.The first one was the symmetry of FEA model in delamination analysis and the second one was the prediction of delamination initiation.To use a half or quarter model to simulate delamination in composite laminates, a strict symmetry condition should be satisfied, so quasi-static lay-up laminates cannot be modeled with a simplified half or quarter model.For the prediction of delamination initiation, strain energy components based criteria should be used instead of average strain energy based criteria.
Frizzell and McCarthy [9] developed a 3D FEA model to simulate the failure process of bolted joint in FML laminates.Intralaminar damage was described by CDM model and interlaminar failure was captured by cohesive zone model.A damage regularization method was used to avoid the difficulty of converge.
Although there already have been many significant achievements in this area, a lot of unsolved problems still exist [10] [11] [12].For example, the interaction of different damage mechanisms were not taken into full consideration, this will cause an overestimation of bolted joints strength in composite laminates.The damages in thickness direction were rarely studied in contrast with their important roles in the failure of composite laminates.Hence there is still a lot work to be done.Herein, an elasto-plastic damage model is developed to consider more factors in the failure analysis of bearing failure in composite laminates.

Modified Sun-Chen One Parameter Plasticity Model
Wang et al. [13] proposed a generalized Hill yield criterion which considers the effect of hydrostatic pressure for anisotropic materials and non-symmetry of yield strength under tensile and compressive load by referring to Ducker-Prager yield criterion.The simplified yield function is as follows: where 2 T Y is transverse tensile yield strength, 66 a is a parameter in Sun-Chen model, m is the ratio of tensile yield strength and compressive yield strength.Γ is a parameter which characterizes the non-symmetry of material behavior under tension and compression and ( ) The effective stress can be expressed as follows: ( )( ) following the incremental derivation method of plastic work per unit volume, the relationship between plastic strain increment and effective plastic strain increment can be deduced: Substituting Equation (2) into Equation ( 4), the expression of plastic strain components in incremental form is derived: From Equation ( 4), it can be easily found that plasticity flows in all material directions except the longitudinal direction which is supposed to have a linear mechanical behavior and the flow of plasticity is closely related to current stress state.
The material is assumed to satisfy isotropic hardening rule and the effective stress is related to the effective plastic strain with a power law: where A and n are both parameters which describe the plasticity of material.
According to the initial yield criterion, subsequent yield surface can be obtained: where k is hardening parameter and it associates with plastic deformation which can be expressed as a function of effective plastic strain, then the function of k on effective plastic strain can be deduced: ( ) ( ) ( ) The expression of plastic compliance is as follows: Combining the plastic compliance with the elastic compliance, incremental stress-strain relationship is derived:    11) where: The determination method of the parameters in this model can be obtaining by referring to Wang et al. [14].A series of off-axis compression experiments are needed.The plasticity parameters A and n in Equation ( 6) can be determined by curve fitting method according to the plot of effective stress and effective plastic strain.Then 66 a is obtained by comparing the response of two off-axis com-

Matrix Failure Criteria
Matrix failure is also called inter-fiber failure which means a crack parallel to fiber direction has developed through the whole lamina, including matrix cracks and fiber-matrix interface cracks.LaRC05 criteria [15] is used here to predict the initiation of matrix failure as in Equation ( 14), when 1 m f = , matrix failure will initiate.
( ) ( ) where T τ and L τ are the transverse shear stress and longitudinal shear stress on the potential fracture plane respectively, N σ is the normal stress on the po- tential fracture plane.The stress components on the potential fracture plane can be deduced according to the transformation matrices in Appendix 1.The direction of fracture plane is defined in Figure 1, where α is fracture angle.T S and L S are transverse fracture resistance and longitudinal resistance on poten- tial fracture plane respectively.T µ and L µ are inclination coefficient or fric- tional coefficient, which represent the influence of normal stress on the fracture resistance.
The parameters can be determined using unidirectional laminate off-axis compression experiment.The off-axis compressive strengths of different off-axis angle are notated as

Fiber Failure Criteria
Fiber kinking (Figure 2) can be defined as localized shear deformation in a band (kinking band) when the material subject to compressive load along fiber direction.As reported by literatures [2] [16] [17] [18], fiber kinking plays a key role in the failure of composite laminate during compression, so in this paper kinking is primarily considered as the pattern of fiber failure in composite laminate subject to compressive load.
The process of kinking band formation is very complicated, phenomenons like fiber inclination, matrix shear deformation, inter-fiber failure in kinking band and fiber breakage at the boundary of kinking band can be observed [16].
It's not possible to develop a model which can represent all the factors during kinking band formation, so a simplified model in Figure 3 is built up to analyze the primary mechanisms of failure in the fiber kinking failure process.
Before failure initiation, the material is assumed to be continuum, so stress keeps balanced everywhere in the material.Then stress components can be rotated to local framework of kinking band and used to evaluate if failure will happen.Hence the local framework of kinking band needs to be determined firstly.
It's supposed that kinking plane is always parallel to fiber direction, so the coordinate of kinking plane ψ Ω can be obtained by rotating material coordi- nate by an angle ψ around axe 1 (Figure 3(a)), and the coordinate of kinking band is to be acquired by rotating kinking plane coordinate ψ Ω by an angle ϕ around axe 3 ψ (Figure 3(b)).The angle ψ depends on transverse stress state and is calculated as follows: ( ) Fiber misalignment angle ϕ is the sum of initial misalignment angle i ϕ and shear strain γ , which means   where the superscripts of stress components represent the coordinate of stress.
Applying the matrix failure criteria in Equation ( 14), the fiber misalignment angle when matrix failure happens in kinking band is derived: It should be noted that, i ϕ is not an initial misalignment angle in reality, which means it's not possible to get the value of i ϕ by experiment.i ϕ is ac- tually an effective initial misalignment angle derived from the longitudinal compressive strength C X according to the material's constitutive model.If differ- ent constitutive models are chosen to represent the mechanical behavior of material, different values of i ϕ will be obtained.
After i ϕ is determined, fiber misalignment angle ϕ can be calculated from Equation ( 22), which means the coordinate of kinking band is determined. ( Then stress components in material coordinate are rotated to kinking band coordinate, LaRC05 criteria is applying to judge the initiation of fiber kinking failure.When the value of fc f exceeds one, failure will happen. where For fiber tensile failure, maximum stress criterion is used to predict failure initiation.When the value of ft f exceeds one, failure will happen.

Damage Evolution
If failure initiates, material stiffness is to be degraded to consider the effect of damage evolution.According to literatures [19], degraded material stiffness matrix is in the form as in Equation ( 25): 22 23 33 ed 44 55 66 1 1 where L d and T d are damage variable for fiber failure mode and matrix failure mode respectively. , where , , , α α λ λ are parameter which define the mechanical response of material during the process of damage evolution.α define the form of stress-strain relationship and λ define the intense of degradation as shown in Figure 4.

Strength Prediction
The elasto-plastic damage model developed in Section 2 is firstly validated by two cases of strength prediction.These two cases are provided by literatures [20] [21].The properties and parameters used for numerical analysis are in Table 1.
Case 1: The material system is T300/PR319 (carbon fiber/epoxy).In the test, hydrostatic pressure of 600 MPa is applied to the material and a shear load is applied to the material at the same time.The comparison between the model prediction curve and the test results is shown in Figure 5.The failure stress prediction error is 9.3% and the failure strain prediction error is 16.1%.Considering the dispersiveness of the test results under the complex load condition, it can be considered that the analytical results of the failure analysis method are in good agreement with the experimental results, which proves the effectiveness of the analytical method for the failure strength analysis of the matrix.
Case 2: The material used in the test was S2-glass/epoxy (high-strength glass fiber/epoxy) system; the material received varying levels of lateral hydrostatic pressure while measuring the axial compressive strength of the material.
The results of the comparison between the model prediction results and the test results are shown in Figure 6.It can be seen that the predicted results of the nonlinear strength analysis using the modified Sun-Chen plasticity model agree well with the experimental values.When using the linear model analysis, as the hydrostatic pressure increases, the prediction results will gradually become larger than reality.Both linear and nonlinear analyses can predict the increase of axial compressive strength with increasing lateral hydrostatic pressure.This is determined by the characteristics of the failure criterion, but it can be predicted by nonlinear analysis.As the load increases, the local stress caused by the nonlinear increase of shear deformation in the kink region increases.Linear analysis cannot catch this phenomenon.

Mechanical Response of Off-Axis Compression Specimen
To verify the validity of the model, the results of off-axis tensile and compression tests of a set of continuous fiber reinforced composite unidirectional laminates were selected [13].The test material used was a carbon fiber/epoxy system (Toho

Finite Element Modeling
In this paper, the delamination damage is analyzed by the VCCT method based on fracture mechanics in ABAQUS software.This method is based on linear elastic fracture mechanics (LEFM) to evaluate the strain energy release rate (SERR) at the crack tip.The virtual crack closure technique is based on the crack closure integration method.The basic idea is to assume that the energy released by the crack propagation Δa is equal to the energy required to close the crack.
The specific method has been built into ABAQUS software.
When using VCCT method in ABAQUS for numerical simulation analysis, appropriate failure criteria should be selected, among which BK criterion is a criterion that is mostly used currently.The BK (Benzeggagh-Kenane) criterion is a commonly used failure criterion for judging the delamination of mixed modes.The expression is as follows: ( ) In order to facilitate the comparison with the existing experimental results, the finite element modeling in this paper is based on the testing program of double lap joint proposed in literature [2].The test sample format and specific dimensions are shown in Figure 10 below.The information of FE model is also shown in Figure 10.
The material used in the experiment is IM600/Q133, the mechanical properties of this material are already listed in Table 2 and the strength properties are

Finite Element Analysis Result
Comparing the results of the numerical analysis with the experimental results, the following Figure 11 is obtained.From the experimental data, it can be seen that the load and the displacement are in a linear relationship at the beginning of the loading, indicating that the stiffness of the laminate has not changed significantly.When the displacement reaches 0.4 mm, the response curve showed a turning point, indicating that there was more damage to some materials.Some fluctuations happen on the curve afterwards, and the fluctuation was large when the displacement reached 0.7 mm, indicating that the damage has progressed, but the overall stiffness of the slab does not decrease during this stage, indicating that the main mode of damage has not changed.When the displacement reaches around 1.0 mm, the curve turns again and the stiffness decreases again.It indicates that a new damage pattern has occurred.
In addition to the overall response, the finite element model developed in this paper can also simulate the intralaminar and interlaminar damage evolution progress as in Figure 12. 2) The model is implemented in commercial FEA software ABAQUS through UMAT subroutine interface.

Conclusions
3) The model is validated for strength prediction and mechanical response prediction of unidirectional laminate by experiment results.And the transformation matrix rotating the stress and strain from original coordinate to new coordinate is expressed as follows. 2 respectively.The stress in loading direction is x σ , and the stress state of unidirectional laminate is as fol- lows:

Figure 1 .
Figure 1.Coordinate system of material and fracture plane.
ϕ is deduced from the longitudinal compressive strength of the material.When unidirectional laminate only subject to longitudinal compressive load 11 C X σ = − , the stress state in kinking band is:

Figure 3 .
Figure 3. Schematic of fiber kinking failure analysis model.

Figure 5 .
Figure 5.The predicted shear response and the experiment result under hydrostatic pressure.

Figure 6 .
Figure 6.Predicted axial compressive strength under different transverse hydrostatic pressure and experiment result.
The analysis model developed is implemented by using the UMAT subroutine interface provided by Abaqus.The flow chart of the subroutine is shown is Figure 9.

Figure 7 .
Figure 7. Off-axis loading of unidirectional composites (x is loading direction; 1 and 2 are material axes).

Figure 8 .
Figure 8.(a) Predicted stress-strain response of IM600/Q133 in loading direction in comparison with experiment results under different off-axis loadings including compression and tension with a off-axis angle of 30˚; (b) Predicted stress-strain response of IM600/Q133 in loading direction in comparison with experiment results under different off-axis loadings including compression and tension with a off-axis angle of 90˚.

1 )
In this paper, an elasto-plastic damage model considering the non-symmetry of composite material behavior under tension and compression, failure judgement and damage evolution is developed to describe the mechanical behavior of composite laminates under both tensile and compressive load.

4 )
A progressive failure of composite laminate under bearing load is proceeded using the elasto-plastic damage model.Delamination is taken into account by a fracture mechanics method implemented using the Virtual Crack Closure Technique (VCCT) available in ABAQUS.The numerical simulation results for joint's progressive damage and mechanical response were compared with the existing experimental data, and the reliability of this model is proved.

Figure 11 .
Figure 11.Mechanical response predicted by FEA in comparison with experiment result.
(mm) Materials Sciences and Applications

m n m n l n l n l m l m m n m n l n l n l m l m m n m n l n l n l m l
)
The properties and parameters used for numerical analysis are in Table 2. K. Xue DOI: 10.4236/msa.2018.97042586 Materials Sciences and Applications