Guidelines and Parameter Selection for the

Simulation of Progressive Delamination

Kyongchan Song
Swales Aerospace/ATK Space, Hampton VA 23669
Carlos G. Dvila, Cheryl A. Rose
NASA Langley Research Center, Hampton VA 23669
Abstract: Turons methodology for determining optimal analysis parameters for the simulation of
progressive delamination is reviewed. Recommended procedures for determining analysis
parameters for efficient delamination growth predictions using the Abaqus/Standard cohesive
element and relatively coarse meshes are provided for single and mixed-mode loading. The
Abaqus cohesive element, COH3D8, and a user-defined cohesive element are used to develop
finite element models of the double cantilever beam specimen, the end-notched flexure specimen,
and the mixed-mode bending specimen to simulate progressive delamination growth in Mode I,
Mode II, and mixed-mode fracture, respectively. The predicted responses are compared with their
analytical solutions. The results show that for single-mode fracture, the predicted responses
obtained with the Abaqus cohesive element correlate well with the analytical solutions. For
mixed-mode fracture, it was found that the response predicted using COH3D8 elements depends
on the damage evolution criterion that is used. The energy-based criterion overpredicts the peak
loads and load-deflection response. The results predicted using a tabulated form of the BK
criterion correlate well with the analytical solution and with the results predicted with the userwritten element.
Keywords: Damage, Composites, Progressive Delamination, Virtual Crack Closure Technique
(VCCT), Cohesive Elements, Fracture Mechanics.

1. Introduction
Delaminations in laminated composite structures usually initiate from discontinuities such as
matrix cracks and free edges or from embedded defects due to the manufacturing processes.
Delaminations can reduce the stiffness and strength of a composite structure. In other cases,
delamination can provide stress relief and delay the final failure of the composite structure.
Therefore, it is important to analyze the progressive growth of delamination in order to predict the
performance of a composite structure and to develop reliable and safe designs.
2008 Abaqus Users Conference

Abaqus/Standard provides the Virtual Crack Closure Technique (VCCT) and a cohesive element
for the simulation of progressive delamination growth. The formulation of the cohesive element is
based on the Cohesive Zone Model (CZM) approach for modeling complex fracture mechanisms
at the crack tip (Dugdale, 1960 and Barenblatt, 1962). The cohesive zone approach models an
extended cohesive zone, or process zone, ahead of the crack-tip using traction-separation laws that
relate the opening displacements in the process zone to the resisting tractions. The tractionseparation laws are defined in each fracture mode by an initial elastic stiffness, the peak traction,
or interfacial strength, and the area under the traction-separation law that is equal to the critical
energy release rate.
Fracture simulations using cohesive elements require significant experience for determining the
mesh requirements and the selection of properties that characterize the traction-separation
response. For typical graphite-epoxy materials, the length of the cohesive zone is less than 1 mm.
For an accurate representation of the distribution of tractions ahead of the crack-tip and an
accurate simulation of delamination propagation, at least 3 to 5 cohesive elements are required in
the cohesive zone. These mesh requirements may not be practical for the simulation of
delamination propagation in structural analyses.
Corigliano, 2003, used two-dimensional models of the Double Cantilever Beam specimen (DCB)
and identified the pathological dependence of the solution on mesh refinement and attributed the
errors in the predictions to the difficulties in calculating the interlaminar stresses accurately. Barua
et al, 2007, conducted parametric studies of the effect of interfacial strength and mesh refinement
on the predicted load-deflection response. Barua proposed establishing mesh requirements by
performing mesh calibration studies using a simple DCB specimen and then applying those
calibration results to more complex components made of the same material.
Turon, 2007, proposed a procedure for relaxing the requirement for extremely fine meshes. The
procedure consists of increasing artificially the length of the cohesive zone by reducing the
interfacial strength. A longer cohesive zone allows for the use of larger elements while
maintaining sufficient accuracy in the calculation of the energy release rate.
Turons procedure for using relatively coarse meshes has been demonstrated with a user-defined
cohesive element (UEL) developed by Turon, 2006. In the present paper, Turons procedure is
reviewed and is applied with the Abaqus cohesive element (COH3D8) to provide basic guidelines
regarding optimal choices for mesh density and cohesive zone parameters for efficient single- and
mixed-mode fracture simulations. Predictions for the load-displacement response for single and
mixed-mode fracture specimens obtained using the Abaqus/Standard cohesive element are
compared to the analytical solutions.

2. Selection of Constitutive Parameters and Mesh Requirements

The constitutive response of the cohesive elements used in this paper is defined by bi-linear
traction-separation law as shown in Figure 1, where T represents the interfacial strength, c is the
critical separation, fail is the separation at failure, and the area under the curve, Gc, is the critical
2

2007 Abaqus Users Conference

strain energy release rate. Under mixed-mode fracture, the properties required to define the
bilinear traction-separation law are the three critical fracture energies GIC, GIIC, GIIIC, the penalty
stiffnesses K1,K2,K3, and the interfacial strengths T, S1, and S2. In the present work, the shear
strengths in the two orthogonal directions, S1 and S2 are assumed to be the same and are referred to
as S.
The initial response of the cohesive element is assumed to be linear until a damage initiation
criterion is met. The penalty stiffness, Ki, of the bi-linear traction-separation law is defined as
Ki =

Ti
ci

(1)

where ci is the critical separation for damage initiation. The value of the penalty stiffness must be
high enough to prevent interpenetration of the crack faces and to prevent artificial compliance
from being introduced into the model by the cohesive elements. However, an overly high value
can lead to numerical problems. Several guidelines have been proposed for obtaining the penalty
stiffness of a cohesive element. Turon, et al., 2007, proposed an equation for the penalty stiffness
for Mode I that ensures that the compliance of the bulk material is much larger than the initial
compliance of the cohesive element:
K3

E3
t

(2)

where is a parameter much larger than 1, E3 is the transverse Youngs modulus of the material,
and t is the thickness of an adjacent sub-laminate. Turon recommends using = 50 . This value is
used for the penalty stiffness in all of the analyses presented herein. In the present paper, the
penalty stiffnesses for all modes are assumed to be the same: K1=K2=K3=K.

Figure 1. Bilinear tractionseparation response of the cohesive element.

To predict the propagation of delamination accurately, it is necessary to have a sufficiently fine
finite element discretization of the cohesive zone to ensure a good representation of the
interlaminar stress fields. The required maximum size of the cohesive elements can be determined
by evaluating the size of the cohesive zone and the number Ne of elements that are needed in the
zone.
2008 Abaqus Users Conference

The length of the cohesive zone, lcz, is the distance from the crack tip to the point when the
maximum cohesive traction is attained. Figure 2 describes the length of the cohesive zone. It has
been shown that, for infinite-sized geometries, the length of the cohesive zone, lcz, is a material
property. For single Mode I and Mode II fracture, the cohesive lengths are approximately (Turon,
2007)
lcz_I = ME2

GIC
T

and lcz_II = ME2

GIIC
S2

(3)

where E2 is the transverse Youngs modulus and M is a parameter that depends on the cohesive
zone theory used to derive the expression for the cohesive zone length. In this study, Hillerborgs
model in which M is equal to 1.0 is used (Hillerborg, 1976). Therefore, the length le of the
cohesive elements must satisfy the following conditions for Mode I and Mode II:
le

lcz _ I
Ne

and le

lcz _ II

(4)

Ne

Figure 2. Length of the cohesive zone.

Turon suggests that a minimum of three elements in the cohesive zone is required to accurately
represent the fracture energy in the cohesive zone (Turon, 2007). However, the cohesive zone
length for typical polymeric composites is less than 1 mm, and such small mesh requirements may
render most structural analyses computationally intractable.
When the length of the cohesive zone is negligible compared to the crack length, the propagation
of damage is governed by linear elastic fracture mechanics (LEFM) principles that state that the
energy release rate for creating a new fracture area must be equal to Gc. Therefore, when LEFM is
valid, the effect of the interfacial strengths can be neglected. It can be observed in Equation 3 that
the length of the cohesive zone increases as the inverse of the square of the interfacial strength.
This suggests that the length of the cohesive zone can be artificially lengthened to ensure that it
spans enough elements of a given size. The interfacial normal and shear strengths that are required
can be obtained using Equations 3 and 4, which give

2007 Abaqus Users Conference

Ta =

E2 GIC
and S a =
N e le

E2 GIIC
N e le

(5)

Therefore, the normal and shear strengths that ensure that enough elements span the cohesive zone
are

T = Min T a , T and S = Min S a , S

(6)

Abaqus provides four damage initiation criteria: maximum stress, maximum strain, quadratic
stress, and quadratic strain criterion. The quadratic stress criterion, shown in Equation 7, is used
throughout this paper.
z

xz 2 yz
+ + =1
S S

(7)

where the Macaulay bracket, < >, signifies that the compressive stress does not contribute to
damage initiation.
Once the damage initiation criterion of the cohesive element is satisfied, the stiffness of the
element is degraded. Equation 8 represents the softening response of the cohesive element

i = (1 d ) K i i , i = 1,2,3

(8)

where d is the damage variable, which has the value d = 0 when the interface is undamaged, and
the value d = 1 when the interface is fully fractured.
In the analyses presented herein, we use the energy-based Benzeggagh and Kenane (BK) damage
evolution criterion shown in Equation 9 (Benzeggagh and Kenane, 1996)

G
GC = GIC + (GIIC GIc ) Shear
GT

(9)

where is the BK material parameter, G Shear = G II + G III , and GT = G I + G II + G III .

In the UEL used herein (Turon, 2006), the BK fracture criterion can be expressed in terms of
separations:

2
GC = G IC + (G IIC G Ic )
1 + 2 2 2

(10)

where is the mode mixity ratio defined in terms of separation.

2008 Abaqus Users Conference

shear
shear + normal

(11)

Abaqus provides two methods to implement the BK criterion. In the standard BK option, the
fracture toughness Gc given by Equation 9 is determined in terms of the accumulated energy
release rates GShear and GT. This option is specified using the following command lines
*DAMAGE EVLOUTION, TYPE=ENERGY,
MIXED MODE BEHAVIOR=BK, POWER=

Alternatively, the BK criterion can also be specified with Gc supplied in tabular form. Using
Equation 10, a table of Gc can be specified as a function of the mode mixity ratio. This option is
specified using the following command lines
*DAMAGE EVLOUTION, TYPE=ENERGY,
MIXED MODE BEHAVIOR=Tabular
Gc, 1

where 1 is defined as

1 =

tan 1 ( / (1 ))

(12)

In the following simulations, both options are used and their corresponding results are compared
to the results obtained with a user-defined element and with analytical solutions.

3. Analyses of progressive delamination growth

The Abaqus cohesive element COH3D8 and the user-defined cohesive element COH_UEL
developed by Turon, 2006 were used to develop finite element (FE) models of the double
cantilever beam (DCB) specimen and the end-notched flexure (ENF) specimen to simulate
progressive delamination growth under Mode I and Mode II fracture, respectively. In Sections 3.1
and 3.2, the results of nonlinear analyses performed to investigate the effect of mesh density and
interfacial strength on the predicted load-deflection response of the DCB and ENF specimens are
presented. From the progressive delamination growth predictions of the FE models for singlemode fracture, the mesh density required for solution convergence and the corresponding normal
and shear strengths for a given material are obtained. Progressive delamination growth predictions
of single-mode fracture obtained with the Abaqus cohesive element and COH_UEL correlate
closely with each other. Therefore, only the progressive delamination growth predictions obtained
with the Abaqus cohesive element are presented with the analytical solutions.
FE models of the mixed-mode bending (MMB) specimen were developed using the mesh density
and the interfacial strengths that resulted in converged results for the single-mode analyses. In
Section 3.3, the predicted progressive delamination results for the MMB specimen with the
Abaqus cohesive element and the COH_UEL are compared with analytical solutions and with
Abaqus VCCT predictions.

2007 Abaqus Users Conference

FE models of each specimen are developed using the material properties of AS4/3501-6 that are
provided in Table 1. All specimens are unidirectional and the fibers are aligned with the direction
of fracture propagation.

Table 1. Mechanical and interface material properties of AS4/3501-6.

E1

E2

E3

12

G12

G13

G23

148 GPa

10.50 GPa

10.50 GPa

0.27

5.61 GPa

5.61 GPa

3.17 GPa

GIIC

GIIIC

53.78 MPa

86.88 MPa

1.8

GIC
0.08 kJ/m

3.1

0.55 kJ/m

0.55 kJ/m

Double Cantilever Beam Specimen

FE models of a DCB specimen are developed to investigate the influence of mesh refinement and
the normal interfacial strength on the predicted responses for Mode I simulation. The length of the
specimen, L, is 101.6 mm long, and the width is 7.62 mm. The specimen is composed of two sublaminates of thickness t=1.524 mm. The initial crack length, ao, is 29.21 mm. Each sub-laminate
has a stacking sequence of [03]. The configuration of the specimen is shown in Figure 3.
The DCB specimen is modeled with two layers of 4-node shell elements S4 that are connected
together with cohesive elements. The reference surfaces of the sub-laminates are moved from the
mid-surface of the sub-laminates to the contact surface between the two sub-laminates using the
shell property option OFFSET
Three levels of mesh refinement were considered and the corresponding models are labeled DCB1, -2, and -3. All three models have a uniform mesh in a 10 mm-long damage propagation zone
ahead of the crack tip, which is labeled c0 in Figure 3. The element sizes in zone co of models
DCB-1, DCB-2, and DCB-3 are 1.0 mm, 0.5 mm, and 0.32 mm, respectively.
F

Crack tip

2t

a0

c0
L

Figure 3. Configuration of double cantilever beam (DCB) specimen.

The cohesive zone length for Mode I fracture calculated from Equation 3 is lcz_I = 0.29 mm. The
penalty stiffness determined from Equation 2 with = 50 is K = 344 GPa/mm. The nominal normal
strength provided in Table 1 and the adjusted normal strength obtained using Equation 5 are used
to study the influence of mesh size and normal strength on the prediction of fracture propagation
2008 Abaqus Users Conference

in Mode I. Three different values for the number of elements within the cohesive zone, Ne, are
substituted into Equation 5 to calculate the adjusted normal strength. The number of elements
within the cohesive zone and the adjusted normal strength corresponding to the different mesh
refinements, DCB-1, -2, and -3, are shown in Table 2.
Table 2. Adjusted normal strengths corresponding to the three mesh refinements.
Number of Elements in
Cohesive Zone
Ne =1
Ne =3
Ne =5

DCB-1
le_ModeI = 1mm
29.6 MPa
17.1 MPa
13.2 MPa

DCB-2
le_mode I = 0.5 mm
41.9 MPa
24.2 MPa
18.7 MPa

DCB-3
le_mode I = 0.32mm
53.8 MPa
30.2 MPa
23.4 MPa

The predicted load-deflection responses obtained using model DCB-1 are compared to the
analytical solution (Reeder, 2004) in Figure 4. The model with the nominal normal strength and an
adjusted normal strength with Ne =1 significantly overpredicts the strength of the DCB specimen
and the response does not follow the softening response of the analytical solution. The responses
predicted using an adjusted strength with Ne =3 and 5 compare well with the analytical solution.
However, these models exhibit over-softening before the peak load and they display spurious
oscillations during delamination propagation. Therefore, model DCB-1 is too coarse to predict
accurately the propagation of delamination, even when the strengths are lowered using Equation 5.
The predicted load-deflection responses obtained using model DCB-2 are shown in Figure 5. The
response predicted with the nominal normal strength and the response predicted with an adjusted
normal strength with Ne = 1 significantly overpredict the peak load. Both responses exhibit minor
oscillations during fracture propagation. The responses and the peak loads obtained with Ne = 3
and Ne = 5 correlate well with the analytical solution.
The load-deflection responses obtained using the most refined model, DCB-3, are shown in Figure
6. For this level of mesh refinement, the adjusted strength with Ne =1 is approximately equal to the
nominal strength. The solution obtained with nominal strength overpredicts the peak load and the
softening response of the analytical solution by approximately 7%. The peak loads and the loaddeflection responses predicted using normal strengths adjusted with Ne = 3 and Ne = 5 correlate
well with the analytical solution.
A comparison of the results obtained with models DCB-1, DCB-2, and DCB-3 indicates that,
when used in conjunction with an adjusted strength with Ne = 3, model DCB-2 is the coarsest
model that provides accurate results.

2007 Abaqus Users Conference

20
18
16

Force, F (N.)

14
12
10
8

Analytical Solution

Nominal Strength
Adjusted Strength (Ne=1)

Adjusted Strength (Ne=3)

2

Adjusted Strength (Ne=5)

0
0

0.2

0.4

0.6

0.8

1.2

1.4

Displacement, (mm.)

Figure 4. Effect of the normal interfacial strength on the predicted load-deflection

response for FE model DCB-1.
20
18
16

Force, F (N.)

14
12
10
8

Analytical Solution

Nominal Strength
Adjusted Strength (Ne=1)

Adjusted Strength (Ne=3)

Adjusted Strength (Ne=5)

0
0

0.2

0.4

0.6

0.8

1.2

1.4

Displacement, (mm.)

Figure 5. Effect of the normal interfacial strength on the predicted load-deflection

response for FE model DCB-2.

2008 Abaqus Users Conference

16
14

Force, F (N.)

12
10
8

Analytical Solution

Adjusted Strength (Ne =1)

Adjusted Strength (Ne=3)

Adjusted Strength (Ne=5)

0
0

0.2

0.4

0.6

0.8

1.2

1.4

Displacement, (mm.)

Figure 6. Effect of the normal interfacial strength on the predicted load-deflection

response for FE model DCB-3.

3.2

End Notched Flexure (ENF) Specimen

FE models of the ENF specimen were developed to investigate the influence of mesh refinement
and the interfacial strength on the predicted response for Mode II fracture simulation. The ENF
specimen has the same dimensions and stacking sequence as the DCB specimen analyzed in
Section 3.1. For the ENF specimen, the force, F, is applied at the middle of the specimen. The
ENF specimen configuration and the boundary conditions are shown in Figure 7.
FE models of the ENF specimens with two levels of mesh refinement, indicated as model ENF-1
and ENF-2, are developed to investigate the sensitivity of the predicted responses to the mesh
refinement. These models have a uniform mesh in a damage propagation zone of length co equal to
21.6 mm ahead of the crack tip. The span-wise element length in zone co of Models ENF-1 and
ENF-2 is 1.0 mm and 0.5 mm, respectively. The cohesive zone length for Mode II fracture
calculated from Equation 3 is lcz_ II = 0.77 mm.
F

Crack tip

0.5L

2t
w =0

a0

c0
L

u, w =0

Figure 7. Configuration and boundary conditions of the End Notched Flexure

specimen.
10

2007 Abaqus Users Conference

A comparison of the predicted load-deflection responses of the FE models and the analytical
solution (Mi, 1998) is shown in Figure 8. The results of the FE models for the ENF specimen are
consistent with the analytical solutions and indicate that both models are sufficiently refined and
that the difference in element size between the two models has a minimal influence on the
predicted responses. In addition, when using an element size of 1 mm or smaller, no adjustment to
the nominal shear strength is necessary. These results indicate that the simulation of progressive
delamination of a composite structure subjected to Mode II l does not require a mesh as fine as a
mesh required for Mode I fracture.
200
180
160

Force, F (N.)

140
120
100
80
60

Analytical Solution

40

Nominal Strength (ENF-1, 1mm)

20

Nominal Strength (ENF-2, 0.5mm)

0
0

0.5

1.5

2.5

3.5

Displacement, (mm.)

Figure 8. Effect of mesh refinement on the predicted load-deflection response of

ENF FE models.

3.3

Mixed-Mode Bending (MMB) Specimen

The Abaqus cohesive element and the COH_UEL are used to develop FE models of the MMB
specimen. The specimen is 100 mm-long and 25.5 mm-wide with two 2.3 mm-thick sublaminates. The initial crack length, ao, is 27 mm. The experimental set-up of Reeder, et al.,1990,
shown in Figure 9 is modeled. The level of mode mixity in the MMB is applied by varying the
length of the applied load lever arm, d, of the test apparatus.
In the FE simulation of the MMB specimen, the lever arm length, d, is equal to 52.8 mm and the
corresponding mixed-mode ratio, GII/GT, is equal to 0.2. The penalty stiffness determined from
Equation 2 with =50 and t=2.3 mm is equal to K=228GPa/mm.

2008 Abaqus Users Conference

11

The FE model of the MMB specimen has a uniform mesh in a damage propagation zone ahead of
the crack tip of length co equal to 23.0 mm. The element size in the propagation zone was selected
to be le = 0.5 mm, and the number of elements in the cohesive zone is chosen as Ne=3. These
values were found to provide converged results in the analyses of the DCB and ENF specimens
described in the previous sections. Models of the MMB were created with the Abaqus cohesive
elements and with user-defined elements. Two damage evolution options were evaluated with the
Abaqus cohesive elements: the standard BK and the tabular options described in Section 2. To
improve the convergence rate, a viscoelastic regularization factor of 10-4 was used in the analyses
for both the Abaqus cohesive element and the user-defined element.
Using Equation 5, the reduction factor of the normal strength, T / T , is 45 %. The ENF analyses
indicate that the nominal shear strength of the material does not need to be reduced. However, to
investigate the influence of the shear strength on the mixed-mode fracture predictions, three
different shear strengths are used: the nominal shear strength, S, the shear strength adjusted by
using Equation 5, S , and the shear strength reduced by the same reduction factor as the normal
~
strength, S . The adjusted normal and shear strengths are reported in Table 3.
Table 3. The shear strengths of the cohesive element.
AS4/3501-6

T
S

: Adjusted Normal Strength by Equation 5

: Nominal Shear Strength

24.2 MPa

S
~
S

: Shear Strength Adjusted by Equation 5

63.7 MPa

: Shear Strength Adjusted by Mode I reduction factor

39.2 MPa

86.9 MPa

L
F(d+L)/L
F(d/L)

0.5 L

Crack tip

2t
w =0

a0

c0
L

u, w =0

Figure 9. Configuration, experimental set-up, and boundary conditions of the MMB

specimen.
The predicted load-deflection responses obtained using Abaqus cohesive elements, the VCCT
tool, and the analytical solution are shown in Figure 10. The VCCT response and the analytical
12

2007 Abaqus Users Conference

solution correlate well with each other. The responses predicted with the standard BK option and
the Abaqus COH3D8 elements overpredict the force required for crack propagation by
approximately 35%. However, when the tabulated form of the BK option is used, the predicted
results correlate well with the analytical solution and with the VCCT prediction. The predicted
solutions are also sensitive to the normal-to-shear strengths used in the analyses. For the ratios
used, the three load-deflection curves obtained with the tabulated BK option differ by
approximately 6%.
The results predicted with the standard BK option and with the tabulated BK option should be
identical to each other if the mode-mixity does not change during damage evolution. The fact that
the results shown in Figure 10 do not correlate well with each other indicates that the mode-mixity
changes during damage evolution. This unexpected result will require further investigation of the
evolution of damage in the MMB test.
The predicted load-deflection responses obtained using the user-defined cohesive element (Turon,
2006), COH_UEL, are shown in Figure 11. When the adjusted normal and nominal shear
strengths, T , S , are used, the FE model overpredicts the peak load and the load-deflection response
compared to the analytical solution and the VCCT prediction. The results also indicate that
reducing the shear strength influences the predicted load-deflection response of the FE model.
When the adjusted normal and shear strengths, T , S , obtained by Equation 5 are used for the
interfacial strengths, the predicted load-deflection responses obtained with the COH_UEL
correlate well with the analytical solution and with Abaqus VCCT predictions. Although the
nominal shear strength can be used as shear strength of COH_UEL for single Mode II fracture, the
predicted response of FE models with COH_UEL indicates that both normal and shear strengths
needed to be adjusted using Turons guideline to obtain accurate fracture simulations for the
mixed-mode fracture.
200
180
160

Force, F (N.)

140
120
Analytical Solution

100

VCCT
BK_Standar
d (T and S )
T
&S

80

BK_Standar
d (T and S )
T
&S
~
BK_Standar
d T and S
T
&S

60

BK_Table
T
& S_t (T and S )

40

T
& S_t (T and S )
BK_Table
T
& S_t T and S~
BK_Table

20

0
0

0.2

0.4

0.6

0.8

Displacement, (mm.)

Figure 10. Load-deflection response of MMB specimen with Abaqus COH3D8

cohesive element.

2008 Abaqus Users Conference

13

200
180
160

Force, F (N.)

140
120
100
Analytical Solution

80

VCCT

(
(

60

T & S BK_Table T and S

40

T & S BK_Table T and S

T & S BK_Table T and S

20

)
)

0
0

0.2

0.4

0.6

0.8

Displacement, (mm.)

Figure 11. Load-deflection response of MMB specimen with user-defined COH_UEL

cohesive element.

4. Concluding Remarks
Turons methodology for selecting analysis parameters for the simulation of delamination
propagation using relatively coarse meshes is reviewed and used to determine analysis parameters
for use with the Abaqus/Standard cohesive element. The Abaqus cohesive element, COH3D8, and
a user-defined cohesive element are used to model a DCB specimen, an ENF specimen, and a
MMB specimen to simulate progressive delamination in Mode I, Mode II, and mixed-mode
fracture, respectively. For single-mode fracture, the predicted responses obtained with the Abaqus
cohesive elements correlate well with the analytical solutions when the interfacial normal and
shear strengths are determined using Turons methodology. For mixed-mode fracture, the
predicted response obtained with Abaqus cohesive element depends on the damage evolution
criterion selected. The results using the standard BK overpredict the analytical results by as much
as 35%, but when the BK criterion is supplied in tabulated form, the results correlate well with the
VCCT results and with the analytical solution.
The selection of cohesive element parameters and mesh requirements for fracture propagation in
structural subcomponents should be established by performing single-mode fracture simulations in
simple specimens of the same material. In the present paper, the optimal mesh density and analysis
parameters for the simulation of fracture of a mixed-mode bending test were determined from
parametric simulations of DCB and ENF specimens. This procedure should also be applied for the
simulation of fracture in more complex structures.

14

2007 Abaqus Users Conference

5. References
1. Barenblatt, G., The Mathematical Theory of Equilibrium Cracks in Brittle Fracture,
Advances in Applied Mechanics, Vol. 7, pp. 55-129, 1962.
2. Barua, A., and Bose, K., "On the Optimum Choices of Cohesive-Zone Parameters Describing
Initiation and Propagation of Cracks," in ECCOMAS Thematic Conference on Mechanical
Response of Composites, Ed. P. P. Camanho et al., Porto, Portugal 2007.
3. Benzeggagh, M. L. and Kenane, M., Measurement of Mixed-mode Delamination Fracture
Toughness of Unidirectional Glass/epoxy Composites with Mixed-mode Bending Apparatus,
Compos. Sci. Technol., Vol. 56, No. 4, pp. 439-449, 1996.
4. Corigliano, A., Milne, I., Ritchie, R. O., and Karihaloo, B., "Damage and Fracture Mechanics
Techniques for Composite Structures," Comprehensive Structural Integrity, Pergamon,
Oxford, 2003, pp. 459-539.
5. Dugdale, D. S., Yielding of Steel Sheets Containing Slits, Journal of the Mechanics and
Physics of Solids, Vol. 8, pp. 100-104, 1960.
6. Hillerborg, A., Moder, M., and Petersson, P. E., "Analysis of Crack Formation and Crack
Growth in Concrete by Means of Fracture Mechanics and Finite Elements," Cement and
Concrete Research, Vol. 6, No. 6, pp. 773-781, 1976.
7. Mi, Y., Crisfield M. A., Davis G. A. O., Hellweg, H. B. Progressive Delamination Using
Interface Elements, J Composites, Vol. 32, pp. 1246-72, 1998.
8. Reeder, J.R., Demarco, K., and Whitley, K. S., "The Use of Doubler Reinforcement in
Delamination Toughness Testing," Composites Part A: Applied Science and Manufacturing,
Vol. 35, No. 11, pp. 1337-1344, 2004.
9. Reeder, J. R., and Crews Jr., J. R. Mixed-Mode Bending Method for Delamination Testing,
AIAA Journal, Vol. 28, No. 7, pp. 1270-1276, 1990.
10. Turon, A. Camanho, P. P., and Dvila, C. G., A Damage Model for the Simulation of
Delamination in Advanced Composites under Variable-Mode Loading, Mechanics of
Materials, Vol. 38, No. 11, pp.1072-1089, 2006.
11. Turon, A., Dvila, C. G., Camanho, P. P., and Costa, J., An Engineering Solution for Mesh
Size Effects in the Simulation of Delamination Using Cohesive Zone Models, Engineering
Fracture Mechanics, Vol. 74, No. 10, pp. 1665-1682, 2007.

2008 Abaqus Users Conference

15