Progressive Damage Numerical Modelling and Simulation of Aircraft Composite Bolted Joints Bearing Response

Finite element numerical progressive damage modelling and simulations applied to the strength prediction of airframe bolted joints on composite laminates can lead to shorter and more efficient product cycles in terms of design, analysis and certification, while benefiting the economic manufacturing of composite structures. In the study herein, experimental bolted joint bearing tests were carried out to study the strength and failure modes of fastened composite plates under static tensile loads. The experimental results were subsequently benchmarked against various progressive damage numerical modelling simulations where the effects of different failure criteria, damage variables and subroutines were considered. Evidence was produced that indicated that both the accuracy of the simulation results and the speed of calculation were affected by the choice of user input and numerical scheme.


Introduction
The use of resin-based carbon fiber composite materials in modern large civil aircraft has been increasing significantly in recent decades. Bolted joints are still widely used on composite airframe structures due to their ease of installation and disassembly and their damage tolerance characteristics. The strength estimation in terms of bolt bearing against composite laminated plates is important for the design optimization of aircraft composite structures. Over the past 30 years, a significant amount of research has been published regarding the numerical and experimental behaviors of composite bolted joints.
Analytical models have been widely used to predict the behavior of single-lap composite bolted joints, and the effects of geometrical parameters, material properties, the stacking sequences of composites, bolt torque and the friction coefficient on the stiffness of composite bolted joints have been included in these analytical models. These models are regarded as valuable preliminary tools for the analysis of the stiffness of composite bolted joints [1][2][3][4][5].
Finite element technology has been applied in the current study to determine the structural behaviors of composite bolted joints where not only could strain and stress be calculated for the whole

Description of Specimen
The single lap joint is a common means of joining airframe components on plate structures. In the experimental study, single-lap composite bolted joints were tested that consisted of IMS-977-2 carbon fiber/epoxy matrix composite lap plates joined by HST12 Hi-Lite fasteners. The geometry and size of specimens is shown in Figure 1, following the ASTM D5961 standard [25]. The stacking sequence of the laminated plate was [45 • /90 • /-45 • /0 • /90 • /0 • /-45 • /90 • /45 • /-45 • ]s, for a total of 20 layers with single-layer nominal thicknesses of 0.188 mm. The HST12 Hi-Lite fastener with self-locking characteristics consisted of titanium alloy Ti-6Al-4V pins and stainless steel CRES347 nuts.

Description of Specimen
The single lap joint is a common means of joining airframe components on plate structures. In the experimental study, single-lap composite bolted joints were tested that consisted of IMS-977-2 carbon fiber/epoxy matrix composite lap plates joined by HST12 Hi-Lite fasteners. The geometry and size of specimens is shown in Figure 1, following the ASTM D5961 standard [25]. The stacking sequence of the laminated plate was [45°/90°/-45°/0°/90°/0°/-45°/90°/45°/-45°]s, for a total of 20 layers with single-layer nominal thicknesses of 0.188 mm. The HST12 Hi-Lite fastener with self-locking characteristics consisted of titanium alloy Ti-6Al-4V pins and stainless steel CRES347 nuts.

Experimental Process
Experimental tests were carried out at room temperature (24 ± 3 °C ) and humidity (55 ± 5%), and the specimen rested in these room/laboratory conditions for three hours. As shown in Figure 2, the tensile experiment was performed on a CMT5504 Electronic universal testing machine (MTS System Corporation, Eden Prairie, MN, USA) with a 100-kN load capability. The relative error of force and displacement indication of the machine was ±0.5%. The grip holder moving speed was adjusted by step less speed regulation, and the accuracy of moving speed could be controlled within ±0.5%. Each specimen contained two support plates bonded to composite components to minimize the eccentricity in applied force from the grip holder of the test machine, and the thickness of support plates was equal to the sum of the composite plate thicknesses.

Experimental Process
Experimental tests were carried out at room temperature (24 ± 3 • C) and humidity (55 ± 5%), and the specimen rested in these room/laboratory conditions for three hours. As shown in Figure 2, the tensile experiment was performed on a CMT5504 Electronic universal testing machine (MTS System Corporation, Eden Prairie, MN, USA) with a 100-kN load capability. The relative error of force and displacement indication of the machine was ±0.5%. The grip holder moving speed was adjusted by step less speed regulation, and the accuracy of moving speed could be controlled within ±0.5%. Each specimen contained two support plates bonded to composite components to minimize the eccentricity in applied force from the grip holder of the test machine, and the thickness of support plates was equal to the sum of the composite plate thicknesses.

Description of Specimen
The single lap joint is a common means of joining airframe components on plate structures. In the experimental study, single-lap composite bolted joints were tested that consisted of IMS-977-2 carbon fiber/epoxy matrix composite lap plates joined by HST12 Hi-Lite fasteners. The geometry and size of specimens is shown in Figure 1, following the ASTM D5961 standard [25]. The stacking sequence of the laminated plate was [45°/90°/-45°/0°/90°/0°/-45°/90°/45°/-45°]s, for a total of 20 layers with single-layer nominal thicknesses of 0.188 mm. The HST12 Hi-Lite fastener with self-locking characteristics consisted of titanium alloy Ti-6Al-4V pins and stainless steel CRES347 nuts.

Experimental Process
Experimental tests were carried out at room temperature (24 ± 3 °C ) and humidity (55 ± 5%), and the specimen rested in these room/laboratory conditions for three hours. As shown in Figure 2, the tensile experiment was performed on a CMT5504 Electronic universal testing machine (MTS System Corporation, Eden Prairie, MN, USA) with a 100-kN load capability. The relative error of force and displacement indication of the machine was ±0.5%. The grip holder moving speed was adjusted by step less speed regulation, and the accuracy of moving speed could be controlled within ±0.5%. Each specimen contained two support plates bonded to composite components to minimize the eccentricity in applied force from the grip holder of the test machine, and the thickness of support plates was equal to the sum of the composite plate thicknesses.  After the machine was preheated for more than 15 min, the grip holder was pulled at a speed of 2 mm/min to simulate quasi-static loading, then stopped after the load dropped by approximately 30% from the maximum value. The numerical values of applied load and grip holder displacement were recorded automatically. The test procedure corresponds to the guidelines given in ASTM D5961 [25].

Material Properties
The basic mechanical property parameters of the IMS-977-2 composite lamina are given in Table 1 [26,27]. The stress-strain curves of titanium alloy Ti-6Al-4V and stainless steel CRES347 are shown in Figure 3 [28], and the corresponding relationship between stress and strain ( Table 2) can be obtained according to the stress-strain curves.

Progressive Damage Model
Many failure criteria have been proposed to predict fiber reinforced composite material failure. The early failure criteria did not consider the actual failure mechanisms, whereas later phenomenological failure criteria considered how failure would occur as long as the stress reached its ultimate strength (i.e., maximum-stress criteria), as shown in Equations (1)-(4). Equations (5)-(8) depict four Hashin criteria failure modes that are considered to be the earliest three-dimensional stress states and failure mechanisms. Micromechanical behavior was considered according to Puck criteria, where the degree of fiber failure is not only related to the stress state of the composite but also to the volume fraction of the fiber and matrices, as shown in Equations (9)- (15). For matrix failure, Puck thought that the matrix fracture occurred on the plane parallel to the fiber direction, while the fracture angle was used to describe the deflection of the fracture plane. The range of the fracture angle was from 0 • to 180 • , meeting the requirement of the maximum damage factor. In LaRc05 criteria, as shown in Equations (16)-(26), Pinho introduced the effects of in-situ strength to predict matrix criteria based on Puck criteria, and as a result fiber kinking prediction was built based on plasticity theory, which is different from the current damage criteria.
To search for the fracture angle of a matrix, Puck proposed the solution of a stepwise calculation of the maximum damage factor using an interval of 1 • within its limits. For each stress state, 180 calculations are required, which is the same problem as for LaRc05 criteria when searching for the fiber misalignment angle. This may be a computationally expensive calculation, and severe convergence problems may occur in implicit finite element. For these reasons, a certain method is proposed in this paper based on the derivatives of continuous functions. As shown in Figure 4, several angles were selected that could derive a damage factor equal to zero, and the damage factors corresponding to these angles were calculated. The maximum damage factor and required matrix fracture angle or fiber misalignment angle were selected more quickly than the conventional method.
Materials 2020, 13, x FOR PEER REVIEW 6 of 16 while the fracture angle was used to describe the deflection of the fracture plane. The range of the fracture angle was from 0° to 180°, meeting the requirement of the maximum damage factor. In LaRc05 criteria, as shown in Equations (16)-(26), Pinho introduced the effects of in-situ strength to predict matrix criteria based on Puck criteria, and as a result fiber kinking prediction was built based on plasticity theory, which is different from the current damage criteria.
To search for the fracture angle of a matrix, Puck proposed the solution of a stepwise calculation of the maximum damage factor using an interval of 1° within its limits. For each stress state, 180 calculations are required, which is the same problem as for LaRc05 criteria when searching for the fiber misalignment angle. This may be a computationally expensive calculation, and severe convergence problems may occur in implicit finite element. For these reasons, a certain method is proposed in this paper based on the derivatives of continuous functions. As shown in Figure 4, several angles were selected that could derive a damage factor equal to zero, and the damage factors corresponding to these angles were calculated. The maximum damage factor and required matrix fracture angle or fiber misalignment angle were selected more quickly than the conventional method.  is the normal stress on the potential fracture plane; α is fracture angle; and are transverse fracture resistance and longitudinal resistance on the potential fracture plane, respectively; and are the inclination coefficient or frictional coefficient, respectively, which represent the influence of normal stress on the fracture resistance [22].
The angle of the kink band, ѱ, is found numerically in the range 0° and 180° so as to maximize τ 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 potential fracture plane; α is fracture angle; S T and S L are transverse fracture resistance and longitudinal resistance on the potential fracture plane, respectively; µ T and µ L are the inclination coefficient or frictional coefficient, respectively, which represent the influence of normal stress on the fracture resistance [22].
The angle of the kink band, ψ, is found numerically in the range 0 • and 180 • so as to maximize the failure index in Equation (17). The misalignment angle, ϕ, is the sum of the initial misalignment angle ϕ0 (manufacturing defect) and the shear strain expressed in a coordinate system aligned with the manufacturing defect. This is calculated based on the linear shear response assumption. The McCauley brackets are defined as <x>+ = max {0,x}.
The above failure criteria in Table 3 were combined with the exponential damage variable, as shown in Equation (27) and sudden drop damage variable as shown in Equations (28) and (29) together, and implemented in the ABAQUS user-defined material subroutine UMAT to study the effect of progressive damage on the numerical simulation results.
Hashin criteria Fiber Tensile Failure Fiber Compressive Failure Matrix Tensile Failure Matrix Compressive Failure Puck criteria Fiber Compressive Failure Matrix Tensile Failure Matrix Compressive Failure Matrix Tensile Failure The exponential damage variable is written as where F is the value of the failure criterion, T is the coefficient in the stiffness matrix (e.g., C 11 , C 22 , C 33 ), ε t is ultimate failure strain, L is element characteristic length, and G is fracture toughness [29]. The sudden drop damage variable is written as

The Effect of Failure Criteria
The load-displacement curves obtained from the experiment and finite element method are presented in Figure 5 as different failure criteria but using the same exponential damage variables. The ultimate failure loads calculated by the finite element method were close to the experimental results. The maximum difference between the numerical simulation and the average value (17612N) of the experiment test was 4.8%, as calculated by the maximum-stress and Hashin criteria. The minimum difference was 1.1%, as calculated by the LaRC05 criteria. On the other hand, the difference in joint stiffness in the linear elastic stage was relatively larger (the maximum and minimum difference were 29.4% and 17.7%, respectively), and the curves obtained by the finite element method did not have a marked decline around the end points when compared with the experiment results. Possible reasons for these differences are as follows: (1) to achieve a good convergence when simulating the failure of a composite in ABAQUS, the viscous regularization of damage variables is often introduced to delay the degradation of the composite; (2) the material properties of fasteners are defined as elastic-plastic, while the protruding heads of fasteners fractured during the experimental tests, as shown in Figure 6.

Fracture of screw head
Tensile direction Figure 6. Experiment specimens.
The microscope images of damaged composite plates in the experimental tests and the prediction of composite failure in the finite element are presented in Table 4. The bolt hole of the composite plate was squeezed seriously by the bolt shank, causing several damage modes to occur such as matrix crushing, fibers pulling out and fiber fracture. It can be seen that while the maximum stress, Hashin and Puck criteria presented the same accuracy of prediction for fiber compression, the LaRc05 criteria could predict the boundary of the damaged region well. For the matrix failure prediction, Larc05 and Puck criteria were able to predict the damaged area of the matrix accurately; however, the damaged area predicted by Hashin criteria was obviously larger than the experimental result, while the damaged area was much smaller when predicted using maximum-stress criteria.

Fracture of screw head
Tensile direction Figure 6. Experiment specimens.
The microscope images of damaged composite plates in the experimental tests and the prediction of composite failure in the finite element are presented in Table 4. The bolt hole of the composite plate was squeezed seriously by the bolt shank, causing several damage modes to occur such as matrix crushing, fibers pulling out and fiber fracture. It can be seen that while the maximum stress, Hashin and Puck criteria presented the same accuracy of prediction for fiber compression, the LaRc05 criteria could predict the boundary of the damaged region well. For the matrix failure prediction, Larc05 and Puck criteria were able to predict the damaged area of the matrix accurately; however, the damaged area predicted by Hashin criteria was obviously larger than the experimental result, while the damaged area was much smaller when predicted using maximum-stress criteria. The microscope images of damaged composite plates in the experimental tests and the prediction of composite failure in the finite element are presented in Table 4. The bolt hole of the composite plate was squeezed seriously by the bolt shank, causing several damage modes to occur such as matrix crushing, fibers pulling out and fiber fracture. It can be seen that while the maximum stress, Hashin and Puck criteria presented the same accuracy of prediction for fiber compression, the LaRc05 criteria could predict the boundary of the damaged region well. For the matrix failure prediction, Larc05 and Puck criteria were able to predict the damaged area of the matrix accurately; however, the damaged area predicted by Hashin criteria was obviously larger than the experimental result, while the damaged area was much smaller when predicted using maximum-stress criteria.

The Effect of Damage Variables
When the same failure criteria (LaRC05) were introduced with different damage variables in UMAT, the numerical simulation results showed an obvious difference. Compared to the use of exponential damage variables, sudden decrease damage variables reduced the bearing strength by 2.1%, as shown in Figure 7. When the sudden decrease damage variables were used, the loaddisplacement curve began to show a significant downward trend after the ultimate failure load was reached, which is close to the experimental results. Unlike sudden decrease damage variables, which are directly equal to 1 when the value of the failure criteria expression is >1, exponential damage variables increase from zero to 1 gradually with increases in the value of the criterion expression. The matrix stiffness of the composite presented faster degradation when the sudden decrease damage variable was used, and the damage propagation on the composite plates was faster as well.

The Effect of Damage Variables
When the same failure criteria (LaRC05) were introduced with different damage variables in UMAT, the numerical simulation results showed an obvious difference. Compared to the use of exponential damage variables, sudden decrease damage variables reduced the bearing strength by 2.1%, as shown in Figure 7. When the sudden decrease damage variables were used, the load-displacement curve began to show a significant downward trend after the ultimate failure load was reached, which is close to the experimental results. Unlike sudden decrease damage variables, which are directly equal to 1 when the value of the failure criteria expression is >1, exponential damage variables increase from zero to 1 gradually with increases in the value of the criterion expression. The matrix stiffness of the composite presented faster degradation when the sudden decrease damage variable was used, and the damage propagation on the composite plates was faster as well.
Moreover, there was plastic deformation on the bolt head, screw and nut when the sudden decrease damage variable was used in UMAT, which is closer to the experimental phenomenon shown in Figure 8b,c. When the exponential damage variable was used, the plastic deformation of the fastener occurred only in the middle of the screw, as shown in Figure 8a.  Moreover, there was plastic deformation on the bolt head, screw and nut when the sudden decrease damage variable was used in UMAT, which is closer to the experimental phenomenon shown in Figure 8b,c. When the exponential damage variable was used, the plastic deformation of the fastener occurred only in the middle of the screw, as shown in Figure 8a.

The Effect of Subroutine
As the same failure criteria, damage variables and damage stiffness matrix were used in subroutine USDFLD and UMAT, respectively, the difference in the ultimate failure load was 14.2%, as shown in Figure 9. There were severe non-convergence in the USDFLD calculation process when   Moreover, there was plastic deformation on the bolt head, screw and nut when the sudden decrease damage variable was used in UMAT, which is closer to the experimental phenomenon shown in Figure 8b,c. When the exponential damage variable was used, the plastic deformation of the fastener occurred only in the middle of the screw, as shown in Figure 8a.

The Effect of Subroutine
As the same failure criteria, damage variables and damage stiffness matrix were used in subroutine USDFLD and UMAT, respectively, the difference in the ultimate failure load was 14.2%, as shown in Figure 9. There were severe non-convergence in the USDFLD calculation process when the material property of the bolt was defined as elastic-plastic, but this was avoided in the UMAT calculation process. Unlike subroutine USDFLD, which can call strain and stress directly, subroutine

The Effect of Subroutine
As the same failure criteria, damage variables and damage stiffness matrix were used in subroutine USDFLD and UMAT, respectively, the difference in the ultimate failure load was 14.2%, as shown in Figure 9. There were severe non-convergence in the USDFLD calculation process when the material property of the bolt was defined as elastic-plastic, but this was avoided in the UMAT calculation process. Unlike subroutine USDFLD, which can call strain and stress directly, subroutine UMAT can only call strains directly, then calculate stress according to the damage stiffness matrix and strain of the composite. The damage variable is the regularized viscosity before the calculation of stress, making the damage variable smaller than the real value and always less than 1, as shown in Equation (30). In other words, the degradation of composite stiffness was delayed by the viscous regularization of the damage variable, improving the convergence of the calculation.
where d V is damage variable of the current incremental step with regularized viscosity, η is the viscosity coefficient, ∆t is the time increment, and d is the damage variable of the previous incremental step with regularized viscosity.
where is damage variable of the current incremental step with regularized viscosity, is the viscosity coefficient, Δ is the time increment, and ′ is the damage variable of the previous incremental step with regularized viscosity.

Conclusions
In order to improve the accuracy and efficiency of simulating the bearing behaviors of composite bolted joints, the effects of failure criteria, damage variables and user subroutines were studied in this paper, and experimental tests were carried out to compare numerical results. An approach based on a derivative method was used to find the fiber misalignment angle and the matrix fracture angle by applying LaRc05 criteria, and it was found this method had better efficiency than the conventional ergodic method.
When combined with the same damage variable, the maximum stress, Hashin and Puck criteria all presented the same accuracy at predicting fiber compression, and the LaRc05 criteria were able to predict the boundary of the damaged region well. LaRc05 and Puck criteria presented more accurate results than maximum-stress and Hashin criteria at predicting the matrix failure.
Compared to exponential damage variables, using the sudden decrease damage variables in UMAT could more accurately replicate experimental results when combined with LaRc05 criteria. As the same criteria and damage variables were incorporated in different user subroutines, there was better convergence in UMAT calculation than USDFLD calculation.

Conclusions
In order to improve the accuracy and efficiency of simulating the bearing behaviors of composite bolted joints, the effects of failure criteria, damage variables and user subroutines were studied in this paper, and experimental tests were carried out to compare numerical results. An approach based on a derivative method was used to find the fiber misalignment angle and the matrix fracture angle by applying LaRc05 criteria, and it was found this method had better efficiency than the conventional ergodic method.
When combined with the same damage variable, the maximum stress, Hashin and Puck criteria all presented the same accuracy at predicting fiber compression, and the LaRc05 criteria were able to predict the boundary of the damaged region well. LaRc05 and Puck criteria presented more accurate results than maximum-stress and Hashin criteria at predicting the matrix failure.
Compared to exponential damage variables, using the sudden decrease damage variables in UMAT could more accurately replicate experimental results when combined with LaRc05 criteria. As the same criteria and damage variables were incorporated in different user subroutines, there was better convergence in UMAT calculation than USDFLD calculation.