October 2006 

NASA/CR-2006-214523 

NIA Report No. 2006-06 

Influence of Finite Element Software on Energy 

Release Rates Computed Using the Virtual 

Crack Closure Technique 

Ronald Krueger 

National Institute of Aerospace (NIA), Hampton, Virginia 

Dirk Goetze 

University of Stuttgart, Stuttgart, Germany 
The NASA STI Program Office . . . 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 Office plays a key part in helping NASA 

maintain this important role. 

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for NASA’s 

scientific and technical information. The NASA STI 

Program Office provides access to the NASA STI 

Database, the largest collection of aeronautical and 

space science STI in the world. The Program Office is 

also NASA’s institutional mechanism for 

disseminating the results of its research and 

development activities. These results are published 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 that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results ... even 

providing videos. 

For more information about the NASA STI Program 

Office, 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 (301) 621-0134 

• Phone the NASA STI Help Desk at 

(301) 621-0390 

• Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7121 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 NAS1-02117

 

October 2006 

NASA/CR-2006-214523 

NIA Report No. 2006-06 

Influence of Finite Element Software on Energy 

Release Rates Computed Using the Virtual 

Crack Closure Technique 

Ronald Krueger 

National Institute of Aerospace (NIA), Hampton, Virginia 

Dirk Goetze 

University of Stuttgart, Stuttgart, Germany 
Available from: 

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS) 

7121 Standard Drive 5285 Port Royal Road 

Hanover, MD 21076-1320 Springfield, VA 22161-2171 

(301) 621-0390 (703) 605-6000 

The use of trademarks or names of manufacturers in the report is for accurate reporting and does not 

constitute an official endorsement, either expressed or implied, of such products or manufacturers by the 

National Aeronautics and Space Administration. 
INFLUENCE OF FINITE ELEMENT SOFTWARE ON ENERGY RELEASE RATES

COMPUTED USING THE VIRTUAL CRACK CLOSURE TECHNIQUE

Ronald Krueger

* and Dirk Goetze+

ABSTRACT

Strain energy release rates were computed along straight delamination fronts of Double 

Cantilever Beam, End-Notched Flexure and Single Leg Bending specimens using the 

Virtual Crack Closure Technique (VCCT). The results were based on finite element

analyses using ABAQUS® and ANSYS® and were calculated from the finite element

results using the same post-processing routine to assure a consistent procedure. Mixedmode strain energy release rates obtained from post-processing finite element results

were in good agreement for all element types used and all specimens modeled.

Compared to previous studies, the models made of solid twenty-node hexahedral

elements and solid eight-node incompatible mode elements yielded excellent results. For 

both codes, models made of standard brick elements and elements with reduced

integration did not correctly capture the distribution of the energy release rate across the 

width of the specimens for the models chosen. The results suggested that element types

with similar formulation yield matching results independent of the finite element

software used. For comparison, mixed-mode strain energy release rates were also

calculated within ABAQUS®/Standard using the VCCT for ABAQUS® add on. For all

specimens modeled, mixed-mode strain energy release rates obtained from ABAQUS®

finite element results using post-processing were almost identical to results calculated

using the VCCT for ABAQUS® add on.

1. INTRODUCTION

One of the most common failure modes for composite structures is delamination [1-4]. The

remote loadings applied to composite components are typically resolved into interlaminar

tension and shear stresses at discontinuities that create mixed-mode I, II and III delaminations.

To characterize the onset and growth of these delaminations, the use of fracture mechanics has

become common practice over the past two decades [4-6]. The total strain energy release rate,

GT, the mode I component due to interlaminar tension, GI, the mode II component due to 

interlaminar sliding shear, GII, and the mode III component, GIII, due to interlaminar scissoring 

shear, as shown in Figure 1, have to be calculated when a methodology based on fracture

mechanics is used. In order to predict delamination onset or growth for two-dimensional

problems, these calculated G components are compared to interlaminar fracture toughness

properties measured over a range from pure mode I loading to pure mode II loading [7-9]. A

quasi static mixed-mode fracture criterion is determined by plotting the interlaminar fracture

toughness, Gc, versus the mixed-mode ratio, GII/GT, determined from data generated using 

Double Cantilever Beam (DCB) (GII/GT=0), End Notched Flexure (ENF) (GII/GT=1), and Mixed

Mode Bending (MMB) tests of varying ratios, as shown in Figure 2 for T300/914C [10]. A curve

fit of these data is performed to determine a mathematical relationship between Gc and GII/GT.

 

* Senior Research Scientist, National Institute of Aerospace (NIA), 100 Exploration Way, Hampton, VA 23666 + University of Stuttgart, Stuttgart, Germany

1
[11]. Failure is expected when, for a given mixed-mode ratio GII /GT, the calculated total energy

release rate, GT, exceeds the interlaminar fracture toughness, Gc. An interaction criterion

incorporating the scissoring shear (mode III), however, has not yet been established. The edgecracked torsion test (ECT) is being considered for standardization [12, 13].

The virtual crack closure technique (VCCT) is widely [14, 15] used for computing energy

release rates based on results from continuum (2D) and solid (3D) finite element analyses to

supply the mode separation required when using the mixed-mode fracture criterion. Although the

original publication on VCCT dates back a quarter century [14], the virtual crack closure

technique has only recently been implemented into ABAQUS® [16]. Currently the virtual crack 

closure technique (VCCT) is used mainly by scientists in universities, research institutions and

government laboratories and is usually implemented in specialized codes or used in postprocessing routines in conjunction with general-purpose finite element software.

Lately, an increased interest in using a fracture mechanics based approach to assess the

damage tolerance of composite structures in the design phase and during certification has also

renewed the interest in the virtual crack closure technique. Reliability and repeatability of results

are required to ensure that methodologies based on fracture mechanics will be incorporated into 

the design and certification process. Computed energy release rates as shown in Figure 3 for the

two-dimensional finite element analysis of a hybrid composite tapered flex-beam have shown

that energy release rates computed from results obtained from different finite element codes are

in good agreement [17].

The first objective of this study was to compare computed strain energy release rates using 

VCCT based on results obtained from three-dimensional analyses using different finite element

codes. For comparison, the finite element codes ANSYS® and ABAQUS®were chosen which

were run on different platforms. For this investigation, simple specimens like the double

cantilever beam (DCB) specimen and the End-Notched Flexure (ENF) specimen with 

unidirectional layup and the Single Leg Bending (SLB) specimen with multi-directional layup 

were chosen for full three-dimensional finite element simulations. The specimens were modeled

with different types of solid brick elements. Mixed-mode strain energy release rates were

calculated using a user written post-processing routine. Computed strain energy release rates

across the width of the specimens were compared with results from previous analyses [18, 19].

The second objective of this study was to test and study the implementation of the virtual

crack closure technique in ABAQUS® [16]. Analyses using models of DCB and ENF specimens

with unidirectional layup and models of a SLB specimen with multi-directional layup were

performed. Computed mixed-mode strain energy release rates across the width of the specimens

were compared with results obtained from a user written post-processing routine and with results

from previous analyses [18, 19].

2. SPECIMEN DESCRIPTION

For the current numerical investigation simple specimens such as the Double Cantilever

Beam (DCB), the End-Notched Flexure (ENF) and the Single Leg Bending (SLB) specimens, as

shown in Figure 4, were chosen. The DCB specimen is used to determine the mode I

interlaminar fracture toughness, GIC (GII/GT=0), the ENF specimen is used for mode II

2
(GII/GT=1) fracture toughness testing. The SLB specimen was introduced for the determination

of fracture toughness as a function of mixed-mode I/II ratio [20, 21]. This test may be performed 

in a standard three point bending fixture such as that used for the ENF test. By varying the

relative thickness of the delaminated regions (t1 and t2) various mode mixities may be achieved.

This type of specimen was chosen to study mode separation. Previously, a number of combined 

experimental and numerical studies of these specimens had been performed where the critical

strain energy release rates were evaluated under mode I, mode II and mixed-mode conditions

[18-20, 22, 23].

In general, mode I, mode II and mixed-mode tests are performed on unidirectionally 

reinforced laminates, which means that delamination growth occurs at a [0/0] interface and crack

propagation is parallel to the fibers. For the current study, DCB and ENF specimens made of

T300/1076 graphite/epoxy and unidirectional layup [0]24 (UD24) were modeled. One SLBspecimen modeled was made of C12K/R6376 graphite/epoxy and unidirectional layup [0]32

(UD32). Although this unidirectional layup is desired for standard test methods to characterize

fracture toughness, this kind of delamination growth will rarely occur in real structures.

Previously, a number of combined experimental and numerical studies on specimens with

multidirectional layup were performed where the critical strain energy release rates of various

interfaces were evaluated under mode I, mode II and mixed-mode conditions [20, 22, 23].

Therefore, another SLB-specimen made of C12K/R6376 graphite/epoxy with a multidirectional

layup was also modeled. The stacking sequence [±30/0/-30/0/30/04/30/0/-30/0/-30/30/↑-

30/30/0/30/0/ -30/04/30/0/30/0/±30] was designated D±30, where the arrow denotes the location

of the delamination. The material properties are given in Table 1 and the layups modeled are

summarized in Table 2.

3. ANALYSIS TOOLS

3.1 Virtual Crack Closure Technique

For delaminations in laminated composite materials where the failure criterion is highly

dependent on the mixed-mode ratio (as shown in Figure 2), the virtual crack closure technique

(VCCT) has been widely used for computing energy release rates [14, 15]. Results based on

continuum (2-D) and solid (3-D) finite element analyses provide the mode separation required

when using the mixed-mode fracture criterion.

The mode I and mode II components of the strain energy release rate, GI, GII are computed 

as shown in Figure 5(a) for a 2-D four-node element as an example of VCCT. The terms X’i , Z’i

are the forces at the crack tip at nodal point i and u’l, w’l (u’l*, w’l*) are the displacements at the

corresponding nodal points l and l* behind the crack tip. For geometric nonlinear analysis where

large deformations may occur, both forces and displacements obtained in the global coordinate

system need to be transformed into a local coordinate system (x', z') which originates at the crack

tip as shown in Figure 5(a). The local crack tip system defines the tangential (x' or mode II) and

normal (z' or mode I) coordinate directions at the crack tip in the deformed configuration. The

total energy release rate GT is calculated from the individual mode components as GT =GI +GII

+GIII, where GIII =0 for the two-dimensional case shown in Figure 5(a). The extension to 3-D is

3
straight forward as shown in Figure 5(b). Equations and additional details about VCCT for

higher order solid elements and shell elements can be found in the literature [14, 15].

3.2 User Written VCCT Post-Processing Routines

Although the original publication on VCCT dates back a quarter century [14], the virtual

crack closure technique has only recently been implemented into ABAQUS® [16]. Currently, the

virtual crack closure technique (VCCT) is used mainly by scientists in universities, research 

institutions, government laboratories and research groups in industry who use their own user

written subroutines. The subroutines may either interface directly with the finite element code

during its execution, provided this option has been made available for this particular software, or

operate as an entirely separate post-processing step.

3.2.1 Post-Processing Routine for ABAQUS®

The application of the virtual crack closure technique based on results from finite element

analysis requires access to the element forces at nodes, the nodal point displacements and the

nodal point coordinates. The flow chart, depicted in Figure 6, shows an independent postprocessing routine where input data for VCCT is extracted directly from an ABAQUS® binary 

result file (.fil):

1. Establish interface with finite element software ABAQUS®.

2. Provide input interface to read variable problem-specific external data such as

crack length, Young’s modulus, etc., and control parameters such as element type,

etc., which is generally not hard coded into the subroutine.

3. Create output file to store echo print of control parameters, input data, and retrieved 

data as well as intermediate and final results.

4. Read element forces, or forces at constraints, nodal point coordinates, and nodal

displacements from binary result file (.fil).

5. Store the data retrieved in step 4 in internal common blocks and arrays for further

access during the calculation.

6. Calculate area virtually closed, ∆A, using nodal point coordinates.

7. Calculate local coordinate system x’, y’ and z’ at crack tip for geometrically

nonlinear analysis and arbitrarily shaped delamination front as discussed in section

3.1 or related references [14, 15].

8.

a. Obtain the forces at the crack tip and in front of the crack tip from the

forces at element nodes by summing the forces at common nodes from

elements belonging to the upper surface.

4
b. Obtain relative displacements between the surfaces from the nodal point

displacements.

c. Transform forces and displacements to local crack tip coordinate system.

9. Calculate mixed-mode strain energy release rates using the equations given for the

individual element type and correct for different element lengths and widths if

necessary.

10. Repeat steps 6 through 9 for all nodes along the delamination front and write

results to output file.

3.2.2 Post-Processing Routine for ANSYS Mechanics

TM

The basic strategy to implement VCCT in a post-processing routine is the same regardless

of data retrieved from a binary data file such as the ABAQUS® .fil file or a data file in ASCII

format such as a PATRAN® neutral file. The original post-processing routine to calculate energy

release rates was modified to allow more flexibility. It is finite element software independent and 

uses three ASCII files as input which need to be provided by the user. The modified approach 

shown in Figure 7 can be used with any finite element software as long as the input data files

required can be created. It was tested in conjunction with the finite element software ANSYS®,

for which a macro (vcct-dat.mac) was written to create the data files required. In a first step, the

model data was written to the file model.dat. The nodal coordinates for elements at the

delamination front are written to the file topbot.dat, the displacements and the forces are written 

to the file result.dat.

In the post-processing routine, the software specific interface to the binary output file in

step 1 of the procedure became obsolete as shown in Figure 7. Major changes were also made to

allow the routine to interface with the ASCII data files model.dat, topbot.dat and result.dat and 

to read the data. Additional changes were also required to allow for the conversion of special

element designations unique to ANSYS® to a generic designation used in the user written postprocessing routine.

3.3 VCCT for ABAQUS®

Currently, VCCT for ABAQUS® is an add-on capability to ABAQUS®/Standard Versions

6.5 and 6.6 that provides a specific implementation of the virtual crack closure technique within

ABAQUS®. This new capability enhances ABAQUS® to solve delamination and debonding

problems in composite materials. The implementation is compatible with all the features in

ABAQUS® such as large-scale nonlinear, models of composite structures including continuum

shells, composite materials, cohesive elements, buckling, and contact. The bond line in twodimensional analyses or the plane of delamination in three-dimensional analyses is modeled using

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

capability. Additional element definitions are not required, and the underlying finite element mesh 

and model does not have to be modified [16].

Beyond simple calculation of the mixed-mode strain energy release rates at the crack tip or

5
along the delamination front, the implementation also offers a crack 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 advance. To avoid sudden loss of

stability when the crack tip is advanced, 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. The propagation

capability was not tested since it exceeded the scope of the current study.

4 FINITE ELEMENT MODELING AND ANALYSIS

4.1 Analysis Procedure

The goal of this investigation was to compare computed strain energy release rates based

on results obtained from different finite element codes. For comparison, the finite element codes

ANSYS® and ABAQUS®were chosen. The input files for both codes were generated by a user

written FORTRAN routine (genmesh.f) as shown in the flow chart of Figure 8. This approach 

was used to avoid user interference with input data and ensure that nodal point coordinates,

element topology and material data as well as load and boundary conditions of the models were

identical for both codes. Results from ABAQUS®were saved in the binary .fil results file, results

from ANSYS®were written to three output files (model.dat, results.dat and

topbot.dat) using a user written ANSYS® macro as mentioned above. For both codes, mixedmode strain energy release rates were calculated from finite element results using the same user

written post-processing routine for VCCT to assure a consistent procedure. For comparison,

mixed-mode strain energy release rates were also calculated within ABAQUS®/Standard using

the VCCT for ABAQUS® add on mentioned above.

4.2 Finite Element Selection

To study the influence of element selection on the computed mixed-mode strain energy 

release rates, several three-dimensional solid element types were used in the past to model the

specimens [18, 19]. In ABAQUS®, the use of standard solid eight-node brick element C3D8,

incompatible mode element C3D8I, reduced integration element C3D8R as well as solid twentynode hexahedral elements C3D20 and reduced integration element C3D20R was studied [19].

Shear locking is common in first-order, fully integrated elements, such as C3D8, that are subject

to bending. The numerical formulation of this element gives rise to shear strains that do not

really exist. Therefore, these elements are too stiff in bending, and many elements over the

thickness are required to obtain acceptable results. Elements where a lower-order, reduced 

integration is used to form the element stiffness such as the C3D8R and C3D20R elements

usually provide more accurate results in bending and reduce running time. Incompatible mode

elements, such as C3D8I are also recommended for bending and contact problems. In these

elements, internal deformation modes are added to the standard displacement modes of the

element in order to eliminate the parasitic shear stresses that occur in bending [24]. The use of

6
continuum shell elements was also investigated. The continuum shell elements SC8R provided in

ABAQUS® allow the user to discretize the entire three-dimensional body. The elements look like

three-dimensional solids from a modeling standpoint and have only displacement degrees of

freedom. The continuum shell elements can be stacked to provide more refined throughthickness response. Due to the kinematic and constitutive formulation of the element and the

first-order layer-wise composite theory used for layered elements, the continuum shell elements

appear attractive for modeling [24].

Three different eight-node elements were chosen from the ANSYS® element library [25,

26]. The solid elements SOLID45, SOLID46, and SOLID185 were used for three-dimensional

modeling of the specimen. Analyses with elements SOLID45 and SOLID46 can be performed 

with extra displacement shapes. The element formulation introduces extra internal degrees of

freedom (inaccessible to ANSYS® users) to modify the displacement/velocity gradient tensor.

All internal degrees of freedom are introduced automatically at the element level and are

condensed. Element SOLID46 additionally allows up to 250 different material layers. Element

SOLID185 does not have extra shapes and therefore can exhibit shear locking in bending types

of applications. In such situations, the use of a uniform reduced integration technique with

hourglass stiffness for controlling hourglass modes is recommended (SOLID185R). which is also 

available for element SOLID45 (SOLID45R). The analysis, however, cannot capture the bending 

behavior with a single layer of elements. For example, in the case of a fixed-end cantilever with

a lateral point load four elements over the thickness are usually recommended. Element 185 also 

has an enhanced strain mode (185E). This mode prevents shear locking in bending-dominated

problems and volumetric locking in nearly incompressible cases.

Two twenty-node elements were chosen from the ANSYS® element library. The solid

elements used in this study are 95 and 186 which are higher order versions of the 45 and 185

types mentioned above and exhibit quadratic displacement behavior. The elements can tolerate

irregular shapes without as much loss of accuracy and are well suited to model curved 

boundaries. Element 186 is available as a structural solid and layered solid [25, 26].

4.3 Finite Element Models

Typical three-dimensional finite element models of double cantilever beam (DCB), EndNotched Flexure (ENF) and Single Leg Bending (SLB) specimens are shown in Figures 9 

through 14. For the DCB and SLB specimen made with unidirectional layup only half of the

specimen width (B/2) was modeled as shown in Figure 9, and symmetry conditions were applied

in the xz-plane (y=0). Along the length, all models were divided into different sections with

different mesh refinement as shown in Figure 9(a). A refined mesh of length d=5 mm with 20

elements as used for the UD24 layup is shown in the detail of Figure 9(b), This section length, d,

had been selected in previous studies [18, 19] and was also used during the current

investigations. Across the width, all models were divided into a center section of width, i, and a

refined edge section, j, to capture local edge effects and steep gradients. These sections appear as

dark areas in the full view of the specimen as shown in Figure 9. The specimens with 

unidirectional layup were modeled with six elements through the specimen thickness (2h) as

shown in the detail of Figure 9(b). The delamination was modeled as a discrete discontinuity in

7
the center of the specimen, with separate, unconnected nodes (with identical coordinates) on the

upper and lower surfaces of the delaminated section. For the DCB with unidirectional layup,

finite element meshes were also developed to model the full width B as shown in Figures 10 and

11. A model with a refined mesh near the edges is shown in Figure 10(a) (j=1.0 mm, 8 elements)

and a model with uniform mesh across the width is shown in Figure 10(b). The difference in the

model shown in Figure 11 is a refined zone along the length near the delamination front (d=10 

mm, 20 elements).

For the analysis of the ENF specimen, a model with refined edges was used as shown in

Figure 12. Interpenetration of the cantilever arms was prevented by introducing multi-point

constraints in the plane of delamination along a line of nodes above the support as discussed in a

previous study [19]. The use of multi-point constraints appears advantageous as less modeling

effort is required, and a computationally expensive contact analysis is avoided.

For the SLB specimen, a length of d=3.0 mm with 12 elements in the vicinity of the

delamination front was chosen as shown in Figure 13. This section length, d, had been selected

in previous studies [18] and remained unchanged during the current investigations. For the

specimen with D±30 layup, two plies on each side of the delamination were modeled

individually using one element for each ply as shown in the detail of Figure 13. The adjacent

four plies were modeled by one element with material properties smeared using the rule of

mixtures [27, 28]. This procedure did not calculate the full A-B-D stiffness matrix contributions

of the plies [27, 28]; however, it appeared suitable to enforce a reasonable model size. The

adjacent element extended over the four 0 ̊ plies. The six outermost plies were modeled by one

element with smeared material properties. Additionally, two models with a finer mesh across the

width were created a shown in Figure 14. A model with a refined mesh near the edges is shown

in Figure 14(a) (j=2.0 mm, 8 elements) and a model with uniform mesh across the width is

shown in Figure 14(b).

5. COMPUTED MIXED-MODE ENERGY RELEASE RATES

5.1 Results Obtained from ABAQUS® and ANSYS® Using Post-Processing

The first objective of this study was to compare computed strain energy release rates using 

VCCT based on results obtained from three-dimensional analyses using different finite element

codes. For comparison, the finite element codes ANSYS® and ABAQUS®were chosen. For this

investigation, double cantilever beam (DCB) specimens and the Single Leg Bending (SLB)

specimens with unidirectional and multi-directional layup were chosen for full three-dimensional

finite element simulations [26]. Mixed-mode strain energy release rates were calculated using a

user written post-processing routine as discussed above. Computed strain energy release rates

across the width of the specimens were compared with results from previous analyses [18, 19].

5.1.1 Mode I Energy Release Rates Distribution Across the Width of a DCB Specimen

The computed mode I strain energy release rate values across the width of a DCB

specimen with unidirectional layup were normalized with respect to the value from beam theory,

8
 

!

GI,beam(a) =

12 " a

2

" P2

B2" h3" E1

where a denotes the delamination length, P the external load, B the specimen width, h the

thickness of the cantilever arms as shown in Figure 4 and E1 the modulus of elasticity [8]. The

normalized values were plotted versus the normalized width, y/B, of the specimen as shown in

Figure 15 for results obtained from ABAQUS®. Strain energy release rates obtained from the

model shown in Figure 9 using C3D8I, C3D20 and C3D20R elements are in excellent agreement

with the values from previous analyses, which were added for comparison [18]. The mode I

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

decreasing towards the edges as explained earlier which will cause the initial straight front to 

grow into a curved front as explained in detail in the literature [29-32]. As expected, the mode II

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

shown. The model made of C3D8 elements yielded results, which did not correctly capture the

decrease of the mode I strain energy release rate towards the edge. The mode I strain energy

release rate computed from the model made of C3D8R elements was noticeably higher across the

entire width compared to the other results.

The results of a brief mesh refinement study for the standard eight-node brick element

(C3D8) and eight node brick elements with reduced integration and hourglass control (C3D8R)

are shown in Figures 16 to 18 for a uniform mesh across the width of the specimen. Two

elements over the thickness of one arm were used to compute the distributions shown in

Figure 16. Three elements over the thickness of one arm yielded the results in Figure 17, and 

four elements over the thickness of one arm yielded the results shown in Figure 18. With

increasing number of elements, the strain energy release rates distributions approach the

reference results. Hence, a minimum number of four elements across the thickness is

recommended when using these elements. Results obtained from models using eight node

incompatible mode elements were included in the plots of Figures 16 to 18. The distributions

suggest that as few as two elements across the thickness are sufficient to obtain a good

agreement with reference results.

Mode I energy release rates obtained from the model shown in Figure 9 using the finite

element software ANSYS® are shown in Figure 19 for all element types used in this study. Strain

energy release rates obtained using SOLID186, SOLID95, SOLID45, SOLID46 and 

SOLID185E type elements were in excellent agreement with the values from previous analyses,

which were added for comparison [18]. The model made of SOLID185 elements yielded results,

which did not correctly capture the decrease of the mode I strain energy release rate towards the

edge. The mode I strain energy release rate computed from the model made of SOLID45R and

SOLID185R elements was noticeably higher across the entire width compared to the other

results.

Strain energy release rates obtained from analyses using ABAQUS® and ANSYS® are

plotted in Figures 20 to 23 for corresponding element types and compared to results obtained

from previous analyses [18]. For both finite element codes, mode I energy release rates obtained

from the models using higher order twenty-node elements are in excellent agreement with the

9
values from previous analyses as shown in Figure 20. The models made of standard eight node

brick elements yielded results that did not correctly capture the decrease of the mode I strain

energy release rate towards the edge as shown in Figure 21. The mode I strain energy release rate

computed from models made of eight node brick elements with reduced integration were

noticeably higher across the entire width compared to the other results obtained from previous

analyses [18] as shown in Figure 22. For both finite element codes, mode I energy release rates

obtained from the models using solid eight node incompatible mode elements are in excellent

agreement with the values from previous analyses as shown in Figure 23. The results obtained 

from the model shown in Figure 9 suggest that element types with similar formulation yield 

matching results independent of the software used. Further, mode I energy release rates obtained

from the models using higher order twenty-node elements and eight node incompatible mode

elements are in excellent agreement with the values from previous analyses.

5.1.2 Mixed-Mode Energy Release Rate Distributions Across the Width of a Unidirectional

SLB Specimen

For this investigation, the symmetry of the SLB specimen was taken into account, and

only one half of the specimen width B/2 was modeled as shown in Figure 9. The finite element

model is basically identical to the one used for the DCB specimens discussed in the previous

section, except for the boundary conditions. The influence of element selection on the computed

mixed-mode strain energy release rates was studied using ABAQUS® C3D8, C3D8I, C3D8R,

C3D20 and C3D20R type elements as well as ANSYS® SOLID45 and 45R, SOLID185, 185R

and 185E and SOLI46, SOLID95 and SOLID186 type elements [26].

In Figures 24 to 26, the computed mode I, II and mode III strain energy release rates from

ABAQUS® analyses are plotted versus the normalized width, y/B, of the specimen. Computed

strain energy release rates obtained from models using C3D8I, C3D20 and C3D20R elements are

in excellent agreement with the values from previous analyses [18]. As shown in Figure 24, the

mode I strain energy release rate is fairly constant in the center part of the specimen and

progressively decreases towards the edges as previously discussed for the DCB specimen. The

model made of C3D8 elements yields results that do not correctly capture the decrease of the

mode I strain energy release rate towards the edge. The values computed from the model made

of C3D8R elements appear excessively high. The mode II strain energy release rate as shown in 

Figure 25 is fairly constant across almost the entire width of the specimen, peaking in the

immediate vicinity of the edges. As shown in Figure 26, the mode III contribution is zero in the

center of the specimen peaking to about only 8% of GII at the edges. The model made of C3D8

elements yields results that do not correctly capture the mode II and III distribution. The mode II

values computed from the model made of C3D8R elements appear excessively high across the

entire width of the specimen. These observations support the results obtained from the study of

the UD24 DCB specimen discussed above. Compared to the model made of C3D20 elements,

the models made of C3D8I elements yield nearly the same results and provide a reduced model

size.

Computed mode I, II and mode III strain energy release rates from ANSYS® analyses are

plotted in Figures 27 to 29 versus the normalized width, y/B, of the specimen. Computed strain

energy release rates obtained from models using SOLID186, SOLID95, elements are in excellent

10
agreement with the values from a previous analysis models [18] as shown in Figure 27. The

model made of SOLID185 elements yields results that do not correctly capture the decrease of

the mode I strain energy release rate towards the edge. The values computed from the model

made of SOLID45R and SOLID185R elements appear excessively high. The mode II strain 

energy release rate, as shown in Figure 28, is fairly constant across almost the entire width of the

specimen, peaking in the immediate vicinity of the edges. As shown in Figure 29, the mode III

contribution is zero in the center of the specimen peaking to about only 8% of GII at the edges.

The model made of SOLID185 elements yields results that do not correctly capture the mode II

and III distribution. The mode II values computed from the model made of SOLID45R and

SOLID185R elements appear excessively high across the entire width of the specimen. These

observations support the results obtained from the study of the UD24 DCB specimen discussed

above. Compared to the models made of SOLID186, SOLID95 type elements, the models made

of SOLID45, SOLID46 and SOLID185E type elements yield nearly the same results that are also

in excellent agreement with reference results [18].

Mixed-mode strain energy release rates obtained from analyses using ABAQUS® and

ANSYS® are plotted in Figures 30 to 41 for corresponding element types and compared to

results obtained from previous analyses [18]. For both finite element codes, mode I, II and III

energy release rates obtained from the models using higher order twenty-node elements are in

excellent agreement with the values from previous analyses as shown in Figures 30 to 32. The

models made of standard eight node brick elements yielded results that did not correctly capture

the distribution of the energy release rates across the width as shown in Figure 33 to 35 for the

individual fracture modes. The strain energy release rates computed from models made of eight

node brick elements with reduced integration were noticeably higher across the entire width 

compared to the other results obtained from previous analyses [18] as shown in Figure 36 to 38

for the individual fracture modes. For both finite element codes, mixed-mode energy release

rates obtained from the models using solid eight node incompatible mode elements are in

excellent agreement with the values from previous analyses as shown for the individual models

in Figures 39 to 41. The results suggest that element types with similar formulation yield

matching results independent of the finite element software used. Further energy release rates

obtained from the models using higher order twenty-node elements and eight node incompatible

mode elements are in excellent agreement with the values from previous analyses. Therefore,

solid eight node incompatible mode type elements (ABAQUS® C3D8I and ANSYS® SOLID46)

were chosen for the following study.

5.1.3 Mixed-Mode Energy Release Rate Distributions Across the Width of a

Multidirectional SLB Specimen

In this study, SLB specimens with a multidirectional layup (D±30) were modeled entirely

with solid elements as shown in Figure 13. The model was made of eight node ABAQUS®

C3D8I or ANSYS® SOLID46 elements. It was shown above that compared to the models made

of twenty-node elements, the models made of incompatible mode elements yielded nearly the

same results and provide a reduced model size. Therefore, these element types were chosen for

this and the following studies. Section widths and mesh sizes used in the current investigation 

11
were taken from a previous study where the effect of mesh size across the width on computed

mixed-mode energy release had been investigated [18, 19].

The computed mode I, II and mode III strain energy release rates, as shown in Figures 42 

to 44, were compared with values from previous analyses [18]. As shown in Figure 42, the zone

with a constant GI distribution in the center becomes smaller compared to the UD32 case and the

decrease towards the edges is more pronounced. The decrease is caused by increased anticlastic

bending due to the lower values of bending rigidities in the individual arms. Results obtained 

from ABAQUS® and ANSYS® are in excellent agreement and also agree well with the reference

results. The mode II and mode III distributions are shown in Figures 43 and 44. The mode II

strain energy release rate is fairly constant across almost the entire width of the specimen and

peaks near the edges accompanied by local mode III contribution. Compared to the UD32 layup

these peaks become more visible for specimens with the D±30 layup caused by increased

anticlastic bending. Results obtained from ABAQUS® follow closely the distributions from the

reference results. Results obtained from ANSYS® are in good agreement in the center of the

specimen and also peak towards the edges, however values for y/B< 0 are below and for y/B>0

above values from ABAQUS® and published the reference values [18].

5.2 Results Obtained from ABAQUS® Using Post-Processing and VCCT for ABAQUS®

The second objective of this study was to test and study the implementation of the virtual

crack closure technique in ABAQUS® [16]. Analyses using models of DCB and ENF specimens

with unidirectional layup and models of a SLB specimen with multi-directional layup were

performed. Computed mixed-mode strain energy release rates across the width of the specimens

were compared with results obtained from a user written post-processing routine and with results

from previous analyses [18, 19].

5.2.1 Mode I Energy Release Rates Distribution Across the Width of a DCB Specimen

Strain energy release rates were computed from finite element results obtained from the

model of a full specimen as shown in Figure 10(a) using C3D8I elements. The incompatible

mode elements (C3D8I) were chosen since models consistently yielded the best results in the

studies described above and were least sensitive to mesh refinement. The computed mode I strain

energy release rate values across the width of a DCB specimen with unidirectional layup were

normalized with respect to the value from beam theory as above and plotted versus the

normalized width, y/B, of the specimen as shown in Figure 45. Plotted are results obtained for

different load increments from VCCT for ABAQUS® and from the user written post-processing

routine as shown in the flow chart of Figure 8. Strain energy release rates obtained from both 

methods are practically identical across the entire width and for all load increments.

Additional analyses were performed using models shown in Figures 10 and 11 which 

allowed a discretization study of the effect of different mesh densities as well as mesh transition

from coarse in the center to fine at edge. The results plotted in Figure 46 indicate that

irrespective of mesh density, energy release rates obtained from VCCT for ABAQUS® are almost

identical to results obtained through post-processing.

12
Analyses were repeated with models where the new continuum shell element SC8R was

used. The computed mode I strain energy release rates are shown in Figure 47 for a uniform

mesh across the width and a mesh where the edge region had been refined as shown in Figure 10.

For both models, energy release rates obtained from VCCT for ABAQUS® are almost identical to 

results obtained through post-processing. The results obtained from the model with a uniform

mesh across the width are in good agreement with previous results, which were included for

comparison [18]. The model made of SC8R elements however yields different results in the

refined region near the edge of the specimen when a non-uniform mesh is used. The cause is

related to the element formulation and the associated problem has been fixed by ABAQUS®.

After solving the problem, ABAQUS® provided a computed mode I strain energy release rate

distribution obtained from the non-uniform mesh, which is in good agreement with previous

results. The problem has been fixed for ABAQUS®/Standard 6.6.2 and higher.

5.2.2 Mode II Energy Release Rates Distribution Across the Width of an ENF Specimen

For this investigation, the entire ENF specimen was modeled across the width as shown in

Figure 12. Incompatible mode elements (C3D8I) were chosen since models consistently yielded

the best results in the studies described above and were least sensitive to mesh refinement. In 

Figure 48, the computed mode II and mode III strain energy release rates are normalized with the

reference value, GII,beam, from classical beam theory (not accounting for transverse shear)

 

!

GII,beam(a) =

9" a

2

" P2

16" B2" h3" E1

.

Here, a denotes the delamination length, P the external load, B the specimen width, h the

thickness of the cantilever arms as shown in Figure 4 and E1 the modulus of elasticity. The

results are plotted versus the normalized width, y/B, of the specimen. The mode II strain energy 

release rate is fairly constant across almost the entire width of the specimen, peaking in the

immediate vicinity of the edges. The mode III contribution is zero in the center of the specimen

peaking to about only 5% of GII,beam at the edges. The computed GI values are nearly zero and

therefore are not shown. The results confirm that energy release rates obtained from VCCT for

ABAQUS® are almost identical to results obtained through post-processing.

5.2.3 Mixed-Mode Energy Release Rate Distributions Across the Width of a

Multidirectional SLB Specimen

For the investigation of the SLB specimen with multidirectional layup (D±30), the finite

element models shown in Figures 13 and 14 were used. The model was made of eight node

ABAQUS® incompatible mode C3D8I elements. The computed mode I, II and III strain energy 

release rate distributions across the width of the specimen are shown in Figures 49 to 51. The

zone with a constant GI distribution in the center is smaller compared to the SLB specimen with 

unidirectional layup discussed earlier, and the decrease towards the edges is more pronounced as

shown in Figure 49. The decrease is caused by increased anticlastic bending due to the lower

values of bending rigidities in the individual arms. The mode II distribution, as shown in

Figures 50, is fairly constant across almost the entire width of the specimen and peaks near the

13
edges. The mode III contribution is zero in the center of the specimen peaking at the edges as

shown in Figure 51. Compared to the unidirectional layup (UD32), these peaks become more

visible for specimens with the multidirectional layup (D±30) caused by increased anticlastic

bending. The results plotted in Figure 49 to 51 indicate that irrespective of mesh size mixedmode energy release rates obtained from VCCT for ABAQUS® are almost identical to results

obtained through post-processing.

6. SUMMARY AND CONCLUDING REMARKS

Simple fracture toughness specimens such as the Double Cantilever Beam (DCB)

specimen and the End-Notched Flexure (ENF) specimen with unidirectional layup and the Single

Leg Bending (SLB) specimen with unidirectional and multi-directional layup were modeled with

full three-dimensional solid finite elements and analyzed using the finite element codes ANSYS®

and ABAQUS®. Model data such as nodal point coordinates, element topology and material data

as well as load and boundary conditions were generated using a user written routine to guarantee

that the models were identical for both codes and to avoid interference with input data. Different

element types were used to model the specimens. Delaminations were modeled as discrete

discontinuities. Strain energy release rates were calculated along a straight front using the virtual

crack closure technique. For finite element results based on ABAQUS® and ANSYS®, mixedmode strain energy release rates were calculated from finite element results using the same postprocessing routine to assure a consistent procedure. For comparison, mixed-mode strain energy

release rates were also calculated within ABAQUS®/Standard using the VCCT for ABAQUS®

add on.

Mixed-mode strain energy release rates obtained from post processing ANSYS® and 

ABAQUS® finite element results were in good agreement with results from previous analyses for

all specimens modeled. Compared to previous studies, the models made of solid twenty-node

hexahedral elements (ABAQUS® types C3D20 and C3D20R with reduced integration, ANSYS®

types SOLID95 and SOLID186) and solid eight-node incompatible mode elements (ABAQUS®

type C3D8I, ANSYS® types SOLID45, SOLID46 and SOLID185E) yielded excellent results.

The models made of standard brick elements (ABAQUS® type C3D8, ANSYS® type SOLID185)

yielded results, which were in agreement but compared to reference results did not correctly

capture the distribution of the energy release rate across the width of the specimen for the model

chosen. The values computed from the model made of elements with reduced integration

(ABAQUS® type C3D8R, ANSYS® types SOLID45R and SOLID185R) were also in good

agreement however appeared consistently higher than the reference results. The results suggest

that element types with similar formulation yield matching results independent of the finite

element software used.

For all specimens modeled, mixed-mode strain energy release rates obtained from

ABAQUS® finite element results using post processing were almost identical to results

calculated within ABAQUS®/Standard using the VCCT for ABAQUS® add on. For identical

models, almost identical results were computed which indicates that the post-processing routine

is redundant and may be entirely replaced by the VCCT for ABAQUS® routine. The results also

14
imply that results computed by the VCCT for ABAQUS® routine can readily be compared to 

existing results computed previously using a post-processing routine.

ACKNOWLEDGEMENTS

The authors gratefully acknowledge Mike Sasdelli and Cheryl Liu of ABAQUS® East

and Satish Gandhi of ABAQUS® Headquarters for seriously investigating the problem found 

with the continuum shell element and providing a solution for upcoming versions of ABAQUS®

Standard.

REFERENCES

[1] A.C. Garg, Delamination - A Damage Mode in Composite Structures, Eng. Fracture

Mech., vol. 29, pp. 557-584, 1988.

[2] V.V. Bolotin, Delaminations in Composite Structures: Its origin, Buckling, Growth and 

Stability, Composites Part B: Engineering, vol. 27B, pp. 129-145, 1996.

[3] T.E. Tay, Characterization and Analysis of Delamination Fracture in Composites - An

Overview of Developments from 1990 to 2001, Applied Mechanics Reviews, vol. 56, pp.

1-32, 2003.

[4] T.K. O'Brien, Characterization of Delamination Onset and Growth in a Composite

Laminate, in Damage in Composite Materials, ASTM STP 775,: American Society for

Testing and Materials, 1982, pp. 140-167.

[5] T.K. O'Brien, Fracture Mechanics of Composite Delamination, in ASM Handbook,

Volume 21, Composites: ASM International, 2001, pp. 241-245.

[6] R.H. Martin, Incorporating Interlaminar Fracture Mechanics into Design, in

International Conference on Designing Cost-Effective Composites: IMechE Conference

Transactions, London, U.K., 1998, pp. 83-92.

[7] ASTM D 6671-01, Standard Test Method for Mixed Mode I-Mode II 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, 2000.

[8] ASTM D 5528-94a, 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, 2000.

[9] T.K. O'Brien, Composite Interlaminar Shear Fracture Toughness, GIIc: Shear

Measurement or Sheer Myth ?, in Composite Materials: Fatigue and Fracture, Seventh 

Volume, ASTM STP 1330, 1998, pp. 3-18.

[10] 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, 1997, pp. 60-81.

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

15
[12] S.M. Lee, An Edge Crack Torsion Method for Mode III Delamination Fracture Testing,

J. of Composite Technology and Research., pp. 193-201, 1993.

[13] J.G. Ratcliffe, Characterization of the Edge Crack Torsion (ECT) Test for Mode III

Fracture Toughness Measurement of Laminated Composites, NASA/TM-2004-213269

September 2004 2004.

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

[15] R. Krueger, Virtual Crack Closure Technique: History, Approach and Applications,

Applied Mechanics Reviews, vol. 57, pp. 109-143, 2004.

[16] VCCT for ABAQUS - User's Manual, ABAQUS 2005.

[17] G.B. Murri, J.R. Schaff, and A.L. Dobyns, Fatigue Life Analysis of Tapered Hybrid

Composite Flexbeams, presented at ICCM-14, San Diego, 2003.

[18] R. Krüger, Three Dimensional Finite Element Analysis of Multidirectional Composite

DCB, SLB and ENF Specimens, Institute for Statics and Dynamics of Aerospace

Structures, University of Stuttgart ISD-Report No. 94/2, 1994.

[19] R. Krueger and T.K. O'Brien, A Shell/3D Modeling Technique for the Analysis of

Delaminated Composite Laminates, NASA/TM-2000-210287, ARL-TR-2207, June

2000.

[20] B.D. Davidson, R. Krüger, and M. König, Three Dimensional Analysis of Center

Delaminated Unidirectional and Multidirectional Single Leg Bending Specimens,

Composites Science and Technology, vol. 54, pp. 385-394, 1995.

[21] A. Pieracci, B.D. Davidson, and V. Sundararaman, Nonlinear Analyses of Homogeneous,

Symmetrically Delaminated Single Leg Bending Specimens, Journal Composite Tech.

Res., vol. 20, pp. 170-178, 1998.

[22] B.D. Davidson, R. Krüger, and M. König, Effect of Stacking Sequence on Energy Release

Rate Distributions in Multidirectional DCB and ENF Specimens, Eng. Fracture Mech.,

vol. 55, pp. 557-569, 1996.

[23] B.D. Davidson, R. Krüger, and M. König, Three Dimensional Analysis and Resulting

Design Recommendations for Unidirectional and Multidirectional End-Notched Flexure

Tests, J. Composite Materials, vol. 29, pp. 2108-2133, 1995.

[24] ABAQUS 6.5, Analysis User's Manual, Volume IV: Elements, ABAQUS 2005.

[25] ANSYS Elements Reference, ANSYS Inc., Canonsburg, PA 15317 2005.

[26] D. Goetze, Bruchanalyse von Faserverbundwerkstoff-Proben mit Finiten Elementen und 

der Virtuellen Rißschließungsmethode (Computational Fracture Analysis of Composite

Specimens Based on Three-Dimensional Finite Element Analysis and the Virtual Crack

Closure Technique), Institut fuer Statik und Dynamik der Luft- und

Raumfahrtkonsturktionen, Universitaet Stuttgart, Germany, Diplomarbeit 2003.

[27] S.W. Tsai, Theory of Composite Design: Think Composites, 1992.

[28] S.W. Tsai and H.T. Hahn, Introduction to Composite Materials: Technomic Publishing

Co., Inc., 1980.

[29] I.S. Raju, K.N. Shivakumar, and J.H. Crews Jr., Three-Dimensional Elastic Analysis of a

Composite Double Cantilever Beam Specimen, AIAA J., vol. 26, pp. 1493-1498, 1988.

[30] R. Krüger, M. König, and T. Schneider, Computation of Local Energy Release Rates

Along Straight and Curved Delamination Fronts of Unidirectionally Laminated DCBand ENF - Specimens, in Proceedings of the 34th AIAA/ASME/ASCE/AHS/ASC SSDM

16
Conference, La Jolla, CA: American Institute of Aeronautics and Astronautics,

Washington, 1993, pp. 1332-1342.

[31] B.D. Davidson and R.A. Schapery, Effect of Finite Width on Deflection and Energy

Release Rate of an Orthotropic Double Cantilever Specimen, J. Composite Mat., vol.

22, pp. 640--656, 1988.

[32] B.D. Davidson, An Analytical Investigation of Delamination Front Curvature in Double

Cantilever Beam Specimens, J. Composite Mat., vol. 24, pp. 1124--1137, 1990.

17
Table 1.

Material Properties.

T300/1076 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

C12K/R6376 Unidirectional Graphite/Epoxy Prepreg

E11 = 146.9 GPa E22 = 10.6 GPa E33 = 10.6 GPa

ν12 = 0.33 ν13 = 0.33 ν23 = 0.33

G12 = 5.45 GPa G13 = 5.45 GPa G23 = 3.99 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 2.

Stacking Sequence.

Layup-ID Stacking Sequence Material

UD24 [0]24 T300/1076

UD32 [0]32 C12K/R6376

D±30 [±30/0/-30/0/30/04/30/0/-30/0/-30/30/↑-30/30/30/0/30/0/-30/04/-30/0/30/0/±30] C12K/R6376

18
Figure 1: Fracture Modes. 

Interlaminar sliding shear

Mode II

Interlaminar scissoring shear

Mode II

Interlaminar tension

Mode I

0

100

200

300

400

500

600

700

800

0 0.2 0.4 0.6 0.8 1

Mixed Mode Ratio GII/GT

GC, 

J/m2

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

curve fit: Gc = G

Ic

 + (GIIc-GIc)(GII/GT)

Figure 2. Mixed-mode fracture criterion for T300/914C

η = 1.62

 G

Ic

 =170.3

 GIIc = 493.6

experimental data

mean values

19
0

100

200

300

400

500

600

700

800

0 0.02 0.04 0.06 0.08 0.1 0.12

G

II

, J/m2

delamination length along interface bl, mm

ANSYS®

ABAQUS®

0° graphite

±45° glass

0° graphite

0° graphite

311 312

bl

b

a

resin pocket resin pocket

(a) Flexbeam specimen with loading and micrograph of observed damage

(b) Finite element mesh of the damaged zone

Figure 3. Investigation of hybrid flex-beam specimen [17].

(b) Computed energy release rate GII along delaminated interface

resin pocket

ply ending

312

ply ending

311

delamination along a resin pocket delamination along bl

delamination along b

constant axial load P

cyclic displacement v

detail of damaged area

20
Figure 4. Specimen configurations.

(c) Single Leg Bending Specimen (SLB)

(b) End Notched Flexure Specimen (ENF)

(a) Double Cantilever Beam Specimen (DCB)

+θ layup UD24 

B 25.0 mm

2h 3.0 mm

2L 150.0 mm

a 111.5 mm

P 12.66 N

x

y

z

2h

 p

p

 a 2L

 B

h

x

y

z

+θ

 2L

 B a

 L

 p

2h

z

x

y

 2L

 B a

 L

 t1

 t2

 p

+θ

 layup UD32 layup D±30

B 25.0 mm 25.4 mm

t1 2.03 mm 2.03 mm

t2 2.03 mm 2.03 mm

2L 177.8 mm 177.8 mm

a 34.3 mm 34.3 mm

P 100.0 N 100.0 N

 layup UD24

B 25.0 mm

2h 3.0 mm

2L 150.0 mm

a 30.0 mm

P 503.0 N

21
(a): VCCT for geometrically nonlinear analysis.

x',u',X'

x,u,X

z,w,Z

global system

z',w',Z'

local crack tip system

undeformed state

(outline)

deformed state

Z'i

X'i

u' *

u' 

w' i

k

a

crack closed

a

a

 *

 w' *

GI = - Z'i ( w'l - w'l * ) / ( 2 a )

GII = - X'i ( u'l - u'l * ) / ( 2 a ) 

X,Y,Z: forces

u,v,w: displacements

i

*

L

 b

a

 ZLi

x',u',X'

z',w',Z'

x,u,X

y,v,Y z,w,Z

global

system

local crack tip 

system

a

 XLi

 YLi

(b). VCCT for eight noded solid elements.

 uL *

 wL *

 uL 

 vL wL 

a

Figure 5: Virtual Crack Closure Technique (VCCT).

GI = - Z'i ( w'l - w'l * ) / ( 2 ab )

GII = - X'i ( u'l - u'l * ) / ( 2 ab ) 

GIII = - Y'i ( v'l - v'l * ) / ( 2 ab ) 

X,Y,Z: forces

u,v,w: displacements

22
 output file

grate.out

contains echo input

and results

Figure 6. Flow chart of routine to calculate strain energy release rates from ABAQUS® results

program grate

start program

LRUNIT(1,1)=8

...

CALL INITPF(FNAME,NRU,...

CALL DBRNU(JUNIT)

ABAQUS result file initialization

user specified

input file

grate.data

contains parameters

to chose options

and input/output

file names

read(nunit,*)word,b1,...

read data from user input file

call gstart(kunit,output)

open file for output

ABAQUS® result file

binary format

specimen.fil

INCLUDE ‘aba_param.inc’

enable access to ABAQUS® files

CALL DBFILE(LOP,ARRAY,JRCD)

read each record of the ABAQUS® result file

...

IF(KEY.EQ.15)THEN

call datin(head,1,nnum,...)

...

identify forces at element nodes, nodal

displacements and nodal point coordinates

record and store data in arrays for later use

call gend(iunit)

end

close output file and terminate program

...

call varea(iunit,ntype,...)

call local(iunit,ntype,...)

call trans(iunit,ntype)

call grate(ntype,iunit,...)

calculate area and local coordinate system,

transform forces and displacements,

compute mixed mode energy release rate

n

l

k

j

s

p

q

r

o

m

j

23
call eltyp(eltype,elnum,...)

convert ANSYS® element type to internally

used convention

call datin(head,1,nnum,...)

read model topology from model. dat

 read elements above the delamination in

group ELCRTOP, elements below the

delamination in group ELCRBOT and

nodes at the delamination front in group

CRFRONT from topbot. dat

read displacement at forces at element

nodes from result.dat

store data in arrays and common blocks for

later use

model.dat

topological

data

ANSYS

FE-analysis

 output file

grate.out

contains echo input

and results

program grate

start program

user specified

input file

grate.data

contains parameters

to choose options

and input/output

file names

read(nunit,*)word,b1,...

read data from user input file

call gstart(kunit,output)

open file for output

call gend(iunit)

end

close output file and terminate program

call datout(...)

call orgtopo(nopel,...)

call varea(iunit,ntype,...)

call local(iunit,ntype,...)

call trans(iunit,ntype)

call grate(ntype,iunit,...)

calculate area and local coordinate system,

transform forces and displacements,

compute mixed mode energy release rate

n

l

k

s

p

q

r

o

m

vcct-dat.mac

ANSYS® macro to write

topology, forces and

displacements

topbot.dat

elements and

nodes around

crack

result.dat

displacements

and forces

PARAMETER (NRU=1)

...

Figure 7. Flow chart of modified routine to calculate strain energy release rates.

24
genmesh.f

read input data

calculate nodal point

coordinates, element topology,

position of loads and

boundaries

write input file

Figure 8. Flow chart of entire analysis procedure.

user specified

input file

genmesh.data

contains loads,

dimensions and

material data

ABAQUS® input file

model.inp

ABAQUS® Standard

read input data

perform finite element analysis

calculate GI, GII, GIII

write output

user specified

input file

extract.data

contains loads

case, increment

and

file names

extract.f

read data from result file(s)

calculate GI, GII, GIII

write energy release rates to

output file

GI, GII, GIII

mixed mode ratio

failure index

grate.txt

ANSYS® input file

model.inp

ANSYS MechanicalTM

read input data

perform finite element analysis

write output

vcct-dat.mac

ANSYS macro to write

topology, forces and

displacements

model.dat result.dat topbot.dat

ABAQUS® output

database

model.odb

ABAQUS® result file

model.fil

ABAQUS® Viewer

read output database

plot results

write GI, GII, GIII

GI, GII, GIII

grate.rpt

25
B/2

length element 

b 70.0 mm 14

c 36.5 mm 9

d 10.0 mm 20

e 3.5 mm 2

f 10.0 mm 2

g 20.0 mm 4

i 11.5 mm 12

detail of local

model around

delamination front

2L

c

b

d

e

f

length

B/2 12.5 mm

2L 150.0 mm

delamination front

2h

d

j

length element 

2h 3.0 mm 6

 j 1.0 mm 8

i

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

(b) Deformed model of DCB specimen and detail of refined region around delamination front

(a) Dimensions

g

26
B

length element 

b 70.0 mm 14

c 36.5 mm 9

d 10.0 mm 8

e 3.5 mm 2

f 10.0 mm 3

g 20.0 mm 3

length element 

B 25.0 mm

i 23.0 mm 21

2h 3.0 mm 6

j 1.0 mm 8

i

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

(b) Deformed model of DCB specimen with uniform mesh

(a) Deformed model of DCB specimen with refined edge

length element 

B=i 25.0 mm 25

2h 3.0 mm 6

length element 

b 70.0 mm 14

c 36.5 mm 9

d 10.0 mm 8

e 3.5 mm 2

f 10.0 mm 3

g 20.0 mm 3

B = i

27
length element 

b 70.0 mm 14

c 36.5 mm 9

d 10.0 mm 20

e 3.5 mm 2

f 10.0 mm 3

g 20.0 mm 3

length element 

B 25.0 mm

i 22.5 mm 10

2h 3.0 mm 6

j 1.25 mm 8

Figure 12. Full three-dimensional finite element model of an ENF specimen

B

i

length element 

B=i 25.0 mm 25

2h 3.0 mm 6

length element 

b 21.0 mm 2

c 25.0 mm 6

d 10.0 mm 20

e 13.0 mm 3

f 48.0 mm 9

g 33.0 mm 3

B = i

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

28
Figure 13. Finite element model of SLB specimen with D±30 layup 

i

delamination front

d

±30/0/-30/0/30

04

30/0/-30/0

-30

30

wc

outline of undeformed configuration

wc: deflection in the center of the specimen

detail of local

model around

delamination front

B

i

t1

t2

length element 

B 25.4 mm

i 21.4 mm 10

length element 

b 1.0 mm 1

c 32.79 mm 5

d 3.0 mm 12

e 53.11 mm 5

f 88.9 mm 10

g 1.2 mm 1

length element 

d 3.0 mm 12

j 2.0 mm 8

t1=t2 2.032 mm 5

29
B

i

Figure 14. Full three-dimensional finite element model of a SLB specimen

(b) Deformed model of SLB specimen with uniform mesh 

(a) Deformed model of SLB specimen with refined mesh at the edge 

B = i

length element 

B 25.4 mm

i 23.4 mm 23

j 1.0 mm 1

t1=t2 2.032 mm 5

length element 

b 1.0 mm 1

c 32.79 mm 5

d 3.0 mm 12

e 53.11 mm 5

f 88.9 mm 10

g 1.2 mm 1

length element 

b 1.0 mm 1

c 32.79 mm 5

d 3.0 mm 12

e 53.11 mm 5

f 88.9 mm 10

g 1.2 mm 1

length element 

B 25.4 mm

i 21.4 mm 21

j 2.0 mm 8

t1=t2 2.032 mm 5

30
0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

Krüger [18]

C3D8

C3D8I

C3D8R

y/B

Figure 16. Influence of element selection on computed strain energy release rate

distribution across the width of a DCB specimen with UD24 layup.

c=10 mm

n=20

t=2

Normalized

Mode I

Energy

Release

Rate

0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

C3D8I

C3D8R

C3D20

C3D20R

Krüger [18]

Normalized

Mode I

Energy

Release

Rate

y/B

Figure 15. Influence of ABAQUS® element selection on computed strain energy release rate

distribution across the width of a DCB specimen with UD24 layup. 

31
0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

Krüger [18]

C3D8

C3D8I

C3D8R

Normalized

Mode I

Energy

Release

Rate

y/B

Figure 18. Influence of element selection on computed strain energy release rate

distribution across the width of a DCB specimen with UD24 layup.

c=10 mm

n=20

t=4

0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

C3D8I

C3D8R

Krüger [18]

y/B

Figure 17. Influence of element selection on computed strain energy release rate

distribution across the width of a DCB specimen with UD24 layup. 

c=10 mm

n=20

t=3

Normalized

Mode I

Energy

Release

Rate

32
0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D20

C3D20R

Krüger [18]

186

95

Normalized

Mode I

Energy

Release

Rate

y/B

Figure 20. Computed strain energy release rate distribution across the width of a 

DCB specimen with UD24 layup modeled with twenty-noded solid elements.

0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

Krüger [18]

45

45R

186

46

185

185R

185E

95

Normalized

Mode I

Energy

Release

Rate

y/B

Figure 19. Influence of ANSYS®element selection on computed strain energy release rate

distribution across the width of a DCB specimen with UD24 layup. 

33
0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8R

Krüger [18]

45R

185R

Normalized

Mode I

Energy

Release

Rate

y/B

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

specimen with UD24 layup modeled with reduced integration eight-noded brick elements.

c=10 mm

n=20

0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

Krüger [18]

185

Normalized

Mode I

Energy

Release

Rate

y/B

Figure 21. Computed strain energy release rate distribution across the width of a 

DCB specimen with UD24 layup modeled with eight-noded brick elements. 

c=10 mm

n=20

34
0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

C3D8I

C3D8R

C3D20

C3D20R

Krüger [18]

Energy

Release

Rate GI,

kJ/m2

y/B

Figure 24. Influence of ABAQUS® element selection on computed mode I strain energy release rate

distribution across the width of a SLB specimen with UD32 layup.

0.0

0.5

1.0

1.5

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8I

Krüger [18]

45

46

185E

Normalized

Mode I

Energy

Release

Rate

y/B

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

specimen with UD24 layup modeled with incompatible mode eight-noded brick elements.

35
0.0000

0.0005

0.0010

0.0015

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

C3D8I

C3D8R

C3D20

C3D20R

Krüger [18] Energy

Release

Rate GIII,

kJ/m2

y/B

Figure 26. Influence of ABAQUS® element selection on computed mode III strain energy release

rate distribution across the width of a SLB specimen with UD32 layup.

0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

C3D8I

C3D8R

C3D20

C3D20R

Krüger [18]

Energy

Release

Rate GII,

kJ/m2

y/B

Figure 25. Influence of ABAQUS® element selection on computed mode II strain energy release 

rate distribution across the width of a SLB specimen with UD32 layup.

36
0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

Krüger [18]

45

45R

185

185R

185E

46

95

186

Energy

Release

Rate GII,

kJ/m2

y/B

Figure 28. Influence of ANSYS® element selection on computed mode II strain energy release rate

distribution across the width of a SLB specimen with UD32 layup.

0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

Krüger [18]

45

45R

185

185R

185E

46

95

186

Energy

Release

Rate GI,

kJ/m2

y/B

Figure 27. Influence of ANSYS® element selection on computed mode I strain energy release rate

distribution across the width of a SLB specimen with UD32 layup.

37
0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D20

C3D20R

Krüger [18]

95

186

Energy

Release

Rate GI,

kJ/m2

y/B

Figure 30. Computed mode I strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with twenty-node solid elements.

0.0000

0.0005

0.0010

0.0015

-0.5 -0.4 -0.3 -0.2 -0.1 0

Krüger [18]

45

45R

185

185R

185E

46

95

186

Energy

Release

Rate GIII,

kJ/m2

y/B

Figure 29. Influence of ANSYS® element selection on computed mode III strain energy release rate

distribution across the width of a SLB specimen with UD32 layup.

38
0.0000

0.0005

0.0010

0.0015

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D20

C3D20R

Krüger [18]

95

186

Energy

Release

Rate GIII,

kJ/m2

y/B

Figure 32. Computed mode III strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with twenty-node solid elements.

0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D20

C3D20R

Krüger [18]

95

186

Energy

Release

Rate GII,

kJ/m2

y/B

Figure 31. Computed mode II strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with twenty-node solid elements.

39
0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

Krüger [18]

185

Energy

Release

Rate GII,

kJ/m2

y/B

Figure 34. Computed mode II strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with solid eight-node brick elements.

0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

Krüger [18]

185

Energy

Release

Rate GI,

kJ/m2

y/B

Figure 33. Computed mode I strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with solid eight-node brick elements.

40
0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8R

Krüger [18]

45R

185R

Energy

Release

Rate GI,

kJ/m2

y/B

Figure 36. Computed mode I strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with reduced integration eight-node brick elements.

0.0000

0.0005

0.0010

0.0015

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8

Krüger [18]

185

Energy

Release

Rate GIII,

kJ/m2

y/B

Figure 35. Computed mode III strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with solid eight-node brick elements.

41
0.0000

0.0005

0.0010

0.0015

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8R

Krüger [18]

45R

185R

Energy

Release

Rate GIII,

kJ/m2

y/B

Figure 38. Computed mode III strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with reduced integration eight-node brick elements.

0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8R

Krüger [18]

45R

185R

Energy

Release

Rate GII,

kJ/m2

y/B

Figure 37. Computed mode II strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with reduced integration eight-node brick elements.

42
0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8I

Krüger [18]

45

185E

46

Energy

Release

Rate GII,

kJ/m2

y/B

Figure 40. Computed mode II strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with incompatible mode eight-node brick elements.

0.000

0.005

0.010

0.015

0.020

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8I

Krüger [18]

45

185E

46

Energy

Release

Rate GI,

kJ/m2

y/B

Figure 39. Computed mode I strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with incompatible mode eight-node brick elements.

43
0.000

0.005

0.010

0.015

0.020

0.025

0.030

-0.4 -0.2 0 0.2 0.4

C3D8I

Krüger [18]

46

Energy

Release

Rate GI,

kJ/m2

y/B

Figure 42. Influence of FE software on computed mode I strain energy

release rate distribution across the width of SLB specimen with D±30 layup.

0.0000

0.0005

0.0010

0.0015

-0.5 -0.4 -0.3 -0.2 -0.1 0

C3D8I

Krüger [18]

45

185E

46

Energy

Release

Rate GIII,

kJ/m2

y/B

Figure 41. Computed mode I strain energy release rate distribution across the width of a 

SLB specimen with UD32 layup modeled with incompatible mode eight-node brick elements.

44
0.000

0.005

0.010

0.015

0.020

-0.4 -0.2 0 0.2 0.4

C3D8I

Krüger [18]

46

Energy

Release

Rate GIII,

kJ/m2

y/B

Figure 44. Influence of FE software on computed mode III strain energy

release rate distribution across the width of SLB specimen with D±30 layup.

0.000

0.005

0.010

0.015

0.020

0.025

0.030

-0.4 -0.2 0 0.2 0.4

C3D8I

Krüger [18]

46

Energy

Release

Rate GII,

kJ/m2

y/B

Figure 43. Influence of FE software on computed mode II strain energy

release rate distribution across the width of SLB specimen with D±30 layup.

45
0.00

0.20

0.40

0.60

0.80

1.00

1.20

-0.4 -0.2 0 0.2 0.4

mesh with refined fine edges

VFA

uniform mesh 

VFA

uniform fine mesh

VFA

Energy

Release

Rate,

N/mm

y/B

Figure 46. Computed strain energy release rate distribution 

across the width of a DCB specimen for different models.

0.00

0.20

0.40

0.60

0.80

1.00

1.20

-0.4 -0.2 0 0.2 0.4

33.2 %

50.3 %

75.9 %

100.0 %

33.2 % VFA

50.3 % VFA

75.9 % VFA

100.0 % VFA

Mode I

Energy

Release

Rate,

N/mm

y/B

Figure 45. Computed strain energy release rate distribution 

across the width of a DCB specimen for different load increments.

46
0.000

0.200

0.400

0.600

0.800

1.000

1.200

1.400

-0.4 -0.2 0 0.2 0.4

GII

GIII

GII VFA

GIII VFA

Energy

Release

Rate,

kJ/m2

y/B

Figure 48. Computed mixed mode strain energy release rate

distributions across the width of an ENF specimen.

0.00

0.50

1.00

1.50

-0.4 -0.2 0 0.2 0.4

SC8R uniform mesh 

SC8R fine mesh at edges 

SC8R uniform mesh VFA

SC8R fine mesh at edges VFA

Krüger [18]

SC8R fine mesh at edges VFA - fixed

y/B

Figure 47. Computed strain energy release rate distribution 

across the width of a DCB specimen for continuum shell elements.

Normalized

Mode I

Energy

Release

Rate

47
0.000

0.005

0.010

0.015

0.020

0.025

0.030

-0.4 -0.2 0 0.2 0.4

uniform mesh

VFA 

refined mesh at the edges

VFA 

refined uniform mesh

VFA 

Mode II

Energy

Release

Rate,

N/mm

y/B

Figure 50. Computed mode II strain energy release rate

distribution across the width of a SLB specimen for different models.

0.000

0.005

0.010

0.015

0.020

0.025

0.030

-0.4 -0.2 0 0.2 0.4

uniform mesh

VFA

refined mesh at the edges

VFA

refined uniform mesh

VFA 

Mode I

Energy

Release

Rate,

N/mm

y/B

Figure 49. Computed mode I strain energy release rate

distribution across the width of a SLB specimen for different models.

48
0.000

0.005

0.010

0.015

0.020

0.025

0.030

-0.4 -0.2 0 0.2 0.4

uniform mesh

VFA 

refined mesh at the edges

VFA 

refined uniform mesh

VFA 

Mode III

Energy

Release

Rate,

N/mm

y/B

Figure 51. Computed moe III strain energy release rate

distribution across the width of a SLB specimen for different models.

49
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Contractor Report

 4. TITLE AND SUBTITLE

Influence of Finite Element Software on Energy Release Rates Computed

Using the Virtual Crack Closure Technique

5a. CONTRACT NUMBER

NAS1-02117

 6. AUTHOR(S)

Krueger, Ronald

Goetze, Dirk

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

NASA Langley Research Center National Institute of Aerospace (NIA)

Hampton, VA 23681-2199 100 Exploration Way

 Hampton, VA 23666

 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. 2006-06

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

Langley Technical Monitor: Jonathan Ransom

An electronic version can be found at http://ntrs.nasa.gov

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 24

Availability: NASA CASI (301) 621-0390

19a. NAME OF RESPONSIBLE PERSON

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

14. ABSTRACT

Strain energy release rates were computed along straight delamination fronts of Double Cantilever Beam, End-Notched Flexure and Single Leg Bending

specimens using the Virtual Crack Closure Technique (VCCT). The results were based on finite element analyses using ABAQUS® and ANSYS® and were

calculated from the finite element results using the same post-processing routine to assure a consistent procedure. Mixed-mode strain energy release rates

obtained from post-processing finite element results were in good agreement for all element types used and all specimens modeled. Compared to previous

studies, the models made of solid twenty-node hexahedral elements and solid eight-node incompatible mode elements yielded excellent results. For both

codes, models made of standard brick elements and elements with reduced integration did not correctly capture the distribution of the energy release rate

across the width of the specimens for the models chosen. The results suggested that element types with similar formulation yield matching results independent

of the finite element software used. For comparison, mixed-mode strain energy release rates were also calculated within ABAQUS®/Standard using the

VCCT for ABAQUS® add on. For all specimens modeled, mixed-mode strain energy release rates obtained from ABAQUS® finite element results using

post-processing were almost identical to results calculated using the VCCT for ABAQUS® add on.

15. SUBJECT TERMS

Finite Element Analysis; Computational Fracture Mechanics; Virtual Crack Closure Technique; Composite Material; Fracture

Toughness Specimens

18. NUMBER

 OF 

 PAGES

54

19b. TELEPHONE NUMBER (Include area code)

(301) 621-0390

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

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

1002

5f. WORK UNIT NUMBER

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/CR-2006-214523

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