1656_C012.

fm Page 563 Monday, May 23, 2005 6:07 PM

Computational Fracture Mechanics 563

(a) (b)

FIGURE 12.6 Virtual crack extension along a three-dimensional crack front: (a) uniform crack advance and
(b) advance over an increment of crack front.

In a three-dimensional problem, G typically varies along the crack front. In computing G, one
can consider a uniform virtual crack advance over the entire crack front or a crack advance over
a small increment, as Figure 12.6 illustrates. In the former case, ∆Ac = ∆aL, and the computed
energy release rate would be a weighted average. Defining Ac incrementally along the crack front
would result in a local measure of G.
For two-dimensional problems, the virtual crack extension formulation of G requires an area
integration, while three-dimensional problems require a volume integration. Such an approach is
easier to implement numerically and is more accurate than contour and surface integrations for
two- and three-dimensional problems, respectively.
Numerical implementation of the virtual crack extension method entails applying a virtual
displacement to nodes within a specified contour. Since the domain integral formulation is very
similar to the above method, further discussion on numerical implementation is deferred to
Section 12.3.3.

12.3 THE ENERGY DOMAIN INTEGRAL

Shih et al. [21,22] formulated the energy domain integral methodology, which is a general frame-
work for numerical analysis of the J integral. This approach is extremely versatile, as it can be
applied to both quasistatic and dynamic problems with elastic, plastic, or viscoplastic material
response, as well as thermal loading. Moreover, the domain integral formulation is relatively simple
to implement numerically, and it is very efficient. This approach is very similar to the virtual crack
extension method.

12.3.1 THEORETICAL BACKGROUND

Appendix 4.2 presents a derivation of a general expression for the J integral that includes the effects
of inertia as well as inelastic material behavior. The generalized definition of J requires that the
contour surrounding the crack tip be vanishingly small:

 ∂u j 
J = lim
Γo →0 ∫ ( w + T )δ1i − σ ij
Γο 
 n dΓ
∂x1  i
(12.23)

where T is the kinetic energy density. Various material behavior can be taken into account through
the definition of w, the stress work.
1656_C012.fm Page 564 Monday, May 23, 2005 6:07 PM

564 Fracture Mechanics: Fundamentals and Applications

Consider an elastic-plastic material loaded under quasistatic conditions (T = 0). If thermal

strains are present, the total strain is given by

ε ijtotal = ε ije + ε ijp + α Θδ ij = ε ijm + ε kkt (12.24)

where a is the coefficient of thermal expansion and Θ is the temperature relative to a strain-free
condition. The superscripts e, p, m, and t denote elastic, plastic, mechanical, and thermal strains,
respectively. The mechanical strain is equal to the sum of elastic and plastic components. The stress
work is given by
m
ε kl
w=
∫ 0
σ ij dε ijm (12.25)

The form of Equation (12.23) is not suitable for numerical analysis, since it is not feasible to
evaluate stresses and strains along a vanishingly small contour. Let us construct a closed contour
by connecting inner and outer contours, as Figure 12.7 illustrates. The outer contour Γ1 is finite,
while Γo is vanishingly small. For a linear or nonlinear elastic material under quasistatic conditions,
J could be evaluated along either Γ1 or Γo, but only the inner contour gives the correct value of J
in the general case. For quasistatic conditions, where T = 0, Equation (12.23) can be written in
terms of the following integral around the closed contour, Γ * = Γ1 + Γ+ + Γ− − Γo [21,22]:

 ∂u j  ∂u j
J=
∫ * σ
Γ
ij
∂x1
− wδ1i  qmi d Γ −

∫
Γ+ + Γ−
σ2j
∂x1
qd Γ (12.26)

where
Γ+ and Γ− = upper and lower crack faces, respectively
mi = outward normal on Γ *
q = arbitrary but smooth function that is equal to unity on Γo and zero on Γ1

Note that mi = − ni on Γo; also, m1 = 0 and m2 = ±1 on Γ+ and Γ−. In the absence of crack-face
tractions, the second integral in Equation (12.26) vanishes.

FIGURE 12.7 Inner and outer contours, which form a closed contour around the crack tip when connected
by Γ+ and Γ−.
1656_C012.fm Page 565 Monday, May 23, 2005 6:07 PM

Computational Fracture Mechanics 565

For the moment, assume that the crack faces are traction free. Applying the divergence theorem
to Equation (12.26) gives

∂  ∂u j  
J=
∫A* ∂xi
σ ij
 ∂x1
− wδ1i  q  dA
 
(12.27)
 ∂u j  ∂q  ∂  ∂u j  ∂w 
=
∫ σ ij
A*  ∂x1
− wδ1i 
 ∂xi
dA + 
A*  ∂xi

∫  σ ij ∂x  − ∂x  qdA
 1 1


where A* is the area enclosed by Γ*. Referring to Appendix 3.2, we see that

∂  ∂u j  ∂w
σ − =0
∂xi  ij ∂x1  ∂x1
(12.28)

when there are no body forces and w exhibits the properties of an elastic potential:

∂w
σ ij = (12.29)
∂ε ij

It is convenient at this point to divide w into elastic and plastic components:

e
ε kl ε klp
w = we + w p =
∫ 0
σ ij dε ije +
∫
0
Sij dε ijp (12.30)

where Sij is the deviatoric stress, defined in Equation (A3.62). While the elastic components of w
and eij satisfy Equation (12.29), plastic deformation does not, in general, exhibit the properties of
a potential. (Equation (12.29) may be approximately valid for plastic deformation when there is
no unloading.) Moreover, thermal strains would cause the left side of Equation (12.28) to be
nonzero. Thus the second integrand in Equation (12.27) vanishes in limited circumstances, but not
in general. Taking account of plastic strain, thermal strain, body forces, and crack face-tractions
leads to the following general expression for J in two dimensions:

 ∂u j  ∂q  ∂ε ijp ∂w p ∂Θ ∂u  
J=
∫ σ ij
A*  
 ∂x1
− wδ1i  + σ ij −
 ∂xi  ∂x1 ∂x1
+ ασ ii
∂x1
− Fi j  q dA
∂x1  
∂u
−
∫ Γ+ + Γ−
σ 2 j j qdΓ
∂x1
(12.31)

where the body force contribution is inferred from the equilibrium equations, and the contribution
from thermal loading is obtained by substituting Equation (12.24) and Equation (12.30) into
Equation (12.27). Inertia can be taken into account by incorporating T, the kinetic energy density,
into the group of terms that are multiplied by q. For a linear or nonlinear elastic material under
quasistatic conditions, in the absence of body forces, thermal strains, and crack-face tractions,
Equation (12.31) reduces to

 ∂u j  ∂q
J=
∫ σ ij
A*  ∂x1
− wδ1i 
 ∂xi
dA (12.32)

Equation (12.32) is equivalent to Rice’s path-independent J integral (Chapter 3). When the sum of
the additional terms in the more general expression (Equation 12.31) is nonzero, J is path dependent.
1656_C012.fm Page 566 Monday, May 23, 2005 6:07 PM

566 Fracture Mechanics: Fundamentals and Applications

FIGURE 12.8 Surface enclosing an increment of a three-dimensional crack front.

Comparing Equation (12.21) and Equation (12.32) we see that the two expressions are identical
if q = ∆x1/∆a. Thus q can be interpreted as a normalized virtual displacement, although the above
derivation does not require such an interpretation. The q function is merely a mathematical device
that enables the generation of an area integral, which is better suited to numerical calculations.
Section 12.3.3 provides guidelines for defining q.

12.3.2 GENERALIZATION TO THREE DIMENSIONS

Equation (12.23) defines the J integral in both two and three dimensions, but the form of this
equation is poorly suited to numerical analysis. In the previous section, J was expressed in terms
of an area integral in order to facilitate numerical evaluation. For three-dimensional problems, it
is necessary to convert Equation (12.23) into a volume integral.
Figure 12.8 illustrates a planar crack in a three-dimensional body; h corresponds to the position
along the crack front. Suppose that we wish to evaluate J at a particular h on the crack front. It is
convenient to define a local coordinate system at h, with x1 normal to the crack front, x2 normal
to the crack plane, and x3 tangent to the crack front. The J integral at h is defined by Equation
(12.23), where the contour Γo lies in the x1–x2 plane.
Let us now construct a tube of length ∆L and radius ro that surrounds a segment of the crack
front, as Figure 12.8 illustrates. Assuming quasistatic conditions, we can define a weighted average
J over the crack front segment ∆L as follows:

J∆ L =
∫∆L
J (η)qdη

 ∂u j 
= lim
rο →0
∫ wδ1i − σ ij
Sο 
 qn dS
∂x1  i
(12.33)

where
J(h) = point-wise value of J
So = surface area of the tube in Figure 12.8
q = weighting function that was introduced in the previous section
1656_C012.fm Page 567 Monday, May 23, 2005 6:07 PM

Computational Fracture Mechanics 567

FIGURE 12.9 Interpretation of q in terms of a virtual crack advance along ∆L.

Note that the integrand in Equation (12.33) is evaluated in terms of the local coordinate system,
where x3 is tangent to the crack front at each point along ∆L.
Recall from the previous section that q can be interpreted as a virtual crack advance. For
example, Figure 12.9 illustrates an incremental crack advance over ∆L, where q is defined by

∆a(η) = q(η)∆amax (12.34)

and the incremental area of the virtual crack advance is given by

∆ Ac = ∆amax
∫∆L
q(η)dη (12.35)

The q function need not be defined in terms of a virtual crack extension, but attaching a physical
significance to this parameter may aid in understanding.
If we construct a second tube of radius r1 around the crack front (Figure 12.10), it is possible
to define the weighted average J in terms of a closed surface, analogous to the two-dimensional
case (Figure 12.7 and Equation 12.26):

 ∂u j  ∂u j
J ∆L =
∫
S*
σ ij
 ∂x1
− wδ1i  qmi dS −

∫
S+ + S−
σ2j
∂x1
qdS (12.36)

FIGURE 12.10 Inner and outer surfaces, So and S1, which enclose V *.
1656_C012.fm Page 568 Monday, May 23, 2005 6:07 PM

568 Fracture Mechanics: Fundamentals and Applications

where the closed surface S* = S1 + S+ + S− − So , and S+ and S− are the upper and lower crack faces,
respectively, that are enclosed by S1. From this point, the derivation of the domain integral formu-
lation is essentially identical to the two-dimensional case, except that Equation (12.34) becomes a
volume integral:

 ∂u  ∂q  ∂ε ij p ∂w p ∂Θ ∂u  
∫
j
J ∆L = σ ij − wδ1i  + σ ij − + ασ ii − Fi j  q  dV
V *  ∂x1  ∂xi  ∂x1 ∂x1 ∂x1 ∂x1  
 
(12.37)
∂u j
−
∫ S+ + S
−
σ2j
∂x1
qd Γ

Equation (12.37) requires that q = 0 at either end of ∆L; otherwise, there may be a contribution to
J from the end surfaces of the cylinder. The virtual crack advance interpretation of q (Figure 12.9)
fulfills this requirement.
If the point-wise value of the J integral does not vary appreciably over ∆L, to a first approxi-
mation, J(h) is given by

J ∆L
J (η) ≈ (12.38)
∫ ∆L q(η, ro )dη

Equation (12.38) is a reasonable approximation if the q gradient along the crack front is steep relative
to the variation in J(h).
Recall that Equation (12.22) was defined in terms of a fixed coordinate system, while Equation
(12.37) assumes a local coordinate system. The domain integral formulation can be expressed in
terms of a fixed coordinate system by replacing q with a vector quantity qi, and evaluating the
partial derivatives in the integrand with respect to xi rather than x1, where the vectors qi and xi are
parallel to the direction of crack growth. Several commercial codes that incorporate the domain
integral definition of J require that the q function be defined with respect to a fixed origin.

12.3.3 FINITE ELEMENT IMPLEMENTATION

Shih et al. [21] and Dodds and Vargas [23] give detailed instructions for implementing the domain
integral approach. Their recommendations are summarized briefly.
In two-dimensional problems, one must define the area over which the integration is to be
performed. The inner contour Γo is often taken as the crack tip, in which case A* corresponds to
the area inside of Γ1. The boundary of Γ1 should coincide with element boundaries. An analogous
situation applies in three dimensions, where it is necessary to define the volume of integration. The
latter situation is somewhat more complicated, however, since J(h) is usually evaluated at a number
of locations along the crack front.
The q function must be specified at all nodes within the area or volume of integration. The shape
of the q function is arbitrary, as long as q has the correct values on the domain boundaries. In a plane
stress or plane strain problem, for example, q = 1 at Γo , which is usually the crack tip, and q = 0 at
the outer boundary. Figure 12.11 illustrates two common examples of q functions for two-dimensional
problems, with the corresponding virtual nodal displacements. This example shows 4-node square
elements and rectangular domains for the sake of simplicity. The pyramid function (Figure 12.11(a))
is equal to 1 at the crack tip but varies linearly to zero in all directions, while the plateau function
(Figure 12.11(b)) equals 1 in all regions except the outer ring of elements. Shih et al. [21] have shown
that the computed value of J is insensitive to the assumed shape of the q function.
1656_C012.fm Page 569 Monday, May 23, 2005 6:07 PM

Computational Fracture Mechanics 569

(a)

(b)

FIGURE 12.11 Examples of q functions in two dimensions, with the corresponding virtual nodal displace-
ments: (a) the pyramid function and (b) the plateau function. Taken from Shih, C.F., Moran, B., and Nakamura,
T., ‘‘Energy Release Rate Along a Three-Dimensional Crack Front in a Thermally Stressed Body.” International
Journal of Fracture, Vol. 30, 1986, pp. 79–102.

Figure 12.12 illustrates the pyramid function along a three-dimensional crack front, where the
crack-tip node of interest is displaced a unit amount, and all other nodes are fixed. If desired, J(h)
can be evaluated at each node along the crack front.
The value of q within an element can be interpolated as follows:
n

q( xi ) = ∑N q
I =1
I I (12.39)

FIGURE 12.12 Definition of q in terms of a virtual nodal displacement along a three-dimensional crack front.
1656_C012.fm Page 570 Monday, May 23, 2005 6:07 PM

570 Fracture Mechanics: Fundamentals and Applications

where
n = number of nodes per element
qI = nodal values of q
NI = element shape functions, which were introduced in Section 12.1.1

The spatial derivatives of q are given by

n 2 or 3
∂q ∂N I ∂ξk
∂xi
= ∑ ∑ ∂ξ
I =1 k =1 k ∂x j
qI (12.40)

where xi are the parametric coordinates for the element.

In the absence of thermal strains, path-dependent plastic strains, and body forces within the
integration volume or area, the discretized form of the domain integral is as follows:

m
 ∂u j  ∂q   ∂x j  
J= ∑ ∑ σ
A* or V * p =1
ij
∂x1
− wδ1i   det 
 ∂xi 
  wp
 ∂ξk  
p

 ∂u j 
− ∑
crack faces
 σ 2 j ∂x
 1
q w

(12.41)

where m is the number of Gaussian points per element, and wp and w are weighting factors. The
quantities within { }p are evaluated at the Gaussian points. Note that the integration over crack
faces is necessary only when there are nonzero tractions.

12.4 MESH DESIGN

The design of a finite element mesh is as much an art as it is a science. Although many commercial
codes have automatic mesh generation capabilities, construction of a properly designed finite
element model invariably requires some human intervention. Crack problems, in particular, require
a certain amount of judgment on the part of the user.
This section gives a brief overview of some of the considerations that should govern the
construction of a mesh for analysis of crack problems. It is not possible to address in a few pages
all of the situations that may arise. Readers with limited experience in this area should consult
the published literature, which contains numerous examples of finite element meshes for crack
problems.
Figure 12.13 shows examples of the arrangement of nodes on the crack plane of two-dimensional
models. When both crack faces are modeled, there are usually matching nodes along each crack
face. If these matching node pairs have identical coordinates, many commercial mesh generation
tools will merge these node pairs into single nodes, which will cause the crack to heal. One technique
for avoiding this unwanted intervention by commercial meshing tools is to create a small gap
between the crack faces. The width of this gap should be much less then the crack length. When
the crack is on a symmetry plane, which is the case for most Mode I problems, only one crack face
needs to be modeled. A symmetry constraint is applied to the uncracked cross section, and the crack
face is normally unconstrained.2

2 For some problems, the crack face may have a traction boundary condition or an imposed displacement field.