April 2009 

NASA/TM-2009-215715 

The Development of Directional Decohesion 

Finite Elements for Multiscale Failure 

Analysis of Metallic Polycrystals 

Erik Saether and Edward H. Glaessgen 

Langley Research Center, 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 

Hampton, Virginia 23681-2199 

April 2009 

NASA/TM-2009-215715 

The Development of Directional Decohesion 

Finite Elements for Multiscale Failure 

Analysis of Metallic Polycrystals 

Erik Saether and Edward H. Glaessgen 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

The use of trademarks or names of manufacturers in this 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. 
Table of Contents 

1. Introduction ...................................................................................................................2

2. Cohesive Zone Models .................................................................................................5 

2.1 CZM Parameters Dependent on Direction of Crack Growth ................................11 

3. Decohesion Finite Elements .......................................................................................14

3.1 Directional Decohesion Element Formulation ......................................................19 

4. Simulation of Fracture in 2-D and 3-D Aluminum Polycrystals ................................24

4.1 Study of CZM Parameter Variation on Fracture Development in a 

 2-D Polycrystal .....................................................................................................24 

4.2 Study of Directional Decohesion Elements on Fracture Development 

 in 3-D Polycrystals ................................................................................................28 

5. Concluding Remarks ...................................................................................................34 

6. References ...................................................................................................................36

i
 Abstract

 Atomistic simulations of intergranular fracture have indicated that grain-scale crack 

growth in polycrystalline metals can be direction dependent. At these material length 

scales, the atomic environment greatly influences the nature of intergranular crack 

propagation, through either brittle or ductile mechanisms, that are a function of adjacent 

grain orientation and direction of crack propagation. Methods have been developed to 

obtain cohesive zone models (CZM) directly from molecular dynamics simulations. 

These CZMs may be incorporated into decohesion finite element formulations to 

simulate fracture at larger length scales. A new directional decohesion element is 

presented that calculates the direction of Mode I opening and incorporates a material 

criterion for dislocation emission based on the local crystallographic environment to 

automatically select the CZM that best represents crack growth. The simulation of 

fracture in 2-D and 3-D aluminum polycrystals is used to illustrate the effect of 

parameterized CZMs and the effectiveness of directional decohesion finite elements. 

1
1. Introduction 

At nanometer length scales, the local atomic environment greatly influences the nature 

of intergranular crack propagation. Crack tip process zones can be characterized by brittle 

or ductile mechanisms that are initiated as a function of adjacent grain orientation and 

direction of crack propagation. Experimental evidence for the crystallographic 

dependence of crack tip processes date back to the work of Bilby and Bullough (1954). 

The atomistic response at the crack front can be obtained through molecular dynamics 

(MD) simulation and cast into an aggregate traction-displacement relationship and 

subsequently used as part of a sequential multiscale analysis (Yamakov et al., 2006; 

Saether, 2008).

Sequential multiscale modeling typically involves some form of averaging of physical 

parameters that can serve as initial conditions or provide material parameters to another 

model which is analyzed separately. A desirable aspect of sequential methods is that 

length and time scales between independent material models do not have to be coupled. 

The averaging or homogenization of information across length scales that is inherent to 

sequential multiscale methods is depicted in Figure 1 where a notional coupling is shown 

across domains representing characteristic features at the sub-atomic through structural 

length scales. 

Cohesive zone models (CZM) are based on traction-displacement relationships and are 

a viable tool for sequential multiscale modeling (Saether, 2008). Within a sequential 

analysis framework, CZMs based on averaged atomistic deformation processes for 

intergranular and transgranular fracture can be used to simulate material behavior at 

larger length scales. 

The notional flow of sequential coupling that utilizes CZMs to carry information of 

microscopic failure mechanisms to predict damage progression at larger length scales is 

depicted in Figure 2. Hitherto, the properties of the CZM have been determined 

empirically or heuristically, and mixed-mode fracture has been limited to empirical 

relationships that have been applied at both microscopic and macroscopic length scales. 

2
In Figure 2, CZMs provide the critical transition between inherently atomistic and 

inherently continuum representations. 

QM

(Atoms)

MD

(Lattice)

Polycrystral

Heterogeneous

Continuum

Homegeneous

Continuum

Homogenization

QM

(Atoms)

MD

(Lattice)

Polycrystral

Heterogeneous

Continuum

Homegeneous

Continuum

Homogenization

Polycrystal

(Grains)

Molecular 

Dynamics

Quantum (Lattice)

Mechanics

(Atoms)

Heterogeneous

Continuum

Homogeneous

Continuum

QM

(Atoms)

MD

(Lattice)

Polycrystral

Heterogeneous

Continuum

Homegeneous

Continuum

Homogenization

QM

(Atoms)

MD

(Lattice)

Polycrystral

Heterogeneous

Continuum

Homegeneous

Continuum

Homogenization

Polycrystal

(Grains)

Molecular 

Dynamics

Quantum (Lattice)

Mechanics

(Atoms)

Heterogeneous

Continuum

Homogeneous

Continuum

Figure 1. Hierarchy of models over length scales from quantum mechanical 

through homogeneous continuum. (Adapted from Oden et al., 2006) 

Local Opening (n) and 

Sliding (t) Modes 

Grain Decohesion and 

Dislocation Motion

c) Embedding Cohesive Zones in 

Finite Element Models

1

1’

2

2’

3

3’

d) Idealization of Metallic Microstructures

Vxx Vxx

Vyy

Vyy

Wxy

Wxy

Gn

Gt

a) Atomistic simulations b) Cohesive Zone Representation of 

Atomistic Results

Wp

Op Ocr

Local Opening (n) and 

Sliding (t) Modes 

Grain Decohesion and 

Dislocation Motion

Gn

Gt

c) Embedding Cohesive Zones in 

Finite Element Models

1

1’

2

2’

3

3’

d) Idealization of Metallic Microstructures

Vxx Vxx

Vyy

Vyy

Wxy

Wxy

Gn

Gt

a) Atomistic simulations b) Cohesive Zone Representation of 

Atomistic Results

Wp

Op Ocr

c) Embedding Cohesive Zones in 

Finite Element Models

1

1’

2

2’

3

3’

1

1’

2

2’

3

3’

1

1’

2

2’

3

3’

d) Idealization of Metallic Microstructures

Vxx Vxx

Vyy

Vyy

Wxy

Wxy

Vxx Vxx

Vyy

Vyy

Wxy

Wxy

Gn

Gt

Wp

Op Ocr

Wp

Op Ocr

b) Cohesive Zone Representation of 

Atomistic Results

a) Atomistic simulations

3

Figure 2. Multiscale analysis with cohesive zone models. 

c) Cohesive Zones embedded in 

Finite Element Models 
Characterization of intergranular fracture for different grain boundary (GB) 

orientations is currently very limited. Yamakov et al. (2006) determined CZM parameters 

for intergranular fracture of a symmetric tilt 699 GB, which, depending on the direction 

of propagation, exhibited either ductile or brittle fracture characteristics in Mode I 

opening. With the lack of predictions of CZM parameters for other GB orientations, 

scaling the results obtained by Yamakov et al. (2006) for fracture along a single GB is 

used in the current study to represent other GB orientations.

 CZMs may be incorporated into decohesion finite element formulations to simulate 

fracture at continuum length scales. The CZMs provide constitutive relations that govern 

the separation of material planes due to relative displacement measures that correspond to 

the three fundamental modes of fracture. A new feature to decohesion finite element 

formulations is presented that allows the element to automatically apply appropriate 

brittle or ductile CZMs for crack growth depending on the direction of fracture 

propagation. This new feature thus modifies the element behavior as the finite element 

analysis proceeds. 

Six-node and 12-node directional decohesion finite elements are formulated that 

incorporate Rice’s criterion (Rice, 1992) for dislocation emission to automatically select 

the appropriate CZM to simulate brittle or ductile crack propagation characteristics. For 

intergranular fracture, element input requires the crystallographic orientation of the 

adjacent grains in the microstructure and additional material parameters including 

surface, GB, stable and unstable stacking-fault, and unstable twinning energy in order to 

predict brittle or ductile Mode I fracture.

Section 2 presents a general discussion of cohesive zone models, and Section 3 

presents an overview of decohesion element formulations. Specific formulations of 

directional decohesion elements are detailed in Section 3.1. Finally, Section 4 presents 

illustrative simulations of intergranular fracture in 2-D and 3-D aluminum polycrystal 

models. A 2-D model is used to demonstrate the effect of CZM parameterization, and two 

4
3-D models are used to assess the unique modeling capability of directional decohesion 

finite elements. 

2. Cohesive Zone Models 

 Cohesive zone models were originally developed to represent complicated nonlinear 

fracture processes in ductile and quasi-brittle materials (Dugdale, 1960; Barenblatt, 

1962). CZMs were later developed to describe general adhesion and frictional slip along 

an interface (Maugis, 1992; Kem et al., 1998). The CZM approach is formulated from a 

constitutive relationship based on applied tractions and relative displacements to 

represent separation in various fracture modes of two initially coincident and bonded 

internal surfaces. The relative displacements associated with the creation of a new 

fracture surface for the three fundamental fracture modes are depicted in Figure 3. 

MODE I MODE II MODE III

Opening 

motion 

Sliding 

motion 

Tearing motion

MODE I MODE II MODE III

Opening 

motion 

Sliding 

motion 

Tearing motion

Figure 3. Fundamental fracture modes in solids 

In a CZM-based approach, the area under the CZM traction-displacement curve 

represents the work of separation to open a crack in a particular fracture mode, i.e. the 

fracture toughness, and is given by Tvergaard and Hutchinson (1992) as 

 (1) 3

'

 W '

f

0

c G d

5
where Gc is the work of separation, W is the applied traction, ' is a general form of the 

displacement, and 'f is the critical displacement at which complete separation has 

occurred and the tractions are zero. 

In general, the traction-displacement relationships W   ' are obtained through 

differentiation of a potential I I  ' , which represents the free energy of decohesion. 

The selection of a potential function is typically based on recovering the assumed 

traction-displacement relationship, and particular forms are generally selected for 

analytical convenience. In practice, various forms have been used as shown in Table 1. 

The existence of a work potential yields the work of separation as independent of the 

loading history and of the specific shape of the traction-displacement law. 

Cohesive properties along an interface have typically been approximated using 

empirical data to define the CZMs (Tvergaard and Hutchinson, 1992; Costanzo and 

Allen, 1995; Camacho and Ortiz, 1996; Klein and Gao, 1998; Zavattieri et al., 2001; 

Zavattieri and Espinosa, 2003; Turon et al., 2004). These models are frequently used in 

conjunction with the finite element method (FEM) to study fracture at macroscopic 

length scales in a wide variety of materials. 

The shape of the CZM law represents the behavior of the material near the crack 

tip under load. Various attempts have been made to determine the shape based on 

fundamental bonding characteristics in metals (Rose et al., 1983; Nguyen and Ortiz, 

2002). The most commonly assumed forms of the traction-displacement law have been 

expressed as exponential, bilinear, and trapezoidal functions. A general review of 

various forms of CZMs is presented by Chandra et al. (2002) and summarized in Table 1. 

The table illustrates the basic form of the CZMs, their key parameters, and important 

features.

6
Table 1. Various cohesive zone models and their parameters. (From Chandra et al., 2002)

7
Despite all of the forms of CZMs that have been proposed, a commonly used 

mathematical form not shown in Table 1, is a bilinear constitutive relation. The bilinear 

form is often chosen because of its mathematical simplicity and its suitability for 

representing brittle and ductile fracture in metals (Yamakov et al., 2006), and purely 

brittle fracture in polymeric and ceramic composite materials (Camacho and Ortiz, 1996). 

Bilinear types of cohesive zone models have been used in several recent finite element 

simulations of brittle fracture during multi-axial dynamic loading of ceramic 

microstructures (Zavattieri et al., 2001, 2003). This form of a CZM is shown in Figure 4.

Gc

'o 'f '

Wo

Damage onset

Complete failure

of cohesive zone

W

Gc

'o 'f '

Wo

Gc

'o 'f '

Wo

Damage onset

Complete failure

of cohesive zone

W

Figure 4. Bilinear cohesive zone model.

 In Figure 4, ' represents relative discontinuous displacement jumps between the 

upper and lower nodes, ' 'top - '

bot. 'R indicates the displacement jump corresponding 

to the maximum traction and indicates the onset of damage. Final failure of the cohesive 

zone is assumed after the relative displacement jump yields zero traction. This maximum 

relative displacement jump is indicated as '

f

.

8
 An interpolation procedure using single-mode CZMs has been developed by Turon et 

al. (2004) to generate an effective CZM for mixed-mode applications. This interpolation 

or coupling scheme is depicted in Figure 5. 

Figure 5. Interpolation of CZM components. 

In Figure 5, the relative displacement jumps for single-mode CZMs are denoted by 

G Gtop - Gbot while 'refers to the interpolated displacement jumps. The interpolation of 

“normal” and “shear” CZMs to obtain a single effective CZM is made to represent mixed 

normal and tangential fracture modes. The Mode II and III displacements, GII and GIII, are 

combined by the root sum of the squares to make a combined “shear” displacement, 

'f

'o

Mode I

Shear Mode

GS '

GI E=0

E=1

T

2

III

2

II G (G )  (G ) S (2) 

and the coupled displacement across the interface is defined by 

  2 2

I S ' G  G (3) 

GS

f

GS

o

f GI

o GI

tS

o

o tI

To

'f

'o

Mode I

Shear Mode

GS '

GI E=0

E=1

T

f GI

GS

f

GS

o

o GI

tS

o

o tI

To

9
where is the MacAuley bracket function, x 1⁄2   x  x . Mode mixity parameters 

are defined by 

I G  G

G E 

S

S

 E  E

E 

1 2 2

B 2

2

 (4) 

where E = 0 for pure Mode I loading, and E = 1 for pure shear loading. The fracture 

onset criterion is defined by 

> @ K ' (G )  (G )  (G ) B o 2

I

o 2

S

o 2

I

o (5) 

and the criterion for complete fracture representing crack propagation is derived from the 

B-K critical energy release rate expression (Benzeggagh and Kenane, 1996) given by 

 

K

 ̧

 ̧

1

·

 ̈

 ̈

©

§  

tot

shear

c Ic shear Ic G

G G G G G (6)

where Gtot = GI + Gshear, Gshear = GII + GIII, and K is an empirical factor used to correlate 

with experimental data. This leads to an expression for the final opening displacement 

jump corresponding to crack propagation given by 

  o

I

o

S I

o

I I S B

'

G G  G G  G G ' 

K R I I I

I

 . (7)

 In the previous expressions,

I

, and 

I

 are the individual nodal displacement 

jumps corresponding to opening in Mode I, II and III, respectively, and are userspecified. The , , and are the individual nodal displacements corresponding to 

the onset of damage and are computed from the user-specified initial stiffness values, ,

, and , and peak tractions, , , and , of the single-mode CZMs. The 

I G ,

I

II G

o

tI

III G

t

o G Io G IIo G III

o kI

o kII

o kIIIo

tII

o

III

10
interpolated parameters, '

R and 'f

, define the mixed-mode CZM used in decohesion 

finite elements for the continuum simulation of fracture. 

2.1 CZM Parameters Dependent on Direction of Crack Growth 

Important damage and failure behavior in polycrystalline metals can be simulated 

including intergranular fracture along grain boundaries in metals, transgranular fracture 

along slip planes, and fracture through interfaces between grains and small second phase 

particles. At atomistic length scales, it has been shown that for a common 699 symmetric 

tilt GB, Mode I fracture can propagate by distinct ductile or brittle modes depending on 

the direction of propagation (Yamakov et al., 2006). A synopsis of the distinct directional 

difference in crack propagation behavior and the atomistic determination of 

corresponding CZM parameters as presented by Yamakov et al. (2006) is shown in 

Figure 6. The system was subjected to a hydrostatic tensile load in order to avoid shear 

stresses away from the crack that could induce undesired dislocation formation. To 

ensure crack nucleation and growth at the nanometer scale, the prestress in the range of 

3.75 to 4.25 GPa was considered. 

Figure 6 shows characteristics of the propagating internal crack at different times 

during a MD simulation. At t = 12 ps, the center crack shown at the top of the figure has 

well-developed differences in the crack tip behavior at the two ends. The left crack tip, 

highlighted at t = 68 ps, shows a large degree of blunting due to several formations that 

are identified as the development of twinning (d), cores of partial or twinning 

dislocations (e) and secondary slip (i). The right crack tip, shown at t = 55 ps, 

demonstrates cleavage associated with brittle fracture. The calculated CZM for ductile 

opening shows a peak traction of ~4.0 GPa and an expected large opening displacement 

for final fracture. For brittle fracture, the CZM is characterized by a peak traction of ~5.0 

GPa and a shorter softening region prior to complete opening in Mode I. The CZM 

curves are computed as moving averages of traction-displacement measurements 

obtained from MD simulation (See Yamakov et al., 2006). In the CZM curves presented 

11
in Figure 6, the individual data points and a least-squares curve fit for a prestress of 4.25 

GPa are shown, together with curve fits for the other magnitudes of applied prestress. 

Ductile crack tip Brittle crack tip

Ductile CZM Brittle CZM 

Figure 6. Direction-dependant crack growth behavior incorporated into 

decohesion finite element formulations (From Yamakov et al., 2006). 

12
 The brittle and ductile response of the crack tips shown in Figure 6 is a function of 

several parameters that include the atomic configuration (e.g. the orientation of slip 

planes on either side of the grain boundary), and the magnitude of applied external loads. 

 Figure 7 shows a configuration of a crack plane and slip plane with the angle between 

the two specified by TThe orientation of the Burgers vector, b, lying in the slip plane, 

showing the direction of slip is given by the angle I For a slip direction aligned with the 

crack opening such that I = 0, the Rice criterion (Rice, 1992) for dislocation emission 

associated with Mode I opening along two propagation directions can be written as 

 

0 for ductile fracture

0 for brittle fracture where

8

1 cos sin

2

2

s GB

us

D 

r T T D t  J  J

J D (8) 

where  cosT is associated with the -x direction and cos  S T cosT is associated 

with the +x direction, and Dis the numerical value of the criterion used to predict brittle 

or ductile fracture depending on its sign. For a 699 symmetric tilt GB in aluminum, 

Yamakov et al. (2006) identifies an unstable stacking-fault energy of Jus = 0.168 J/m2, a 

z

y

x

b

T

I

Slip plane 

Crack plane 

Crack front 

Figure 7. Geometric configuration of a crack with an intersecting slip plane. 

13
surface energy of the GB as Js = 0.952 r 0.010 J/m2

, and an excess GB energy given by 

JGB = 0.60 r 0.05 J/m2. Substituting in the energy quantities and multiplying through by 

the denominator in the second term of Equation (8) (because only the sign of D is 

important), the resulting Rice criterion for this system can be expressed as 

 

D r    r T T2 1.0307 0.01 1 cos sin . (9) 

 3. Decohesion Finite Elements 

Decohesion finite elements (also referred to in the literature as “interface elements”, 

“cohesive elements,” or “CZM elements”) were initially introduced by Hilleborg et al. 

(1976) to study fracture in concrete. Needleman (1987, 1990a, 1990b) later formulated a 

decohesion element to study dynamic fracture in isotropic solids. The formulations of the 

elements presented here are based on the work of Beer (1985). The basic topology of 

decohesion elements is formed by separable surfaces that can open according to 

constitutive relationships that relate tractions and displacements. The element 

formulations have found broad application for modeling brittle and ductile fracture 

whether or not crack paths are known a priori. 

When the location of fracture is not known a priori, decohesion elements can be 

placed between all continuum elements within a finite element mesh. While decohesion 

elements can add artificial compliance to the model (Diehl, 2008), more complex 

simulations using adaptive meshing can limit the number of decohesion elements by 

placing these elements only in regions of high stress (Zhang et al., 2007). Figure 8 depicts 

a micromechanical simulation of intergranular fracture with decohesion elements placed 

along all grain boundaries.

The formulation of decohesion finite elements is based on the dimensionality, order of 

element interpolants, and the particular choice of CZMs that are incorporated into the 

formulation. Decohesion elements are formulated with two sets of initially coincident 

14
nodes defining two superposed surfaces forming a cohesive interface. In Figure 9, :R

represents the middle surface of the element, and : and : represent the top and bottom 

surfaces of the element after separation. The particular configuration depicted in Figure 9 

represents an 8-node linear decohesion element used to simulate the cohesive behavior 

between connected 8-node solid continuum elements. Shape functions appropriate to the 

order of assumed variation over the element domain govern the interpolation of quantities 

over the superposed surfaces. 

:+

:O

:-

x1

x2

x3

 Figure 9. Decohesion element cohesive surfaces. 

G3, e

G2, e2

G1, e1

.

. .

.

.

.

.

1 3

4

6

7

8

Figure 8. Embedding decohesion elements along GBs to study microstructural fracture. 

5

2

15
 Fundamental to the response of decohesion elements is the form of the assumed 

cohesive zone models (CZMs) that govern the prediction of separation – fracture – of 

element surfaces (Rabinovitch, 2008). As discussed in Section 2, CZMs are defined by 

displacement jumps, denoted by 'i, which measure the discontinuity of the i

th

displacement components between initially coincident top (  and bottom points on the 

surfaces ( . Thus, the resulting finite element is initially of zero thickness across the 

assumed fracture plane see Beer, 1985; Camanho et al., 2003; and Cook et al., 2001). 

Discretizing each surface into k finite element nodes and assuming that the displacement 

field varies over the surface defined by standard Lagrangian shape functions Nk yields the 

relationship 

 G i

 G i

 (10) 

°

¿

°

3⁄4

1⁄2

°

 ̄

°

®

­

G

G

G

 

°

¿

°

3⁄4

1⁄2

°

 ̄

°

®

­

G  G

G  G

G  G

 

°

¿

°

3⁄4

1⁄2

°

 ̄

°

®

­

'

'

'

 

 

 

k 3k

k 2k

k 1k

3 3

2 2

1 1

3

2

1

N

N

N

where Gik is the i

th relative displacement at the kth node referenced to the middle surface of 

the element :2 Transforming between local (e1,e2,e3) and global (x1,x2,x3) coordinate 

systems using the transformation tensor Tij yields 

^ `u > @ T ^ `' > @> @ T N ^ `G > @ B ^ `G (11) 

Relating the element tractions, Wi, to the relative element displacements, ui, is performed 

using the decohesion element constitutive operator, Dij, as 

^ `W > @ D ^ `u (12) 

The governing equation can be obtained from the principle of virtual work 

d^` ^` ^` ^` u d f d u 0 o

T o T 3 W :  

:

 (13) 

16
which yields, for a geometrically linear problem, the relationship 

> @ ^` ^ ` 3  o

T o T B d f

:

W : 0 (14) 

where [B] is constant. Substituting Equations (11) and (12) into Equation (14) yields the 

standard governing equation for linear elastostatics 

> @ K ^` ^` G f (15) 

 where 

> @ > @ > @> @ 3 :

:o

T o K (16) B D B d

The constitutive matrix, [D], is formulated based on the assumed traction-displacement 

relationship.

 Different features can be incorporated into the element formulation. For example, 

Turon et al. (2004) utilizes a tangent constitutive relation based on bilinear CZMs that 

incorporates a damage parameter as a state variable for mixed-mode failure and a 

variable term that accounts for interpenetration due to possible closing in Mode I fracture. 

Various additional details of decohesion element formulations can be found in Chen et 

al., 1999; Foulk et al., 2000; Alfano and Crisfield, 2001; Goyal et al., 2002; Camanho et 

al., 2003; Segurado and Llorca, 2004; Gustafson et al., 2008; Hamitouche et al., 2008; 

and Van den Bosch et al., 2007. However, most of the differences in these formulations 

pertain to specific definitions of CZM relationships incorporated into the decohesion 

element. 

An additional issue relevant to the development of CZMs regards the initial stiffness 

associated with the constitutive relationship for each fracture mode. For macroscopic 

analysis, this initial stiffness assumes the role of a penalty constraint. A formal penalty 

parameter would assume the highest value that would not cause ill-conditioning in the 

17
global element stiffness matrix. However, in practice, this constraint is set to a high value 

to minimize the relative displacements between the top and bottom surfaces before the 

onset of softening, but does allow numerically small displacements to be calculated that 

are compared with the onset relative displacement at which softening begins. This onset 

displacement, denoted by 'R in Figure 4, is equal to the peak traction divided by the 

selected penalty stiffness and is, thus, a small quantity. However, a consequence of this is 

the introduction of a spurious compliance into the model that can alter results, 

particularly in models where decohesion elements are placed between many continuum 

elements. However, at nanoscopic length scales, this initial stiffness is typically 

calculated based on the elastic properties of the grain boundary or between lattice planes 

because the thickness of the cohesive surface can no longer be generally regarded as 

infinitesimal. In this case, the initial stiffness of the CZM no longer approximates a 

mathematical penalty parameter but is reduced to represent the actual stiffness of a finite 

thickness cohesive zone. Under arbitrarily small loads, the displacement field is not 

homogeneous, and the cohesive surface exhibits an elastic separation due to straining 

below the peak traction. To avoid adding spurious compliance to the model due to a 

numerically finite initial stiffness, either an initial thickness may be assigned to the 

decohesion elements or the elastic properties of adjacent continuum elements may be 

adjusted.

Convergence difficulties have been encountered and have been related to numerical 

problems caused by abrupt changes in the slope of the CZM relationship – such as at the 

vertices of triangular and trapezoidal CZM relations - during loading and unloading 

cycles in which large variations can occur in computing the tangent stiffness matrix (Gao 

and Bower, 2004). Another effect that can hinder convergence involves the type of 

numerical quadrature used to evaluate the element integrals. It has been found that 

Gaussian quadrature tends to couple kinematic degrees of freedom across the element 

that can be seen in the eigenmodes of the element deformation states. Conversely, 

Newton-Cotes quadrature methods act to uncouple the eigenmodes and have been 

associated with an improvement in convergence behavior (Schellekens and de Borst, 

1993). Additional numerical aspects, such as bounds on element size, have been 

18
presented (Allen and Searcy, 2000; Tomar et al., 2004; and Turon et al., 2007), and 

augmented solution schemes that include numerical relaxation to improve convergence 

characteristics have been advanced (Gao and Bower, 2004).

3.1 Directional Decohesion Element Formulation 

The directional dependence of Mode I fracture presented in Section 2.1 motivated the 

development of a directional decohesion finite element formulation for simulating 

intergranular crack propagation. The unique feature of this element is the adaptive 

application of different CZMs depending on the grain orientation and the internal 

computation of the direction of Mode I propagation. This element incorporates the 

direction-dependant CZMs for Mode I opening and utilizes a single Mode II and Mode 

III CZM for computing an effective mixed-mode CZM as explained in Section 2.0. The 

appropriate brittle or ductile Mode I CZM is selected based on Rice’s criterion that 

accounts for the orientation of the surrounding slip systems in the direction of 

propagation (see Figure 6). Both 1-D and 2-D element configurations are presented 

which incorporate Rice’s criterion (Rice, 1992) to apply the appropriate CZM during an 

incremental-iterative solution procedure. The directional decohesion element 

formulations have been encoded into user-defined element (UEL) subroutines to enable 

simulations to be performed using the commercial finite element code 

ABAQUS®1/Standard (2008).

 Until the final effective opening displacement has been reached and the decohesion 

element stiffness terms are zero, calculations are performed using the relative u, v, and w

nodal displacements to approximate the displacement gradients over the element surface 

and determine the direction of propagation. For cases in which various opening 

displacements are equal or zero, a default direction is assumed. A propagation vector, d

*

,

is used to specify the direction of Mode I fracture as is required in the application of 

Rice’s criterion. A second propagation vector,s

* , gives the direction of relative tangential 

 

1 ABAQUS® is manufactured by Dassault Systèmes Simulia Corp. (DSS), Providence, RI, USA 

19
motion of the cohesive surfaces. At microscopic length scales where the behavior of 

discrete atoms and slip planes are significant, there is a distinction between Mode II and 

Mode III sliding and fracture. Fracture by definition generates new free surface area, 

while sliding involves bond breakage which reform between the upper and lower 

surfaces. The determination of relative amounts of tangential sliding and fracture is an 

issue that is currently unresolved and will not be elaborated here; a further discussion can 

be found in Saether (2008).

 The directional decohesion elements developed in the current work consist of a 1-D 6-

node decohesion element for joining 2-D quadratic triangular continuum elements and a 

2-D 12-node decohesion element for joining 3-D quadratic tetrahedral elements. For both 

elements, the orientation of slip planes with respect to the cohesive surface is required 

input. To apply Rice’s criterion, the direction of Mode I opening for each decohesion 

element must be determined. 

 For the 1-D decohesion element, the direction of opening or sliding is simply given 

by the sign of the slope of the relative normal displacements or relative tangential 

displacements along the element. Figure 10 shows the 6-node decohesion element 

configuration. Although the element is quadratic, the Mode I opening direction is 

determined by computing an approximate linear gradient of the relative opening 

displacements, 'vi, over the element using only the end nodes. The degree of sliding in 

Mode II is computed by a summation of all relative x-displacements, 'ui, over the 

element. 

Figure 10. Relative opening displacements for Mode I and Mode II fracture in 

a 6-node quadratic decohesion element.

d

*

y’, v 

'v1 'v2 'v3

, s

*

X1 X2 X3

x’, u 

1e

*

2e

*

'u1 'u2 'u3

20
The kinematics and integration scheme used in the new 6-node 1-D directional 

decohesion element is similar to the element presented by Alfano and Crisfield (2001). 

The sign of the difference between relative displacements determines the directionality of 

Mode I opening and Mode II sliding across the element. In determining the sense of 

relative motion, only the end nodes of the upper or lower surface need to be considered. 

The propagation vectors are then computed as 

  3 1 1 d Sgn v v e

* *

 '  ' (17) 

  1 3 1 s Sgn u u e * * '  ' (18) 

where Sgn indicates the Sign function. 

 For the 2-D 12-node decohesion element shown in Figure 11, the surface may be 

defined by the local mapped x’ and y’ coordinates corresponding to the node positions, 

and can be defined using a general scalar field representation as I = F(x’,y’), where I

represents the variation of a relative displacement component over this region. Using this 

surface, the components of the propagation vectors can be expressed in terms of the 

element x’ and y’ coordinates. The equation of a surface is given by 

 Ax’ + By’ + Cz’ = D (19) 

where x’ and y’ are the mapped coordinate axes, and z’ is the value of the scalar field at 

each (x’,y’) point over the domain. The normal to this surface is given by 

x' y' z' n A e Be Ce * & * *   (20) 

where are unit basis vectors in the local element coordinate system, and the 

components of interest are calculated as 

De

*

21
y ( - ) y ( - ) y ( - )

1 y

1 y

1 y

A 1 2

'

3 1 3

'

2 3 2

'

1

3

'

3

2

'

2

1

'

1

 I I  I I  I I

I

I

I

 

 (21) 

(x - x ) (x - x ) (x - x )

x 1

x 1

x 1

B '

2

'

3 1

'

1

'

2 3

'

3

'

1 2

3

'

3

2

'

2

1

'

1

 I  I  I

I

I

I

 .

 All global displacement quantities are mapped to a local element coordinate system for 

the determination of opening behavior. Failure propagation vectors are computed in the 

local frame and mapped back to the global coordinate system to compare with grain 

orientation angles to select which CZM law is to be applied. 

The 12-node decohesion element shown in Figure 11 is defined with respect to a local 

(x’,y’,z’) element coordinate system, and the kinematics are similar to an element 

presented by Segurado and Llorca (2004). 

d

*

x’, u 

z’, w 

'v3

'u3

'w3

'v2

'u2

'w2

'v1

'u1

s

*

'w1

y’, v 

1 e

*

3 e

* 2 e*

Figure 11. Estimated direction of GB opening and sliding in

mixed-mode fracture. 

22
For Mode I, the expected opening direction is given by the propagation vector, d

*

, as 

opening I 1 I 2 d A e B e * * *  (22) 

where the AI and BI are given in Equations (21) with Ii = 'wi. Because the current 

development assumes a symmetric tilt GB, the orientation angle T specifies the relative 

angle between a slip plane and the cohesive zone crack plane. The direction of opening in 

Mode I is given by BI AI

1 tan  E .

This element can be enhanced to include a full description of grain boundary 

orientation with respect to surrounding slip planes. For the present analysis, it is assumed 

that the crack surface is offset from the nearest slip plane by a single tilt angle, T as 

shown in Figure 9. Thus, defining a local coordinate system, (x’,y’), aligned with the GB, 

the direction of opening is simply given by the sign of the AI component of the 

propagation vector d

*

where Sgn(AI) < 0 and Sgn(AI) > 0 indicates opening in the –x and 

+x directions, respectively. 

 The component directions of Mode II and Mode III fracture are given by Equations 

(21) with Ii = 'ui and Ii = 'vi, respectively. The propagation vectors for these two 

modes are given by 

v III 1 III 2

u II 1 II 2

s A e B e

s A e B e * * *

* * *

 

 

'

' (23) 

 The direction of the in-plane propagation vectors may then be used to quantitatively 

determine the degree of mode mixity for fracture tangential to the crack plane. 

Directional dependence of Mode II and Mode III fracture at atomistic length scales have 

yet to be developed, however, once defined with suitable selection criteria, these may be 

included into the element formulation in an identical manner as with the automated 

selection of appropriate Mode I CZMs.

23
4.0 Simulation of Fracture in 2-D and 3-D Aluminum Polycrystals 

Investigations of fracture in 2-D and 3-D aluminum polycrystals have been 

performed. One study involves a 2-D polycrystal model in which CZM parameters 

describing brittle and ductile intergranular fracture propagation are varied to simulate the 

effect of different GB fracture behavior. This simulation does not include directionality 

of fracture but is instead focused on the effect of CZM parameterization on the overall 

stress-strain behavior of the model and final fracture path. Another study examines two 3-

D polycrystal models to assess the effectiveness of directional decohesion finite elements 

that automatically select a brittle or ductile CZM using Rice’s criterion based on GB 

orientation with respect to adjacent slip planes and the opening direction of Mode I 

fracture.

4.1 Study of CZM Parameter Variation on Fracture Development

in a 2-D Polycrystal 

Figure 12 shows a tessellated 2-D polycrystalline microstructure having many grain 

boundaries. The model has dimensions of lx = ly = 47.5 nm, and consists of 50 grains. 

Six-node triangular finite elements are used to represent the elastic continuum having the 

elastic properties of aluminum (E=72

GPa, Q=0.34). Parametric studies were 

performed by linearly combining the 

CZM parameters using the MD 

extracted cohesive zone curves 

describing brittle and ductile Mode I 

fracture shown in Figure 8 (Glaessgen 

et al., 2006). The polycrystalline 

microstructure is simulated under a 

uniform tensile strain Uy loaded 

incrementally to 22% strain. CZM 

Figure 12. Idealized microstructural model. properties for shear sliding under 

24
Mode II loading were arbitrarily selected to simulate an elastic-perfectly plastic behavior. 

All parameters are shown in Table 2. 

 Using brittle Mode I CZM properties for all the GBs, Figure 13 shows the sequence 

of decohesion processes under normal applied displacements. The CZM parameters for 

Mode I opening obtained from MD analysis yielded an initial GB stiffness that is similar 

to the stiffness of the adjoining grains, and thus does not act as a numerical penalty 

constraint. Therefore, these interfaces begin to exhibit relative displacement at load 

levels below the peak traction in the elastic region of the bilinear CZM. This is depicted 

in Figure 13a where no magnification has been applied to the relative displacements of 

the cohesive zones along the GBs and the black regions indicate significant relative 

displacements along the grain interfaces prior to full opening of the CZMs. At higher 

applied strains, the lowest energy associated with GB opening dictates the location at 

which a dominant microcrack begins to form (Figure 13b). At a critical applied strain, 

the local coalescence of GB opening forms a dominant microcrack as shown in Figure 

13c. It is also seen that due to unloading after the dominant microcrack has fractured the 

model, the stiffness in other decohesion elements acts to pull the opening cohesive zones 

back to an unloaded state. 

(a) Initial opening of GBs, (b) Local crack formation, (c) Dominant microcrack, 

 Hy = 16% Hy = 17% Hy = 22%

 Figure 13. Dominant microcrack formation in a polycrystal model using cohesive 

 zone models. 

25
Prior to demonstrating the behavior of directional decohesion elements, the 

qualitative influence of CZM parameters will be illustrated through comparative 

simulations in which properties for Mode I fracture are assigned as a parametric series of 

proportions of the Mode I brittle and ductile CZM parameters presented in Table 3. These 

parameters, I, were obtained as

Icombined = D Ibrittle + (1 - D) Iductile (24)

where D ranged from 1.0 to 0.0. The Mode II shear parameters were arbitrarily selected 

to simulate an elastic-perfectly plastic behavior. The values chosen were W

o

mx = 0.8 GPa 

and '

c

mx = 100 nm which were kept constant in all simulations. Material properties for 

the grains were assumed as linear isotropic with a Young’s modulus of 100.0 GPa and 

Poisson’s ratio of 0.34.

 Table 2. CZM parameters defining Mode I behavior at the 2-D crack tip.

CZM Brittle Mode I Ductile Mode I 

W

o

mx (GPa) 5.00 3.95

'

o

 (nm) 1.07 0.85

'

c

mx (nm) 2.17 5.65

G (J/m2) 5.43 11.2

Table 3. CZM parameters used to simulate different proportions of elastic

 and plastic Mode I opening behavior at the 2-D crack tip.

Mode I CZM parameters.

G (J/m2 ' ) c CZM D mx (nm) 

C1 1.00 2.17 5.43

C2 0.75 3.04 7.20

C3 0.50 3.91 8.75

C4 0.25 4.78 10.1

C5 0.00 5.65 11.2

Figure 14 shows the resulting stress-strain response of the polycrystal model using 

variations in the proportion of brittle and ductile Mode I cohesive zone behavior. The 

global stress versus strain results demonstrate the imposed trend of decreasing maximum 

26
stress with increasing ductility in the CZM law. Because the relative tangential sliding 

motion in shear is assumed to be independent of Mode I separation, all the parameterized 

models follow a similar extended softening response that continues until ultimate failure. 

Due to the different parameters used to define Mode I cohesive properties of the GBs, 

three different crack paths where exhibited by the models and are shown in Figure 14. 

0.0E+00

5.0E-10

1.0E-09

1.5E-09

2.0E-09

2.5E-09

0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0

Strain

Stress (N/mm2

)

CZM_C1

CZM_C2

CZM_C3

CZM-C4

CZM_C5

Developed crack paths 

C1

C2, C3, C4 

C5

Figure 14. Parametric study of Mode I stress/strain response in a polycrystal 

with brittle vs ductile GBs. 

(%)

4.2 Study of Directional Decohesion Elements on Fracture

Development in 3-D Polycrystals 

 Two 3-D polycrystal models have been developed to illustrate the influence of CZM 

parameters and the functionality of directional decohesion elements that account for grain 

27
orientation and Mode I fracture opening direction. The first model focuses on the 

preferred propagation of fracture along GBs forming a junction between multiple grains. 

In a 3-D model, grain boundaries typically form triple lines and quadripole junctions. 

However, to best illustrate the directional behavior of fracture propagation, the model is 

limited to 3 grains forming a triple junction. Modeling the cohesive surfaces placed along 

each GB using the 12-node directional decohesion element discussed in Section 3.1 will, 

in general, require five independent angles to fully define the orientation of a flat GB in a 

defect-free material and the availability of accurate CZMs to describe directional fracture 

characteristics. This information will be required as input to each decohesion element 

defining a particular GB. However, because of the limited characterization of GB fracture 

behavior, for the present analysis, each GB is assumed as a symmetric tilt GB that can be 

characterized by a single misorientation angle. In order to apply Rice’s criterion, the 

local direction of opening must be computed and the angle, T between the opening 

direction and the local slip plane orientation must be determined. 

 The grains are assumed to be isotropic and linear elastic (E = 100 GPa, Q 0.3), and 

the CZM parameters used in this simulation are presented in Table 4. The directional 

decohesion element detailed in Section 3.1 was coded into a user-defined element for use 

in ABAQUS. This element routine implemented Rice’s criterion which was 

automatically applied during the simulation to determine whether brittle or ductile Mode 

I CZM parameters were to be applied. The model geometry, loading, and results of 

changing the tilt angle for each GB is shown in Figure 15. 

Table 4. CZM parameters defining GB properties in the triple junction model. 

CZM Brittle Mode I Ductile Mode I Shear

W

o

mx (GPa) 5.00 8.00 0.80 

'

o

 (nm) 1.07 1.71 0.40

'

c

mx (nm) 2.17 8.00 1.00

As shown in Figure 15, simulation of a 3-D grain junction under a constant biaxial 

load state shows that different fracture paths are predicted using the formulated 

directional decohesion element by changing the tilt angle associated with each GB.

28
Figure 15. Simulation of a triple junction showing different opening behavior 

with different assumed GB tilt angles using directional decohesion elements. 

The second model is a 3-D polycrystal model that is used to illustrate the effect of 

changing CZM parameters associated with GBs on the global stress-strain response and 

fracture surface creation. This model was constructed by extruding the 2-D model used 

in Section 4.1 in the thickness (z-) direction and is depicted in Figure 16. The same 

loading and boundary conditions shown in Figure 12 were used, however, to fix the 

model against motion in the thickness direction, roller supports were placed on the z = 0 

surface of the model. The model incorporates directional decohesion elements along each 

GB and assigns an initial random material orientation for each grain. CZM parameters are 

systematically varied to simulate the fracture behavior of other types of symmetric tilt 

GBs.

29
 To simplify the geometry of the model, each grain was assumed to have a principal slip 

system normal to the global x-y coordinate system. Thus, grain orientation of the n

th grain 

can be assigned using a single angle, En, ranging from (-S,+Sand defined with respect to 

the x-axis as shown in Figure 17. 

 Figure 16. Model of a 3-D 50-grain microstructure. 

z

y

x

x

y

Ej

Ei

Dk

v

w

Tb1 Tb2

x

y

Ej

Ei

Dk

v

w

Tb1 Tb2

Grain i Grain j 

Figure 17. Geometric description of grain orientation and grain boundaries. 

30
 For a simple tilt grain geometry, the grain boundary orientation is defined by the angle, 

Dk, which is measured counterclockwise to the global y-axis. This angle is computed as 

 ̧

 ̧

1

·

 ̈

 ̈

©

§



 D 

y y

1 x x k

v w

v w Tan (25)

where v and w denote vertex positions defining the ends of the GB. It follows that the tilt 

angles, Tb1 and Tb2, of the GB with respect to the grain slip system are given by 

b2 j k

b1 i k

2

2

T S  E  D

T E  S  D (26) 

A further simplification is required to allow the parameterized MD results for a 

symmetric tilt grain boundary to be used. Therefore, each tilt angle is replaced by an 

average value, T , obtained as 

  2   2 T Tb1  Tb2 Ei  Ej (27) 

The 3-D 25-grain microstructure model, shown in Figure 16, has x and y dimensions 

of 47.5 nm and a thickness of 4.75 nm. The grain properties were modeled as linear 

isotropic with a modulus of 100 GPa and a Poisson’s ratio of 0.34. Although similar 

geometry and grain properties where used for the 2-D and 3-D polycrystal models, CZM 

parameters were not made to coincide because these parameters were used differently in 

the two models to generate effective properties for the decohesion elements in the two 

models that precludes a direct comparison of deformation response. The Mode I CZM 

parameters used = 5 GPa for brittle fracture and = 4 GPa for ductile fracture in 

the current simulation. Both CZMs used an initial stiffness of 4.69 GPa. For Mode II and 

Mode III deformation fields, the maximum traction was held constant at 0.8 GPa, and the 

maximum opening displacement was held constant at 2.0 nm. The varied CZM 

o

mx W omx W

31
parameters used to represent other GB configurations with differing degrees of elastic 

and plastic response are presented in Table 5. These parameters were derived from CZM 

properties calculated by Yamakov et al. (2006) for a symmetric tilt 699 GB in aluminum. 

Table 5. CZM parameters used to simulate different proportions of elastic and

 plastic Mode I opening behavior at the 3-D crack tip.

Brittle Mode I CZM. Peak traction held constant at W

o

mx = 5.0 GPa. 

'

c

mx (nm) G (J/m2 CZM )

D1 2.17 5.42

D2 2.89 7.22

D3 4.33 10.84

D4 6.50 16.25

Du a.

CZM '

c ) G (J/m2

)

ctile Mode I CZM. Peak traction held constant at W

o

mx = 4.0 GP

mx (nm

D1 3.75 7.50

D2 5.00 10.0

D3 6.25 12.5

D4 10.0 20.0

Figure 18 shows a representative evolution of a dominant microcrack in the 3-D 

polycrystal model. As directional decohesion elements along various GBs exhibit 

softening due to damage accumulation, local load redistribution influences which GBs 

will ultimately fracture and create part of the dominant microcrack. For simplicity in the 

CZM formulation and finite element implementation, unloading is assumed to follow a 

path back to the origin to avoid storing and enforcing as a nodal constraint any evolving 

inelastic displacements during simulation. Thus, as various regions away from the 

dominant microcrack are unloaded, the directional decohesion elements along GBs in 

these regions pull the opening GBs back to a closed state. 

32
The global stress-strain response of the model is shown in Figure 19. For each set of 

CZM parameters listed in Table 5, the initial random assignment of grain orientation 

causes the slight difference in the initial elastic stiffness between ductile and brittle 

CZMs such that each model exhibits a similar global stress-strain response until local 

softening due to damage accumulation commences at ~9.35% strain. The variation in the 

maximum ultimate displacement for each model shows the expected increase in the 

degree of softening with increasing in the nonlinear portion of the global stressstrain curves. The difference in predicted crack paths is due to changes in the internal 

load transfer due to variations in the Mode I CZM properties of the GBs. 

c 'mx

H = 0.0% H = 4.40% H = 9.68% 

H = 10.1% H = 11.1% H = 11.7% 

Figure 18. Fracture sequence in the 3-D polycrystal model. 

33
0.0E+00

1.0E-12

2.0E-12

3.0E-12

4.0E-12

5.0E-12

0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 25.0

Strain 

Stress (N/nm2

)

Model D1

Model D2

Model D3

Model D4

Developed crack paths

5. Concluding Remarks 

 The directional decohesion element developed herein adds an additional feature to 

decohesion element formulations that allows the element to alter its cohesive behavior 

during the course of nonlinear finite element analysis. For simulation of microstructural 

failure within polycrystals, the directional decohesion element incorporates Rice’s 

criterion for dislocation emission to automatically apply an appropriate brittle or ductile 

CZM. The selected CZM governs intergranular fracture based on surrounding grain 

orientation and opening direction in various fracture/sliding modes. Two simplistic 

 Figure 19. Stress-strain curves for the 3-D polycrystal model. 

 D1, D2 D3 D4

34
illustrative examples of how results from atomistic simulations might be incorporated 

within the direction decohesion elements were presented. 

 For arbitrary GB orientations, the directional decohesion element will require 

additional rotation angles to be input to define the GB and perform more involved 

calculations of the relative orientation of the cohesive plane and surrounding slip 

systems. Some additional derivations might be required to modify the opening criterion 

for the case when the crack plane does not lie along an atomic plane and intersects a slip 

plane in a skewed orientation. However, the most important requirement for the realistic 

simulation of failure in metallic polycrystals is the determination of CZMs to quantify 

GB fracture behavior for the full multi-dimensional space of GB configurational 

parameters. 

35
6. References 

ABAQUS/Standard User’s Manual, 2008, Dassault Systèmes Simulia Corp. (DSS), 

Providence, RI, USA. 

Alfano, G., and Crisfield, M.A., 2001, “Finite Element Interface Models for the Delamination 

Analysis of Laminated Composites: Mechanical and Computational Issues,” Int. J. Num. 

Meth. Engrg., 50, pp. 1701-1736. 

Allen, D.H., and Searcy, C.R., 2000, “Numerical Aspects of a Micromechanical Model of a 

Cohesive Zone,” J. Reinf. Plast. and Comp., 19, pp. 240-248. 

Barenblatt, G.I., 1959, “The Formation of Equilibrium Cracks during Brittle Fracture. General 

Ideas and Hypothesis. Axially-Symmetric Cracks,” Prikl. Matem. I mekham, 23, pp. 434–

444.

Barenblatt, G.I., 1962, “Mathematical Theory of Equilibrium Cracks in Brittle Fracture,” In: 

Dryden, H.L., von Karman, T. (Eds.), Advances in Applied Mechanics, vol. VII. Academic 

Press, New York, pp. 55–125. 

Beer, G., 1985, “An Isoparametric Joint/Interface Element for the Finite Element Analysis,” 

Int. J. Numer. Meth. Engrg., 21, pp. 585-600. 

Benzeggagh, M.L. and Kenane, M., 1996, “Measurement of Mixed-Mode Delamination 

Fracture Toughnesss of Unidirectional Glass/Epoxy Composites with Mixed-Mode 

Bending Apparatus,” Comp. Sci. and Tech., 49, pp. 439-439. 

Bilby, B.A., and Bullough, R., 1954. “The Formation of Twins by a Moving Crack,” Phil. 

Mag., 45, pp. 631-646. 

Camacho, G.T., and Ortiz, M., 1996, “Computational Modeling of Impact Damage in Brittle 

Materials,” Int. J. Solids and Struct., 33, pp. 2899–2938. 

Camanho, P.P., Davila, C.G., and De Moura, M.F., 2003, “Numerical Simulation of MixedMode progressive Delamination in Composite Materials,” J. Comp. Mat., 37, pp. 1415-

1438.

Chandra, N., Li, H., Shet, C., and Ghonem, H., 2002, ”Some Issues in the Application of 

Cohesive Zone Models for Metal-Ceramic Interfaces,” Int. J. Solids & Struct., 39, pp. 

2827-2855.

Chen, J., Crisfield, M., Kinloch, A.J., Busso, E.P., Matthews, F.L., and Qui, Y., 1999, 

“Predicting Progressive Delamination of Composite Material Specimens via Interface 

Elements,” Mech. Comp. Mat. Struct., 6, pp. 301-317. 

Cook, R.D, Malkus, D.S, and Plesha, M.E, 2001, Concepts and Applications of Finite Element 

Analysis, Third Edition, John Wiley & Sons. 

36
Diehl, T., 2008, “On Using a Penalty-Based Cohesive-Zone Finite Element Approach, Part I: 

Elastic Solution Benchmarks,” Int. J. Adhesion & Adhesives, 28, pp.237-255. 

Dugdale, D.S., 1960, “Yielding of Steel Sheets Containing Slits,” J. Mech. Phys. Solids, 8, pp. 

100–104. 

Foulk, J.W., Allen, D.H., and Helms, K.L.E., 2000, “Formulation of a Three-Dimensional 

Cohesive Zone Model for Application to a Finite Element Algorithm,” Comput. Methods 

Appl. Mech. Engrg., 183, pp. 51-66. 

Gao, Y.F., and Bower, A.F., 2004, “A Simple Technique for Avoiding Convergence Problems 

in Finite Element Simulations of Crack Nucleation and Growth on Cohesive Interfaces,” 

Mod. Sim. Mater. Sci. Engrg., 12, pp. 453-463. 

Glaessgen, E.H., Saether, E., Phillips, D.R., and Yamakov, V., 2006, “Multiscale Modeling of 

Grain-Boundary Fracture: Cohesive Zone Models Parameterized from Atomistic 

Simulations,” Proceedings of the 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural 

Dynamics, and Materials Conference, AIAA-2006-1674-CP, Newport, RI, May 1-4. 

Goyal, V.K., Johnson, E.R., Davila, C.G., and Jaunky, N., 2002, “An Irreversible Constitutive 

Law for Modeling the Delamination Process Using Interface Elements,” NASA/CR-2002-

211758. 

Gustafson, P.A., and Waas, A.M., 2008, “Efficient and Robust Traction Laws for the Modeling 

of Adhesively Bonded Joints,” Proceedings of the 49th AIAA/ASME/ASCE/AHS/ASC 

Structures, Structural Dynamics, and Materials Conference, AIAA 2008-1847, 

Schaumburg, IL. 

Hamitouche, L., Tarfaoui, M., and Vautrin, A., 2008, “An Interface Debonding Law Subject to 

Viscous Regularization for Avoiding Instability: Application to the Delamination 

Problems,” Eng. Fract. Mech., 75, pp. 3084-3100. 

Hilleborg, A., Modeer, M., and P., I.E., 1976, “Analysis of Crack Formation and Crack Growth 

in Concrete by Means of Fracture Mechanics and Finite Elements,” Cement & Concrete 

Res., 6, pp. 773-792. 

 Kem, K.S., McMeeking, R.M., and Johnson, K.I., 1998, “Adhesion, Slip, Cohesive Zones and 

Energy Fluxes for Elastic Spheres in Contact,” J. Mech. Phys. Solids, 46, pp. 243-266. 

 Klein, P. and Gao, H., 1998, “Crack Nucleation and Growth as Strain Localization in a 

Virtual-Bond Continuum,” Engr. Fract. Mech, 61, pp. 21-48. 

 Maugis, E., 1992, “Adhesion of Spheres – the JKR-DMT Transition using a Dugdale Model,”

J. Colloid Interface Sci., 150, pp. 243-269. 

37
Needleman, A., 1987, “A Continuum Model for Void Nucleation by Inclusion Debonding,” J.

Appld. Mech., 54, pp. 525-531. 

Needleman, A., 1990a, “An Analysis of Tensile Decohesion Along an Interface,” J. Mech. 

Phys. Solids, 38, pp. 289–324. 

Needleman, A., 1990b, “An Analysis of Decohesion Along an Imperfect Interface,” Int. J. 

Fract., 42, pp. 21-40. 

 Nguyen, O. and Ortiz, M., 2002, “Coarse-Graining and Renormalization of Atomistic Binding 

Relations and Universal Macroscopic Cohesive Behavior,” J. Mech. Phys. Solids, 50, pp. 

1727-1741. 

Oden, J.T., Prudhomme, S., Romkes, A., and Bauman, P., 2006, "Multi-scale Modeling of 

Physical Phenomena: Adaptive Control of Models," SIAM Journal on Scientific 

Computing, 28, pp. 2359-2389.

Rabinovitch, O., 2008, “Debonding Analysis of Fiber-Reinforced-Polymer Strengthened 

Beams: Cohesive Zone Modeling versus a Linear Elastic Fracture Mechanics Approach,” 

Eng. Fract. Mech., 75, pp. 2842-2859.

Rice, J.R., 1992, “Dislocation Nucleation from a Crack Tip: An Analysis Based on the Peierls 

Concept,” J. Mech. Phys. Solids, 40, pp. 239-271. 

Rice, J.R., and Wang, J.S., 1989, “Embrittlement of Interfaces by Solute Segregation,” Mat. 

Sci. Eng. A, 107, pp. 23–40. 

Saether, E., 2008, “A Multiscale Methodology for Simulating Fracture in Polycrystalline 

Metals,” PhD Thesis, Virginia Polytechnic Institute and State University. 

Schellekens, J.C.J., and de Borst, R., 1993, “On the Numerical Integration of Interface 

Elements,” Int. J. Num. Meth. Engrg., 36, pp. 43-66. 

Segurado, J., and LLorca, J., 2004, “A New Three-Dimensional Interface Finite Element to 

Simulate Fracture in Composites,” Int. J. Solids Struct., 41, pp. 2977-2993. 

Tomar, V., Zhai, J., and Zhou, M., 2004, “Bounds for Element Size in a Variable Stiffness 

Cohesive Finite Element Model,” Int. J. Num. Meth. Engrg., 61, pp. 1894-1920. 

Turon, A., Camanho, P., Costa, J., and Davila, C., 2004, “An Interface Damage Model for the 

Simulation of Delamination Under Variable-Mode Ratio in Composite Materials,” 

NASA/TM-2004-213277. 

Turon, A., Davila, C.G., Camanho, P.P., and Costa, J., 2007, “An Engineering Solution for 

Mesh Size Effects in the Simulation of Delamination using Cohesive Zone Models,” Eng. 

Fract. Mech., 74, pp. 1665-1682. 

38
39

Tvergaard, V., 1990, “Effect of Fibre Debonding in a Whisker-Reinforced Metal,” Mat. Sci. 

and Eng. A, 125, pp. 203–213. 

Tvergaard, V., and Hutchinson, J.W., 1992, “The Relation between Crack Growth Resistance 

and Fracture Process Parameters in Elastic–Plastic Solids,” J. Mech. Phys. Solids, 40, pp. 

1377–1397. 

Van den Bosch, M.J., Schreurs, P.J.G., and Geers, M.G.D, 2007, “A Cohesive Zone Model 

with a Large Displacement Formulation Accounting for Interfacial Fibrilation,” European

J. Mech. A/Solids, 26, pp.1-19. 

Xu, X.P., and Needleman, A., 1993, “Void Nucleation by Inclusion Debonding in a Crystal 

Matrix,” Mod. Sim. Mat. Sci. and Eng., 1, pp. 111–132. 

Yamakov, V., Saether, E., Phillips, D.R., and Glaessgen, E.H., 2006, “Molecular-Dynamics 

Simulation-Based Cohesive Zone Representation of Intergranular Fracture Processes in 

Aluminum,” J. Mech. and Phys. Solids, 54, pp. 1899-1928. 

 Zavattieri, P.D., Raghuram, P.V., and Espinosa, H.D., 2001, “A Computational Model of 

Ceramic Microstructures Subjected to Multi-Axial Dynamic Loading,” J. Mech. Phys. 

Solids, 49, pp. 27-68. 

 Zavattieri, P.D. and Espinosa, H.D., 2003, “An Examination of the Competition Between Bulk 

Behavior and Interfacial Behavior of Ceramics Subjected to Dynamic Pressure-Shear 

Loading,” J. Mech. Phys. Solids, 51, pp. 607-635. 

Zhang, Z., Paulino, G.H., and Celes, W., 2007, “Extrinsic Cohesive Modelling of Dynamic 

Fracture and Microbranching Instability in Brittle Materials,” Int. J. Num. Meth. Engrg.,

72, pp. 893-923. 
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

19a. 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-04-2009 Technical Memorandum

The Development of Directional Decohesion Finite Elements for

Mulstiscale Failure Analysis of Metallic Polycrystals

698259.02.07.07.03.01

Saether, Erik; Glaessgen, Edward H.

NASA Langley Research Center

Hampton, VA 23681-2199

L-19610

National Aeronautics and Space Administration

Washington, DC 20546-0001

NASA

NASA/TM-2009-215715

Unclassified-Unlimited

Subject Category 26 Availability: NASA CASI (443) 757-5802

 Atomistic simulations of intergranular fracture have indicated that grain-scale crack growth in polycrystalline metals can be

direction dependent. At these material length scales, the atomic environment greatly influences the nature of intergranular crack

propagation, through either brittle or ductile mechanisms, that are a function of adjacent grain orientation and direction of crack

propagation. Methods have been developed to obtain cohesive zone models (CZM) directly from molecular dynamics simulation.

These CZMs may be incorporated into decohesion finite element formulations to simulate fracture at larger length scales. A new

directional decohesion element is presented that calculates the direction of Mode I opening and incorporates a material criterion for

dislocation emission based on the local crystallographic environment to automatically select the CZM that best represents crack

growth. The simulation of fracture in 2-D and 3-D aluminum polycrystals is used to illustrate the effectiveness of directional

decohesion finite elements.

Finite element method; Polycrystalline metals; Multiscale analysis

U U U UU 45

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

(443) 757-5802