Progressive Failure Analysis of Laminates with Embedded Wrinkle Defects Based on an Elastoplastic Damage Model

Out-of-plane wrinkling has a significant influence on the mechanical performance of composite laminates. Numerical simulations were conducted to investigate the progressive failure behavior of fiber-reinforced composite laminates with out-of-plane wrinkle defects subjected to axial compression. To describe the material degradation, a three-dimensional elastoplastic damage model with four damage modes (i.e., fiber tensile failure, matrix failure, fiber kinking/splitting, and delamination) was developed based on the LaRC05 criterion. To improve the computational efficiency in searching for the fracture angle in the matrix failure analysis, a high-efficiency and robust modified algorithm that combines the golden section search method with an inverse interpolation based on an existing study is proposed. The elastoplastic damage model was implemented in the finite-element code Abaqus using a user-defined material subroutine in Abaqus/Explicit. The model was applied to the progressive failure analysis of IM7/8552 composite laminates with out-of-plane wrinkles subjected to axial compressive loading. The numerical results showed that the compressive strength prediction obtained by the elastoplastic damage model is more accurate than that derived with an elastic damage model. The present model can describe the nonlinearity of the laminate during the damage evolution and determine the correct damage locations, which are in good agreement with experimental observations. Furthermore, it was discovered that the plasticity effects should not be neglected in laminates with low wrinkle levels.


Introduction
Fiber waviness is a type of manufacturing defect that occurs mostly during filament winding. Ply level out-of-plane waviness can result in severe degradation of mechanical properties, in particular, such as the compressive strength of composites. Hsiao and Daniel [1] conducted theoretical and experimental studies on unidirectional composites with out-of-plane wrinkles under compressive loading. They discovered that the stiffness and strength of the laminates decreases significantly with increasing fiber waviness. A similar study was conducted by Davidson and Waas [2], who found that for thick unidirectional carbon fiber polymer matrix composites, there exists a fiber misalignment angle at which the compressive strength is the global minimum. In a subsequent work, Davidson and Waas [3] developed a fiber waviness tolerance and criticality assessment framework employing surrogate modeling and Monte Carlo methods to predict the compressive strength and failure mode of a composite structure with fiber waviness. Adams et al. [4,5] investigated the effect of fiber wrinkles on the compressive strength of symmetric cross-ply laminates. According to their results, the compressive strength of those laminates was up to 36% lower than that of pristine laminates. Furthermore, various studies have shown that composites exhibit significant nonlinear behaviors before the final collapse of their structures, especially for laminates with wrinkles. Wisnom and Atkinson [6] performed finite-element (FE) analyses and experiments on T800/924 unidirectional laminates. Their results showed that fiber misalignments could induce shear nonlinearities. Chun et al. [7] discovered experimentally that DMS 2224 carbon/epoxy composite laminates with wrinkles exhibit both material and geometric nonlinearities under tension and compression. In addition, Makeev et al. [8] observed a nonlinear shear behavior in IM7/8552 laminates with fiber wrinkles using digital image correlation, and they proposed a nonlinear shear stress-strain relation obtained from numerical simulations. Davidson and Waas [9] used an odd polynomial series fitting obtained from experimental results to model the nonlinear shear response of the matrix in fiber-reinforced composites. Mukhopadhyay et al. [10] captured the progressive damage in IM7/8552 laminates with wrinkles using a high-speed video camera during compression experiments, and then, adopted a nonlinear shear stress-strain relation to describe the material behavior in a subsequent numerical analysis. However, the above studies focused only on shear nonlinearity, and disregarded plasticity and nonlinearity in the transverse direction.
In many progressive damage models developed for composite laminates, plasticity has been introduced to describe material nonlinearity. In a microscale, Prabhakar and Waas [11][12][13], Sun et al. [14], and Yuan et al. [15] applied a micromechanics model to predict the compressive failure behavior of unidirectional fiber-reinforced laminated composites using plasticity to approximate the nonlinearity of matrices. At the ply level, Lemanski et al. [16] presented a perfectly plastic model based on the Hill yield criterion and the kinematic hardening rule to simulate the nonlinear behavior of AS4/8552 composite laminates fabricated by Wang et al. [17], although only delamination damage was considered in their analysis. Wang et al. [18][19][20] performed off-axis tension and compression tests of IM600/Q133 unidirectional laminates and discovered that a material performs differently under tension and compression. Therefore, they proposed an elastoplastic constitutive model with distinguished tension and compression performances as an improvement on the Hill yield criterion. Xue et al. [21,22] applied this model to analyze the progressive elastoplastic failure behavior of IM600/Q133 and IM7/8552 composite laminates. Chen et al. [23] developed a combined elastoplastic damage model that accounted for plasticity on both the transverse and in-plane shear directions; however, the model could only be applied for two-dimensional damage analysis. According to an experiment, in which composite laminates were subjected to traverse compression [24,25], the fracture plane is not parallel to the loading direction. Moreover, fibers exhibit kinking or splitting failure under axial compression [26,27]. Therefore, existing damage models cannot be applied for the compression failure analysis of laminates considering plasticity effects.
To the best of our knowledge, no studies have been conducted on the longitudinal compression failure of multi-ply laminates with wrinkles that consider plastic damage in the matrix. Therefore, we extended the elastoplastic damage model described in [23] to a three-dimensional (3D) damage analysis and conducted a progressive failure simulation of composite laminates with out-of-plane wrinkle defects subjected to axial compression. To rapidly search for the fracture angle in the matrix failure analysis, we propose a highly accurate approach that combines the golden section search method and an inverse interpolation. We implemented the damage model based on the LaRC05 criterion [28] in a modified algorithm using a user-defined material subroutine in Abaqus/Explicit (VUMAT). To evaluate the effectiveness of the proposed method, we compared the predicted results with experimental observations [10].
The contributions of this study are as follows: (1) We conducted a progressive failure analysis of multidirectional fiber-reinforced polymer laminates with embedded wrinkle defects based on an elastoplastic damage model. (2) We demonstrated the nonlinearity of the laminate during damage evolution and correctly determined the damage location through numerical simulations. The remainder of the paper is organized as follows: Section 2 describes the elastoplastic damage model adopted in the present study. In Section 3, an approach is proposed to determine the orientation of the fracture plane in the matrix failure analysis. The implementation of the elastoplastic model using the user-defined subroutine VUMAT is presented in Section 4. Section 5 shows numerical simulation examples and discusses the predictions by comparing them with test results [10]. Finally, conclusions are presented in Section 6.

Stress-Strain Relationships
For composite materials exhibiting a plasticity response, the total strain tensor ε is expressed as the sum of the elastic and plastic strain parts, ε e and ε p , respectively, as follows: According to damage mechanics theories, the Cauchy nominal stress tensor σ and the effective stress tensorσ obey the following relationship: , where d represents the damage variable and d 1 , d 2 , and d 3 denote the damage in the fiber and in the transverse and shear directions, respectively.
To ensure the irreversibility of the damage, d 1 , d 2 , and d 3 are expressed as: where d ft is the damage caused by tension in the fiber direction; d kink and d split denote the fiber kinking and splitting, respectively, caused by compression in the fiber direction; d mt represents both the damage caused by tension in the matrix and the degradation of the interface between the fibers and the matrix due to decohesion; and d mc denotes the damage caused by compression in the matrix. The relationship between stress and strain for undamaged orthotropic anisotropy composites is as follows:σ =C 0 ε e (4) where C 0 is the stiffness tensor of the undamaged unidirectional laminated composite. By substituting Equation (4) into Equation (2), the relationship between the elastic strain tensor ε e and the Cauchy nominal stress tensor can be expressed as follows: where S 0 is the undamaged flexibility tensor. The stress-strain relationship for damaged composite materials can be expressed in the form: where S d is the damage flexibility tensor. By comparing Equations (5) and (6), the relationship between S 0 and S d can be expressed as follows: Materials 2020, 13, 2422 4 of 26 According to Matzenmiller's assumption [29], the compliance tensor of the damaged lamina can be obtained by adjusting the Poisson ratios. In the present study, the compliance tensor was degraded similar to [29]. Therefore, the following form of S d was adopted:

Plastic Model
The plastic yield function can be expressed in terms of effective stresses as follows [23]: where F p is the plastic potential and k is the hardening parameter. The power law proposed by Sun and Chen [30] is expressed as follows: where A and n are coefficients that can be determined using an approach based on the linear regression analysis of off-axis tensile tests performed on unidirectional composite laminate specimens. σ and ε p are the equivalent stress and the equivalent plastic strain, respectively. For simplicity, Chen et al. [23] converted the isotropic hardening law to an equivalent form in which the equivalent plastic strain ε p is used as an internal variable: where the coefficients β and m are related to A and n, respectively, through the relationships β = A −1/n and m = 1/n. Here, the 3D general plastic potential proposed by Sun and Chen [30] is employed. Hence, F p can be rewritten as: where a 66 is a material parameter that describes the level of plastic deformation developed under shear loading compared with that under transverse loading and σ ij represents the effective stress component. Assuming the associated plastic-flow rule and the associated hardening rule for composite materials, the plastic strain rate dε p ij can be expressed as: where λ is a plastic consistency parameter. The plastic work per unit volume dW p is defined as: By substituting Equations (12) and (13) into Equation (14), the equivalent strain rate is expressed as follows:

Damage Model
In general, the damage modes in composite structures include fiber failure, matrix cracking, and delamination. Fiber failure and matrix cracking occur in the plane and can be further categorized into tensile and compressive modes. Pinho [28] classified the fiber failure caused by compression into fiber splitting and fiber kinking.

Fiber Tensile Failure
Typically, fiber failure occurs when the longitudinal stress reaches the longitudinal tensile strength. Therefore, the maximum stress criterion is considered: where σ 11 is the normal stress in the fiber direction; X T is the axial tensile strength of the composite; and f ft is the exposure factor corresponding to the tension-induced fiber failure.

Matrix Failure
According to experimental observations [24,25], a fracture plane appears in the longitudinal direction of fibers under a transverse stress or/and an in-plane stress acts. Matrix failure occurs in the fracture plane. Puck introduced a criterion using stress components in the fracture plane. Pinho [28] adopted the advantages of the Puck criterion and improved it in the second World-Wide Failure Exercise (WWFE-II). The stress components on the fracture plane, as shown in Figure 1, can be obtained through a transformation of coordinates using the following formulas: where σ, τ, and τ L are the normal, longitudinal shear, and transverse shear stresses in the crack plane, respectively. φ denotes the angle of the fracture plane.
By substituting Equations (12) and (13) into Equation (14), the equivalent strain rate is expressed as follows:

Damage Model
In general, the damage modes in composite structures include fiber failure, matrix cracking, and delamination. Fiber failure and matrix cracking occur in the plane and can be further categorized into tensile and compressive modes. Pinho [28] classified the fiber failure caused by compression into fiber splitting and fiber kinking.

Fiber Tensile Failure
Typically, fiber failure occurs when the longitudinal stress reaches the longitudinal tensile strength. Therefore, the maximum stress criterion is considered: 11 11 1, 0 where 11  is the normal stress in the fiber direction; XT is the axial tensile strength of the composite; and fft is the exposure factor corresponding to the tension-induced fiber failure.

Matrix Failure
According to experimental observations [24,25], a fracture plane appears in the longitudinal direction of fibers under a transverse stress or/and an in-plane stress acts. Matrix failure occurs in the fracture plane. Puck introduced a criterion using stress components in the fracture plane. Pinho [28] adopted the advantages of the Puck criterion and improved it in the second World-Wide Failure Exercise (WWFE-II). The stress components on the fracture plane, as shown in Figure 1, can be obtained through a transformation of coordinates using the following formulas: where , , and L are the normal, longitudinal shear, and transverse shear stresses in the crack plane, respectively. ϕ denotes the angle of the fracture plane.  If the normal stress in the fracture plane is a tensile stress, i.e., σ N; ≥ 0, then matrix tensile failure will occur; otherwise, matrix compression failure will occur. The matrix failure criteria are expressed as follows: where Y T is the transverse tensile strength; µ L and µ T are the friction coefficients in the longitudinal and transverse directions, respectively; and S L and S T are the in situ longitudinal and transverse shear Materials 2020, 13, 2422 6 of 26 strengths, respectively. It is worth noting that, although S L and S T are not the same as the parameters S 12 and S 23 , respectively, used in physics, Puck et al. [31] found that the value of S L can be set to that of S 12 . The other parameters can be determined by analyzing the pure transverse compression of a composite laminate, which is calculated as follows [32]: where φ 0 is the angle of the fracture plane for pure compression [24,25], i.e., φ 0 = 53 ± 2 • , and Y C is the transverse compressive strength of the composite.

Fiber Compression Failure
Schultheisz and Waas [26] observed a local matrix deformation accompanied by fiber fracture (i.e., "fiber kinking") that differs from fiber microbuckling. Furthermore, Argon [27] assumed that initial microbuckling could result in fiber rotation, matrix shearing, fiber-matrix debonding, fiber kinking, and splitting in kink bands.
Although the fiber kinking mechanism is similar to that of matrix failure, a major difference is that in the matrix failure model, the associated stresses are calculated with respect to the fracture plane, whereas in fiber kinking analysis, stresses are calculated with respect to the kink plane, which is aligned with the fiber rotation. As shown in Figure 2, system 1-2-3 describes the material coordinates, where ψ is the angle between the kink plane and axis 2. Fiber kinking occurs on the plane with local coordinate system 1 m -2 m -3 m , which is obtained by rotating the coordinate system 1-2 ψ -3 ψ by an angle θ.
where YT is the transverse tensile strength; L and T are the friction coefficients in the longitudinal and transverse directions, respectively; and SL and ST are the in situ longitudinal and transverse shear strengths, respectively. It is worth noting that, although SL and ST are not the same as the parameters S12 and S23, respectively, used in physics, Puck et al. [31] found that the value of SL can be set to that of S12. The other parameters can be determined by analyzing the pure transverse compression of a composite laminate, which is calculated as follows [32]: where ϕ0 is the angle of the fracture plane for pure compression [24,25], i.e., ϕ0 = 53 ± 2°, and YC is the transverse compressive strength of the composite.

Fiber Compression Failure
Schultheisz and Waas [26] observed a local matrix deformation accompanied by fiber fracture (i.e., "fiber kinking") that differs from fiber microbuckling. Furthermore, Argon [27] assumed that initial microbuckling could result in fiber rotation, matrix shearing, fiber-matrix debonding, fiber kinking, and splitting in kink bands.
Although the fiber kinking mechanism is similar to that of matrix failure, a major difference is that in the matrix failure model, the associated stresses are calculated with respect to the fracture plane, whereas in fiber kinking analysis, stresses are calculated with respect to the kink plane, which is aligned with the fiber rotation. As shown in Figure 2, system 1-2-3 describes the material coordinates, where is the angle between the kink plane and axis 2. Fiber kinking occurs on the plane with local coordinate system 1 m -2 m -3 m , which is obtained by rotating the coordinate system 1-2 -3 by an angle .  The associated stress transformations in the transformation of the above coordinates are expressed as follows [28]:  According to the concept of matrix compression failure, and taking the effect of matrix transverse tension due to fiber misalignment into account, the fiber kinking/splitting criterion can be expressed as [28]: where σ m 22 , τ m 12 , and τ m 23 are the transverse normal stress and the in-plane and out-of-plane shear stresses, respectively, in the coordinate system 1 m -2 m -3 m , and <·> denotes the Macauley symbol, i.e., <x> = (x + |x|)/2. Pinho et al. [28] determined experimentally that fiber kinking takes place only for compressive stress σ 11 ≤ −X C /2; otherwise, fiber splitting occurs.
Angle ψ is expressed by [32]: The misalignment frame orientation θ is the sum of the initial misalignment angle θ i and the additional shear strain γ due to loading: The shear strain γ in the initial misalignment frame is defined as follows: where f CL is the shear function (i.e., τ = f CL (γ)). For linear shear, Equation (25) can be simplified to:

Damage Propagation Criterion
Damage evolution is accompanied by release of strain energy and degradation of the material properties. The loading/unloading and softening stress-strain curves for a combined elastoplastic damage model are shown in Figure 3, where σ 0 and ε 0 are the initial values of the failure stress and strain in all directions, respectively.  Without considering plastic deformation in the fiber direction, the material exhibits a linear elastic behavior before damage. After damage occurs, the stiffness decreases gradually. Before damage is initiated, although irreversible plastic deformations in both the shear and transverse directions are observed, the material stiffness is not degraded. Therefore, a nonlinear behavior is shown during loading. Unlike in the fiber direction, both plastic deformation and stiffness softening occur beyond the point of damage initiation. Stiffness softening is expressed by an exponential damage parameter as follows [33]: Without considering plastic deformation in the fiber direction, the material exhibits a linear elastic behavior before damage. After damage occurs, the stiffness decreases gradually. Before damage is initiated, although irreversible plastic deformations in both the shear and transverse directions are observed, the material stiffness is not degraded. Therefore, a nonlinear behavior is shown during loading. Unlike in the fiber direction, both plastic deformation and stiffness softening occur beyond the point of damage initiation. Stiffness softening is expressed by an exponential damage parameter as follows [33]: A discrete element is the basic unit of the FE method. Although elements with different sizes obey the same stress-strain relationship, the energy release rates of different elements are unequal and proportional to the element size. To alleviate the dependence of the energy release rate on the element size, an element characteristic length L C is introduced with a critical strain energy release rate of G I,C . The damage energy dissipated per unit volume, g I,C , can be defined with the exponential factor A I as an internal variable [23] as follows: where G I,C contains G kink , G split , G IC , and G IIC , which are identical to the fiber kinking and splitting fracture toughness and modes I and II of the matrix fracture toughness, respectively. g I,C (A I ) can be obtained from the following integration: By substituting Equation (29) into Equation (28), the following equation is derived, which is used to solve A I by numerical iterations: Details on the approach used to determine A I are provided in Chen et al. [23].

Cohesive Model
A modified cohesive model [34] was adopted to simulate interfacial failure, as shown in Figure 4a. In this model, mode II and mode III fractures, which were considered as the combined resulting transverse shear modes, were used with mode I to compose the mixed mode, as shown in Figure 4b. The total mixed-mode relative displacement is expressed as: where δ I and δ II are defined as: Here, δ 1 is the normal opening relative displacement, and δ 2, δ 3 are the resulting transverse shear relative displacements.
A quadratic damage initiation criterion under a multiaxial stress state was used to predict the onset of delamination in the cohesive zone [34]: where σ max I and σ max II denote the interlaminar tensile and shear strengths (see in Figure 4), respectively, and σ I and σ II are the resulting normal interlaminar stress and shear stress of the interface, correspondingly. Assuming linear softening in the interface elements after the delamination onset, the fracture energy under mixed-mode loading with the power law criterion [35] is expressed as: where α ∈[1.0, 2.0] is an empirical parameter derived from mixed-mode tests, and G IC and G IIC are the critical energy release rates for pure modes I (opening) and II (shear), respectively. The superscripts "0" and "f " indicate the initial and final values of the effective displacement δ m , respectively, and the max function represents the irreversible damage.
where max I  and max II  denote the interlaminar tensile and shear strengths (see in Figure 4), respectively, and Ι and ΙΙ are the resulting normal interlaminar stress and shear stress of the interface, correspondingly. Assuming linear softening in the interface elements after the delamination onset, the fracture energy under mixed-mode loading with the power law criterion [35] is expressed as: where α ∈[1.0, 2.0] is an empirical parameter derived from mixed-mode tests, and GIC and GIIC are the critical energy release rates for pure modes I (opening) and II (shear), respectively. The superscripts "0" and "f" indicate the initial and final values of the effective displacement m, respectively, and the max function represents the irreversible damage.

Calculation of the Angle of the Fracture Plane
Knops [36] proposed an algorithm to determine the fracture plane in composite laminates. As shown in Figure 5, the fiber orientation of each ply in the laminate is assumed to be in direction 1.

Calculation of the Angle of the Fracture Plane
Knops [36] proposed an algorithm to determine the fracture plane in composite laminates. As shown in Figure 5, the fiber orientation of each ply in the laminate is assumed to be in direction 1. Subsequently, each ply is divided into discrete equal angular intervals of 1 • from direction 1 (0 • ) to 180 • . The stress components σ N (φ), τ T (φ), and τ L (φ) on the planes oriented at each angle are determined for a designated stress state and substituted into Equation (18) to derive the value of f mat . Then, the angle corresponding to the maximal f mat gives a potential orientation of the fracture plane. If f mat reaches 1.0 with increasing external loading, which indicates matrix failure initiation, then the associated angle of a potential fracture plane is the real orientation φ fp of the fracture plane.
Although the procedure developed by Knops can be used to obtain a precise fracture angle, its efficiency is low. For a model with N elements and M increments, M × N × 180 iteration steps are required. Here, we propose a modified algorithm to determine the fracture angle with high efficiency and robustness based on the studies of Wiegand et al. [37] and Schirmaier et al. [38]. According to the research of Schirmaier et al. [38], the number of f mat thresholds does not exceed three, and the distance between two local maxima is always greater than 25 • . Therefore, the fracture angle can be searched for in [−90 • , 90 • ] with a step size of 10 • , as shown in Figure 6. (If the curve increases or decreases monotonously, which is not shown in Figure 6, the fracture angle is 90 • or −90 • ). The local extrema intervals (denoted as "range 1" and "range 2") are determined by applying the golden section search method and an inverse interpolation.
Subsequently, each ply is divided into discrete equal angular intervals of 1° from direction 1 (0°) to 180°. The stress components N (ϕ), T(ϕ), and L(ϕ) on the planes oriented at each angle are determined for a designated stress state and substituted into Equation (18) to derive the value of fmat. Then, the angle corresponding to the maximal fmat gives a potential orientation of the fracture plane. If fmat reaches 1.0 with increasing external loading, which indicates matrix failure initiation, then the associated angle of a potential fracture plane is the real orientation ϕfp of the fracture plane. Although the procedure developed by Knops can be used to obtain a precise fracture angle, its efficiency is low. For a model with N elements and M increments, M × N × 180 iteration steps are required. Here, we propose a modified algorithm to determine the fracture angle with high efficiency and robustness based on the studies of Wiegand et al. [37] and Schirmaier et al. [38]. According to the research of Schirmaier et al. [38], the number of fmat thresholds does not exceed three, and the distance between two local maxima is always greater than 25°. Therefore, the fracture angle can be searched for in [−90°, 90°] with a step size of 10°, as shown in Figure 6. (If the curve increases or decreases monotonously, which is not shown in Figure 6, the fracture angle is 90° or −90°). The local extrema intervals (denoted as "range 1" and "range 2") are determined by applying the golden section search method and an inverse interpolation.  For range 1 (see the magnified image at the top left corner of Figure 6), i.e., [ϕ1, ϕ4], the golden section points ϕ2 and ϕ3 can be determined as follows:  Although the procedure developed by Knops can be used to obtain a precise fracture angle, its efficiency is low. For a model with N elements and M increments, M × N × 180 iteration steps are required. Here, we propose a modified algorithm to determine the fracture angle with high efficiency and robustness based on the studies of Wiegand et al. [37] and Schirmaier et al. [38]. According to the research of Schirmaier et al. [38], the number of fmat thresholds does not exceed three, and the distance between two local maxima is always greater than 25°. Therefore, the fracture angle can be searched for in [−90°, 90°] with a step size of 10°, as shown in Figure 6. (If the curve increases or decreases monotonously, which is not shown in Figure 6, the fracture angle is 90° or −90°). The local extrema intervals (denoted as "range 1" and "range 2") are determined by applying the golden section search method and an inverse interpolation.  For range 1 (see the magnified image at the top left corner of Figure 6), i.e., [ϕ1, ϕ4], the golden section points ϕ2 and ϕ3 can be determined as follows: For range 1 (see the magnified image at the top left corner of Figure 6), i.e., [φ 1 , φ 4 ], the golden section points φ 2 and φ 3 can be determined as follows: Subsequently, the values of f mat at φ 2 and φ 3 must be compared. If f mat (φ 2 ) < f mat (φ 3 ), the local extremum interval is substituted by [φ 2 , φ 4 ]; otherwise, the local extremum interval is updated by [φ 1 , φ 3 ]. For f mat (φ 2 ) < f mat (φ 3 ) in the current case, let φ a = φ 2 , φ c = φ 4 , and φ b = φ 3 . By applying an interpolation in interval [φ a , φ c ], the function f mat (φ) is approximated as follows: Therefore, the maximum point of range 1, φ fp1 , is expressed as: The maximum point in range 2, denoted as φ fp2 , is obtained by applying a similar procedure on range 2. Finally, the fracture angle in interval [−90 • , 90 • ] is determined as the larger one of φ fp1 and φ fp2 .
To verify the effectiveness of the above algorithm, four typical stress states were selected to determine the angle of the fracture plane: pure shear in-plane and out-of-plane, uniaxial compression, and arbitrary 3D stress, which are listed in Table 1.   A comparison between the proposed algorithm and Knops' algorithm [36] in terms of efficiency and accuracy is shown in Table 2, where N is the number of calculation points and "time" represents the real calculation time. As shown in Table 2, the angles calculated using the proposed algorithm are similar to those derived with Knops' algorithm with step size 0.1°, whereas the associated calculation time is reduced to between 1/10 and 1/15 of the latter. It can be concluded that the algorithm proposed here is highly accurate and efficient.  A comparison between the proposed algorithm and Knops' algorithm [36] in terms of efficiency and accuracy is shown in Table 2, where N is the number of calculation points and "time" represents the real calculation time. As shown in Table 2, the angles calculated using the proposed algorithm are similar to those derived with Knops' algorithm with step size 0.1 • , whereas the associated calculation time is reduced to between 1/10 and 1/15 of the latter. It can be concluded that the algorithm proposed here is highly accurate and efficient.

VUMAT Subroutine
Discontinuity problems do not easily converge in implicit solvers, such as damage problems and nonlinear problems. Hence, we developed a 3D elastoplastic damage algorithm implemented in the user-defined subroutine VUMAT as an extension to that of Chen et al. [23], which is suitable for a plane stress state. To accommodate compressive failure analysis, we replaced the LaRC criterion with the Hashin criterion [39]. In addition, the cohesive model was implemented in Abaqus/Explicit [40] using the user-defined subroutine VUMAT. The flow chart of the subroutine VUMAT is shown in Figure 8, which includes VUMAT-1 and VUMAT-2 for the elastoplastic damage and cohesive zone models, respectively.

Elastoplastic Damage Algorithm
Similar to the study by Chen [23], the user-defined subroutine VUMAT of the elastoplastic damage algorithm was driven by a strain increment, in which the loading history was discretized into a sequence of time intervals, i.e., [tn, tn+1](n = 1, 2, 3). Therefore, the discrete backward Euler algorithm was applied to update the effective stress and strain components, and this process can be described as follows: the Abaqus/Explicit main program yields  

Elastoplastic Damage Algorithm
Similar to the study by Chen [23], the user-defined subroutine VUMAT of the elastoplastic damage algorithm was driven by a strain increment, in which the loading history was discretized into a sequence of time intervals, i.e., [t n , t n+1 ](n = 1, 2, 3). Therefore, the discrete backward Euler algorithm was applied to update the effective stress and strain components, and this process can be described as follows: the Abaqus/Explicit main program yields ∆ε n+1 , ε n , ε e n , ε p n , ε p n , σ n ,σ n , σ n , d 1,n , d 2,n , d 3,n , φ fp,n , a result of the n th increment, used as the initial condition to the n+1 th increment, which is updated to a novel set variable ε n+1 , ε e n+1 , ε p n+1 , ε p n+1 , σ n+1 ,σ n+1 , σ n+1 , d 1,n+1 , d 2,n+1 , d 3,n+1 , φ fp,n+1 at the end of the increment. Therefore, the incremental elastoplastic constitutive algorithm using the backward Euler explicit integration procedure is formulated as follows: The nonlinearity in Equation (40) can be solved by applying the Newton-Raphson method.
Calculate the effective stress on the fracture plane σ N , τ L , and τ T (Equation (17)) and the effective stress in the fiber misalignment frame σ m 22 , τ m 12 , and τ m 23 (Equation (21)). c.
Damage judgment: (1) Check the fiber failure criterion: If σ 11,n+1 ≥ 0 then Calculate f f t (Equation (16)) and d f t Else Calculate f kink or f split Equation (22) and d f c End. (2) Check the matrix failure criterion: If σ N,n+1 ≥ 0 then Calculate f mt (Equation (18)) and d mt Else Calculate f mc (Equation (18)) and d mc End.
(5) Correct the nominal Cauchy stress: a. Update the state dependent variables: b.

Cohesive Zone Algorithm
The numerical algorithm for the cohesive zone implemented in the user-defined subroutine VUMAT is based on [34]. According to the flow chart of VUMAT-2 in Figure 8, the algorithm primarily comprises the calculation of the relative displacement, evaluation of initial and ultimate failure, and derivation of the damage parameter d, which reflects the failure evolution. Therefore, the details of FE implementation in the cohesive zone are as follows: (1) Obtain material properties: Calculate the relative displacement in the local orthogonal coordinate system: b. Calculate the mixed-mode relative displacement δ m and its components: c.

Numerical Examples
To verify the elastoplastic damage model developed in this study, numerical simulations of the progressive failure of composite laminates with fiber waviness were performed. Mukhopadhyay et al. [10] conducted an experimental and numerical study of the compressive failure of IM7/8552 laminates with wrinkles. For comparison, the IM7/8552 laminates with a layup of [45 2 /90 2 /−45 2 /0 2 ] 3s used by Mukhopadhyay et al. [10] were selected in the present simulation. Thus, the numerical and experimental results in [10,41] can be used for comparison with the results obtained from the proposed elastoplastic damage model. The geometry of the laminates and the profile of the embedded wrinkles are described in the following.
The elastic properties [10] and plastic parameters [22] are listed in Table 3, where the plastic parameters a 66 , β and m were converted using Equation (11). Table 3. Material properties of IM7/855 and model parameters [10,22]. First, the profile of the embedded wrinkles was plotted in MATLAB; subsequently, the associated FE model was generated using a Python script run by Abaqus GUI, shown in Figure 9. The profile of the wrinkles can be described with a cosine function as follows [10]: where h w is the nodal coordinate of the through-thickness in the wrinkled configuration; h 0 is the nodal coordinate in a wrinkle-free flat laminate; L is the wavelength of the wrinkle; and δ is the wave amplitude. The value of B is unity on the centerline and decreases linearly with the thickness, with a ratio of 1.0:0.63:0.39:0.0 [42]. First, the profile of the embedded wrinkles was plotted in MATLAB; subsequently, the associated FE model was generated using a Python script run by Abaqus GUI, shown in Figure 9. The profile of the wrinkles can be described with a cosine function as follows [10]: where hw is the nodal coordinate of the through-thickness in the wrinkled configuration; h0 is the nodal coordinate in a wrinkle-free flat laminate; L is the wavelength of the wrinkle; and is the wave amplitude. The value of B is unity on the centerline and decreases linearly with the thickness, with a ratio of 1.0:0.63:0.39:0.0 [42]. A 3D quasi-static progressive failure simulation was implemented in Abaqus/Explicit using the user-defined subroutine VUMAT. To reduce the computational cost, two adjacent plies with the same orientation were modeled using continuum 3D eight-noded reduced integration (C3D8R) elements, with one element through the thickness of the laminate, without considering delamination failure between those two plies. Layers of eight-noded 3D cohesive (COH3D8) interface elements with zero A 3D quasi-static progressive failure simulation was implemented in Abaqus/Explicit using the user-defined subroutine VUMAT. To reduce the computational cost, two adjacent plies with the same orientation were modeled using continuum 3D eight-noded reduced integration (C3D8R) elements, with one element through the thickness of the laminate, without considering delamination failure between those two plies. Layers of eight-noded 3D cohesive (COH3D8) interface elements with zero thickness were inserted between the plies at different orientations to model delamination failure, with the cohesive material properties shown in Table 4. Usually, there are cohesive elements with zero and finite thickness, as shown in Figure 10. The roles of the two types of cohesive element are identical. However, a finite-thickness cohesive element must be meshed using extremely small thickness intervals, which greatly increases the computation time in Abaqus/Explicit. To avoid extremely dense grids in thickness, which reduce the computing efficiency, zero-thickness elements were used in the present work.  Table 4. Material properties of cohesive layers [10].   Table 4. Material properties of cohesive layers [10].   The laminate dimensions were 30 × 30 × 6 mm, with a nominal ply thickness of 0.125 mm, as shown in Figure 11a. A fine mesh with in-plane dimensions of 0.25 × 0.25 mm was applied at the location of the wrinkle and towards the laminate edges, whereas comparatively coarser meshes of dimensions 0.5 × 0.25 mm and 0.75 × 0.25 mm were used elsewhere for computational efficiency. A typical model of the compression specimen comprised approximately 270,000 C3D8R elements and 260,000 COH3D8 elements. The boundary conditions shown in Figure 11b are as follows: fully constrained boundary conditions were applied on the left end (fixed) through reference point 1, and displacement boundary conditions were applied to the load direction by constraining the other directions to reference point 2 on the right end.

Compressive Failure Stress
Three levels of wrinkle severity, with a maximum waviness angle of 5.6 • , 9.9 • , and 11.4 • , were investigated. The predicted compressive stresses versus the displacement curves calculated using three models are illustrated in Figure 12. As shown, both the predicted stiffness and the strength decrease with increasing wrinkle level for the three models. The stiffness predicted by the elastic model is the highest, while that obtained with the model described in [10] is the lowest. The strength calculated with the model of [10] is higher than that derived with the elastic model and the proposed elastoplastic model for defect angle 5.6 • , but an opposite trend is observed for the higher angle of 11.4 • . The model in [10] considers the effects of both shear nonlinearity and residual thermal stress. The above opposite trends might be an effect of the residual thermal stress on the stiffness and strength of laminates with different wrinkle levels. It can be seen that the displacement at final failure predicted by the elastoplastic model is significantly greater than those of the other two models, owing to the plastic effects. Three levels of wrinkle severity, with a maximum waviness angle of 5.6°, 9.9°, and 11.4°, were investigated. The predicted compressive stresses versus the displacement curves calculated using three models are illustrated in Figure 12. As shown, both the predicted stiffness and the strength decrease with increasing wrinkle level for the three models. The stiffness predicted by the elastic model is the highest, while that obtained with the model described in [10] is the lowest. The strength calculated with the model of [10] is higher than that derived with the elastic model and the proposed elastoplastic model for defect angle 5.6°, but an opposite trend is observed for the higher angle of 11.4°. The model in [10] considers the effects of both shear nonlinearity and residual thermal stress. The above opposite trends might be an effect of the residual thermal stress on the stiffness and strength of laminates with different wrinkle levels. It can be seen that the displacement at final failure predicted by the elastoplastic model is significantly greater than those of the other two models, owing to the plastic effects. To demonstrate the effectiveness and predictive accuracy of the elastoplastic model, the compressive strengths predicted by the three models for different wrinkle severities are compared with test data [10,41] in Figure 13. As shown, the discrepancy between the results of the elastoplastic model and the test data for the three levels of wrinkle severity are 1.07%, 8.45%, and 1.97%. The errors of the elastic model are 6.93%, 12.04%, and 5.28%, and the prediction errors of the model in [10] are 10.73%, 3.64%, and −2.79%. The elastoplastic model provided better predictions than the elastic model. Compared with the model proposed in [10], the elastoplastic model achieved more accurate results for wrinkle levels 5.6° and 11.4°. Theoretically, the compressive strength of the laminate with wrinkle level 9.9° should be greater than that of wrinkle level 11.4°, but the opposite was observed in the experiment reported in [10] (as shown in Figure 13). The authors of [10] explained in [41] that this test result was a statistical average value of the compressive strength. They found that the test error for wrinkle level 9.9° was greater than those for 5.6° and 11.4°, owing to the strong fluctuations in the dataset obtained at wrinkle level 9.9°. This may be the reason that the model in [10] seems superior to the elastoplastic model in the compressive strength comparison for wrinkle level 9.9°. To demonstrate the effectiveness and predictive accuracy of the elastoplastic model, the compressive strengths predicted by the three models for different wrinkle severities are compared with test data [10,41] in Figure 13. As shown, the discrepancy between the results of the elastoplastic model and the test data for the three levels of wrinkle severity are 1.07%, 8.45%, and 1.97%. The errors of the elastic model are 6.93%, 12.04%, and 5.28%, and the prediction errors of the model in [10] are 10.73%, 3.64%, and −2.79%. The elastoplastic model provided better predictions than the elastic model. Compared with the model proposed in [10], the elastoplastic model achieved more accurate results for wrinkle levels 5.6 • and 11.4 • . Theoretically, the compressive strength of the laminate with wrinkle level 9.9 • should be greater than that of wrinkle level 11.4 • , but the opposite was observed in the experiment reported in [10] (as shown in Figure 13). The authors of [10] explained in [41] that this test result was a statistical average value of the compressive strength. They found that the test error for wrinkle level 9.9 • was greater than those for 5.6 • and 11.4 • , owing to the strong fluctuations in the dataset obtained at wrinkle level 9.9 • . This may be the reason that the model in [10] seems superior to the elastoplastic model in the compressive strength comparison for wrinkle level 9.9 • .
The predicted curves of compressive strength versus the wrinkle severity by the three models were shown in Figure 14. As can be seen, the computed compressive strengths by the elastoplastic model match better with the measured data [10] than those by the other two models. It should be noted that the discrepancy between the prediction by the elastoplastic model and the elastic model decreases with the increase of wrinkle severity.

Failure Mechanism
To further demonstrate the advantages of considering the elastoplastic effect, the damage evolution is examined in this section. Mukhopadhyay et al. [10] used high-speed imaging to record the failure behavior of IM7/8552 composite laminates with fiber waviness subjected to compressive loading. The predicted curves of compressive strength versus the wrinkle severity by the three models were shown in Figure 14. As can be seen, the computed compressive strengths by the elastoplastic model match better with the measured data [10] than those by the other two models. It should be noted that the discrepancy between the prediction by the elastoplastic model and the elastic model decreases with the increase of wrinkle severity.
The predicted curves of compressive strength versus the wrinkle severity by the three models were shown in Figure 14. As can be seen, the computed compressive strengths by the elastoplastic model match better with the measured data [10] than those by the other two models. It should be noted that the discrepancy between the prediction by the elastoplastic model and the elastic model decreases with the increase of wrinkle severity.

Failure Mechanism
To further demonstrate the advantages of considering the elastoplastic effect, the damage evolution is examined in this section. Mukhopadhyay et al. [10] used high-speed imaging to record the failure behavior of IM7/8552 composite laminates with fiber waviness subjected to compressive loading. Figure 15a,b show the damage evolution at wrinkle level 5.6° predicted by the elastic model and the elastoplastic model, respectively, and the corresponding high-speed camera images of the

Failure Mechanism
To further demonstrate the advantages of considering the elastoplastic effect, the damage evolution is examined in this section. Mukhopadhyay et al. [10] used high-speed imaging to record the failure behavior of IM7/8552 composite laminates with fiber waviness subjected to compressive loading. Figure 15a,b show the damage evolution at wrinkle level 5.6 • predicted by the elastic model and the elastoplastic model, respectively, and the corresponding high-speed camera images of the damage sequence reported in [10] are shown in Figure 15c. It can be seen that the predictions of the elastoplastic model agree better with the test images than those of the elastic model. According to the prediction of the elastoplastic model and the test images, fiber failure first occurred in one of the 0 • plies, followed by delamination in the interface between the 90 • and −45 • plies. Subsequently, further fiber failure and through-thickness delamination occurred when the loading was increased. It is clear that the damage locations predicted by the proposed model match well with those in the test images. damage sequence reported in [10] are shown in Figure 15c. It can be seen that the predictions of the elastoplastic model agree better with the test images than those of the elastic model. According to the prediction of the elastoplastic model and the test images, fiber failure first occurred in one of the 0° plies, followed by delamination in the interface between the 90° and −45° plies. Subsequently, further fiber failure and through-thickness delamination occurred when the loading was increased. It is clear that the damage locations predicted by the proposed model match well with those in the test images. The elastoplastic damage prediction results for laminates of wrinkle levels 11.4° and 9.9° are shown in Figures 16 and 17, respectively. The same damage sequence was observed in these two cases. Inter-ply delamination initiated earlier than fiber failure for the wrinkle severity levels 11.4° and 9.9°, in contrast to the case of wrinkle level 5.6°. This is in accordance with the results of Mukhopadhyay et al. [10]; they reported that fiber compressive failure was the dominant mechanism for lower-severity wrinkles (wrinkle severity lower than 8.7° for the laminates discussed), whereas inter-ply delamination was the driving mode of damage for higher-severity wrinkles (wrinkle severity exceeding 8.7° for the laminates discussed). The elastoplastic damage prediction results for laminates of wrinkle levels 11.4 • and 9.9 • are shown in Figures 16 and 17, respectively. The same damage sequence was observed in these two cases. Inter-ply delamination initiated earlier than fiber failure for the wrinkle severity levels 11.4 • and 9.9 • , in contrast to the case of wrinkle level 5.6 • . This is in accordance with the results of Mukhopadhyay et al. [10]; they reported that fiber compressive failure was the dominant mechanism for lower-severity wrinkles (wrinkle severity lower than 8.7 • for the laminates discussed), whereas inter-ply delamination was the driving mode of damage for higher-severity wrinkles (wrinkle severity exceeding 8.7 • for the laminates discussed). Figure 18 shows the damage location prediction of a laminate with wrinkle level 11.4 • . As shown, similar to the case of wrinkle level 5.6 • , the fiber failure and delamination locations obtained with the elastoplastic model match better with those in the test image than the elastic model. and 9.9°, in contrast to the case of wrinkle level 5.6°. This is in accordance with the results of Mukhopadhyay et al. [10]; they reported that fiber compressive failure was the dominant mechanism for lower-severity wrinkles (wrinkle severity lower than 8.7° for the laminates discussed), whereas inter-ply delamination was the driving mode of damage for higher-severity wrinkles (wrinkle severity exceeding 8.7° for the laminates discussed).   Figure 18 shows the damage location prediction of a laminate with wrinkle level 11.4°. As shown, similar to the case of wrinkle level 5.6°, the fiber failure and delamination locations obtained with the elastoplastic model match better with those in the test image than the elastic model.

Plastic Effects
The curves of equivalent stress versus equivalent strain for the location of the matrix damage initiation are shown in Figure 19, with clusters of the damage evolution sequence shown together. As shown in Figure 19a,b, matrix damage was initiated in phase ①, where the associated equivalent plastic strains were 2.032 × 10 −4 and 1.5 × 10 −4 . This indicates that plasticity occurred before damage initiation. The matrix damage was initiated at the location of the maximum misalignment angle for these two cases, as shown in the damage cloud in phase ①. As loading proceeded, regions of fiber wrinkles underwent deformation owing to the large strain in the matrix and then triggered the next damage mode, i.e., fiber failure (see Figure 19a in phase ②) or delamination (see Figure 19b in phase ③). At the later stage of damage development, the strains experienced within the kind band were large (4.43 × 10 −4 and 3.42 × 10 −4 for wrinkle levels 5.6° and 11.4°, respectively), as shown in phase ⑤. This is similar to the compressive failure mechanism in unidirectional composite laminates [11]. The nonlinear deformation characteristics of the matrix were described well by the proposed elastoplastic damage model, which could not be achieved with the elastic damage model.

Plastic Effects
The curves of equivalent stress versus equivalent strain for the location of the matrix damage initiation are shown in Figure 19, with clusters of the damage evolution sequence shown together. As shown in Figure 19a,b, matrix damage was initiated in phase ①, where the associated equivalent plastic strains were 2.032 × 10 −4 and 1.5 × 10 −4 . This indicates that plasticity occurred before damage initiation. The matrix damage was initiated at the location of the maximum misalignment angle for these two cases, as shown in the damage cloud in phase ①. As loading proceeded, regions of fiber wrinkles underwent deformation owing to the large strain in the matrix and then triggered the next damage mode, i.e., fiber failure (see Figure 19a in phase ②) or delamination (see Figure 19b in phase ③). At the later stage of damage development, the strains experienced within the kind band were large (4.43 × 10 −4 and 3.42 × 10 −4 for wrinkle levels 5.6° and 11.4°, respectively), as shown in phase ⑤. This is similar to the compressive failure mechanism in unidirectional composite laminates [11]. The nonlinear deformation characteristics of the matrix were described well by the proposed elastoplastic damage model, which could not be achieved with the elastic damage model.

Plastic Effects
The curves of equivalent stress versus equivalent strain for the location of the matrix damage initiation are shown in Figure 19, with clusters of the damage evolution sequence shown together. As shown in Figure 19a,b, matrix damage was initiated in phase 1 , where the associated equivalent plastic strains were 2.032 × 10 −4 and 1.5 × 10 −4 . This indicates that plasticity occurred before damage initiation. The matrix damage was initiated at the location of the maximum misalignment angle for these two cases, as shown in the damage cloud in phase 1 . As loading proceeded, regions of fiber wrinkles underwent deformation owing to the large strain in the matrix and then triggered the next damage mode, i.e., fiber failure (see Figure 19a in phase 2 ) or delamination (see Figure 19b in phase Materials 2020, 13, x FOR PEER REVIEW 23 of 27 (a) Wrinkle level 5.6°; numbers ①, ② and ③ indicate the initial damage of the matrix, fiber, and delamination, respectively; numbers ④and ⑤ correspond to the damage evolution of the matrix.
(b) Wrinkle level 11.4°; numbers ①, ③ and ④ indicate the initial damage of the matrix, delamination, and fiber, respectively, whereas numbers ② and ⑤ represent the damage evolution of the matrix. For comparison, the predicted compressive stress versus displacement curves obtained using the elastic model and the elastoplastic model are plotted in Figure 20 for wrinkle levels 5.6° and 11.4°, with clouds of damage development sequences shown together. For the two wrinkle levels, the For comparison, the predicted compressive stress versus displacement curves obtained using the elastic model and the elastoplastic model are plotted in Figure 20 for wrinkle levels 5.6 • and 11.4 • , with clouds of damage development sequences shown together. For the two wrinkle levels, the initiations of matrix damage, fiber failure, and delamination were delayed in the elastoplastic model compared with the elastic model. It is worth noting that the ultimate displacement at final failure for wrinkle level 5.6 • was larger than that for 11.4 • . This indicates that the laminate with wrinkle level 5.6 • is relatively ductile, whereas that with wrinkle level 11.4 • exhibits a brittle behavior. Moreover, the discrepancy between the ultimate displacements obtained from the elastic model and the elastoplastic model for wrinkle level 5.6 • are much greater than those for wrinkle level 11.4 • . Combined with the predicted curves of compressive strength versus wrinkle severity in Figure 14, it may imply that the plastic effect should not be neglected for laminates with lower wrinkle levels, at least for the laminate with the ply sequence discussed herein. initiations of matrix damage, fiber failure, and delamination were delayed in the elastoplastic model compared with the elastic model. It is worth noting that the ultimate displacement at final failure for wrinkle level 5.6° was larger than that for 11.4°. This indicates that the laminate with wrinkle level 5.6° is relatively ductile, whereas that with wrinkle level 11.4° exhibits a brittle behavior. Moreover, the discrepancy between the ultimate displacements obtained from the elastic model and the elastoplastic model for wrinkle level 5.6° are much greater than those for wrinkle level 11.4°. Combined with the predicted curves of compressive strength versus wrinkle severity in Figure 14, it may imply that the plastic effect should not be neglected for laminates with lower wrinkle levels, at least for the laminate with the ply sequence discussed herein.

Conclusions
The progressive failure of multidirectional fiber-reinforced polymer laminates with embedded wrinkle defects was numerically simulated using an elastoplastic damage model. Damage evolution analysis was performed according to the LaRC05 criterion with four damage modes (i.e., fiber tensile failure, matrix failure, fiber kinking/splitting, and delamination). A modified algorithm with high efficiency and robustness was proposed based on a previous study to rapidly search for the fracture angle in the matrix failure analysis. This algorithm achieves high accuracy and significantly improved computational stability by combining the golden section search method and an inverse interpolation. The elastoplastic damage model was applied to simulate the compressive failure behavior of IM7/8552 [45 2 /90 2 /−45 2 /0 2 ] 3s composite laminates with out-of-plane wrinkles. The results showed that the model can reproduce the nonlinearity of the laminate during the evolution of the damage and provide more accurate compressive strength predictions than the elastic model and a previous model [10]. In addition, the proposed model could determine the damage locations during the progressive failure process by comparing with test images. Based on the comparison of the compressive stress-displacement curves predicted by the elastoplastic and elastic models, it can be concluded that plasticity effects should not be neglected for laminates with low wrinkle levels.