

Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

49, 4, pp. 1115-1133, Warsaw 2011

A NEW FINITE ELEMENT TECHNIQUE FOR A PHASE

FIELD MODEL OF BRITTLE FRACTURE

Charlotte Kuhn

Ralf Müller

University of Kaiserslautern, Institute of Applied Mechanics, Kaiserslautern, Germany

e-mail: chakuhn@rhrk.uni-kl.de; ram@rhrk.uni-kl.de

Phase field models for fracture employ a continuous field variable to
indicate cracks. The width of the transition zone between cracked and
uncracked areas is controlled by a regularization parameter. The nu-
merical implementation of such models is sensible to the choice of this
parameter in conjunction with the mesh size, as the mesh has to be fi-
ne enough to resolve high gradients of the crack field appearing in the
transition zones. This is the main computational limit and challenge of
the implementation. To overcome this limitation, a finite element me-
thod using exponential shape functions is introduced. Numerical exam-
ples show that these new shape functions perform better than standard
Lagrange shape functions.

Key words: phase field, fracture, finite elements, exponential shape
functions

1. Introduction

In phase field models for fracture, such as introduced by e.g. Aranson et al.
(2000), Karma et al. (2001), Eastgate et al. (2002), Brener and Spatschek
(2003) and Spatschek et al. (2006), cracks are represented by an additional
field variable which is 0 if the material is cracked and 1 if it is undamaged,
and cracking is addressed as a phase transition problem. All these phase field
fracture models differ in detail, but commonly they feature a regularization
parameter which has the dimension of length and controls the width of the
transition zone where the crack field interpolates between 1 and 0. Finite
element discretizations of those models are faced with the difficulty that the
mesh has to be sufficiently fine in relation to the regularization length, as it
has to be fine enough to resolve high gradients of the crack field appearing


1116 C. Kuhn, R. Müller

around the cracks. On the other hand, the regularization parameter must be
chosen sufficiently small in order to obtain reasonable results. Amor et al.
(2009) suggest 1% of the global geometric dimension of the sample to be an
appropriate value for the regularization parameter.

There are different approaches to meet the requirements for a sufficiently
fine resolution on the one hand and to keep the computation time within
bounds on the other hand. Eastgate et al. (2002) use Fourier transforms to
solve the linear terms of their phase field model in order to increase the effi-
ciency of the computations.However, this restricts the simulations toproblems
with periodic boundary conditions.Amore commonapproach are adaptive re-
meshing strategies as can be found in Provatas et al. (1998) for a phase field
solidification problem or in Bourdin andChambolle (2000) for an approxima-
tion of theMumford-Shah functional. Exploiting the fact that the phase field
order parameter varies significantly only near an interface, themesh is refined
only where it is needed. This also applies for the material force based h-type
mesh refinement algorithm suggested by Welschinger et al. (2010).

In the presentwork, we follow a different approach inspired byLaZghab et
al. (2002). Special shape functions that qualitatively capture the shape of the
crackfieldare constructedand implemented into thefiniteelement codeFEAP.
These shape functions allow a coarser discretization without compromise on
the accuracy of the results.

2. Material model

2.1. Governing equations

The core of the present phase field model of fracture is a regularized ver-
sion of the variational formulation of brittle fracture by Francfort andMarigo
(1998), whichwas introduced inBourdin (2007). In the regularizedmodel, the
energy density of a broken linear elastic material with the stiffness tensor C
and the cracking resistance Gc is approximated by

ψ(ε,s)= ψel(ε,s)+ψsurf (s)

=
1

2
(s2+η)ε · (Cε)+Gc

( 1

4ǫ
(1−s)2+ ǫ|∇s|2

) (2.1)

This regularized energy density ψ is a function of the linearized strain tensor
ε = 1

2
(∇u+(∇u)⊤), i.e. the symmetric part of the gradient of the displa-

cements u and the continuous scalar field s, which describes the condition


A new finite element technique for a phase field model... 1117

of the material. Cracks are indicated by s = 0, while s = 1 is assigned to
sound material. The parameter ǫ, appearing twice in the surface part ψsurf ,
has the dimension of length and regulates the width of the transition zone be-
tween broken and unbrokenmaterial, where s interpolates between 0 and 1.
The factor (s2+η) in the elastic part of the energy density ψel accounts for
the change in stiffness between broken and unbroken material. Consequently,
material law (2.2) for the Cauchy stress tensor σ is modified to

σ=
∂ψ

∂ε
=(s2+η)Cε (2.2)

Thus, the parameter η ≪ 1 can be seen as the remaining stiffness if s = 0
and is needed to secure positive definiteness of the elastic energy. The local
balance law for the Cauchy stress tensor

divσ=0 (2.3)

where body forces and inertia terms are neglected, remains unaltered upon
the regularization.

Interpreting s as the order parameter of a phase fieldmodel, ṡ is assumed
to be proportional to the variational derivative of the energy density ψ with
respect to s, i.e.

ṡ =−M
δψ

δs
=−M

[

sε · (Cε)−Gc
(

2ǫ∆s+
1−s

2ǫ

)]

(2.4)

Themobility factor M is a constant, which describes the kinetics of the pro-
cess. For sufficiently large values of M, the solution of the evolution equation
can be considered as stationary. In order to take into consideration the irre-
versible character of cracking, s(x, t) is fixed to 0 for all future times t > t∗

if it becomes 0 at any time t∗.

2.2. 1D solution

The study of evolution equation (2.4) in one dimension is a very simple yet
instructive way to understand the impact of the regularization parameter ǫ
on the solution. Neglecting the elastic contribution and assuming the solution
to be stationary (ṡ =0), the one dimensional evolution equation reduces to

s′′−
s

4ǫ2
=−

1

4ǫ2
(2.5)


1118 C. Kuhn, R. Müller

If boundary conditions s(0) = 0 and s′(±L) = 0 apply, Eq. (2.5) yields the
(piecewise defined) solution

s(x)= 1−
cosh
(
L−|x|
2ǫ

)

cosh
(
L
2ǫ

) (2.6)

which converges to

s(x)= 1− exp
(

−
|x|

2ǫ

)

(2.7)

for L ≫ ǫ. Figure 1 shows plots of this function for different values of ǫ. The
smaller ǫ gets, the higher gradients and curvatures of the solution s(x) appear
in the vicinity of the crack at x =0. The limit ǫ → 0 yields a discontinuous
function, which is 0 at x =0 and 1 elsewhere.

Fig. 1. 1D stationary crack field

3. Numerical implementation

3.1. Finite element formulation

The set of equations formed by balance equation (2.3) and evolution equ-
ation (2.4) together with the respective boundary conditions is solved using
the finite elementmethod.For the 2D case, 4 nodequadrilateral elementswith
3 degrees of freedom (ux,uy,s) per node are used in the discretization. The
starting point of the FE formulation is the weak form of the field equations.
With virtual displacements δu and δs, they read

∫

Ω

∇δu ·σ dV =

∫

∂Ωt

δu · t∗n dA (3.1)


A new finite element technique for a phase field model... 1119

and
∫

Ω

[

δs
ṡ

M
−∇δs ·q+ δs

(

sε · (Cε)+
Gc
2ǫ
(s−1)

)]

dV =

∫

∂Ωq

δsq∗n dA (3.2)

with q = −2Gcǫ∇s. The boundary conditions for the stresses σ and for q
are prescribed by the traction t∗n and the normal flux q

∗
n, which is usually

assumed to be zero.
In the discretization, the displacements u, the crack field s, as well as

their virtual counterparts δu and δs are approximated by shape functions
NuI , N

s
I , N

δu
I , and N

δs
I , which interpolate the respective nodal values ûI, ŝI,

δûI, and δŝI.UsingVoigt-notation –denotedbyanunderbar in the following–
the approximations read

u=
N
∑

I=1

NuI ûI s =
N
∑

I=1

NsI ŝI

δu=
N
∑

I=1

NδuI δûI δs =
N
∑

I=1

NδsI δŝI

(3.3)

With thematrices

[BuI ] =







NuI,x 0

0 NuI,y
NuI,y N

u
I,x







[BsI] =

[

NsI,x
NsI,y

]

[BδuI ] =







NδuI,x 0

0 NδuI,y
NδuI,y N

δu
I,x







[BδsI ] =

[

NδsI,x
NδsI,y

]

(3.4)

defined by the derivatives of the shape functions, the approximations of the
gradient expressions yield

ε=
N
∑

I=1

[BuI ]ûI ∇s =
N
∑

I=1

[BsI]ŝI

δε=
N
∑

I=1

[BδuI ]δûI ∇δs =
∑N
I=1[B

δs
I ]δŝI

(3.5)

AsEq. (3.1) andEq. (3.2)must hold for any choice of the virtual quantities
δu and δs, the respective nodal values cancel out of the system of equations.
Thus, the left hand sides of Eq. (3.1) and Eq. (3.2) form the nodal residuals


1120 C. Kuhn, R. Müller

[RI] =

[

RuI
RsI

]

=

∫

Ω





[BδuI ]
⊤σ

NδsI
ṡ

M
− [BδsI ]

⊤q+NδsI

(

sε⊤ · [Cε]+
Gc
2ǫ
(s−1)

)



 dV

(3.6)
Thetime integration of the transient terms isperformedwith the implicitEuler
method, and the overall nonlinear set of equations is solved iteratively with a
Newton-Raphson algorithm.For this algorithm, the stiffnessmatrix [KIJ] and
the dampingmatrix [DIJ] must be provided. They are obtained by derivation
of the nodal residuals [RI] with respect to the nodal values (ûJ, ŝJ)

[KIJ]=

∫

Ω





[BδuI ]
⊤(s2+η)C[BuJ] [B

δu
I ]
⊤2sCεNsJ

NδsI 2s(Cε)
⊤[BuJ] 2Gcǫ[B

δs
I ]
⊤[BsJ]+N

δs
I

(

ε⊤·Cε+
Gc
2ǫ

)

NsJ



dV

(3.7)

DIJ =

∫

Ω





0 0

0
1

M
NδsI N

s
J



dV

If the same shape functions are chosen for the approximation of the actual
values and thevirtual quantities, i.e. NuI = N

δu
I and N

s
I = N

δs
I , the assembled

system matrix becomes symmetric; differently chosen shape functions render
a non-symmetric system matrix.

4. Exponential shape functions

Usually, the linear Lagrangian shape functions

NlI(ξ,η)=
1

4
(1+ ξIξ)(1+ηIη) I =1, . . . ,4 (4.1)

with (ξI,ηI) according to Fig.2 are used for all the shape functions N
u
I , N

δu
I ,

NsI , and N
δs
I as well as for the approximation of the geometry within the

isoparametric concept

x=
N
∑

I=1

NlIx̂I (4.2)

The computational limit of this method becomes apparent looking at the
3D representation of the phase field of an initial crack in an unloaded sample
in Fig.3. The mesh has to be chosen fine enough to resolve the high gra-
dients and curvatures of the phase field in the transition zone between cracked
and uncracked areas. In this example, the regularization parameter ǫ was set


A new finite element technique for a phase field model... 1121

Fig. 2. Node and edge numbering of the quadrilateral element in global (left) and
natural coordinates (right)

Fig. 3. 3D representation of the phase field of a crack

to 0.01L, where L is the edge length of the square sample. A uniform mesh
with 200×200 elements was used for the discretization, thus the edge length
of the elements is h =0.005L. Bourdin et al. (2008) show that linear triangu-
lar elements overestimate the surface energy by a factor f(h/ǫ) = 1+h/4ǫ,
and hence different authors (Amor et al., 2009;Miehe et al., 2010) empirically
found h ≈ ǫ as an upper bound for the element size in a 2D setting.

Figure 4 compares the computed phase field values to analytical solu-
tion (2.7) at the specimen edge far behind the crack tip along the segment
{0}× [0,0.5L] and in front of the crack tip along the segment [0.5L,1L]×{0}
for two different values of ǫ. At the specimen edge, the analytical solution
captures the computed values quite perfectly, yet in front of the crack tip the
computed crack field is even steeper than the analytical solution predicts. Ho-
wever, the shape still resembles an exponential function, and the impact of ǫ
remains the same, i.e. the smaller ǫ, the smaller the transition zone between


1122 C. Kuhn, R. Müller

s =0 and s =1. Hence, particularly along the crack lips but also at the crack
tip, the shape functions derived from the 1D solution promise a very accurate
interpolation of the crack field.

Fig. 4. Comparison of the computed 2D phase field to 1D solution (2.7) at the
specimen edge (top) and at the crack tip (bottom)

4.1. Construction of 1D exponential shape functions

The simulation of extrusion processes, where shear boundary layers cha-
racterized by an exponential velocity profile need to be resolved, is faced with
a similar problem.As an alternative to a finenumerical resolution, LaZghab et
al. (2002) introduce exponential shape functions that capture the shape of the
sharp velocity field. Since the numerical 2D solution shows the same characte-
ristic exponential shape thatwas found for the 1D case, this ansatz seems also
very promising for the discretization of the present phase fieldmodel.Adapted
to the present phase field model, these 1D shape functions read

N̄e1(ξ,δ) = 1−
exp
(

−
δ(1+ξ)
4

)

−1

exp
(

−δ
2

)

−1

N̄e2(ξ,δ) =
exp
(

−
δ(1+ξ)
4

)

−1

exp
(

−δ
2

)

−1

(4.3)

in natural coordinates on the interval [−1,1]. The coefficient δ = h/ǫ in
the exponential shape functions is the ratio of the element size h and the


A new finite element technique for a phase field model... 1123

regularizationparameter ǫ. In the limit case δ → 0or equivalently h → 0 these
shape functions converge to the one dimensional linear shape functions, i.e.

lim
δ→0

N̄e1(ξ,δ) =
1−ξ

2
= Nl1(ξ) lim

δ→0
N̄e2(ξ,δ) =

1+ ξ

2
= Nl2(ξ) (4.4)

These shape functions are designed to perfectly match analytical solu-
tion (2.7), if ŝ1 ¬ ŝ2 holds for the according nodal values. However, they do
not match in the other case, see Fig.5. This problem arises due to the fact
that the exponential shape functions are unsymmetric with respect to ξ, i.e.
N̄e1(−ξ) 6= N̄

e
2(ξ). In order to resolve this deficiency, the orientation of the

exponential shape functions needs to be switched according to the nodal va-
lues of s

Ne1(ξ)=

{

N̄e1(ξ,δ) if ŝ1 ¬ ŝ2

N̄e2(−ξ,δ) if ŝ1 > ŝ2

Ne2(ξ)=

{

N̄e2(ξ,δ) if ŝ1 ¬ ŝ2

N̄e1(−ξ,δ) if ŝ1 > ŝ2

(4.5)

Fig. 5. Approximation with unswitched shape functions (left), and switched shape
functions (right)

The respective derivatives which are necessary for the approximation of
the gradient ∇s in the finite element discretization are

∂Ne1
∂ξ
=
δ

4

exp
(

−
δ(1±ξ)
4

)

exp
(

−δ
2

)

−1

∂Ne2
∂ξ
=−

δ

4

exp
(

−
δ(1±ξ)
4

)

exp
(

−δ
2

)

−1
(4.6)

where + applies for ŝ1 ¬ ŝ2 and − for ŝ1 > ŝ2.


1124 C. Kuhn, R. Müller

4.2. Extension to 2D

Thevelocity field forwhich the exponential shape functionswere construc-
ted by LaZghab et al. (2002) is exponential only in one direction. Thus, it is
sufficient to combine the 1D exponential shape functionswith linear 1D shape
functions in the second direction in order to obtain shape functions for a 2D
quadrilateral element. Inour case, thiswouldbesufficient for thediscretization
of the crack lips but not at the crack tip. Thereforewe use the 1D exponential
shape functions in both directions to construct 2D shape functions.
The 2D linear shape function of each single element node can be obtained

bymultiplying the 1D linear shape functions belonging to the adjacent edges
of the respective node. The adaptation of this strategy for the construction of
2D exponential shape functions yields

N̄e1(ξ,η,δi)= N
e
1(ξ,δ1)N

e
1(η,δ4) N̄

e
2(ξ,η,δi)= N

e
2(ξ,δ1)N

e
1(η,δ2)

N̄e3(ξ,η,δi)= N
e
2(ξ,δ3)N

e
2(η,δ2) N̄

e
4(ξ,η,δi)= N

e
1(ξ,δ3)N

e
2(η,δ4)
(4.7)

where the element nodes and the element edges are numbered according to
Fig.2. Each shape function depends on the ratio δi = hi/ǫ of both adjacent
element edges. For an appropriate approximation behavior, it is postulated
that the orientation of the 1D shape functions of opposite edges must be
the same. The so constructed shape functions possess the Kronecker delta
property, i.e. N̄eI(ξJ,ηJ,δi) = δIJ. Continuity across element borders holds,
if the orientation of the shared edge of two neighbor elements is the same.
However, summation of the 4 shape functions gives

4
∑

I=1

N̄eI(ξ,η,δi)= 1− [N
e
1(ξ,δ1)−N

e
1(ξ,δ3)][N

e
1(η,δ2)−N

e
1(η,δ4)]

︸ ︷︷ ︸

=R(ξ,η,δi)

(4.8)

Thus, partition of unity does not hold in every case. If arbitrary quadrilateral
elements are used, the shape functions can bemodified, for example, to

NeI(ξ,η,δi)= N̄
e
I(ξ,η,δi)+

1

4
R(ξ,η,δi) for I =1, . . . ,4. (4.9)

Under the constraint that the orientation of the 1Dshape functions of opposite
edges must be the same, the correction term R(ξ,η,δi) vanishes whenever

δ1 = δ3 or δ2 = δ4 (4.10)

holds. This is especially true for square and rectangular elements as used in
Section 5.


A new finite element technique for a phase field model... 1125

As R(ξ,η,δi) vanishes on the entire boundary of the unit square
[−1,1] × [−1,1], this correction does not disturb the continuity and the
Kronecker delta property of the shape functions. Furthermore, one can show
that limδi→0R(ξ,η,δi)= 0.

4.3. Derivatives of the 2D exponential shape functions

The computation of thederivatives NeI,ξ and N
e
I,η of the exponential shape

functions with respect to the natural coordinates ξ and η is straightforward,
using the 1D derivatives in Eq. (4.6). The derivatives NeI,x and N

e
I,y with

respect to the global coordinates x and y follow from the relation

∂NeI
∂ξ
=
∂NeI
∂x

∂x

∂ξ
+
∂NeI
∂y

∂y

∂ξ

∂NeI
∂η
=
∂NeI
∂x

∂x

∂η
+
∂NeI
∂y

∂y

∂η

⇔

[

NeI,ξ
NeI,η

]

=

[

x,ξ y,ξ
x,η y,η

]

︸ ︷︷ ︸

=J

[

NeI,x
NeI,y

]

(4.11)

where

J=
4
∑

I=1

[

NlI,ξx̂I N
l
I,ξŷI

NlI,ηx̂I N
l
I,ηŷI

]

(4.12)

if the geometry is approximated with linear shape functions according to
Eq. (4.2).

5. Results

The performance of the 2D exponential shape functions is tested in this sec-
tion. In all simulations, linear shape functions (4.1) were used for the ap-
proximation of the geometry and the actual and virtual displacements, i.e.
NuI = N

δu
I = N

l
I. Three different versions of approximating the crack field s

and its virtual counterpart δs are compared to each other: The standard ap-
proximationwith linear shape functions NsI = N

δs
I = N

l
I (labeled lin/lin), the

complete approximation with exponential shape functions NsI = N
δs
I = N

e
I

(labeled exp/exp), and amixed formulationwith NsI = N
e
I but N

δs
I = N

l
I (la-

beled lin/exp). As regular meshes with solely rectangular or square elements
were used, condition (4.10) holds and the correction term R of the exponential
shape function vanishes.


1126 C. Kuhn, R. Müller

5.1. Stationary evolution equation

For the first numerical assessment of the 2D exponential shape functions,
we stick to the example of Section 4.2. No mechanical loads are applied and
the problem reduces to solving evolution equation (2.4) under the constraint
s(x,y) = 0 if (x,y)∈ [0,L/2]×{0}. A regular mesh with square elements is
used for thediscretization. Figure 6 shows the crack field for ǫ =0.01L compu-
tedwith only 4×4 elements (δi =25) and 5Gauss points per direction.With
such a coarsemesh, the standard approach using linear shape functions fails to
give a reasonable solution. However, both computations using the exponential
shape functions yield already qualitatively very good results, which are at first
sight almost identical. Only directly in front of the crack tip the result from
themixed formulation (lin/exp) is even a bit more accurate.

Fig. 6. Comparison of the crack field computed with 4×4 elements

Figure 7 shows the results from a further investigation of the performan-
ce of the exponential shape functions in this example. For ǫ = 0.01L (left)
and ǫ =0.002L (right), the stationary solution of the evolution equation was
computed with meshes within the range of 2×2 to 400×400 elements. The
surface energy

Esurf =

∫

Ω

ψsurf (s) dV =

∫

Ω

Gc
( 1

4ǫ
(1−s)2+ ǫ|∇s|2

)

dV (5.1)


A new finite element technique for a phase field model... 1127

Fig. 7. Evaluation of the surface energy

associatedwith the computed crack field is plotted in the left column of Fig.7.
The results obtainedwith both versions using the exponential shape functions
are very similar. Even if a very coarsemesh is used, the surface energy is only
slightly underestimated, while the solution with linear shape functions ove-
restimates it by far. This is even more significant for the smaller value of ǫ.


1128 C. Kuhn, R. Müller

In standard 2D finite element analysis, the integrals in residuals (3.6), stiff-
ness matrix (3.7)1, and dampingmatrix (3.7)2 are computed using the Gauss
quadrature formula with 2 integration points per direction. This quadrature
methodwas used for the computation of the surface energy in the first rowand
the relative error of the surface energy in the second row of Fig.7. The usage
of the exponential shape functions yields a smaller relative error, yet their full
potential only reveals itself if a more precise quadrature method is employed.
The plots in the last row of Fig.7 show the relative error in the surface energy,
if 10× 10 Gauss points are used for the integration. Especially for large va-
lues of the ratio δ = h/ǫ, this measure drastically decreases the relative error
compared to the calculations with 2×2 quadrature points. Thus, in the latter
case, the largest part of the error is due to the quadrature method. Standard
linear shape functions and a non-uniform mesh with square elements of edge
length h = 7.1429 · 10−4L and h = 4.8828 · 10−4L in the vicinity the crack,
were used to obtain the reference solutions Esurf = 0.51017344300GcL for
ǫ =0.01L and Esurf =0.50241252899GcL for ǫ =0.002L, respectively.

5.2. Crack growth

In this section, we test the performance of the exponential shape func-
tions, when mechanical loads are applied to the system, i.e. the whole set of
coupled equations has to be solved. Therefore, the sample is now loaded by
a linear increasing displacement load u∗(t) =

√

GcL/(2µ)t along the edge
[0,L]×{0.5L}.With this scaling of the displacements, the geometric length L
and the cracking resistance Gc can be factored out from the energy functional.
If themobility M is chosen large enough to assume quasistatic cracking, then
the solution of the coupled problem only depends on the ratio of the Lamé
constants λ/µ and the regularization parameter ǫ in conjunctionwith L. The
simulations presented here, refer to the case of equal Lamé constants λ = µ
with plane strain assumptions, corresponding to a Poisson’s ratio of ν =0.25
and the regularization parameter ǫ =0.0005L. Exploiting the symmetryof the
problem, only half of the area is discretized. The discretization in x-direction
is done with 150 elements. A varying number of n elements plus one row of
elements of fixed height, to model the initial crack, discretize the structure in
the y-direction, see Fig.8. TheFEMsolutions based on the exponential shape
functions in Section 5.1 were very similar, no matter if linear or exponential
shape functions were used for the discretization of δs. Themixed formulation
yields an unsymmetric system matrix, which is computationally more expen-
sive, and thus this option is not treated in the following. Gauss quadrature
with 5×5 integration points was used for the integration.


A new finite element technique for a phase field model... 1129

Fig. 8. Simulation setup: contour plot of the initial crack field (left) and finite
elementmesh (right)

The two plots in Fig.9 show the evolution of the elastic energy with re-
spect to the load factor t for different numbers of n. At the beginning, the
elastic energy increases with the loading until rupture occurs and it drops to
zero. Impressively, the simulation with only n =2 elements in the y-direction
already gives a qualitatively good result not that far away from the simula-
tionswithmore elements, when the exponential shape functions are employed.
Using the standard linear shape functions, no rupture can be observed in the
simulation with n = 2 elements up to a load factor of t = 3, which is about
twice the actual critical loading. Also the simulation with n = 16 elements
still overestimates the critical loading by far. Only the simulations with more
elements produce as accurate results as the simulations with the exponential
shape functions. Seeking for a reliable prediction of the stability of a structure,
another advantage of the exponential shape functions is that they underesti-
mate the true critical load value. The overestimation of the critical load value
of the linear shape functions stems from the overestimation of the surface
energy associated with the initial crack.

Fig. 9. Elastic energy


1130 C. Kuhn, R. Müller

6. Conclusions and outlook

The aim of this work was to provide an alternative to an expensive mesh
refinement in finite element simulations of a phase field model for fracture in
cases where the regularization parameter is very small. To this end, special
shape functions, which capture the analytical stationary solution of the 1D
crack field,were derived and implemented into a 2Delement of a finite element
code. As the exponential terms of these shape functions are controlled by the
regularization length of the phase field model, they are able to adjust to the
shape of the crack field, which depends on this parameter in a similar way.
This adaptivepropertyof the exponential shape functionsallows computations
with virtually arbitrarily small values of the regularization parameter, which
would require an extensivemesh refinement if standard linear shape functions
are used.

Theperformanceof the 2Dexponential shape functions has been evaluated
in two examples, the first considering only the stationary evolution equation,
the second considering the whole coupled problem of the mechanical force
balance and the evolution equation. In both cases, usage of the exponential
shape functions allowed a considerable reduction of the level of refinement
without compromise on the accuracy of the results.

Two different approaches of incorporating the exponential shape functions
in the finite element scheme have been tested. The mixed formulation with
exponential shape functions only for the crack field itself and linear shape
functions for the test functions seems to yield slightly more accurate results
than the purely exponential formulation. Yet the prize to pay is an unsymme-
tric systemmatrix causing a significant increase in the computation time and
required memory. In this regard, the purely exponential formulation is to be
preferred.

Another issue that comes up when the exponential shape functions are
used is the choice of an appropriate quadrature method for computation of
the integrals in the residual and the systemmatrix. The standard integration
scheme using 2Gauss quadrature points per direction fails to approximate the
integrals sufficiently well if there are large gradients in the crack field, and the
ratio δ = h/ǫ of element length and regularization parameter becomes too
large. In this case, a more exact quadrature method with more integration
points has to be employed. Yet, this again increases the computation time.
Thus, an automatic choice of a reasonable number of quadrature points in
conjunction with δ is useful.


A new finite element technique for a phase field model... 1131

So far, the simulations with the elements with exponential shape functions
are restricted to simple examples, where the crack path is known in advan-
ce. This is due to the fact that the adequate orientation of those elements
with respect to the crack position has to be defined a priori. The considered
phase fieldmodel, however, naturally contains the possibility to determine al-
so complicated crack paths including crack initiations and crack branching.
To this end, the development of a stable algorithm that, if necessary, re-
defines the orientation of the exponential elements after every time step is
essential.

The exponential shape functions are especially adequate to approximate
the crack field in fractured zones. In undamaged zones, where the crack field
is almost constant, the linear shape functions do as well. Due to the expo-
nential terms involved, the evaluation of the exponential shape functions is
computationallymore expensive than the evaluation of linear shape functions.
Therefore, they should only be employed where they are needed, i.e. in the
vicinity of cracks. A combination of both raises the problem of constructing
blending elements which blend elements with exponential shape functions to
those with linear shape functions. The difficulty here is to preserve continuity
and partition of unity properties of the shape functions.

In this work, the focus was laid on improving the shape functions for the
crack field.However, also thedisplacement field features high gradients around
the cracks, which also cannot be captured by the linear shape functions if the
mesh is too coarse. Thus, further room for improvement lies in enhancing
the shape functions used for the approximation of the displacement field. An
ansatz going in this direction would be the construction of a 9 node element,
combining the exponential shape functions with 3 nodes per direction, which
were also introduced in LaZghab et al. (2002), with quadratic Lagrange shape
functions for the displacement field.

References

1. Amor H., Marigo J.-J., Maurini C., 2009, Regularized formulation of the
variational brittle fracture with unilateral contact: Numerical experiments, Jo-
urnal of the Mechanics and Physics of Solids, 57, 8, 1209-1229

2. Aranson I.S, Kalatsky V.A., Vinokur V.M., 2000, Continuum field de-
scription of crack propagation,Physical Review Letters, 85, 1, 118-121


1132 C. Kuhn, R. Müller

3. Bourdin B., 2007, Numerical implementation of the variational formulation
of quasi-static brittle fracture, Interfaces and Free Boundaries, 9, 411-430

4. Bourdin B., Chambolle A., 2000, Implementation of an adaptive finite-
element approximation of the mumford-shah functional, Numerische Mathe-
matik, 85, 4, 609-646

5. Bourdin B., Francfort G., Marigo J.-J., 2008, The variational approach
to fracture, Journal of Elasticity, 91, 1, 5-148

6. Brener E.A., Spatschek R., 2003, Fast crack propagation by surface diffu-
sion,Physical Review E, 67, 1, 016112

7. Eastgate L.O., Sethna J.P., Rauscher M., Cretegny T., ChenC.-S.,
Myers C.R., 2002, Fracture in mode I using a conserved phase-field model,
Physical Review E, 65, 3, 036117

8. FrancfortG.A.,MarigoJ.-J., 1998,Revisitingbrittle fractureasanenergy
minimization problem, Journal of the Mechanics and Physics of Solids, 46, 8,
1319-1342

9. Karma A., Kessler D.A., Levine H., 2001, Phase-field model of mode iii
dynamic fracture,Physical Review Letters, 87, 4, 45501

10. LaZghab S., Aukrust T., Holthe K., 2002, Adaptive exponential finite
elements for the shear boundary layer in the bearing channel during extru-
sion, Computer Methods in Applied Mechanics and Engineering, 191, 11/12,
1113-1128

11. Miehe C.,Welschinger F., HofackerM., 2010,Thermodynamically con-
sistent phase-fieldmodels for fracture: Variational principles andmulti-field fe
implementations, International Journal for NumericalMethods in Engineering,
83, 10, 1273-1311

12. Provatas N., Goldenfeld N., Dantzig J., 1998, Efficient computation
of dendritic microstructures using adaptive mesh refinement, Physical Review
Letters, 80, 15, 3308-3311

13. SpatschekR.,HartmannM.,BrenerE.,Müller-KrumbhaarH.,Kas-
snerK., 2006,Phase fieldmodeling of fast crackpropagation,Physical Review
Letters, 96, 1, 015502

14. Welschinger F., Hofacker M., Miehe C., 2010, Configurational-force-
based adaptive fe solver for a phase fieldmodel of fracture, [In:]Proceedings in
Applied Mathematics and Mechanics, 10, 689-692


A new finite element technique for a phase field model... 1133

Nowa metoda elementów skończonych dla modelu pól fazowych przy

opisie kruchego pękania

Streszczenie

Modele pól fazowychw opisie procesu pękania wykorzystują zmienne ciągłe pola
do wykrywania pęknięć. Szerokość strefy przejściowej pomiędzy obszarem pęknięcia
a nieuszkodzonym jest opisana parametrem regularyzacji. Numeryczna implementa-
cja takich modeli jest wrażliwa na dobór tego parametru w połączeniu z rozmiarem
siatki elementów skończonych, któramusi być odpowiednio gęsta, by uwzględnić du-
że gradienty pola z pęknięciem w strefie przejściowej. Jest to główne ograniczenie
w przeprowadzaniu obliczeń i duże wyzwanie symulacyjne. W pracy zaproponowano
użycie wykładniczych funkcji kształtu do metody elementów skończonych w celu eli-
minacji tego ograniczenia.Przedstawioneprzykładypokazały, że zastosowanie funkcji
wykładniczych zamiast standardowych funkcji Lagrange’awyraźniepoprawiłowydaj-
ność numerycznąmodeli.

Manuscript received December 1, 2010; accepted for print January 31, 2011