Multi-Scale Nonlinear Progressive Damage and Failure Analysis for Open-Hole Composite Laminates

: In order to study the nonlinear behaviors and interactions among the constituents for the composite material structure under the tensile load, multiscale damage model using generalized method of cells (GMC) and a lamina-level progressive damage model were established, respectively, for ﬁber reinforced composite laminates with a central hole, which were based on the thermodynamic Schapery Theory (ST) at either the micro-level or the lamina level. Once the nonlinear progressive degradation of the matrix material reached the lower limit value for the ST method, matrix failures naturally occurred, the failure of the ﬁber was determined by the maximum stress failure criterion. For the multiscale progressive damage model, the GMC model consisting of a ﬁber subcell and three matrix subcells was imposed at each integral point of FEM elements, and the three matrix subcells undergo independent damage evolution. The load versus displacement curves and failure modes of the open-hole laminates were predicted by using the two progressive failure models, and the results were compared with that obtained by the Hashin-Rotem progressive failure model and the experimental results. The results show that the ST based method can obtain the nonlinear progressive damage evolution states and failure states of the composite at both the lamina level and the multiscale level. Finally, the damage contours and failure paths obtained are also presented.


Introduction
Polymer matrix composites (PMCs) have been widely used in aircraft structures owing to their low weight, high strength, high resistance to fatigue, and many other superior advantages. Progressive damage analysis of composite laminates is regarded as an important and complex subject, which is highlighted by the World-Wide Failure Exercise (WWFE) [1][2][3]. It is mainly based on macroscopic failure analysis methods within the WWFE to predict the final failure behavior of fiber reinforced composites. The associated parameters in these theories generally rely on extensive mechanical experiments, which result in myriad costs of time and expense. Although many criteria take different failure modes into account and incorporate progressive failure modeling [4][5][6], they are essentially phenomenological methods that cannot capture the interaction between the constituents at the microscale level. In order to more accurately characterize the mechanical response of the fiber reinforced composite structures, damage and failure mechanisms must be treated separately [7][8][9]. The gradual expansion of matrix micro-damage or micro-crack originates from the distributed micro-pores or other micro-defects produced in the manufacturing of composite structures. The damage expansion leads to stiffness reduction in local areas of the structure, resulting in redistribution of stress and strain fields, which is the main cause of nonlinear phenomena before the catastrophic failure of fiber reinforced composite materials occurs [10].

ST Model
The Schapery Theory can be used for the progressive damage analysis of the fiber reinforced composites, in which the total applied potential W T of the material is divided into the recoverable elastic part W and the dissipative part W S that can cause structural deformation [25].
A portion of the total applied work potential causes the microcracking and any other structural changes when the material is loaded, which affects the elastic properties of the material, and at the same time, the other portion of the applied work potential can be recovered when the material is unloaded. This process is shown in Figure 1, where the shaded area represents the dissipated potential W S and the area below the straight line represents the elastic potential W.
According to the ST method [22,23], both W and Ws are functions of a group of internal state variables (ISVs) S m (m = 1, 2, . . . , M). The total applied work potential is a constant with respect to each of the internal state variables [23], According to the ST method [22,23], both W and Ws are functions of a group of internal state variables (ISVs) Sm (m = 1, 2, …, M). The total applied work potential is a constant with respect to each of the internal state variables [23], The dissipated potential WS is allowed to be a function of any number of internal state variables, and different internal state variables can be used to explain different damage and failure mechanisms. In order to describe the nonlinear progressive damage of the matrix material of the fiber reinforced composites, Ws is assumed to be a function of a single internal variable S [25]. For convenience, Ws = S can be selected, Calculating the derivative of Equation (3) with respect to S yields,

Constitutive Relation
The constitutive relation of a unidirectional lamina under the plane stress state can be written as, 11 11 11 12 22 where γ12 is the engineering shear strain, and The dissipated potential W S is allowed to be a function of any number of internal state variables, and different internal state variables can be used to explain different damage and failure mechanisms. In order to describe the nonlinear progressive damage of the matrix material of the fiber reinforced composites, Ws is assumed to be a function of a single internal variable S [25]. For convenience, Ws = S can be selected, Calculating the derivative of Equation (3) with respect to S yields,

Constitutive Relation
The constitutive relation of a unidirectional lamina under the plane stress state can be written as, where γ 12 is the engineering shear strain, and where E 11 , E 22 , υ 12 , υ 21, and G 12 denote axial elastic modulus, transverse elastic modulus, Poisson's ratio, transverse Poisson's ratio, and shear modulus respectively. It can be assumed that υ 12 υ 21 1, then Equation (6) becomes Aerospace 2022, 9, 59 4 of 15

Calculating the Damage State
The damage mechanism associated with the internal state variable S only is limited to the matrix phase of composite materials, so it is considered that the damage is limited to E 22 and G 12 , or only degrades the matrix material in the RUC. Both E 22 and G 12 are functions of the internal state variable S, which can be expressed as, where, E 220 and G 120 represent the original transverse elastic modulus and shear modulus, and e s (S) and g s (S) represent the polynomial function of the transverse and shear modulus with respect to S. According to the plane stress constitutive relation, the expression of elastic strain energy density W can be expressed as, By assuming Q 12 is constant, and differentiating Equation (10) with respect to S yields, Experiments show that variable S is a function of ε 3 [25], so a reduced internal state variable S γ can be introduced, Then, Equation (11) becomes Once S γ is solved by Equation (13), the elastic moduli E 22 and G 12 can undergo nonlinear degradation according to Equations (8) and (9).

Generalized Method of Cells
The generalized method of cells (GMC) is employed to predict the strengths of continuous fiber composites in this work. The fiber reinforced composite is modeled as a rectangular repeating unit cell that contains a preset number of rectangular subcells as shown in Figure 2, and each subcell (βγ) may be occupied by a distinct homogeneous material. The constitutive equation in each subcell for the micromechanical model is denoted below, where σ (βγ) , C (βγ) , and ε (βγ) represent the average stress components, the elastic stiffness matrix, and the average total strain components in subcell (βγ), respectively. Inelastic strain and thermal strain are not took into considered herein. The basic displacement assumptions in GMC is a linear polynomial, as shown in Formula (15), with the local coordinates (x 3 ) located at the center of each subcell, where W (βγ) i(00) represents the displacements at the center of subcell (βγ), microvariables W The basic displacement assumptions in GMC is a linear polynomial, as shown mula (15), with the local coordinates (     The formulation of GMC involves imposition of several governing conditions, including application of continuity of displacements and tractions at each of the subcell interfaces and the periodic boundaries of the repeating unit cell (RUC) in an average sense. By regarding the stress components in the subcells as the basic unknowns, a system of equations relating the unknown stress variables to the macroscopic strains can be represented as follows, where the G matrix consists of information on the geometry and mechanical properties of the material in the individual subcells (βγ), the T vector represents the subcell stress components that need to be solved, and the f m vector incorporates information of geometrical dimensions and the global strains on the RUC. The expression that relates the average strain in the subcell to the externally global strain in virtue of the solutions of Equation (16), then the final form of the global constitutive equation that relates the average stress σ and strain ε is determined as follow, where C* is the average or effective stiffness matrix. The global stress can be established in the GMC repeating unit cell as follow, where h and l represent the height and length of the RUC, h β and l γ represent the height and length of the subcell (βγ), as shown in Figure 2. The GMC model is applied to each integration point of each finite element, and the strain states of each integral point are taken as the load input to the RUC, then the stress and strain states of each constitutive materials can be obtained according to Equation (16). The damage evolution and failure are evaluated at the micro-level, and after updating the stiffness matrix using the homogenization method, the multi-scale analysis will entry the next finite element iterative calculation. Figure 3 is a schematic diagram of the analysis process. Aerospace 2022, 9, x FOR PEER REVIEW 6 of 16

Analysis Method
Three analysis methods are adopted in this paper, one is a multi-scale analysis method based on GMC, the second method is based on ST at the lamina level, and the last one employs an empirical method based on the 2D Hashin-Rotem failure criteria and implicit solver of static analysis is used. The maximum stress failure criterion for fiber subcells and ST theory for matrix subcells are employed for the first method, the maximum stress failure criterion together with a degradation limit are used for the second method. Once the reduced internal state variable Sγ is determined, the transverse and shear stiffness of matrix properties degrades according to the damage functions es and gs. These damage and failure models for the static analysis are all implemented in ABAQUS/Standard employing the user-defined subroutine UMAT.

Lamina Level Damage and Failure Model Based on ST
According to the experimental data of transverse elastic modulus E22 and shear modulus G12 corresponding to internal variable Sγ in literature [26], damage functions es and gs are obtained by cubic polynomial fitting, as shown in equations below. The coefficients are shown in Table 1. Figure 4 shows the polynomial es and gs obtained by fitting and the experimental data, respectively. The fitting effect is very ideal. The damage evolution rate of G12 is obviously higher than that of E22.

Coefficients
Values G12 Coefficients Values

Analysis Method
Three analysis methods are adopted in this paper, one is a multi-scale analysis method based on GMC, the second method is based on ST at the lamina level, and the last one employs an empirical method based on the 2D Hashin-Rotem failure criteria and implicit solver of static analysis is used. The maximum stress failure criterion for fiber subcells and ST theory for matrix subcells are employed for the first method, the maximum stress failure criterion together with a degradation limit are used for the second method. Once the reduced internal state variable S γ is determined, the transverse and shear stiffness of matrix properties degrades according to the damage functions e s and g s . These damage and failure models for the static analysis are all implemented in ABAQUS/Standard employing the user-defined subroutine UMAT.

Lamina Level Damage and Failure Model Based on ST
According to the experimental data of transverse elastic modulus E 22 and shear modulus G 12 corresponding to internal variable S γ in literature [26], damage functions e s and g s are obtained by cubic polynomial fitting, as shown in equations below. The coefficients are shown in Table 1. Figure 4 shows the polynomial e s and g s obtained by fitting and the experimental data, respectively. The fitting effect is very ideal. The damage evolution rate of G 12 is obviously higher than that of E 22 .  According to Equation (13), the value of Sγ depends on the strain state of the repeating unit cell, that is, it is determined by the transverse strain ε22 and shear strain γ12. Figure 5 show the changing curves of Sγ with ε22 in the interval [0, 0.035] when γ12 is selected to be 0, 0.0025, 0.005, 0.0075, and 0.01. It can be seen that the positive solution of Sγ is unique and monotonically increasing. The growth rate increases first, then decreases, and finally flattens out. The values of Sγ are between 0 and 1.2. It is generally believed that fiber fracture occurs instantaneously and results in immediate local failure of laminates. Therefore, the properties of fiber is deemed to be linear before it reaches the ultimate failure, so the strength prediction can be satisfied by using a simple maximum stress failure criterion. , 0 where σ11 is the stress in the fiber direction, XT and XC are the fiber direction tensile and compressive strengths. Fiber failure will initiate when df exceed 1.
For the matrix damage of monolayer plate, the degradation of E22 and G12 for a laminate or the RUC is reduced by Equations (8) and (9). Since material softening may cause According to Equation (13), the value of S γ depends on the strain state of the repeating unit cell, that is, it is determined by the transverse strain ε 22 and shear strain γ 12 . Figure 5 show the changing curves of S γ with ε 22 in the interval [0, 0.035] when γ 12 is selected to be 0, 0.0025, 0.005, 0.0075, and 0.01. It can be seen that the positive solution of S γ is unique and monotonically increasing. The growth rate increases first, then decreases, and finally flattens out. The values of S γ are between 0 and 1.2.  According to Equation (13), the value of Sγ depends on the strain state of the repeating unit cell, that is, it is determined by the transverse strain ε22 and shear strain γ12. Figure 5 show the changing curves of Sγ with ε22 in the interval [0, 0.035] when γ12 is selected to be 0, 0.0025, 0.005, 0.0075, and 0.01. It can be seen that the positive solution of Sγ is unique and monotonically increasing. The growth rate increases first, then decreases, and finally flattens out. The values of Sγ are between 0 and 1.2. It is generally believed that fiber fracture occurs instantaneously and results in immediate local failure of laminates. Therefore, the properties of fiber is deemed to be linear before it reaches the ultimate failure, so the strength prediction can be satisfied by using a simple maximum stress failure criterion. , 0 where σ11 is the stress in the fiber direction, XT and XC are the fiber direction tensile and compressive strengths. Fiber failure will initiate when df exceed 1.
For the matrix damage of monolayer plate, the degradation of E22 and G12 for a laminate or the RUC is reduced by Equations (8) and (9). Since material softening may cause It is generally believed that fiber fracture occurs instantaneously and results in immediate local failure of laminates. Therefore, the properties of fiber is deemed to be linear before it reaches the ultimate failure, so the strength prediction can be satisfied by using a simple maximum stress failure criterion.
where σ 11 is the stress in the fiber direction, X T and X C are the fiber direction tensile and compressive strengths. Fiber failure will initiate when d f exceed 1.
For the matrix damage of monolayer plate, the degradation of E 22 and G 12 for a laminate or the RUC is reduced by Equations (8) and (9). Since material softening may cause convergence problems for the finite element model using the implicit algorithm, the degradation ratio of E 22 and G 12 was limited to 0.4, after which the damage evolution was terminated and E 22 and G 12 would remain constant, and the matrix material is considered to be failed.

Micro-Scale Damage and Failure Model
The micromechanical GMC model can capture the interaction between the constituents of composite materials, and can easily analyze the damage or failure of fiber and matrix phases independently. Each material point in the finite element model corresponds to a 2 × 2 RUC [21], as shown in Figure 6, containing one fiber subcell and three matrix subcells. The fiber volume fraction of T800/3900-2 lamina is 0.6, and its properties of fiber and matrix constituents are shown in Tables 2 and 3. convergence problems for the finite element model using the implicit algorithm, the degradation ratio of E22 and G12 was limited to 0.4, after which the damage evolution was terminated and E22 and G12 would remain constant, and the matrix material is considered to be failed.

Micro-Scale Damage and Failure Model
The micromechanical GMC model can capture the interaction between the constituents of composite materials, and can easily analyze the damage or failure of fiber and matrix phases independently. Each material point in the finite element model corresponds to a 2 × 2 RUC [21], as shown in Figure 6, containing one fiber subcell and three matrix subcells. The fiber volume fraction of T800/3900-2 lamina is 0.6, and its properties of fiber and matrix constituents are shown in Tables 2 and 3.   Fiber subcells are considered to be linear before the failure occurs, and the fiber subcell failure is determined by the maximum stress failure criterion, which is shown in Equations (23) and (24). The material properties of fiber constituent after failure are immediately reduced to 0.01% of its initial value.  Fiber subcells are considered to be linear before the failure occurs, and the fiber subcell failure is determined by the maximum stress failure criterion, which is shown in Equations (23) and (24). The material properties of fiber constituent after failure are immediately reduced to 0.01% of its initial value. σ (11) 11 where σ (11) 11 is the stress of the fiber subcell in the longitudinal direction, X fc and X fc are its tensile and compressive strengths. When damage variable d f is greater than or equal to 1, fiber subcells would fail.
In each time step of the multiscale analysis, the GMC micromechanics model calculates the overall elastic moduli should be consistent with the laminate level moduli obtained from Equations (8) and (9). Therefore, the matrix subcells should be degraded in a certain way in order to produce the same E 22 and G 12 . Matrix constituent properties are also functions of the reduced internal state variable S γ , The value of S γ is determined by the following formula, where ε (βγ) 22 and γ (βγ) 12 represent transverse strain and shear strain of matrix subcell (βγ). Matrix subcell stress vector σ (βγ) can be obtained by the following formula, A is a third-order matrix connecting the macroscopic strain components and the stress components of each subcell under the plane stress state.
Using the constitutive relation of the subcell material, there is where S (βγ) is the flexibility matrix of the subcell (βγ) material. Combined with Equations (28) and (29), the expressions of strain of each subcell are as follows: Due to the different stress-strain states of each matrix subcell at the same material point in the analysis process, the damage evolution rates of different matrix subcells are different. In order to keep the same as the damage limit 0.4 of lamina level E 22 and G 12 , the matrix subcell material with damage evolution dropping to 0.4 earlier stops, degrading until the last matrix subcell reaches the degradation limit, which represents the final matrix failure of the RUC. The flow chart of multi-scale analysis algorithm based on a four-subcells GMC model is as shown in Figure 7.

Lamina-Scale Failure Model-Based Hashin-Rotem Theory
The Hashin-Rotem failure theory include separate criteria for fiber and matrix. The fiber failure criteria are the same as the maximum stress criteria as shown in Equations (21) and (22). Matrix failure are contributed by the transverse and shear stresses in a lamina, which are dictated by where σ 22

Lamina-Scale Failure Model-Based Hashin-Rotem Theory
The Hashin-Rotem failure theory include separate criteria for fiber and matrix. Th fiber failure criteria are the same as the maximum stress criteria as shown in Equation (21) and (22). Matrix failure are contributed by the transverse and shear stresses in a lam ina, which are dictated by where σ22 and τ12 are the transverse and shear stresses in a lamina, YT, YC, and S are re spectively the transverse compressive strength, transverse tensile strength, and shea strength in the composite lamina.    Table 4.  Figure 8.

Experimental Result
The tensile test process of open-hole laminates is conducted in accordance with ASTM D5766-standard test method for open-hole tensile strength of polymer matrix composite laminates. The open-hole tensile test of laminates was carried out on MTS 370 electro-hydraulic servo test machine, which has a maximum load-carrying capacity of 250 kN with a accuracy of 0.5%. The displacement load curves of the test pieces are shown in Figure 8. The average failure load of all specimens is 38.34 kN, and the dispersion coefficient is 2.13%. It can be seen from the displacement versus load curve that all the specimens are obviously brittle before the ultimate failure determined by the tensile strength of 0° layer onset. The photograph of the failed specimen is present in Figure 9, it can be observed that the catastrophic failures are concentrated around the open hole which causes the stress concentrations and finally span the width of the specimen, the failure modes of the specimen are mainly as follows: tensile fracture of fiber for 0° layer and tensile fracture of the matrix for ±45° layer.  The average failure load of all specimens is 38.34 kN, and the dispersion coefficient is 2.13%. It can be seen from the displacement versus load curve that all the specimens are obviously brittle before the ultimate failure determined by the tensile strength of 0 • layer onset. The photograph of the failed specimen is present in Figure 9, it can be observed that the catastrophic failures are concentrated around the open hole which causes the stress concentrations and finally span the width of the specimen, the failure modes of the specimen are mainly as follows: tensile fracture of fiber for 0 • layer and tensile fracture of the matrix for ±45 • layer.
The tensile test process of open-hole laminates is conducted in accordance with ASTM D5766-standard test method for open-hole tensile strength of polymer matrix composite laminates. The open-hole tensile test of laminates was carried out on MTS 370 electro-hydraulic servo test machine, which has a maximum load-carrying capacity of 250 kN with a accuracy of 0.5%. The displacement load curves of the test pieces are shown in Figure 8. The average failure load of all specimens is 38.34 kN, and the dispersion coefficient is 2.13%. It can be seen from the displacement versus load curve that all the specimens are obviously brittle before the ultimate failure determined by the tensile strength of 0° layer onset. The photograph of the failed specimen is present in Figure 9, it can be observed that the catastrophic failures are concentrated around the open hole which causes the stress concentrations and finally span the width of the specimen, the failure modes of the specimen are mainly as follows: tensile fracture of fiber for 0° layer and tensile fracture of the matrix for ±45° layer.

Results of Finite Element Analysis
The three methods implemented using the subroutine UMAT introduced in Section 2 above are used to conduct the finite element analysis of the open-hole laminate specimens. The S4R mesh type is used, and the mesh is refined around the open hole. The composite constituents properties can be seen in Tables 2 and 3, and the laminate properties are shown in Table 4. Additionally, the tensile and compressive strengths of the composite laminate are 2690 MPa and 1380 MPa, respectively, which are provided by vendors and used to perform FEA at the micro-level. The finite element model is restricted from moving in the thickness direction and all the in-plane rotational degrees of freedom to account for plane stress simulation. The left end of the model is fixed and the right edge is constrained in the transverse displacement and all rotations, then a tensile displacement of 5 mm is applied to the right end. The load versus displacement curves obtained from the three progressive damage methods are shown in Figure 10. It can be observed that the ultimate load prediction for ST/FF (38.76 kN), Hashin-Rotem (39.49 kN) and ST/GMC/FF (38.93 kN) agree well with the experimental value (38.45 kN). The load displacement curve obtained by Hashin-Rotem progressive damage method is linear, but both ST/FF and ST/GMC methods achieve satisfactory nonlinear reduction before the displacement-load curve reaches its peak. The load displacement curve obtained in the experiment has shown, to some extent, nonlinear behavior, especially near the end of the experimental curve. In the experiment, the crack initiates around the open hole and propagates rapidly to the edge of the specimen. It is not possible to capture the rapid progression of the crack by using the implicit solver, so the gradual failure of fiber and matrix yield the exhibited softening in the ST or ST/GMC model.
for plane stress simulation. The left end of the model is fixed and the right edge is constrained in the transverse displacement and all rotations, then a tensile displacement of 5 mm is applied to the right end. The load versus displacement curves obtained from the three progressive damage methods are shown in Figure 10. It can be observed that the ultimate load prediction for ST/FF (38.76 kN), Hashin-Rotem (39.49 kN) and ST/GMC/FF (38.93 kN) agree well with the experimental value (38.45 kN). The load displacement curve obtained by Hashin-Rotem progressive damage method is linear, but both ST/FF and ST/GMC methods achieve satisfactory nonlinear reduction before the displacementload curve reaches its peak. The load displacement curve obtained in the experiment has shown, to some extent, nonlinear behavior, especially near the end of the experimental curve. In the experiment, the crack initiates around the open hole and propagates rapidly to the edge of the specimen. It is not possible to capture the rapid progression of the crack by using the implicit solver, so the gradual failure of fiber and matrix yield the exhibited softening in the ST or ST/GMC model.  Figure 11 shows the fiber damage and matrix nonlinear damage state of 0° layer and 45° layer under different loads by ST/FF method. As the damage states of the −45° layer and 45° layer are basically symmetric with respect to the 0° axis, the damage contour of 0° layer and 45° layer are selected to be shown below. At 80% ultimate load, tensile and shear damage begin to fail at the upper and lower edges of the open hole, while no element in the transverse direction reaches the maximum failure value for 0° layer. Under this load, a small number of elements have failed in the transverse direction and more shear failure are presented around the open hole for the 45° layer, but no fiber failure occurs in this layer. At 90% ultimate load, the fiber failure, the progressive degradation of transverse and shear moduli for 0° layer and 45° layer all continue to extend. When reaching the ultimate load, all damage contours of 0° or 45° layer have spanned the width of the specimen except for the fiber failure of 45° layer, after which the crack propagates rapidly until    Figure 12. The shear moduli of first and third matrix subcells at 0 • layer evolve to the limit value 0.4 for most elements, as shown in the third row of Figure 12, while the damage of transverse elastic moduli evolve along the direction of x 2 axis, which does not span the finite element model. It can be seen that the elastic modulus of second matrix subcell degrade faster than the other two matrix subcells, because the second matrix subcell locates at the same column with the fiber subcell, which have a large value of stiffness component according to the stress continuity assumption of GMC, i.e., σ 22 in the same column of RUC are equal on the x 2 axis shown in Figure 2. This phenomenon is very obvious at 0 • layer, but not at 45 • layer, indicating that the reduction evolution of matrix subcells is sensitive to the layer angle. It can be observed the degree of damage for shear modulus are greater than that of transverse modulus in both 0° and 45° layers, which is consistent with the theoretical model shown in Figure 4 above. The final damage patterns of 0° or 45° layer are in agreement with the experimental results, which proves that the prediction of failure mode is reliable.   Figure 12 shows the nonlinear damage contours of three matrix subcells for 0° layer and 45° layer under both 16.7 kN and the ultimate loads using ST/GMC/FF method. Under the loading of 16.7 kN, the shear modulus of the second matrix subcell at 0° layer begins to fail for a few elements at the upper and lower edges of the open hole, while none of the damage of transverse modulus of any element reaches its failure value. The other two matrix subcells for 0° layer have the same damage evolution and do not reach the lower limit of damage 0.4. Under this load, the shear elastic moduli gs of the three matrix subcells for 45° layer have been extended into a cross shape of ±45°, however, just a few transverse failure appears at the upper and lower edges of the open hole for the matrix subcells. The transverse and shear modulus progressive damage patterns are basically the same. Under the ultimate load, the transverse elastic modulus and shear modulus of all subcells in 45° layer have evolved to the failure value in almost all finite elements around the hole, which is not shown Figure 12. The shear moduli of first and third matrix subcells at 0° layer evolve to the limit value 0.4 for most elements, as shown in the third row of Figure 12, while the damage of transverse elastic moduli evolve along the direction of x2 axis, which does not span the finite element model. It can be seen that the elastic modulus of second matrix subcell degrade faster than the other two matrix subcells, because the second matrix subcell locates at the same column with the fiber subcell, which have a large value of stiffness component according to the stress continuity assumption of GMC, i.e., σ22 in the same column of RUC are equal on the x2 axis shown in Figure 2. This phenomenon is very obvious at 0° layer, but not at 45° layer, indicating that the reduction evolution of matrix subcells is sensitive to the layer angle.   As shown in Figure 13, the damage state of 0° layered fiber subcells is similar to the damage contour obtained by ST/FF method, in which the failure starts from the elements at the upper and lower edges of the hole and expands along the oblique crossing direction, until it runs horizontally through the finite element model. Because the RUC model used contains only one fiber subcell, the failure of the fiber subcell means the fiber failure for the whole RUC at the finite element scale. As shown in Figure 13, the damage state of 0 • layered fiber subcells is similar to the damage contour obtained by ST/FF method, in which the failure starts from the elements at the upper and lower edges of the hole and expands along the oblique crossing direction, until it runs horizontally through the finite element model. Because the RUC model used contains only one fiber subcell, the failure of the fiber subcell means the fiber failure for the whole RUC at the finite element scale.
As shown in Figure 13, the damage state of 0° layered fiber subcells is similar to the damage contour obtained by ST/FF method, in which the failure starts from the elements at the upper and lower edges of the hole and expands along the oblique crossing direction, until it runs horizontally through the finite element model. Because the RUC model used contains only one fiber subcell, the failure of the fiber subcell means the fiber failure for the whole RUC at the finite element scale.

Conclusions
By where the ST method well captured the nonlinear behavior of the laminates and reasonable failure mode prediction was obtained. The ST method makes the matrix damage become a natural nonlinear damage evolution and replaces the traditional matrix failure criterion model. Fiber tensile failure shows obvious brittleness behavior, so the failure criterion of fiber is still needed in both macroscopic and microscopic models. 2. Based on the multi-scale progressive damage model of GMC model, ST/GMC/FF methods will resolve stresses or strains of a laminate to the repeating unit cell at the microscopic level, which make the matrix and fiber constituents to be distinguishable, and different matrix subcells also have independent damage evolvement and

Conclusions
By using ST/FF, Hashin progressive damage model and ST/GMC/FF multi-scale progressive damage model, the ultimate load, failure mode and nonlinear damage evolution of fiber reinforced composite laminates with a central open hole were analysed, and the tensile experiment results of 8 specimens were compared. The following conclusions can be drawn:

1.
ST/FF nonlinear damage model, Hashin-Rotem progressive damage analysis method, and the ST/GMC/FF model based on the multi-scale progressive damage model are ideal to predict the damage state, failure pattern and the ultimate load, where the ST method well captured the nonlinear behavior of the laminates and reasonable failure mode prediction was obtained. The ST method makes the matrix damage become a natural nonlinear damage evolution and replaces the traditional matrix failure criterion model. Fiber tensile failure shows obvious brittleness behavior, so the failure criterion of fiber is still needed in both macroscopic and microscopic models.

2.
Based on the multi-scale progressive damage model of GMC model, ST/GMC/FF methods will resolve stresses or strains of a laminate to the repeating unit cell at the microscopic level, which make the matrix and fiber constituents to be distinguishable, and different matrix subcells also have independent damage evolvement and failure states. It shows that the multi-scale method based on the GMC model can better explain the crack path and failure mechanism for the matrix constituent of a laminate.