November 2011 

NASA/CR–2011-217305 

NIA Report No. 2011-02 

 

Development and Application of Benchmark 

Examples for Mode II Static Delamination 

Propagation and Fatigue Growth Predictions 

Ronald Krueger 

National Institute of Aerospace, Hampton,Virginia 
NASA STI Program . . . in Profile

 Since its founding, NASA has been dedicated to 

the advancement of aeronautics and space science. 

The NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

 The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. It 

collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI program 

provides access to the NASA Aeronautics and Space 

Database and its public interface, the NASA Technical 

Report Server, thus providing one of the largest 

collections of aeronautical and space science STI in 

the world. Results are published in both non-NASA 

channels and by NASA in the NASA STI Report 

Series, which includes the following report types: 

 TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but having 

less stringent limitations on manuscript length 

and extent of graphic presentations. 

 TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

 CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

 CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA. 

 SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

 TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission. 

 Specialized services also include creating custom 

thesauri, building customized databases, and 

organizing and publishing research results. 

 For more information about the NASA STI 

program, see the following: 

 Access the NASA STI program home page at 

http://www.sti.nasa.gov

 E-mail your question via the Internet to 

help@sti.nasa.gov

 Fax your question to the NASA STI Help Desk 

at 443-757-5803 

 Phone the NASA STI Help Desk at 

443-757-5802 

 Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center Prepared for Langley Research Center 

Hampton, Virginia 23681-2199 under Cooperative Agreement NNL09AA00A 

 

November 2011 

NASA/CR–2011-217305 

NIA Report No. 2011-02 

 

Development and Application of Benchmark 

Examples for Mode II Static Delamination 

Propagation and Fatigue Growth Predictions 

Ronald Krueger 

National Institute of Aerospace, Hampton,Virginia 
Available from:

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076-1320

443-757-5802

Trade names and trademarks are used in this report for identification only. Their usage does not 

constitute an official endorsement, either expressed or implied, by the National Aeronautics and Space 

Administration.
DEVELOPMENT AND APPLICATION OF BENCHMARK EXAMPLES FOR MODE II 

STATIC DELAMINATION PROPAGATION AND FATIGUE GROWTH PREDICTIONS

Ronald Krueger*

ABSTRACT

The development of benchmark examples for static delamination propagation and 

cyclic delamination onset and growth prediction is presented and demonstrated for a 

commercial code. The example is based on a finite element model of an End-Notched 

Flexure (ENF) specimen. The example is independent of the analysis software used and 

allows the assessment of the automated delamination propagation, onset and growth 

prediction capabilities in commercial finite element codes based on the virtual crack closure 

technique (VCCT). First, static benchmark examples were created for the specimen. 

Second, based on the static results, benchmark examples for cyclic delamination growth 

were created. Third, the load-displacement relationship from a propagation analysis and the 

benchmark results were compared, and good agreement could be achieved by selecting the 

appropriate input parameters. Fourth, starting from an initially straight front, the 

delamination was allowed to grow under cyclic loading. The number of cycles to 

delamination onset and the number of cycles during delamination growth for each growth 

increment were obtained from the automated analysis and compared to the benchmark 

examples. Again, good agreement between the results obtained from the growth analysis 

and the benchmark results could be achieved by selecting the appropriate input parameters. 

The benchmarking procedure proved valuable by highlighting the issues associated with 

choosing the input parameters of the particular implementation. Selecting the appropriate 

input parameters, however, was not straightforward and often required an iterative 

procedure. Overall the results are encouraging, but further assessment for mixed-mode 

delamination is required.

1. INTRODUCTION

Over the past two decades, the use of fracture mechanics has become common practice to 

characterize the onset and growth of delaminations. In order to predict delamination onset or 

growth, the calculated strain energy release rate components are compared to interlaminar fracture 

toughness properties measured over a range from pure mode I loading to pure mode II loading. 

The virtual crack closure technique (VCCT) is widely used for computing energy release rates 

based on results from continuum (2D) and solid (3D) finite element (FE) analyses and to supply the 

mode separation required when using the mixed-mode fracture criterion [1, 2]. The virtual crack 

closure technique was recently implemented into several commercial finite element codes. As new 

methods for analyzing composite delamination are incorporated into finite element codes, the need 

for comparison and benchmarking becomes important since each code requires specific input 

parameters unique to its implementation. These parameters are unique to the numerical approach 

chosen and do not reflect real physical differences in delamination behavior.

An approach for assessing the mode I, and mixed-mode I and II, delamination propagation 

capabilities in commercial finite element codes under static loading was recently presented and 

 

*

R. Krueger, National Institute of Aerospace, 100 Exploration Way, Hampton, VA, 23666, resident at Durability, 

Damage Tolerance and Reliability Branch, MS 188E, NASA Langley Research Center, Hampton, VA, 23681, USA.

1
demonstrated for VCCT for ABAQUS®1 [3] as well as MD NastranTM and MarcTM2 [4]. First, 

benchmark results were created manually for full three-dimensional finite element models of the 

Double Cantilever Beam (DCB) and the Single Leg Bending (SLB) specimen. Second, starting 

from an initially straight front, the delamination was allowed to propagate using the automated 

procedure implemented in the finite element software. The approach was then extended to allow the 

assessment of the delamination fatigue growth prediction capabilities in commercial finite element 

codes [5]. As for the static case, benchmark results were created manually first. Second, the 

delamination was allowed to grow under cyclic loading in a finite element model of a commercial 

code. In general, good agreement between the results obtained from the propagation and growth 

analysis and the benchmark results could be achieved by selecting the appropriate input 

parameters. Overall, the results were encouraging but showed that additional assessment for mode 

II and other mixed-mode delamination cases are required.

The objective of the present study was to create additional benchmark examples, independent of 

the analysis software used, which allows the assessment of the static delamination propagation as 

well as the onset and growth prediction capabilities in commercial finite element codes. For the 

simulation of mode II fracture, the three-point bending End Notched Flexure (ENF) specimen was 

selected as shown in Figure 1. Dimensions, layup and material properties were taken directly from a 

related experimental study [6]. 

Static benchmark results were created based on the approach developed earlier [3], using twodimensional finite element models for simulating the ENF specimens with different delamination 

lengths a0. For each delamination length modeled, the load, Q, and the displacement, w, (center 

deflection) were monitored. The mode II strain energy release rate, GII, was calculated for a fixed 

applied load. It is assumed that the delamination propagates when computed energy release rate, GII, 

reaches the fracture toughness GIIc. Thus, critical loads and critical displacements for delamination 

propagation were calculated for each delamination length modeled. From these critical 

load/displacement results, benchmark solutions were created. It is assumed that the 

load/displacement relationship computed during automatic propagation should closely match the 

benchmark cases. 

Benchmark cases to assess the growth prediction capabilities were created based on the finite 

element models of the ENF specimen used for the static benchmark case. First, the number of cycles 

to delamination onset, ND, was calculated from the mode II fatigue delamination growth onset data 

of the material [6]. Second, the number of cycles during delamination growth, ΔNG, was obtained 

incrementally from the material data for mode II fatigue delamination propagation [6] by using 

growth increments of Δa=0.1 mm. Third, the total number of growth cycles, NG, was calculated by 

summing over the increments ΔNG. Fourth, the corresponding delamination length, a, was 

calculated by summing over the growth increments Δa. Finally, for the benchmark case where 

results for delamination onset and growth were combined, the delamination length, a, was 

calculated and plotted versus an increasing total number of load cycles NT=ND+NG. 

After creating the benchmark cases, the approach was demonstrated for the commercial finite 

element code ABAQUS®. Starting from an initially straight front, the delamination was allowed to 

propagate under static loading or grow under cyclic loading based on the algorithms implemented 

into the software. Input control parameters were varied to study the effect on the computed 

delamination propagation and growth. The benchmark enabled the selection of the appropriate 

 

1

 ABAQUS® is a product of Dassault Systèmes Simulia Corp. (DSS), Providence, RI, USA 2

 MD NastranTM and MarcTM are manufactured by MSC.Software Corp., Santa Ana, CA, USA. NASTRAN® is a 

registered trademark of NASA.

2
input parameters that yielded good agreement between the results obtained from the growth 

analysis and the benchmark results. Once the parameters have been identified, they may then be 

used with confidence to model delamination growth for more complex configurations.

In this paper, the development of the benchmark cases for the assessment of the static 

delamination propagation as well as the onset and growth prediction capabilities are presented. 

Examples of automated propagation and growth analyses are shown, and the selection of the 

required code specific input parameters are discussed.

2. METHODOLOGY

For the current numerical investigation, the three-point End-Notched Flexure (ENF) specimen, 

as shown in Figure 1, was chosen since it is simple, only exhibits the mode II opening fracture mode 

and had been used previously to develop an approach to assess the quasi-static delamination 

propagation simulation capabilities in commercial finite element codes [4]. The methodology for 

delamination propagation, onset and growth was applied to the ENF specimen to create the 

benchmark example [7, 8]. For the current study, an ENF specimen made of IM7/8552 

graphite/epoxy with an unidirectional layup, [0]24, was modeled. The material, layup, overall 

specimen dimensions including initial crack length, a, were identical to specimens used in a related 

experimental study [6]. The material properties are given in Tables I and II.

2.1 Static fracture toughness

A quasi-static mixed-mode fracture criterion is determined first, since the parameters are 

required input to the ABAQUS® VCCT solution as discussed in the appendix. The mixed-mode 

fracture criterion for a material is determined by plotting the interlaminar fracture toughness, Gc, 

versus the mixed-mode ratio, GII/GT as shown in Figure 2. The fracture criterion is generated 

experimentally using pure Mode I (GII/GT=0) Double Cantilever Beam (DCB) tests, pure Mode II 

(GII/GT=1) End-Notched Flexure (ENF) test (as shown in Figure 1), and Mixed Mode Bending 

(MMB) tests of varying ratios of GI and GII. For the material used in this study, the experimental 

data (open symbols) and mean values (filled symbols) are shown in Figure 2. The different symbols 

(circles and squares) denote different data sources for the material [6, 9]. A 2D fracture criterion 

was suggested by Benzeggah and Kenane [10] using a simple mathematical relationship between Gc

and GII/GT

€

Gc = GIc + ( ) GIIc − GIc ⋅ GII

GT













η

. (1)

In this expression, GIc and GIIc are the experimentally-determined fracture toughness data for mode I 

and II as shown in Figure 2. The factor 

€

η was determined by a curve fit using the LevenbergMarquardt algorithm in the KaleidaGraphTM graphing and data analysis software [11]. 

2.2 Fatigue delamination growth onset

The number of cycles to delamination onset, ND, can be obtained from the delamination onset 

curve plotted in Figure 3 [6, 12]. The onset curve (solid green line) is a power law fit 

€

G = m0 ⋅ ND

m1 (2)

3
of the experimental data (open, green circles) obtained from an ENF test using a draft standard for 

delamination growth onset [6].

2.3 Fatigue delamination growth

The number of cycles during delamination growth, NG, can be obtained from the fatigue 

delamination propagation relationship (Paris Law) plotted in Figure 4 [6, 12]. The delamination 

growth rate (solid black line) can be expressed as a power law function

€

da

dN = c ⋅Gmax

n (3)

where da/dN is the increase in delamination length per cycle and Gmax is the maximum energy 

release rate at the front at peak loading. The factor c and exponent n were obtained by fitting the 

curve to the experimental data obtained during the ENF tests [6]. The critical energy release rate or 

fracture toughness, GIIc, was included in the plot of Figure 4 (blue solid vertical line). Since 

composites do not exhibit the same threshold behavior commonly observed in metals, a cutoff 

value, Gth, was chosen below which delamination growth was assumed to stop (green solid vertical 

line) [6].

3. FINITE ELEMENT MODELING

3.1 Model description

An example of a two-dimensional finite element model of an End-Notched Flexure (ENF) 

specimen with boundary conditions is shown in Figure 5. The specimen was modeled with solid 

plane strain elements (CPE4, CPE4I) and solid plane stress elements (CPS4) in ABAQUS®

Standard 6.9EF and 6.10. Along the length, all models were divided into different sections with 

different mesh refinement as shown in Figure 5a. The ENF specimen was modeled with six 

elements through the specimen thickness (2h) as shown in the detail of Figure 5b. The resulting 

element length at the delamination tip was Δa=1.0 mm. Finer meshes, resulting in Δa=0.5 mm and 

Δa=0.25 mm were also generated as shown in Figures 5c and d, respectively.

The plane of delamination was modeled as a discrete discontinuity in the center of the specimen. 

For the analysis with ABAQUS® 6.9EF and 6.10, the models were created as separate meshes for 

the upper and lower part of the specimens with identical nodal point coordinates in the plane of 

delamination [13]. Two surfaces (top and bottom surface) were defined to identify the contact area 

in the plane of delamination as shown in Figure 5b. Additionally, a node set was created to define 

the intact (bonded nodes) region. 

Examples of three-dimensional finite element models of the ENF specimen are shown in 

Figures 6 and 7. Along the length, all models were divided into different sections with different 

mesh refinement. A refined mesh was used in the center of the ENF specimen as shown in the detail 

of Figure 6b. Across the width, a uniform mesh (25 elements) was used to avoid potential problems 

at the transition between a coarse and finer mesh [3-5]. Through the specimen thickness (2h), six 

elements were used as shown in the detail of Figure 6b. The resulting element length at the 

delamination tip was Δa=1.0 mm. The specimen was modeled with solid brick elements (C3D8I), 

which had yielded excellent results in a previous studies [3.4]. A refined mesh with an increased 

number of elements in the length direction, resulting in Δa=0.5 mm, was also generated (not 

4
shown). For an additional model, the width direction was also refined (50 elements) resulting in the 

model shown in Figure 7a.

An example of a continuum shell element (SC8R) finite element model of the ENF specimen is 

shown in Figure 7b. Typically, the continuum shell elements in ABAQUS® are used to model an 

entire three-dimensional body. Unlike conventional shells, which model a reference surface, the 

SC8R elements have displacement degrees of freedom only, use linear interpolation, and allow 

finite membrane deformation and large rotations and, therefore, are suitable for geometric nonlinear 

analysis. The continuum shell elements are based on first-order layer-wise composite theory and 

include the effects of transverse shear deformation and thickness change [13]. In the x-y plane, the 

models had the same fidelity as the models made of solid brick elements C3D8I shown in Figures 6 

and 7a, resulting in an element length at the delamination tip Δa=1.0 mm, and Δa=0.5 mm, 

respectively. In the z-direction, only one element was used to model the thickness of each specimen 

arm. These less-refined models were used to study the effect on performance (CPU time), computed 

load/displacement behavior and growth prediction in comparison with the more refined solid brick 

(C3D8I) models shown in Figures 6 and 7a.

3.2 Static delamination propagation analysis

For the automated delamination propagation analysis, the VCCT implementation in ABAQUS®

Standard 6.9EF and 6.10 were used. The plane of delamination in three-dimensional analyses is 

modeled using the existing ABAQUS®/Standard crack propagation capability based on the contact 

pair capability [13]. Additional element definitions are not required, and the underlying finite 

element mesh and model does not have to be modified [13]. The implementation offers a crack and 

delamination propagation capability in ABAQUS®. It is implied that the energy release rate at the 

crack tip is calculated at the end of a converged increment. Once the energy release rate exceeds the 

critical strain energy release rate (including the user-specified mixed-mode criteria as shown in 

Figure 2), the node at the crack tip is released in the following increment, which allows the crack to 

propagate. To avoid sudden loss of stability when the crack tip is propagated, the force at the crack 

tip before advance is released gradually during succeeding increments in such a way that the force is 

brought to zero no later than the time at which the next node along the crack path begins to open 

[13]. 

In addition to the mixed-mode fracture criterion, VCCT for ABAQUS® requires additional input 

for the propagation analysis. If a user specified release tolerance is exceeded in an increment 

€

(G −Gc )/Gc > release tolerance, a cutback operation is performed which reduces the time 

increment. In the new smaller increment, the strain energy release rates are recalculated and 

compared to the user specified release tolerance. The cutback reduces the degree of overshoot and 

improves the accuracy of the local solution [13]. A release tolerance of 0.2 is suggested in the 

handbook [13]. To help overcome convergence issues during the propagation analysis, ABAQUS®

provides:

• contact&stabilization,"which"is"applied"across"only"selected"contact"pairs"and"used"to"

control" the"motion" of" two" contact" pairs"while" they"approach"each" other"in"multi5

body"contact."The"damping"is"applied"when"bonded"contact"pairs"debond"and"move"

away"from"each"other"[13]

• automatic&or&static&stabilization"which"is"applied"to"the"motion"of"the"entire"model"

and"is"commonly"used"in"models" that"exhibit" statically"unstable"behavior" such"as"

buckling"[13]

• viscous®ularization,"which"is"applied"only"to"nodes"on"contact"pairs"that"have"just"

5
debonded. The"viscous"regularization"damping"causes"the"tangent"stiffness"matrix"

of" the" softening" material" to" be" positive" for" sufficiently" small" time" increments."

Viscous" regularization" damping" in" VCCT" for" ABAQUS®" is" similar" to" the" viscous"

regularization"damping"provided" for"cohesive"elements"and" the"concrete"material"

model"in"ABAQUS®/Standard"[13]. Further details about the required input parameters 

are discussed in the appendix.

For automated propagation analysis, it was assumed that the computed behavior should closely 

match the benchmark results created below. For all analyses, the elastic constants (given in Table I) 

and the input to define the fracture criterion (given in Table II) were kept constant. The following 

items were varied to study the effect on the automated delamination propagation behavior during 

the analysis:

• The release tolerance (reltol) was varied. 

• The input for contact stabilization (cs) was varied. 

• The input for viscous®ularization (damv) was varied. 

• Two- and three-dimensional models with different types of elements were used.

• Models with different element length at the crack tip/delamination front, Δa, were used.

3.3 Fatigue delamination onset and growth analysis

For the automated delamination onset and growth analysis, the low-cycle fatigue analysis in 

ABAQUS® Standard 6.9EF and 6.10 was used to model delamination growth at the interfaces in 

laminated composites [13, 14]. A direct cyclic approach is part of the implementation and provides 

a computationally effective modeling technique to obtain the stabilized response of a structure 

subjected to constant amplitude cyclic loading. The theory and algorithm to obtain a stabilized 

response using the direct cyclic approach are described in detail in reference 14. Delamination onset 

and growth predictions are based on the calculation of the strain energy release rate at the 

delamination front using VCCT. To determine propagation, computed energy release rates are 

compared to the input data for onset and growth from experiments as discussed in the methodology 

section. The implementation is set up to release at least one element at the interface after the loading 

cycle is stabilized [13, 14].

For automated delamination onset and growth analyses, it was assumed that the computed 

behavior should closely match the benchmark results created below. For all analyses, the elastic 

constants (given in Table I), the input to define the fracture criterion (given in Table II), and the 

parameters for delamination onset and delamination growth (Paris Law) (given in Table II) were 

kept constant. The parameters to define the load frequency (f=5 Hz), the load ratio (R=0.1) as well 

as the minimum and maximum applied displacement (wmin and wmax) were also kept constant 

during all analyses. A Fourier series is used during the execution of the ABAQUS® Standard to 

approximate the periodic cyclic loading. Based on previous results, it was decided to use 50 

Fourier terms to approximate the periodic cyclic loading [5]. Note, that number of Fourier terms 

is required input in addition and independent of the definition of the applied cyclic loading. Further 

details about the required input parameters are discussed in the appendix where a sample input file 

is also provided. The following items were varied to study the effect on the predicted delamination 

onset and growth behavior during the analysis:

• The input to define the cyclic loading was studied. The starting time, t0, which causes a 

phase-shift was varied. 

• Two- and three-dimensional models with different types of elements were used.

6
• For three-dimensional models, the number of terms used to define a Fourier series was 

reduced to decrease computation time. 

• Models with different element length at the crack tip/delamination front, Δa, were used.

4. DEVELOPMENT OF THE STATIC BENCHMARK CASE

The static benchmark case was created based on the approach developed earlier [3]. Twodimensional finite element models simulating ENF specimens with 15 different delamination 

lengths a0 were created (25.4 mm≤ a0≤76.2 mm). For each delamination length modeled, the load,

Q, and displacement, w, were monitored as shown in Figure 8 (colored lines). Using VCCT, the 

mixed-mode strain energy release rate components were computed for the applied load Q=2000 N 

as shown in Figure 8. As expected, the results were predominantly mode II. Therefore, a failure 

index GII/GIIc was calculated by correlating the results with the mode II fracture toughness, GIIc, of 

the graphite/epoxy material. It is assumed that the delamination propagates when the failure index 

reaches unity. Therefore, the critical load, Qcrit, can be calculated based on the relationship between 

load, Q, and the energy release rate, G [15], 

€

G = Q2

2 ⋅

∂CP

∂A . (4)

In equation (4), CP is the compliance of the specimen, and ∂A is the increase in surface area 

corresponding to an incremental increase in load or displacement at fracture. The critical load, Qcrit, 

and critical displacement, wcrit, were calculated for each delamination length modeled 

€

GII

GIIc

= Q2

Qcrit

2 ⇒ Qcrit = Q GIIc

GII

 , wcrit = w

GIIc

GII

 (5)

and the results were included in the load/displacement plots as shown in Figure 9 (solid red circles). 

These critical load/displacement results indicated that, with increasing delamination length, less 

load is required to extend the delamination. For the first ten delamination lengths, a0, investigated, 

the values of the critical displacements also decreased at the same time. This means that the ENF 

specimen exhibits unstable delamination propagation under load as well as displacement control in 

this region. The remaining critical load/displacement results pointed to stable propagation. 

From these critical load/displacement results (solid red line), two benchmark solutions can be 

created as shown in Figure 10. During the analysis, either prescribed displacements, w, or nodal 

point loads, Q, are applied. For the case of prescribed displacements, w, (dashed blue line), the 

applied displacement must be held constant over several increments once the critical point (Qcrit, 

vcrit) is reached, and the delamination front is advanced during these increments. Once the critical 

path (solid black line) is reached, the applied displacement is increased again incrementally. For the 

case of applied nodal point loads (dashed red line), the applied load must be held constant while the 

delamination front is advanced during these increments. Once the critical path (solid black line) is 

reached, the applied load is increased again incrementally. 

The benchmark result for prescribed displacements may also be visualized by plotting the 

prescribed displacements, w, at delamination growth onset versus the increase in delamination 

length, a*, as illustrated in Figure 11. For the case of applied nodal point loads, the benchmark 

7
result may be visualized by plotting the applied loads, Q, at delamination growth onset versus the 

increase in delamination length, a*, as illustrated in Figure 12. It is assumed that the 

load/displacement, load/delamination-length or applied displacement/delamination-length 

relationship computed during automatic propagation should closely match the benchmark cases.

5. CREATING A BENCHMARK EXAMPLE FOR GROWTH PREDICTION

For the development of the benchmark case for delamination onset from an initial flaw and 

subsequent growth under cyclic loading, guidance was taken from test results for mode II 

interlaminar fracture toughness and fatigue characterization [6]. In the report, a series of tests are 

documented which were conducted under load control with a maximum load, Qmax, which caused 

the energy release rate at the front, GIImax, to reach values of 50% of GIIc

€

GIImax

GIIc

= 0.5 (6)

The maximum load, Q50,max, and the maximum displacement, w50,max, were calculated using the 

known quadratic relationship between energy release rate and applied load or displacement

€

GIImax

GIIc

= Qmax

2

Qcrit

2 ⇒ Qmax = Qcrit

GIImax

GIIc

 , wmax = wcrit

GIImax

GIIc

 (7)

€

Qmax = Qcrit 0.5 , wmax = wcrit 0.5 (8)

where Qcrit and wcrit are the critical values. For the current study, a critical energy release rate 

GIIc=0.78 kJ/m2

 was used, and the critical values, Qcrit , (solid red horizontal line) and, wcrit, (solid 

red vertical line) were obtained from the benchmark for static delamination propagation (discussed 

above) shown in the load-displacement plot in Figure 13. The calculated maximum load, Q50,max, 

and calculated maximum displacement, w50,max, are shown in Figure 13 (dashed red line) in 

relationship to the static benchmark case (solid grey circles and dashed grey line) mentioned above. 

During constant amplitude cyclic loading of an ENF specimen under load control, the applied 

maximum load, Q50,max=1079 N, is kept constant while the displacement increases with increasing 

delamination length (horizontal dashed red line). For simulations performed under displacement 

control, the applied maximum displacement, w50,max=1.0 mm, is kept constant while the load 

decreases as the delamination length increases (vertical dashed red line). In the report on mode II 

interlaminar fracture toughness and fatigue characterization [6], tests are also documented that were 

performed at 40%, 30% and 20% of GIIc. The corresponding maximum loads and the maximum 

displacements were calculated using equation (8) and were included in the plot of Figure 13 

(Q40,max, Q30,max, Q20,max - horizontal lines; w40,max, w30,max, w20,max – vertical lines).

For each of the 15 finite element models representing 15 delamination lengths (25.4 mm ≤

a0≤76.2 mm) as shown in Figure 8, the energy release rate corresponding to an applied load, 

Q50,max=1079 N, was calculated using equation (7). The energy release rate first increased with an 

increase in delamination length, a*, as shown in Figure 14 (open red circles and dashed red line). 

After reaching a peak, the energy release rate decreased with increasing delamination length. 

Delamination growth was assumed to become unstable once the calculated energy release rate 

passed the fracture toughness value GIIc (solid blue line in Figure 14). Later, delamination growth 

8
was assumed to become stable again after the calculated energy release rate dropped below the 

fracture toughness value GIIc. Additionally, the energy release rate dependence on the crack length 

was calculated for Q40,max, Q30,max, and Q20,max, and the results were included in the plot of Figure 14 

(dashed lines with open symbols). The static benchmark case (solid grey circles and dashed grey 

line in Figure 14), where the delamination propagates at constant GIIc (solid blue line in Figure 14) 

was included for comparison. Also included was the cutoff value, Gth, (green solid horizontal line).

For each of the 15 finite element models representing 15 delamination lengths (25.4 mm ≤

a0≤76.2 mm) as shown in Figure 8, the energy release rate corresponding to an applied 

displacement w50,max=1.0 mm was also calculated using equation (7). The energy release rate first 

increased with an increase in delamination length, a*, as shown in Figure 15 (solid red circles and

dashed red line). After reaching a peak, the energy release rate decreased with increasing 

delamination length. Delamination growth was assumed to stop once the calculated energy release 

rate dropped below the cutoff value, Gth, (green solid horizontal line). Additionally, the energy 

release rate dependence on the crack length was calculated for w40,max, w30,max, and w20,max, and the 

results were included in the plot of Figure 15 (dashed lines with solid symbols). The static 

benchmark case (solid grey circles and dashed grey line in Figure 15), where the delamination 

propagates at constant GIIc (solid blue line in Figure 15) was included for comparison.

The ENF fatigue characterization tests were performed at a load ratio R=0.1 [6]. The 

corresponding minimum load, Qmin, and minimum displacement, wmin,, were calculated

€

R = Qmin

Qmax

= wmin

wmax

= 0.1 ⇒ Qmin = 0.1⋅Qmax and wmin = 0.1⋅wmax (9)

Further, the tests were performed at a frequency f=5 Hz [6]. A graphical representation of the cyclic 

fatigue loading is plotted in Figure 16 for the example of displacement control w50,max=1.0 mm. The 

applied displacement, w, is represented as a function of time, t

€

w = A0 + B1 ⋅sinω t − t [ ( )0 ]⋅wmax (10)

where wmax is the maximum displacement. The constants A0=0.55, B1=0.45, the circular frequency 

ω=10π=31.42 and the starting time t0=0.05 are calculated from load ratio R=0.1 and the frequency 

f=5 Hz for testing. The resulting equation to calculate the applied displacement, w, is shown in 

Figure 16.

5.1 Fatigue delamination growth onset

The number of cycles to delamination onset, ND, may be obtained by solving equation (2) for 

ND. 

 

€

G = m0 ⋅ ND

m1 ⇒ ND = 1

m0













1

m1

c1



⋅G

1

m1

 ⇒ ND = c1 ⋅Gc2 (11)

 

€

where c2 = 1

m1

. Values for the constants m0 and m1 as well as c1 and c2 are shown in Table II.

9
At the beginning of the test, the specimen is loaded initially so that the energy release rate at the 

front, GII, reaches about 50% of GIIc corresponding to G50=0.392 kJ/m2

. The initial energy release 

rate is shown in the delamination onset plot of Figure 17 as a horizontal dashed red line. From the 

delamination onset curve, the number of cycles to delamination onset is determined as ND50=80, 

shown as a vertical dashed red line. The values corresponding to 40%, 30% and 20% of GIIc were 

also calculated using equation (11) and were included in the plot of Figure 17 (G40, G30, G20 - 

horizontal dashed lines; ND40, ND30, ND20 – vertical dashed lines). 

 5.2 Fatigue delamination growth 

The number of cycles during delamination growth can be obtained by solving equation (3) for 

NG

€

NG = ∫ dN = 1

c

Gmax

-n ∫ ⋅ da (12)

As mentioned above, the specimen is loaded initially so that the energy release rate at the front, GII, 

reaches about 50% of GIIc corresponding to a GIImax= G50=0.392 kJ/m2 in the current study as shown 

in the Paris Law plot of Figure 18 (solid red square).

For practical applications, equation (3) can be replaced by an incremental equivalent expression

€

Δa

ΔN = c ⋅Gmax

n (13)

where for the current study, increments of Δa=0.1 mm were chosen. Starting at the initial 

delamination length, a0=25.4 mm, the energy release rates, Gi,max, were obtained for each increment, 

i, from the curve fit (open red circles and dashed red line) plotted in Figure 14. These energy release 

rate values were then used to obtain the increase in delamination length per cycle or growth rate 

Δa/ΔN from the Paris Law in Figure 18. The growth rate first increased before reaching a peak and 

then decreased with increasing a* as shown in Figure 19 (Q50,max, dashed red line), where a* is the 

increase in delamination length. The values corresponding to 40%, 30% and 20% of GIIc were also 

determined using the same approach and were included in the plot of Figure 19 (Q40,max, Q30,max, and 

Q20,max – dashed lines). For all load levels, the number of cycles during delamination growth, NG, 

was calculated by summing the increments ΔNi

€

 NG = ΔNi

i=1

k

∑ = 1

c

Gi,max

-n ⋅ Δa

i=1

k

∑ (14)

where k is the number of increments. The corresponding delamination length, a, was calculated by 

adding the incremental lengths, Δa, to the initial length, a0, 

€

 a = a0 + a* = a0 + Δa

i=1

k

∑ = a0 + k ⋅ Δa . (15)

For the delamination growth phase, the increase in delamination length, a*, is plotted in Figure 20 

for an increasing number of load cycles NG. 

10
The analyses were repeated for the case of applied displacements. The calculated growth rates 

Δa/ΔN for w50,max, w40,max, w30,max, and w20,max are plotted in Figure 21. After an initial increase the 

growth rates reached a flat peak and then decreased rapidly with increasing a*, as shown in 

Figure 21. The number of cycles during delamination growth, NG, and the corresponding increase in 

delamination length, a*, were also calculated using equations (14) and (15), respectively. The 

results are plotted in Figure 22.

5.3 Combined fatigue delamination onset and growth 

For the combined case of delamination onset and growth, the total life, NT, may be expressed as 

 

€

NT = ND +NG (16)

where, ND, is the number of cycles to delamination onset and NG, is the number of cycles during 

delamination growth [8]. For this combined case, the increase in delamination length, a*, is plotted 

in Figure 23 for an increasing number of load cycles, NT, for all load levels (Q50,max, Q40,max, Q30,max, 

and Q20,max – dashed lines). For the first ND cycles, the delamination length remains constant 

(horizontal dashed red line), followed by a growth section where - over NG cycles - the delamination

length increases following the Paris Law (dashed red line). The total life for applied cyclic 

displacements is plotted in Figure 24 for all levels (w50,max, w40,max, w30,max, and w20,max). For the first 

ND cycles, the delamination length remains constant (horizontal dashed red line), followed by a 

growth section where - over NG cycles - the delamination length increases following the Paris Law 

(dashed red line). Once a delamination length is reached where the energy release rate drops below 

the assumed cutoff value, Gth (as shown in Figure 15), the delamination growth no longer follows 

the Paris Law (dashed grey line) and stops (horizontal dashed red line) as shown in Figure 24. For 

the applied cyclic displacements, w50,max, two cutoff values Gth=0.08 kJ/m2 and Gth=0.05 kJ/m2 were 

assumed resulting in a shift of the cutoff as shown in Figure 24.

5.4 Using the benchmark example to assess an automated analysis in a commercial FE code

The computed energy release rate may serve as an initial verification step before the 

delamination growth analysis is initiated. The computed values at the initial crack tip in a 2D model 

or along the delamination front in a three-dimensional model should reach but not exceed the target 

value of GIImax=0.78 kJ/m2. Once delamination growth starts in the model, the computed energy 

release rate should increase with increasing delamination length, a, as shown in Figures 14 and 15. 

The curve fits could therefore be used to check the computed energy release rate during 

delamination growth as demonstrated later. 

The growth rate Δa/ΔN, shown in Figures 18 and 21, may be used as a check for the correct 

implementation of the Paris Law provided this output is available. For the delamination growth 

analysis, the increase in delamination length, a*, is expected with increasing number of cycles, N, as 

shown in Figures 20 and 22. The curve fits can therefore be used as a benchmark.

A delamination length prediction analysis that accounts for delamination fatigue onset as well 

as growth should yield results that closely resemble the plots in Figures 23 and 24. The curve fits 

can therefore be used as a benchmark as demonstrated later.

11
6. STATIC ANALYSIS BENCHMARKING

6.1 Computed load/displacement behavior to verify correct input data and model response

First, static analyses were performed with propagation disabled to ensure that the input was 

correct and that the models responded as expected. For this static analysis, an applied center load of 

Q=2000 N was chosen for all models. To ensure that the model geometry and material input data 

for all models produced consistent results, the computed load-displacement behavior as shown in 

Figures 25 was evaluated. In comparison to the benchmark (solid grey circles and solid grey line) 

from Figure 8, which is based on analyses using enhanced plane strain elements (CPE4I), the results 

are in good agreement. The model where classical plane strain elements (CPE4) were used exhibits 

a slightly more compliant behavior (solid blue line) as shown in Figure 25. Results from models 

using enhanced and classical plane stress elements (CPS4I, CPS4) and solid elements (C3D8I) 

(solid green, light blue and black lines) are basically identical to the benchmark results. The results 

from the model with continuum shell elements (SC8R) exhibit a slightly stiffer response (solid red 

line) compared to the benchmark results. Based on the results, it was assumed that the geometry, 

maximum applied displacements and material input were defined correctly and all models 

adequately represented the benchmark case.

6.2 Computed energy release rate distribution across the specimen width

For all the models, the mode II strain energy release rate values were also computed for an 

applied center load Q=2000 N and plotted versus the normalized width, y/B, of the specimen as 

shown in Figure 26. The results were obtained from models shown in Figures 5 through 7. The 

results obtained from two-dimensional finite element models are plotted at the center of the 

specimen y/B =0. In comparison to the benchmark (solid grey circle) from Figure 8, which is based 

on analyses using enhanced plane strain elements (CPE4I), the results are in good agreement. Only 

the model where classical plane strain elements (CPE4) were used yielded results that were about 

11% higher (solid blue square) as shown in Figure 26. As expected, the mode I and mode III strain 

energy release rates were computed to be nearly zero and hence are not shown. The results from 

three-dimensional models (C3D8I, SC8R) indicate that qualitatively, the mode II strain energy 

release rate is fairly constant over almost the entire width of the specimen and only increases 

sharply near the edges. The results obtained from simulations where continuum shell elements 

(SC8R) had been used to model the specimen (see Figure 7b) were about 6.6% lower at the center 

of the specimen (solid red diamonds) compared to the results from solid elements (C3D8I) (open 

black triangles). These low results had not been observed previously when the energy release rate 

distribution across the width of a DCB had been computed [5]. Doubling the number of elements 

across the width did not alter the results (red Xs). Only when the number of elements used to model 

the thickness of the specimen was increased to three as in the original solid element model 

(Figure 6a) did the results obtained from continuum shell element models (red crosses) resemble the 

results obtained from solid models. However, if a model of the same fidelity through the thickness 

is required when using continuum shell elements, the advantage of simpler models and faster 

analysis is lost. Based on the results, it was assumed that, with respect to the computed energy 

release rates, the models adequately represented the benchmark case.

12
6.3 Results from automated delamination propagation analysis

6.3.1 Computed delamination propagation for applied static displacements

The propagation analysis was performed in two steps using the models shown in Figures 5 

through 7 for an initial delamination length, a0=25.4 mm. In the first step, a displacement of w=1.3 

mm was applied which nearly equaled the critical displacement, wcrit=1.42 mm, determined earlier. 

In the second step, the applied displacement was increased to w=5.0 mm. For this second step, 

automatic incrementation was used in ABAQUS® and a small increment size (0.5% of the total 

step) was chosen at the beginning of the step. To reduce the risk of numerical instability and early 

termination of the analysis, a minimum allowed increment size (10-18 of the total step) was also 

chosen. The analysis was limited to 1000 increments. 

Initially, analyses were performed using two-dimensional planar models without stabilization or 

viscous regularization. Release tolerance values were varied between the default value (reltol=0.2) 

and 0.5. Using reltol=0.2, 0.3 and 0.4, the load dropped at the critical point, but the displacement 

kept increasing with decreasing load as shown in Figure 27 (solid red, blue and green lines, 

respectively) where the computed resultant force (load Q) in the center of the ENF specimen is 

plotted versus the applied displacement w. Then, the analysis terminated early due to convergence 

problems. By increasing the release tolerance - as suggested in the ABAQUS® error in the message 

(.msg) file – it was possible to complete the analysis without an error message as shown in 

Figure 27, using a release tolerance value of 0.5 (solid purple line). As desired, the load dropped at 

the critical point, however, the displacement kept increasing with decreasing load. As the 

displacement continued to increase over 2.0 mm, the computed load/displacement path converged 

to the stable propagation branch of the benchmark result. For the stable path, a saw-tooth pattern is 

observed where the top results are in good agreement with the benchmark result. Based on the 

results, it was decided to introduce additional stabilization to obtain better agreement with the 

benchmark case.

Based on problems identified during previous analyses, automated or static stabilization was not 

used in this study [5]. The results computed when contact stabilization"(cs)"was added are plotted in 

Figure 28. For a small stabilization factor (cs=1x10-6) and a release tolerance suggested in the 

handbook (reltol=0.2) [13], the load decreased, and delamination propagation started shortly after 

reaching the critical point of the benchmark solution (thin solid red line). The load/displacement 

path then ran parallel to the constant displacement branch of the benchmark result and followed the 

stable propagation branch of the benchmark result. To reduce the observed overshoot, the release 

tolerance was reduced. For a stabilization factor of cs=1x10-6 and release tolerance reltol=0.1 

(dashed blue line), the overshoot was reduced, and the computed load/displacement path then ran 

closer to the constant displacement branch of the benchmark result. Further reducing the release 

tolerance (reltol=0.01) yielded results that were in excellent agreement with the benchmark for a 

wide range of stabilization factors (1x10-4≤ cs ≤1x10-8) (solid black, light blue, orange and green 

lines – on top of each other). For the stable path, a saw-tooth pattern was observed where the top 

results were in good agreement with the benchmark result. 

When viscous"regularization"was added to help overcome convergence issues, a value of 0.2 

was used initially for the release tolerance as suggested in the handbook [13]. For a viscosity 

coefficient damv=0.1 and a release tolerance of 0.2, the load and displacement overshot the critical 

point as shown in Figure 29 (solid red line). The load/displacement path then ran parallel to the 

constant displacement branch of the benchmark result and followed the stable propagation branch of 

the benchmark result. Lowering the viscosity coefficient (damv=1x10-3) did not improve the results 

(dashed blue line). Further lowering the viscosity coefficient (damv=1x10-4) reduced the overshoot 

13
at the critical point, however, the computed displacement started to increase before the transition 

between the constant displacement branch and the stable propagation was reached (solid black line). 

Additional reduction in the viscosity coefficient alone (reltol=0.2 and 1x10-5≤ damv ≤1x10-6) 

reduced the overshoot and shifted the results closer to the constant displacement branch of the 

benchmark (dashed light blue and green lines). Only simultaneously reducing the release tolerance 

(reltol=0.1) yielded results that were in excellent agreement with the benchmark including the 

constant displacement branch and the stable propagation branch of the benchmark result (solid 

orange line). For the stable path, a saw-tooth pattern was observed where the top results were in 

good agreement with the benchmark result. Results did not change when plane stress elements 

(CPS4I) were used for the model (solid purple line). Additionally reducing the release tolerance 

(reltol=0.01) did not yield any improvement (solid blue line) but increased the analysis time.

An alternative way to plot the benchmark is shown in Figure 30 where the applied displacement 

w is plotted versus the increase in delamination length a*. This way of presenting the results is 

shown, since it may be of advantage for large structures where local delamination propagation may 

have little effect on the global stiffness of the structure and may therefore not be visible in a global 

load/displacement plot. However, extracting the delamination length a from the finite element 

results required more manual, time consuming post-processing of the results compared to the 

relatively simple and readily available output of nodal displacements and forces. The results plotted 

in Figure 30 are selected examples that were discussed above and were shown in the global 

load/displacement plots of Figures 27 through 29. The conclusions that can be drawn from this plot 

are identical to those discussed above. 

Results obtained from three-dimensional models are shown in Figures 31 through 33. Based on 

the results from two-dimensional planar models shown above, contact stabilization or viscous 

regularization was added for all analyses to help overcome convergence issues. Analyses where 

viscous regularization was used (1x10-2≤ damv ≤1x10-6) did not converge and terminated early for a 

range of release tolerance values (0.5≤ reltol ≤0.2) that had yielded converged results earlier for 

two-dimensional planar finite element models. The results computed when a small stabilization 

factor (cs=1x10-6) was added are plotted in Figure 31. For an initial release tolerance value 

(reltol=0.2) as suggested in the user’s manual [13], the load dropped (solid blue line) and 

delamination propagation started shortly before reaching the critical point of the benchmark solution 

(solid grey line). The load/displacement path then ran parallel to the constant displacement branch 

of the benchmark result. At the transition between the constant displacement branch and the 

stable propagation branch of the benchmark result, the computed load was about 13% lower 

compared to the benchmark. As displacements increased, the results closely followed the stable 

propagation branch of the benchmark result. To attempt to improve the results, the release tolerance 

was varied. Neither increasing (reltol=0.5) nor reducing (reltol=0.1) the release tolerance had a 

significant effect on the results (red solid line and green solid line, respectively). Additionally, the 

mesh was refined, first by dividing the original element length (Δa=1.0 mm) in half (Δa=0.5 mm) 

and keeping the number of elements constant across the width. Second, the number of elements 

across the width was doubled resulting in the FE-model shown earlier in Figure 7a. Changing the 

mesh also did not have a significant effect on the results (black solid line and orange solid line, 

respectively). Deviation from the benchmark may be explained by the fact that the benchmark 

results were created using two-dimensional planar finite element models of the ENF specimen. As 

shown in Figure 26, the three-dimensional models yield an energy release rate distribution where 

the peak values near the edges are somewhat higher than in the center and higher than the results 

from two-dimensional planar finite element models. Therefore, delamination propagation is 

14
expected to start before the peak obtained from two-dimensional planar finite element models is 

reached which shifts the entire results plot towards lower loads. 

Results obtained from models where continuum shell elements (SC8R) were used to model the 

ENF specimen (Figure 7b) are plotted in Figure 32. For an initial release tolerance value of 

reltol=0.1 (solid light-blue line), the load dropped and delamination propagation started 

significantly before reaching the critical point of the benchmark solution (solid grey line). In an 

attempt to improve the results, the mesh was refined by doubling the nodes in the width direction. 

Changing the mesh increased the critical load, but overall did not significantly improve the results 

(solid orange line) compared to the benchmark or to results using solid brick elements (C3D8I) as 

shown in Figure 31. Selected examples that were discussed above are plotted in Figure 33, where 

the applied displacement w is plotted versus the increase in delamination length a*. The discrepancy 

between the results obtained from models where solid brick elements (C3D8I) had been used (solid 

red, blue and green lines) and results obtained from models where continuum shell elements (SC8R) 

had been used (solid light blue and orange lines) is significant. This type of discrepancy had not 

been observed previously when propagation for a DCB specimen had been computed [3]. The 

energy release rate distributions, as shown in Figure 26, indicate that the values obtained from 

models where continuum shell elements (SC8R) had been used are lower in the center and the same 

at the edge compared to values obtained from models where solid brick elements (C3D8I) had been 

used. Based upon the values at the edge, propagation should start in continuum shell and solid

models at the same displacement. Based upon the values in the center, propagation should start in 

the shell-element models for higher displacements. The reason for the shift of the results towards 

smaller displacements as shown in Figure 33 (solid light blue and orange lines) remains unclear and 

requires further investigation and discussion with the finite element code developers. These results 

highlight the importance of benchmarking to detect specific issues in a finite element 

implementation. 

6.3.2 Computed delamination propagation for applied static center load

The propagation analysis was performed in two steps using the models shown in Figures 5 

through 7 for an initial delamination length, a0=25.4 mm. In the first step, a center load Q=1500 N 

was applied which equaled nearly the critical load, Qcrit=1526 N, determined earlier. In the second 

step, the total load was increased (Q=1800 N). Automatic incrementation was used with a small 

increment size at the beginning (0.5% of the total step) and a very small minimum allowed 

increment (10-18 of the total step) to reduce the risk of numerical instability and early termination of 

the analysis. The analysis was limited to 5000 increments. 

The same steps discussed in the section on applied displacement were followed. Initially, 

analyses were performed using two-dimensional planar models without stabilization or viscous 

regularization. Release tolerance values were varied between the default value (reltol=0.2) and 0.7. 

For the default value (reltol=0.2), the load increase stopped at the critical point, but the analysis 

terminated immediately due to convergence problems. Even by increasing the release tolerance - as 

suggested in the ABAQUS® error in the message (.msg) file – it was not possible to complete the 

analysis without termination and reach the stable propagation path of the benchmark case. 

Therefore, additional stabilization had to be introduced in order to obtain agreement with the 

benchmark case.

The results computed when contact stabilization" (cs)"was added are plotted in Figure 34. For 

small stabilization factors (cs=1x10-4 and cs=1x10-6) and a small release tolerance (reltol=0.01), the 

load increased up to the critical point and delamination propagation started while the load remained 

15
constant (solid black and green lines). For the stable propagation path, the results were in good 

agreement with the benchmark result. Computed results did not change when plane stress elements 

(CPS4I) were used for the model and the stabilization factor was kept constant at a low value 

(cs=1x10-6). Initially, the release tolerance value was set at the default value (reltol=0.2) (dashed 

orange line). Then, the release tolerance was reduced to reltol=0.1 (dashed purple line) and further 

to reltol=0.01 (dashed green line). The results were in good agreement with the benchmark results 

over the entire load/displacement range.

When viscous"regularization"was added to help overcome convergence issues, a value of 0.2 

was used initially for the release tolerance as suggested in the handbook [13] in combination with a 

viscosity coefficient damv=0.1. The computed load/displacement behavior for this parameter 

combination was in good agreement with the benchmark case as shown in Figure 35 (solid red line). 

Reducing the viscosity coefficient to values that previously yielded good results (1x10-4≤ damv 

≤1x10-6) caused the analysis to terminate. Reducing the release tolerance to reltol=0.1 and using 

damv=0.1 yielded results (solid blue line) that were in excellent agreement with the benchmark 

(solid grey line). Further reducing the release tolerance to reltol=0.01 however, caused the analysis 

to terminate. Results did not change when plane stress elements (CPS4I) were used for the model 

(reltol=0.2 and damv=0.1) as shown by the dashed red line. 

An alternative way to plot the benchmark is shown in Figure 36 where the applied center load Q

is plotted versus the increase in delamination length a*. This way of presenting the results may be of 

advantage for large structures where local delamination propagation may have little effect on the 

global stiffness of the structure and may therefore not be visible in a global load/displacement plot. 

However, extracting the delamination length a from the finite element results required more 

manual, time consuming, post-processing of the results compared to the relatively simple and 

readily available output of nodal displacements and forces. The results plotted in Figure 36 are 

selected examples that were discussed above and were shown in the global load/displacement plots 

of Figures 34 and 35. The conclusions that can be drawn from the plot in Figure 36 are identical to 

those discussed in the two paragraphs above for Figures 34 and 35. 

Results obtained from three-dimensional models are shown in Figures 37 and 38. Based on the 

results from two-dimensional planar models shown above, only contact stabilization was added to 

help overcome convergence issues. The load/displacement result computed when a small 

stabilization factor (cs=1x10-6) was added is plotted in Figure 37. A release tolerance value 

(reltol=0.1) was used which had yielded good results in the case where two-dimensional planar 

models had been used for the analyses. The load increase stopped (solid red line) and delamination 

propagation started shortly before reaching the critical point of the benchmark solution (solid grey 

line). As mentioned earlier, deviation from the benchmark may be explained by the fact that the 

benchmark results were created using two-dimensional planar finite element models of the ENF 

specimen. As shown in Figure 26, the three-dimensional models yield energy release rate 

distributions where the peak values near the edges are somewhat higher than in the center and 

higher than the results from two-dimensional planar finite element models. Therefore, delamination 

propagation is expected to start before the peak obtained from two-dimensional planar finite 

element models is reached which shifts the entire results plot towards lower loads. The result 

discussed was also plotted in Figure 38, where the applied center load Q is plotted versus the 

increase in delamination length a*. The conclusions that can be drawn from this plot are identical to 

those discussed above. 

In summary, good agreement between analysis results and the benchmark could be achieved for 

different release tolerance values in combination with contact stabilization or viscous 

16
regularization. Selecting the appropriate input parameters, however, was not straightforward and 

often required several iterations in which the parameters had to be changed. Increasing the 

release tolerance, as suggested in the user's manual [13], may help to obtain a converged 

solution, however, leads to an undesired overshoot of the computed result compared to the 

benchmark. Therefore, a combination of release tolerance and contact stabilization or viscous 

damping is suggested to obtain more accurate results. A gradual reduction of the release 

tolerance and contact stabilization or viscous damping over several analyses is suggested. 

6.3.3 Computed delamination front shape

Besides matching the load displacement behavior of benchmark results, a delamination 

propagation analysis should also yield a delamination front shape that is representative of the actual 

failure. An example of delamination front shapes observed by opening a tested ENF specimen are 

shown in Figure 39 [6]. Starting from the initial straight delamination front which is formed by the 

edge of the polytetrafluoroethylene (PTFE) insert, the delamination appears to grow uniformly 

across the width and remains straight during the static loading used to pre-crack the specimen [6]. 

The front observed after cyclic loading does not appear uniform across the width. The front is 

somewhat jagged, suggesting that growth happens in one location then stops and continues at 

another location across the width. On average, however, the front remains straight compared to the 

pronounced thumbnail shaped fronts observed in a Double Cantilever Beam (DCB) Specimen [3]. 

A straight front across almost the entire width of the specimen is to be expected when looking at 

the computed mode II strain energy release rate distribution plotted in Figure 26. Due to peak Gvalues near the edges, delamination propagation is expected to start in those areas first. 

Delamination propagation computed using the model with a uniform mesh across the width 

(Figure 6) is shown in Figure 40. Plotted on the bottom surface (defined in Figure 6) are the 

contours of the bond state, where the delaminated section appears in blue and the intact (bonded) 

section in red. The transition between the colors (in green/yellow) indicates the location of the 

delamination front. The initial straight front is shown in Figure 40a and was added for clarification 

in Figure 40c. The first propagation was observed at the edges of the specimen (Figure 40b) as 

expected from the distribution of the energy release rate (Figure 26). Later, the delamination 

propagated across the width of the specimen with a jagged front shape as shown in Figure 40c. 

Propagation occurred in one location then stopped and continued at another location across the 

width. Subsequently, the entire front moved forward and, on average, remained straight. To improve 

the results, the element length, Δa, was divided in half (Δa=0.5 mm). The refined model and results 

are shown in Figure 41. Note that the number of elements was kept constant across the width. As 

before, delamination propagation started at the edges of the specimen (Figure 41b). Later, when the 

delamination propagated across the width of the specimen, the front shape still appeared somewhat 

jagged as shown in Figure 41c. To further improve the results, the number of elements across the 

width was doubled, resulting in the FE-model shown earlier in Figure 7a. As for the previous two 

models, delamination propagation started at the edges of the specimen (Figure 42a). Later, when the 

delamination propagated across the width of the specimen, the front was jagged locally as shown in 

Figure 42b but on average moved forward as a straight front across the entire width, which is in 

agreement with the front shape observed after testing. 

17
7. ASSESSMENT OF AUTOMATED DELAMINATION ONSET AND GROWTH 

ANALYSIS UNDER CYCLIC LOADING

7.1 Static analysis to verify correct input data and model response

First, static analyses were performed with growth disabled to ensure that the input was correct 

and that the models responded as expected. For this static analysis, only a single load cycle was 

analyzed, starting at the previously applied maximum displacement w=1.0 mm. This step was 

performed to check that the amplitude input was defined correctly and resulted in the desired 

periodic cyclic loading during the analysis. Therefore, the displacement, w, obtained as analysis 

output, was plotted versus the step time, t, as shown in Figure 43. The models using plane strain 

elements (CPE4I) (open red circles) yielded the same output as the models using plane stress 

elements (CPS4I) (blue Xs) and solid elements (C3D8I) (green crosses). Based on a comparison 

with the desired fatigue loading (grey line) it is assumed that the amplitude input was defined 

correctly and the increments are small enough to adequately represent the desired periodic fatigue 

loading shown in Figure 16. 

Second, the computed mode II energy release rate – also obtained as analysis output - was 

plotted versus the step time, t, as shown in Figure 44. As desired, the energy release rate cycled 

between the expected maximum value G50,max= 0.5 GIIc= 0.392 kJ/m2

 and minimum value 

G50,min=0.1 G50,max = 0.00392 kJ/m2

. The model using plane strain elements (CPE4I) (open red 

circles and red line) yielded the same output as the models using plane stress elements (CPS4I) 

(blue Xs and solid blue line), which is consistent with the observations made above. The results 

from the solid model (C3D8I), taken in the center of the specimen at y=0.0 (green crosses and solid 

green line) are, as expected, in agreement with the results obtained from the two-dimensional planar 

models.

7.2 Assessment of automated delamination onset and growth analysis

In Figures 45 to 54, the increase in delamination length, a*, is plotted versus the total number of 

cycles, NT, for different input parameters and models. Input parameters were varied to study the 

effect on the computed onset and growth behavior during the analysis. Based on previous results 

[5], it was decided to use an initial time increment of i0=0.001 (1/200 of a single loading cycle, 

ts=0.2s), for all analyses. Also, the solution controls were kept at fixed values. For all results shown, 

the analysis stopped when a cycle limit - used as input to terminate the analysis - was reached. 

Further details about the input parameters are discussed in the appendix where a sample input file is 

also provided.

7.2.1 Computed delamination onset and growth for applied displacement

Initially, analyses were performed for a benchmark case where a maximum cyclic 

displacement, w50, corresponding to G50,max= 0.5 GIIc= 0.392 kJ/m2

, was applied. Two-dimensional 

planar models made of plane strain elements (CPE4I) with an element length, Δa=1.0 mm, at the 

delamination tip were used. The analysis stopped when a cycle limit, 108

, was reached.

First, the input to define the cyclic loading was varied. In order to define the cyclic applied 

displacement, w, a set of parameters needs to be determined as shown in the equation in Figures 16 

and 43. A Fourier series is used in ABAQUS® Standard during the analysis to approximate the 

periodic cyclic loading. Based on previous results [5], it was decided to use 50 Fourier terms for 

the analysis to approximate the periodic cyclic loading. Note, that number of Fourier terms is 

required input in addition and independent of the definition of the applied cyclic displacement 

defined in equation 10 and shown in Figure 16. Further details about the required input parameters 

18
are discussed in the appendix where a sample input file is also provided. The selection of the 

starting time, t0, causes a phase-shift (see Figure A5 in the appendix for details). The results 

obtained for t0 = 0.0 and t0= -0.05 (open red circles and open green diamond) were in excellent 

agreement with the benchmark results (solid grey line), as shown in Figure 45. The results obtained 

for t0=0.05 (open blue squares) qualitatively follows the benchmark but is shifted towards higher 

number of cycles. For all results shown, the predicted onset occurs after the benchmark result onset, 

ND50=80 cycles. The threshold cutoff, where delamination growth is terminated and the 

delamination length remains constant, is predicted close to the number of cycles defined by the 

benchmark.

Second, the effect of mesh size was studied. Therefore, the element length, Δa, at the 

delamination tip was varied as shown in the finite element models depicted in Figure 5. For this 

analysis, the cutoff value was also changed from Gth=0.08 kJ/m2 to Gth=0.05 kJ/m2 (see Figure 24). 

This change served as an additional checkpoint to assure that the automated growth was terminated 

and the delamination length remained constant once the cutoff value was reached. The results 

obtained from the respective models are plotted in Figures 46. Excellent agreement with the 

benchmark curve (grey solid line) could be achieved for all element lengths studied. For smaller 

element length, Δa=0.25 mm (open green diamonds and solid green line), the predicted onset occurs 

at a slightly lower number of cycles compared to the larger element length of Δa=1.0 mm (open red 

circles and solid red line) and Δa=0.5 mm (open blue squares and solid blue line). The observed 

mismatch is largely due to the increased element length, which causes the first growth step to be 

larger compared to the results obtained from shorter elements. For all models, the threshold cutoff, 

where delamination growth is expected to stop and the delamination length remains constant, is 

predicted close to the number of cycles defined by the benchmark.

The analyses were repeated for the other benchmark cases with applied maximum cyclic 

displacement, w40 (dashed red line), w30 (dashed blue line), and w20 (dashed green line), as shown in 

Figure 47. Two-dimensional planar models - made of plane strain elements (CPE4I) with an 

element length, Δa=1.0 mm, at the delamination tip - were used. As before, it was decided to use 50 

Fourier terms to approximate the periodic cyclic loading. A starting time, t0=0.0, was selected 

based on the results from above. The computed results for w40 (open red circles and solid red line), 

w30 (open blue squares and solid blue line), and w20 (open green diamond and solid green line) are in 

excellent agreement with the respective benchmark cases as shown in Figure 47.

The results obtained for models made of three-dimensional solid elements (shown in Figures 6 

and 7) are presented in Figure 48. The observed delamination fronts during growth, which will be 

discussed later, were not as uniform across the width of the specimen as would be expected from the 

static propagation (see Figures 38-42) or from the distribution of the energy release rate across the 

width (see Figure 26). Therefore, the increase in delamination length, a*, was plotted for three 

selected locations: the center of the specimen at y=0.0 (open circles); edge 1 at y=0.5 (open 

squares); and edge 2 at y= -0.5 (open diamonds), as shown in Figure 48. To reduce the 

computational effort, the analysis was initially performed for an element length, Δa=1.0 mm and 

only one Fourier term to approximate the periodic cyclic loading (results in red – CPU time 8.23 

105 s or 9.5 days*). Good agreement with the benchmark curve (grey solid line), however, could be 

achieved only when the number of Fourier terms to approximate the periodic cyclic loading was 

 

*

 The analyses were performed on a Quad-Core AMD Opteron(tm) Processor 8378 running openSUSE 11.3 (x86_64) 

using ABAQUS® Standard 6.9 and 6.10 on a single CPU.

19
increased to 50 as before (results in blue), which increased the computation time somewhat to 8.60 

105

 s (4.3 % increase). 

To improve the results, the element length, Δa, was divided in half (Δa=0.5 mm). Note that the 

number of elements was kept constant across the width. The results obtained when one Fourier term 

was used deviated from the benchmark as shown in Figure 49 (in green). Increasing the number of 

Fourier terms to approximate the periodic cyclic loading to 50 improved the results somewhat 

(purple lines and open symbols). The onset and initial growth were captured sufficiently accurately. 

Results for all locations (center and both edges), however, started to deviate from the benchmark 

after the initial growth phase of the analysis. To further improve the results, the number of elements 

across the width was doubled resulting in the FE-model shown earlier in Figure 7a. The results 

started to deviate from the benchmark very early in the analysis (in black). The observed 

performance for the refined meshes remains unclear and further discussion with the software 

developers is required.

7.2.2 Computed delamination onset and growth for applied center load

Analyses were also performed for benchmark cases where a maximum cyclic center load, Q, 

was applied to the model of the specimen. Two-dimensional planar models made of plane strain 

elements (CPE4I) with an element length, Δa=1.0 mm, at the delamination tip were used first and a 

center load, Q50, corresponding to G50,max= 0.5 GIIc= 0.392 kJ/m2 was applied. The analysis stopped 

when a cycle limit of 6000 was reached.

First, the input to define the cyclic loading was varied. In order to define the cyclic applied 

load, Q, a set of parameters needs to be determined as shown in the equation in Figures 16 and 43. 

As mentioned above, it was decided to use 50 Fourier terms to approximate the periodic cyclic 

loading. The selection of the starting time, t0, causes a phase-shift (see Figure A5 in the appendix 

for details). The results obtained for t0=-0.05 and t0=0.0 (open red circles and open green diamond) 

were in excellent agreement with the benchmark results (solid grey line), as shown in Figure 50. 

The results obtained for t0=0.05 (open blue squares) qualitatively follows the benchmark but is 

shifted towards a higher number of cycles. For all results shown, the predicted onset occurs after the 

benchmark result onset, ND50=80 cycles. 

Second, the effect of mesh size was studied. Therefore, the element length, Δa, at the 

delamination tip was varied as shown in the finite element models depicted in Figure 5. The results 

obtained from the respective models are plotted in Figure 51. Excellent agreement with the 

benchmark curve (grey solid line) could be achieved for both element lengths studied. For the 

smaller element length, Δa=0.5 mm, (open blue squares and solid blue line), the predicted onset 

occurs at a slightly lower number of cycles compared to an element length of Δa=1.0 mm (open red 

circles and solid red line) and is closer to the benchmark. The observed mismatch is largely due to

the increased element length, which causes the first growth step to be larger compared to the results 

obtained from shorter elements.

The analyses were repeated for the other benchmark cases with applied maximum cyclic center 

loads, Q40 (dashed red line), Q30 (dashed blue line), and Q20 (dashed green line), as shown in 

Figure 52. Two-dimensional planar models - made of plane strain elements (CPE4I) with an 

element length, Δa=1.0 mm, at the delamination tip - were used. As before, it was decided to use 50 

Fourier terms to approximate the periodic cyclic loading. A starting time, t0=0.0, was selected 

based on the results from above. The computed results for Q40 (open red circles and solid red line), 

Q30 (open blue squares and solid blue line), and Q20 (open green diamond and solid green line) are 

in excellent agreement with the respective benchmark cases as shown in Figure 52.

20
The results obtained for models made of three-dimensional solid elements (shown in Figures 6 

and 7) are plotted in Figure 53. The observed delamination fronts during growth, which will be 

discussed later, were not as uniform across the width of the specimen as would be expected from the 

static propagation (see Figures 38-42) and the distribution of the energy release rate across the width

(see Figure 26). Therefore, the increase in delamination length, a*, was plotted for three locations: 

the center of the specimen at y=0.0 (open circles); edge 1 at y= 0.5 (open squares); and edge 2 at 

y= -0.5 (open diamonds), as shown in Figure 53. As before, the analysis was initially performed for 

an element length, Δa=1.0 mm and only one Fourier term to approximate the periodic cyclic 

loading in order to reduce the computational effort (results in red – CPU time 2.8 105 s). The results 

in the center and at edge 1 (y= 0.5) were in good agreement with the benchmark curve (grey solid 

line). Upon increasing the number of Fourier terms to 50, the results in the center were not 

significantly affected (open circles). However, the results at the edges reversed with edge 2 (y= -0.5) 

edge being in good agreement with the benchmark curve (open blue diamonds) and the results at 

edge 1 (y= 0.5) slowing prematurely (open blue squares). The observed performance remains 

unclear and further discussion with the software developers is required. The corresponding 

delamination front shapes during growth will be discussed later.

To improve the results, the element length, Δa, was divided in half (Δa=0.5 mm). Note that the 

number of elements was kept constant across the width. The results obtained when one Fourier term 

was used deviated from the benchmark as shown in Figure 54 (in green). The results in the center 

and at the edges were shifted towards a lower number of cycles compared to previous results from 

two-dimensional analyses and the benchmark (solid grey line). Increasing the number of Fourier 

terms to 50 did not improve the results (purple line and open symbols). Only the onset and initial 

growth was captured sufficiently accurately. The computational effort was substantial and reached 

1.83 106 s (21 days). For further improvement, the number of elements across the width was 

doubled resulting in the FE-model shown earlier in Figure 7a. This model had yielded results that 

had been in excellent agreement with the benchmark in the static case discussed above. However, 

for the analysis of onset and growth under a cyclic center load the results only followed the 

benchmark during the onset and initial growth phase (black lines and open symbols) and then 

quickly deviated from the benchmark. The analysis was terminated by the user after 1.44 106 s (16.6 

days). The corresponding delamination fronts during growth will be discussed later. The observed 

performance remains unclear and further discussion with the software developers is required. Long 

analysis times may be avoided if the node release is made independent of the complicated cyclic 

fatigue implementation in ABAQUS® Standard. 

7.2.3 Computed delamination front shape

Besides matching the delamination length versus number of cycles behavior of the benchmark 

results, a delamination onset and growth analysis should also yield a delamination front shape that is 

representative of the actual failure observed during testing. An example of delamination front 

shapes observed by opening a tested ENF specimen was shown previously in Figure 39 [6]. The 

front observed after cyclic loading does not appear uniform across the width. The front is somewhat 

jagged, suggesting that growth happens in one location, then stops and continues at another location 

across the width. On average, however, the front remains straight compared to the pronounced 

thumbnail shaped fronts observed in a Double Cantilever Beam (DCB) Specimen [3]. A straight 

front across almost the entire width of the specimen is to be expected when looking at the computed 

mode II strain energy release rate distribution plotted in Figure 26. Due to peaks near the edges 

delamination propagation is expected to start in those areas. 

21
Delamination growth computed using the model with a uniform mesh across the width 

(Figure 6), applied displacement, w50, and one Fourier term used in the analysis is shown in 

Figure 55. Plotted on the bottom surface (defined in Figure 6) are the contours of the bond state, 

where the delaminated section appears in blue and the intact (bonded) section in red. The transition 

between the colors (in green/yellow) indicates the location of the delamination front. The initial 

straight front was added for clarification in Figures 55a, c, d and e. The first growth was observed at 

the edges of the specimen (Figure 55b) as expected from the distribution of the energy release rate 

(Figure 26). The front shown in Figure 55b corresponds to a total number of NT =400 cycles. Later, 

the delamination propagated across the entire width of the specimen and the front started to advance 

in the center first as shown in Figure 55c. Further growth occurred at the center and edge 1 at y=0.5, 

which resulted in a curved front as shown in Figure 55d. Later, growth occurred predominately at 

edge 2 (y= -0.5), resulting in the slightly slanted final front after NT =108 cycles shown in 

Figure 55e. The images of the delamination fronts help visualize the results plotted previously in 

Figure 48.

To improve the results, 50 Fourier terms were used as discussed earlier and the resulting 

computed growth is shown in Figure 56. As before, the first growth was observed at the edges of the 

specimen (Figure 56b) as expected from the distribution of the energy release rate (Figure 26). 

Later, the delamination propagated across the entire width of the specimen and the front started to 

advance in the center first as shown in Figure 56c. As the delamination grew, the front always 

turned into a straight front across the entire width before new growth started from the center as 

shown in Figure 56d. The final front after NT =108

 cycles shown in Figure 56e is in good agreement 

with the experimental finding shown in Figure 39. 

For a smaller element length (Δa=0.5 mm), selected details of the growing fronts are shown in 

Figures 57 through 59. Delamination growth computed using applied displacement, w50, and one 

Fourier term used in the analysis is shown in Figure 57. Note that the number of elements was kept 

constant across the width. As before, delamination growth initially started at the edges of the 

specimen. Later, the delamination grew across the width of the specimen, however, the center and 

both edges lagged behind as shown in Figure 57a. The front retained this double-curved shape 

throughout the analysis, as shown in Figure 57b. To improve the results, 50 Fourier terms were used 

in the analysis and the resulting fronts are shown in Figure 58. Again, delamination growth initially 

started at the edges of the specimen. As the delamination grew across the width of the specimen, 

however, the center and edge 1 (y=0.5) lagged behind as shown in Figure 58a. Later, growth 

occurred predominately at that edge, resulting in the changed curved front shown in Figure 58b. To 

further improve the results, the number of elements across the width was doubled, resulting in the 

FE-model shown earlier in Figure 7a. As for the previous models, delamination propagation started 

initially at the edges of the specimen. Later, when the delamination propagated across the width of 

the specimen the front is jagged locally as shown in Figure 59a. The front remained irregular as 

shown in Figure 59b until the analysis was terminated after 35 days. The observed performance for 

the refined meshes remains unclear and further discussion with the software developers is required.

Delamination growth computed using the model with a uniform mesh across the width (Δa=1.0 

mm, Figure 6); applied load, Q50; and one Fourier term used in the analysis is shown in Figure 60. 

The initial straight front was added for clarification in Figures 60a, c, d and e. The first growth was 

observed at the edges of the specimen (Figure 60b) as expected from the distribution of the energy 

release rate (Figure 26). Later, the delamination propagated across the entire width of the specimen 

and the front started to advance faster at edge 2 (y= -0.5) developing into a somewhat slanted front, 

as shown in Figure 60c. Then growth occurred more in the center and at edge 1 (y=0.5), which 

22
resulted in a curved front, as shown in Figure 60d. The final front after NT = 6000 cycles shown in 

Figure 60e remained curved with the growth at edge 2 (y= -0.5) lagging behind. The images of the 

delamination fronts help visualize the results plotted in Figure 53.

For a smaller element length (Δa=0.5 mm), selected details of the growing fronts are shown in 

Figures 61 through 63. Delamination growth computed using applied load, Q50; and one Fourier 

term used in the analysis is shown in Figure 61. Note that the number of elements was kept constant 

across the width. As before, delamination growth initially started at the edges of the specimen. 

Later, the delamination started to advance faster at edge 1 (y=0.5), developing into a somewhat 

slanted front as shown in Figure 61a. The front shape later advanced at edge 2 (y= -0.5) and the 

center and turned into a double-curved front as shown in Figure 61b at the end of the analysis after 

6000 cycles. To improve the results, 50 Fourier terms were used in the analysis and the resulting 

fronts are shown in Figure 62. Again, delamination growth initially started at the edges of the 

specimen. As the delamination grew across the width of the specimen, however, the delamination 

started to advance faster at edge 1 (y=0.5), developing into a somewhat slanted front, as shown in 

Figure 62a. Later, growth occurred predominately at edge 2 (y= -0.5), resulting in the final curved 

front after NT =6000 cycles shown in Figure 62b. To improve the results, the number of elements 

across the width was doubled resulting in the FE-model shown earlier in Figure 7a. As for the 

previous models, delamination propagation started initially at the edges of the specimen as shown in 

Figure 63a. Later, when the delamination propagated across the width of the specimen, the front 

was jagged locally and remained irregular as shown in Figure 63b until the analysis was terminated 

after 1.44 106 s CPU time (NT = 9540 cycles). The observed performance for the refined meshes 

remains unclear and further discussion with the software developers is required.

8. SUMMARY AND CONCLUSIONS

The development of benchmark examples, which allow the assessment of the static 

delamination propagation as well as the growth prediction capabilities was presented and 

demonstrated for ABAQUS® Standard. The example is based on a finite element model of mode II 

End Notched Flexure (ENF) specimen. The model is independent of the analysis software used and 

allows the assessment of the delamination growth prediction capabilities in commercial finite 

element codes. 

First, static benchmark results were created based on the approach developed in reference 3, 

using two-dimensional finite element models for simulating the ENF specimens with different 

initial delamination lengths, a0. For each delamination length modeled, the load and displacements 

were monitored. The mode II strain energy release rate was calculated for a fixed applied load. It 

was assumed that the delamination propagated when the mode II strain energy release rate reached 

the fracture toughness value. Thus, critical loads and critical displacements for delamination 

propagation were calculated for each initial delamination length modeled. From these critical 

load/displacement results, benchmark solutions were created. It was assumed that the 

load/displacement relationship computed during automatic propagation should closely match the 

benchmark cases. 

Second, the development of a benchmark example for delamination fatigue growth prediction 

was presented step by step. The number of cycles to delamination onset was calculated from the 

material data for mode II fatigue delamination growth onset from reference 6. The number of cycles 

during stable delamination growth was obtained incrementally from the material data for mode II 

fatigue delamination propagation. For the combined benchmark case of delamination onset and 

23
growth, the delamination length was calculated for an increasing total number of load cycles. 

Second, starting from an initially straight front, the delamination was allowed to grow under cyclic 

loading. The number of cycles to delamination onset and the number of cycles during stable 

delamination growth for each growth increment were obtained from the analysis.

After creating the benchmark cases, the approach was demonstrated for the commercial finite 

element code ABAQUS®. Starting from an initially straight front, the delamination was allowed to 

propagate under static loading or grow under cyclic loading based on the algorithms implemented 

into the software. Input control parameters were varied to study the effect on the computed 

delamination propagation and growth. 

The results showed the following:

• The benchmarking procedure was capable of highlighting the issues associated with 

the input parameters of a particular implementation. 

• In general, good agreement between the results obtained from the propagation and 

growth analysis and the benchmark results could be achieved by selecting the 

appropriate input parameters. However, selecting the appropriate input parameters 

was not straightforward and often required an iterative procedure. 

• Different sets of input parameters were identified to be important for the automated 

propagation analysis under static loading compared to the automated onset and 

growth analysis under cyclic loading. 

• The results for automated delamination propagation analysis under static loading 

showed the following: 

o Increasing the release tolerance as suggested in the ABAQUS® handbook 

helped to obtain a converged solution, however, it led to undesired overshoot 

of the result. 

o A combination of release tolerance and contact stabilization or viscous 

damping is required to obtain more accurate results. 

o A gradual reduction of the release tolerance and contact stabilization or 

viscous damping over several analyses is suggested. 

o Results obtained from continuum shell models were offset compared to results 

obtained from three-dimensional solid models. Energy release rates were lower 

in the center of the specimen. Values were higher at the edges causing 

propagation to begin at lower loads and displacements than the solid models or 

the benchmark case.

o Computed delamination fronts obtained from three-dimensional solid models 

were somewhat jagged. Fronts obtained from refined meshes matched the 

experimentally observed shapes.

• The results for automated delamination growth analysis under static loading showed 

the following: 

o Stabilization or viscous damping is not required. The release tolerance has no 

effect.

o The onset prediction appeared much more sensitive to the input parameters 

than the growth prediction. 

o Best results were obtained when a sine curve representation of the cyclic applied 

displacement was selected in combination with the starting time, t0=0.0. 

24
o Consistent results were obtained when input parameters were selected such 

that 50 terms in the Fourier series were used during the execution of 

ABAQUS® Standard to approximate the periodic cyclic loading. 

o Results obtained from three-dimensional solid models appeared mesh 

dependent. Results obtained from refined meshes deviated from the benchmark 

compared to results from coarser meshes. Computed delamination fronts did not 

match the experimentally observed shape. Since the reason for this discrepancy 

is unclear, further discussion with the software developers is required.

o The analyses for three-dimensional models of a simple ENF specimen required 

many days of computation time. Improvements to the implementation by 

making the node release independent of the complicated cyclic fatigue 

implementation may alleviate this problem. 

Overall, the results are promising and the current findings concur with previously published 

conclusions [3,4,5]. Further studies, however, are required which should include the assessment of 

the propagation capabilities in more complex mixed-mode specimens and on a structural level.

Assessing the implementation in one particular finite element code illustrated the value of 

establishing benchmark solutions since each code requires specific input parameters unique to its 

implementation. Once the parameters have been identified, they may then be used as a starting point 

to model delamination growth for more complex configurations.

ACKLOWLEDGEMENTS

This research was partially supported by the Aircraft Aging and Durability Project as part of 

NASA’s Aviation Safety Program, and by the Subsonic Rotary Wing Project as part of NASA’s 

Fundamental Aeronautics Program. 

The analyses were performed at the Durability, Damage Tolerance and Reliability Branch at 

NASA Langley Research Center, Hampton, Virginia, USA. 

REFERENCES

1. E. F. Rybicki and M. F. Kanninen, "A Finite Element Calculation of Stress Intensity Factors by a Modified 

Crack Closure Integral," Eng. Fracture Mech., Vol. 9, pp. 931-938, 1977.

2. R. Krueger, "Virtual Crack Closure Technique: History, Approach and Applications," Applied Mechanics 

Reviews, Vol. 57, pp. 109-143, 2004.

3. R. Krueger, "An Approach to Assess Delamination Propagation Simulation Capabilities in Commercial Finite 

Element Codes," NASA/TM-2008-215123, 2008.

4. A. C. Orifici and R. Krueger, “Assessment of Static Delamination Propagation Capabilities in Commercial 

Finite Element Codes Using Benchmark Analysis,” NASA/CR-2010-216709, NIA report no. 2010-03, 2010.

5. R. Krueger, “Development of a Benchmark Example for Delamination Fatigue Growth Prediction,” 

NASA/CR-2010-216723, NIA report no. 2010-04, 2010.

6. T. K. O'Brien, W. M. Johnston, and G. Toland, “Mode II Interlaminar Fracture Toughness and Fatigue 

Characterization of a Graphite Epoxy Composite Material,” NASA/TM-2010-216838, 2010.

7. I. S. Raju and T. K. O'Brien, "Fracture Mechanics Concepts, Stress Fields, Strain Energy Release Rates, 

Delamination and Growth Criteria," in Delamination Behavior of Composites, S. Sridharan, Ed.: Woodhead 

Publishing in Materials, 2008.

8. "Composite Fatigue Damage Onset," in Composite Materials Handbook CMH-17, Vol. 3, Section 12.6.4, 

draft of revision G, 2009.

25
9. P. Hansen and R. Martin, "DCB, 4ENF and MMB Delamination Characterisation of S2/8552 and IM7/8552," 

Materials Engineering Research Laboratory Ltd. (MERL), Hertford, UK N68171-98-M-5177, 1999.

10. M. L. Benzeggagh and M. Kenane, "Measurement of Mixed-Mode Delamination Fracture Toughness of 

Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending Apparatus," Composites Science and 

Technology, Vol. 56, pp. 439-449, 1996.

11. "KaleidaGraph: Version 4.1, 2009.

12. "Fatigue Fracture Toughness," in Composite Materials Handbook CMH-17. Rev. G, Vol. 1, Section 6.9.4, 

ASTM International, 2010.

13. Abaqus Analysis User’s Manual, ABAQUS® Standard, Version 6.10, DSS Simulia, 2010.

14. Abaqus Theory Manual, ABAQUS® Standard , Version 6.10, DSS Simulia, 2010.

15. D. Broek, The Practical Use of Fracture Mechanics: Kluwer Academic Publishers, 1991.

TABLE I. MATERIAL PROPERTIES [8].

Unidirectional Graphite/Epoxy Prepreg 

E11 = 161 GPa E22 = 11.38 GPa E33 = 11.38 GPa

ν12 = 0.32 ν13 = 0.32 ν23 = 0.45

G12 = 5.2 GPa G13 = 5.2 GPa G23 = 3.9 GPa

The material properties are given with reference to the ply coordinate axes where index 11 denotes the ply principal 

axis that coincides with the direction of maximum in-plane Young’s modulus (fiber direction). Index 22 denotes the 

direction transverse to the fiber in the plane of the lamina and index 33 the direction perpendicular to the plane of the 

lamina.

TABLE II. FRACTURE PARAMETERS.

Fracture Toughness Data [6,9] – Figures 2, A1

GIc = 0.21 kJ/m2 GIIc = GIIIc =0.78 kJ/m2 η= 2.57

Delamination Growth Onset Data [6] – Figures 3, 17, A3,

m0=0.78 m1=-0.16

c1=0.213 c2=-6.329

Delamination Growth Rate Data (Paris Law) [6] – Figures 4, 18, A5

GIIc = 0.78 kJ/m2 Gth = 0.08 kJ/m2

c=0.33 n=5.55

26
APPENDIX

Delamination fatigue growth analysis in ABAQUS® 

Delamination growth at the interfaces in laminated composites subjected to cyclic loadings 

can be simulated in ABAQUS® by specifying the low-cycle fatigue criterion propagation analysis 

using the direct cyclic approach [13]. The interface along which the delamination (or crack) 

propagates must be indicated in the model using a fracture criterion definition. The onset and 

growth of fatigue delamination at the interfaces are characterized by the relative fracture energy 

release rate. The fracture energy release rates at the crack tips in the interface elements are 

calculated based on the virtual crack closure technique (VCCT). The low-cycle fatigue analysis 

in ABAQUS® is a quasi-static analysis on a structure subjected to sub-critical cyclic loading. The 

low-cycle fatigue analysis in ABAQUS® uses the direct cyclic approach to obtain the stabilized 

cyclic response of the structure directly. The direct cyclic analysis uses a combination of Fourier 

series and time integration of the nonlinear material behavior to obtain the stabilized cyclic 

response of the structure iteratively and therefore avoids the numerical expense associated with a 

transient analysis. The direct cyclic analysis in ABAQUS® is therefore suited for very large 

problems in which many load cycles must be applied to obtain the stabilized response. The direct 

cyclic analysis in ABAQUS®, however, is limited to geometrically linear behavior and fixed 

contact conditions. The theory and algorithm to obtain a stabilized response using the direct 

cyclic approach are described in detail in the ABAQUS® Theory Manual [14]. 

Required input for ABAQUS®

The input required to perform a delamination onset and growth analysis in ABAQUS®

Standard is discussed in the following paragraphs. It is assumed that the reader is familiar with 

ABAQUS® Standard and the syntax used in the input file (.inp). The focus is therefore on the 

specific input that relates to delamination propagation, low-cycle fatigue and the direct cyclic 

approach in ABAQUS®. An example input file is given at the end of this appendix to provide an 

overview of an entire analysis and assist the readers in creating their own analyses. 

Input for delamination propagation

The interface along which the delamination (or crack) propagates must be indicated in the 

model using a fracture criterion definition: 

*DEBOND,SLAVE=VCCT_TOP,MASTER=VCCT_BOT,FREQ=1

*FRACTURE CRITERION,TYPE=fatigue,MIXED MODE BEHAVIOR=BK, TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

where VCCT_TOP and VCCT_BOT are interface surfaces as shown in Figure A1 and <tol> is the 

tolerance within which the crack propagation criterion must be satisfied. The input parameters 

for the fracture criterion are obtained from the static mixed-mode failure criterion, the 

delamination growth onset criterion and the growth rate shown in Figures A2 –A5. 

The critical energy release rates <GIc>,<GIIc>,<GIIIc> and the curve fit parameter <eta> are 

obtained from the mixed-mode failure criterion as shown in Figure A2. A quasi static mixed-mode 

fracture criterion is determined for a material by plotting the interlaminar fracture toughness, Gc, 

versus the mixed-mode ratio, GII/GT as discussed earlier and shown in Figure A2.

The parameters <c1>,<c2> are obtained by solving the law for growth onset (shown in Figure 

A3) for the number of cycles N as shown in equation (8). The parameters <c3>,<c4> are obtained 

directly from the Paris Law as shown in Figure A4. The parameter <r1> is calculated from the 

27
energy release rate cutoff value, Gth, and the fracture toughness, GIc, as shown in Figure A4. To 

calculate the parameter <r2>, the user needs to define an energy release rate upper limit, Gpl, 

above which the fatigue crack will grow at an accelerated rate as shown in Figure A4. For the 

current benchmark example, Gpl was chosen to be 90% of the fracture toughness. 

Input for cyclic loading

Defining a low-cycle fatigue analysis using the direct cyclic approach in ABAQUS®

Standard requires the definition of an amplitude curve which describes the relative load 

magnitude: 

*amplitude,name=test, DEFINITION=PERIODIC

<n>,<omega>,<t0>,<A0> 

<A1>,<B1>

where test is the label to be used to refer to the amplitude curve. The parameters defined in the 

first line are the number of terms in the Fourier series, <n>, the circular frequency, <omega>, in 

radians per time, the starting time, <t0>, and the constant term in the Fourier series, <A0>, as 

shown in Figure A5. Using different starting times, <t0>, causes a phase-shift in the amplitude 

curve (dashed lines) as shown in Figure A5. The parameters defined in the second line are the 

first coefficient of the cosine terms, <A1>, and the first coefficient of the sine terms, <B1>, also 

shown in Figure A5. 

The amplitude curve is then referenced in the definition of the cyclic loading. In the current 

example, prescribed displacements were used to simulate the displacements of the ENF specimen:

*BOUNDARY,AMPLITUDE=test

LFRONT, 1, 1, -1.0

where LFRONT is the node set located at the center of the specimen where the displacements are 

applied as shown in Figure A1. The factor 1.0 is used to multiply the relative magnitude defined by 

the amplitude curve (shown in Figure A5) and obtain the applied cyclic displacement, w, as shown 

in Figure A6.

The direct cyclic approach in ABAQUS® Standard is used to obtain the stabilized cyclic 

response of a structure directly: 

*direct cyclic,fatigue

 <i0>,<tS>,,,<ni>,<nmax>,<Δn>, <imax>,

,,<NT>,,

where the parameter fatigue is used to perform a low-cycle fatigue analysis. The parameters 

defined in the first line are the initial time increment, <i0>, the time of a single loading cycle, 

<tS> (as shown in Figure A6), the minimum time increment allowed (not used), the maximum 

time increment allowed (not used), the initial number of terms in the Fourier series, <ni>, the 

maximum number of terms in the Fourier series, <nmax>, the increment in number of terms in the 

Fourier series, <Δn> , and the maximum number of iterations allowed in a step, <imax>. The 

parameters defined in the second line are the minimum increment in number of cycles over 

which the damage is extrapolated forward (default used), the maximum increment in number of 

cycles over which the damage is extrapolated forward (default used), the total number of cycles 

allowed in a step, <NT>, and the damage extrapolation tolerance (default used). The time of a 

single loading cycle was kept constant at ts=0.2 s for all analyses. Most analyses were run up to 

NT=108 cycles in order to reach the threshold after which delamination growth stops as shown in 

28
Figure 24. All other input parameters were varied and the effect on the results was studied as 

discussed in the main part of this report.

Control parameters direct cyclic analysis

Solution controls in ABAQUS® Standard can be reset and modified by using keyword

*controls,type=direct cyclic

<IPI>, <CRn> ,<CUn>,<CR0>,<CU0>,

where the parameter direct cyclic is used to set parameters that will be used to control the 

stabilized state and plastic ratcheting detections and to specify when to impose the periodicity 

condition for direct cyclic analysis. If this keyword is omitted, default parameters are used. The 

parameters defined in the first line are the iteration number at which the periodicity condition is 

first imposed, <IPI>, (default used), the stabilized state detection criterion for the ratio of the 

largest residual coefficient on any terms in the Fourier series to the corresponding average flux 

norm, <CRn>, the stabilized state detection criterion for the ratio of the largest correction to the 

displacement coefficient on any terms in the Fourier series to the largest displacement 

coefficient, <CUn>, plastic ratcheting detection criterion for the ratio of the largest residual 

coefficient on the constant term in the Fourier series to the corresponding average flux norm, 

<CR0>, and the plastic ratcheting detection criterion for the ratio of the largest correction to the 

displacement coefficient on the constant term in the Fourier series to the largest displacement 

coefficient, <CU0> [13]. 

Example input files 

An input file is given to provide an overview of an entire analysis and assist the readers in 

creating their own analyses. The analysis was divided into two steps. In the first step, a small 

static preload step was introduced as a work around to avoid problems discovered with the initial 

contact conditions (ABAQUS® bug v68_1987). The second step performed the desired cyclic 

analysis. It was found that the prescribed displacements in the static preload step had to be small 

(0.001 mm) compared to the prescribed displacements (vmin =0.1, vmax=-1.0), which were used to 

simulate the cyclic opening of the arms of the ENF specimen. In a previous study, analyses using 

larger preload steps did not converge in the second step and terminated prematurely [5].

For all analyses, the input to define the fracture criterion (<GIc>,<GIIc>,<GIIIc>,<eta>), the 

parameters for delamination onset (<c1>,<c2>), and delamination growth (Paris Law) 

(<c3>,<c4>,<r1>,<r2>) were kept constant. The parameters to define the load frequency 

(<omega>,<A0>) as well as the minimum and maximum applied displacement (<tS>, vmax=-1.0)

were also kept constant during all analyses. Based on previous results [5], it was decided to use an 

initial time increment of i0=0.001 (1/200 of a single loading cycle, ts=0.2s), for all analyses. Also 

the solution controls were kept at fixed values. Other parameters required to define the Fourier 

series, which is used to define the cyclic load, were varied. The ABAQUS® keywords shown in 

bold type were discussed in detail in the previous paragraphs. 

29
Input file for fatigue onset and growth analysis

 *HEADING

ENF-UD-IM7/8552, a=10 mm 

units: mm, N, MPa 

*** elements, nodes, material, etc

*ELEMENT, TYPE=CPE4I , ELSET=EALL

... 

... 

*NSET, NSET=BONDED, GENERATE

 205, 845, 8

**** surface and contact definition ****

*SURFACE, TYPE=ELEMENT, NAME=VCCT_BOT

EL_BOT, S3

*SURFACE, TYPE=ELEMENT, NAME=VCCT_TOP

EL_TOP, S1

*INITIAL CONDITION, TYPE=CONTACT

VCCT_TOP, VCCT_BOT, BONDED

*CONTACT PAIR, INTERACTION=VCCT, ADJUST=BONDED

VCCT_TOP, VCCT_BOT

**** VCCT -- STAGE I ****

*SURFACE INTERACTION, NAME=VCCT

 <width>

**** VCCT fatigue input

*parameter

** Damage and tolerance parameters

tol=0.001 

** Fracture toughness:

 GIc = 0.208

 GIIc = 0.78

 GIIIc = 0.78

** B-K parameter:

 eta=2.5713

** width in the plane stress/strain direction

 width =25.4

** fatigue crack growth data **

 c1=0.213

 c2=-6.329

 c3=0.33185

 c4=5.5519

 r1=0.102

 r2=0.9

***

*** amplitude ****

*amplitude,name=test, DEFINITION=PERIODIC

1,31.42,0.05,0.55

0,0.45

*** history data ***

*** static ramp up step ****

*STEP, NLGEOM, INC= 10000

*STATIC

0.001, 0.001

**

*DEBOND,SLAVE=VCCT_TOP,MASTER=VCCT_BOT,FREQ=1

*FRACTURE CRITERION,TYPE=VCCT,MIXED MODE BEHAVIOR=BK

1.0e6,1.0e6,1.0e6,<eta>

30
**

  ENF loading

** * center deflection for ENF, 2D-full model

*BOUNDARY, TYPE=DISPLACEMENT

LFRONT, 1, 1, -0.001

**

*** center load for ENF, 2D-full model

 CLOAD

** * 408,1, -1.0

**

** field and history output **

*OUTPUT, FIELD, VARIABLE=PRESELECT, FREQ=1

*Output, history,VARIABLE=PRESELECT,freq=1

*NODE output,NSET=LFRONT

 RF1

*NODE output,NSET=LFRONT

 U1

*NODE output,NSET=SUPPORT

 RF1

*END STEP

**** cyclic loading step ****

*STEP, INC= 10000

*direct cyclic,fatigue

 0.001,0.2,, , 50, 50, , 20

3,6,100000000,,

*DEBOND,SLAVE=VCCT_TOP,MASTER=VCCT_BOT,FREQ=1

*FRACTURE CRITERION,TYPE=fatigue,MIXED MODE BEHAVIOR=BK, TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

*** run analysis first with default values

*controls,type=direct cyclic

,100,1E-3,100,1E-3 

**** ENF loading

*BOUNDARY,AMPLITUDE=test

LFRONT, 1, 1, -1.0

**

*** center load for ENF, 2D-full model

***CLOAD,AMPLITUDE=test

** 408,1, -682.5

**

** field and history output **

*OUTPUT, FIELD, VARIABLE=PRESELECT, FREQ=25

*ELEMENT OUTPUT

 cycleini,status,sdeg

*CONTACT OUTPUT, MASTER=VCCT_BOT, SLAVE=VCCT_TOP

 dbt,dbsf,dbs,openbc,crsts,enrrt,efenrrtr,bdstat

*Output, history,VARIABLE=PRESELECT,freq=25

*NODE output,NSET=LFRONT

 RF1

*NODE output,NSET=LFRONT

 U1

*NODE output,NSET=SUPPORT

 RF1

*END STEP

31
Q

2h

+θ

y

z

x

dimensions

B 25.4 mm

2h 4.5 mm

2L 101.6 mm

a0 25.4 mm

Figure 1. End-Notched Flexure Specimen (ENF) [6].

a0

2L L

B

w

material and layup: IM7/8552 [0]24

fatigue loading

wmax 1.0 mm

wmin 0.1 mm

R 0.1

f 5.0 Hz

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

0 0.2 0.4 0.6 0.8 1

experimental data [9]

mean values [9]

experimental data [6]

mean values [6]

DCB, Mode I MMB, Mode I and II

Gc ,

kJ/m2

mixed mode ratio GII/GT

Figure 2. Mixed mode fracture criterion.

ENF, Mode II

G

Ic

 =0.21

 GIIc = 0.78

curve fit: Gc = G

Ic

 + (GIIc-GIc)(GII/GT)

η η = 2.57

32
10-8

10-7

10-6

10-5

10-4

10-3

10-2

10-1

100

10-2 10-1 100

da/dN, 

mm/cycle

G

II,max, kJ/m2

Figure 4. Delamination growth rate (Paris Law) for mode II for IM7/8552 [6].

GIIc = 0.78

da/dN = c Gn

c = 0.33185

n = 5.5519

Paris Law [6]

cutoff

Gth = 0.08

0.00

0.20

0.40

0.60

0.80

1.00

100 101 102 103 104 105 106 107

experimental data [6]

N

D , cycle

Gmax ,

 kJ/m2 (N/mm)

Figure 3. Delamination growth onset for mode II for IM7/8552 [6].

curve fit: G= m

0

 * N 

m

0

= 0.78

m

1

= -0.16

m1

33
initial delamination front bonded nodes

2h

initial delamination front bonded nodes

 bottom surface

top surface

Δa Δa=1.0 mm

a. Deformed FE-model of ENF specimen with initial delamination before growth.

Figure 5. Two-dimensional finite element model of an ENF specimen.

d. Delamination front detail of a further refined FE-model of an ENF specimen.

X

Y Z

X

Y Z

c. Delamination front detail of a refined FE-model of an ENF specimen.

b. Delamination front detail of a FE-model of an ENF specimen.

Δa Δa=0.5 mm

Δa Δa=0.25 mm

w

2L

a0

x,u

z,w

w=0 u,w=0

34
X

Y Z

B

2L

top and bottom 

surface

bonded 

nodes

x y

z

a0

delamination front

edge 1

edge 2

a. Deformed three-dimensional FE-model with initial delamination front before growth.

Figure 6. Full three-dimensional finite element model of an ENF specimen.

b. Detail of three-dimensional FE-model around delamination front.

X

Y Z

delamination front

 bottom surface

top surface

bonded 

nodes

2h

Δa=1.0 mm

Δa

35
bonded 

nodes

initial straight delamination front

a. Deformed three-dimensional model of an ENF specimen with a fine mesh before propagation.

Figure 7. Three-dimensional finite element models of an ENF specimen.

initial straight delamination front

bonded 

nodes

Δa=0.5 mm

Δa=1.0 mm

b. Deformed continuum shell model of an ENF specimen after delamination propagation. 

edge 1

edge 2

y x

z

x y

z

36
0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

25.4

26.4

27.4

28.4

29.4

30.4

35.5

40.6

45.7

50.8

55.9

60.9

66.0

71.1

76.2

load 

Q, N

displacement w, mm

G

IIc

 = 0.78 kJ/m2=const.

G

II

=1.3 kJ/m2 2.6 3.3 4.2 4.8 3.9 2.33.0 1.11.7

Figure 9. Calculated critical load-displacement behavior for ENF specimen.

delamination length a

0

=

0.0 1.0 2.0 3.0 4.0 5.0 6.0

25.4

26.4

27.4

28.4

29.4

30.4

35.5

40.6

45.7

50.8

55.9

60.9

66.0

71.1

76.2

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

load 

Q, N

displacement w, mm

G

II

=1.3 kJ/m2 2.6 3.3 4.2 4.8 3.9 3.0 2.3 1.7 1.1

Figure 8. Load-displacement behavior for ENF specimens with 

different initial delamination lengths, a0.

delamination length a

0

=

a

0

= 25.4 mm

a

0

= 76.2 mm

37
0 5 10 15 20 25 30 35 40 45 50 55

0

1

2

3

4

5

displacement

w, mm

increase in delamination length a*, mm

G

IIc

 = 0.78 kJ/m2=const.

benchmark for applied displacement w

Figure 11. Benchmark case for applied displacement, w, plotted versus 

increase in delamination length, a*.

0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

load 

Q, N

displacement w, mm

wcrit=1.42 mm

G

IIc

= 0.78 kJ/m2=const.

Qcrit

=1526 N

displacement w

load Q

benchmark for applied

Figure 10. Benchmark cases for applied load, Q, and displacement, w, for ENF specimen.

38
0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

u

u

u

u

u

u

u

u

u

u

u

u

u

u

u

v

v

load

 Q, N

displacement w, mm

Figure 13. Maxium cyclic loads and applied displacments for an ENF specimen.

wcrit=1.42 mm

G

IIc

= 0.78 kJ/m2 = constant

Qcrit

=1526 N

Q

50,max=1079 N; w50,max=1.0 mm

Q40,max=965 N; w40,max=0.90 mm

Q30,max=835 N; w30,max=0.78 mm

Q20,max=682 N; w20,max=0.64 mm

a

0

=25.4 mm a

0

=76.2 mm

0 5 10 15 20 25 30 35 40 45 50 55

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

load 

Q, N

increase in delamination length a*, mm

G

IIc

 = 0.78 kJ/m2=const.

benchmark for applied load Q

Figure 12. Benchmark case for applied load, Q, plotted versus increase in delamination length, a*.

39
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0 5 10 15 20 25 30 35 40 45 50 55

w50,max=1.0 mm

w40,max=0.9 mm

w30,max=0.78 mm

w20,max=0.64 mm G, kJ/m2

Figure 15. Energy release rate - delamination length behavior 

for ENF specimen with constant displacement.

G

IIc

G= G =0.78

IIc

 =0.78 = const.

Gth = 0.08

increase in delamination length a*, mm

0.00

0.50

1.00

1.50

0 5 10 15 20 25 30 35 40 45 50 55

Q50,max=1079 N

Q40,max=965 N

Q30,max=835 N

Q20,max=682 N

G, kJ/m2

Figure 14. Energy release rate - delamination length behavior 

for ENF specimen with constant load.

G

IIc

 =0.78

Gth = 0.05

increase in delamination length a*, mm

G= G

IIc

 =0.78 = const.

40
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

100 101 102 103 104 105 106 107

experimental data [6]

N

D , cycle

Gmax ,

 kJ/m2

 (N/mm)

Figure 17. Fatigue growth onset for ENF specimen.

curve fit: G= m

0

 * N 

D

m

0

=0.78

m

1

=-0.16

m1

G

50 = 0.5 GIIc

 = 0.39

N

D50

=80

G

40 = 0.4 GIIc

 = 0.31

G

20 = 0.2 GIIc

 = 0.16

G

30 = 0.3 GIIc

 = 0.24

N

D20

N =26550 D30

=2040

N

D40

=330

N

D

 = c

1

*G c2

c

1 0.213

c

2 -6.329

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0

0.2

0.4

0.6

0.8

1

1.2

1.4

time t, s

applied 

displacement 

w, mm

w50,max=1.0

w

50,min

=0.1

t

s

Figure 16. Applied cyclic center deflection for ENF specimen.

R=0.1

f= 5 Hz t

s

= =0.2 s 1

f

_

w=[A0+B1 sinω(t-t

0

)]*wmax

ω=10π=31.42

A

0

=0.55

B

1

=0.45

t

0

=-0.05

41
10-5

10-4

10-3

10-2

10-1

100

101

0 5 10 15 20 25 30 35 40 45 50 55

Q50,max=1079 N

Q40,max=965 N

Q30,max=835 N

Q20,max=682 N

Δa/ΔN,

mm/cycle

increase in delamination length a*, mm

Figure 19. Delamination growth rate behavior for ENF specimen for applied cyclic loading.

10-8

10-7

10-6

10-5

10-4

10-3

10-2

10-1

100

10-2 10-1 100

da/dN, 

mm/cycle

G

II,max, kJ/m2

Figure 18. Delamination growth rate (Paris Law) for ENF specimen.

GIIc = 0.783

da/dN = cGn

c=0.33185

n=5.5519

cutoff

Gth = 0.08

Paris Law [6]

increasing a

G

50 = 0.5 GIIc

 = 0.39

G

40 = 0.4 GIIc

 = 0.31

G

20 = 0.2 GIIc

= 0.16

G

30 = 0.3 GIIc

 = 0.24

42
10-9

10-8

10-7

10-6

10-5

10-4

10-3

10-2

0 5 10 15 20 25 30 35 40 45 50 55

w50,max=1.0 mm

w40,max=0.9 mm

w30,max=0.78 mm

w20,max=0.64 mm

Δa/ΔN,

mm/cycle

increase in delamination length a*, mm

Figure 21. Delamination growth rate behavior for ENF specimen for applied center deflection.

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108

Q50,max=1079 N

Q40,max=965 N

Q30,max=835 N

Q20,max=682 N

increase in

delamination

 length 

a*, mm

number of cycles N

G

Figure 20. Delamination growth behavior for ENF specimen for applied cyclic loading.

43
0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108

Q50,max=1079 N

Q40,max=965 N

Q30,max=835 N

Q20,max=682 N

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 23. Total delamination growth behavior for ENF specimen for applied cyclic loading.

N

G

N (delamination growth) D

(delamination 

onset)

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108 109

w50,max=1.0 mm

w40,max=0.9 mm

w30,max=0.78 mm

w20,max=0.64 mm

increase in

delamination

 length 

a*, mm

number of cycles N

G

Figure 22. Stable delamination growth behavior for ENF specimen for applied center deflection.

44
0.0 0.5 1.0 1.5 2.0 2.5 3.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

CPE4I

C3D8I

CPE4

CPS4

CPS4I

SC8R

load 

Q, N

displacement w, mm

Figure 25. Load-displacement behavior for an ENF specimen obtained from different models.

wcrit=1.42 mm

Qcrit

=1526 N

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108 109

w50,max=1.0 mm

w40,max=0.9 mm

w30,max=0.78 mm

w20,max=0.64 mm

increase in

delamination

 length 

a*, mm

number of cycles N

T

N

G

(delamination growth)

N

D

(delamination 

onset)

cutoff

Figure 24. Total delamination growth behavior for ENF specimen for applied center deflection.

Gth = 0.05

0.08

45
0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

reltol=0.2

reltol=0.3

reltol=0.4

reltol=0.5

load 

Q, N

displacement w, mm

Figure 27. Computed critical load-displacement behavior for ENF specimen obtained 

from two-dimensional planar models (CPE4I) with different release tolerance settings.

wcrit=1.42 mm

Qcrit

=1526 N

benchmark

0.00

0.50

1.00

1.50

2.00

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5

C3D8I

SC8R six elemts through thickness

SC8R refined in width direction

SC8R (Fig. 7b)

CPE4I (benchmark)

CPE4

CPS4

CPS4I

G

II ,

kJ/m2

y/B

Figure 26. Computed strain energy release rate distribution across the width of an ENF specimen.

46
0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

reltol=0.2; damv=0.1

reltol=0.2; damv=1E-3

reltol=0.2; damv=1E-4

reltol=0.2; damv=1E-5

reltol=0.2; damv=1E-6

reltol=0.1; damv=1E-6

reltol=0.1; damv=1E-6 (CPS4I)

reltol=0.01; damv=1E-6 (CPS4I)

load 

Q, N

displacement w, mm

Figure 29. Computed critical load-displacement behavior for ENF specimen obtained 

from two-dimensional planar models with added viscous regularization.

wcrit=1.42 mm

Qcrit

=1526 N

benchmark

0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

reltol=0.2; cs=1E-6

reltol=0.1; cs=1E-6

reltol=0.01; cs=1E-6

reltol=0.01; cs=1E-4

reltol=0.01; cs=1E-7

reltol=0.01; cs=1E-8

load 

Q, N

displacement w, mm

Figure 28. Computed critical load-displacement behavior for ENF specimen obtained 

from two-dimensional planar models (CPE4I) with added contact stabilization.

wcrit=1.42 mm

Qcrit

=1526 N

benchmark

47
0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

benchmark

reltol=0.5; cs=1E-6

reltol=0.2; cs=1E-6

reltol=0.1; cs=1E-6

reltol=0.1; cs=1E-6, Fig.7a

reltol=0.1; cs=1E-6, Δa=0.5 mm

load 

Q, N

displacement w, mm

Figure 31. Computed critical load-displacement behavior for ENF specimen

 obtained from results using three-dimensional solid finite element models.

wcrit=1.42 mm

Qcrit

=1526 N

0 5 10 15 20 25 30 35 40 45 50 55

0

1

2

3

4

5

critical

reltol=0.5

reltol=0.2; damv=1E-4

reltol=0.2; damv=1E-5

reltol=0.2; damv=1E-6

reltol=0.1; damv=1E-6

reltol=0.2; cs=1E-6

reltol=0.01; cs=1E-6

reltol=0.01; cs=1E-8

reltol=0.1; damv=1E-6 (CPS4I)

reltol=0.01; damv=1E-6 (CPS4I)

displacement

w, mm

Figure 30. Critical displacement-propagation behavior for an ENF specimen.

increase in delamination length a*, mm

benchmark

48
0 5 10 15 20 25 30 35 40 45 50 55

0

1

2

3

4

5

critical

reltol=0.5; cs=1E-6

reltol=0.2; cs=1E-6

reltol=0.1; cs=1E-6

SC8R, model Fig. 7b

SC8R,width refined

displacement

w, mm

Figure 33. Computed critical displacement-propagation for ENF specimen

 obtained from results using three-dimensional finite element models.

increase in delamination length a*, mm

benchmark

0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

SC8R, model Fig. 7b

SC8R , width refined

load 

Q, N

displacement w, mm

Figure 32. Computed critical load-displacement behavior for ENF specimen

 obtained from results using continuum shell finite element models.

wcrit=1.42 mm

Qcrit

=1526 N

benchmark

49
0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

reltol=0.2; damv=0.1

reltol=0.1; damv=0.1

benchmark

reltol=0.2; damv=0.1 (CPS4I)

load 

Q, N

displacement w, mm

Figure 35. Computed critical load-displacement behavior for ENF specimen obtained 

from two-dimensional planar models with added viscous regularization.

wcrit=1.42 mm

Qcrit

=1526 N

0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

benchmark

reltol=0.01; cs=1E-4

reltol=0.01; cs=1E-6

reltol=0.2; cs=1E-6 (CPS4I)

reltol=0.1; cs=1E-6 (CPS4I)

reltol=0.01; cs=1E-6 (CPS4I)

load 

Q, N

displacement w, mm

Figure 34. Computed critical load-displacement behavior for ENF specimen obtained 

from two-dimensional planar models with added contact stabilization.

wcrit=1.42 mm

Qcrit

=1526 N

50
0.0 1.0 2.0 3.0 4.0 5.0 6.0

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

reltol=0.1; cs=1E-6

load

Q, N

displacement w, mm

Figure 37. Computed critical load-displacement behavior for an ENF specimen 

obtained from three-dimensional models subjected to an applied center load Q.

wcrit=1.42 mm

Qcrit

=1526 N

benchmark

0 5 10 15 20 25 30 35 40 45 50 55

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

reltol=0.2; damv=0.1

reltol=0.1; damv=0.1

reltol=0.01; cs=1E-4

reltol=0.01; cs=1E-6

reltol=0.2; damv=0.1; CPS4I

reltol=0.01; cs=1E-6; CPS4I

Figure 36. Computed critical load-propagation for ENF specimen obtained from

results using two-dimensional planar finite element models for applied center deflection.

load 

Q, N

increase in delamination length a*, mm

benchmark

51
0 5 10 15 20 25 30 35 40 45 50 55

0

200

400

600

800

1000

1200

1400

1600

1800

2000

2200

critical

reltol=0.1; cs=1E-6

Figure 38. Computed critical load-propagation for ENF specimen obtained from 

results using three-dimensional solid finite element models subjected to an applied center load, Q.

load 

Q, N

increase in delamination length a*, mm

benchmark

52
PTFE film

Figure 40. Computed delamination front shape for ENF specimen (FE model in Figure 6).

static

pre-crack

B

growth under 

cyclic loading

a. Initial delamination front shape (Bottom surface of FE model in Figure 6).

location of initial delamination front

delaminated section

intact section

Figure 39. Fracture surface of typical ENF fatigue specimen [6].

b. Detail of delamination initiation at c. Detail of delamination propagation. 

edges.

location of initial 

delamination

front

edge 1

edge 2

53
a. Initial delamination front shape (Bottom surface of FE model).

Figure 42. Computed delamination front shape for ENF specimen (FE model in Figure 7a).

location of initial delamination front

delaminated section

intact section

Figure 41. Computed delamination front shape for ENF specimen (FE model with Δa=0.5 mm).

c. Detail of delamination propagation. 

a. Detail of delamination initiation at 

edges.

b. Detail of delamination propagation. 

edge 1

edge 2

b. Detail of delamination initiation at 

edges.

54
0.00 0.05 0.10 0.15 0.20 0.25 0.30

0

0.1

0.2

0.3

0.4

0.5

2D plane strain model (CPE4I)-Fig.5

2D plane stress model (CPS4I)-Fig.5

3D solid model (C3D8I)-Fig.6

time t, s

Figure 44. Computed energy release rate for ENF specimen.

R=0.1

f= 5 Hz

GII ,

 kJ/m2

 (N/mm)

Gmax,50 = 0.5 G

IIc =

 0.392 kJ/m2

G

min,50

 = 0.00392 kJ/m2

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0

0.2

0.4

0.6

0.8

1

1.2

1.4

desired fatigue loading from Fig.16

2D plane strain model (CPE4I)-Fig.5

2D plane stress model (CPS4I)-Fig.5

3D solid model (C3D8I)-Fig.6

time t, s

applied 

displacement 

w, mm

t

s

Figure 43. Computed cyclic fatigue loading for ENF specimen.

w=[A0+B1 sinω(t-t

0

)]*wmax

ω=10π=31.42

A

0

=0.55

B

1

=0.45

t

0

=-0.05

w50,max=1.0

w

50,min

=0.1

R=0.1

f= 5 Hz t

s

= =0.2 s 1

f

_

55
0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108 109

benchmark w

50

Δa=1.0 mm 

Δa=0.5 mm 

Δa=0.25 mm 

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 46. Computed delamination onset and growth obtained for different element length, Δa.

N

G

(delamination growth)

N

D

(delamination 

onset)

cutoff

Gth = 0.05

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108 109

benchmark w

50

t

0

=-0.05

t

0

=0.05

t

0

=0.0

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 45. Computed delamination onset and growth obtained for different values of t0.

N

G

(delamination growth)

N

D

(delamination 

onset)

cutoff

Gth = 0.08

56
0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108 109

benchmark w

50

center, one Fourier term

edge 2, one Fourier term

edge 1, one Fourier term

center, 50 Fourier terms

edge 1, 50 Fourier terms

edge 2, 50 Fourier terms

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 48. Computed delamination onset and growth obtained from full three-dimensional models 

(Δa=1.0 mm, Fig. 6).

N

G

(delamination growth)

N

D

(delamination 

onset)

cutoff

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108 109

benchmark w

40

benchmark w

30

benchmark w

20

computed growth

for w

40

computed growth

for w

30

computed growth

for w

20

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 47. Computed delamination onset and growth for different applied cyclic displacements.

N

G

(delamination growth)

N

D

(delamination onset)

cutoff

57
0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108

benchmark Q

50

t

0

=-0.05

t

0

=0.05

t

0

=0.0

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 50. Computed delamination onset and growth obtained for different values of t0.

N

G

N (delamination growth) D

(delamination 

onset)

plus refined width, Fig.7a

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108 109

benchmark w50

center

edge 1

edge 2

center

edge 1

edge 2

center

edge 1

edge 2

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 49. Computed delamination onset and growth obtained from full three-dimensional models.

N

G

N (delamination growth) D (delamination onset) cutoff

Δa=0.5 mm, 

one Fourier term

Δa=0.5 mm, 

50 Fourier terms

58
0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108

benchmark Q

40

benchmark Q

30

benchmark Q

20

computed growth

for Q

40

computed growth

for Q

30

computed growth

for Q

20

increase in 

delamination

 length 

a*, mm

number of cycles N

T

Figure 52. Computed delamination onset and growth for different applied cyclic loads.

Δa=1.0 mm

50 Fourier terms

t

s

=0

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108

benchmark Q

50

Δa=1.0 mm 

Δa=0.5 mm 

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 51. Computed delamination onset and growth obtained for different element lengths, Δa.

N

G

N (delamination growth) D

(delamination 

onset)

59
0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108

benchmark Q

50

center

edge 2

edge 1

center

edge 2

edge 1

center

edge 2

edge 1

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 54. Computed delamination onset and growth obtained from 

full three-dimensional models (Δa=0.5 mm).

Δa=0.5 mm, one Fourier term

Δa=0.5 mm, 50 Fourier terms

N

G

N (delamination growth) D

(delamination 

onset)

plus refined width, Fig.7a

0

5

10

15

20

25

30

35

40

45

50

55

100 101 102 103 104 105 106 107 108

benchmark Q

50

center, one Fourier term

edge 1, one Fourier term

edge 2, one Fourier term

center, 50 Fourier terms

edge 1, 50 Fourier terms

edge 2, 50 Fourier terms

increase in

delamination

 length 

a*, mm

number of cycles N

T

Figure 53. Computed delamination onset and growth obtained from 

full three-dimensional models (Δa=1.0 mm, Fig. 6).

N

G

N (delamination growth) D

(delamination 

onset)

60
location of initial delamination front

delaminated section

intact section

a. Initial delamination front shape (Bottom surface of FE model Δa=1.0 mm).

c. Detail of delamination growth. 

edge 1

edge 2

b. Detail of delamination initiation at edges.

d. Detail of delamination growth. e. Detail of final front shape. 

Figure 55. Computed delamination front shape for ENF specimen with one Fourier term

(FE model in Figure 6).

NT =1.0 108

NT = 400

NT = 5500

NT = 740,000

61
location of initial delamination front

delaminated section

intact section

b. Detail of delamination initiation at edges.

d. Detail of delamination growth. e. Detail of final front shape.

Figure 56. Computed delamination front shape for ENF specimen with 50 Fourier terms 

(FE model in Figure 6).

c. Detail of delamination growth. 

a. Initial delamination front shape (Bottom surface of FE model Δa=1.0 mm).

edge 1

edge 2

NT = 1.0 108

NT = 400 NT = 5300

NT = 9590

62
a. Detail of delamination growth. b. Detail of final front shape.

Figure 57. Computed delamination front shape for ENF specimen with one Fourier term

(Δa=0.5 mm).

a. Detail of delamination growth. b. Detail of final front shape. 

a. Detail of delamination growth. 

Figure 59. Computed delamination front shape for ENF specimen (Δa=0.5 mm, model Fig. 7a).

b. Detail of final front shape. 

Figure 58. Computed delamination front shape for ENF specimen with 50 Fourier terms

(Δa=0.5 mm).

NT = 3.5 106

NT = 2377

NT = 1.0 108

NT = 1.0 108 NT = 16000

NT = 8960

63
location of initial delamination front

delaminated section

intact section

a. Initial delamination front shape (Bottom surface of FE model Δa=1.0 mm).

d. Detail of delamination growth. e. Detail of final front shape. 

Figure 60. Computed delamination front shape for ENF specimen (FE model in Figure 6).

b. Detail of first delamination 

initiation at edges.

c. Detail of delamination growth.

edge 1

edge 2

NT = 6000

NT = 329 NT = 1482

NT = 1504

64
b. Detail of final front shape. 

Figure 61. Computed delamination front shape for ENF specimen with one Fourier term 

(Δa=0.5 mm).

a. Detail of delamination growth. b. Detail of final front shape. 

a. Detail of delamination intitiation at edges. 

Figure 63. Computed delamination front shape for ENF specimen (Δa=0.5 mm, model Fig. 7a).

b. Detail of final front shape. 

Figure 62. Computed delamination front shape for ENF specimen with 50 Fourier terms 

(Δa=0.5 mm).

NT = 6000

NT = 6000

NT = 9540

a. Detail of delamination growth. 

NT = 144

NT = 642

NT = 712

65
initial delamination front bonded nodes

X

Y Z

Figure A1. Deformed 2D finite element model of an ENF specimen.

w

DETAIL

2h

bonded nodes

(BONDED)

node set

(LFRONT)

initial delamination front

 bottom surface

(VCCT_BOT)

top surface

(VCCT_TOP)

w=0 u,w=0

2,u

1,w

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

0 0.2 0.4 0.6 0.8 1

experimental data [9]

mean values [9]

experimental data [6]

mean values [6]

DCB, Mode I MMB, Mode I and II

Gc ,

kJ/m2

mixed mode ratio GII /GT

Figure A2. Mixed mode fracture criterion.

ENF, Mode II

G

Ic

 =0.21

 GIIc = 0.78

curve fit: G η = 2.57 c = G

Ic

 + (GIIc-GIc)(GII/GT)

η

*FRACTURE CRITERION, TYPE=fatigue, MIXED MODE BEHAVIOR=BK, 

TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

66
10-8

10-7

10-6

10-5

10-4

10-3

10-2

10-1

100

10-2 10-1 100

da/dN, 

mm/cycle

G

II,max, kJ/m2

Figure A4. Delamination growth rate (Paris Law) for IM7/8552.

GIIc = 0.78

Gth = 0.08

Paris Law [6]

da/dN = c

3

 Gc4

c

3

=0.33185

c

4

=5.5519

r

2

=Gpl /GIc = 0.9

r

1

=Gth/GIc= 0.102

*FRACTURE CRITERION,TYPE=fatigue,MIXED MODE BEHAVIOR=BK,TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

Gpl = 0.705

0.00

0.20

0.40

0.60

0.80

1.00

100 101 102 103 104 105 106 107

experimental data [6]

N

D , cycle

Gmax ,

 kJ/m2 (N/mm)

Figure A3. Delamination growth onset for mode II for IM7/8552 [6].

curve fit: G= m

0

 * N 

D

m

0

= 0.78

m

1

= -0.16

m1 N

D

 = c

1

*G c2

c

1 0.213

c

2 -6.329

*FRACTURE CRITERION, TYPE=fatigue, MIXED MODE BEHAVIOR=BK, 

TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

67
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

0

0.2

0.4

0.6

0.8

1

1.2

1.4

time t, s

applied 

displacement 

w, mm

Figure A6. Cyclic fatigue loading input.

w50,max=1.0

w

50,min

=0.1

w=[A0+B1 sinω(t-t0)]*1.0

*DIRECT CYCLIC, fatigue

0.001,0.2,,,50,50,,20

,,10000000,,

*BOUNDARY,AMPLITUDE=test

LFRONTM, 1, 1, -1.0 t

s

R=0.1

f= 5 Hz t

s

= =0.2 s 1

f

_

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0

0.2

0.4

0.6

0.8

1

1.2

1.4

time t, s

amplitude a

Figure A5. Amplitude curve (relative load magnitude).

ω=10π=31.42

A

0

=0.55

A1

=0

B

1

=0.45

t

0

=-0.05

t

0

= 0.0

t

0

= 0.05

R=0.1

f= 5 Hz

*AMPLITUDE, NAME=test, DEFINITION=PERIODIC

1,31.42,-0.05,0.55

0,0.45

a=A

0

+A1 cosω(t-t0)+B1 sinω(t-t0)

68
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Contractor Report

 4. TITLE AND SUBTITLE

Development and Application of Benchmark Examples for Mode II Static 

Delamination Propagation and Fatigue Growth Predictions 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Krueger, Ronald

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

NIA Report No. 2011-02

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

Langley Technical Monitor: Jonathan B. Ransom

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 24

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

The development of benchmark examples for static delamination propagation and cyclic delamination onset and growth prediction is presented and demonstrated for a commercial 

code. The example is based on a finite element model of an End-Notched Flexure (ENF) specimen. The example is independent of the analysis software used and allows the 

assessment of the automated delamination propagation, onset and growth prediction capabilities in commercial finite element codes based on the virtual crack closure technique 

(VCCT). First, static benchmark examples were created for the specimen. Second, based on the static results, benchmark examples for cyclic delamination growth were created. 

Third, the load-displacement relationship from a propagation analysis and the benchmark results were compared, and good agreement could be achieved by selecting the appropriate 

input parameters. Fourth, starting from an initially straight front, the delamination was allowed to grow under cyclic loading. The number of cycles to delamination onset and the 

number of cycles during stable delamination growth for each growth increment were obtained from the automated analysis and compared to the benchmark examples. Again, good 

agreement between the results obtained from the growth analysis and the benchmark results could be achieved by selecting the appropriate input parameters. The benchmarking 

procedure proved valuable by highlighting the issues associated with choosing the input parameters of the particular implementation. Selecting the appropriate input parameters, 

however, was not straightforward and often required an iterative procedure. Overall the results are encouraging, but further assessment for mixed-mode delamination is required. 

15. SUBJECT TERMS

benchmarking; composites; delamination; finite element analysis; growth prediction; virtual crack closure technique

18. NUMBER

 OF 

 PAGES

73

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

NNL09AA00A

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

877868.02.07.07.05.01.01 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/CR-2011-217305

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 11 - 2011