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

doi:10.4236/ojcm.2013.34012

Paper Menu >>

Journal Menu >>

Open Journal of Composite Materials, 2013, 3, 113-126

http://dx.doi.org/10.4236/ojcm.2013.34012 Published Online October 2013 (http://www.scirp.org/journal/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

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 three-dimensional 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 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 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 strength-to-

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 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 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 traction-separa-

tion law instead of a conventional stress-strain 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 [3-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 delamina-

tion, however, some disadvantages of large process zone

approach have been noted. These include numerical in-

stability (e.g., elastic snap-back), 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

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

so-called “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 [14-21]. Moёs and Belyts-

chko [22] introduced an analysis framework capable of

considering cohesive cracks and frictional contact be-

tween crack surfaces in two-dimensional (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 traction-separation 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 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 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

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 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 level-set 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

,





, inside

the material domain is described in two parts: the con-

ventional finite element approximation and the XFEM

enriched field representing the discontinuity [12]:









IJ J

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

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

matrix-form equation is obtained:



(3)

3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates:

A Sensitivity Analysis of Modeling Parameters

115

where U denotes a vector containing the nodal parame-

ters including ordinary degrees of freedom “u” and the

enriched degrees of freedom “a”:



,Uua (4)

The stiffness matrix K and the external load vector F

are defined as:

uu ua

ij ij

ij au aa

ij ij

KKK









(5)



iii

FFF (6)

The stiffness components



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:

iYi X

iZ iY

iZ iX

BNN











































(8)









BNH 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

B[22]:

crack

Enriched by Heaviside function

Influence domain of node J

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









, ,

iiii

ii ii

Z X

NH H

NH HNHH

NH HNH H



 























































(10)

3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates:

A Sensitivity Analysis of Modeling Parameters

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:

TMatGeo SS

KK KBDBGMG



  



(11)

,,,

BDGM are the strain gradient matrix, the ma-

terial constitutive matrix, Cartesian gradient matrix and

re-arranged second Piola-Kirchhoff stress matrix, Sij,

using the identity matrix, 33



, defined as: (see Equa-

tions (12)-(16); including the footnotes).











































































(14)





























































































































































(15)

33 33 33

33 33

XX XY XZ

SYYYZ

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

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

117

cohh

coh

= T

u = ū

Ω

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:

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 non-zero thickness cohesive zone model can

be established in the form of tangential stiffness matrix

for 3D problems [25]:



TMat Geo Coh

CohInterfaceCoh

KKK K

BD BGMG

BDB









 









(20)

3.1. Choosing an Interface Traction-Separation

Law

Needleman [26] implemented a potential function to

model cohesive zone behaviour for ductile interfaces as

follows.



00 0

,11expexp

nnt

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 non-con-

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

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



 



  



  







     







 

 

  

iiiiii

NH NH

YYZ

NH NHNH NHNHNH

xyyz z

XXZZXXZZXXZ

 

 



 







 



 

 



     

 

(13)

3D Nonlinear XFEM Simulation of Delamination in Unidirectional Composite Laminates:

A Sensitivity Analysis of Modeling Parameters

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 traction-separation 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,28-30]:





 



elastic part

1softening part

0decohesion part

Pen

Pen f

TKvvv

TDKvvvv

Tvv





 









(22)

T is the interfacial strength of the material and D is the

damage index defined as:

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

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 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

 (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 elastic-snap 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

max

Figure 3. A bilinear traction-separation 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



 (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 traction-separation 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 sub-section. 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 traction-separation 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.

max

 (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

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

traction-separation 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:

max

cz tran

lEt







(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 level-set

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-

tion-separation 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 J-integral

method [35]:

kjijj

JWn



















 (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-

ear-elastic 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 crack-tip coordinates with

the x1-axis parallel to the crack face. Equation (28) may

not be in a well-suited 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]:

kij k

Wq V

















 (29)

where the integral paths





and c

 denote near-field

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

120

tion region (hatched region in Figure 5) which is sur-

rounded by .



 and c

; V



represents the near-

field domain of crack-tip where in LEFM stress singular-

ity 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-inte-

gral is required to be evaluated by the following coordi-

nate transformation:





0llk k



 (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 J-integral corre-

sponds to the rate of change in potential energy per unit

crack extension, namely, the energy release rate G:

1020

cos sin

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 mixed-mode crack propagation where the failure

prior to the complete debonding is evaluated by the fol-

lowing power law [37]:

I IIIII

IC IICIIIC

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 D5528-01. Composite materials

considered in the current study are T300/977-2 carbon

fiber reinforced epoxy and AS4/PEEK carbon fiber-re-

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).

Figure 5. Local crack-tip co-ordinates and the contour Γ

and its interior area, VΓ.

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 com-

posites 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 sub-

routine, UEL (accessible free-of-charge 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 post-processing 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/977-2 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

121

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/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/977-2 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 load-displacement during delamina-

tion. The mesh-free 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/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 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

load-displacement 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 traction-separation

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

122

Figure 7. Load-displacement curve for fine mesh (le = 0.4 mm) simulation with different cohesive zone lengths for T300/977-2

sample.

Figure 8. Load-displacement curve for coarse mesh (le = 1.25 mm) simulation with different cohesive zone length for T300/

977-2 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-

tion-separation 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/977-2 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

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

123

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

Figure 10. The comparison of load-displacement 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/977-2 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 XFEM-cohesive 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

124

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

works.

Figure 12. The comparison between XFEM load-displacement 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 Institute-National 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

125

773-781.

http://dx.doi.org/10.1016/0008-8846(76)90007-7

[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.

1397-1434.

http://dx.doi.org/10.1016/0022-5096(94)90003-5

[3] G. T. Camacho and M. Ortiz, “Computational Modelling

of Impact Damage in Brittle Materials,” International

Journal of Solids and Structures, Vol. 33, No. 20-22,

1996, pp. 2899-2938.

http://dx.doi.org/10.1016/0020-7683(95)00255-3

[4] P. P. Camanho, C. G. Dávila and M. F. De Moura, “Nu-

merical Simulation of Mixed-Mode Progressive Delami-

nation in Composite Materials,” Journal of Composite

Materials, Vol. 37, No. 16, 2003, pp. 1415-1424.

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 Adhesively-Bonded

Joints,” International Journal of Fracture, Vol. 119, No.

1, 2003, pp. 25-46.

[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. 453-463.

http://dx.doi.org/10.1088/0965-0393/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. 11-12, 2004, pp. 2977-2993.

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. 107-137.

http://dx.doi.org/10.1007/s10704-005-4729-6

[9] M. Nishikawa, T. Okabe and N. Takeda, “Numerical Si-

mulation of Interlaminar Damage Propagation in CFRP

Cross-Ply Laminates Under Transverse Loading,” Inter-

national Journal of Solids and Structures, Vol. 44, No. 10,

2007, pp. 3101-3113.

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. 1665-1682.

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. 601-620.

http://dx.doi.org/10.1002/(SICI)1097-0207(19990620)45:

5<601::AID-NME598>3.0.CO;2-S

[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. 131-150.

http://dx.doi.org/10.1002/(SICI)1097-0207(19990910)46:

1<131::AID-NME726>3.0.CO;2-J

[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. 1-4, 1996, pp. 289-314.

http://dx.doi.org/10.1016/S0045-7825(96)01087-0

[14] N. Sukumar, N. Moёs, B. Moran and T. Belytschko, “Ex-

tended Finite Element Method for Three-Dimensional

Crack Modeling,” International Journal for Numerical

Methods in Engineering, Vol. 48, No. 11, 2000, pp. 1549-

1570.

http://dx.doi.org/10.1002/1097-0207(20000820)48:11<15

49::AID-NME955>3.0.CO;2-A

[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. 51-52, 2001, pp.

6825-6846.

http://dx.doi.org/10.1016/S0045-7825(01)00260-2

[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. 1873-1905. 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. 1-15.

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. 760-788. 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,” Thin-Walled Structures, Vol. 44, No. 9,

2006, pp. 1031-1038.

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. 3373-3393.

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. 21-39.

http://dx.doi.org/10.1007/s10704-009-9423-7

[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. 813-833.

http://dx.doi.org/10.1016/S0013-7944(01)00128-X

[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

126

Applied Mechanics and Engineering, Vol. 196, No. 41-44,

2007, pp. 4087-4100.

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. 4367-4378.

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. 384-401.

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. 525-531.

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. 192-200.

[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.

4774-4792.

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. 3866-3880.

http://dx.doi.org/10.1016/j.engfracmech.2008.02.010

[30] E. Barbieri and M. Meo, “A Meshfree Penalty-Based

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. 139-157.

http://dx.doi.org/10.1007/978-94-011-3638-9_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. 2095-2113.

http://dx.doi.org/10.1016/j.biomaterials.2005.09.040

[33] P. P. Camanho and C. G. Dávila, “Mixed-Mode 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

Size-Scale Transition: Dimensional Analysis vs. Renor-

malization Group Theory,” Engineering Fracture Me-

chanics, Vol. 70, No. 14, 2003, pp. 1809-1839.

http://dx.doi.org/10.1016/S0013-7944(03)00126-7

[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. 379-386.

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.

3967-4001.

http://dx.doi.org/10.1016/S0020-7683(03)00176-8

[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.

543-550. 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.