INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2000; 48:1549{1570
Extended nite element method for
three-dimensional crack modelling
N. Sukumarz, N. Moesx, B. Moran{ and T. Belytschko;y;k
Department of Civil and Mechanical Engineering; Northwestern University; 2145 Sheridan Road;
Evanston; IL 60208; U.S.A.
SUMMARY
An extended nite element method (X-FEM) for three-dimensional crack modelling is described. A discontin-
uous function and the two-dimensional asymptotic crack-tip displacement elds are added to the nite element
approximation to account for the crack using the notion of partition of unity. This enables the domain to be
modelled by nite elements with no explicit meshing of the crack surfaces. Computational geometry issues
associated with the representation of the crack and the enrichment of the nite element approximation are
discussed. Stress intensity factors (SIFs) for planar three-dimensional cracks are presented, which are found
to be in good agreement with benchmark solutions. Copyright ? 2000 John Wiley & Sons, Ltd.
KEY WORDS: extended nite element method; partition of unity; local enrichment; elastostatics; planar three-
dimensional cracks
1. INTRODUCTION
The accurate modelling of three-dimensional cracks in nite bodies remains a challenging prob-
lem in computational mechanics. The relevance and importance of the computation of fracture
parameters and the simulation of three-dimensional crack growth stems from the widespread use
of numerical fracture mechanics in fatigue life predictions of safety-critical components such as
aircraft fuselages, pressure vessels, automobile components, and castings. Fatigue failure usually
Correspondence to: T. Belytschko, Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road,
Evanston, IL 60208, U.S.A.
yE-Mail: t-belytschko@nwu.edu
zPost-Doctoral Research Fellow, Theoretical and Applied Mechanics
xResearch Associate, Department of Mechanical Engineering
{Associate Professor of Civil Engineering
kWalter P. Murphy, Professor of Computational Mechanics
Contract=grant sponsor: National Science Foundation; contract=grant number: CMS-9732319
Contract=grant sponsor: Georgia Institute of Technology ONR; contract=grant number: N0014-95-1-0539
Contract=grant sponsor: Federal Aviation Administration; contract=grant number: DTFA03-98-F-IA025
Contract=grant sponsor: Oce of Naval Research
Received 15 September 1999
Copyright
?
2000 John Wiley & Sons, Ltd. Revised 29 November 1999

1550 N. SUKUMAR ET AL.
occurs due to the initiation and propagation of surface or near-surface cracks, which are often
assumed to be elliptical or semi-elliptical in shape for numerical modelling. Closed-form solutions
for the stress intensity factors (SIFs) are available for simple crack geometries in three dimensions;
however, for arbitrary-shaped cracks in nite specimens, numerical methods are the only recourse
to modelling three-dimensional fatigue crack growth.
Currently, nite element methods based on non-singular as well as singular elements are widely
used in linear elastic fracture analysis. Methods based on singular elements [1{3] are able to pro-
vide accurate stress intensity factors. Enriched element methods which incorporate the asymptotic
crack-tip elds in the trial functions [4; 5] provide stress intensity factors directly as part of the
solution. In spite of the successes using nite elements in computational fracture, mesh generation
in three dimensions is time consuming and especially burdensome for multiple crack congurations
and crack growth simulations. It is dicult to explicitly model the crack topology as part of the
nite element since accuracy considerations require signicant renement in the vicinity of the
crack front.
The modelling of growing cracks using nite elements and boundary elements with re-meshing
has been pursued by many researchers [6{8]. The nite element alternating method which is based
on the Schwartz{Neumann alternating method has been applied to the modelling of arbitrary crack
congurations in three dimensions [9; 10]. Numerical techniques such as the boundary integral
equation method [11], body force method [12], and the self-similar crack expansion method [13]
also enable the accurate computation of stress intensity factors for three-dimensional cracks. The
requirement of a Green’s function imposes restriction on the scope of application of the boundary
integral equation based methods. These methods are not readily extendable to non-linear problems
nor to the study of cracks in anisotropic materials. The modelling of cracks as a continuous distri-
bution of innitesimal dislocation loops has also gained prominence. In the dislocation distribution
approach, the crack surface in three dimensions is assumed to be made up of a continuous dis-
tribution of innitesimal dislocation loops with the Burger’s vector corresponding to the jump in
displacement across the crack surface; a hypersingular integral equation is obtained on the crack
surface which is solved in the Hadamard nite-part sense [14; 15]. More recently, the element-free
Galerkin (EFG) method, a mesh-free method, has been successful in modelling static and dynamic
fracture in two dimensions [16; 17] and three dimensions [18; 19].
The extended nite element method (X-FEM) alleviates shortcomings associated with meshing
of the crack surfaces in existing methods. The nite element method is used as the building block
in the extended nite element method, and hence much of the theoretical and numerical devel-
opments in nite elements can be readily extended and applied. In this paper, the nite element
approximation is enriched by additional functions through the notion of partition of unity [20].
Strouboulis et al. [21] used local enrichment functions within the partition of unity framework
for modelling re-entrant corners in two dimensions, and in Reference [22], enrichment functions
for holes are proposed. In X-FEM, the additional functions are used to model the presence of
cracks, voids or inhomogeneities, and also to improve accuracy in problems where some aspects
of the functional behaviour of the solution eld is known a priori. For crack modelling, a dis-
continuous function and the two-dimensional asymptotic crack-tip displacement elds are added
to the displacement-based nite element approximation. Partition of unity enrichment methods for
discontinuities and near-tip crack elds were introduced by Belytschko and Black [23]. Moes
et al. [24] proposed the generalized Heaviside function as a means to model the crack away from
the crack tip, and developed simple rules for the introduction of the discontinuous and crack-tip
enrichments. The numerical results presented in this paper demonstrate that for three-dimensional
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1551
cracked bodies, the accuracy and performance of X-FEM is excellent. This opens up many exciting
possibilities for its further development and advancement.
The outline of this paper follows. In the following section, we introduce the extended nite
element method, and in Section 3, we describe crack modelling in R2 and R3. In Section 4, the
computational methodology adopted in X-FEM and the numerical issues involved are discussed.
The strong and weak forms of the Galerkin method for elastostatics, along with the discrete
equations for X-FEM are given in Section 5. In Section 6, a brief description of three-dimensional
domain integral computations is presented, which is followed by numerical results for several
problems in three-dimensional linear elastic fracture mechanics. The numerical SIF results are
compared to available reference solutions from the literature. Some nal remarks and conclusions
are made in Section 7.
2. EXTENDED FINITE ELEMENT METHOD
In nite element methods, a partition of a domain into sub-domains (elements) forms the basis of
mesh generation. The presence of
aws or inhomogeneities such as cracks, voids, and inclusions
must be taken into account in mesh generation—the mesh must conform to these geometric enti-
ties. In unstructured mesh generation, a Delaunay tessellation of the domain results in triangular
elements in R2 and tetrahedral elements in R3. Alternatively, algorithms for the construction of
quadrilateral (in two dimensions) or hexahedral meshes (in three dimensions) are also available.
In computational fracture, quadrilateral or hexahedral elements are preferred to triangular or tetra-
hedral elements due to the higher-order approximation in the former. However, in comparison
to unstructured mesh generation, hexahedral mesh generation has not reached a mature stage of
development, and hence mesh generation in three dimensions for computational fracture applica-
tions continues to present many challenges.
The extended nite element method alleviates much of the burden associated with mesh gener-
ation by not requiring the nite element mesh to conform to cracks, and in addition, provides a
seamless means to use higher-order elements or special nite elements without signicant changes
in the formulation. The essence of X-FEM lies in sub-dividing a model problem into two distinct
parts: mesh generation for the geometric domain (cracks not included), and enriching the nite
element approximation by additional functions that model the
aw(s) and other geometric entities.
The enrichment of the nite element approximation is described as follows. Consider a point
x of Rd (d = 1{3) that lies inside a nite element e. Denote the nodal set N = {n1; n2; : : : ; nm},
where m is the number of nodes of element e. (m = 2 for a linear one dimensional nite element,
m = 3 for a constant-strain triangle, m = 8 for a trilinear hexahedral element, etc.) The enriched
displacement approximation for a vector-valued function u(x) :Rd →Rd assumes the form
uh(x) =
P
I
nI∈N
I (x)uI +
P
J
nJ∈N
g
J (x) (x)aJ ; (uI ; aJ ∈ R
d) (1)
where the nodal set Ng is dened as
Ng= {nJ : nJ ∈N; !J ∩
g 6= ∅} (2)
In the above equation, !J =supp(nJ ) is the support of the nodal shape function J (x), which
consists of the union of all elements with nJ as one of its vertices; and
g is the domain associated
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1552 N. SUKUMAR ET AL.
Figure 1. Co-ordinate conguration for crack front enrich-
ment functions.
Figure 2. Partitioning algorithm for a crack in
two dimensions.
with a geometric entity such as a hole, crack surface, or crack front. In general, the choice of the
function (x) depends on the geometric entity.
3. CRACK MODELLING
The crack is modelled by enriching the nodes whose nodal shape function support intersects the
interior of the crack by a discontinuous function, and enriching the nodes whose nodal shape
function support intersects the crack front (crack-tip in R2) by the two-dimensional asymptotic
crack-tip elds.
3.1. Enrichment
We restrict the decription of the implementation to planar cracks in three dimensions. The
enrichments concepts that follow are easily extended to non-planar cracks but the implementa-
tion is more dicult.
Consider a single crack in three dimensions, and let c be the crack surface and c the crack
front. Note that for an internal crack, the crack front corresponds to the boundary of the crack:
c = @c whereas for an edge crack, the crack front is only part of the boundary: c⊂ @c. The
interior of a planar crack is modelled by the enrichment function H (x), which we refer to as a
generalized Heaviside function. The function H (x) takes on the value +1 above the crack and −1
below the crack. More precisely, let x be the closest point to x on the crack c, and n be the
normal to the crack plane (Figure 1). The H (x) function is then given by +1 if (x − x) · n>0
and −1 otherwise, i.e.
H (x)=
(
1 if (x− x) · n>0
−1 otherwise
(3)
To model the crack front and also to improve the representation of crack-tip elds in three-
dimensional computations, crack-tip enrichment functions are used in elements which contain the
crack front. In the neighbourhood of the crack front, the asymptotic elds are two-dimensional in
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1553
nature. The crack front enrichment consists of functions which incorporate the radial and angular
behaviour of the two-dimensional asymptotic crack-tip displacement eld:
(x)≡{ 1; 2; 3; 4}=

√
r cos

2
;
√
r sin

2
;
√
r sin  sin

2
;
√
r sin  cos

2

(4)
where r and  are polar co-ordinates in the x1ˆ{x2ˆ plane (Figure. 1). Note that the second function
in the above equation is discontinuous on the crack plane.
If a three-dimensional planar crack coincides with element boundaries, an equivalence exists
between the nite element space with the crack explicitly modelled by the mesh and the X-FEM
space with H (x) used as the enrichment function. This correspondence for two-dimensional cracks
was shown in References [24; 25], which readily extends to three dimensions.
4. COMPUTATIONAL METHODOLOGY
Computational geometry issues associated with the representation of the crack and the enrichment
of the nite element approximation are discussed. The X-FEM implementation is performed in
C++: the object-oriented features of C++ provide for better data management, encapsulation, code
re-usability,
exibility, and maintenance. In the following sub-section, we describe some of the im-
portant computational geometry issues, with an aim towards an accurate and robust implementation.
4.1. Mesh–geometry interactions
The nite element mesh consists of tetrahedral, prismatic, and hexahedral elements. In two dimen-
sions, cracks are represented by line segments and the crack front is a point. In three dimensions,
a crack is represented by a polygon partition into triangles, and the crack front consists of line
segments as in Reference [19].
The notion of geometric predicates which is widely used in computational geometry, is also
an integral part of the X-FEM mesh{geometric computations. In the following sub-section, we
brie
y touch upon the concept of geometric predicates. In crack modelling, dierent functions are
used to enrich the displacement approximation in the interior of the crack and on the crack front.
Therefore, one of the rst tasks is to determine the nite elements that intersect the crack. These
nite elements are then partitioned into simplices so that numerical integration of the weak form
accounts for the discontinuities on either side of the crack surface. These sub-tetrahedrons are also
used to determine whether a node is to be enriched (see Section 4.1.4).
4.1.1. Geometric predicates. Since data is stored and computations performed using nite-precision
arithmetic, it is essential that the robustness of algorithms is maintained even for small perturba-
tions in the data. Eorts have been made to develop robust geometric predicates [26], which are
especially important in the development of algorithms for the Delaunay tessellation and Voronoi
diagram of a point set.
The incircle and orientation tests are widely used in computational geometry. The orientation
test determines whether a point lies to the left of, to the right of, or on a line or plane dened
by other points. The incircle test determines whether a point lies inside, outside, or on a circle
(sphere) dened by other points. Each of these tests is performed by evaluating the sign of a
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1554 N. SUKUMAR ET AL.
determinant. Instead of explicitly computing intersections, geometric predicates provide an easier
and robust means to evaluate queries.
In X-FEM, the orientation test is used to determine the nodal enrichment for a query point x and
in the algorithm for the partitioning of the nite elements into tetrahedrons in three dimensions.
In the evaluation of the predicates associated with a planar crack and for generality to both two
and three dimensions, we assume the crack plane (crack segment in two dimensions) is dened
by a point x0 that lies on it and a unit normal m. In two dimensions, if e1 and e2 are unit vectors
in the plane, then t×m is a vector along e3, where t is the tangent vector to the crack segment,
and in three dimensions m is the unit normal to the plane of the crack. The equation of the crack
segment in R2 or the crack plane in R3 can be written as
f(x)= (x− x0) ·m=0 (5)
For a query point x, the ternary predicates ABOVE (f(x)¿), BELOW (f(x)¡−), and ON
(−6f(x)6) are used, where =10−6 is used as the tolerance for a nite element mesh
with element edge length of O(1).
4.1.2. Crack–mesh intersection. A nave approach to nding all nite elements that intersect the
crack is to test for intersection with each element in the mesh. Clearly, this is of O(N ) complexity,
where N is the number of elements in the mesh, and hence is not a computationally attractive
choice. A more appealing alternative that is computationally feasible is adopted. Since determining
the elements that intersect with the crack is non-trivial, we rst nd a candidate set that is possibly
larger than the desired set. We use a bounding box (BB) of cell dimensions 10×10×10 that
encompasses the entire nite element mesh. The cell dimensions of the BB-mesh is used as a
compromise between computational costs and storage requirements. For a given query point x that
is contained in a BB-element b, all the associated nite elements tb that intersect b are obtained. In
three dimensions we use a map of the vertices of a triangle from x→ [; ], where (0; 0); (1; 0),
and (0; 1) are the vertex co-ordinates of the reference right-angled triangle. The starting search
path for a triangle [xa; xb; xc] is (; )= (0; 0); (; )= (1; 0), and (; )= (0; 1), In order to nd
the set of all intersected elements (may be larger than the exact set), a recursive algorithm is
implemented. The stopping criterion is
1. if all points are in the same nite element; or
2. if the element sets for all adjacent points have at least one element in common; or
3. if the distance between all points in the reference co-ordinate is less than a prescribed tol-
erance. In three dimensions, |1 − 2|¡10−5=
√
V and |1 − 2|¡10−5=
√
V , where V is the
area of the triangle (simplex partition of the crack).
4.1.3. Partitioning of nite elements. The algorithm to compute the intersected nite elements is
outlined in Section 4.1.2. Here we describe the algorithm to partition a given three-dimensional
nite element into tetrahedrons. The partitioning is carried out for all the elements that are selected
on the basis of the algorithm presented in Section 4.1.2; as mentioned earlier, this set of elements
can be larger than the exact set of elements that intersect the crack domain. In Section 4.1.4, an
additional criterion is invoked so that the nodal enrichment is precise.
We rst illustrate the procedure in two dimensions. Each crack segment is dened by a point
x0 on it and the unit normal vector m that is dened in the previous sub-section. The face and
element connectivity of each nite element are also known, with counterclockwise orientation of
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1555
the nodes assumed in the connectivity. Standard Template Library (STL) containers in C++ are
used to create maps for node-to-co-ordinates, node-to-orientation, and edge-to-connectivity for the
nite element. The orientation (ABOVE= + 1; BELOW= − 1; ON=0) of the nodes of the nite
element are set and stored in the node-to-orientation map.
The element is partitioned if and only if there exists two nodes in the connectivity with orien-
tation +1 and −1. A loop over the edges of the nite element is carried out, and if there is an
intersection for an edge, the point of intersection is computed and added to the node-to-co-ordinate
and edge-to-connectivity maps. The orientation of the point of intersection is set to +10 if the
orientation of the rst point of the edge is +1, and −10 if the orientation of the rst point is −1.
This updates the edge map for the nite element. Now, using the edge map of the nite element,
two new edge maps for the domains above and below the crack are created. The two edges (upper
and lower) that correspond to the crack segment are also created such that counterclockwise orien-
tation of the nodes is preserved. The centroid of the two domains is evaluated with the orientation
for the node that belongs to the surface above the crack set to +100, and the orientation set to
−100 for the node that belongs to the surface below the crack. A loop over the two edge maps is
carried out and on using the centroid for each surface, the partitioned elements are obtained. The
partitioned elements for a sample crack are illustrated in Figure 2.
The partitioning of three-dimensional polyhedrons is more involved than the two-dimensional
case. Even for the simplied case of partitioning a three-dimensional nite element that is cut
by a plane, no public-domain package that can readily interface to a C++ program is currently
available. An algorithm was hence developed and implemented for the partitioning of a three-
dimensional nite element (such as tetrahedron, prism, or hexahedron) by a single crack plane.
The crack is dened by a point x0 that lies on it and the unit normal vector m that is dened in the
previous sub-section. The nodal connectivities of the nite element are stored in counterclockwise
orientation, and the face connectivities are stored such that an edge that is common to two faces
is traversed in opposite directions. The nodal and face connectivity of a tetrahedron are such that
if the nodes on a face are traversed in a certain sense, then the thumb, using the right-hand rule,
points in the direction of the fourth node. Standard Template Library (STL) containers are used
to create maps for node-to-co-ordinates, node-to-orientation, and face-to-connectivity for the nite
element. The orientation (ABOVE = +1; BELOW = −1, ON =0) of the nodes of the nite element
are set and stored in the node-to-orientation map.
The element is partitioned if there exists two nodes in the connectivity with orientation +1 and
−1. A loop over the faces of the nite element is carried out, and if the crack plane intersects
an edge, the point of intersection is computed and added to the node-to-co-ordinate and face-to-
connectivity maps. The orientation of the point of intersection is set to +10 if the orientation of
the rst point on the edge is +1, and −10 if the orientation of the rst point is −1. This updates
the face-to-connectivity map for the nite element.
Next, using the face map of the nite element, two new maps for the domains above and
below the crack plane are created. The two faces (upper and lower) that correspond to the crack
plane are also created so that the orientation of the nodes in the face meet the requirement—
counterclockwise for nodes in the above surface and clockwise orientation for the nodes in the
below surface preserve the correct orientation of the partitioned tetrahedrons. The centroid of the
nodes in the face maps is evaluated and the orientation of the new point (node) that belong to the
above face map is set to +100, and −100 for the point (node) below the face map. The centroid
of the two volumes is evaluated with orientation +1000 for the node that belongs to the volume
above the crack and orientation set to −1000 for the node that belongs to the volume below the
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1556 N. SUKUMAR ET AL.
Figure 3. Partitioning algorithm for an edge
crack in three dimensions.
Figure 4. Elastostatic boundary value problem.
crack. A loop over the two face maps is carried out and using the centroid for each face and the
centroid of the two volumes, the partitioned elements are obtained. The partitioned elements for
an edge crack in a 4× 4× 4 hexahedral mesh are shown in Figure 3.
4.1.4. Nodal enrichment. We next describe the enrichment for crack modelling. The enriched
nite element approximation is
uh(x)=
P
I
nI∈N
I (x)uI +
P
J
nJ∈N
c
J (x)H (x)aJ +
P
K
nK∈N
f
K (x)

4
X
l=1
l(x)b
l
K
!
(6)
The second and third terms on the right-hand side of the above equation are the discontinuity and
front enrichments, respectively. The set Nf consists of those nodes for which the closure of the
nodal shape function support intersects the crack front. The set Nc is the set of nodes whose nodal
shape function support is intersected by the crack and which do not belong to Nf :
Nf = {nK : nK ∈N; !K ∩c 6= ∅} (7)
Nc = {nJ : nJ ∈N; !J ∩c 6= ∅; nJ =∈N
f
} (8)
Note that for any node in Nc, the support of the nodal shape function is fully cut into two
disjoint pieces by the crack. If for a certain node nI , one of the two pieces is very small compared
to the other, then the generalized Heaviside function used for the enrichment is almost a constant
over the support, leading to an ill-conditioned stiness matrix. Therefore, in this case, node nI is
removed from the set Nc. The criterion for nodal inclusion in Nc is as follows: the volume above
the crack is Vabove, and the volume below the crack is V below! :V!=V
above
! +V
below
! . If either of the
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1557
two ratios, V above! =V! or V
below
! =V! is below a prescribed tolerance, the node is removed from the
set N c. We use a tolerance =10−4.
5. GOVERNING EQUATIONS
5.1. Strong form
Consider a body
⊂R3, with boundary . The boundary consists of the sets u, t, and ic , such
that =u ∪t ∪mi=1
i
c . All the internal surfaces
i
c are assumed to be traction-free. The boundary
value problem of elastostatics solves for the displacement u(x) of a body
which is xed on
u and subjected to surface forces (tractions) along t (Figure 4). We now consider the boundary
value problem for small displacement elastostatics. The eld equations of elastostatics are
r · b + b= 0 in
(9a)
b=C : U (9b)
U=rsu (9c)
where rs is the symmetric gradient operator and C is the tensor of elastic moduli for a homoge-
neous isotropic material.
The essential and natural boundary conditions are
u= u on u (10a)
b · n= t on t (10b)
b · n= 0 on ic (i=1; 2; : : : ; m) (10c)
where n is the unit outward normal to
, u and t are prescribed displacements and tractions,
respectively, and m is the number of internal surfaces. Note that Equation (10c) imposes the
condition that the internal surfaces ic be traction-free.
5.2. Weak form and discrete system
We consider the weak form for the equilibrium equation of elasticity with the associated boundary
conditions. For the discrete system, the weak form (principle of virtual work) is
Find uh ∈Vh such that
Z

h
b(uh):U(vh)=
Z

h
b · v d
+
Z
ht
t · v d ∀vh ∈Vh0 (11)
where uh(x)∈Vh and vh(x)∈Vh0 are the approximating trial and test functions used in X-FEM.
Since we are not concerned with convergence proofs, we restrict the discussion to those spaces
used in the construction of the discrete approximation. The space Vh is the enriched nite element
space that satisfy the Dirichlet boundary conditions, and which include basis functions that are
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1558 N. SUKUMAR ET AL.
discontinuous across the crack surfaces. The space Vh0 is the corresponding space with homogeneous
Dirichlet boundary conditions.
In a Bubnov{Galerkin procedure, the trial functions uh as well as the test functions vh are
represented as linear combinations of the same shape functions. The trial and test functions, which
are based on Equation (6) are
uh(x) =
P
I
nI∈N
I (x)uI +
P
J
nJ∈N
c
J (x)H (x)aJ +
P
K
nK∈N
f
K (x)

4
X
l=1
l(x)b
l
K
!
(12)
vh(x) =
P
I
nI∈N
I (x)vI +
P
J
nJ∈N
c
J (x)H (x)cJ +
P
K
nK∈N
f
K (x)

4
X
l=1
l(x)e
l
K
!
(13)
where I (x) are the nite element shape functions, and j(x) (j=1{4) are the enriched functions
for the crack front, which are given in Equation (4).
On substituting the trial and test functions from Equation (12) into Equation (11), and using
the arbitrariness of nodal variations, the following discrete system of linear equations is obtained:
Kd= f (14)
where d is the vector of unknowns at the nodes, and the stiness matrix K is
KI J =
2
6
4
KuuI J K
ua
I J K
ub
I J
KauI J K
aa
I J K
ab
I J
KbuI J K
ba
I J K
bb
I J
3
7
5
(15a)
K

I J =
Z

h
(BI )
T
CB

J d
(; = u; a; b) (15b)
In Equation (14), the external force vector f is dened as
fI = {f
u
I ; f
a
I ; f
b
I1; f
b
I2; f
b
I3; f
b
I4} (16a)
fuI =
Z
ht
I t d +
Z

h
Ib d
(16b)
faI =
Z
ht
IH t d +
Z

h
IHb d
(16c)
fbIj =
Z
ht
I jt d +
Z

h
I jb d
(j=1{4) (16d)
The discrete equations have three degrees of freedom for unenriched nodes. Nodes in the set
Nc have each three degrees of freedom, and nodes in the set Nf have each twelve degrees of
freedom—see Section 4.1.4 for details on the nodal sets Nc and Nf .
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1559
In Equation (15), C is the constitutive matrix for an isotropic linear elastic material, and BuI ,
BaI , and B
b
Ij are the matrix of shape function derivatives which are given by
BuI =
2
6
6
6
6
6
6
6
4
I; x 0 0
0 I;y 0
0 0 I; z
0 I; z I;y
I; z 0 I; x
I;y I; x 0
3
7
7
7
7
7
7
7
5
(17a)
BaI =
2
6
6
6
6
6
6
6
4
(IH); x 0 0
0 (IH);y 0
0 0 (IH); z
0 (IH); z (IH);y
(IH); z 0 (IH); x
(IH);y (IH); x 0
3
7
7
7
7
7
7
7
5
(17b)
BbI = [B
b
I1 B
b
I2 B
b
I3 B
b
I4] (17c)
BbIj =
2
6
6
6
6
6
6
6
4
(I j); x 0 0
0 (I j);y 0
0 0 (I j); z
0 (I j); z (I j);y
(I j); z 0 (I j); x
(I j);y (I j); x 0
3
7
7
7
7
7
7
7
5
(j=1{4) (17d)
6. NUMERICAL RESULTS
Several problems are presented to illustrate the accuracy and versatility of the extended nite
element method in three-dimensional elastostatics. We rst solve some benchmark crack problems,
and then examine the performance of X-FEM for cracks in nite specimens. Finite element mesh
generation is carried out by gmsh [27], which is available in the public-domain. In all problems,
numerical integration is carried out using Gauss{Legendre quadrature. In hexahedral elements
associated with only the nite element shape functions, 2×2×2 quadrature is used, and in elements
that also have enriched degrees of freedom, 6 × 6 × 6 quadrature is used. In all problems, the
generalized Heaviside function H (x) and the two-dimensional asymptotic crack-tip displacement
functions (see Section 3) are used to model the crack. A conjugate gradient sparse iterative solver
IML++ [28] is used with the convergence tolerance  set to 10−8. The elastic constants used in the
computations are: Young’s modulus E=105 psi and Poisson’s ratio =0:3. The implementation
of the domain form of the contour J -integral and the parameters used in the computations are
described in the following section.
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1560 N. SUKUMAR ET AL.
6.1. Computation of stress intensity factors
Domain integral methods [29{31] are used to evaluate stress intensity factors along the three-
dimensional crack front. For the mode I crack problems considered here, the stress intensity factor
at a point s on the crack front is given by
KI (s)=
r
J (s)E
1− 2
(18)
With the normal to the crack front (and in the crack plane) oriented along the xˆ1-axis of a local
co-ordinate system, the pointwise J -integral is given by
J (s)= lim
→0
Z
(s)
H1ˆn d (= 1ˆ; 2ˆ) (19)
where
H1ˆj =W1ˆj − ijui;1ˆ (i; j= 1ˆ; 2ˆ; 3ˆ) (20)
Domain integral representations of the crack-tip contour integral provide a convenient and accurate
method for evaluating stress intensity factors in two- or three-dimensional fracture. For linear
elastostatics, in the absence of body forces and material inhomogeneities, and assuming traction-
free crack surfaces, the volume form of the domain integral is given by [29]
J (s)= −
Z
V
(Hkjqk; j + Hkj; jqk) dV
Z
Lc
lknk ds
(21)
where V is a volume enclosing the crack front, nk(s) are components of the in-plane unit outward
normal at s, lk(s) are components of an arbitrary unit vector at s lying in the plane of the crack
and Lc is the perturbed segment (virtual extension) along the crack front (Figure 5). The vector
eld qk is dened in V as
qk =
8
>
<
>
:
lk on St
0 on S0
arbitrary otherwise
(22)
6.1.1. Implementation in X-FEM. The mode I stress intensity factors are computed by the domain
form of the contour J -integral. Consider a point s on the crack front where the stress intensity
factor is to be computed. In order to evaluate Equation (21), a virtual extension domain is required
around the point s on the crack front [20]. Typically, in fracture analysis with nite elements, the
mesh is constructed so that the virtual extension domain is the union of nite elements in the
vicinity of the point s. Since the crack is not modelled as part of the nite element mesh in
X-FEM, this is not a natural choice. In the two-dimensional implementation of X-FEM [24],
elements that were within a characteristic distance of the crack-tip were included in the virtual
extension domain, which does not readily generalize to an easy-to-implement algorithm in three
dimensions. An alternative approach is to use an independent grid of hexahedral cells around the
point s to dene the virtual extension domain. This approach is used in EFG for fracture problems,
and we adopt the same approach in the three-dimensional implementation of X-FEM.
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1561
Figure 5. Volume representation for J -integral
evaluation.
Figure 6. Local crack front co-ordinate axes for
virtual extension domain.
We describe the virtual extension domain for the elastostatic problems presented in Section 6. Let
the perturbed segment (Lc) along the crack front be parameterized by the curvilinear
co-ordinate s :R3→ [0; Lc]. Consider the evaluation of the domain integral at a point s on the
crack front. The dimensions of the virtual extension domain are: L1ˆ; L2ˆ and L3ˆ along the x1ˆ- ; x2ˆ-,
and x3ˆ-co-ordinate directions, respectively. The point s is located at the origin of the local orth-
ogonal x1ˆx2ˆx3ˆ-co-ordinate system (Figure 6). The x2ˆ-axis is contained in the plane of the crack
and is tangential to the crack front at s. The J -domain is sub-divided into 6× 2× 6 hexahedral
cells, and 6× 6× 6 quadrature rule is used in each integration cell. In three-dimensional fracture,
the denominator in Equation (21) is well approximated by the formula (L+
2ˆ
+ L−
2ˆ
)=2 [32], where
L+
2ˆ
and L−
2ˆ
are the element lengths on either side of s in the xˆ2-direction. On using the same
approximation, we obtain
Z
Lc
lknk ds ≈
L2ˆ
2
(23)
which is used in the X-FEM domain integral computations. We dene the vector function qk in
the local crack-front co-ordinate system. The function q1ˆ is chosen to be a trilinear function which
vanishes on the boundary of the virtual extension domain and is unity at the point s. In the local
crack-front co-ordinate system,
q1ˆ(x1ˆ; x2ˆ; x3ˆ) =

1−
2|x1ˆ|
L1ˆ

1−
2|x2ˆ|
L2ˆ

1−
2|x3ˆ|
L3ˆ

; q2ˆ = 0; q3ˆ = 0 (24)
with the aid of the above formulas for various terms that appear in Equation (21), we evaluate
the domain form of the J -integral.
6.2. Benchmark problems
For an innite domain, two well known and widely used benchmarks are the penny crack and the
elliptical crack under pure mode I loading conditions. The problem of an edge crack in a nite-
thickness plate under tension [33] is also used to study the accuracy of new numerical techniques
in three-dimensional SIF computations. We compare the numerical stress intensity factors to the
SIF solutions for the above problems.
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1562 N. SUKUMAR ET AL.
Figure 7. Mesh (surface) for the penny crack problem.
6.2.1. Planar penny crack in an innite domain. Let a be the radius of a penny crack with x3 the
co-ordinate axis normal to the plane of the crack. The closed-form solution for the stress intensity
factor along the crack front of a penny crack in an innite domain under uniaxial tension is given
by [34; 35]
KEI =2
0
33
r
a

(25)
Consider a penny crack of radius a=0:1 inside a bi-unit cube. The specimen is subjected to the
stress 033 = 1 in the x3-direction. Since the eects of the nite size of the model are minimal, we
use the exact solution given in Equation (25). Two dierent meshes are used: (a) Mesh 1 consists
of 24× 24× 24 hexahedral elements, and (b) Mesh 2 has 24× 24× 25 hexahedral elements. In
the former, the crack lies on element faces, whereas in the latter, the crack is located in the centre
of the elements (no intersection with any element faces). Both meshes have graded renement
towards the centre of the cube and in the plane of the crack (Figure 7). The nodes normal to the
crack surface are equi-spaced with length h.
In Figure 8, a contour plot of the crack plane (x3 = 0) nodal enrichment for the 24× 24× 24
mesh is shown. The domain of enrichment extends an element length (h) above and below the
crack plane: identical contour plots to that shown appear on x3 =± h. In Figure 8(a), the nodes
enriched on the crack plane by the generalized Heaviside function are indicated, and in Figure 8(b),
the nodes enriched by the crack-tip functions are shown. In both the contour plots, colour codes
are assigned a number between 0 and 4, where the number indicates the number of nodes on the
face of the element that are enriched. As one can observe, the mesh is fairly coarse in the vicinity
of the crack front. In Table I, the parameters for the system equations are indicated.
The dimensions of the virtual extension domain are: L1ˆ = 2a; L2ˆ = a=2; L3ˆ = 2a, and L1ˆ =L2ˆ =
a=2; L3ˆ = 2:5a for meshes 1 and 2, respectively. The SIF results have four-fold symmetry; hence
results for only 0◦6690◦ are presented. In Table II, the numerical SIF results are presented as
a function of , and in Figure 9, the normalized SIF is shown as a function of . The X-FEM
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1563
Figure 8. Nodal enrichment on the crack plane for the penny crack problem:
(a) generalized Heaviside function; and (b) crack-tip functions.
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1564 N. SUKUMAR ET AL.
Table I. System parameters for the penny crack problem.
Number of Non-zero Sparsity
Mesh unknowns entries (%)
24 24 24 48 948 4 212 194 0.17
24 24 25 42 678 3 531 390 0.19
Table II. Stress intensity factors for the penny crack problem.
 KI KI
(deg) KEI [24 24 24] [24 24 25]
0 0.3568 0.3508 0.3617
10 0.3568 0.3583 0.3503
20 0.3568 0.3592 0.3559
30 0.3568 0.3467 0.3595
40 0.3568 0.3535 0.3383
50 0.3568 0.3535 0.3383
60 0.3568 0.3467 0.3595
70 0.3568 0.3592 0.3559
80 0.3568 0.3583 0.3503
90 0.3568 0.3508 0.3617
Figure 9. Normalized stress intensity factors for
the penny crack problem.
Figure 10. Geometric denitions for an
elliptical crack.
results are observed to be in good agreement with the exact solution: the errors in the SIFs are
between 0.4{2.9 per cent for mesh 1 and 0.3{5.2 per cent for mesh 2.
6.2.2. Planar elliptical crack in an innite domain. We consider an elliptical crack with semi-
major axis a=0:1 and semi-minor axis b=0:05 inside a bi-unit cube. The body is subjected to unit
tractions 033 = 1 on x3 =± 1. The co-ordinate x3-axis is normal to the plane of the crack. Since the
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1565
Table III. Stress intensity factors for the elliptical
crack problem.

(deg) KEI KI
0 0.2314 0.2358
10 0.2365 0.2378
20 0.2495 0.2459
30 0.2662 0.2564
40 0.2830 0.2776
50 0.2983 0.2965
60 0.3107 0.3089
70 0.3198 0.3163
80 0.3253 0.3194
90 0.3273 0.3202
crack dimensions are small compared to the specimen, we use the innite domain solution as the
reference solution. The exact SIF solution for a planar elliptical crack in an innite domain is [35]
KEI =
033
√
b
E(k)

sin2 +
b2
a2
cos2 
1=4
(26)
where  is the elliptic angle (Figure 10), 033 is the far-eld applied stress in the x3-direction, and
E(k) which is the elliptic integral of the second kind is given by
E(k)=
Z =2
0
p
1− k2sin2  d; k2 =
a2 − b2
a2
(27)
The nite element mesh consists of 24× 24× 24 hexahedral elements. The number of unknowns
in the matrix system is 48 324, with 3 993 602 non-zero entries in the stiness matrix{matrix spar-
sity is 0.17 per cent. The dimensions of the virtual extension domain are: L1ˆ = 2b; L2ˆ = b; L3ˆ = 4b.
The SIF results have four-fold symmetry; hence results for only 0◦6690◦ are presented. For
the chosen values of a and b, the value of the elliptic integral E(k)= 1:211096 [36]. In Table III,
the SIF results are presented as a function of the elliptic angle , and in Figure 11, the normalized
SIF is plotted versus . The agreement between the exact solution and the numerical results is
good for the entire range of . The minimum (=10◦) and maximum (=30◦) errors in the
stress intensity factors are 0.6 and 3.7 per cent, respectively.
6.2.3. Single edge-crack tension specimen. A single edge-crack tension specimen (Figure 12)
subjected to unit tractions in the x3-direction is analysed. The dimensions of the specimen are:
h=w=1:75; a=w=1:0, and t=w=3:0. This benchmark problem was studied by Raju and Newman
[33] using singular nite elements, and by Li et al. [37], who used a symmetric weak-form
boundary integral equation method. Two dierent nite element meshes are considered in the X-
FEM model: (a) Mesh 1 consists of 20× 20× 20 hexahedral elements, with each element a cuboid
of dimensions h1 × h2 × h3 (h1 = 0:5, h2 = 0:75, h3 = 0:875); and (b) Mesh 2 has 40 × 40 × 40
hexahedral elements, with each element a cuboid of dimensions h1×h2×h3 (h1 = 0:25, h2 = 0:375,
h3 = 0:4375). Unit tractions 033 = 1 are imposed on the top and bottom surfaces. The solution
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1566 N. SUKUMAR ET AL.
Figure 11. Stress intensity factors for the elliptical
crack problem
Figure 12. Single edge-crack tension specimen.
Table IV. System parameters for the
single edge-crack problem.
Number of Non-zero Sparsity
Mesh unknowns entries (%)
20 20 20 30 612 2 941 268 0.30
40 40 40 213 522 17 955 158 0.04
parameters for the dierent meshes are indicated in Table IV. The dimensions of the virtual
extension domain are: L1ˆ = 6h1, L2ˆ = 2h2, and L3ˆ = 6h2.
The SIF results for the hexahedral meshes are presented in Table V, and in Figure 13, com-
parisons with the results obtained by Raju and Newman [33] and Li et al. [37] are shown. The
arc length co-ordinate s is measured along the crack front, from its centre (s=0) to the point
where it meets the free surface
(
s= t=2


. The stress state at the centre of the specimen is under
near plane strain conditions and the numerical SIF approaches the two-dimensional plane strain
result. Good agreement between the X-FEM results and the reference solutions are obtained in
the interior, away from the free surface. The SIF results due to Li et al. [37] are greater than the
other numerical results towards the centre of the crack front.
The stress state at the intersection of a crack with a free surface has been the subject of many
studies. The order of the singularity at the intersection of a crack and a free surface is r− (¡ 12 ),
where  is dependent on the Poisson’s ratio [38]. Bazant and Estenssoro [39] investigated the stress
singularity index by a special nite element method for mode I{III cracks. Tan and Fenner [40]
used boundary integral equations to study the nature of the singularity at the intersection between
a crack and a free surface and showed numerically that the stresses (strains) near the free surface
are less singular than the inverse square root form, and the eect becomes more pronounced with
increasing Poisson’s ratio. The trend towards a diminishing KI as the free surface is approached
is observed in the X-FEM as well as other reference solution results.
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1567
Table V. Stress intensity factors for the single edge-crack
under uniaxial tension problem.
KI
033
p
a
KI
033
p
a
s
t
[20 20 20] [40 40 40]
0.00 2.8470 2.7873
0.05 2.8481 2.7869
0.10 2.84703 2.7858
0.15 2.8449 2.7836
0.20 2.8412 2.7799
0.25 2.8350 2.7736
0.30 2.8244 2.7629
0.35 2.8056 2.7441
0.40 2.7702 2.7084
0.45 2.6796 2.6283
Figure 13. Comparison of stress intensity factors for the single edge-crack problem.
Table VI. System parameters for the penny crack in a
nite cube problem.
Number of Non-zero Sparsity
Mesh unknowns entries (%)
20 20 20 50 568 4 770 050 0.19
40 40 40 213 252 18 115 070 0.04
6.3. Planar penny crack in a nite cube
As an example of an embedded crack in a nite domain, we consider the problem of a penny
crack of radius a located at the centre of a bi-unit cube. We choose a=0:5, and compare the
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1568 N. SUKUMAR ET AL.
Table VII. Stress intensity factors for the penny crack in a
nite cube problem.
 KI KI
(deg) [24 24 24] [40 40 40]
0 0.8641 0.8685
10 0.8706 0.8792
20 0.8848 0.8933
30 0.8702 0.8863
40 0.8754 0.8920
50 0.8754 0.8920
60 0.8704 0.8864
70 0.8853 0.8935
80 0.8709 0.8799
90 0.8644 0.8689
numerical results to those obtained by Li et al. [37] who used a symmetric weak-form boundary
integral equation formulation. As a reference solution, a rened mesh solution was adopted in
Reference [37] for the stress intensity factor at =0◦:
K refI =
2:213033

√
a; (=0◦) (28)
where 033 is the applied stress in the x3-direction. Two dierent nite element meshes are con-
sidered: (a) Mesh 1 consists of 24 × 24 × 24 hexahedral elements, with each element a cube of
edge length h= 112 ; and (b) Mesh 2 has 40 × 40 × 40 hexahedral elements, with each element
a cube of edge length h=0:05. Unit tractions 033 = 1 are imposed on the surfaces x3 = ± 1. In
Table VI, the solution parameters are indicated. The dimensions of the virtual extension domain
are: L1ˆ = 2h, L2ˆ = 2h, L3ˆ = 2h. In Table VII, the SIF results are presented as a function of . With
renement, improved SIFs are obtained for all , which indicates the correct trend towards con-
vergence. Comparing the numerical solution at =0◦ to the reference solution given in Equation
(28), we observe that the dierence between the numerical result and the reference solution for
meshes 1 and 2 are about 2.0 per cent and 1.6 per cent, respectively.
7. CONCLUSIONS
The formulation and implementation of the extended nite element method (X-FEM) for three-
dimensional crack modelling was described. In X-FEM, the nite element space is enriched by
adding special functions to the approximation using the notion of partition of unity. For three-
dimensional crack modelling, a discontinuous function was used to model the interior of the crack
surface, and functions from the two-dimensional asymptotic crack-tip displacement elds were
used for the crack front enrichment. These enrichment functions were added to the nite element
approximation within the context of a displacement-based Galerkin formulation.
Issues pertaining to accurate and robust geometric computations in X-FEM were addressed, and
some of the algorithms were described. The performance of the extended nite element method for
three-dimensional static cracks was studied. Benchmark mode I problems of penny and elliptical
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

FINITE ELEMENT METHOD FOR CRACK MODELLING 1569
cracks in an innite domain were solved. The numerical stress intensity factors (SIFs) were found
to be in good agreement with the exact solution for these problems.
The SIFs obtained using X-FEM for a single edge-crack in a nite-thickness specimen
approached the plane strain results towards the centre of the specimen and showed a signi-
cant decrease as the free surface was approached, which corresponded well to other reference
solutions. The SIF results at =0◦ for an embedded penny crack in a nite cube also showed
good agreement with the solution due to Li et al. [37].
This study demonstrates the accuracy of the extended nite element method in three-dimensional
SIF computations. The method provides a robust and versatile numerical tool to solve crack prob-
lems in complex structural components without the need to explicitly align the mesh with the
crack. By eliminating the need to include the crack surfaces in the model, mesh generation is
greatly simplied. Furthermore, it facilitates the modelling of crack growth for fatigue studies,
since a single mesh suces.
ACKNOWLEDGEMENTS
The authors are grateful for the research support of the National Science Foundation through contract CMS-
9732319, the Multi-University Center for Integrated Diagnostics, Georgia Institute of Technology ONR Grant
N0014-95-1-0539, the Federal Aviation Administration contract DTFA03-98-F-IA025 and the Oce of Naval
Research. The authors also thank John Dolbow, Jean-Francois Remacle, and Christophe Geuzaine for many
helpful discussions.
REFERENCES
1. Henshell RD, Shaw KG. Crack tip nite elements are unnecessary. International Journal for Numerical Methods in
Engineering 1975; 9:495{507.
2. Barsoum RS. Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements. International Journal
for Numerical Methods in Engineering 1977; 11:85{98.
3. Akin JE. The Generation of elements with singularities. International Journal for Numerical Methods in Engineering
1976; 10:1249{1259.
4. Benzley SE. Representation of singularities with isoparametric nite elements. International Journal for Numerical
Methods in Engineering 1974; 8:537{545.
5. Giord Jr LN, Hilton PD. Stress intensity factors by enriched nite elements. Engineering Fracture Mechanics 1978;
10:485{496.
6. Gerstle WH, Martha L, Ingraea AR. Finite and boundary element modeling of crack propagation in two- and three
dimensions. Engineering with Computers 1987; 2:167{183.
7. Gerstle WH, Ingraea AR, Perucchio R. Three-Dimensional fatigue crack propagation analysis using the boundary
element method. International Journal of Fatigue 1988; 10(3):187{192.
8. Martha LF, Wawrzynek PA, Ingraea AR. Arbitrary crack representation using solid modeling. Engineering with
Computers 1993; 9:63{82.
9. Vijayakumar K, Atluri SN. An embedded elliptical
aw in an innite solid, subject to arbitrary crack-face tractions.
Journal of Applied Mechanics 1981; 48:88{96.
10. Nishioka T, Atluri SN. Analytical solution for embedded cracks, and nite element alternating method for elliptical
surface cracks, subjected to arbitrary loading. Engineering Fracture Mechanics 1983; 17:247{268.
11. Cruse TA. Boundary Element Analysis in Computational Fracture. Kluwer Academic Publishers: Dordrecht, The
Netherlands, 1988.
12. Murakami Y, Nemat-Nasser S. Interacting dissimilar semi-elliptical surface
aws under tension and bending.
Engineering Fracture Mechanics 1982; 16(3):373{386.
13. Xu Y, Moran B, Belytschko T. Self-similar crack expansion method for three-dimensional crack analysis. ASME,
Journal of Applied Mechanics 1997; 64:729{737.
14. Hills DA, Kelly PA, Dai DN, Korsunsky AM. Solution of Crack Problems: The Distributed Dislocation Technique.
Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996.
15. Dai DN, Hills DA, Nowell D. Modelling of growth of three-dimensional cracks by a continuous distribution of
dislocation loops. Computational Mechanics 1997; 19:538{544.
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570

1570 N. SUKUMAR ET AL.
16. Belytschko T, Lu YY, Gu, L. Element-free Galerkin methods. International Journal for Numerical Methods in
Engineering 1994; 37:229{256.
17. Fleming M, Chu YA, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip elds. International
Journal for Numerical Methods in Engineering 1997; 40:1483{1504.
18. Sukumar N, Moran B, Black T, Belytschko T. An element-free Galerkin method for three-dimensional fracture
mechanics. Computational Mechanics 1997; 20:170{175.
19. Krysl P, Belytschko T. The element-free Galerkin method for dynamic propagation of arbitrary 3-d cracks. International
Journal for Numerical Methods in Engineering 1999; 44(6):767{800.
20. Melenk JM, Babuska I. The partition of unity nite element method: basic theory and applications. Computer Methods
in Applied Mechanics and Engineering 1996; 139:289{314.
21. Strouboulis T, Babuska I, Copps K. The design and analysis of the generalized nite element method. Computer
Methods in Applied Mechanics and Engineering 2000; 181:43{71.
22. Strouboulis T, Copps K, Babuska I. The generalized nite element method: an example of its implementation and
illustration of its performance. International Journal for Numerical Methods in Engineering 2000; 47(8):1401{1417.
23. Belytschko T., Black T. Elastic crack growth in nite elements with minimal remeshing. International Journal for
Numerical Methods in Engineering 1999; 45(5):601{620.
24. Moes N, Dolbow J, Belytschko T. A nite element method for crack growth without remeshing. International Journal
for Numerical Methods in Engineering 1999; 46(1):131{150.
25. Dolbow J. An extended nite element method with discontinuous enrichment for applied mechanics. Ph.D. Thesis,
Theoretical and Applied Mechanics, Northwestern University, Evanston, IL, USA, 1999.
26. Shewchuk JR. Adaptive precision
oating-point arithmetic and fast robust geometric predicates. Discrete and
Computational Geometry 1997; 18:305{363.
27. Remacle JF, Geuzaine C. Gmsh nite element grid generator. Available at http://scorec rpi edu/ ˜remacle/
Gmsh Eng.html, 1998.
28. Dongarra J, Lumsdaine A, Pozo R, Remington K. IML version 1.2: Iterative methods library reference guide. Technical
Report, 1998. Available at http://www.math.nist.gov/iml++.
29. Moran B, Shih CF. Crack tip and associated domain integrals from momentum and energy balance. Engineering
Fracture Mechanics 1987; 27(6):615{641.
30. Moran B, Shih C.F. A general treatment of crack tip contour integrals. International Journal of Fracture 1987;
35:295{310.
31. Nikishkov GP, Atluri SN. Calculation of fracture mechanics parameters for an arbitrary three-dimensional crack by
the ‘equivalent domain integral method’. International Journal for Numerical Methods in Engineering 1987;
24:1801{1821.
32. Gosz M, Dolbow J, Moran B. Domain integral formulation for stress intensity factor computation along curved three-
dimentional interface cracks. International Journal of Solids and Structures 1998; 35:1763{1783.
33. Raju IS, Newmann Jr JC. Three-dimensional nite-element analysis of nite-thickness fracture specimens. Technical
Report NASA TND-8414, NASA Langley Research Center, Hampton, VA 23665, 1977.
34. Green AE, Sneddon IN. The distribution of stress in the neighborhood of a
at elliptical crack in an elastic solid.
Proceedings of the Cambridge Philosophical Society 1950; 47:159{164.
35. Irwin GR. The crack extension force for a part-through crack in a plate. ASME, Journal of Applied Mechanics 1962;
29:651{654.
36. Byrd RF, Friedman MD. Handbook of Elliptic Integrals for Engineers & Scientists. Springer: Berlin; 1971.
37. Li S, Mear ME, Xiao L. Symmetric weak-form integral equation method for three-dimensional fracture analysis.
Computer Methods in Applied Mechanics and Engineering 1998; 151:435{459.
38. Benthem JP. The state of stress at the vertex of a quarter-innite crack in a half space. International Journal of Solids
and Structures 1977; 13:479{492.
39. Bazant ZP, Estenssoro, L. Surface singularity and crack propagation. International Journal of Solids and Structures
1979; 15:405{426.
40. Tan CL, Fenner RT. Elastic fracture mechanics analysis by the boundary integral equation method. Proceedings of the
Royal Society of London 1979; (369);A:243{260.
Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2000; 48:1549{1570