INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2003; 56:2015–2037 (DOI: 10.1002/nme.653)
Brittle fracture in polycrystalline microstructures with the
extended nite element method
N. Sukumar
1; 2; ∗; †
, D. J. Srolovitz
2; 3
, T. J. Baker
3
and J.-H. Prevost
2; 4
1
Department of Civil and Environmental Engineering; University of California; Davis; CA 95616; U.S.A.
2
Princeton Materials Institute; Bowen Hall; Princeton University; NJ 08544; U.S.A.
3
Department of Mechanical and Aerospace Engineering; Princeton University; NJ 08544; U.S.A.
4
Department of Civil and Environmental Engineering; Princeton University; NJ 08544; U.S.A.
SUMMARY
A two-dimensional numerical model of microstructural eects in brittle fracture is presented, with an
aim towards the understanding of toughening mechanisms in polycrystalline materials such as ceram-
ics. Quasi-static crack propagation is modelled using the extended nite element method (X-FEM) and
microstructures are simulated within the framework of the Potts model for grain growth. In the X-FEM,
a discontinuous 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.
Hence, crack propagation can be simulated without any user-intervention or the need to remesh as the
crack advances. The microstructural calculations are carried out on a regular lattice using a kinetic
Monte Carlo algorithm for grain growth. We present a novel constrained Delaunay triangulation algo-
rithm with grain boundary smoothing to create a nite element mesh of the microstructure. The fracture
properties of the microstructure are characterized by assuming that the critical fracture energy of the
grain boundary (G
gb
c
) is dierent from that of the grain interior (G
i
c
). Numerical crack propagation sim-
ulations for varying toughness ratios G
gb
c
=G
i
c
are presented, to study the transition from the intergranular
to the transgranular mode of crack growth. This study has demonstrated the capability of modelling
crack propagation through a material microstructure within a nite element framework, which opens-up
exciting possibilities for the fracture analysis of functionally graded material systems. Copyright
? 2003
John Wiley & Sons, Ltd.
KEY WORDS: polycrystalline microstructure; grain boundary; strong discontinuity; two-dimensional
crack propagation; Potts grain growth model; Delaunay triangulation; meshing micro-
structures; partition of unity; extended nite element method
∗
Correspondence to: N. Sukumar, Department of Civil and Environmental Engineering, University of California,
Davis, CA 95616, U.S.A.
†
E-mail: nsukumar@ucdavis.edu
Contract=grant sponsor: Idaho National Engineering Laboratory; contract=grant number: KOO-182412
Contract=grant sponsor: National Science Foundation; contract=grant number: NSF-9988788
Received 31 August 2001
Revised 31 January 2002
Copyright
?
2003 John Wiley & Sons, Ltd. Accepted 10 June 2002
2016 N. SUKUMAR ET AL.
1. INTRODUCTION
Understanding deformation and failure mechanisms in brittle polycrystalline materials such
as ceramics is critical for improvements in the development and application of advanced
structural materials. The material microstructure plays a pivotal role in dictating the modes of
fracture and failure, and the macroscopic response of real materials. The grain morphology,
elastic modulus, and the toughness of the individual microstructural constituents and interfaces
are key parameters that control the failure mechanisms in polycrystalline materials. Concepts
such as grain boundary design and control and grain boundary engineering to improve the
fracture resistance of polycrystalline materials are well-recognized [1–3]. Quasi-static crack
propagation through a material microstructure depends on the mechanical state in the vicinity
of the crack-tip, and hence local dierences in toughness (grain interior versus grain bound-
aries) signicantly inuence the crack path. In light of the above, it is clear that any numerical
fracture model that is able to model crack propagation by incorporating these microstructural
features has the potential to describe toughening mechanisms in polycrystals and provide a
framework for microstructural design.
Brittle fracture simulations in disordered (heterogeneous) materials using spring networks
has received wide attention since the late 1980s [4–6]. A detailed study on spring-networks and
nite element methods for crack propagation simulations was conducted by Jagota and Ben-
nison [7, 8], who pointed out the dilemma associated with spring networks—regular meshes
can model uniform strain but show strong anisotropy in crack propagation whereas random
networks cannot in general represent uniform strain. Schlangen and Garboczi [9] considered
the appropriate selection of cross-sectional area and moment of inertia for beam elements in
a random lattice to simulate a homogeneous medium with showed no mesh dependency. An
elegant solution and partial resolution to the above shortcomings was provided by Bolander
and Saito [10] who used rigid-body spring networks to model brittle fracture in homogeneous
isotropic materials such as cement and concrete. They proposed a model for random networks
based on the Voronoi tessellation which was able to produce homogeneous deformation on
uniform straining and showed little bias towards crack propagation directions.
Lattice spring models have been used to study brittle fracture and damage in polycrystalline
materials. Yang and co-workers [11] used the Potts grain growth model [12–16] to generate
a polycrystalline microstructure which was mapped onto a triangular lattice. The mechanics
of this structure was represented by a spring network on the lattice, where a spring fails
if the stored elastic energy in the spring exceeded a critical value. In Reference [11], the
transition from intergranular (growth along the grain boundary) to transgranular (growth in the
grain interior) fracture with increasing grain boundary toughness was observed. In Reference
[17, 18], the eects of thermal-mismatch on microcracking was studied. Holm [19] considered
surface formation energies to study intergranular fracture in polycrystals; the inuence of
low/high-angle grain boundaries and grain boundary microcracking on the fracture path was
investigated. Kim and co-workers [20] analysed two-dimensional crack propagation through
a polycrystal as a function of the grain boundary toughness, focusing on the competition
between intergranular and transgranular mode of crack propagation.
An alternative approach to modelling fracture phenomena for arbitrary microstructure is the
use of cohesive surfaces within a nite element formulation. Zhai and Zhou [21] proposed
a micromechanical model in which the cohesive surface formulation of Xu and Needleman
[22] is used to study failure modes in composite microstructures, whereas in Reference [23],
Copyright ? 2003 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2003; 56:2015–2037
BRITTLE FRACTURE IN POLYCRYSTALLINE MICROSTRUCTURES WITH THE X-FEM 2017
intergranular microcrack initiation and evolution is assessed in ceramic microstructures under
dynamic loading conditions. The cohesive surface formulation is a phenomenological frame-
work in which the fracture characteristics of the material are embedded in a cohesive surface
traction–displacement relation. Using this approach, an inherent length scale is introduced into
the model, and in addition no fracture criterion (K-dominant eld) is required; crack growth
and the crack path are outcomes of the analysis. However, the constraint that cohesive sur-
faces can only lie along element edges does tend to render crack paths that are biased by the
mesh spacing and orientation [24].
Early numerical models for treating discontinuities in nite elements can be traced to the
work of Ortiz and co-workers [25] and Belytschko and co-workers [26]. They modelled shear
bands as ‘weak’ (strain) discontinuities that could pass through nite elements using a multi-
eld variational principle. Dvorkin and co-workers [27] considered ‘strong’ (displacement)
discontinuities by modifying the principle of virtual work statement. A unied framework
for analysing strong discontinuities by taking into account the softening constitutive law and
the interface traction–displacement relation was put forth by Simo and co-workers [28, 29].
Applications and extensions of this approach have been proposed by many researchers todate;
for a few related studies, see References [30–33]. In the strong discontinuity approach, the
displacement consists of regular and enhanced components, where the enhanced component
yields a jump across the discontinuity surface. An assumed enhanced strain variational for-
mulation is used, and the enriched degrees of freedom are statically condensed on an element
level to obtain the tangent stiness matrix for the element. It is to be noted that, in this
approach, the discontinuity surface can only terminate on element boundaries. A comprehen-
sive review and comparison of various embedded discontinuity approaches is provided by
Jirasek [34].
The use of nite elements to study discrete crack propagation through a polycrystalline
microstructure has not been fully explored todate. This stems from the fact that in order
to accurately capture the microstructural features and the crack path, rened meshes with
continuous adaptive remeshing techniques are required for discrete crack growth simulations.
The computational eort and complexity involved is signicant even in two-dimensions, and
hence the above approach has not attracted wide attention.
A signicant improvement in discrete crack modelling has been realized with the devel-
opment of the extended nite element method (X-FEM) [35–37]. In this approach, the do-
main is modelled by nite elements with no explicit meshing of the crack surfaces. The
location of the crack discontinuity can be arbitrary with respect to the underlying nite
element mesh, and quasi-static or fatigue crack propagation simulations can be performed
without the need to remesh as the crack advances. In the X-FEM, a discontinuous func-
tion (generalized Heaviside step function) and the two-dimensional asymptotic crack-tip dis-
placement elds are added to a standard displacement-based nite element approximation to
account for the presence of the crack using the notion of partition of unity [38, 39]. A par-
ticularly appealing feature is that the nite element framework and its properties (sparsity
and symmetry) are retained, and a single-eld (displacement) variational principle is used
to obtain the discrete equations. The classical nite element degrees of freedom as well as
the enriched degrees of freedom are found simultaneously by solving the discrete system.
This technique provides an accurate and robust numerical method that removes the need
to mesh the crack geometry in both two-dimensional [35, 36] and three-dimensional crack
modelling [40].
Copyright ? 2003 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2003; 56:2015–2037
2018 N. SUKUMAR ET AL.
In this paper, we present a numerical technique to carry out crack propagation simulations
through arbitrary material microstructures. Quasi-static crack propagation is modelled using the
X-FEM (Section 2) and microstructures are simulated within the framework of the Potts model
for grain growth (Section 3). The microstructural calculations are carried out on a regular
lattice, and a novel constrained Delaunay triangulation algorithm is developed to construct
the initial nite element mesh of the microstructure. A detailed description of the Delaunay
algorithm is presented in Section 4. A standard displacement-based Galerkin method is used
to obtain the discrete equations for linear elastostatics (Section 5). The fracture properties of
the microstructure are characterized by assuming that the critical fracture energy of the grain
boundary (G
gb
c
) is dierent from that of the grain interior (G
i
c
). The crack growth criterion
in the grain interior is governed by the maximum circumferential stress criterion, whereas
along the grain boundary, the growth direction is determined by selecting the one that has the
maximum value for G=G
k
c
(k is either gb or i). In Section 6, numerical simulations of crack
propagation through a microstructure are performed for dierent toughness ratios G
gb
c
=G
i
c
,to
study the transition from the intergranular to the transgranular mode of crack growth. The
main results and conclusions from this study are discussed in Section 7.
2. EXTENDED FINITE ELEMENT METHOD
The X-FEM [35, 36] is a numerical method to model internal (or external) boundaries such
as holes, inclusions, or multiple cracks, without the need for the mesh to conform to these
boundaries. The X-FEM is based on a standard Galerkin procedure, and uses the concept of
partition of unity [38, 39] to accommodate the internal boundaries in the discrete model. The
partition of unity method generalized nite element approximations by presenting a means to
embed local solutions of boundary-value problems into the nite element approximation. This
idea was exploited in References [41, 42] for problems with cracks and holes—the numerical
technique was referred to as the generalized nite element method (GFEM).
Partition of unity enrichment for discontinuities and near-tip crack elds was introduced
by Belytschko and Black [43]. The displacement enrichment functions for crack problems are
functions that span the asymptotic near-tip displacement eld. A signicant improvement in
discrete crack growth modelling without the need for any remeshing strategy was conceived
in Reference [35], with further extensions of the technique for modelling holes and branched
cracks presented in Reference [36]. The generalized Heaviside step function was proposed as
a means to model the crack away from the crack-tip, with simple rules for the introduction
of the discontinuous and crack-tip enrichments. This advance has clearly provided an accurate
and robust computational tool for modelling discontinuities independent of the mesh geometry
in two-dimensions [35, 36] and three-dimensions [40]. In addition, recent studies have explored
the use of fast marching and level sets for evolving crack discontinuities in three-dimensions
within the X-FEM framework: growth of multiple planar cracks are handled using the fast
marching method [44, 45], whereas non-planar crack growth is carried out using level sets
[46, 47]. The ideas and developments in the X-FEM have had an impact in other related areas.
For example, Wells and Sluys [48] proposed a cohesive crack model using nite elements
that adopts the notion of partition of unity and the use of the Heaviside step function as an
enrichment function to model the displacement discontinuity.
Copyright ? 2003 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2003; 56:2015–2037
BRITTLE FRACTURE IN POLYCRYSTALLINE MICROSTRUCTURES WITH THE X-FEM 2019
2.1. Displacement approximation
The enrichment of the nite element approximation is described as follows. Consider a point
x of
R
d
(d= 1–3) that lies inside a nite element e. Denote the nodal set N = {n
1
;n
2
;:::;n
m
},
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):
R
d
→ R
d
assumes
the form:
u
h
(x)=
I
n
I∈N
I
(x)u
I
classical
+
J
n
J ∈N
g
J
(x) (x)a
J
enriched
(u
I
; a
J
∈R
d
) (1)
where the nodal set N
g
is dened as
N
g
= {n
J
: n
J
∈N;!
J
∩
g
= ∅} (2)
In the above equation, !
J
= supp(n
J
) is the support of the nodal shape function
J
(x), which
consists of the union of all elements with n
J
as one of its vertices, or in other words the
union of elements in which
J
(x) is non-zero. In addition,
g
is the domain associated with
a geometric entity such as crack-tip [35], crack surface in three-dimensions [40], or material
interface [49]. In general, the choice of the enrichment function (x) that appears in Equation
(1) is dependent on the geometric entity.
2.2. Two-dimensional crack modelling
The crack is assumed to consist of linear one-dimensional segments. 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 supports
intersect the crack tip by the two-dimensional asymptotic crack-tip elds. A short description
of the selection of nodes for enrichment as well as the computation of enrichment functions
follows; for further details, the interested reader can refer to Reference [35]. In addition to the
above, partitioning algorithms are also implemented if the crack intersects the nite elements;
in this study, the algorithm described in Reference [40] is used with some minor additions
for improvement.
2.2.1. Enrichment functions. Consider a single crack in two-dimensions, and let
c
be the
crack surface and
c
the crack tip. The interior of a 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 segment that contains
x
∗
. The H (x) function is then given by
H (x)=
1 if (x − x
∗
) · n¿0
−1 otherwise
(3)
To model the crack tip and also to improve the representation of crack-tip elds in fracture
computations, crack-tip enrichment functions are used in the element which contains the crack
Copyright ? 2003 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2003; 56:2015–2037
