3 D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates : A Sensitivity Analysis of Modeling Parameters

This article presents a three-dimensional extended finite element (XFEM) approach for numerical simulation of delamination in unidirectional composites under fracture mode I. A cohesive zone model in front of the crack tip is used to include interface material nonlinearities. To avoid instability during simulations, a critical cohesive zone length is defined such that user-defined XFEM elements are only activated along the crack tip inside this zone. To demonstrate the accuracy of the new approach, XFEM results are compared to a set of benchmark experimental data from the literature as well as conventional FEM, mesh free, and interface element approaches. To evaluate the effect of modeling parameters, a set of sensitivity analyses have also been performed on the penalty stiffness factor, critical cohesive zone length, and mesh size. It has been discussed how the same model can be used for other fracture modes when both opening and contact mechanisms are active.


Introduction
Today, composite structures are widely used in high tech engineering applications including aeronautical, marine and automotive industries.They have high strength-toweight ratios, good corrosion resistance, and superior fracture toughness.In addition, they can be engineered based on required strength or performance objectives in each given design.Although fiber-reinforced composites have been proven to provide numerous advantages, they are still prone to cracking, interlaminar delamination, fiber breakage and fiber pull-out failure modes.Among these, delamination is known to be the most common mode that often occurs because of a weak bonding between composite layers, an existing crack in the matrix, broken fibres, fatigue or severe impact.
For modeling delamination, numerous investigations have been performed over the past few decades.Hillerborg et al. [1] introduced a combination of the finite element method (FEM) and an analytical solution to simu-late crack growth.The approach is frequently referred to as "fictitious crack modeling" where a traction-separation law instead of a conventional stress-strain relationship is utilized in the crack tip zone to capture degradation of material properties due to the damage.Xu and Needleman [2] applied an energy potential function to implement a cohesive zone model (CZM) concept during the analysis of interface debonding.Further investigations on improving cohesive interface models were performed in [3][4][5][6][7][8][9].Based on these reports, CZM has been proven to be capable of modeling "large process zones"in the present case the composite delamination interface.When utilized in the simulation of progressive delamination, however, some disadvantages of large process zone approach have been noted.These include numerical instability (e.g., elastic snap-back), reduction of stress intensity upon the delamination initiation, and the softening of the original body in the process zone [10].
In other investigations, a relatively new feature of FEM, known as the extended finite element method (XFEM), has been implemented for numerical modeling of discontinuities.The original XFEM approach was introduced by Belytschko and Black [11] and enhanced by Moёs et al. [12].They implemented the concept of the partition of unity method (PUM), which was earlier employed to develop a method for modeling discontinuity in materials [13].In the basic XFEM, a modified Heaviside step function is implemented to model the crack surface by adding extra degrees of freedom to each node of the so-called "enriched elements" [12].Further improvements of XFEM were presented in different applications with plastic/elastic material domains, fluid/solid phases, and static/dynamic loadings [14][15][16][17][18][19][20][21].Moёs and Belytschko [22] introduced an analysis framework capable of considering cohesive cracks and frictional contact between crack surfaces in two-dimensional (2D) problems.Later on, a similar approach was implemented to model cohesive cracks in concrete specimens [23].The application of the 2D model was also extended to composite materials in [24], by means of utilizing XFEM including CZM with a linear traction-separation law to predict delamination.
In the present article, the above cohesive crack modeling approach is applied to 3D domains and used for predicting mode I fracture behaviour of unidirectional laminates.To this end, an ABAQUS user-element subroutine has been developed to model the nonlinear behaviour of composite samples (T300/977-2 carbon fiber reinforced epoxy and AS4/PEEK carbon fiber reinforced polyether ether ketone) under the standard double cantilever beam (DCB) test.Cohesive zone was added to enrich elements in the crack front under a bilinear traction-separation law.In addition, a technique is introduced for simple implementation of the cohesive zone by avoiding material softening due to the application of large process zone.Namely, to decrease the computational time and to avoid instability during simulations, a critical length of cohesive zone in vicinity of the crack tip is defined such that the user-defined XFEM elements are only assigned inside this region.It is shown that the new technique avoids predefining a complete delamination path along the specimen length, and leads to more realistic prediction of experimental data.Finally, a set of sensitivity analyses have been performed to identify effects of different modeling parameters, while comparing the results to conventional FEM, the mesh free method, and the interface element approaches.

Nonlinear Extended Finite Element (XFEM): A Review of Fundamentals
There are several numerical techniques available for analyzing stress and displacement fields in engineering struc-tures, including FEM, finite difference method (FDM) and meshless methods.Among these, FEM has shown popularity in terms of modeling material nonlinearity effects as well as different complex geometries and boundary conditions.As addressed earlier, the FEM capability in modeling discontinuities was first realized by introducing the partition of unity [13] into the approximating functions-later known as the extended finite element [11,22].In some complex crack problems, XFEM has demonstrated more accurate and stable solutions while the conventional finite element results were rough or highly oscillatory [22].
In XFEM, as an advantage, the finite element mesh is generated regardless of the location of discontinuities.Subsequently, search algorithms such as the level-set or fast marching methods can be utilized to identify the location of any discontinuity with respect to the existing mesh, and also to distinguish between different types of required enrichments for affected elements.Finally, additional auxiliary degrees of freedom are added to the conventional FEM approximation functions in the selected nodes around the discontinuity.To review the method mathematically, let us assume a discontinuity (a crack) within an arbitrary finite element mesh, as depicted in Figure 1.
The displacement field of a point x ,   g u x , inside the material domain is described in two parts: the conventional finite element approximation and the XFEM enriched field representing the discontinuity [12]:   x   is the conventional shape function,   x   is the enrichment function, N is the finite element mesh nodes and f N is the number of enriched nodes of the mesh, I u is the classic degrees of freedom at each node and J a are the additional enriched degrees of freedom at the J th node.The displacement approximation in Equation ( 1) can be implemented in numerical solutions of Linear Elastic Fracture Mechanics (LEFM) to predict displacement fields.Considering the total potential energy governing the problem, one can write: where  ,  , u , b f and t f are the stress tensor, strain tensor, displacement vector, body forces and traction forces, respectively. is the traction boundary and  is the integration domain.Discretizing Equation (2) and applying the variational formulation, the following matrix-form equation is obtained: where U denotes a vector containing the nodal parameters including ordinary degrees of freedom " u " and the enriched degrees of freedom " a": The stiffness matrix K and the external load vector F are defined as: The stiffness components   , , in Equation (5) include the classical ( uu ), enriched ( aa ) and coupled ( ua ) arrays of XFEM approximation: where C is the material constitutive relationship arrays and i B is the shape functions derivatives matrix defined for each degree of freedom in 3D problems as: , N i is the conventional FEM shape functions, H is the Heaviside step function value; X, Y and Z are the reference coordinates.For numerical implementation, integration points can be included by the following shifting amendment in a i B [22]: where i  is the numerical integration (e.g., Gauss quadrature) coordinates in the local system.In order to include cohesive properties to XFEM formulations, Khoei et al. [25] introduced an approach based on contact modeling (to be discussed further in Section 3).Their approach also considered a 3D modeling of nonlinear (large deformation) formulation by forming a total tangential matrix based on both material and geometrical stiffness matrices: B D G M are the strain gradient matrix, the material constitutive matrix, Cartesian gradient matrix and re-arranged second Piola-Kirchhoff stress matrix, S ij , using the identity matrix, 3 3 I  , defined as: (see Equations ( 12)-( 16); including the footnotes). .
In Equations ( 12) to ( 16), x, y and z represent the spatial material coordinate system and also the Heaviside function value is interpolated at each integration point.

Cohesive Interface Modeling
Before In composite materials when the crack propagation initiates, a damage zone appears in front of the crack tip and dissipates the high stress intensity expected in LEFM.This damage zone can be interpreted as a cohesive zone (e.g., due to fiber bridging) which complicates the identification of crack tip (Figure 2).In the presence of a cohesive crack, the total potential energy formulation can be rewritten to consider the traction forces over the crack surface.Neglecting body forces, the related governing equation is written as: where t f , t T and v are the external forces on domain boundary, cohesive traction vector on the crack surface and the crack tip opening displacement vector, respectively.Based on constitutive behaviour of the interface material, the traction forces are directly related to the opening and sliding displacements of the crack faces.A can be implemented to rearrange the displacement nodal vector and rewrite the crack opening displacement and cohesive traction as follows [23].
where Interface D represents the interface material properties matrix.
Conventionally, in order to add traction forces into XFEM, a "contact" bond region within enriched elements is assumed and the gauss integration points within this bond are considered for formulation of the contact stiffness matrix.Applying the same idea for the cohesive traction, a non-zero thickness cohesive zone model can be established in the form of tangential stiffness matrix for 3D problems [25]:

Choosing an Interface Traction-Separation Law
Needleman [26] implemented a potential function to model cohesive zone behaviour for ductile interfaces as follows.
  where φ 0 , v n , v t , v n0 , v t0 are the material fracture energy, normal and tangential crack opening displacements, normal and tangential model parameters, respectively.This potential function has been utilized in many research works to extract the mechanical behavior of interface layers; however it has shown some disadvantages such as introducing softening and numerical instability to FEM, especially during crack propagation steps [10].Also, it neglects the material behavior dependency on different fracture modes which can result in a non-conservative estimation of critical forces.To overcome the aforementioned problems, a rigid cohesive model was proposed in [27].In this model, a high initial stiffness is applied to the interface elements and the degradation of material is assessed by damage indices to reduce the penalty stiffness, K Pen , continuously.Despite reliable results acquired by this approach, numerical instabilities were observed when the material degradation was commenced in front of the cracked region and the convergence became dependent of the penalty stiffness.In the present work, a bilinear traction-separation law (Figure 3) was implemented to model the interface material degradation in the opening mode in both normal and tangential directions, which has been also widely used in earlier investigations [4,[28][29][30]: T is the interfacial strength of the material and D is the damage index defined as: Remark: For modeling contact/closing modes, a similar approach is used to consider a diagonal matrix for where a penalty stiffness value is assigned to each diagonal component depending on the normal or tangential (sliding) behaviour of the interface [25].For general-purpose fracture simulations, our model through a user-defined code in ABAQUS recognizes whether the crack faces are in the opening or closing mode and applies the corresponding Interface D for each enriched element.
In the cohesive opening mode all the three diagonal terms of Interface D may be considered identical and the model assumes the element degradation begins when the crack relative displacement exceeds the critical crack tip opening value, v 0 .In turn, v 0 can be defined in terms of the penalty stiffness, K Pen , and the maximum interfacial strength, T max , as follows.
Therefore, selecting an appropriate value of K Pen becomes critical to establish a stable cohesive finite element model.While choosing a large value of K Pen may help true estimation of the elements stiffness before the crack initiation, it will reduce the required critical relative displacement value for the crack initiation and, hence, cause numerical instability upon the crack initiation, known as the elastic-snap back.Using a reasonable value for K Pen can lead to accurate results while attaining a low computational cost.Earlier works have been undertaken to formulate K Pen based on different types of material properties.Turon et al. [10] proposed a simple relationship between the transverse modulus of elasticity, E tran , the specimen thickness, t, and the penalty stiffness, K Pen : where α was proposed to be equal to 50 to prevent the stiffness loss.Nonetheless, in several other cases, a comparison between numerical analysis and experimental results is required to identify an optimum value of penalty stiffness.

Length of the Cohesive Zone
Another important factor in the numerical simulation of delamination is the length of cohesive zone.As opening or sliding displacement increases, elements in the cohesive zone gradually reach the maximum interfacial strength and the maximum stress rises up to the critical interfacial stress ahead of the crack tip.Upon this point, the affected elements' stiffness moves into the softening region of traction-separation law and experiences an irreversible degradation.The maximum length of cohesive zone occurs when the crack tip elements are considered debonded completely (Figure 4), and then the simulation should move to the next stage by expanding the crack length.The criterion for such complete opening is discussed in the subsequent sub-section.Considering a correct value of the length of cohesive zone is essential in numerical modeling of delamination, to prevent difficulties due to implementing the traction-separation law instead of a conventional (continuous) constitutive relationship.Other researchers have studied extensively this topic for general crack problems.Hillerborg et al. [1] proposed a characteristic length parameter for isotropic materials as follows.
2 max where l cz , G C and E are the cohesive zone length, the critical energy release rate and the Young's modulus of the material, respectively.Other equations can be found for estimating the cohesive zone length [10].For various traction-separation laws, Planas and Elices [31] introduced a different equation for isotropic materials.For orthotropic materials, like composite laminates, Yang et al. [32] discussed the possible effects of longitudinal, transverse and shear moduli as well as the laminate thickness, t, on the cohesive zone length.Subsequently, they suggested a modified formulation for measuring the cohesive zone length in slender composite laminates: For FEM implementations, the number of elements attaining a cohesive constitutive behavior is directly related to the cohesive zone length.Accordingly, a range of values for proper mesh size in cohesive zone has been proposed by other researchers [33] and [34], yet it is deemed difficult to estimate an exact value that can be optimum for all fracture simulations.

Crack Initiation and Growth Criteria
In conventional application of cohesive zone models, a predefined crack path is utilized to model the cracking behavior in the structure using a separate layer of cohesive elements.
The crack evolution (opening and propagation) can only occur by failure of these elements under a given loading condition and failure criteria.Damage indices are normally employed to reduce the affected elements' stiffness in subsequent simulation steps.In the present work, however, a separate set of cohesive elements have not been employed; instead, by applying the level-set method, nonlinear XFEM elements are embedded with a cohesive behavior in front of the crack.This in turn reduces the cost of computations and prevents the traction-separation law from imposing unnecessary softening into simulations, especially at the stage of crack evolution.As a threshold criterion, the crack growth is directly related to the energy release rate utilizing the J-integral method [35]: where Г is an "arbitrary" contour surrounding the cracktip with no intersection to other discontinuities, W is the strain energy density defined as W = (1/2)σ ij ε ij for a linear-elastic material, n j is the j th component of the outward unit normal to Г, δ 1j is Kronecker delta, and the coordinates are taken to be the local crack-tip coordinates with the x 1 -axis parallel to the crack face.Equation (28) may not be in a well-suited form for finite element implementations; hence an equivalent form of this equation has been proposed by exploiting the divergence theorem and additional assumptions for homogeneous materials [36]: where the integral paths   and c  denote near-field and crack surface paths, respectively.V  is the integra-tion region (hatched region in Figure 5) which is surrounded by  .  and c  ; V  represents the nearfield domain of crack-tip where in LEFM stress singularity is expected (Figure 5).q is a function varying linearly from 1 (near the crack-tip) to 0 (towards the exterior boundary,  ).To extract the energy release rate from the J-integral, the crack-axis components of the J-integral is required to be evaluated by the following coordinate transformation: where α lk is the coordinate transformation tensor and θ 0 is the crack angle with respect to the global coordinate system.The tangential component of the J-integral corresponds to the rate of change in potential energy per unit crack extension, namely, the energy release rate G: In Mode I and Mode II fracture analyses, the crack propagation occurs when the measured energy release rate exceeds its critical value.This can also be depicted in the mixed-mode crack propagation where the failure prior to the complete debonding is evaluated by the following power law [37]: where G I , G II and G III are the energy release rates for Mode I, Mode II and Mode III; G IC , G IIC and G IIIC are the critical energy release rates for Mode I, Mode II and Mode III which can be found from standard fracture tests; α, β and γ are empirical fracture critical surface parameters fitted using experimental data.
If crack propagation happens in any step of the numerical analysis, the crack tip extends by the length of cohesive zone (e.g., Equation ( 27)) and presents elements debonding in the simulation.It is also expected that the energy release rate per crack extension reaches its critical value when each element in the cohesive zone exceeds the critical opening displacement, v 0 .

Illustrative Example: Numerical Simulation of DCB Test
Double Cantilever Beam (DCB) is a standard test used to evaluate the mode I fracture toughness and failure properties of materials.The DCB samples are normally fabricated based on ASTM D5528-01.Composite materials considered in the current study are T300/977-2 carbon fiber reinforced epoxy and AS4/PEEK carbon fiber-reinforced polyether ether ketone which are used, e.g., in the aerospace industries to manufacture airframe structures with reduced weight (instead of steel components).The T300/977-2 specimens have a 150 mm length, 20 mm width, and 1.98 mm thickness for each arm, with an initial crack length of 54 mm as shown in Figure 6.For AS4/PEEK samples, specimens were 105 mm long, 25.4 mm wide and 1.56 mm thick for each arm, with a 33 mm initial crack.Material properties of each specimen are summarized in Table 1.
Previous numerical works on both types of these composites have been performed using cohesive interface layers via conventional finite elements and the mesh-free methods [4,10,30].In the present study, ABAQUS finite element package was employed to implement the new 3D nonlinear XFEM approach (Sections 2 and 3) along with the CZM properties via a user-defined element subroutine, UEL (accessible free-of-charge for research purposes via contacting the corresponding author).A MAT-LAB code was also developed and linked to the FEM package to undertake the analysis framework by performing post-processing of numerical results and evaluating the stability of the crack propagation as well as updating ABAQUS input files and elements properties for each step of the analysis.Results of the XFEM are also compared to other standard numerical approaches to provide further understanding of the XFEM performance in terms of the prediction accuracy and numerical stability.Next, effects of different important modeling variables such as interface stiffness (the penalty factor) and the cohesive region length are studied.

Effects of Different Modeling Approaches
Turon et al. [10] investigated the effective cohesive zone length for T300/977-2 specimens.They suggested a cohesive zone length of 0.9 mm from numerical simulations on a very fine mesh (with element length, l e , of 0.125 mm).Based on their work, the size of elements in the cohesive zone region should not exceed 0.5 mm and a minimum of two elements is required in this region for acceptable modeling results.In Figure 6, the present XFEM results of the fine mesh simulation with the pen-   alty stiffness of 1 × 10 6 N/mm 3 and the cohesive zone length of 1.5 mm are compared to the results available in the literature by means of different types of numerical approaches [4,10] and [30].
Figure 6 shows that all the models predict a similar trend of the global load-displacement during delamination.The mesh-free method [30] overestimates the stiffness of the material and leads to a higher peak opening force by 5%, while cohesive finite approach [10] underestimates the resisting force by 10% in comparison to the experimental data [4].The XFEM estimates the peak opening force by 3% difference from the experiment, and similar to [4] provides more conservative estimation of the fracture behavior of the DCB samples.

Effects of Mesh Size and Cohesive Zone Length
The DCB test of T300/977-2 specimens was simulated using two different mesh sizes, namely the element lengths of 0.4 mm and 1.25 mm, to demonstrate the effect of coarse and fine meshes on the XFEM results.In addition, in each case to present the influence of cohesive zone length, simulations were rerun (Figures 7 and 8) with different l cz values and with a fixed penalty stiffness of 1 × 10 6 N/mm 3 following [10].Recalling Figure 7, in the fine mesh (l e = 0.4 mm) models, it is observed that using 3 (l cz = 1.5 mm) to 6 (l cz = 2.5 mm) elements within the cohesive zone would lead to an accurate estimation of the experimental data, while increasing this critical value to 8 (l cz = 3.5 mm) elements would introduce an unrealistic global softening behavior to the model which regards the previous investigation in term of mesh sensitivity [38].In the coarse mesh (l e = 1.25 mm) runs (Figure 8), only for the case with 3 (l cz = 3.5 mm) elements, the simulation result became relatively agreeable with the experimental values.It is worth adding that in an earlier work, Harper and Hallet [28] had also obtained accurate load-displacement results using different mesh sizes in interface elements.Namely, for smoother numerical results, they decreased the elements size to prevent dynamic effects of larger elements failure such as the sudden drop of fracture energy release rate.They also introduced a global damping factor of 5% into simulations to dissipate the oscillation caused by the cohesive element debonding and the loss of stiffness in each step of crack propagation.In the present study, the enriched elements in the cohesive zone have the aggregation of stiffness from XFEM approximation and the traction-separation law (Equation ( 20)).Hence, when the complete debonding occurs at each delamination step, the affected elements' stiffness will not completely disappear by elimination of the cohesive zone stiffness, and hence the XFEM approximation would inherently prevent the os-  cillations to a certain degree without adapting a damping ratio to the model.This feature can be especially beneficial regarding computational time in explicit finite element analysis.

Effect of Different Penalty Stiffness Factors
As discussed in Section 3.1, accuracy of the bilinear traction-separation law in modeling the process zone can directly depend on the penalty stiffness value, the optimum value of which may change from one crack problem to another.In this section, in order to evaluate the accuracy of XFEM predictions against different penalty stiffness values, a set of simulations with fine mesh were The XFEM results were less sensitive to the larger order of penalty stiffness values (from 10 3 to 10 5 N/mm 3 ) in comparison to the conventional finite element method results [10].Finally, within the latter recommended Pen K range, two sets of complimentary simulations on AS4/PEEK samples were run to see the effect of interaction between the mesh size and the penalty stiffness.According to results in Figures 11 and 12, the mesh sensitivity decreases using lower values of the penalty stiff-  ness, and vice versa.As AS4/PEEK has a higher critical energy release rate in comparison to T300/977-2 samples (Table 1), a larger cohesive zone region should be expected and, hence, the sensitivity of simulations to the element size is reduced.Conversely, the crack simulation of a material with low fracture toughness would necessitate a smaller cohesive zone length and subsequently, would show more sensitivity to the mesh size.

Conclusions
The present work demonstrated the effectiveness of a combined XFEM-cohesive zone model (CZM) approach in 3D numerical simulation of Mode I fracture (delamination) in fiber reinforced composites in the presence of large deformation effects and interface material nonlinearity.Sets of sensitivity analyses were performed to evaluate the effect of modeling parameters on the variation of numerical predictions.For reliable simulations, a minimum of two elements is required within the cohesive zone region (regardless of critical length value) in front of the crack tip.On the other hand, considering a very long cohesive zone would introduce a global softening to simulations and can lead to the underestimation of the peak opening force.A maximum of six elements with a fine mesh was recommended as the limit within the cohesive zone region for Mode I fracture analysis of the studied unidirectional composites, which was in close agreement with previous reports.It was also observed that reducing the penalty stiffness value in the traction-  separation law improves the convergence of numerical simulations and reduces the mesh size sensitivity; however, using conventional FEM this can again cause a softening problem and reduce the peak opening force prediction.The XFEM approach with embedded CZM is found to be less sensitive to the aforementioned effects, particularly when the penalty stiffness value is chosen arbitrarily within the range of transverse and longitudinal moduli of the composite.
Although the present work relied on deterministic fracture behavior of the material, there is no question that in practice mechanical properties of the same composite can vary from one sample/manufacturing process to another, and hence "stochastic" XFEM modeling of composites is worthwhile [39].

Figure 1 .
Figure 1.The influence domain of node J in an arbitrary finite element mesh.

Figure 2 .
Figure 2. (a) A general body in the state of equilibrium under different tractions and boundary conditions; (b) Example of fiber bridging from an actual tested sample under fracture mode I.

Figure 3 .
Figure 3.A bilinear traction-separation law for modeling the interface material degradation.

Figure 4 .
Figure 4. Formation steps of the cohesive zone in front of crack tip during numerical simulation.

Figure 5 .
Figure 5. Local crack-tip co-ordinates and the contour Γ and its interior area, V Γ .

Figure 6 .
Figure 6.The comparison between the global load-displacement (F-Δ) results using different analysis methods on T300/977-2 sample; For XFEM, fine mesh along with the penalty stiffness of 1 × 106 N/mm 3 and the cohesive zone length of 1.5 mm were used.
performed with a wide range of values varying from 10 2 N/mm 3 to 10 5 N/mm 3 , and with an identical cohesive zone length of l cz = 2.5 mm.Results are shown in Figures 9 and 10 and compared to experimental data [4] for both T300/977-2 and AS4/PEEK specimens, respectively.

Figure 9 .
Figure 9.The comparison of load-displacement curves of T300/977-2 sample using different penalty stiffness values.

Figure 10 .
Figure 10.The comparison of load-displacement curves of AS4/PEEK sample using different penalty stiffness values.

Figure 11 .
Figure 11.The comparison between XFEM load-displacement curves for fine mesh analysis of AS4/PEEK and previous works.

Figure 12 .
Figure 12.The comparison between XFEM load-displacement curves for coarse mesh analysis of AS4/PEEK and previous works.