Progressive Failure Analysis in Open-Hole Tensile Composite Laminates of Airplane Stringers Based on Tests and Simulations

: The failure types and ultimate loads for eight carbon-epoxy laminate specimens with a central circular hole subjected to tensile load were tested experimentally and simulated using two different progressive failure analysis (PFA) methodologies. The ﬁrst model used a lamina level modeling based on the Hashin criterion and the Camanho stiffness degradation theory to predict the damage of the ﬁber and matrix. The second model implemented a micromechanical analysis technique coined the generalized method of cells (GMC), where the 3D Tsai–Hill failure criterion was used to govern matrix failure, and the ﬁber failure was dictated by the maximum stress criterion. The progressive failure methodology was implemented using the UMAT subroutine within the ABAQUS/implicit solver. Results of load versus displacement and failure types from the two different models were compared against experimental data for the open hole laminates subjected to tensile displacement load. The results obtained from the numerical simulation and experiments showed good agreement. Failure paths and accurate damage contours for the tested specimens were also predicted.


Introduction
Utilizing composite laminates in aircraft structures plays quite an important role owing to their low weight, high strength, the capability to design varieties of geometries and other advantages. The ultimate load and failure type of composite structures are significant parameters to satisfy the engineering requirement, especially those discontinuous structural components such as laminates with notches or open holes that may cause complex stress states and damages [1][2][3][4].
Failure analysis of composite panels containing notches or holes is regarded as an important subject in structural design [5,6]. An open hole existing in a laminate can lead to stress concentration that produces complicated damage processes and failure mechanisms [7,8]. The breakage occurring around the open hole induces the redistribution of the stresses in the composite laminates that contributes to causing the catastrophic failure [9]. Matrix cracking, fiber breakage, delamination, etc., observed in an open hole laminate subjected to tensile loads are common failure mechanisms [10][11][12]. Many analytical and numerical models have been used to predict the failure strength of composite laminates with open holes [13,14]. In most cases, linear elasticity and a failure criterion are utilized to forecast the onset of failure in composite structures [15,16]; some of these models are examined in the world wide failure exercise (WWFE) [17][18][19]. Although many criteria distinguish the failure of fiber and matrix and take progressive failure modeling into account [20][21][22], they are not available to capture the interaction between the constituents at the microscale level.
Employing micromechanics to acquire the damage states of composite structures is of great significance [23,24]. Variations on constituent material properties along with the microstructural characteristics, for instance, the fiber volume fraction, orientation, and microlevel architecture can be taken into account by using microscopic mechanical methods to better research the failure mechanism of composite structures. In addition, the micromechanical method can also be easy to estimate the interactive effects among the different constituents in the composites.
The generalized method of cells (GMC), initiated by Paley and Aboudi [25,26] and then developed for improving the computational efficiency [27], is actually analytical, which has been proved to be an extremely effective micromechanics method. Its formulation contains the application of a few governing equations using a homogenization method, by which traction and displacement continuity equations are imposed at the whole subcell interfaces and the periodic boundaries of the repeating unit cell (RUC) in an average sense. A rectangular RUC divided into desired quantities of rectangular subcells is used to account for composites properties, and each subcell is available to occupy a specified material. These continuity conditions directly contribute to formulating a strain concentration matrix that is related to the local subcell strains to the global average applied strains. Given the local constitutive law, the subcell stresses can be represented by the global average applied strains. Ultimately, the entire RUC stiffness is obtained by employing a homogenization strategy. It is robust enough to implement a multiscale framework using the micromechanics model GMC together with the FEM, wherein the integration points of the finite element model call the GMC as the material point response, and damage in the constituents in the composite laminates are captured locally via the GMC [28].
The main objective of the current work is to use the GMC to model damage and failure mechanisms in the composite laminates with an open hole under tensile loading, and the results obtained are compared to those obtained from the tests and the popular Hashin failure criteria. The failure loads and the failure types of angle-ply laminates with a centrally located open hole were studied through the tensile tests and the finite element method. Numerical analyses were modeled in two ways. The first model utilized the 3D Tsai-Hill criterion for the matrix constituent and maximum stress criteria for the fiber constituent by using the multiscale analysis method of GMC. The second model employed the Hashin criterion in conjunction with the Camanho stiffness degradation method to perform the macrolevel analysis [29,30]. The results obtained were arranged in tabular and graphical form, which exhibit damage modes and their corresponding ultimate loads for the composite panels. Numerical results, such as ultimate loads and failure modes, were compared with the experimental data. The analysis model with failure criteria for PFA was accomplished using a user-defined material subroutine (UMAT) in the ABAQUS finite element analysis code, and its implicit solver was utilized in the simulation models.

Constitutive Models
It is deemed that fiber-reinforced composites can be justified as a collection of RUCs which are divided into N β × N γ (β, γ = 1, 2, . . . , N β , N γ ) subcells as displayed in Figure 1. The number of subcells is arbitrary, and each of them can contain a specified material representing the constituents in the composite; the constituents are matrices and fibers in general, and the interfaces between the subcells can be ignored in the perfect bonding conditions [26].
3 ) is set for each subcell of the RUCs, where the x 1 axis shown in Figure 1 represents the fiber direction, and the size of each subcell is defined as h β and l γ along the direction of the x 2 and the x 3 . It is assumed employing a linear displacement approximation is suitable to obtain the average behavior of the RUCs, where ω ij , the local strain tensors can be represented as follows, where Herein, only elastic subcell strains, ε ij (βγ), will be taken into account; Hooke's law is utilized to express subcell stresses in terms of its corresponding strains is written as, where C (βγ) ij denotes the elastic stiffness matrix of the subcell βγ, and D max represents the stiffness deduction coefficient of each failed subcell.

Continuity Conditions
The interfaces between adjoining subcells along with the boundaries between neighboring RUCs are required to be imposed the displacement and traction continuity conditions in an average sense. Hill's strain concentration matrix is constructed using these continuity conditions to make the unknown average subcell strains be related to the applied global average strains. Hence, the subcell strains are utilized to provide the global average strains of RUCs by defining the equation as follows: where h, l denote the sizes of RUC, as well as h β and l γ represent the length of subcell βγ, respectively. Next, imposing interfacial displacement continuity between adjacent subcells produces the following equations: Equations of (6) are available to be represented in the matrix form, where ε = [ε 11 , ε 22 , ε 33 , 2ε 23 , 2ε 13 , 2ε 12 ], ε s = ε (11) , ε (12) , . . . , ε (N β N γ ) , J = [1, h, l, h, l, hl] (8) and where ε (βγ) are components of the average subcell strains, and the matrix A G contains geometrical parameters of each subcell. In order to determine the unknown subcell strains successfully, the interfacial traction continuity is necessary to be imposed to ensure enough equations, which can be written as follows, whereβ andγ are given by These traction conditions can be reformulated in matrix form in terms of the average subcell strains utilizing the subcell constitutive relationships, where the matrix A M possesses the material properties in each subcell.

Determining Subcell Strains
After A G , A M and J have been determined, the subcell strains can be expressed in compact matrix form, where After the subcell strains are acquired, it is convenient to obtain the subcell stresses employing the constitutive law.

Finite Element Models
The open-hole laminated composite plates containing 24 plies were obtained from a new stringer that was cut to a specified length, and the width and length of the specimens were 24 mm and 250 mm, respectively. An open hole with a diameter of 4.76 mm was centrally located at the laminates. The number of these laminate coupons was 8, and each specimen is numbered from 1 to 8. The composite laminate material used was T800S/3900-2B, which was widely applied to the load-bearing structure of both military and civil aircraft, and the adhesion between the fiber and matrix was good. The lamina properties, such as the elastic properties and the strength parameters, as shown in Tables 1 and 2, were provided by vendors. Each ply thickness was 0.19 mm, and the number of plies, fiber orientation and total thicknesses of the laminates are shown in Table 3. The matrix material was considered to be isotropic, having Young's modulus E m and a Poisson's ratio ν m . Transversely isotropic material was applied to the fiber constituent, where five independent elastic parameters existed. E f11 represented the longitudinal modulus, E f22 denoted the transverse modulus, ν f12 was the longitudinal Poisson's ratio, and G f12 was the longitudinal shear modulus. The elastic properties of fiber and matrix, as shown in Table 4, were obtained from vendor data or were partially backed out from elastic lamina properties measured for T800S/3900-2B. A fiber volume fraction of 50.4 percent was used for the current simulation models. The strengths of the composite constituents are shown in Table 5, where X ft and X fc represented the tensile and compressive strength of the fibers, respectively, and Y mt , Y mc and T m denoted the tensile, compressive and shear strength of the matrixes, respectively. Three-dimensional eight-node reduced integration elements (C3D8R) were used in the region around the open hole of the FEA models, and the rest region utilized C3D8 in order to guarantee the accuracy of solutions, though it was also available to choose a 2D shell element for the analysis of laminate specimens. The finite element model established is shown in Figure 2. The meshes around the open hole regions were refined in size to precisely capture the stress fields and damage states. The X-Y-Z global coordinate system was introduced with its origin located in the center of the coupons, X-axis was aligned with the direction of the 0 • fiber, and Y-axis was the direction perpendicular to the 0 • fiber so Z represented the direction of the ply thickness. The left edge of the simulation model was fixed so that the whole translational degrees in each node were restricted, and the horizontal edges of the model had no constraints enforced. While the right edge was only allowed to move in the X direction, and it was imposed a displacement load. Static analysis was implemented using ABAQUS/Standard, in which the coupled damage and failure models were implemented using the user subroutine UMAT, and the microlevel GMC along with the failure criteria was also included. Once the failure criterion was met, the elastic performance of each constituent material would be reduced to 0.1% of the virgin value.

Microlevel Constituent Failure Criteria
Microlevel failure based on GMC was implemented using proper failure criterion within the constituents of the composite. Maximum stress criteria were utilized to detect the failure of fiber, which was deemed to be brittle before it broke.
where β f and γ f denoted the fiber subcell indices, and X ft and X fc represented the tensile and compressive fiber strengths, finally σ 11 was the applied axial stress. A 3D Tsai-Hill criterion was chosen to predict the failure arising in the matrix subcells, which was shown as below [31].
where Y mt , Y mc, and T m belonging to the properties of the matrix subcells represented the transverse tensile strengths, transverse compressive and shear strength, respectively. Once d m was greater than or equal to 1 in a matrix subcell or d f was less than or equal to in a fiber subcell, failure occurred, and then all the elastic properties of those failed subcells were degraded appropriately. All the stiffness matrices of the constituents were reduced by multiplying 1-D max , and D max was set to 99.9% herein.
At each material point, a micromechanical analysis was performed using GMC by calling the UMAT subroutine, in which maximum stress criteria and Tsai-Hill failure criteria were implemented to predict the failure of subcells. The square-packed RUC consisting of 3 matrix subcells colored green and 1 fiber subcells colored blue was chosen to analyze the properties of composite materials, though arbitrary architecture defined by the users was available. The number of subcells embedded in the RUCs is shown in Figure 3. The flow chart showing the hierarchy of this multiscale algorithm is shown in Figure 4.
ABAQUS/standard provided the total strain tensor fields at every element integration point of the FEM model, and the global strains were imposed as basic known variables on the RUC to calculate the global stress at each iteration based on the GMC. Additionally, the criteria applied on the constituent material reduced and updated the global stiffness calculated with the GMC, by which the multiscale method of analysis was implemented.

Lamina Level Hashin Criteria
Lamina level Hashin failure criteria [29] were also used to predict the mechanical response of the test coupons to compare with micro-level Tsai-Hill in combination with maximum stress criteria based on GMC. Hashin criteria utilized take four different failure modes into account, which are matrix cracking, matrix crushing, fiber buckling and fiber fracture. The expressions of lamina level Hashin criteria are represented, as shown below.

4.
Matrix crushing under transverse compression and shearing (σ 22 < 0): where d ft , d fc , d mt and d mc are damage variables corresponding to each damage mode.
The value of each damage variable is 0 before damage occurs. X T , X C are, respectively, the tensile and compression strength in the fiber direction of each ply, and Y T and Y C are, respectively, the transverse tensile and transverse compression strength, which are perpendicular to the fiber direction. The last parameter S ij, is the shear strength. Once any equation above is met, the damage at that integration point occurs. When the stress states of any element satisfy multiple Equations in (20)-(23) at the same time, multiple damage modes corresponding to that element are considered to have occurred, and thus the material properties embedded in the damaged elements will degrade accordingly. Herein the local degradation model, Camanho [30], is used to describe the stiffness reduction after the element is damaged, as shown in Table 6.  Where Q denotes elastic material parameters required to reduce according to the different failure modes. Specifically, Q represents the elastic modulus E 11 to be reduced when fiber rupture or fiber buckling occurs, and in the case that matrix-cracking or matrixcrushing happens, Q is set to E 11 , G 12 , G 23 . Q d indicates the corresponding elastic properties that will be degraded. If more than one failure mode has occurred simultaneously at an integration point, the stiffness of the material would be degraded cumulatively in terms of the degradation scheme.

Results and Discussion
In this chapter, the ultimate load values, load versus displacement curves and failure modes were obtained by testing and simulating the laminate coupons, which contained only ±45 • and 0 • layers, as displayed in Table 3 above. The FEA was accomplished for the laminates using both models (lamina level Hashin and micro-level GMC/Tsai-Hill), simulation results acquired were compared against experimental data, and the damage states and the final failure mode of each laminate are also described.

Experiment
All eight specimens were tested according to the ASTM D5766 standard (standard test method for open-hole tensile strength of polymer matrix composite laminates), and the MTS electrohydraulic servo loaded system having a maximum load of 250 kN with a precision of 0.5% was used in the experiment. One end of the specimen was fixed, and the other end was gradually imposed with a displacement load until the specimen tested was completely destroyed, and then the test was stopped immediately.
The experimental ultimate loads of these laminates were shown in Table 7. The maximum and minimum loads were 63.16 and 58.09, respectively. The average ultimate load for the whole 8 plates was 60.74 kN with a discrete coefficient of 2.126%, and the specimen "1-5" had the closest ultimate load of 60.99 kN to the average values. Test results showed that the load versus displacement curve of the specimen was mainly linear before reaching the ultimate load, which may be due to the fact that 0 • fiber lamination was the main load-bearing object under tensile load, and the fiber had typical brittleness characteristics before the failures occurred. Stress concentrations were present in the specimens due to the existence of the open hole, around which microcracks, such as matrix cracking, may appear. The crack expanded as the displacement load continued to increase; after the specimen reached the ultimate load, the crack rapidly propagated until the specimen was completely destroyed, and the tests were stopped. During this process, the failure modes of the specimen were mainly as follows: tensile fracture of fiber at 0 • and tensile fracture of the matrix at ±45 • , as shown in Figure 5.

Simulations
Two methods mentioned above were used to carry out the finite element analysis, and the load versus displacement curve of specimen numbering "1-5" was represented as shown in Figure 6 due to its closest failure load compared with the average experimental ultimate load. It could be seen the numerical results matched well with the experimental data, as shown in Figure 7. The ultimate load obtained using the multiscale method GMC/Tsai-Hill was 60.80 kN, which was slightly higher than 60.03 kN acquired by using the lamina-level Hashin method. The reaction force versus displacement response obtained using both methods was also desirable. The final failure load of the specimen was determined by the 0 • layer fiber, and the difference between the strength values obtained by the two analysis modes was quite small. On the whole, the simulation results were in good agreement with the experimental results.  It was observed that the failure displacement of the simulation model was smaller than that of the experimental data, which may be the reason that the effect of nonlinearity and plasticity was not taken into account in the analysis model [32]. In addition, the absence of the interface phase of composite material embedded in the GMC may be another factor contributing to the difference [33], which would be considered in our next work.
Because the specimens tested were subjected to tensile load, tensile failure of fiber in 0 • layers may be the main failure pattern, while the ±45 • layers mainly produced matrix tensile failure. Therefore, the fiber tensile failure mode of the 0 • layer and the matrix tensile failure mode of the ±45 • layer were of significant concern. The failure modes and damage processes between −45 • and 45 • for the damage contours plots were roughly the same, so the damage states of only 45 • layers were represented here. It should be emphasized that a fiber subcell failure utilizing GMC/Tsai-Hill method represented the fiber failure for the whole RUC, and the matrix subcell whose failure was dominant in the three matrix subcells would be displayed for comparison and illustration. The simulation results showed that the failure of the third subcell was dominant, so in this paper, the failure of the No. 3 matrix subcell represented the failure of the matrix constituent. Figures 8 and 9 showed contour plots of matrix damage state in 45 • layers of laminate analyzed at different percentages of elements or subcells failure using both methods. The matrix damage measured around the open hole showed a progressive outward growth pattern for these layups, with damage initiating at 50 kN of the load. Simulations using both methods generated similar matrix failure patterns, which progressed in the same manner.    Figure 10a, and continued to expand as the tensile load gradually increased. Figure 10b displayed the damage situation at 90 percent failure under the tensile load of 59 kN, after which laminated plates rapidly reached the ultimate load, and then fiber damage passed through the specimen quickly. The two methods produced different fiber failure patterns, as shown in Figure 11. A gradient of material properties due to the matrix or fiber failure could result in a bifurcated state around the open hole, but it not occurred in the experiment and the simulation by GMC/Tsai-Hill model. A similar situation existed in some of the literature [3,34]. Hence the GMC/Tsai-Hill method could produce reasonable mechanism and physical failure, and the fiber failure path produced by GMC/Tsai-Hill was considered to be more reasonable.

Conclusions
Failure loads and the failure types of angle-ply laminates with a centrally located open hole were investigated through the tensile tests and the finite element method. All eight specimens were tested, and the average ultimate load was calculated. The specimen 1-5 that had the closest ultimate load to the average value was chosen, and then its load versus displacement date and failure modes were used to compare with results obtained from two different simulation models. Tensile failure of fiber in 0 • layers was the main failure pattern, while the ±45 • layers mainly produced matrix tensile failure mode.
The analysis results of the two PFA models were all in good agreement with the test results, and the error was controlled within 7%. On the whole, the GMC/Tsai-Hill multiscale analysis method based on the microlevel failure prediction of constituent material was more accurate than the macrolevel Hashin methodology. The failure modes were basically the same between experimental data and simulation results. The load versus displacement curve was basically linear before the ultimate load appeared, and then the curve fell steeply, which was in good agreement with the test data. Both analysis models predicted good results, including failure load, failure modes, failure position and load versus displacement curves.
The matrix damage propagation predicted by the two analysis methods was similar, but the difference was present in the propagation of fiber damage, bifurcation state occurred in the macro Hashin analysis, which not appeared in the multiscale GMC analysis method. The comparison with the experimental results showed that the GMC/Tsai-Hill method could produce a desirable mechanism, physical failure and failure path.