July 2010 

NASA/CR-2010-216723

NIA Report No. 2010-04 

Development of a Benchmark Example for 

Delamination Fatigue Growth Prediction 

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

 

July 2010 

NASA/CR-2010-216723

NIA Report No. 2010-04 

Development of a Benchmark Example for 

Delamination Fatigue Growth Prediction 

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 OF A BENCHMARK EXAMPLE FOR DELAMINATION FATIGUE 

GROWTH PREDICTION 

Ronald Krueger*

ABSTRACT 

The development of a benchmark example for cyclic delamination growth prediction is 

presented and demonstrated for a commercial code. The example is based on a finite 

element model of a Double Cantilever Beam (DCB) specimen, which is independent of the 

analysis software used and allows the assessment of the delamination growth prediction 

capabilities in commercial finite element codes. First, the benchmark result was created for 

the specimen: The number of cycles to delamination onset was calculated from the material 

data for mode I fatigue delamination growth onset. Then, the number of cycles during stable 

delamination growth was obtained incrementally from the material data for mode I fatigue 

delamination propagation. For the benchmark case, where the results for delamination onset 

and growth were combined, 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 in a finite element model of a commercial code. 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. In general, good 

agreement between the results obtained from the growth analysis and the benchmark results 

could be achieved by selecting the appropriate input parameters. 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. 

An approach for assessing the delamination propagation capabilities in commercial finite 

element codes under static loading was recently presented and demonstrated for VCCT for 

ABAQUS®1 [3, 4]. In these recent papers, benchmark results were created for full three-dimensional 

finite element models of the Double Cantilever Beam (DCB) and the Single Leg Bending (SLB) 

specimen. Then, starting from an initially straight front, the delamination was allowed to propagate 

in the finite element model. The load-displacement relationship and the total strain energy obtained 

from the propagation analysis results and the benchmark results were compared and good 

 

*

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

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

1
agreements could be achieved by selecting the appropriate input parameters. Overall, the results 

were encouraging but it was determined that further assessment for mixed-mode delamination is 

required [3, 4]. 

The objective of the present study was to create a benchmark example, independent of the 

analysis software used, which allows the assessment of the delamination fatigue growth prediction 

capabilities in commercial finite element codes. At the beginning, a benchmark example is created 

based on incremental finite element models of a DCB specimen. To avoid unnecessary 

complications, experimental anomalies such as fiber bridging were not addressed. First, a sample 

material and cyclic loading was selected. Second, the number of cycles to delamination onset, ND, 

was calculated from the mode I fatigue delamination growth onset data of the material. Third, the 

number of cycles during stable delamination growth, NG, was obtained incrementally from the 

material data for mode I fatigue delamination propagation by using growth increments of a=0.1 

mm. Fourth, the total number of growth cycles, NG, was calculated by summing over the increments 

NG. Fifth, 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, the approach was demonstrated 

for the commercial finite element code ABAQUS®. Starting from an initially straight front, the 

delamination was allowed to grow based on the algorithms implemented into the commercial finite 

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

delamination increase during cyclic loading. It was assumed that delamination length increase 

during cyclic loading obtained from finite element analysis should closely match the growth shown 

in the benchmark example. The benchmark enabled the selection of the appropriate 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. 

2. METHODOLOGY 

The methodology for delamination propagation, onset and growth was applied to the DCB 

specimen to create the benchmark example [5, 6]. Since the required material property input data 

were not readily available in the open literature for a single material, a fictitious set of properties 

was constructed for this benchmarking exercise. Individual properties for commonly used 

graphite/epoxy tape materials were obtained from the open literature to create this set to represent a 

typical graphite epoxy composite. The material properties are given in Tables I and II. 

2.1 Static fracture toughness 

The mode I fracture toughness (mixed-mode ratio GII/GT=0) is generated experimentally using 

the Double Cantilever Beam (DCB) tests (as shown in Figure 1) [7]. A fracture toughness GIc= 

0.17 kJ/m2 was used in this benchmarking exercise [8]. 

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 2 [9, 10]. The onset curve (solid green line) is a power law fit 

G = m0 ND

m1 (1) 

2
of the experimental data (open, green circles) obtained from a DCB test using the respective 

standard for delamination growth onset [9]. 

2.3 Fatigue delamination growth 

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

delamination propagation relationship (Paris Law) plotted in Figure 3 [8]. The delamination growth 

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

da

dN = c Gmax

n (2) 

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 are obtained by fitting the curve 

to the experimental data (open black circles) obtained from DCB tests [8]. The critical energy 

release rate or fracture toughness, GIc, was included in the plot of Figure 3 (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) [8]. It has to be noted that this benchmarking exercise ignores branching or fiber bridging and 

hence the Paris Law was not normalized with the static R-curve as recently suggested [11, 12]. 

3. SPECIMEN AND FATIGUE TEST DESCRIPTION 

For the current numerical investigation, the Double Cantilever Beam (DCB) specimen, as 

shown in Figure 1, was chosen since it is simple, only exhibits the mode I 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 [3, 4]. To avoid unnecessary 

complications, experimental anomalies such as fiber bridging were not addressed. For the current 

study, a DCB specimen made of graphite/epoxy with an unidirectional layup, [0]24, was modeled. 

The material properties are given in Table I [8]. The material, layup, overall specimen dimensions 

including initial crack length, a, were identical to the specimen used earlier [3, 4]. 

For the cyclic loading of the specimen, guidance was taken from a draft standard designed to 

determine mode I fatigue delamination propagation [11]. In the draft document, it is recommended 

to start the test at a maximum displacement, max, which causes the energy release rate at the front, 

GImax, to reach initially about 80% of GIc

GImax

GIc

= 0.8 (3) 

The maximum load, Pmax, and maximum displacement, max/2, were calculated using the known 

quadratic relationship between energy release rate and applied load or displacement 

GImax

GIc

= Pmax

2

Pcrit

2 Pmax = Pcrit

GImax

GIc

 , max = crit

GImax

GIc

 (4) 

3
 max = crit 0.8 (5) 

where Pcrit and crit are the critical values. For the current study, a critical energy release rate 

GIc=0.17 kJ/m2

 was used and the critical values Pcrit and crit (grey dashed lines) were obtained from 

the benchmark for static delamination propagation [3, 4] shown in the load-displacement plot in 

Figure 4. The calculated maximum load, Pmax, and calculated maximum displacement, max/2, are 

shown in Figure 4 (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 a DCB 

specimen under displacement control, the applied maximum displacement, max/2=0.67 mm, is kept 

constant while the load drops as the delamination length increases (solid red circles and solid red 

line). The energy release rate corresponding to an applied maximum displacement max/2=0.67 mm 

was calculated for different delamination lengths a using equation (5). The energy release rate 

decreases with increasing delamination length, a, as shown in Figure 5 (solid red circles and solid 

red line). Delamination growth was assumed to stop once the calculated energy release rate drops 

below the cutoff value, Gth, (green solid horizontal line). The static benchmark case (solid grey 

circles and dashed grey line in Figure 5), where the delamination propagates at constant GIc (solid 

blue line in Figure 5) was included for comparison. 

In the ASTM draft document, it is suggested to use a load ratio R=0.1 for testing [11]. The 

corresponding minimum load, Pmin, and minimum displacement, min/2, were calculated 

 

R = Pmin

Pmax

= min

max

= 0.1  min = 0.1max (6) 

Further, it was suggested in the ASTM draft document to use a frequency f=10 Hz for testing [11]. 

A graphical representation of the cyclic fatigue loading is plotted in Figure 6. The applied 

displacement /2 is represented as a function of time, t

 /2 = a0 + b1 sin t t [ ] ( )0 max /2 (7) 

where max/2=0.67 mm is the maximum displacement. The constants a0=0.55, b1=0.45, the circular 

frequency =20=62.832 and the starting time t0=0.025 are calculated from load ratio R=0.1 and 

the frequency f=10 Hz for testing. The resulting equation to calculate the applied displacement /2 

is shown in Figure 6. 

4. CREATING A BENCHMARK EXAMPLE FOR GROWTH PREDICTION 

4.1 Fatigue delamination growth onset 

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

ND. 

 

G = m0 ND

m1 ND = 1

m0





 



	 
 

1

m1

c1

  

G

1

m1

 ND = c1 Gc2 (8) 

4
 

where c2 = 1

m1

. Values for the constants c1 and c2 are shown in Figure 7. 

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

front, GImax, reaches about 80% of GIc corresponding to GImax=0.1362 kJ/m2

. The initial energy 

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

the delamination onset curve, the number of cycles to delamination onset is determined as, ND=150 

shown as a vertical dashed red line. 

4.2 Fatigue delamination growth 

The number of cycles during stable delamination growth can be obtained by solving equation 

(2) for NG

NG = dN = 1

c

Gmax

-n  da (9) 

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

GImax, reaches about 80% of GIc corresponding to GImax=0.1362 kJ/m2

 in the current study as shown 

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

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

a

N = c Gmax

n (10) 

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

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

from the curve fit (solid red circles and solid red line) plotted in Figure 5. 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 8. The growth rate rapidly decreases with increasing delamination 

length a as shown in Figure 9 (solid red circles and solid red line). The number of cycles during 

stable delamination growth, NG, was calculated by summing the increments Ni

 NG = Ni

i=1

k

 = 1

c

Gi,max

-n  a

i=1

k

 (11) 

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

i=1

k

 = a0 + k  a (12) 

For stable delamination growth, the delamination length, a, is plotted in Figure 10 for an 

increasing number of load cycles NG (crosses and solid red line). 

5
4.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 (13) 

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

delamination growth [12]. For this combined case, the delamination length, a, is plotted in Figure 11 

for an increasing number of load cycles NT. For the first ND cycles, the delamination length remains 

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

delamination length increases following the Paris Law (crosses and solid red line). Once a 

delamination length is reached where the energy release rate drops below the assumed cutoff value, 

Gth, (as shown in Figure 5) the delamination growth no longer follows the Paris Law (dashed grey 

line) and stops (horizontal solid red line). 

Applying the relationship 

 t = N

f

, (14) 

where t is the time and f is the frequency, the development of the delamination length a can be 

plotted on a time scale assuming a frequency of 10 Hz as shown in Figure 12 (solid red line). 

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

The load/displacement behavior of the DCB specimen - as shown in Figure 4 (solid red line) - 

can serve as an initial check for the finite element model. The correct input of model dimensions, 

material, layup and load application can thus be verified. During initial loading, the load and 

displacement should increase in a linear fashion and follow the initial slope until the maximum load, 

Pmax, and maximum displacement, max, are reached. To minimize problems with numerical stability 

of the analysis, it is suggested that prescribed displacements, max/2 and min/2 are applied in the 

analysis instead of nodal point loads, Pmax and Pmin. The same approach was used to create the static 

benchmark case mentioned above [3, 4]. Once the delamination starts to grow, the load is expected 

to drop while the applied maximum displacement max/2=0.67 is expected to remain constant. 

The energy release rate may serve as an additional 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 3D model should reach but not exceed the target value of GImax=0.136 kJ/m2

. 

Once delamination growth starts in the model, the computed energy release rate should decrease 

with increasing delamination length, a, as shown in Figure 5. The curve fit (solid red line) can 

therefore be used to check the computed energy release rate during delamination growth. 

The growth rate a/N, shown in Figure 9, decreases and the curve fit (solid red line) 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 delamination length, a, should increase with the number 

of cycles, N, as shown in Figure 10. The curve fit (solid red line) can therefore be used as a 

benchmark. 

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

as stable growth should yield results that closely resemble the plot in Figure 11. The curve fit (solid 

red line) can therefore be used as a benchmark. 

6
5. FINITE ELEMENT MODELING 

A typical two-dimensional finite element model of a Double Cantilever Beam (DCB) specimen 

is shown in Figure 13. The specimen was modeled with solid plane strain elements (CPE4) and 

solid plane stress elements (CPS4) in ABAQUS® Standard 6.8, 6.9 and 6.9EF. Along the length, all 

models were divided into different sections with different mesh refinement. The DCB specimen was 

modeled with six elements through the specimen thickness (2h) as shown in the detail of Figure 

13b. The resulting element length at the delamination tip was a=0.5 mm. A finer mesh, resulting in 

a=0.25 mm, was also generated as is shown in Figure 13c. Additionally, three coarser meshes with 

a reduced number of elements in the length direction were also generated as shown in Figures 14, 

resulting in a=1.0 mm (Figure 14b), a=1.25 mm (Figure 14c) and a=1.67 mm (Figure 14a). 

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

specimen. For the analysis with ABAQUS® 6.8, 6.9 and 6.9EF, 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 13b. Additionally, a node set was 

created to define the intact (bonded nodes) region. 

Typical three-dimensional finite element models of the DCB specimen are shown in Figures 15 

and 16. Along the length, all models were divided into different sections with different mesh 

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

Figure 15b. Across the width, a uniform mesh was used to avoid potential problems at the transition 

between a coarse and finer mesh [3,4]. Through the specimen thickness (2h), six elements were 

used as shown in the detail of Figure 15b. The resulting element length at the delamination tip was 

a=0.5 mm. The specimen was modeled with solid brick elements (C3D8I) which had yielded 

excellent results in a previous studies [3,4]. Two coarser meshes with a reduced number of elements 

in the width and length directions, resulting in a=1.0 mm and a=2.0 mm, were also generated as 

shown in Figures 16a and b. 

Three models of the DCB specimen were generated with continuum shell elements (SC8R) as 

shown in Figures 17a to c. 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 

4b, 5a and 5b, resulting in an element length at the delamination tip a=0.5 mm, a=1.0 mm and 

a=2.0 mm, respectively. In the z-direction, only one element was used to model the thickness of 

the specimen. 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 

models discussed above. 

For all the analyses performed, the low-cycle fatigue analysis in ABAQUS® Standard 6.8, 6.9 

and 6.9EF 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 

7
determine propagation, computed energy release rates are compared to the input data for onset and 

growth from experiments as discussed in the methodology section. During the analysis, at least one 

element length at the crack tip is released along the interface after each stabilized cycle. 

For all analyses, the elastic constants, the input to define the fracture criterion, and the 

parameters for delamination onset and delamination growth (Paris Law) were kept constant. The 

elastic constants are given in Table I. The fracture toughness values and the parameters required for 

delamination onset and growth are given in Table II. The parameters to define the load frequency 

(f=10 Hz), the load ratio (R=0.1) as well as the minimum and maximum applied displacement 

(min/2=0.067 mm and max/2=0.67 mm) were also kept constant during all analyses. To study the 

effect on the computed onset and growth behavior during the analysis: 

• The number of terms used to define a Fourier series was varied. A Fourier series is 

used during the execution of the ABAQUS® Standard to approximate the periodic 

cyclic loading. 

• The size of the initial time increment used in the analysis was varied. 

• The input required to define the cyclic loading was altered. 

• The release tolerance was varied. Once a user specified release tolerance 

((G Gc )/Gc > release tolerance) is exceeded during the analysis in ABAQUS®

Standard, a cutback operation is performed which reduces the time increment. The 

cutback reduces the degree of overshoot and improves the accuracy of the local 

solution. 

• The solution controls were varied. 

It was assumed that the computed onset and growth behavior should closely match the 

benchmark results established earlier. Setting the value of the input parameters correctly is often 

an iterative procedure, which will be discussed later. Further details about the required input 

parameters are discussed in the appendix where sample input files are also provided. 

6. ANALYSIS 

6.1 Static analysis to verify correct input data and model response 

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

correct and that the models responded as expected. For this static analysis, the applied displacement 

at the tip of the arms was ramped up to the maximum displacement max/2=0.67 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 18 and 19 was evaluated. In 

comparison to the fatigue benchmark (solid grey circles and solid grey line) from Figure 4, which is 

based on analyses using solid 3D volume elements (C3D8I), the results are in good agreement. The 

model where plane strain elements (CPE4) were used exhibits a slightly more compliant behavior 

(dashed blue line) as shown in Figure 18. The model where plane stress elements (CPS4) were used 

exhibits a slightly stiffer behavior (dash-dot green line) whereas the results from the solid model 

(C3D8I) are as expected identical (solid black line) to the benchmark results. The results obtained 

from solid models with different mesh densities (see Figures 15 and 16) and from the continuum 

shell models (see Figure 17) are in good agreement with the benchmark as shown in Figure 19.

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. 

8
In a second analysis step, only a single load cycle was analyzed, starting at the previously 

applied maximum displacement max/2=0.67. 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, /2, - obtained as analysis output - was plotted versus the step time, t, 

as shown in Figure 20. The models where plane strain elements (CPE4) were used (blue crosses) 

yielded the same output as the models where plane stress elements (CPS4) were used (green x’s) 

and solid model (C3D8I) were used (open black circles). 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 6. 

Additionally, the computed mode I energy release rate – also obtained as analysis output - was 

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

between the expected maximum value GImax=0.136 kJ/m2

 and minimum value GImin=0.00136 kJ/m2. 

The model where plane strain elements (CPE4) were used yields lower energy release rates (blue 

crosses and solid blue line) while the model where plane stress elements (CPS4) were used yields 

higher energy release rates (green x’s and solid green line) which is consistent with the observations 

made above with respect to the slight variation in model stiffness. The results from the solid model 

(C3D8I), taken in the center of the specimen at y=0.0, are as expected identical to the benchmark 

results (open black circles and solid black line). 

For the models made of 3D continuum elements (C3D8I and SC8R), the computed mode I 

strain energy release rate values were also computed for an applied maximum displacement 

max/2=0.67 and plotted versus the normalized width, y/B, of the specimen as shown in Figure 22. 

The results were obtained from models shown in Figures 15 through 17. Qualitatively, the mode I 

strain energy release rate is fairly constant in the center part of the specimen and drops progressively 

towards the edges. As desired, the energy release rate in the center of specimen reached the 

expected maximum value GImax=0.136 kJ/m2

 for all models used. Additionally, these results served 

as a verification that the VCCT procedure implemented in ABAQUS® 6.8 yielded the same results 

(open symbols) as an external post-processing routine (crosses) [4]. As expected, the mode II and 

mode III strain energy release rates were computed to be nearly zero and hence are not shown. 

Based on results, it was assumed that with respect to the computed energy release rates all models 

adequately represented the benchmark case. 

6.2 Results from fatigue onset and growth analysis 

In Figures 23 to 32, the delamination length, a, is plotted versus the number of cycles, N, for 

different input parameters and models. For all results shown, the analysis stopped when a 

10,000,000 cycles limit - used as input to terminate the analysis - was reached. 

6.2.1 Initial results

Initial results, as plotted in Figures 23 and 24, were obtained using the specified default values 

as input parameters (see appendix). The results obtained from two-dimensional plane strain (open 

blue circles and solid blue line) and plane stress (open red squares and solid red line) models as well 

as full three-dimensional models (open green diamonds and solid green line) were within 4% of the 

benchmark results (grey crosses and solid grey line) as shown in Figure 23. For better visualization 

of the results and to be able to identify the differences in the results, the scale on the vertical axis 

was expanded as shown in Figure 24. For all results shown, the predicted onset occurs prior to the 

benchmark result onset, ND=150 cycles. The onset value is lowest for the plane stress results, 

9
followed by the results obtained from 3D solid models and the plane strain results. This sequence 

can be explained by the computed energy release rate shown in Figure 21, where the values 

obtained from the plane stress model are slightly higher compared to the results obtained from 3D 

solid and plane strain models. During the growth phase of the analysis, the results lie on curves with 

nearly the same slope parallel to the benchmark which suggests that the Paris Law was implemented 

correctly and is - as expected - independent of the model. For all models, 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. 

6.2.2 Variation of input parameters 

Input parameters were varied to study the effect on the computed onset and growth behavior 

during the analysis and results are plotted in Figures 25 to 27. For this parametric study, only 

models with a refined mesh (a=0.25 mm, see Figure 13c) made of plane strain elements (CPE4) 

were used. Once a set of parameters was established that yielded good results, the effects of mesh 

size and element type on the results were studied. 

First, the effect of the initial time increment used in the analysis settings was studied as 

shown in Figure 25 (see appendix for details). The initial time increment was varied between 

i0=0.01 (one tenth of a single loading cycle, ts=0.1s, open blue circles and solid blue line) and 

i0=0.0001 (one thousandth of a single loading cycle, open green diamonds and solid green line). 

For larger initial time increments, the onset of delamination shifted towards a lower number of 

cycles. Reducing the initial time increments, however, significantly increased the computation 

time. Based on the results, it was therefore decided to use an initial time increment of i0=0.001 

(open red squares and solid red line) for the remainder of the study to save computation time. This 

step is justified by the fact that the results obtained for i0=0.001 were almost identical to the values 

obtained from the analysis where a smaller initial time increment was used (i0=0.0001). 

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

displacement, /2, a set of parameters need to be determined as shown in the equation in Figures 26 

and 27 (see appendix for details). Selecting A1=0 results in a sine curve representation of the cyclic 

load (analysis results are shown in Figure 26). Selecting B1=0 results in a cosine curve 

representation of the cyclic load (analysis results are shown in Figure 27). The selection of the 

starting time, t0, causes a phase-shift. Further details about the input parameters are discussed in 

detail in the appendix where the corresponding cyclic applied displacements /2 are plotted in 

Figure A8. As shown in Figures 26 and 27, the number of cycles to delamination onset is the most 

sensitive to input variations. The results obtained for the growth phase however, lie on curves with 

nearly the same slope, parallel to the benchmark. Also 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. A narrow band of results (green solid lines), was 

obtained when the number of terms used to define the Fourier series was increased to 50. A 

Fourier series is used in ABAQUS® Standard during the analysis to approximate the periodic 

cyclic loading (see appendix). The narrow band of results was in good agreement with the 

benchmark curve (grey solid line) compared to the results obtained for the default setting of 11 

Fourier terms which showed more scatter (black dashed lines). It should be noted that one term in 

a Fourier series is sufficient to exactly represent a simple periodic function such as the simple sine 

function used in the current example. At the same time, a large number of terms should improve the 

approximation of a more complicated cyclic load as demonstrated. It should also be noted that good 

results could only be obtained when the analyses were performed with ABAQUS® Standard 6.9EF. 

10
Earlier versions showed unexplainable wide variations in results. Based on the current results, it was 

decided to use a sine curve representation (A1=0) in combination with the starting time, t0=0.0 for 

the remainder of the study. 

Third, the release tolerance was varied. Once a user specified release tolerance ((G Gc )/Gc >

release tolerance) is exceeded during the analysis in ABAQUS® Standard, a cutback operation is 

performed which reduces the time increment The cutback reduces the degree of overshoot and 

improves the accuracy of the local solution [13]. In the current study, varying the input between 0.2 

(the default value) and 0.01 did not have any effect on the computed onset and growth behavior 

during the analysis. It is assumed that the release tolerance only affects static delamination 

propagation and not growth under cyclic loading studied here. 

Fourth, the solution controls were varied. It is generally not required nor recommended to 

modify the solution controls in ABAQUS® Standard. Since some of the solution controls, however, 

were modified in the DCB example problem provided by ABAQUS® [15], it was decided to briefly 

address the effect on the computed onset and growth behavior during the analysis. The first input 

parameter that defines the iteration number at which the periodicity condition is first imposed, 

was not modified and the default value was kept. A series of analyses were performed where the 

other four input parameters were varied between 100 and 10-4. Changing the input did not have any 

effect on the computed onset and growth behavior during the analysis, however it significantly 

influenced the computation time. The analysis required only about a tenth of the computation time 

when the second input parameter was modified from the default value (510-3) to values larger than 

10-1. Changing the other remaining three parameters did not have any effect for this example. 

Further details about the input parameters are discussed in the appendix where also sample input 

files are provided. 

6.2.3 Variation of mesh size and element type 

The results obtained for models with different mesh sizes and different element types are shown 

in Figures 28 to 33. For models made of plane strain elements (CPE4), the element length, a, at the 

delamination tip was varied as shown in Figures 13 and 14. The results obtained from the respective 

models are plotted in Figures 28 and 29. Excellent agreement with the benchmark curve (grey 

crosses and grey solid line) could be achieved for element lengths up to a=1.25 mm, as shown in 

Figure 28. For element length a=1.67 mm (purple triangles and solid purple line), the predicted 

onset occurs for a slightly higher number of cycles. The observed mismatch is largely due to the 

increased element length, which causes the first growth step to be larger. The following growth 

increments are relatively large due to the increased element length (a=1.67 mm). However, during 

the growth phase of the analysis, the results from all models lie on curves with the same slope 

parallel to the benchmark, which suggests that the Paris Law was implemented correctly and is, as 

expected independent, of the mesh. The analyses were repeated with the goal to eliminate the onset 

part of the analysis and focus on the stable growth section. Input parameters were chosen so that the 

onset part was extremely short and stable growth basically started immediately. Further details 

about the input parameters are discussed in detail in the appendix. The results are shown in 

Figure 29. As mentioned above, excellent agreement with the benchmark curve (grey crosses and 

grey solid line) could be achieved for all element lengths except for element length a=1.67 mm 

(purple triangles and solid purple line). For all models, the threshold cutoff, where delamination 

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

11
number of cycles defined by the benchmark. The total computation time was between 70 s for the 

coarsest mesh and 1030 s for the finest mesh as shown in Figure 301. 

The analyses were repeated for models made of plane stress elements (CPS4), where the 

element length, a, at the delamination tip was varied as before. The results obtained from the 

respective models are plotted in Figure 31. As before for the plane strain models, excellent 

agreement with the benchmark curve (grey crosses and grey solid line) could be achieved for all 

element lengths except for element length a=1.67 mm (purple triangles and solid purple line). The 

total computation time took between 80 s for the coarsest mesh and 1050 s for the finest mesh1

. The 

computation times are included in Figure 30 for comparison. 

The results obtained for models made of 3D solid and continuum shell elements (shown in 

Figures 15 to 17) are shown in Figure 32. Good agreement with the benchmark curve (grey crosses 

and grey solid line) could be achieved for small element lengths (a=0.5 mm). For larger element 

lengths, the results started to deviate from the benchmark. However, during the growth phase of the 

analysis, the results from all models lie on curves with the same slope parallel to the benchmark as 

observed before. Results obtained from the continuum shell element models (blue dashed lines) 

were almost identical when compared to results obtained from corresponding full 3D solid models 

(green solid line) with the same element length, a. The computational effort, however, was 

reduced by a factor of 2.5 to 2.9 for analyses performed using the continuum shell element models. 

For the models made of 3D solid elements, the total computation time was between 5800 s for the 

coarsest mesh and about 7 days for the finest mesh. For the corresponding models made of 

continuum shell elements, the total computation time was between 2000 s for the coarsest mesh and 

about 2.8 days for the finest mesh1

. The computation times are also included in Figure 30 for 

comparison. 

For final comparison, finite element analyses were repeated with two-dimensional and threedimensional models with the same element length, a=0.5 mm, at the delamination tip. The results 

obtained from two-dimensional plane strain (open blue circles and solid blue line) and plane stress 

(open red squares and solid red line) models as well as full 3D solid models (open green diamonds 

and solid green line) and continuum shell elements (orange x’s and solid orange line) were within 

1% of the benchmark curve (grey crosses and solid grey line) as shown in Figure 33. The results 

obtained from the continuum shell element model (orange x’s and solid orange line) were almost 

identical compared to results obtained from the full 3D solid model (open green diamonds and solid 

green line) as mentioned above. The results obtained from plane stress models (open red squares 

and solid red line) are close to the results obtained from three-dimensional models. For all results 

shown, the predicted onset occurs first for the plane stress results and last for the plane strain results. 

Also, the delamination length obtained from plane strain models are slightly lower. This observation 

can be explained by looking at the computed energy release rate in Figure 21. The plane stress 

model yields a slightly higher energy release rate compared to the 3D solid and the plane strain 

model. Delamination onset would therefore occur first in plane stress models as indicated by the 

results plotted in Figure 33. During the stable growth phase, however, the results for all models lie 

on curves with nearly the same slope parallel to the benchmark as mentioned earlier. For all models, 

the threshold cutoff, where delamination growth is assumed to stop and the delamination length 

remains constant, is predicted close to the number of cycles defined by the benchmark. 

 

1

 CPU time on Dual-Core AMD Opteron(tm) Processor 8220 SE 

12
7. SUMMARY AND CONCLUSIONS 

The development of a benchmark example for cyclic delamination growth prediction is 

presented and demonstrated for the commercial finite element code ABAQUS® Standard. The 

example is based on a finite element model of a Double Cantilever Beam (DCB) specimen, which is 

independent of the analysis software used and allows the assessment of the delamination growth 

prediction capabilities in commercial finite element codes. First, 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 I fatigue delamination 

growth onset. The number of cycles during stable delamination growth was obtained incrementally 

from the material data for mode I fatigue delamination propagation. For the combined benchmark 

case of delamination onset and 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. 

The results showed the following: 

• In general, good agreement between the results obtained from the 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. 

• The onset prediction appeared much more sensitive to the input parameters than the 

growth prediction. 

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

• Good agreement between analysis results and the benchmark could be achieved when 

the initial time increment used in the analysis was about one tenth of a single loading 

cycle. 

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

• The release tolerance did not have an effect on the analysis or the computed results. 

• Accurately computing the onset and growth required fine meshes with an element 

length at the tip a1.0 mm. 

• The solution controls in ABAQUS® Standard had to be modified in order to reduce 

computation time. Even with carefully selected input parameters, the analyses for threedimensional models of a simple DCB specimen required days. 

• Although implemented in ABAQUS® Standard 6.8, version 6.9EF was required to 

obtain consistently good results. 

• Improvements are needed to make this analysis applicable to real case scenarios such as 

more complex specimens or structural components. 

Overall, the results are promising. In a real case scenario, however, where the results are 

unknown, obtaining the right solution will remain challenging. Further studies are required which 

should include the assessment of the propagation capabilities in more complex mixed-mode 

specimens and on a structural level.

13
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 with confidence 

to model delamination growth for more complex configurations. 

ACKLOWLEDGEMENTS 

This research was supported by the Aircraft Aging and Durability Project as part of NASA’s 

Aviation Safety Program. 

The analyses were performed at the Durability, Damage Tolerance and Reliability Branch at 

NASA Langley Research Center, Hampton, Virginia, USA. 

The author would like to thank Satish Gadhi of SIMULIA Headquarters and Babu 

Aminjikarai and Brian Baillargeon of the SIMULIA Eastern Region office for their support and 

advice in getting this analysis started. 

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 for Assessing Delamination Propagation Capabilities in Commercial Finite 

Element Codes," in proceeding of the 22nd Annual Technical Conference of the American Society for 

Composites, Seattle, Washington, 2007. 

4. R. Krueger, "An Approach to Assess Delamination Propagation Simulation Capabilities in Commercial Finite 

Element Codes," NASA/TM-2008-215123, 2008. 

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

6. "Composite Fatigue Damage Onset," in Composite Materials Handbook CMH-17. Vol. 3, section 12.6.4, 

draft of revision G, 2009. 

7. "ASTM D 5528-01, Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional 

Fiber-Reinforced Polymer Matrix Composites," in Annual Book of ASTM Standards. vol. 15.03: American 

Society for Testing and Materials, 2008. 

8. M. König, R. Krüger, K. Kussmaul, M. v. Alberti, and M. Gädke, "Characterizing Static and Fatigue 

Interlaminar Fracture Behaviour of a First Generation Graphite/Epoxy Composite," in Composite Materials: 

Testing and Design - (13th Vol.), ASTM STP 1242, S. J. Hooper, Ed.: American Society for Testing and 

Materials, pp. 60-81, 1997. 

9. "ASTM D 6115-97, Standard Test Method for Mode I Fatigue Delamination Growth Onset of Unidirectional 

Fiber-Reinforced Polymer Matrix Composites," in Annual Book of ASTM Standards. vol. 15.03: American 

Society for Testing and Materials, 2008. 

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

11. "Standard Test Method for Mode I Fatigue Delamination Propagation of Unidirectional Fiber-Reinforced 

Polymer Matrix Composites," draft standard, ASTM International, Committee D30 on Composites, 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.9, DSS Simulia, 2009. 

14. Abaqus Theory Manual, ABAQUS® Standard , Version 6.9, DSS Simulia, 2009. 

15. Abaqus Example Problems Manual, ABAQUS® Standard , Version 6.9, DSS Simulia, 2009. 

16. M. L. Benzeggagh and M. Kenane, "Measurement of mixed-mode delamination fracture toughness of 

14
unidirectional glass/epoxy composites with mixed-mode bending apparatus," Composites Science and 

Technology, vol. 56, pp. 439-449, 1996. 

17. "KaleidaGraph: Data Analysis/Graphing Application for Macintosh and Windows Operating Systems," 

Synergy Software, 1996. 

TABLE I. MATERIAL PROPERTIES [8]. 

Unidirectional Graphite/Epoxy Prepreg 

E11 = 139.4 GPa E22 = 10.16 GPa E33 = 10.16 GPa 

12 = 0.30 13 = 0.30 23 = 0.436 

G12 = 4.6 GPa G13 = 4.6 GPa G23 = 3.54 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 [8] – Figure A1 

GIc = 0.17 kJ/m2 GIIc = GIIIc =0.49 kJ/m2 = 1.62 

Delamination Growth Onset Data [10] – Figures 2, A3,

m0=0.2023 m1=-0.078924

Delamination Growth Rate Data (Paris Law) [8] – Figures 3, A5

GIc = 0.17 kJ/m2 Gth = 0.06 kJ/m2

c=2.44 106 n=10.61

15
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®. Two example input files are 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 shown in Figure A2. The fracture criterion is generated 

experimentally using pure Mode I (GII/GT=0) Double Cantilever Beam (DCB) tests (as shown in 

Figure 1) [7], pure Mode II (GII/GT=1) End-Notched Flexure (ENF) test, and Mixed Mode Bending 

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

16
data (open, blue circles) and mean values (filled, blue circles) are shown in Figure A2 [8]. A 2D 

fracture criterion was suggested by Benzeggah and Kenane [16] using a simple mathematical 

relationship between Gc and GII/GT

Gc = GIc + ( ) GIIc GIc  GII

GT





 

	






. (A1) 

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

and II as shown in Figure A2. The factor  was determined by a curve fit using the LevenbergMarquardt algorithm in the KaleidaGraphTM graphing and data analysis software [17]. 

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) and illustrated in Figure A4 (black curve). 

In the case where a law for growth onset is not available or immediate onset is desirable in the 

analysis, parameters such as suggested in the ABAQUS® example problem [15] may be chosen as 

shown in Figure A3 (red box and curve). The parameters <c3>,<c4> are obtained directly from the 

Paris Law as shown in Figure A5. The parameter <r1> is calculated from the energy release rate 

cutoff value, Gth, and the fracture toughness, GIc, as shown in Figure A5. 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 A5. 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 A6. 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 A6. 

The amplitude curve is then referenced in the definition of the cyclic loading. In the current 

example, prescribed displacements were used to simulate the cyclic opening of the arms of the DCB 

specimen: 

*BOUNDARY,AMPLITUDE=test 

LFRONTP, 1, 1, 0.67 

LFRONTM, 1, 1, -0.67 

where LFRONTP and LFRONTM are node sets located at the tip of each arm where the displacements 

are applied as shown in Figure A1. The factor 0.67 is used to multiply the relative magnitude 

defined by the amplitude curve (shown in Figure A6) and obtain the applied cyclic displacement, 

/2, as shown in Figure A7. 

The direct cyclic approach in ABAQUS® Standard is used to obtain the stabilized cyclic 

response of a structure directly: 

17
*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 A7), 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.1 s for all analyses. Most analyses were run up to 

NT=107 cycles in order to reach the threshold after which delamination growth stops as shown in 

Figures 11 and 12. All other input parameters were varied and the effect on the results 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 ratchetting 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 ratchetting 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]. These control parameters were varied and the effect on the results studied 

as discussed in the main part of this report. 

Example input files 

Two example input files are 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 was set-up to perform the 

desired cyclic analysis. It was found that the prescribed displacements in the static preload step 

had to be small (0.0067 mm) compared to the prescribed displacements (min/2=0.067, 

max/2=0.67) which were used to simulate the cyclic opening of the arms of the DCB specimen. 

18
Analyses where larger preload steps were chosen did not converge in the second step and 

terminated prematurely. 

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>, max/2=0.67)

were also kept constant during all analyses. Other parameters required to define the Fourier 

series which is used to define the cyclic load, the initial time increment and the optional input 

parameters to control the solution in ABAQUS® Standard were varied. The ABAQUS® keywords 

shown in bold type were discussed in detail in the previous paragraphs. 

Input file for fatigue onset and growth analysis 

*HEADING 

DCB-UD-T300/1076, a=30.5 mm 

units: mm, N, MPa 

*** elements, nodes, material, etc 

... 

... 

*NSET, NSET=BONDED, GENERATE 

 253, 1693, 8 

**** surface and contact definition for VCCT **** 

*SURFACE, TYPE=ELEMENT, NAME=VCCT_BOT 

 EL_BOT, S3 

*SURFACE, TYPE=ELEMENT, NAME=VCCT_TOP 

 EL_TOP, S1 

*CONTACT PAIR, INTERACTION=VCCT, ADJUST=BONDED, small sliding 

 VCCT_TOP, VCCT_BOT 

*SURFACE INTERACTION, NAME=VCCT 

 <width> 

*INITIAL CONDITION, TYPE=CONTACT 

 VCCT_TOP, VCCT_BOT, BONDED 

** 

*NSET, NSET=LFRONTP 

 8 

*NSET, NSET=LFRONTM 

 1 

**** VCCT fatigue input 

*parameter 

** Damage and tolerance parameters 

tol=0.001 

** Fracture toughness: 

 GIc = 0.17030 

 GIIc = 0.49360 

 GIIIc = 0.5 

** B-K parameter: 

 eta=1.62 

** width in the plane stress/strain direction 

 width =25.0 

** fatigue crack growth data ** 

 c1=2.8E-09 

 c2=-12.415 

 c3=2.44E+06 

 c4=10.61 

19
 r1=0.353 

 r2=0.9 

*** 

*** amplitude **** 

*amplitude,name=test, DEFINITION=PERIODIC 

1,62.832,0.,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> 

** 

*BOUNDARY, TYPE=DISPLACEMENT 

LFRONTP, 1, 1, 0.00067 

LFRONTM, 1, 1, -0.00067 

** field and history output ** 

*OUTPUT, FIELD, VARIABLE=PRESELECT, FREQ=1 

*Output, history,VARIABLE=PRESELECT,freq=1 

*NODE output,NSET=LFRONTP 

 RF1 

*NODE output,NSET=LFRONTP 

 U1 

*END STEP 

*STEP, INC= 10000 

*direct cyclic,fatigue 

 0.001,0.1,,,11,11,5,10, 

,,10000000,, 

*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,5.E-3,5.E-3,5.E-3

*BOUNDARY,AMPLITUDE=test 

LFRONTP, 1, 1, 0.67 

LFRONTM, 1, 1, -0.67 

** 

** field and history output ** 

*OUTPUT, FIELD, VARIABLE= 

*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=LFRONTP 

 RF1 

*NODE output,NSET=LFRONTP 

 U1 

*END STEP 

20
Input file to verify correct input data and model response 

It is recommended to perform a static analysis first, to verify the correct input parameters. A 

single cycle is performed by replacing the input for direct cyclic analysis 

*direct cyclic,fatigue 

 0.001,0.1,,,11,11,5,10, 

,,10000000,, 

in the second step with 

*STATIC 

0.001, 0.1, 1.E-10, 0.001 

By monitoring the applied displacements, /2, during the simulation, the input data for the cyclic 

analysis can be checked as shown in the example in Figure A8. For all analysis shown, the first step 

consisted of a static preload step up to /2=0.00067. In the second step, the number of terms in the 

Fourier series, <n>, the circular frequency, <omega>, the constant term in the Fourier series, <A0>, 

and the time of a single loading cycle, <tS> were kept the same for all analyses. Further, the factor 

0.67 used to multiply the relative magnitude defined by the amplitude curve to obtain the applied 

displacement, /2, was kept constant. The applied displacement, /2, therefore varied between 

max/2=0.67 mm and min/2=0.067 mm at a frequency of 10 Hz as desired. As an example, the 

starting time, <t0>, the first coefficient of the cosine terms, <A1>, and the first coefficient of the 

sine terms, <B1>, were varied to create different sine waves (in red) and cosine waves (in blue) 

with the same displacement maximum, minimum and frequency as shown in Figure A8. 

*HEADING 

DCB-UD-T300/1076, a=30.5 mm 

units: mm, N, MPa 

*** elements, nodes, material, etc 

... 

... 

*NSET, NSET=BONDED, GENERATE 

 253, 1693, 8 

**** surface and contact definition for VCCT **** 

*SURFACE, TYPE=ELEMENT, NAME=VCCT_BOT 

 EL_BOT, S3 

*SURFACE, TYPE=ELEMENT, NAME=VCCT_TOP 

 EL_TOP, S1 

*CONTACT PAIR, INTERACTION=VCCT, ADJUST=BONDED, small sliding 

 VCCT_TOP, VCCT_BOT 

*SURFACE INTERACTION, NAME=VCCT 

 <width> 

*INITIAL CONDITION, TYPE=CONTACT 

 VCCT_TOP, VCCT_BOT, BONDED 

** 

*NSET, NSET=LFRONTP 

 8 

*NSET, NSET=LFRONTM 

 1 

**** VCCT fatigue input 

*parameter 

** Damage and tolerance parameters 

tol=0.001 

** Fracture toughness: 

21
 GIc = 0.17030 

 GIIc = 0.49360 

 GIIIc = 0.5 

** B-K parameter: 

 eta=1.62 

** width in the plane stress/strain direction 

 width =25.0 

** fatigue crack growth data ** 

 c1=2.8E-09 

 c2=-12.415 

 c3=2.44E+06 

 c4=10.61 

 r1=0.353 

 r2=0.9 

*** 

*** amplitude **** 

*amplitude,name=test, DEFINITION=PERIODIC 

1,62.832,0.,0.55 

0,0.45

*** history data *** 

*** static ramp up step **** 

*STEP, NLGEOM, INC= 10000 

*STATIC 

0.001, 0.1 

** 

*DEBOND,SLAVE=VCCT_TOP,MASTER=VCCT_BOT,FREQ=1 

*FRACTURE CRITERION,TYPE=VCCT,MIXED MODE BEHAVIOR=BK 

1.0e6,1.0e6,1.0e6,<eta> 

** 

*BOUNDARY, TYPE=DISPLACEMENT 

LFRONTP, 1, 1, 0.067 

LFRONTM, 1, 1, -0.067 

** field and history output ** 

*NODE PRINT, NSET=BONDED, GLOBAL=YES, FREQ=1 

 COORD 

*OUTPUT, FIELD, VARIABLE=PRESELECT, FREQ=1 

*Output, history,VARIABLE=PRESELECT,freq=1 

*NODE output,NSET=LFRONTP 

 RF1 

*NODE output,NSET=LFRONTP 

 U1 

*CONTACT OUTPUT, MASTER=VCCT_BOT, SLAVE=VCCT_TOP,NSET=BONDED 

 dbt,enrrt,bdstat 

*END STEP 

*STEP, NLGEOM, INC= 10000 

*** run static analysis with one cycle first to verify input 

*STATIC 

0.001, 0.1, 1.E-10, 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> 

*BOUNDARY,AMPLITUDE=test 

LFRONTP, 1, 1, 0.67 

LFRONTM, 1, 1, -0.67 

** 

** field and history output ** 

*OUTPUT, FIELD, VARIABLE=PRESELECT, FREQ=1 

*ELEMENT OUTPUT 

22
 status 

*CONTACT OUTPUT, MASTER=VCCT_BOT, SLAVE=VCCT_TOP 

 dbt,dbsf,dbs,openbc,crsts,enrrt,efenrrtr,bdstat 

*Output, history,VARIABLE=PRESELECT,freq=1 

*NODE output,NSET=LFRONTP 

 RF1 

*NODE output,NSET=LFRONTP 

 U1 

*CONTACT OUTPUT, MASTER=VCCT_BOT, SLAVE=VCCT_TOP,NSET=BONDED 

 dbt,enrrt,bdstat 

*END STEP 

23
P

P a0

h

2L

2h

B

+

y

z

x

dimensions

B 25.0 mm

2h 3.0 mm

2L 150.0 mm

a0 30.5 mm



layup: [0]24

Figure 1. Double Cantilever Beam Specimen (DCB)

fatigue loading

max/2 0.67 mm

min/2 0.067 mm

R 0.1

f 10.0 Hz

0.00

0.05

0.10

0.15

0.20

0.25

100 101 102 103 104 105 106 107

N

D , cycles

Gmax ,

kJ/m2

 (N/mm)

Figure 2. Delamination growth onset for DCB specimen.

experimental data [10]

curve fit: G= m

0

 * N m

m

0

=0.2023

m

1

=-0.078924

1

24
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

0

10

20

30

40

50

60

70

80 a=30 mm

a=31 mm

a=32 mm

a=33 mm

a=34 mm

a=35 mm

a=40 mm

static [3,4]

fatigue

load P, N

applied opening displacement /2, mm

Figure 4. Critical load-displacement behavior for a DCB specimen.

Pmax =54.3

 max /2=0.67 = const.

G= G

Ic

 =

0.17 kJ/m2

= const.

Pcrit =60.7

 crit /2=0.75

benchmark

10-7

10-6

10-5

10-4

10-3

10-2

10-1

100

10-2 10-1 100

experimental data [8]

Paris Law fit

da/dN, 

mm/cycle

G, kJ/m2

Figure 3. Delamination growth rate (Paris Law) [8].

GIc = 0.17

da/dN = c Gn

c=2.44 *106

n=10.61

R=0.1

f= 10 Hz

cutoff

Gth = 0.060

25
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

time t, seconds

applied 

displacement 

/2, mm

Figure 6. Cyclic fatigue loading for DCB specimen.

 max/2=0.67



min

/2=0.067

R=0.1

f= 10 Hz

applied displacement:

/2 = [0.55+0.45*sin(62.832*(t+0.025))]*0.67

0.00

0.05

0.10

0.15

0.20

0.0 10.0 20.0 30.0 40.0

static benchmark [3,4]

fatigue G= f(a)

G, kJ/m2

delamination length a, mm

Figure 5. Energy release rate - delamination length behavior for DCB specimen.

G

Ic

 =0.17

 crit /2=0.67 = const.

Gmax = 0.8 G

Ic

= 0.1362 

Gth = 0.06

G= G

Ic

 =0.17 = const.

26
10-7

10-6

10-5

10-4

10-3

10-2

10-1

100

10-2 10-1 100

experimental data

Paris Law fit

da/dN, 

mm/cycle

G, kJ/m2

Figure 8. Delamination growth rate (Paris Law).

GIc = 0.17

da/dN = c Gn

c=2.44*106

n=10.61

R=0.1

f= 10 Hz

Gmax = 0.8 G

Ic

increasing a

cutoff

Gth = 0.060

0.00

0.05

0.10

0.15

0.20

0.25

100 101 102 103 104 105 106 107

N

D , cycles

Gmax ,

kJ/m2

 (N/mm)

Figure 7. Delamination growth onset for DCB specimen.

Gmax = 0.8 G

Ic

= 0.1362 

N

D

=150

N = c

1

*Gc

2.8461*10 c -091

c

2 -12.415

curve fit: G= m

0

 * N m

m

0

=0.2023

m

1

=-0.078924

1

2

27
0

10

20

30

40

50

100 101 102 103 104 105 106 107 108

calculated values

delamination

 length 

a, mm

number of cycles NG

Figure 10. Stable delamination growth behavior for DCB specimen.

0.0 100

5.0 10-4

1.0 10-3

1.5 10-3

2.0 10-3

0.0 10.0 20.0 30.0 40.0

calculated values

a/N,

mm/cycle

delamination length a, mm

Figure 9. Delamination growth rate behavior for DCB specimen.

28
0

10

20

30

40

50

100 101 102 103 104 105 106

delamination

 length 

a, mm

time t, seconds

Figure 12. Delamination onset and growth behavior for DCB specimen.

t=N/f

f= 10 Hz= 10 1/s

N

G

(delamination growth)

N cutoff D

(onset)

15 s 3.2*105

 s

0

10

20

30

40

50

100 101 102 103 104 105 106 107

delamination

 length 

a, mm

number of cycles NT

Figure 11. Delamination onset and growth behavior for DCB specimen.

N

G

(delamination growth)

N

D

(delamination onset)

cutoff

150 cycles 3.2*106

 cycles

29
X

Y Z

2L

2h

initial delamination front bonded nodes

X

Z

a. Deformed FE-model of DCB specimen with initial delamination before growth

a

initial delamination front

top surface

 bottom surface

bonded nodes

b. Detail of a FE-model of a DCB specimen

Figure 13. Two-dimensional finite element model of a DCB specimen.

c. Detail of a refined FE-model of a DCB specimen

initial delamination front bonded nodes

a

a=0.5 mm

x a=0.25 mm

z

x

z

30
X

Y Z

X

Y Z

X

Y Z

initial delamination front bonded nodes

b. Detail of a FE-model of a DCB specimen

Figure 14. Details of different two-dimensional finite element models of a DCB specimen.

c. Detail of a refined FE-model of a DCB specimen

initial delamination front bonded nodes

a. Detail of a FE-model of a DCB specimen

initial delamination front bonded nodes

a=1.67 mm

a=1.0 mm

a=1.25 mm

a

x

z

x

z

x

z

31
1

2 3

B

2L

Figure 15. Full three-dimensional finite element model of a DCB specimen.

top and bottom 

surface

a

delamination front

bonded 

nodes

 bottom surface

top surface

bonded 

nodes

x y

z

2h

b. Detail of three-dimensional FE-model around delamination front

a. Deformed three-dimensional FE-model with initial delamination front before growth

a=0.5 mm

a

32
bonded 

nodes

initial straight delamination front

x y

z

a. Deformed three-dimensional model of DCB specimen with a coarse mesh

Figure 16. Coarse full three-dimensional finite element models of a DCB specimen.

b. Deformed three-dimensional model of a DCB specimen with a coarse mesh

initial straight delamination front

bonded 

nodes

x y

z

a=1.0 mm

a=2.0 mm

33
y x

z

bonded 

nodes

initial straight delamination front

a. Deformed continuum-shell model of DCB specimen with a fine mesh

b. Deformed continuum-shell model of a DCB specimen with a coarse mesh

Figure 17. Continuum-shell finite element models of a DCB specimen.

c. Deformed continuum-shell model of a DCB specimen with a coarse mesh

y x

z

y x

z

initial straight delamination front

initial straight delamination front

bonded 

nodes

bonded 

nodes

a=0.5 mm

a=1.0 mm

a=2.0 mm

34
0.0 0.2 0.4 0.6 0.8 1.0

0

10

20

30

40

50

60

model Fig.15

model Fig.16a

model Fig.16b

model Fig.17b

model Fig.17c

model Fig.17a

load P, N

applied opening displacement /2, mm

Figure 19. Critical load-displacement behavior for a DCB specimen.

Pmax =54.3

 max /2=0.67 = const.

3D solid model (C3D8I)

continuum shell model 

(SC8R)

fatigue benchmark

 from Fig.4

0.0 0.2 0.4 0.6 0.8 1.0

0

10

20

30

40

50

60

fatigue benchmark

from Figure 4

2D plane strain 

model (CPE4)-Fig.13

2D plane stress

model (CPS4)-Fig.13

3D solid model

 (C3D8I) - Fig.15

load P, N

applied opening displacement /2, mm

Figure 18. Critical load-displacement behavior for a DCB specimen.

Pmax =54.3

 max /2=0.67 = const.

35
0.00 0.05 0.10 0.15 0.20 0.25 0.30

0

0.05

0.1

0.15

2D plane strain

model (CPE4)-Fig.13

2D plane stress

model (CPS4)-Fig.13

3D solid

model (C3D8I)-Fig.15

time t, s

Figure 21. Computed energy release rate.

R=0.1

f= 10 Hz

Gmax ,

kJ/m2

 (N/mm)

Gmax = 0.8 G

Ic

 = 0.1362 kJ/m2

G

min

 = 0.00136 kJ/m2

0.00 0.05 0.10 0.15 0.20 0.25 0.30

0

0.2

0.4

0.6

0.8

desired fatigue loading from Figure 6

2D plane strain model (CPE4)-Fig.13

2D plane stress model (CPS4)-Fig.13

3D solid model (C3D8I)-Fig.15

time t, s

applied 

displacement 

/2, mm

Figure 20. Cyclic fatigue loading for DCB specimen.

 max/2=0.67



min

/2=0.067

t

s

R=0.1

f= 10 Hz 

t

s

=1/f=0.1

36
N

D

(delamination onset)

delamination

 length 

a, mm

0

10

20

30

40

50

100 101 102 103 104 105 106 107

benchmark

plane strain (CPE4)

plane stress (CPS4)

3D solid (C3D8I)

number of cycles NT

Figure 23. Computed delamination onset and growth: Initial results.

N

G

(delamination growth) cutoff

0.00

0.05

0.10

0.15

-0.4 -0.2 0 0.2 0.4

model Fig.15 - VCCT post-processing

model Fig.15 - ABAQUS

model Fig.16a - VCCT post-processing

model Fig.16a - ABAQUS

model Fig.16b - VCCT post-processing

model Fig.16b - ABAQUS

model Fig.17b - VCCT post-processing

model Fig.17b - ABAQUS

model Fig.17c - VCCT post-processing

model Fig.17c - ABAQUS

model Fig.17a - VCCT post-processing

model Fig.17a - ABAQUS

G

I

,

J/m2

y/B

Figure 22. Computed strain energy release rate distribution across the width of a DCB specimen.

Gmax = 0.8 G

Ic

 = 0.1362 kJ/m2

3D solid model (C3D8I)

continuum shell model (SC8R)

37
N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

benchmark

initial time increment

i

0

=0.01

i

0

=0.0001

i

0

=0.001

number of cycles NT

Figure 25. Computed delamination onset and growth obtained for different initial time increments.

N

G

(delamination growth) cutoff

N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

benchmark

plane strain (CPE4)

plane stress (CPS4)

3D solid (C3D8I)

number of cycles NT

Figure 24. Computed delamination onset and growth: Detail of initial results.

N

G

(delamination growth) cutoff

38
N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

t

0

=0.0

t

0

=-0.025

t

0

=0.025

t

0

=0.0

t

0

=-0.025

t

0

=0.025

number of cycles NT

Figure 27. Computed delamination onset and growth obtained for 

cosine representation of the cyclic loading.

N

G

(delamination growth)

/2=[0.55+A1 cos62.8(t-t0)+B1 sin62.8(t-t0)]*0.67 A1=0.45, B1=0

cutoff

benchmark

50 Fourier terms 11 Fourier terms

N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

t

0

=0.0

t

0

=-0.025

t

0

=0.025

t

0

=0.0

t

0

=-0.025

t

0

=0.025

number of cycles NT

Figure 26. Computed delamination onset and growth obtained for 

sine representation of the cyclic loading.

N

G

(delamination growth)

/2=[0.55+A1 cos62.8(t-t0)+B1 sin62.8(t-t0)]*0.67 B1=0.45, A1=0

cutoff

benchmark

50 Fourier terms 11 Fourier terms

39
25

30

35

40

100 101 102 103 104 105 106 107

crack tip element

a=1.67 mm

a=0.25 mm

a=0.5 mm

a=1.0 mm

a=1.25 mm

benchmark

delamination

 length 

a, mm

number of cycles NG

Figure 29. Computed stable delamination growth behavior for different plane strain models.

N

G

(delamination growth) cutoff

N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

benchmark

crack tip element

a=0.5 mm

a=1.67 mm

a=1.0 mm

a=1.25 mm

a=0.25 mm

number of cycles NT

Figure 28. Computed delamination onset and growth behavior for different plane strain models.

N

G

(delamination growth) cutoff

40
N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

benchmark

crack tip element

a=0.5 mm

a=1.67 mm

a=1.0 mm

a=1.25 mm

a=0.25 mm

number of cycles NT

Figure 31. Computed delamination onset and growth behavior for different plane stress models.

N

G

(delamination growth) cutoff

1

10

100

1000

10000

100000

1000000

0.3 0.5 1 1.3 1.7 2

plane strain (CPE4)

plane stess (CPS4)

3D solid (C3D8I)

continuum shell (SC8R)

time, 

seconds

crack tip element length a, mm

Figure 30. Required analysis time for models with different crack tip element lengths.

41
N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

benchmark

crack tip element a=0.5 mm

plane strain (CPE4)

plane stress (CPS4)

3D solid (C3D8I)

continuum shell (SC8R)

number of cycles NT

Figure 33. Computed delamination onset and growth behavior for different 2D and 3D models.

N

G

(delamination growth) cutoff

N

D

(delamination onset)

delamination

 length 

a, mm

25

30

35

40

100 101 102 103 104 105 106 107

benchmark

crack tip element

C3D8I, a=0.5 mm

SC8R, a=1.0 mm

SC8R, a=2.0 mm

SC8R, a=0.5 mm

C3D8I, a=1.0 mm

C3D8I, a=2.0 mm

number of cycles NT

Figure 32. Computed delamination onset and growth behavior for different continuum models.

N

G

(delamination growth) cutoff

42
2h

initial delamination front bonded nodes

X

Z

initial delamination front

top surface

(VCCT_TOP)

 bottom surface

(VCCT_BOT)

bonded nodes

(BONDED)

Figure A1. Two-dimensional finite element model of a DCB specimen.

node set

(LFRONTP)

node set

(LFRONTM)

DETAIL

/2=0.67

/2=0.67

2

1

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0 0.2 0.4 0.6 0.8 1

mixed mode ratio GII/GT

Gc ,

kJ/m2

DCB, Mode I MMB, Mode I and II ENF, Mode II

Figure A2. Mixed mode fracture criterion.

 = 1.62

G

Ic

 =0.17

experimental data [6]

mean values

curve fit: Gc = G

Ic

 + (GIIc-GIc)(GII/GT)



 GIIc = 0.49

*FRACTURE CRITERION, TYPE=fatigue, MIXED MODE BEHAVIOR=BK, 

TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

43
10-3

10-1

101

103

105

107

0.00 0.05 0.10 0.15 0.20 0.25

Gmax, kJ/m2

 (N/mm)

Figure A4. Delamination growth onset input data.

N, cycles

N = c

1

*Gc2

c

1 0.5

c

2 -0.1

N = c

1

*Gc2

c

1 2.8461e-09

c

2 -12.415

*FRACTURE CRITERION, TYPE=fatigue, MIXED MODE BEHAVIOR=BK, 

TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

immediate onset

curve fit from 

DCB test data

c

2

c

2

0.00

0.05

0.10

0.15

0.20

0.25

100 101 102 103 104 105 106 107

N

D , cycle

Gmax ,

kJ/m2

 (N/mm)

Figure A3. Delamination growth onset for DCB specimen.

curve fit: G= m

0

 * N 

m

0

=0.2023

m

1

=-0.078924

N = c

1

*Gc2

c

1 2.8461e-09

c

2 -12.415

*FRACTURE CRITERION, TYPE=fatigue, MIXED MODE BEHAVIOR=BK, 

TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

m1 c

2

44
0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

amplitude a

time t, seconds

Figure A6. Amplitude curve (relative load magnitude).

=20=62.832

A

0

=0.55

A1

=0

B

1

=0.45

t

0

=-0.025

R=0.1

f= 10 Hz

*AMPLITUDE, NAME=test, DEFINITION=PERIODIC

1,62.832,-0.025,0.55

0,0.45

a=A

0

+A1 cos(t-t0)+B1 sin(t-t0)

10-7

10-6

10-5

10-4

10-3

10-2

10-1

100

10-2 10-1 100

experimental data

Paris Law fit

da/dN, 

mm/cycle

G, kJ/m2

Figure A5. Delamination growth rate (Paris Law) for T300/914C.

Gth = 0.060

GIc = 0.170

da/dN = c

3

 Gc4

c

3

=2.44 106

c

4

=10.61

Gpl = 0.153

r

2

=Gpl/GIc = 0.9

r

1

=Gth/GIc=0.353

*FRACTURE CRITERION, TYPE=fatigue, MIXED MODE BEHAVIOR=BK, 

TOLERANCE=<tol>

<c1>,<c2>,<c3>,<c4>,<r1>,<r2>,<GIc>,<GIIc>,

<GIIIc>,<eta>

c

4

45
0

0.2

0.4

0.6

0.8

1

0 0.05 0.1 0.15 0.2 0.25 0.3

B

1

=0.45, A1=0, t0=0.0

B1=0.45, A1=0, t0=-0.025

B1=0.45, A1=0, t0=0.025

B1=0, A1=0.45, t0=0

B1=0, A1=0.45, t0=-0.025

B1=0, A1=0.45, t0=0.025

time t, s

applied 

displacement 

/2, mm

Figure A8. Cyclic fatigue loading for DCB specimen.



min

/2=0.067

/2=[A0+A1 cos(t-t0)+B1 sin(t-t0)]*0.67 =20=62.832, A0=0.55

static preload step

 max/2=0.67

R=0.1

f= 10 Hz -> t

s

=1/f=0.1

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

time t, s

applied 

displacement 

/2, mm

Figure A7. Cyclic fatigue loading input.

 max/2=0.67



min

/2=0.067

R=0.1

f= 10 Hz -> t

s

=1/f=0.1

/2 =[A0

+B1 sin(t-t0)]*0.67

*DIRECT CYCLIC, fatigue

0.001,0.1,,,11,11,5,200

,,10000000,,

*BOUNDARY,AMPLITUDE=test

LFRONTP, 1, 1, 0.67

LFRONTM, 1, 1, -0.67 t

s

46
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

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) 2. REPORT TYPE 3. DATES COVERED (From - To) 

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 

REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITOR'S ACRONYM(S) 

11. SPONSORING/MONITORING

REPORT NUMBER 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 

ABSTRACT 

18. NUMBER 

OF 

PAGES 

19b. NAME OF RESPONSIBLE PERSON 

a. REPORT b. ABSTRACT c. THIS PAGE 

19b. TELEPHONE NUMBER (Include area code) 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18 

01-07-2010 Contractor Report

Development of a Benchmark Example for Delamination Fatigue Growth

Prediction

NNL09AA00A

814358.02.01.07

Krueger, Ronald

NASA Langley Research Center National Institute of Aerospace (NIA)

Hampton, VA 23681-2199 100 Exploration Way

Hampton, VA 23666

NIA Report No. 2010-04

National Aeronautics and Space Administration

Washington, DC 20546-0001

NASA

NASA/CR-2010-216723

Unclassified - Unlimited

Subject Category 24

Availability: NASA CASI (443) 757-5802

Langley Technical Monitor: Jonathan Ransom

The development of a benchmark example for cyclic delamination growth prediction is presented and demonstrated for a

commercial code. The example is based on a finite element model of a Double Cantilever Beam (DCB) specimen, which is

independent of the analysis software used and allows the assessment of the delamination growth prediction capabilities in

commercial finite element codes. First, the benchmark result was created for the specimen. Second, starting from an initially

straight front, the delamination was allowed to grow under cyclic loading in a finite element model of a commercial code. 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. In general, good agreement between the results obtained from the growth analysis and the

benchmark results could be achieved by selecting the appropriate input parameters. Overall, the results are encouraging but further

assessment for mixed-mode delamination is required.

Composites; Delamination propagation; Finite element analysis; Fracture mechanics; Virtual crack closure technique

U U U

UU

51

STI Help Desk (email: help@sti.nasa.gov)

(443) 757-5802