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 joints, so the strength of mechanically fastened joints in composite laminates has great influence on the strength and integrity of aircraft structures.
The failure of mechanically fastened joints in composite materials can be divided into three patterns macroscopically: shear-out, net-tension and bearing, in which bearing failure shows the best carrying capacity  . So, bearing failure is an ideal failure pattern of mechanically fastened joints in composite materials and it’s of great importance to research on it.
Xiao and Ishikawa   studied bearing strength of bolted composite joints and bearing failure behavior of composite laminates by experiment. A FEA model was also built up using commercial FE software Abaqus which considered the contact condition between bolt hole and bolt, accumulated damage model, finite deformation and non-linear behavior of composite materials. Non-linear shear theory was incorporated with continuum damage mechanics (CDM) in their study. Failure initiation was predicted by Hashin failure criteria and Yamada-Sun criteria; stiffness degradation model was employed to consider the redistribution of stress during failure process. Despite these advantages, an assumption of plane stress state was made in this research, so delamination and damage in thickness direction which played key roles in bearing failure of composite laminates could not be taken into consideration.
Gamble et al.  developed a 3D analysis model based on shell element provided by Abaqus in which modified Hill criteria was utilized to predict fiber breakage, matrix crack and delamination. A resin layer which had the same modulus and failure strength with pure resin was assumed to exist between laminas; failure of resin layer represented the initiation of delamination. The predicted failure strength had only 5% error with the experiment result in a case of [0/90]s lay-up laminate tensile test.
Camanho et al.  studied delamination in bolted composite joints by a 3D FEA model. They concluded that pressure applied by bolt and friction between washer and laminate influenced the initiation and accumulation of damage a lot. Simultaneously, the failure patterns of bolted composite joints-bearing, net tension and shear out were also analyzed with the 3D model.
Nguyen  analyzed the failure of a composite laminate with an open-hole under tensile loading by introducing a 3D FEA model into SAMCEF. It was concluded that mesh density had little impact on the prediction of damage evolution and FEA model based on CDM could predict transverse matrix cracks and damage accumulation well.
Atas et al.  used cohesive element to research delamination in pin-loaded composite laminates in detail. A bilinear traction-separation law was introduced to describe interlaminar behavior during the propagation of delamination and full integration element was employed around pin-hole to capture the stress gradient more accurately. Bolted joint in a [0/90]s lay-up composite laminate was simulated to find a non self-similar delamination propagation with an eclipse crack tip distribution. Delamination initiated at interface of 0 and 90 layers where a maximum average shear stress exists. But no progressive damage model was incorporated into this analysis which means that the effect of intralaminar damage was ignored.
Shen et al.  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  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    . 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.
2. Elasto-Plastic Damage Model
2.1. Modified Sun-Chen One Parameter Plasticity Model
Wang et al.  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 is transverse tensile yield strength, 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:
Generalized Hill yield function is taken as the plastic potential function and 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:
Using Voigt tensor notation, the non-zero terms in is expressed explicitly in Equation (11):
The determination method of the parameters in this model can be obtaining by referring to Wang et al.  . 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 is obtained by comparing the response of two off-axis compression experiment of different angles.
2.2. Failure Criteria
2.2.1. 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  is used here to predict the initiation of matrix failure as in Equation (14), when , matrix failure will initiate.
where and are the transverse shear stress and longitudinal shear stress on the potential fracture plane respectively, is the normal stress on the potential 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. and are transverse fracture resistance and longitudinal resistance on potential fracture plane respectively. and are inclination coefficient or frictional 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 , the stress components on fracture plane are notated as , and respectively. The stress in loading direction is , and the stress state of unidirectional laminate is as follows:
Firstly, a small angle (typically ) is taken to determine :
Then two other angles ( ) and (typically ) are used to determinate and :
Figure 1. Coordinate system of material and fracture plane.
2.2.2. 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     , 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  . 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 coordinate 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 (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 and shear strain , which means . is deduced from the longitudinal compressive strength of the material. When unidirectional laminate only subject to longitudinal compressive load , the stress state in kinking band is:
Figure 2. Typical fiber kinking failure mode  .
Figure 3. Schematic of fiber kinking failure analysis model.
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, is not an initial misalignment angle in reality, which means it’s not possible to get the value of by experiment. is actually an effective initial misalignment angle derived from the longitudinal compressive strength according to the material’s constitutive model. If different constitutive models are chosen to represent the mechanical behavior of material, different values of will be obtained.
After 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 exceeds one, failure will happen.
where is Macauley operator, represent .
For fiber tensile failure, maximum stress criterion is used to predict failure initiation. When the value of exceeds one, failure will happen.
2.3. Damage Evolution
If failure initiates, material stiffness is to be degraded to consider the effect of damage evolution. According to literatures  , degraded material stiffness matrix is in the form as in Equation (25):
where and 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.
3. Model Validation
3.1. 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   . 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.
3.2. 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  . The test material used was a carbon fiber/epoxy system (Toho Rayon IM600/Q133, ). The properties and parameters used for numerical analysis are in Table 2.
Table 1. Material properties and parameters for T300/PR319 and S2-glass/Epoxy.
Table 2. Material properties for IM600/Q133.
Figure 4. Influence of parameter on degradation process.
Figure 5. The predicted shear response and the experiment result under hydrostatic pressure.
Figure 6. Predicted axial compressive strength under different transverse hydrostatic pressure and experiment result.
The modified Sun-Chen plasticity model considering the asymmetry of tension and compression was used to calculate the off-axis tensile and compression tests for 30° and 90°. The stress-strain response and test results obtained in the loading direction (shown in Figure 7) were obtained. For comparison, the verification results obtained are shown in Figure 8. It is easy to see that this model can describe the nonlinear response of fiber-reinforced composites under off-axis loads and the difference in response under tensile and compressive loads.
4. Bearing Failure in Composite Laminates
4.1. UMAT Implementation
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. Off-axis loading of unidirectional composites (x is loading direction; 1 and 2 are material axes).
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˚.
Figure 9. Flow chart of UMAT subroutine.
4.2. 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:
where , , are fracture toughness for mode I, II, III respectively. is fracture toughness for mix-mode. is the mix parameter. When exceeds 1.0, crack will initiate and is equivalent strain energy release rate of nodes.
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  . 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 as follows: XT = 2700 MPa, XC = 1037 MPa, YT = 63.7 MPa, YC = 235 MPa, SC = 140 MPa, GIC = 0.33 N/mm, GIIC = 0.8 N/mm, GIIIC = 0.8 N/mm. The thickness of single lamina is 0.125 mm and the lay-up of specimen is [45/-45/90/0]2S and total thickness of specimen is 2 mm.
4.3. 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.
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.
2) The model is implemented in commercial FEA software ABAQUS through UMAT subroutine interface.
3) The model is validated for strength prediction and mechanical response prediction of unidirectional laminate by experiment results.
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 10. Flow chart of UMAT subroutine.
Figure 11. Mechanical response predicted by FEA in comparison with experiment result.
Figure 12. Intralaminar and interlaminar damage state at different stages.
5) For further study of this topic, development of a model which can reveal the complex mechanisms of interactions among different damage patterns is highly recommended by the author.
No acknowledgement is declared by the author.
Appendix 1: Stress and Strain Rotation Relationship
The stress and strain in original coordinate ( ) are notated as
The stress and strain in new coordinate ( ) are notated as
The relationship between and as well as and is as follows.
, , ,
And the transformation matrix rotating the stress and strain from original coordinate to new coordinate is expressed as follows.
, , ,
where represent the included angle cosines between the axis of original coordinate and new coordinate.