A Geometrically Nonlinear Phase Field Theory of Brittle

Fracture

by JD Clayton and J Knap

ARL-RP-0511 October 2014

A reprint from the Int J Fract, 2014;189(2):139–148. 

Published online in the Springer Science+Business Media Dordrecht [accessed 2014 Aug 19].

http://link.springer.com/article/10.1007/s10704-014-9965-1.

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the Army position unless 

so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the 

use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
Army Research Laboratory

Aberdeen Proving Ground, MD 21005-5069

ARL-RP-0511 October 2014

A Geometrically Nonlinear Phase Field Theory of Brittle

Fracture

JD Clayton and J Knap

Weapons and Materials Research Directorate, ARL

A reprint from the Int J Fract, 2014;189(2):139–148. 

Published online in the Springer Science+Business Media Dordrecht [accessed 2014 Aug 19].

http://link.springer.com/article/10.1007/s10704-014-9965-1.

Approved for public release; distribution is unlimited.
ii

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

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 information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the 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)

October 2014

2. REPORT TYPE

Reprint

3. DATES COVERED (From - To)

October 2013–October 2014

4. TITLE AND SUBTITLE

A Geometrically Nonlinear Phase Field Theory of Brittle Fracture

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

JD Clayton and J Knap

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

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

US Army Research Laboratory

ATTN: RDRL-WMP-C

Aberdeen Proving Ground, MD 21005-5069

8. PERFORMING ORGANIZATION

 REPORT NUMBER

ARL-RP-0511

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

Approved for public release; distribution is unlimited.

13. SUPPLEMENTARY NOTES

A reprint from the Int J Fract, 2014;189(2):139–148. Published online in the Springer Science+Business Media Dordrecht 

[accessed 2014 Aug 19]. http://link.springer.com/article/10.1007/s10704-014-9965-1.

14. ABSTRACT

Phase field theory is developed for solids undergoing potentially large deformation and fracture. The elastic potential depends 

on a finite measure of elastic strain. Surface energy associated with fracture can be anisotropic, enabling description of preferred

cleavage planes in single crystals, or isotropic, applicable to amorphous solids such as glass. Incremental solution of the 

Euler–Lagrange equations corresponds to local minimization of an energy functional for the solid, enabling prediction of 

equilibrium crack morphologies. Predictions are in close agreement with analytical solutions for pure mode I or pure mode II 

loading, including the driving force for a crack to extend from a pre-existing plane onto a misoriented cleavage plane. In an 

isotropic matrix, the tendency for a crack to penetrate or deflect around an inclusion is shown to depend moderately on the ratio 

of elastic stiffness in matrix and inclusion and strongly on their ratio of surface energy. Cracks are attracted to (shielded by)

inclusions softer (stiffer) than the surrounding matrix. The theory and results apparently report the first fully three-dimensional 

implementation of phase field theory of fracture accounting for simultaneous geometric nonlinearity, nonlinear elasticity, and 

surface energy anisotropy.

15. SUBJECT TERMS

brittle fracture, phase field, nonlinear elasticity, elastic inclusion, ceramics

 16. SECURITY CLASSIFICATION OF: 

17. LIMITATION

OF ABSTRACT

UU

18. NUMBER

OF PAGES

16

19a. NAME OF RESPONSIBLE PERSON

JD Clayton

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

410-278-6146

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
Int J Fract (2014) 189:139–148

DOI 10.1007/s10704-014-9965-1

OR IG INAL PAPER

A geometrically nonlinear phase field theory of brittle

fracture

J. D. Clayton · J. Knap

Received: 8 May 2014 / Accepted: 29 July 2014 / Published online: 19 August 2014

© Springer Science+Business Media Dordrecht (outside the USA) 2014

Abstract Phase field theory is developed for solids

undergoing potentially large deformation and fracture.

The elastic potential depends on a finite measure of

elastic strain. Surface energy associated with fracture

can be anisotropic, enabling description of preferred

cleavage planes in single crystals, or isotropic, applicable to amorphous solids such as glass. Incremental

solution of the Euler–Lagrange equations corresponds

to local minimization of an energy functional for the

solid, enabling prediction of equilibrium crack morphologies. Predictions are in close agreement with analytical solutions for pure mode I or pure mode II loading, including the driving force for a crack to extend

from a pre-existing plane onto a misoriented cleavage

plane. In an isotropic matrix, the tendency for a crack

to penetrate or deflect around an inclusion is shown

to depend moderately on the ratio of elastic stiffness

in matrix and inclusion and strongly on their ratio of

surface energy. Cracks are attracted to (shielded by)

inclusions softer (stiffer) than the surrounding matrix.

The theory and results apparently report the first fully

three-dimensional implementation of phase field theory of fracture accounting for simultaneous geometric

J. D. Clayton (B)

Impact Physics, US ARL,

Aberdeen, MD 21005-5066, USA

e-mail: jclayton@arl.army.mil

J. Knap

Computational Science, US ARL,

Aberdeen, MD 21005-5066, USA

nonlinearity, nonlinear elasticity, and surface energy

anisotropy.

Keywords Brittle fracture · Phase field ·

Nonlinear elasticity · Elastic inclusion

1 Introduction

Brittle solids—which include ceramics, typical rocks

and minerals, glass, and some metals—demonstrate a

tendency to fracture rather than deform plastically (e.g.,

by dislocation glide) when subjected to stresses exceeding their elastic limit. Brittle solids often demonstrate

strong directional bonding at the atomic scale and a

relatively large ratio of shear to bulk modulus (Gilman

2003), though exceptions are possible. In crystalline

solids, fractures may be transgranular (i.e., cleavage

Lawn 1968; Schultz et al. 1994) or intergranular (i.e., at

grain boundaries), while in amorphous or glassy solids,

fractures are independent of intrinsic microstructure

but highly dependent on pre-existing flaws, especially

surface flaws (Wilshaw 1971). In ceramics and glass,

depending on loading conditions, fractures may consist

of dominant crack(s) or distributed micro-cracks (Lawn

et al. 1994). Performance of ceramics in structural

applications is most often inhibited by their low fracture toughness. Toughness can often be improved by

modifying the microstructure to promote crack bridging (i.e., deflection or pull-out at second-phase particles, grain boundaries, or other heterogeneities), thus

123
140 J. D. Clayton, J. Knap

shielding stresses experienced by a propagating crack

and causing it to meander or arrest. Efforts to improve

ductility of ceramics and ceramic nanocomposites in

this manner include controlled introduction of weak

interfaces, coarsening and elongating the grain structure, and embedding the polycrystal with strong inclusions (Faber and Evans 1983; Lawn et al. 1994; Ohji et

al. 1998).

Resistance to fracture is often the controlling factor

limiting applications in which a brittle material can be

used; grain size and morphology can be engineered

to optimize effective ductility, impact resistance, or

fatigue life for a given application (Lawn et al. 1994;

Wiederhorn 1984). Selection or design of a material

for such an application requires understanding afforded

by a predictive model of fracture. Computer simulations enable descriptions of fracture in brittle solids

under complex loading conditions and for nonlinear

and anisotropic material behaviors, difficult, if not

impossible, to address using existing analytical methods (e.g., solutions available from linear elastic fracture

mechanics). Engineering finite element (FE) simulations often invoke continuum damage mechanics theories, wherein the tangent stiffness of a material element degrades as “damage” accumulates. Conventional

continuum damage mechanics theories (Clayton and

McDowell 2003, 2004; Sun and Khaleel 2004; Clayton 2006, 2008) require prescription of phenomenological kinetic equations specifying the rate of damage

accumulation; associated parameters in such equations

must be tuned to experiments similar to those simulated

numerically. Furthermore, solutions can depend on

mesh size. Simple models based on the notion of theoretical strength (Gilman 1960;Clayton 2009, 2010) can

provide insight into directionality of fracture resistance

but do not enable prediction of morphological crack

evolution for complex stress states. Cohesive models

of fracture (Xu and Needleman 1993; Espinosa and

Zavattieri 2003; Clayton 2005; Arias et al. 2007; Foulk

and Vogler 2010) enable simulation of discrete cracks

and branching; however, mesh construction should not

be arbitrary since crack paths are constrained to follow

element boundaries. Extended FE methods (Moes et

al. 1999) alleviate the latter problem by permitting discontinuities to traverse elements, but predictive physics

is compromised since the crack propagation direction

must be specified by a user-defined criterion. Similarly, numerical methods involving incremental crack

growth laws (Kim et al. 1996) require specification of

criteria for crack extension such as maximum energy

release (Nuismer 1975). Multi-scale techniques (Knap

and Ortiz 2003; Zhang et al. 2007) blending atomic

and continuum theory appear promising for describing problems wherein localized damage or defects are

contained within a relatively small volume of the entire

body, but are inhibited by difficulties with linking discrete and continuous regions and attendant numerical

complexity; purely atomic methods (Knap and Sieradzki 1999) are necessarily restricted to relatively small

system sizes and to short time scales for dynamic simulations.

Phase field theories of fracture (Jin et al. 2001; Eastgate et al. 2002; Del Piero et al. 2007; Hakim and

Karma 2009; Kuhn and Muller 2010; Abdollahi and

Arias 2012; Alber 2012; Borden et al. 2012; Hofacker

and Miehe 2012; Spatschek et al. 2011; Voyiadjis and

Mozaffari 2013) enable prediction of complex fracture

processes without introduction of spurious physics.

Apart from the usual elastic constants, essential material parameters entering such theories can usually be

directly related to surface energy and diffuse interfacial

width, the latter associated with regularization inherent in the phase field approach that renders solutions

mesh size independent. In essence, fractures evolve

naturally with applied loading as the body seeks minimum energy configurations. Mathematical analysis

(Del Piero et al. 2007) has demonstrated aspects of Γ -

type convergence of certain variational models towards

sharp interface models of Griffith type as interfacial

width is reduced.

This paper develops a phase field theory for fracture

of nonlinear elastic materials based on a variational

approach, with numerical implementation involving

incremental minimization of a suitable free energy

functional. The present novel model, whose governing

(Euler–Lagrange) equations are derived in an analogous way to an existing theory for deformation twinning (Clayton and Knap 2011a), is similar to a finite

strain variational model implemented in two dimensions (Del Piero et al. 2007), but is new in its incorporation of anisotropy and fully three-dimensional (3D)

FE implementation. Previous phase field approaches

have considered anisotropic fracture energy (Jin et al.

2001; Hakim and Karma 2009), albeit in the context

of linear elasticity and 2D numerical simulations. Furthermore, these and most other prior works (Eastgate et

al. 2002; Kuhn and Muller 2010; Abdollahi and Arias

2012; Alber 2012; Borden et al. 2012; Hofacker and

123
Phase field theory of brittle fracture 141

Miehe 2012; Voyiadjis and Mozaffari 2013) invoked a

dynamic approach, requiring specification of one or

more parameters controlling the time scale of fracture kinetics. In contrast, the present model invokes

incremental energy minimization, enabling the solution

of quasi-static fracture problems without the need for

kinetic parameter(s). This work represents the first fully

3D implementation of a phase field theory of fracture

accounting for geometric nonlinearity, nonlinear elasticity, and possible surface energy anisotropy. Results

of numerical simulations of mode I and mode II loading

validate the model. Simulation results on crack deflection around or penetration through a spherical inclusion then follow, providing new insight into crack–

inclusion interactions in nonlinear elastic solids not

available from 2D linear analytical solutions (Tamate

1968; Atkinson 1972; Erdogan et al. 1974; Evans

1974).

2 Phase field theory

2.1 Governing equations

Finite deformation theory (Clayton 2011; Clayton and

Knap 2011a, b, 2013) is used, where x = x(X) =

X + u(X) are spatial Cartesian coordinates of a material particle initially at X. Let ∇ denote the material

gradient operator. The deformation gradient, symmetric finite deformation tensor, and ratio of deformed to

initial volume are, respectively,

F = ∇x = 1 + ∇u, C = FTF, J = √

det C. (1)

Degrees of freedom are displacement u(X) and order

parameter η(X), where η = 0 for perfect material,

η = 1 for fully fractured material, and η ∈ (0, 1) within

diffuse interfaces between intact and failed domains,

for example. The total free energy functional for the

body Ω is

(u, η) =



Ω

[W(∇u, η) +

Γ

l η2

+κ : ∇η ⊗ ∇η]dΩ. (2)

Elastic strain energy per unit reference volume is W

and may degrade with increasing η, surface energy per

unit reference area is Γ , constant l is proportional to

the equilibrium crack width or thickness of the diffuse

interface, and κ is a symmetric second-order material

property tensor, in what follows of the form

κ = Γ l[1 + β(1 − m ⊗ m)], (3)

with m a unit normal vector to a preferred cleavage

plane in a crystal, for example. Parameter β penalizes

fractures on planes with orientations different than m.

For isotropic surface energy, e.g., as in a glass or a representation of a macroscopically homogeneous brittle

isotropic solid, β = 0. From functional (2), respective

Euler–Lagrange equations in Ω and boundary conditions on its external surface ∂Ω are derived using standard variational methods and the divergence theorem

(Clayton and Knap 2011a):

∇ · P = 0, ∂W/∂η + 2Γ η/l = 2∇ · κ∇η; (4)

t = Pn, h = 2κ : ∇η ⊗ n. (5)

First Piola–Kirchhoff stress is P = ∂W/∂ F =

∂W/∂∇u and is generally not symmetric, traction on a

surface element with reference normal n is t, and h is a

conjugate force to η on boundary ∂Ω (note that h = 0

along a free surface).

2.2 Elasticity

Similarly to previous phase field models for twinning

(Clayton and Knap 2011b, 2013), compressible neoHookean elasticity is considered with the following

strain energy density function:

W(C, η) = 1

2

[λ(ln J )

2 − μ(2 ln J − trC + 3)]. (6)

Lamé coefficients μ and λ depend on order parameter

η and possibly volume change measure J :

μ(η)=μ0[ζ +(1−ζ )(1−η)2], λ(η, J )=k− 2

3

μ.

(7)

Shear modulus μ degrades from its reference value μ0

to a minimum value ζμ0 in fully fractured domains,

where 0 < ζ  1. Bulk modulus k degrades in tension

but not in compression, in order to prohibit interpenetration:

k = k0{[ζ + (1 − ζ )(1 − η)2]J − 1+1 − J 

∗}.

(8)

The initial (not degraded) bulk modulus is k0 = λ0 +

2

3μ0. The following notation applies: x = 1∀x >

0, x = 0∀x ≤ 0, x∗ = 1∀x ≥ 0, and x∗ =

0∀x < 0. Poisson’s ratio in the undamaged solid is

ν = λ0/(2λ0+2μ0). The stress tensor is then, invoking

(6),

123
142 J. D. Clayton, J. Knap

P = μF + (λ ln J − μ)F−T

+

1

2

J (ln J )

2(∂k/∂ J )F−T. (9)

The rightmost term accounts for tension-compression

asymmetry in the degradation of the bulk stiffness.

Anisotropic elastic constants are not included here but

may be considered in the future by adapting methods

in Clayton and Knap (2011a); however, in many single

crystals, surface energy anisotropy is thought to dominate elastic anisotropy (Lawn 1968).

For purposes of comparison with known analytical

solutions, linear elasticity is also sometimes considered

in this work, whereby the strain energy density function

is

W(∇u, η)= 1

2

λ(∇ · u)

2 +

1

4

μtr{[(∇u + (∇u)

T)]2},

(10)

and elastic coefficients degrade in the same way as in

(7) and (8), with volume change in the linear approximation J → 1 +∇· u. The stress tensor [derivative

of (10) with respect to ∇u] in the linear theory is symmetric:

P = λ(∇ · u)1+2μ(∇u)sym+

1

2

(∇ · u)

2(∂k/∂∇u)1.

(11)

Comparison of nonlinear and linear results also enables

evaluation of possible importance of nonlinear aspects

of the model and thereby may suggest domains of valid

loading regimes for the linear theory.

The present model permits degradation of the bulk

modulus only when volume change is tensile and degradation of the shear modulus regardless of whether loading is tensile or compressive. This approach, which is

analogous to that implemented in Amor et al. (2009)

and one described in Spatschek et al. (2011), enables

the material to lose shear strength even when all three

principal strain components are negative. A different

theory implemented in Miehe et al. (2010), Hofacker

and Miehe (2012), Borden et al. (2012) in the small

strain setting decomposes the strain tensor into positive and negative parts following diagonalization. The

elastic strain energy dependent on the positive (tensile)

components degrades upon fracture, while the compressive energy is unaffected by damage (i.e., does

not vary with changing values of the order parameter). In contrast to the present approach, shear strength

loss does not occur when all three principal strains are

compressive (Borden et al. 2012). Suitability of one

approach over the other likely depends on the material

and loading regime, and could be evaluated by considering the failure behavior of a material subjected to

triaxial compression (Murrell 1965) or shock compression under lateral pre-stress (Clayton 2014). In brittle

solids such as ceramics and rocks, when confining or

lateral stress is not too large, shear fractures may occur

under compression, but extreme lateral stress/pressure

may suppress propagation of mode II cracks due to frictional resistance that tends to increase with increasing

confinement (Murrell 1965; Clayton 2010).

3 Numerical implementation

The theory is implemented numerically using Lagrangian finite elements, following general procedures outlined in Clayton and Knap (2011a). Incremental solutions to (4) are sought using conjugate gradient minimization of total energy functional  of (2), subject

to constraints associated with external boundary conditions on ∂Ω. Let δη(X) be a local change in order

parameter η at material point X induced by incremental

loads. The additional internal constraint δη(X) ≥ 0 if

η(X) ≥ 0.9 is used to render fracture irreversible (Del

Piero et al. 2007).

In the current implementation of the phase field theory, cracks represented by positive values of the order

parameter are predicted to follow paths dictated by

incremental total energy minimization, subject to the

irreversibility constraint described above. When this

constraint is active, the incremental energy minimization problem can be viewed as minimization of energy

of an alternative system with time dependent boundary

conditions associated with introduction of new surfaces

along which η ≥ 0.9 is prescribed; equilibrium equations (4) remain satisfied in solutions thus obtained for

this alternative system. If the irreversibility constraint

on δη is not enforced, then damage is reversible and

cracks will heal fully upon unloading, a feature noted

in other phase field models (Hakim and Karma 2009).

If damage is regarded as completely irreversible, then

dissipated energy can be defined as the contribution

to  of the final two terms in the integrand of (2)

corresponding to the surface energy of fracture associated with the order parameter and its spatial gradient.

This is consistent with definitions of dissipated energy

used elsewhere (Miehe et al. 2010; Borden et al. 2012;

Spatschek et al. 2011) when viscous effects are omit123
Phase field theory of brittle fracture 143

ted. Another definition of dissipated energy, valid for

partial reversibility of damage, could be constructed for

a loading–unloading cycle: net dissipated energy could

be measured by the difference between work expended

and then recovered by the boundary conditions, minus

any residual energy remaining in the system upon external unloading. Further review of capabilities of phase

field approaches to predict crack paths can be found in

Spatschek et al. (2011).

In order to make the results broadly accessible and

applicable to a number of classes of brittle materials,

numerical analyses are performed later in the context

of dimensionless parameters rather than property sets

peculiar to any specific solid. Energy density can be

normalized by μ0, such that W/μ0 depends only on

ν and other terms in the integrand of (2) depend on

dimensionless constants Γ ̄ = Γ /μ0l, β, and ζ . Taking ζ = 0.01 and l/R0 = constant, where R0 is a

fixed characteristic length associated with the problem

of interest, dimensionless solutions depend only the

choice of {ν, Γ,β  ̄ }, where β = 0 for isotropy. Properties representative of brittle solids are later chosen as

ν = 0.25 and Γ ̄ = 0.01, while β and m are explored

parametrically.

In order to evaluate accuracy and validity of the

model, simulations on a block of elastic material with

a pre-existing straight crack/notch are performed. The

block is of reference dimensions 50R0 ×50R0 ×25R0,

where R0 = 2l is the finite radius of the initial notch

tip, and the initial crack/notch is of length a = 25R0.

The plane of this pre-crack is X2 = 0 in local Cartesian

reference coordinates, and let (r, θ) be local reference

polar coordinates with origin at the notch tip. Along the

outer boundary, displacement boundary conditions are

imposed for either pure mode I or pure mode II loading

(Zhang et al. 2007; Rice 1968); e.g., for mode I:

u1 = Δ

4ar

2π

cos

θ

2



1 − 2ν + sin2 θ

2



, (12)

u2 = Δ

4ar

2π

sin

θ

2



2 − 2ν − cos2 θ

2



, (13)

where the load parameter is Δ = KI/2μ

√a. Similar

equations apply for mode II, where Δ = KII/2μ

√a

(Clayton and Knap 2013). For a sharp crack in an infinite medium, linear elastic fracture mechanics predicts

crack extension will occur when Δ ≥ KC/2μ

√a,

where KC = √2μG/(1 − ν), and G = 2Γ for an

isotropic material. Because of the finite domain and

crack radius, values of KI/II entering Δ are corrected

herein by computing the J-integral (Rice 1968) numerically. Linear elastic fracture concepts (e.g., J = G)

are used only to define boundary conditions and aid in

analysis and interpretation of results; the phase field

simulations do not directly incorporate linear elastic

fracture mechanics. To mimic plane strain u3 = 0 is

imposed along the boundary, and h = 0 is imposed on

all external surfaces. Finite element meshes consist of

≈ 4M tetrahedral elements of size small relative to l to

resolve order parameter gradients at fracture surfaces.

4 Numerical results

4.1 Cracking for isotropic and anisotropic surface

energies

Propagation of a pure mode I crack in a notched body

with isotropic surface energy (β = 0) is shown in

Fig. 1a, b for nonlinear elasticity for two different load

increments, and in Fig. 1c for linear elasticity. Initiation at the notch tip occurs at KI/KC  1.03 for both

nonlinear and linear elasticity, validating the model for

isotropic surface energy, and neo-Hookean elasticity

predicts a 10–15% longer crack than linear elasticity

at the larger applied displacement. Figure 1d shows the

effect of anisotropic surface energy with β = 100 for

a cleavage plane twisted at angle φ = π/4 relative to

the pre-existing crack plane; i.e., m = [cos φ, 0,sin φ],

where normalization of the loading for twist is (Wiederhorn 1984)

KC = 

2μG/(1 − ν)sec2 φ (twist). (14)

Figure 1e shows a similar result for a cleavage plane

tilted at angle φ = π/4, i.e., m = [cos φ,sin φ, 0],

where normalization for pure tilt misorientation is by

the factor (Wiederhorn 1984; Xu et al. 2003)

KC = 

2μG/(1 − ν)sec2(φ/2) (tilt). (15)

Shown in Fig. 2a, b are components of first Piola–

Kirchhoff stress tensor P of (9), normalized by the initial shear modulus. Stress concentrations are of large

magnitude (e.g., locally on the order of 10% of the

modulus) but remain bounded in part due to the finite

radius of the crack tip resolved numerically.

As validated in Fig. 3, the phase field model correctly predicts that twist misorientation resists crack

123
144 J. D. Clayton, J. Knap

Fig. 1 Phase field predictions for pure mode I or mode II loading,

with fractured material (η > 0.7) removed for visualization: a

nonlinear elastic, β = 0, KI/KC = 1.1 b nonlinear elastic,

β = 0, KI/KC = 2.4 c linear elastic, β = 0, KI/KC = 2.4

d nonlinear elastic, φ = π/4 (twist), β = 100, KI/KC = 1.1

e nonlinear elastic, φ = π/4 (tilt), β = 100, KI/KC = 1.5 (f)

linear elastic, β = 100, KII/KC = 1.5

propagation more than the same angle of tilt misorientation (Wiederhorn 1984). Figure 1f depicts a result for

pure mode II loading with β = 100 set to favor crack

propagation on m = [1, 0, 0]. In all cases shown in

Figs. 1 and 2, when imposed boundary displacements

are large enough that cracks become overdriven (e.g.,

K KC), they do not necessarily propagate completely through the specimen in an unstable manner

because of the boundary conditions that tend to repel

the crack as it approaches the lower edge and possible effects of nonzero ζ that serves to maintain some

residual stiffness in the fully damaged material.

In the examples involving misorientation (Fig. 1d,

e), cracks demonstrate non-uniform fronts suggestive

Fig. 2 Phase field predictions for pure mode I loading, nonlinear

elastic, β = 0, KI/KC = 2.4, with fractured material (η > 0.7)

removed for visualization: a normal stress component P11/μ0 =

Px X /μ0 b normal stress component P22/μ0 = PyY /μ0 [referred

here to global coordinates with (X, Y ) oriented (normal, parallel)

to the crack face]

Fig. 3 Predictions of phase field model (nonlinear elastic, β =

100) for extension of a mode I crack on plane misoriented by

an angle φ of twist or tilt from the original cleavage surface,

compared to linear elastic fracture mechanics solution (Gell and

Smith 1967; Wiederhorn 1984; Xu et al. 2003). Kφ/K0 is the

ratio of driving force for extension on the misoriented plane relative to that for φ = 0◦

of possible contributions of mode III loading. Presumably, such effects arise from the finite thickness of the

domain in the X3 direction, since far-field in-plane

boundary conditions (12) and (13) apply for pure mode

I loading of the initial notched configuration under

plane strain. For the case shown in Fig 1d, dependence

of the solution on X3 is expected regardless since the

cleavage plane normal m has a nonzero out-of-plane

component, though the predicted cusps in the crack

profile may be affected by or attributed to mode IIItype contributions. Analytical solutions to which the

simulation results are compared in Fig. 3, as stated

in Wiederhorn (1984) and derived in Gell and Smith

(1967), Xu et al. (2003) do not account for finite thickness effects. Their agreement with initiation predicted

by the phase field simulations is reasonable, however,

and it is noted that non-uniformity appears less pronounced at the onset of crack growth. The present

123
Phase field theory of brittle fracture 145

phase field theory, like the cited analytical solutions,

prescribes a single fracture energy for a given material

or interface, and does not distinguish among possibly

different work of separation (toughness) values under

modes I, II, and/or III. In contrast, some cohesive models of fracture such as those described in Espinosa and

Zavattieri (2003) permit the freedom of different values of work of separation under different pure modes

of loading as well as under mixed mode conditions.

4.2 Crack interaction with a spherical inclusion

A second set of phase field simulations considers the

presence of an elastic inclusion embedded in an otherwise homogeneous matrix material of the same geometry (i.e., with notch or pre-crack of length a = 50l), and

subjected to the same mode I displacement boundary

conditions as considered in Figs. 1 and 2. The spherical inclusion has initial radius RI, is centered at a distance 25l from the initial notch tip, and here obeys the

same general constitutive laws as the matrix material,

though inclusion and matrix may have different elastic stiffness and different dimensionless fracture energy

Γ ̄. Both matrix and inclusion are assigned ν = 0.25

(i.e., λ0 = μ0) and β = 0 (isotropy). Depending on the

mismatch in properties of the two materials, a propagating mode I crack will either penetrate the inclusion or

deflect around it. Representative examples of the latter

are shown in Fig. 4 (order parameter η) and Fig. 5 (normal and shear stress components for a planar slice of the

three-dimensional domain). As demonstrated in Fig. 5,

both normal and shear stresses are large but bounded

near the crack tip as it propagates around the stiff and

strong inclusion. Recall, on the other hand, that in linear elastic fracture mechanics, singularities arise at the

crack tip (with shear stress vanishing at θ = 0) for pure

mode I loading (Rice 1968). The present phase field

simulations, which omit anisotropy of fracture toughness and elastic constants, are physically representative of a brittle glassy material with perfect bonding

between inclusion and matrix, for example.

Analytical solutions for crack–inclusion interactions (Tamate 1968; Atkinson 1972; Erdogan et al.

1974) idealize the problem as 2D and linear elastic, but

do enable prediction of effects of elastic properties and

geometry on local stress intensity factors. These works

all predict a decrease in crack driving force as the stiffness of the inclusion increases. A similar effect is preFig. 4 Phase field prediction of crack deflection around a strong

second-phase inclusion of radius RI = 7.5l for far-field mode I

loading (Γ /Γ0 = 5, λ/λ0 = 1) with failed material (η > 0.7)

removed to visualize crack propagation

dicted by the phase field model, as shown in Fig. 6: for

the same far-field boundary conditions, the crack moves

closer to the soft inclusion than the stiff one. The present

phase field approach predicts evolving crack geometry

(i.e., growth) as the specimen deforms to possibly large

local strains, whereas analytical results (Tamate 1968;

Atkinson 1972; Erdogan et al. 1974) only indicate a

tendency for crack extension from a pre-existing flaw

of fixed geometry.

Numerous simulations of this inclusion–crack interaction problem are performed, wherein elastic stiffness

λ and fracture energyΓ of the inclusion are varied, with

analogous properties λ0 = μ0 and Γ0 of the matrix

held fixed, and with the same far-field mode I boundary conditions imposed in each case. Results enable

quantification of effects of properties on propagation

behavior (i.e., deflection/bridging versus particle frac123
146 J. D. Clayton, J. Knap

Fig. 5 Phase field predictions for stresses along mid-plane X3 =

0 for crack deflection around a strong and stiff second-phase

inclusion (RI/l = Γ /Γ0 = λ/λ0 = 5) with failed material

(η > 0.7) removed to visualize crack propagation: a normal

stress component P11/μ0 = Px X /μ0 b shear stress component

P12/μ0 = PxY /μ0 [referred here to global coordinates with

(X, Y ) oriented (normal, parallel) to the initial crack face]

ture) shown in Fig. 7 for RI = 5l. Such trends were also

found to be nearly identical for a larger inclusion with

RI = 7.5l (not shown). The tendency for crack deflection to become more favorable as Γ /Γ0 increases is

understandable since surface energy increases linearly

with this ratio for a planar crack splitting the inclusion. Regarding the effect of λ/λ0, phase field predictions in Fig. 7 agree with trends derived analytically via

linear elastic fracture mechanics (He and Hutchinson

1989) for penetration or deflection of a straight crack

by a planar interface separating two isotropic linear

elastic materials: deflection becomes more favorable

as the modulus of the potentially penetrated material

increases.

5 Conclusions

A new phase field theory of fracture has been implemented in 3D simulations, with results validated for

isotropic and anisotropic cracking in a homogeneous

material under distinct mode I and mode II far-field

loading conditions. Potential crack deflection around a

second-phase inclusion has been studied, demostrating

effects of relative stiffness and strength of the inclusion to that of the surrounding matrix material. These

results are thought profound with regards to the possible design of tailored nanocomposites with increased

macroscopic toughness associated with crack bridging: second-phase particles should be selected with the

appropriate properties to induce deflection, noting that

a trade-off usually exists between stiffness and toughFig. 6 Phase field prediction of crack attraction a towards soft

(λ/λ0 = 0.5) second-phase inclusion and shielding b by stiff

(λ/λ0 = 10) second-phase inclusion. In each case, the same farfield mode I boundary conditions are applied, RI = 5l, Γ /Γ0 =

1, and failed material (η > 0.7) is removed to visualize crack

propagation

ness of candidate particulate materials (e.g., stiff but

weak ceramic inclusions versus compliant yet tough

polymer inclusions). In other words, the present study

of crack penetration versus deflection through inclusions illustrates how the phase field method might be

applied to gain understanding of how nano-composites

may be tailored for improved overall toughness. Future

work shall consider additional geometries and property sets, the latter allowing for interfacial decohesion

(Evans 1974; Xu and Needleman 1993) and cleavage

anisotropy (Schultz et al. 1994), both of which can be

accommodated by the present theory.

123
Phase field theory of brittle fracture 147

Fig. 7 Predicted tendency for mode I crack to penetrate (contour

value 0) versus deflect (contour value 1) around spherical inclusion (radius RI = 5l) with fracture energy Γ and initial stiffness

λ, where Γ0 and λ0 are corresponding fixed properties of the

matrix. Partial deflection is denoted by contour values between

0 and 1

The phase field method is not restricted to simple

geometries and elastic linearity typical of known analytical solutions; crack propagation is insensitive to

mesh construction; and no spurious parameters controlling crack propagation direction or opening displacement are required. Since ceramic crystals of interest

such as alumina and silicon carbide exhibit competition

among fracture, twinning, and glide of partial dislocations (Clayton 2009, 2010), future work will address

such mechanisms simultaneously by assigning multiple order parameters, effectively merging the present

theory with that in Clayton and Knap (2011a).

References

Abdollahi A, Arias I (2012) Numerical simulation of intergranular and transgranular crack propagation in ferroelectric polycrystals. Int J Fract 174:3–15

Alber HD (2012) A model for brittle fracture based on the hybrid

phase field model. Contin Mech Thermodyn 24:391–402

Amor H, Marigo JJ, Maurini C (2009) Regularized formulation of the variational brittle fracture with unilateral contact:

numerical experiments. J Mech Phys Solids 57:1209–1229

Arias I, Knap J, Chalivendra V, Hong S, Ortiz M, Rosakis A

(2007) Numerical modelling and experimental validation of

dynamic fracture events along weak planes. Comput Methods

Appl Mech Eng 196:3833–3840

Atkinson C (1972) The interaction between a crack and an inclusion. Int J Eng Sci 10:127–136

Borden M, Verhoosel C, Scott M, Hughes T, Landis C (2012)

A phase-field description of dynamic brittle fracture. Comput

Methods Appl Mech Eng 217:77–95

Clayton J (2005) Dynamic plasticity and fracture in high density

polycrystals: constitutive modeling and numerical simulation.

J Mech Phys Solids 53:261–301

Clayton J (2006) Continuum multiscale modeling of finite deformation plasticity and anisotropic damage in polycrystals. Theoret Appl Fract Mech 45:163–185

Clayton J (2008) A model for deformation and fragmentation in

crushable brittle solids. Int J Impact Eng 35:269–289

Clayton J (2009) A continuum description of nonlinear elasticity,

slip and twinning, with application to sapphire. Proc R Soc

Lond A 465:307–334

Clayton J (2010) Deformation, fracture, and fragmentation in

brittle geologic solids. Int J Fract 163:151–172

Clayton J (2010) Modeling nonlinear electromechanical behavior of shocked silicon carbide. J Appl Phys 107:013,520

Clayton J (2011) Nonlinear Mechanics of Crystals. Springer,

Dordrecht

Clayton J (2014) Finite strain analysis of shock compression of

brittle solids applied to titanium diboride. Int J Impact Eng

73:56–65

Clayton J, Knap J (2011) A phase field model of deformation

twinning: nonlinear theory and numerical simulation. Physica

D 240:841–858

Clayton J, Knap J (2011) Phase field modeling of twinning in

indentation of transparent crystals. Model Simul Mater Sci

Eng 19:085,005

Clayton J, Knap J (2013) Phase-field analysis of fracture-induced

twinning in single crystals. Acta Mater 61:5341–5353

Clayton J, McDowell D (2003) Finite polycrystalline elastoplasticity and damage: multiscale kinematics. Int J Solids Struct

40:5669–5688

Clayton J, McDowell D (2004) Homogenized finite elastoplasticity and damage: theory and computations. Mech Mater

36:799–824

Del Piero G, Lancioni G, March R (2007) A variational model

for fracture mechanics: numerical experiments. J Mech Phys

Solids 55:2513–2537

Eastgate L, Sethna J, Rauscher M, Cretegny T (2002) Fracture

in mode I using a conserved phase-field model. Phys Rev E

65:036,117

Erdogan F, Gupta G, Ratwani M (1974) Interaction between a

circular inclusion and an arbitrarily oriented crack. J Appl

Mech 41:1007–1013

Espinosa H, Zavattieri P (2003) A grain level model for the study

of failure initiation and evolution in polycrystalline brittle

materials. Part I: theory and numerical implementation. Mech

Mater 35:333–364

Evans A (1974) The role of inclusions in the fracture of ceramic

materials. J Mater Sci 9:1145–1152

Faber K, Evans A (1983) Intergranular crack-deflection toughening in silicon carbide. J Am Ceram Soc 66:C94–C96

Foulk J, Vogler T (2010) A grain-scale study of spall in brittle

materials. Int J Fract 163:225–242

Gell M, Smith E (1967) The propagation of cracks through grain

boundaries in polycrystalline 3% silicon–iron. Acta Metal

15:253–258

Gilman J (1960) Direct measurements of the surface energies of

crystals. J Appl Phys 31:2208–2218

Gilman J (2003) Electronic basis of the strength of materials.

Cambridge University Press, Cambridge

123
148 J. D. Clayton, J. Knap

Hakim V, Karma A (2009) Laws of crack motion and phase-field

models of fracture. J Mech Phys Solids 57:342–368

He MY, Hutchinson J (1989) Crack deflection at an interface between dissimilar elastic materials. Int J Solids Struct

25:1053–1067

Hofacker M, Miehe C (2012) Continuum phase field modeling

of dynamic fracture: variational principles and staggered FE

implementation. Int J Fract 178:113–129

Jin Y, Wang Y, Khachaturyan A (2001) Three-dimensional phase

field microelasticity theory and modeling of multiple cracks

and voids. Appl Phys Lett 79:3071–3073

Kim BY, Wakayama S, Kawahara M (1996) Characterization of

2-dimensional crack propagation behavior by simulation and

analysis. Int J Fract 75:247–259

Knap J, Ortiz M (2003) Effect of indenter-radius size on Au(001)

nanoindentation. Phys Rev Lett 90:226,102

Knap J, Sieradzki K (1999) Crack tip dislocation nucleation in

FCC solids. Phys Rev Lett 82:1700–1703

Kuhn C, Muller R (2010) A continuum phase field model for

fracture. Eng Fract Mech 77:3625–3634

Lawn B (1968) Hertzian fracture in single crystals with the diamond structure. J Appl Phys 39:4828–4836

Lawn B, Padture N, Cai H, Guiberteau F (1994) Making ceramics

“ductile”. Science 263:1114–1116

Miehe C, Welschinger F, Hofacker M (2010) Thermodynamically consistent phase-field model of fracture: variational principles and multi-field FE implementations. Int J Numer Methods Eng 83:1273–1311

Moes N, Dolbow J, Belytschko T (1999) A finite element method

for crack growth without remeshing. Int J NumerMethods Eng

46:131–150

Murrell S (1965) The effect of triaxial stress systems on the

strength of rocks at atmospheric temperatures. Geophys J R

Astr Soc 10:231–281

Nuismer R (1975) An energy release rate criterion for mixed

mode fracture. Int J Fract 11:245–250

Ohji T, Jeong YK, Choa YH, Niihara K (1998) Strengthening

and toughening mechanisms of ceramic nanocomposites. J

Am Ceram Soc 81:1453–1460

Rice J (1968) Mathematical analysis in the mechanics of fracture.

In: Liebowitz H (ed) Fracture: an advanced treatise. Academic

Press, New York, pp 191–311

Schultz R, Jensen M, Bradt R (1994) Single crystal cleavage of

brittle materials. Int J Fract 65:291–312

Spatschek R, Brener E, Karma A (2011) Phase field modeling of

crack propagation. Philos Mag 91:75–95

Sun X, Khaleel M (2004) Modeling of glass fracture damage

using continuum damage mechanics-static spherical indentation. Int J Damage Mech 13:263–285

Tamate O (1968) The effect of a circular inclusion on the stresses

around a line crack in a sheet under tension. Int J Fract Mech

4:257–266

Voyiadjis G, Mozaffari N (2013) Nonlocal damage model using

the phase field method: theory and applications. Int J Solids

Struct 50:3136–3151

Wiederhorn S (1984) Brittle fracture and toughening mechanisms in ceramics. Ann Rev Mater Sci 14:373–403

Wilshaw T (1971) The Hertzian fracture test. J Phys D Appl Phys

4:1567–1581

Xu L, Huang Y, Rosakis A (2003) Dynamic crack deflection and

penetration at interfaces in homogeneous materials: experimental studies and model predictions. J Mech Phys Solids

51:461–486

Xu XP, Needleman A (1993) Void nucleation by inclusion

debonding in a crystal matrix. Model Simul Mater Sci Eng

1:111–132

Zhang S, Zhu T, Belytschko T (2007) Atomistic and multiscale

analyses of brittle fracture in crystal lattices. Phys Rev B

76:094,114

123
11

1 DEFENSE TECHNICAL

(PDF) INFORMATION CTR

DTIC OCA

2 DIRECTOR

(PDF) US ARMY RESEARCH LAB

RDRL CIO LL

IMAL HRA MAIL & RECORDS MGMT

1 GOVT PRINTG OFC

(PDF) A MALHOTRA

37 DIR USARL

(PDF) RDRL CIH C

J KNAP

RDRL WM

B FORCH

J MCCAULEY

RDRL WML B

I BATYREV

B RICE

D TAYLOR

N WEINGARTEN

RDRL WML H 

C MEYER

B SCHUSTER

RDRL WMM

J BEATTY

RDRL WMM B

G GAZONAS

D HOPKINS

B POWERS

C RANDOW

T SANO

RDRL WMM E

J SWAB

RDRL WMM F

M TSCHOPP

RDRL WMM G

J ANDZELM

RDRL WMP

S SCHOENFELD

RDRL WMP B

C HOPPEL

S SATAPATHY

M SCHEIDLER

A SOKOLOW

T WEERASOORIYA

RDRL WMP C

R BECKER

S BILYK

T BJERKE

D CASEM

J CLAYTON 

D DANDEKAR

M GREENFIELD

R LEAVY

J LLOYD

S SEGLETES

A TONGE

C WILLIAMS

RDRL WMP D

R DONEY
12

INTENTIONALLY LEFT BLANK.