A Continuum Damage Model for Intralaminar Progressive Failure Analysis of CFRP Laminates Based on the Modified Puck’s Theory

A continuum damage model is proposed to predict the intralaminar progressive failure of CFRP laminates based on the modified Puck’s theory. Puck’s failure criteria, with consideration of the in situ strength effect, are employed to evaluate the onset of intralaminar failure including fiber fracture and inter-fiber fracture. After damage initiation, a bilinear constitutive relation is used to describe the damage evolution process. In strict accordance with Puck’s concept of action plane, the extent of damage is quantified by the damage variables defined in the fracture plane coordinate system, rather than the traditional material principal coordinate system. Theoretical and experimental evaluation of CFRP laminates under different loading conditions demonstrates the rationality and effectiveness of the proposed numerical model. The model has been successfully implemented in a finite element (FE) software to simulate the intralaminar progressive failure process of CFRP laminates. A good agreement between the experimental and numerical results demonstrates that the present model is capable of predicting the intralaminar failure of CFRP laminates.


Introduction
Carbon fiber reinforced polymers (CFRPs) are being increasingly used in industry due to their advantageous properties such as high specific strength and stiffness, good resistance to fatigue and corrosion, as well as flexibility in design. CFRP laminates are now widely applied in various fields, including aerospace structures, windmill blades, and pressure vessels. Nevertheless, the design of CFRP composite structures is still rather conservative in engineering practice. A large quantity of time-consuming and expensive tests must be carried out to ensure structure safety. To efficiently reduce the time and cost as well as fully exploit the advantages of CFRPs, there is an imperative need for the accurate theoretical prediction of the failure of CFRP laminates. However, the failure mechanisms of CFRP laminates are very complex due to their inherent anisotropy and the variety of failure modes. Hence, developing a reliable failure theory for CFRP laminates is a very challenging task.
The failure process of CFRP laminates can be divided into two stages. First the damage initiates in a ply, and then it gradually propagates through the laminate until the structure reaches the ultimate failure load [1]. The continuum damage mechanics (CDM) approach, which was originally developed by Kachanov [2], has been successfully employed in the progressive failure analysis of CFRP laminates Despite its good performance for unidirectional laminae, Puck's theory predicts an initial failure stress of laminates that is significantly lower than experimental results [16]. The underlying reason is that Puck's theory does not consider the "in situ" effect properly. When confined to a multidirectional laminate, the lamina's IFF strengths are observed to be higher than those measured in isolated unidirectional laminae due to the constraint of neighboring plies with different fiber orientation [17]. However, directly replacing the strength values in Puck's IFF criteria with the corresponding in situ strengths may result in incorrect predictions for both the first ply failure load and the fracture angle [18]. Hence Puck's IFF criteria should be re-examined and modified in order to take into account the in situ effect properly.
Once damage initiates, a stiffness degradation law is required to characterize damage propagation. The direct stiffness degradation method can be easily implemented in a finite element code, but it is purely empirical and lacks generality [19]. The energy-based CDM approach progressively degrades material properties until enough energy is dissipated for complete failure Figure 1. Definition of the material principal coordinate system, 1-2-3, and the fracture plane coordinate system, l-n-t.
Despite its good performance for unidirectional laminae, Puck's theory predicts an initial failure stress of laminates that is significantly lower than experimental results [16]. The underlying reason is that Puck's theory does not consider the "in situ" effect properly. When confined to a multidirectional laminate, the lamina's IFF strengths are observed to be higher than those measured in isolated unidirectional laminae due to the constraint of neighboring plies with different fiber orientation [17]. However, directly replacing the strength values in Puck's IFF criteria with the corresponding in situ strengths may result in incorrect predictions for both the first ply failure load and the fracture angle [18]. Hence Puck's IFF criteria should be re-examined and modified in order to take into account the in situ effect properly.
Once damage initiates, a stiffness degradation law is required to characterize damage propagation. The direct stiffness degradation method can be easily implemented in a finite element code, but it is purely empirical and lacks generality [19]. The energy-based CDM approach progressively degrades material properties until enough energy is dissipated for complete failure [20]. It associates the damage variables which represent the possible damage modes, with their respective fracture energies. A bilinear or an exponential constitutive relation is often established to describe the failure process. While many progressive damage models of CFRP laminates reduce the stiffness coefficients directly in the material principal coordinate system [21][22][23], a few researchers have noticed the influence of the fracture plane orientation on material degradation [5,24]. As shown in Figure 1, the external stress action plane may not coincide with the fracture plane in the case of inter-fiber failure. Hence, the material principal coordinate system (coordinate 1-2-3) generally differs from the fracture plane coordinate system (coordinate l-n-t). According to Puck's failure theory, inter-fiber fracture is exclusively determined by the stresses acting on the fracture plane [25]. After damage occurs, the load carrying capacity on the fracture plane will be reduced directly. Therefore, it is more reasonable to define the damage variables in the fracture plane coordinate system instead of the traditional material principal coordinate system.
In the present study, a continuum damage model for intralaminar progressive failure analysis of CFRP laminates is developed based on Puck's fracture plane theory. The in situ strength effect and the shear nonlinear behavior of CFRPs are considered in the model. The modified Puck's failure criteria are employed to predict the initiation of fiber and inter-fiber fracture. In strict accordance with Puck's failure hypothesis, the material stiffness matrix is degraded using the damage variables defined in the fracture plane coordinate system rather than the traditional material principal coordinate system. Theoretical and experimental evaluation validates the rationality and effectiveness of the proposed model.

Stress and Strain Analysis
The intralaminar damage modes of CFRP laminates normally can be divided into two categories: fiber fracture (FF) and inter-fiber fracture (IFF). For the FF mode, the failure plane is approximately the plane on which the longitudinal normal stress σ 1 acts [8], while the IFF mode occurs in an inclined plane parallel to the fiber direction; see Figure 1. The stress and strain components in the fracture plane coordinate system l-n-t can be transformed from those in the material principal coordinate system using the coordinate transformation matrix T: where σ = [σ 1 σ 2 σ 3 τ 23 τ 31 τ 12 ], ε = [ε 1 ε 2 ε 3 γ 23 γ 31 γ 12 ], σ fp = [σ l σ n σ t τ nt τ lt τ nl ], ε fp = [ε l ε n ε t γ nt γ lt γ nl ] are the stress and strain components defined in the coordinate system 1-2-3 and l-n-t respectively. The inverse matrix of T is given by: where θ (−90 • ≤ θ ≤ 90 • ) is the potential inter-fiber fracture angle. The constitutive equations for the undamaged unidirectional lamina in the coordinate 1-2-3 and l-n-t are respectively given by: where S 0 and S f p 0 are the initial compliance matrices defined in the coordinates 1-2-3 and l-n-t, respectively.  (1) and (2) into Equation (5), we obtain: Comparing Equation (4) and Equation (6) results in: Once damage initiates, the CDM approach introduces a damage variable matrix D to establish the relationship between the effective stressσ and the actual stress σ. As is discussed in Section 1, it is more reasonable to define the effective stress in the fracture plane coordinate system l-n-t instead of the material principal coordinate system 1-2-3:σ where D fp is the damage variable matrix defined in the coordinate l-n-t, whose form is given by: where d i (i= l, n, t) and d ij (i, j= n, t, l) are the damage variables corresponding to different damage modes. Based on the assumption of energy equivalence, the form of the strain energy of the damaged material is identical to that of the undamaged material [26]: Substituting Equation (8) into Equation (10) results in: Comparing the left and right sides of Equation (11) gives: Following the similar deriving process of Equation (7), we obtain: Substituting Equations (7) and (13) into Equation (12), the compliance matrix of the damaged material in the coordinate 1-2-3, S d , is expressed as: where D is defined as the damage variable matrix in the coordinate 1-2-3: Finally, the constitutive equations for the damaged unidirectional lamina in the coordinate 1-2-3 is given by: CFRPs usually exhibit nonlinear shear response in the in-plane shear directions. The nonlinear shear constitutive model proposed by Hahn and Tsai [27] is adopted in this paper: where G 12 and G 13 are the initial shear moduli, and β is the shear nonlinearity factor. In the present study, the nonlinear terms in Equation (17) are linearized using the stable algorithm provided in Abaqus documentation [28]. Although the linearization algorithm does not consider the unloading process, it is adopted due to its simplicity.
Puck's failure criteria for inter-fiber fracture are as follows [13] (p. 74): Inter-fiber tensile fracture (IFFT): Inter-fiber compressive fracture (IFFC): Y T and Y C are the transverse tensile and compressive strengths respectively, while S L is the longitudinal shear strength. p t ⊥// , p c ⊥// , p t ⊥⊥ , p c ⊥⊥ are four inclination parameters of contour lines of Substituting Equations (18)- (20) into Equations (21) and (22), the inter-fiber fracture function, F IFF , can be written as the function of all stress components except σ 1 as well as the potential fracture angle θ: The θ-F IFF relation can be obtained under any given three-dimensional stress sate. Inter-fiber fracture is considered to occur once the maximum value of F IFF is equal to 1, and the actual fracture angle θ fp is the corresponding potential fracture angle: The maximum value of F IFF and the corresponding fracture angle θ fp can be searched numerically within the interval −90 • ≤ θ ≤ 90 • [17].
Given the in situ effect of multidirectional laminates, the in situ strength values, Y is T , Y is C , and S is L , should be used in Puck's IFF criteria. However, simply replacing the strength values in Equations (25)- (28) with the corresponding in situ strengths may result in incorrect predictions of both the failure stress and the corresponding fracture angle [18]. For example, the strength values of CFRP IM7/8552 are listed in Table 1. The following simple stress states are considered: transverse tensile , and longitudinal shear stress (τ 21 = S L or S is L ). As shown in Figure 2, for the ply in a unidirectional laminate, the maximum values of the failure function F IFF and the corresponding fracture angles are well predicted under all these simple stress states. For the ply embedded in a multidirectional laminate, Figure 3 shows that the predicted results are reasonable under transverse compression (see Appendix A for details) and longitudinal shear. Nevertheless, the maximum value of F IFF exceeds 1 under pure transverse tension, indicating that material failure has already occurred before σ 2 = Y is T . Besides, the corresponding tensile fracture angle is no longer 0 • , which is obviously wrong [18]. To reasonably combine Puck's IFF criteria with the in situ strength theory, the procedure to determine the seven unknown parameters, namely, the three strength parameters (R At ⊥ , R A ⊥⊥ and R A ⊥// ) and the four inclination parameters (p t ⊥// , p c ⊥// , p t ⊥⊥ and p c ⊥⊥ ), should be re-examined in detail. Table 1. Strength values of CFRP IM7/8552 [29].
Under the longitudinal shear stress state, the corresponding in situ strength is denoted as S is L , and the fracture angle θ ls f p = 0 • . According to Equations (18)- (20), the only stress component on the fracture plane is τ nl = S is L . Substituting (τ nt = 0, τ nl = S is L , σ n = 0) into Equation (21) results in: Under transverse compression, the corresponding in situ strength is denoted as is C Y . The compressive fracture angle is slightly above 50° for CFRPs [25]; thus, 51 c fp θ  is taken in the present study. According to Equations (18)- (20), the stress components on the fracture plane are τnt = sin cos Under transverse compression, the corresponding in situ strength is denoted as is C Y . The compressive fracture angle is slightly above 50° for CFRPs [25]; thus, 51 c fp θ  is taken in the present study. According to Equations (18)- (20), the stress components on the fracture plane are τnt = sin cos c fp θ is the fracture angle corresponding to the maximum value of the failure function FIFFC: Under transverse compression, the corresponding in situ strength is denoted as Y is C . The compressive fracture angle is slightly above 50 • for CFRPs [25]; thus, θ c f p = 51 • is taken in the present study. According to Equations (18)- (20), the stress components on the fracture θ c f p is the fracture angle corresponding to the maximum value of the failure function F IFFC : Solving Equations (32) and (33) together, we obtain: Under transverse tension, the corresponding in situ strength is denoted as Y is T , and the fracture angle θ t f p = 0 • . According to Equations (18)- (20), the only stress component on the fracture plane is The inclination parameters p t ⊥// and p c ⊥// can be derived experimentally [25]. Nevertheless, the inclination parameters p t ⊥⊥ and p c ⊥⊥ in Puck's IFF criteria are not determined by experiments. Their recommended values, Equation (28), are given based on mathematical reasons, rather than physical ones [25]. If the in situ strengths are taken into account, Equation (28) will lead to incorrect predictions under pure transverse tension (as shown in Figure 3). In the present study, the rigorously derived Equation (35) is employed to determine p c ⊥⊥ , while p t ⊥⊥ is obtained using the following method. According to Equations (18)-(20), the stress components on the potential fracture plane under pure transverse tensile stress (σ 2 > 0) are given by: The failure stress is Y is T , and the fracture angle θ t f p = 0 • . This implies the following conditions: It is easy to prove that Equation (38) can be satisfied. From Equation (39) we obtain: Substituting Equation (36) into Equation (40) results in: Puck et al. believed that p t ⊥⊥ and p c ⊥⊥ ought to be approximately of the same magnitude, and setting p t ⊥⊥ = p c ⊥⊥ will not lead to any unacceptable contradictions for unidirectional laminates [25]. Since the modified IFF criteria should also be applicable in the case of unidirectional laminates, we propose to minimize the difference between p t ⊥⊥ and p c ⊥⊥ without violation of Equation (41). In other words, p t ⊥⊥ can be obtained by solving the following problem: Puck suggested the use of a simple maximum stress formulation to predict fiber fracture in 1969, and believed it was sufficient for a preliminary analysis [13] (p. 37). More sophisticated FF criteria were developed afterwards to account for the transverse effect, requiring the measurements for the elastic modulus E 1f and Poisson's ratio υ 12f of the fibers [17]. Since these values are seldom provided in references, the simple maximum stress criteria are adopted in the present study: Fiber tensile fracture (FFT): Fiber compressive fracture (FFC): where X T and X C are the longitudinal tensile and compressive strengths, respectively.

Damage Evolution Law
Damage propagation is a process accompanied by the dissipation of energy. A fracture energy-based approach is employed in the present study to characterize damage evolution. Once damage initiation is predicted by the failure criteria, the material stiffness will be progressively degraded until enough energy is dissipated for complete failure.

Damage Variables
As shown in Figure 4, a bilinear constitutive relation is used to characterize the damage evolution process. The internal damage variable corresponding to each intralaminar failure mode is defined as: where ε eq,I , ε 0 eq,I and ε f eq,I are the equivalent strains in the current state, the damage initiation state, and the final failure state, respectively. All these equivalent strains will be defined precisely later. d I = 0 means the material is undamaged, while d I = 1 represents the complete failure of material.  The normal tensile stress on the inter-fiber fracture plane (σn > 0) tends to open the cracks, so no forces can be transmitted as a result of there being no contact between the crack faces. Nevertheless, cracks are closed under the normal compressive stress (σn < 0), and thus forces can be transmitted across the cracks [5]. Consequently, the damage variable IFFC d only has an effect on the shear moduli, and the damage variable with respect to the n-direction is given by: . n IFFT dd  (46) The normal tensile stress on the inter-fiber fracture plane (σ n > 0) tends to open the cracks, so no forces can be transmitted as a result of there being no contact between the crack faces. Nevertheless, cracks are closed under the normal compressive stress (σ n < 0), and thus forces can be transmitted across the cracks [5]. Consequently, the damage variable d IFFC only has an effect on the shear moduli, and the damage variable with respect to the n-direction is given by: The damage variable with respect to the l-direction is given by [21]: The damage variable d t represents the damage with respect to the t-direction. Since the t-direction is always perpendicular to the damage axes of IFF and FF failure (see Figure 1), no damage occurs along the t-axis [31], i.e., d t = 0.
The damage variables associated with the shear moduli G ij (i, j = n, t, l) can be expressed as: with

Equivalent Stress and Strain
According to Puck's fracture hypothesis, material failure is exclusively determined by the stress components acting on the fracture plane. For the IFF mode, the stress components on the fracture plane are τ nt , τ nl and σ n . The tensile normal stress (σ n > 0) promotes fracture in combination with the shear stresses τ nt and τ nl , while the compressive normal stress (σ n < 0) impedes material failure [25]. Therefore, the equivalent stress and strain for the IFF mode are defined as: where is the McCauley operator defined as x = (x + |x|)/2 for x ∈ . For the FF mode, the stress components on the fracture plane are τ 12 , τ 13 and σ 1 [8]. However, the contribution of the shear stress components τ 12 and τ 13 to fiber fracture is very small, and hence can be neglected [32]. Therefore, the equivalent stress and strain for the FF mode are defined as: The equivalent stress and strain in the damage initiation state, σ 0 eq,I and ε 0 eq,I , correspond to the equivalent stress and strain when the failure function F I equals 1.

Equivalent Strain in the Final Failure State
The crack band theory introduces the characteristic length L c to correlate the energy release rate (i.e., energy dissipated per unit area) with the energy dissipated per unit volume [33]. For the IFF failure mode, a quadratic interaction criterion is established under mixed-mode loading: where g n , g nt , and g nl are the strain energy densities associated with the corresponding stress components. G Ic and G IIc are the critical energy release rates for fracture mode I and mode II respectively. L c IFF is the characteristic length for IFF.
The strain energy densities in the final IFF failure state are given by [12]: where σ 0 i is the stress at the initiation of IFF failure. β i denotes the mixed-mode ratio, which can be expressed as: Substituting Equation (58) into Equation (57), the equivalent strain in the final IFF failure state is given by: For the FF mode, the equivalent strain in the final failure state is determined by the bilinear constitutive relation shown in Figure 4: where G t(c) f c is the critical energy release rate corresponding to longitudinal tension (compression). X T(C) is the longitudinal tensile (compressive) strength. L c FF is the characteristic length for FF.

Ply Failure Analysis
Two plies, one in a unidirectional laminate and the other embedded in a multidirectional laminate, are taken as the case study. The material system is IM7/8552, whose strength values are listed in Table 1.
As shown in Figure 5, the results predicted by the original and modified Puck's IFF criteria are almost identical for the ply in a unidirectional laminate. There are only negligible discrepancies between the curves due to the small differences between the parameters used (p t ⊥⊥ = p c ⊥⊥ = 0.258 and 0.262 respectively). The failure stresses and the fracture angles are correctly predicted under these simple stress states.
For the embedded lamina, the predicted results are reasonable under transverse compression (see Appendix A for details) and longitudinal shear; see Figure 6b,c. Under pure transverse tension, the failure stress is underestimated by the original IFF criteria, and the predicted fracture angle is unreasonable; see Figure 6a. This problem is solved by using the modified IFF criteria. Figure 5, the results predicted by the original and modified Puck's IFF criteria are almost identical for the ply in a unidirectional laminate. There are only negligible discrepancies between the curves due to the small differences between the parameters used ( = =0.258 tc pp   and 0.262 respectively). The failure stresses and the fracture angles are correctly predicted under these simple stress states.

As shown in
(a) (b) (c) For the ply in a unidirectional laminate, again there is a good agreement between the predicted σ2-τ21 failure envelopes; see Figure 7a. For the embedded lamina, Puck's original IFF criteria underestimate the failure stress in the high tensile stress (σ2 > 0) region (see Figure 7b), while the modified criteria do not have this serious defect.
For the embedded lamina, the predicted results are reasonable under transverse compression (see Appendix A for details) and longitudinal shear; see Figure 6b,c. Under pure transverse tension, the failure stress is underestimated by the original IFF criteria, and the predicted fracture angle is unreasonable; see Figure 6a. This problem is solved by using the modified IFF criteria.
For the ply in a unidirectional laminate, again there is a good agreement between the predicted σ2-τ21 failure envelopes; see Figure 7a. For the embedded lamina, Puck's original IFF criteria underestimate the failure stress in the high tensile stress (σ2 > 0) region (see Figure 7b), while the modified criteria do not have this serious defect.
For the ply in a unidirectional laminate, again there is a good agreement between the predicted σ 2 -τ 21 failure envelopes; see Figure 7a. For the embedded lamina, Puck's original IFF criteria underestimate the failure stress in the high tensile stress (σ 2 > 0) region (see Figure 7b), while the modified criteria do not have this serious defect.

Progressive Failure Analysis of CFRP Laminates
The proposed model focuses on the intralaminar damage of CFRP laminates, yet ignores the interlaminar damage between adjacent layers (i.e., delamination). Since the influence of delamination on the mechanical behavior of composite laminates with dispersed plies is very small under in-plane loading conditions [29,34], notched CFRP laminates under uniaxial tension and compression are selected as the validation cases. To avoid convergence problems, quasi-static analysis was performed in the FE software Abaqus/Explicit using a user-defined material subroutine (VUMAT). The FE model of the specimen is shown in Figure 8. The continuum 3D 8-node reduced integration element (C3D8R) was used per ply thickness. The mesh around the circular notch was refined (0.5 mm × 0.5 mm × t), while a relatively coarse mesh was used in the remaining regions (1.5 mm × 0.5 mm × t). Ideally, the characteristic length is a function of the fracture angle as well as the direction of crack propagation. For the sake of simplicity, the characteristic lengths for IFF and FF are here taken as the

Progressive Failure Analysis of CFRP Laminates
The proposed model focuses on the intralaminar damage of CFRP laminates, yet ignores the interlaminar damage between adjacent layers (i.e., delamination). Since the influence of delamination on the mechanical behavior of composite laminates with dispersed plies is very small under in-plane loading conditions [29,34], notched CFRP laminates under uniaxial tension and compression are selected as the validation cases. To avoid convergence problems, quasi-static analysis was performed in the FE software Abaqus/Explicit using a user-defined material subroutine (VUMAT). The FE model of the specimen is shown in Figure 8. The continuum 3D 8-node reduced integration element (C3D8R) was used per ply thickness. The mesh around the circular notch was refined (0.5 mm × 0.5 mm × t), while a relatively coarse mesh was used in the remaining regions (1.5 mm × 0.5 mm × t). Ideally, the characteristic length is a function of the fracture angle as well as the direction of crack propagation. For the sake of simplicity, the characteristic lengths for IFF and FF are here taken as the cubic root of the volume of the elements. This method has been proven to be efficient and effective for 3D elements if the element size is sufficiently small [24,35]. selected as the validation cases. To avoid convergence problems, quasi-static analysis was performed in the FE software Abaqus/Explicit using a user-defined material subroutine (VUMAT). The FE model of the specimen is shown in Figure 8. The continuum 3D 8-node reduced integration element (C3D8R) was used per ply thickness. The mesh around the circular notch was refined (0.5 mm × 0.5 mm × t), while a relatively coarse mesh was used in the remaining regions (1.5 mm × 0.5 mm × t). Ideally, the characteristic length is a function of the fracture angle as well as the direction of crack propagation. For the sake of simplicity, the characteristic lengths for IFF and FF are here taken as the cubic root of the volume of the elements. This method has been proven to be efficient and effective for 3D elements if the element size is sufficiently small [24,35].
Only half of the laminate was modelled due to the stacking symmetry in the z direction. The loading direction is along the longitudinal direction of the specimen, which coincides with the direction of the 0° ply. The laminate was clamped on the left end, while a uniform displacement was applied at the right edge. The z-symmetry boundary conditions were applied at the mid-plane. Two specimens were modelled in the present study. One is a quasi-isotropic laminate with the stacking sequence [90/0/±45]3s under unidirectional tensile loading [29], while the other is an angleply laminate containing six ±45° sub-laminates under uniaxial compression [36]. A sketch of the laminate is shown in Figure 9, and the geometric dimensions are reported in Table 2. A strain gauge was placed on the outer surface (lgau = 50 mm) to monitor the axial strain ε of the quasi-isotropic laminate, while an extensometer was used to measure the relative displacement Δ of the angle-ply laminate (lext = 25.4 mm). The laminates are manufactured from CFRP IM7/8552 and T300/976 respectively, whose mechanical properties are listed in Tables 3 and 4. The in situ strengths are calculated using the analytical formulas proposed by Camanho and his co-workers [30,37]; see Table  5. Only half of the laminate was modelled due to the stacking symmetry in the z direction. The loading direction is along the longitudinal direction of the specimen, which coincides with the direction of the 0 • ply. The laminate was clamped on the left end, while a uniform displacement was applied at the right edge. The z-symmetry boundary conditions were applied at the mid-plane.
Two specimens were modelled in the present study. One is a quasi-isotropic laminate with the stacking sequence [90/0/±45] 3s under unidirectional tensile loading [29], while the other is an angle-ply laminate containing six ±45 • sub-laminates under uniaxial compression [36]. A sketch of the laminate is shown in Figure 9, and the geometric dimensions are reported in Table 2. A strain gauge was placed on the outer surface (l gau = 50 mm) to monitor the axial strain ε of the quasi-isotropic laminate, while an extensometer was used to measure the relative displacement ∆ of the angle-ply laminate (l ext = 25.4 mm). The laminates are manufactured from CFRP IM7/8552 and T300/976 respectively, whose mechanical properties are listed in Tables 3 and 4. The in situ strengths are calculated using the analytical formulas proposed by Camanho and his co-workers [30,37]; see Table 5.        Table 5. Analytical formulas to calculate the in situ strengths.

Type of Ply
For the laminate under uniaxial tensile loading, a total of 110,784 elements were used in the model. The axial average stress is defined as the external load per unit cross-sectional area: σ = F/wt L , where w and t L are the width and thickness of the specimen. The axial average stress-strain curves are shown and compared in Figure 10. Good correlation between the experimental and numerical results is observed, and the predicted ultimate average stress (394.8 MPa) is slightly higher than the experimental value (375.7 MPa). Not considering interlaminar damage might be another possible reason for the overprediction apart from the uncertainty of experiments. As shown in Figure 11a, inter-fiber damage first occurs at the notch edge of the outer 90 • ply. Subsequently, fiber tensile fracture initiates in the 0 • plies; see Figure 11b. As the load increases, damage propagates perpendicular to the loading direction. In the final failure state of the laminate, fiber damage extends along the transverse direction in the 0 • plies, while inter-fiber damage in the 90 • plies fully extends across the width of the specimen; see Figure 11c,d. Inter-fiber tensile fracture also occurs in the ±45 • plies (see Figure 11e,f), but the damage zones are smaller than those in the 90 • plies. The predicted damage pattern is in line with the net-section failure mode observed in the experiment [29]. fracture initiates in the 0° plies; see Figure 11b. As the load increases, damage propagates perpendicular to the loading direction. In the final failure state of the laminate, fiber damage extends along the transverse direction in the 0° plies, while inter-fiber damage in the 90° plies fully extends across the width of the specimen; see Figure 11c,d. Inter-fiber tensile fracture also occurs in the ±45° plies (see Figure 11e,f), but the damage zones are smaller than those in the 90° plies. The predicted damage pattern is in line with the net-section failure mode observed in the experiment [29]. perpendicular to the loading direction. In the final failure state of the laminate, fiber damage extends along the transverse direction in the 0° plies, while inter-fiber damage in the 90° plies fully extends across the width of the specimen; see Figure 11c,d. Inter-fiber tensile fracture also occurs in the ±45° plies (see Figure 11e,f), but the damage zones are smaller than those in the 90° plies. The predicted damage pattern is in line with the net-section failure mode observed in the experiment [29]. For the laminate under uniaxial compressive loading, a total of 38,892 elements were used in the model. The F-Δ curve obtained from the experimental data is plotted in Figure 12, where F is the external load measured in the test, and Δ is the relative displacement. As shown in Figure 12, the curve exhibits a very pronounced nonlinear behavior, and the ultimate failure load is approximately 13.5 kN. The simulated curve agrees well with the experimental result, and the ultimate failure load is well predicted (12.6 kN) with a relative error of 6.7%. No fiber damage is observed in the whole failure process of the laminate. Only inter-fiber compressive damage occurs in the vicinity of the hole, and then propagates along the ±45° direction; see Figure 13. The predicted damage pattern is consistent with the experimental observation [36]. For the laminate under uniaxial compressive loading, a total of 38,892 elements were used in the model. The F-∆ curve obtained from the experimental data is plotted in Figure 12, where F is the external load measured in the test, and ∆ is the relative displacement. As shown in Figure 12, the curve exhibits a very pronounced nonlinear behavior, and the ultimate failure load is approximately 13.5 kN. The simulated curve agrees well with the experimental result, and the ultimate failure load is well predicted (12.6 kN) with a relative error of 6.7%. No fiber damage is observed in the whole failure process of the laminate. Only inter-fiber compressive damage occurs in the vicinity of the hole, and then propagates along the ±45 • direction; see Figure 13. The predicted damage pattern is consistent with the experimental observation [36].

Conclusions
In the present study, a continuum damage model based on the modified Puck's theory is developed to simulate the intralaminar progressive failure of CFRP laminates. The in situ strength effect and the nonlinear shear behavior of CFRPs are considered in the model. The modified Puck's failure criteria are adopted to determine damage initiation, while a bilinear constitutive relation is used to describe damage evolution. In strict accordance with Puck's concept of action plane, the equivalent stress/strain and the damage variables are defined in the fracture plane coordinate system rather than the traditional material principal coordinate system. Theoretical and experimental evaluation of CFRP laminates validates the rationality and effectiveness of the proposed model. The numerical model has been implemented in an FE software to simulate the progressive failure of CFRP laminates, and good correlation between the numerical and experimental results is observed. Future research will combine the proposed model with the interface fracture modeling techniques to simulate both intralaminar and interlaminar damage of CFRP laminates.

Conclusions
In the present study, a continuum damage model based on the modified Puck's theory is developed to simulate the intralaminar progressive failure of CFRP laminates. The in situ strength effect and the nonlinear shear behavior of CFRPs are considered in the model. The modified Puck's failure criteria are adopted to determine damage initiation, while a bilinear constitutive relation is used to describe damage evolution. In strict accordance with Puck's concept of action plane, the equivalent stress/strain and the damage variables are defined in the fracture plane coordinate system rather than the traditional material principal coordinate system. Theoretical and experimental evaluation of CFRP laminates validates the rationality and effectiveness of the proposed model. The numerical model has been implemented in an FE software to simulate the progressive failure of CFRP laminates, and good correlation between the numerical and experimental results is observed. Future research will combine the proposed model with the interface fracture modeling techniques to simulate both intralaminar and interlaminar damage of CFRP laminates.

Conclusions
In the present study, a continuum damage model based on the modified Puck's theory is developed to simulate the intralaminar progressive failure of CFRP laminates. The in situ strength effect and the nonlinear shear behavior of CFRPs are considered in the model. The modified Puck's failure criteria are adopted to determine damage initiation, while a bilinear constitutive relation is used to describe damage evolution. In strict accordance with Puck's concept of action plane, the equivalent stress/strain and the damage variables are defined in the fracture plane coordinate system rather than the traditional material principal coordinate system. Theoretical and experimental evaluation of CFRP laminates validates the rationality and effectiveness of the proposed model. The numerical model has been implemented in an FE software to simulate the progressive failure of CFRP laminates, and good correlation between the numerical and experimental results is observed. Future research will combine the proposed model with the interface fracture modeling techniques to simulate both intralaminar and interlaminar damage of CFRP laminates.

Appendix A
As shown in Figure A1, the in situ transverse compressive strength is C Y in Table 1 is increased from the unidirectional transverse compressive strength YC to 3YC, but the maximum value of FIFF always equals 1, and the fracture angle is also within the reasonable range (50°~54°). The reason is as follows: Under transverse compression, the corresponding in situ strength is denoted as Y is C . According to Equations (18)- (20), the stress components on the potential fracture plane are τ nt = Y is C sin θ cos θ, σ n = −Y is C cos 2 θ. Substituting (τ nt = Y is C sin θ cos θ, τ nl = 0, σ n = −Y is C cos 2 θ) into Equation (22) results in: holds in both the original and modified criteria. Substituting it into Equation (A1), we obtain: It can be seen from the above equation that under transverse compressive stress σ 2 = − Y is C , the F IFF -θ curve solely depends on the value of p c ⊥⊥ . The rigorously derived Equation (35) shows the relation between p c ⊥⊥ and θ c f p : In Puck's original criteria, p c ⊥⊥ is determined by Substituting the parameters of CFRP IM7/8552 (with the variance of Y is C ) into Equation (A5), the fracture angles predicted by Puck's original IFFC criterion can be obtained (50 •~5 4 • ).
The above analysis shows that even if the in situ transverse compressive strength increases significantly, the predicted result does not change much. Hence the results predicted by Puck's original IFFC criterion are reasonable under transverse compression.