Open Journal of Composite Materials, 2013, 3, 113126 http://dx.doi.org/10.4236/ojcm.2013.34012 Published Online October 2013 (http://www.scirp.org/journal/ojcm) Copyright © 2013 SciRes. OJCM 113 3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Damoon Motamedi, Abbas S. Milani* School of Engineering, University of British Columbia, Kelowna, Canada. Email: *abbas.milani@ubc.ca Received August 9th, 2013; revised September 9th, 2013; accepted September 20th, 2013 Copyright © 2013 Damoon Motamedi, Abbas S. Milani. This is an open access article distributed under the Creative Commons At tribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is prop erly cited. ABSTRACT This article presents a threedimensional extended finite element (XFEM) approach for numerical simulation of de lamination 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 de fined such that userdefined 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 parame ters, 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. Keywords: Composite Materials; Fracture Properties; Double Cantilever Beam; Extended Finite Element; Cohesive Zone Model 1. Introduction Today, composite structures are widely used in high tech engineering applications including aeronautical, marine and automotive industries. They have high strengthto weight 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 fiberreinforced composites have been proven to provide numerous advantages, they are still prone to cracking, interlaminar delamination, fiber breakage and fiber pullout failure modes. Among these, delamination is known to be the most common mode that often occurs because of a weak bonding be tween 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. Hiller borg et al. [1] introduced a combination of the finite ele ment method (FEM) and an analytical solution to simu late crack growth. The approach is frequently referred to as “fictitious crack modeling” where a tractionsepara tion law instead of a conventional stressstrain relation ship is utilized in the crack tip zone to capture degrada tion 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 investiga tions on improving cohesive interface models were per formed in [39]. 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 delamina tion, however, some disadvantages of large process zone approach have been noted. These include numerical in stability (e.g., elastic snapback), reduction of stress in tensity upon the delamination initiation, and the soften ing 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 *Corresponding author.
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 114 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 em ployed 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 socalled “enriched elements” [12]. Further improve ments of XFEM were presented in different applications with plastic/elastic material domains, fluid/solid phases, and static/dynamic loadings [1421]. Moёs and Belyts chko [22] introduced an analysis framework capable of considering cohesive cracks and frictional contact be tween crack surfaces in twodimensional (2D) problems. Later on, a similar approach was implemented to model cohesive cracks in concrete specimens [23]. The applica tion of the 2D model was also extended to composite materials in [24], by means of utilizing XFEM including CZM with a linear tractionseparation law to predict de lamination. In the present article, the above cohesive crack model ing approach is applied to 3D domains and used for pre dicting mode I fracture behaviour of unidirectional lami nates. To this end, an ABAQUS userelement subroutine has been developed to model the nonlinear behaviour of composite samples (T300/9772 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 tractionseparation law. In addition, a technique is introduced for simple imple mentation of the cohesive zone by avoiding material sof tening due to the application of large process zone. Name ly, to decrease the computational time and to avoid insta bility during simulations, a critical length of cohesive zone in vicinity of the crack tip is defined such that the userdefined XFEM elements are only assigned inside this region. It is shown that the new technique avoids predefining a complete delamination path along the spe cimen length, and leads to more realistic prediction of experimental data. Finally, a set of sensitivity analyses have been performed to identify effects of different mod eling parameters, while comparing the results to conven tional FEM, the mesh free method, and the interface element approaches. 2. Nonlinear Extended Finite Element (XFEM): A Review of Fundamentals There are several numerical techniques available for ana lyzing 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 ca pability in modeling discontinuities was first realized by introducing the partition of unity [13] into the approxi mating 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 levelset or fast marching methods can be utilized to identify the lo cation of any discontinuity with respect to the existing mesh, and also to distinguish between different types of required enrichments for affected elements. Finally, addi tional auxiliary degrees of freedom are added to the con ventional 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 , g ux , inside the material domain is described in two parts: the con ventional finite element approximation and the XFEM enriched field representing the discontinuity [12]: f IJ g IJ J IJ nN nN uxxux xa (1) is the conventional shape function, is the enrichment function, N is the finite element mesh nodes and N is the number of enriched nodes of the mesh, u is the classic degrees of freedom at each node and a are the additional enriched degrees of freedom at the Jth node. The displacement approximation in Equation (1) can be implemented in numerical solutions of Linear Elastic Fracture Mechanics (LEFM) to predict displace ment fields. Considering the total potential energy gov erning the problem, one can write: ddd bt ufu (2) where , , u, b and t are the stress tensor, strain tensor, displacement vector, body forces and trac tion forces, respectively. is the traction boundary and is the integration domain. Discretizing Equation (2) and applying the variational formulation, the following matrixform equation is obtained: UF (3)
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 115 where U denotes a vector containing the nodal parame ters including ordinary degrees of freedom “u” and the enriched degrees of freedom “a”: T ,Uua (4) The stiffness matrix K and the external load vector F are defined as: uu ua ij ij e ij au aa ij ij KK KKK (5) T , ua iii FFF (6) The stiffness components ,, rs ij rs ua in Equa tion (5) include the classical (uu), enriched (aa) and cou pled (ua) arrays of XFEM approximation: Td,, rsr s iji j BCB rsua (7) 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: , , , ,, ,, ,, 00 00 00 0 0 0 iX iY iZ u i iYi X iZ iY iZ iX N N N BNN NN NN (8) , , , ,, ,, ,, 00 00 00 0 0 0 iX iY iZ a i ii YX ii ZY ii ZX NH NH NH BNH NH NH NH NH NH (9) Ni is the conventional FEM shape functions, H is the Heaviside step function value; X, Y and Z are the refer ence coordinates. For numerical implementation, integra tion points can be included by the following shifting amendment in a i B[22]: crack Enriched by Heaviside function Influence domain of node J J Figure 1. The influence domain of node J in an arbitrary finite element mesh. , , , ,, ,, , , 00 00 00 0 0 0 ii X ii Y ii Z a i iiii YX iiii ZY ii ii Z X NH H NH H NH H B NH HNHH NH HNH H NH HNH H (10)
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 116 where i is the numerical integration (e.g., Gauss quad rature) coordinates in the local system. In order to in clude cohesive properties to XFEM formulations, Khoei et al. [25] introduced an approach based on contact mod eling (to be discussed further in Section 3). Their approach also considered a 3D modeling of nonlinear (large deforma tion) formulation by forming a total tangential matrix based on both material and geometrical stiffness matrices: TT dd ep TMatGeo SS KK KBDBGMG (11) ,,, ep SS BDGM are the strain gradient matrix, the ma terial constitutive matrix, Cartesian gradient matrix and rearranged second PiolaKirchhoff stress matrix, Sij, using the identity matrix, 33 , defined as: (see Equa tions (12)(16); including the footnotes). 00 00 00 00 00 00 00 00 00 i i i i ui i i i i i N X N X N N Y N GY N Y N Z N Z N (14) 00 00 00 00 00 00 00 00 00 i i i i ai i i i i i NH X NH X NH X NH Y NH G Y NH Y NH Z NH Z NH Z (15) 33 33 33 33 33 33 . XX XY XZ SYYYZ ZZ SI SI SI MSISI symS I (16) In Equations (12) to (16), x, y and z represent the spa tial material coordinate system and also the Heaviside function value is interpolated at each integration point. 3. Cohesive Interface Modeling Before In composite materials when the crack propaga tion 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 cohe sive zone (e.g., due to fiber bridging) which complicates the identification of crack tip (Figure 2). In the presence iii iii iii u i iiiiii ii ii ii NNN xyz XX XX XX NNN xyz YY YY YY NNN xyz ZZ ZZ ZZ BNNNNNN xyyzz YX XY YX XY YX XY NN NN NN xyyzz YYZZYYZZYYZ iiiiii NNNNNN xyyzz XXZZXXZZXXZ (12)
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 117 Г u f Г cohh f coh = T u = ū Г t Ω Г Fiber bridging Bridged Fibers (a) (b) Figure 2. (a) A general body in the state of equilibrium un der different tractions and boundary conditions; (b) Exam ple of fiber bridging from an actual tested sample under fracture mode I. of a cohesive crack, the total potential energy formula tion can be rewritten to consider the traction forces over the crack surface. Neglecting body forces, the related governing equation is written as: tt dd d coh coh fu Tv (17) where t , 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, respec tively. 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 transformation matrix Coh B can be implemented to rearrange the displacement nodal vector and rewrite the crack opening displacement and cohesive traction as fol lows [23]. Coh vBu (18) nterface Coh TD Bu (19) where nter f ace D represents the interface material proper ties matrix. Conventionally, in order to add traction forces into XFEM, a “contact” bond region within enriched ele ments 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 nonzero thickness cohesive zone model can be established in the form of tangential stiffness matrix for 3D problems [25]: TT T dd d TMat Geo Coh ep SS CohInterfaceCoh KKK K BD BGMG BDB (20) 3.1. Choosing an Interface TractionSeparation Law Needleman [26] implemented a potential function to model cohesive zone behaviour for ductile interfaces as follows. 2 02 00 0 ,11expexp nnt nt nn t vvv vv vv v (21) where φ0, vn, vt, vn0, vt0 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 re search works to extract the mechanical behavior of inter face 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 noncon servative 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 iii iii iii a i iiiiii ii ii NH NH NH xyz XX XX XX NH NH NH xyz YY YY YY NH NH NH xyz ZZ ZZ ZZ BNH NHNH NHNH NH xyyzz YX XY YX XY YX XY NH NHNHNH xx yy ZY YZZY YZ ii iiiiii NH NH zz YYZ NH NHNH NHNHNH xyyz z XXZZXXZZXXZ (13)
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 118 material is assessed by damage indices to reduce the pen alty stiffness, KPen, 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 be came dependent of the penalty stiffness. In the present work, a bilinear tractionseparation law (Figure 3) was implemented to model the interface material degradation in the opening mode in both normal and tangential direc tions, which has been also widely used in earlier investi gations [4,2830]: 0 0 elastic part 1softening part 0decohesion part Pen Pen f f TKvvv TDKvvvv Tvv (22) T is the interfacial strength of the material and D is the damage index defined as: 0 0 f f vvv Dvv v (23) Remark: For modeling contact/closing modes, a simi lar approach is used to consider a diagonal matrix for nter f ace D where a penalty stiffness value is assigned to each diagonal component depending on the normal or tangential (sliding) behaviour of the interface [25]. For generalpurpose fracture simulations, our model through a userdefined code in ABAQUS recognizes whether the crack faces are in the opening or closing mode and ap plies the corresponding nter f ace D for each enriched ele ment. In the cohesive opening mode all the three diagonal terms of nter f ace 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, v0. In turn, v0 can be defined in terms of the penalty stiffness, KPen, and the maximum interfacial strength, Tmax, as follows. max 0 en T vK (24) Therefore, selecting an appropriate value of KPen be comes critical to establish a stable cohesive finite ele ment model. While choosing a large value of KPen may help true estimation of the elements stiffness before the crack initiation, it will reduce the required critical rela tive displacement value for the crack initiation and, hence, cause numerical instability upon the crack initia tion, known as the elasticsnap back. Using a reasonable value for KPen can lead to accurate results while attaining a low computational cost. Earlier works have been un dertaken to formulate KPen based on different types of T max T v v 0 v f K pe Figure 3. A bilinear tractionseparation law for modeling the interface material degradation. material properties. Turon et al. [10] proposed a simple relationship between the transverse modulus of elasticity, Etran, the specimen thickness, t, and the penalty stiffness, KPen: tran Pen E Kt (25) where α was proposed to be equal to 50 to prevent the stiffness loss. Nonetheless, in several other cases, a com parison between numerical analysis and experimental results is required to identify an optimum value of pen alty stiffness. 3.2. 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 cohe sive 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 tractionseparation law and experiences an ir reversible 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 dis cussed in the subsequent subsection. Considering a cor rect value of the length of cohesive zone is essential in numerical modeling of delamination, to prevent difficul ties due to implementing the tractionseparation law in stead of a conventional (continuous) constitutive rela tionship. 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 C cz G lE T (26) where lcz, GC and E are the cohesive zone length, the critical energy release rate and the Young’s modulus of
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 119 Figure 4. Formation steps of the cohesive zone in front of crack tip during numerical simulation. the material, respectively. Other equations can be found for estimating the cohesive zone length [10]. For various tractionseparation laws, Planas and Elices [31] intro duced 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 thick ness, t, on the cohesive zone length. Subsequently, they suggested a modified formulation for measuring the co hesive zone length in slender composite laminates: 14 34 2 max C cz tran G lEt T (27) For FEM implementations, the number of elements at taining a cohesive constitutive behavior is directly related to the cohesive zone length. Accordingly, a range of val ues for proper mesh size in cohesive zone has been pro posed by other researchers [33] and [34], yet it is deemed difficult to estimate an exact value that can be optimum for all fracture simulations. 3.3. 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 cohe sive 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 levelset method, nonlinear XFEM elements are embedded with a cohesive behavior in front of the crack. This in turn re duces the cost of computations and prevents the trac tionseparation law from imposing unnecessary softening into simulations, especially at the stage of crack evolu tion. As a threshold criterion, the crack growth is directly related to the energy release rate utilizing the Jintegral method [35]: 1 1 d j kjijj u JWn x (28) where Г is an “arbitrary” contour surrounding the crack tip with no intersection to other discontinuities, W is the strain energy density defined as W = (1/2)σijεij for a lin earelastic material, nj is the jth component of the outward unit normal to Г, δ1j is Kronecker delta, and the coordi nates are taken to be the local cracktip coordinates with the x1axis parallel to the crack face. Equation (28) may not be in a wellsuited form for finite element implemen tations; hence an equivalent form of this equation has been proposed by exploiting the divergence theorem and additional assumptions for homogeneous materials [36]: ,d j kij k k V u Wq V x (29) where the integral paths and c denote nearfield and crack surface paths, respectively. V is the integra
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 120 tion region (hatched region in Figure 5) which is sur rounded by . and c ; V represents the near field domain of cracktip where in LEFM stress singular ity is expected (Figure 5). q is a function varying linearly from 1 (near the cracktip) to 0 (towards the exterior boundary, ). To extract the energy release rate from the Jintegral, the crackaxis components of the Jinte gral is required to be evaluated by the following coordi nate transformation: 0 0llk k J (30) where αlk is the coordinate transformation tensor and θ0 is the crack angle with respect to the global coordinate sys tem. The tangential component of the Jintegral corre sponds to the rate of change in potential energy per unit crack extension, namely, the energy release rate G: 0 1020 cos sin l GJ JJ (31) 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 mixedmode crack propagation where the failure prior to the complete debonding is evaluated by the fol lowing power law [37]: 1 I IIIII IC IICIIIC GG G GG G (32) where GI, GII and GIII are the energy release rates for Mode I, Mode II and Mode III; GIC, GIIC and GIIIC 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 parame ters fitted using experimental data. If crack propagation happens in any step of the nu merical 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, v0. 4. 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 prop erties of materials. The DCB samples are normally fab ricated based on ASTM D552801. Composite materials considered in the current study are T300/9772 carbon fiber reinforced epoxy and AS4/PEEK carbon fiberre inforced polyether ether ketone which are used, e.g., in the aerospace industries to manufacture airframe struc tures with reduced weight (instead of steel components). V Г x 1 Г ε Г c V ε Θ 0 X 2 X 1 x 2 Figure 5. Local cracktip coordinates and the contour Γ and its interior area, VΓ. The T300/9772 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 com posites have been performed using cohesive interface layers via conventional finite elements and the meshfree 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 userdefined element sub routine, UEL (accessible freeofcharge for research pur poses via contacting the corresponding author). A MAT LAB code was also developed and linked to the FEM package to undertake the analysis framework by per forming postprocessing of numerical results and evalu ating 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 stabil ity. Next, effects of different important modeling vari ables such as interface stiffness (the penalty factor) and the cohesive region length are studied. 4.1. Effects of Different Modeling Approaches Turon et al. [10] investigated the effective cohesive zone length for T300/9772 specimens. They suggested a co hesive zone length of 0.9 mm from numerical simula tions on a very fine mesh (with element length, le, 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
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 121 Figure 6. The comparison between the global loaddisplacement (FΔ) results using different analysis methods on T300/9772 sample; For XFEM, fine mesh along with the penalty stiffness of 1 × 106 N/mm3 and the cohesive zone length of 1.5 mm were used. Table 1. Material properties of the test samples [4,10,30] used in numerical simulations. T300/9772 AS4/PEEK Elastic Properties Fracture Properties Elastic Properties Fracture Properties E11 = 150 GPa Tmax = 45 MPa E11 = 122.7 GPa Tmax = 80 MPa E22 = E33 = 11 GPa GIC = 268 J/m2 E22 = E33 = 10.1 GPa GIC = 969 J/m2 G12 = G13 = 6 GPa G12 = G13 = 5.5 GPa G23 = 3 GPa G23 = 2.2 GPa v12 = v13 = 0.25 v12 = v13 = 0.25 v23 = 0.5 v23 = 0.48 alty stiffness of 1 × 106 N/mm3 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 loaddisplacement during delamina tion. The meshfree method [30] overestimates the stiff ness of the material and leads to a higher peak opening force by 5%, while cohesive finite element 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 experi ment, and similar to [4] provides more conservative es timation of the fracture behavior of the DCB samples. 4.2. Effects of Mesh Size and Cohesive Zone Length The DCB test of T300/9772 specimens was simulated using two different mesh sizes, namely the element lengths of 0.4 mm and 1.25 mm, to demonstrate the ef fect 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 lcz values and with a fixed penalty stiffness of 1 × 106 N/mm3 following [10]. Recalling Figure 7, in the fine mesh (le = 0.4 mm) models, it is observed that using 3 (lcz = 1.5 mm) to 6 (lcz = 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 (lcz = 3.5 mm) elements would introduce an unrealis tic global softening behavior to the model which regards the previous investigation in term of mesh sensitivity [38]. In the coarse mesh (le = 1.25 mm) runs (Figure 8), only for the case with 3 (lcz = 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 loaddisplacement results using different mesh sizes in interface elements. Namely, for smoother numerical re sults, they decreased the elements size to prevent dy namic effects of larger elements failure such as the sud den drop of fracture energy release rate. They also intro duced 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 tractionseparation law (Equation (20)). Hence, when the complete debond ing occurs at each delamination step, the affected ele ments’ stiffness will not completely disappear by elimi nation of the cohesive zone stiffness, and hence the XFEM approximation would inherently prevent the os
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 122 Figure 7. Loaddisplacement curve for fine mesh (le = 0.4 mm) simulation with different cohesive zone lengths for T300/9772 sample. Figure 8. Loaddisplacement curve for coarse mesh (le = 1.25 mm) simulation with different cohesive zone length for T300/ 9772 sample. cillations to a certain degree without adapting a damping ratio to the model. This feature can be especially benefi cial regarding computational time in explicit finite ele ment analysis. 4.3. Effect of Different Penalty Stiffness Factors As discussed in Section 3.1, accuracy of the bilinear trac tionseparation law in modeling the process zone can directly depend on the penalty stiffness value, the opti mum value of which may change from one crack prob lem 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 performed with a wide range of values varying from 102 N/mm3 to 105 N/mm3, and with an identical cohesive zone length of lcz = 2.5 mm. Results are shown in Fig ures 9 and 10 and compared to experimental data [4] for both T300/9772 and AS4/PEEK specimens, respectively. The XFEM results were less sensitive to the larger order of penalty stiffness values (from 103 to 105 N/mm3) in comparison to the conventional finite element method results [10]. Finally, within the latter recommended en range, two sets of complimentary simulations on AS4/PEEK samples were run to see the effect of interac tion between the mesh size and the penalty stiffness. According to results in Figures 11 and 12, the mesh sen sitivity decreases using lower values of the penalty stiff
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 123 Figure 9. The comparison of loaddisplacement curves of T300/9772 sample using different penalty stiffness values. Figure 10. The comparison of loaddisplacement curves of AS4/PEEK sample using different penalty stiffness values. ness, and vice versa. As AS4/PEEK has a higher critical energy release rate in comparison to T300/9772 samples (Table 1), a larger cohesive zone region should be ex pected and, hence, the sensitivity of simulations to the element size is reduced. Conversely, the crack simulation of a material with low fracture toughness would necessi tate a smaller cohesive zone length and subsequently, would show more sensitivity to the mesh size. 5. Conclusions The present work demonstrated the effectiveness of a combined XFEMcohesive zone model (CZM) approach in 3D numerical simulation of Mode I fracture (delami nation) in fiber reinforced composites in the presence of large deformation effects and interface material nonlin earity. Sets of sensitivity analyses were performed to evaluate the effect of modeling parameters on the varia tion 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 co hesive 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
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 124 Figure 11. The comparison between XFEM loaddisplacement curves for fine mesh analysis of AS4/PEEK and previous works. Figure 12. The comparison between XFEM loaddisplacement curves for coarse mesh analysis of AS4/PEEK and previous works. separation law improves the convergence of numerical simulations and reduces the mesh size sensitivity; how ever, using conventional FEM this can again cause a sof tening problem and reduce the peak opening force pre diction. 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 an other, and hence “stochastic” XFEM modeling of com posites is worthwhile [39]. 6. Acknowledgements Financial support from the Natural Sciences and Engi neering Research Council (NSERC) of Canada is greatly acknowledged. Authors are also grateful to their col leagues Drs. M. Bureau, D. Boucher, and F. Thibault from the Industrial Materials InstituteNational Research Council Canada for their useful feedback and discus sions. REFERENCES [1] A. Hillerborg, M. Modeer and P. E. Petersson, “Analysis of Crack Formation and Crack Growth in Concrete by Means of Fracture Mechanics and Finite Elements,” Ce ment and Concrete Research, Vol. 6, No. 3, 1976, pp.
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 125 773781. http://dx.doi.org/10.1016/00088846(76)900077 [2] X. P. Xu and A. Needleman, “Numerical Simulations of Fast Crack Growth in Brittle Solids,” Journal of the Me chanics and Physics of Solids, Vol. 42, No. 9, 1994, pp. 13971434. http://dx.doi.org/10.1016/00225096(94)900035 [3] G. T. Camacho and M. Ortiz, “Computational Modelling of Impact Damage in Brittle Materials,” International Journal of Solids and Structures, Vol. 33, No. 2022, 1996, pp. 28992938. http://dx.doi.org/10.1016/00207683(95)002553 [4] P. P. Camanho, C. G. Dávila and M. F. De Moura, “Nu merical Simulation of MixedMode Progressive Delami nation in Composite Materials,” Journal of Composite Materials, Vol. 37, No. 16, 2003, pp. 14151424. http://dx.doi.org/10.1177/0021998303034505 [5] B. R. K. Blackman, H. Hadavinia, A. J. Kinloch and J. G. Williams, “The Use of a Cohesive Zone Model to Study the Fracture of Fiber Composites and AdhesivelyBonded Joints,” International Journal of Fracture, Vol. 119, No. 1, 2003, pp. 2546. [6] Y. F. Gao and A. F. Bower, “A Simple Technique for Avoiding Convergence Problems in Finite Element Simu lations of Crack Nucleation and Growth on Cohesive In terfaces,” Modelling and Simulation in Materials Science and Engineering, Vol. 12, No. 3, 2004, pp. 453463. http://dx.doi.org/10.1088/09650393/12/3/007 [7] T. M. J. Segurado and C. T. J. Llorca, “A New Three Dimensional Interface Finite Element to Simulate Frac ture in Composites,” International Journal of Solids and Structures, Vol. 41, No. 1112, 2004, pp. 29772993. http://dx.doi.org/10.1016/j.ijsolstr.2004.01.007 [8] Q. Yang and B. Cox, “Cohesive Models for Damage Evo lution in Laminated Composites,” International Journal of Fracture, Vol. 133, No. 2, 2005, pp. 107137. http://dx.doi.org/10.1007/s1070400547296 [9] M. Nishikawa, T. Okabe and N. Takeda, “Numerical Si mulation of Interlaminar Damage Propagation in CFRP CrossPly Laminates Under Transverse Loading,” Inter national Journal of Solids and Structures, Vol. 44, No. 10, 2007, pp. 31013113. http://dx.doi.org/10.1016/j.ijsolstr.2006.09.007 [10] A. Turon, C. G. Dávila, P. P. Camanho and J. Costa, “An Engineering Solution for Mesh Size Effects in the Simu lation of Delamination Using Cohesive Zone Models,” Engineering Fracture Mechanics, Vol. 74, No. 10, 2007, pp. 16651682. http://dx.doi.org/10.1016/j.engfracmech.2006.08.025 [11] T. Belytschko and T. Black, “Elastic Crack Growth in Finite Elements with Minimal Remeshing,” International Journal for Numerical Methods in Engineering, Vol. 45, No. 5, 1999, pp. 601620. http://dx.doi.org/10.1002/(SICI)10970207(19990620)45: 5<601::AIDNME598>3.0.CO;2S [12] N. Moёs, J. Dolbow and T. Belytschko, “A Finite Ele ment Method for Crack Growth without Remeshing,” In ternational Journal for Numerical Methods in Engineer ing, Vol. 46, No. 1, 1999, pp. 131150. http://dx.doi.org/10.1002/(SICI)10970207(19990910)46: 1<131::AIDNME726>3.0.CO;2J [13] J. M. Melenk and I. Babuška, “The Partition of Unity Finite Element Method: Basic Theory and Applications,” Computer Methods in Applied Mechanics and Engineer ing, Vol. 139, No. 14, 1996, pp. 289314. http://dx.doi.org/10.1016/S00457825(96)010870 [14] N. Sukumar, N. Moёs, B. Moran and T. Belytschko, “Ex tended Finite Element Method for ThreeDimensional Crack Modeling,” International Journal for Numerical Methods in Engineering, Vol. 48, No. 11, 2000, pp. 1549 1570. http://dx.doi.org/10.1002/10970207(20000820)48:11<15 49::AIDNME955>3.0.CO;2A [15] J. Dolbow, N. Moёs and T. Belytschko, “An Extended Finite Element Method for Modeling Crack Growth with Frictional Contact,” Computer Methods in Applied Me chanics and Engineering, Vol. 190, No. 5152, 2001, pp. 68256846. http://dx.doi.org/10.1016/S00457825(01)002602 [16] T. Belytschko, H. Chen, J. Xu and G. Zi, “Dynamic Crack Propagation Based on Loss of Hyperbolicity and a New Discontinuous Enrichment,” International Journal for Numerical Methods in Engineering, Vol. 58, No. 12, 2003, pp. 18731905. http://dx.doi.org/10.1002/nme.941 [17] T. Belytschko and H. Chen, “Singular Enrichment Finite Element Method for Elastodynamic Crack Propagation,” International Journal of Computational Methods, Vol. 1, No. 1, 2004, pp. 115. http://dx.doi.org/10.1142/S0219876204000095 [18] P. M. A. Areias and T. Belytschko, “Analysis of Three Dimensional Crack Initiation and Propagation Using the Extended Finite Element Method,” International Journal for Numerical Methods in Engineering, Vol. 63, No. 1, 2005, pp. 760788. http://dx.doi.org/10.1002/nme.1305 [19] A. Asadpoure, S. Mohammadi and A. Vafai, “Crack Ana lysis in Orthotropic Media Using the Extended Finite Element Method,” ThinWalled Structures, Vol. 44, No. 9, 2006, pp. 10311038. http://dx.doi.org/10.1016/j.tws.2006.07.007 [20] D. Motamedi and S. Mohammadi, “Dynamic Analysis of Fixed Cracks in Composites by the Extended Finite Ele ment Method,” Engineering Fracture Mechanics, Vol. 77, No. 17, 2010, pp. 33733393. http://dx.doi.org/10.1016/j.engfracmech.2010.08.011 [21] D. Motamedi and S. Mohammadi, “Dynamic Crack Pro pagation Analysis of Orthotropic Media by the Extended Finite Element Method,” International Journal of Frac ture, Vol. 161, No. 1, 2010, pp. 2139. http://dx.doi.org/10.1007/s1070400994237 [22] N. Moёs and T. Belytschko, “Extended Finite Element Method for Cohesive Crack Growth,” Engineering Frac ture Mechanics, Vol. 69, No. 1, 2002, pp. 813833. http://dx.doi.org/10.1016/S00137944(01)00128X [23] J. F. Unger, S. Eckardt and C. Könke, “Modelling of Co Hesive Crack Growth in Concrete Structures with the Extended Finite Element Method,” Computer Methods in
3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates: A Sensitivity Analysis of Modeling Parameters Copyright © 2013 SciRes. OJCM 126 Applied Mechanics and Engineering, Vol. 196, No. 4144, 2007, pp. 40874100. http://dx.doi.org/10.1016/j.cma.2007.03.023 [24] E. Benvenuti, “A Regularized XFEM Framework for Em bedded Cohesive Interfaces,” Computer Methods in Ap plied Mechanics and Engineering, Vol. 197, No. 1, 2008, pp. 43674378. http://dx.doi.org/10.1016/j.cma.2008.05.012 [25] A. R. Khoei, S. O. R. Biabanaki and M. Anahid, “A La grangian Extended Finite Element Method in Modeling Large Plasticity Deformations and Contact Problems,” International Journal of Mechanical Sciences, Vol. 51, No. 5, 2009, pp. 384401. http://dx.doi.org/10.1016/j.ijmecsci.2009.03.012 [26] A. Needleman, “A Continuum Model for Void Nuclea tion by Inclusion Debonding,” Journal of Applied Me chanics, Vol. 54, No. 3, 1987, pp. 525531. http://dx.doi.org/10.1115/1.3173064 [27] H. Yuan, G. Lin and A. Cornec, “Applications of Cohe sive Zone Model for Assessment of Ductile Fracture Processes,” Journal of Engineering Materials and Tech nology, Vol. 118, No. 1, 1996, pp. 192200. [28] P. Harper and S. R. Hallett, “Cohesive Zone Length in Numerical Simulations of Composite Delamination,” En gineering Fracture Mechanics, Vol. 75, No. 16, 2008, pp. 47744792. http://dx.doi.org/10.1016/j.engfracmech.2008.06.004 [29] C. Fan, J. P. Y. Ben and J. J. R. Cheng, “Cohesive Zone with Continuum Damage Properties for Simulation of Delamination Development in Fiber Composites and Fail ure of Adhesive Joints,” Engineering Fracture Mechanics, Vol. 75, No. 13, 2008, pp. 38663880. http://dx.doi.org/10.1016/j.engfracmech.2008.02.010 [30] E. Barbieri and M. Meo, “A Meshfree PenaltyBased Approach to Delamination in Composites,” Composites Science and Technology, Vol. 69, No. 13, 2009, pp. 2169 2177. http://dx.doi.org/10.1016/j.compscitech.2009.05.015 [31] J. Planas and M. Elices, “Nonlinear Fracture of Cohesive Materials,” International Journal of Fracture, Vol. 51, No. 2, 1991, pp. 139157. http://dx.doi.org/10.1007/9789401136389_10 [32] Q. D. Yang, B. N. Cox, R. K. Nalla and R. O. Ritchie, “Fracture Length Scales in Human Cortical Bone: The Necessity of Nonlinear Fracture Models,” Biomater, Vol. 27, No. 9, 2006, pp. 20952113. http://dx.doi.org/10.1016/j.biomaterials.2005.09.040 [33] P. P. Camanho and C. G. Dávila, “MixedMode Decohe sion Finite Elements for the Simulation of Delamination in Composite Materials,” NASA/TM, No. 211737, 2002. [34] A. Carpinteri, P. Cornetti, F. Barpi and S. Valente, “Co hesive Crack Model Description of Ductile to Brittle SizeScale Transition: Dimensional Analysis vs. Renor malization Group Theory,” Engineering Fracture Me chanics, Vol. 70, No. 14, 2003, pp. 18091839. http://dx.doi.org/10.1016/S00137944(03)001267 [35] J. R. Rice, “A Path Independent Integral and the Ap proximate Analysis of Strain Concentration by Notches and Cracks,” Journal of Applied Mechanics, Vol. 35, No. 1, 1968, pp. 379386. http://dx.doi.org/10.1115/1.3601206 [36] J. H. Kim and G. H. Paulino, “The Interaction Integral for Fracture of Orthotropic Functionally Graded Materials: Evaluation of Stress Intensity Factors,” International Jour nal of Solids and Structures, Vol. 40, No. 15, 2003, pp. 39674001. http://dx.doi.org/10.1016/S00207683(03)001768 [37] E. M. Wu and R. C. Reuter Jr., “Crack Extension in Fi berglass Reinforced Plastics,” University of Illinois, Cham paign, 1965. [38] M. L. Falk, A. Needleman and J. R. Rice, “A Critical Evaluation of Cohesive Zone Models of Dynamic Frac ture,” Journal de Physique IV, Vol. 11, No. 5, 2001, pp. 543550. http://dx.doi.org/10.1051/jp4:2001506 [39] D. Motamedi, A. S. Milani, M. Komeili, M. N. Bureau, F. Thibault and D. Boucher, “A Stochastic XFEM Model to Study Delamination in PPS/Glass UD Composites: Effect of Uncertain Fracture Properties,” Applied Composite Materials, in Press.
