Biomech Model Mechanobiol (2016) 15:479–496

DOI 10.1007/s10237-015-0702-0

ORIGINAL PAPER

A phase-field model for fracture in biological tissues

Arun Raina1 · Christian Miehe1

Received: 28 March 2015 / Accepted: 1 July 2015 / Published online: 14 July 2015
© Springer-Verlag Berlin Heidelberg 2015

Abstract This work presents a recently developed phas- Keywords Soft biological tissues · Phase-field fracture ·
e-field model of fracture equipped with anisotropic crack Anisotropic failure criterion · Coupled multi-field
driving force to model failure phenomena in soft biolog- problems
ical tissues at finite deformations. The phase-field models
present a promising and innovative approach to thermody-
namically consistent modeling of fracture, applicable to both 1 Introduction
rate-dependent or rate-independent brittle and ductile failure
modes. One key advantage of diffusive crack modeling lies Computational modeling of fracture in soft biological tis-
in predicting the complex crack topologies where methods sues presents a promising research area due to its scope
with sharp crack discontinuities are known to suffer. The to help improve the understanding behind surgical treat-
starting point is the derivation of a regularized crack sur- ments or injuries involving mechanical stresses, e.g., balloon
face functional that Γ -converges to a sharp crack topology angioplasty, ruptured aneurysm or ligament tear. An extra-
for vanishing length-scale parameter. A constitutive balance cellular matrix consisting of networks of elastin and collagen
equation of this functional governs the crack phase-field fibers surrounded by smooth muscle cells and fluids consti-
evolution in a modular format in terms of a crack driving tutes a major part of a tissue; see, e.g., Humphrey (2002).
state function. This allows flexibility to introduce alterna- Elastin networks bear long-range reversible small deforma-
tive stress-based failure criteria, e.g., isotropic or anisotropic, tions and can endure several million deformation cycles
whose maximum value from the deformation history drives without undergoing fatigue as shown in Gundiah et al.
the irreversible crack phase field. The resulting multi-field (2007). The triple-stranded helical structure of collagen fibers
problem is solved by a robust operator split scheme that endows structural integrity to soft tissues at finite deforma-
successively updates a history field, the crack phase field tions due to unraveling of its wavy form (Gelsea et al. 2003).
and finally the displacement field in a typical time step. Due to this microstructural composition, soft biological tis-
For the representative numerical simulations, a hyperelastic sues are nearly incompressible with highly nonlinear and
anisotropic free energy, typical to incompressible soft bio- anisotropic response to mechanical stimulus where more
logical tissues, is used which degrades with evolving phase details can be found in Rhodin (1980) and Chuong and
field as a result of coupled constitutive setup. A quantitative Fung (1984). The constitutive material modeling of bio-
comparison with experimental data is provided for verifica- logical tissues is a constantly evolving research area. See,
tion of the proposed methodology. e.g., Demiray (1972) and Delfino et al. (1997) where an
isotropic strain energy density for soft biological tissues is
formulated or Vaishnav et al. (1973), Lanir (1979), Fung
B Christian Miehe et al. (1979) and Billiar and Sacks (2000) which addi-
cm@mechbau.uni-stuttgart.de tionally account for material anisotropy. Another set of
1 Institute of Applied Mechanics (CE), Chair I,
approaches are based on the structural tensors for incorporat-
University of Stuttgart, Pfaffenwaldring 7, ing anisotropy into polyconvex free energy functions. See,
70569 Stuttgart, Germany e.g., Holzapfel et al. (2000), Menzel and Steinmann (2001)

123
480 A. Raina, C. Miehe

and Gasser et al. (2006) where the latter also accounts for parameters and confirmed the role of different collagen
fiber dispersion obtained from the orientation distribution fiber orientations in different layers yielding anisotropic
function. The conditions of polyconvexity for constructing behavior. Sommer et al. (2008) performed peeling tests on
transversely isotropic or orthotropic free energy functions medial layers of abdominal aortas to investigate phenom-
are shown in Balzani et al. (2006). These material models ena of artery dissection during surgical treatment. Given the
for soft biological tissues represent a class of phenomeno- complexity of material composition and associated deforma-
logical macroscopic constitutive theories. Another class of tion mechanisms at hand, the numerical simulations of soft
models, presumably more physical, takes into account the biological tissues in the post-critical range due to fracture
entropic elasticity of collagen fibers arising due to their are rare. Gasser and Holzapfel (2006) performed numeri-
nanoscale dimensions. These models are based on the well- cal simulations of arterial dissection combining transversely
known theory of network models for rubber elasticity (see isotropic energy density with extended finite element method
James and Guth (1943), Flory and Rehner (1943), Arruda (Belytschko and Black 1999; Moës et al. 1999) where frac-
and Boyce (1993), Miehe et al. (2004)) where a randomly ture is modeled by incorporating strong discontinuities via
connected network of polymer chains is found cross-linked nodal enrichment strategies. Ferrara and Pandolfi (2008,
at the microstructure. The micromechanically motivated net- 2010) used the cohesive damage model (Xu and Needle-
work models for soft biological tissues treat the collagen man 1994) by incorporating cohesive interface elements
fibers as wormlike chains from the statistical point of view; along the edges of bulk finite elements to model tissue
see Flory (1969) for a detailed description. Some notable failure due to balloon angioplasty and medial dissection
applications of network models to soft matter such as bio- under supra-physiological loading conditions. Balzani et al.
logical tissues and fabrics can be found in Bischoff et al. (2012) resorted to continuum damage mechanics formulation
(2002), Kuhl et al. (2005), Alastrué et al. (2009) and Raina (Govindjee and Simo 1991) by introducing internal damage
and Linder (2014). variables along each fiber direction to simulate angioplasty-
The mechanical behavior of soft biological tissues under induced rupture.
normal physiological loading conditions can be qualitatively In this paper, our recent advances in the phase-field model-
simulated with the aforementioned modeling approaches. ing of fracture (see, e.g., Miehe et al. (2010a, b, 2014), Miehe
The successful applications of the same, particularly with and Schänzel (2014)) are used to simulate the complex frac-
respect to human arterial walls, have been demonstrated ture phenomena in soft biological tissues described above.
in the references mentioned above. However, under cer- The usual shortcomings of the classical Griffith-type theory
tain pathological conditions, the soft biological tissues may of brittle fracture, such as inability to determine branching
undergo fracture which is beyond the scope of modeling cracks, can be overcome by variational methods based on
capability of those methods. One such disease is atheroscle- energy minimization as suggested by Francfort and Marigo
rosis which affects the flow of oxygen-carrying blood through (1998); see also Bourdin et al. (2008), Dal Maso and Toader
the arteries due to narrowing of the lumen (stenosis). This (2002), Buliga (1999). The regularized setting of their frame-
is caused by accumulation of the extracellular lipids, cal- work considered in Bourdin et al. (2000, 2008) is obtained
cium, foam cells and necrotic tissue forming a plaque-like by convergence inspired by the work of image segmenta-
substance inside the artery wall as shown in Lendon et al. tion by Mumford and Shah (1989). We refer to Ambrosio
(1991) and Loree et al. (1992). When undiagnosed, it can and Tortorelli (1990) and the reviews of Dal Maso (1993)
lead to myocardial infarction, angina or smoker’s leg due to and Braides (1998, 2002) for details on convergent approx-
plaque rupture which releases thrombogenic materials into imations of free discontinuity problems. The approximation
blood stream (Gertz and Roberts 1990; Loree et al. 1991). regularizes a sharp crack surface topology in the solid by
Upon diagnosis, the treatment of advanced atherosclerosis diffusive crack zones governed by a scalar auxiliary vari-
is generally followed by a stent implantation in the affected able. This variable can be considered as a phase field that
artery with the balloon angioplasty procedure. The mecha- interpolates between the unbroken and the broken states
nism of balloon angioplasty requires the concerned artery of the material. Conceptually similar are recently outlined
to be subjected to supra-physiological loading which can phase-field approaches to brittle fracture based on the clas-
possibly lead to tissue degeneration and fracture (Castaneda- sical Ginzburg–Landau-type evolution equation as reviewed
Zuniga et al. 1980; Gasser and Holzapfel 2007). One refers in Hakim and Karma (2009). These models may be con-
to the reviews by Rhodin (1980) and Silver et al. (1989) sidered as time-dependent viscous regularizations of the
for a detailed description of the histological and biomechan- above-mentioned rate-independent theories of energy min-
ical features of an artery. The experimental investigations imization. The phase-field approach is embedded into the
in Schulze-bauer et al. (2002) and Holzapfel et al. (2004) theory of gradient-type materials with a characteristic length
tested layer-specific mechanical properties of the human scale, as outlined in Capriz (1989), Mariano (2001), Frémond
stenotic iliac arteries up to rupture to identify their strength (2002) and Miehe (2011) and conceptually follows the mod-

123
A phase-field model for fracture in biological tissues 481

els of continuum damage mechanics; see, e.g., Kachanov Bt = ϕt (B0 )

(1986) or Frémond and Nedjar (1996). It is to be empha- ϕ
sized here that finite element modeling of fracture such as d
X x
cohesive zone models (Xu and Needleman (1994), Ortiz
and Pandolfi (1999)), configurational-force-driven adaptive
remeshing techniques (Ortiz and Quigley (1991), Gürses B0
(a) Γ
and Miehe (2009)), element enrichment techniques (Simo
et al. (1993), Oliver (1996), Linder and Armero (2007), Bt = ϕt (B0 )
ϕ
Linder and Raina (2013)) or nodal enrichment techniques
d
(Belytschko and Black (1999), Moës et al. (1999) Wells and X x
Sluys (2001)) suffer in three-dimensional applications with
complex crack topologies. In contrast, the phase-field-type
diffusive crack modeling is a spatially smooth continuum B0
(b) Γc ϕt (Γc )
formulation that avoids the modeling of discontinuities and
can be implemented in a straightforward manner as shown reference configuration current configuration
in Miehe et al. (2010a, b, 2014) to model complex crack pat-
Fig. 1 Finite deformation of a solid with a regularized crack. The
terns. deformation map ϕ maps at time t ∈ T the reference configuration
In Sect. 2, a descriptive motivation of a regularized crack B0 ∈ Rδ onto the current configuration Bt . a The crack phase field
topology based on a phase field is presented in a purely d ∈ [0, 1) defines a regularized crack surface functional Γl (d) that con-
geometric setting. As a result, a crack surface functional verges in the limit l → 0 to the sharp crack surface Γ . b The level set
Γc = {X | d = c} defines for a constant c ≈ 1 the crack faces in the
is defined with a constitutive law that governs the evo- regularized setting. Parts of the continuum with d > c are considered
lution of the crack phase field. This crack functional is to be free space and are not displayed
considered as the crack surface itself which should stay
constant or grow for arbitrary loading processes. For numer-
ical simulations of interest, the crack surface potential is unbroken and for d(X, t) = 1 the fully broken state of the
discretized with finite elements with an operator splitting material at X ∈ B0 and time t ∈ T . It governs the regularized
strategy adopted for the evolution of the phase field. In crack surface
Sect. 3, the local crack driving force, with regard to brittle
failure in isotropic and anisotropic materials is introduced.   
For anisotropy, a general orthotropic failure criterion is 1 2 l
Γl (d) = d + |∇d| d V
2
(1)
developed to model fracture in soft biological tissues. The B0 2l 2
anisotropic contributions to the free energy with comments
on its convexity are given in Sect. 4. Finally, Sect. 5 out- that is formulated in terms of the crack surface density func-
lines representative numerical examples of fracture in soft tion γl (d, ∇d) = d 2 /2l + l|∇d|2 /2 per unit volume of the
biological tissues which demonstrate the features and algo- solid. It is governed by the length-scale parameter l. Assum-
rithmic robustness of the proposed phase-field models of ing a given sharp crack topology by prescribing the Dirichlet
fracture. condition d = 1 on Γ ⊂ B0 , the regularized crack phase field
d in the full domain B0 is obtained by a minimization princi-
ple of diffusive crack topology d = Arg{inf d Γl (d)}, yielding
2 The geometric approach to phase-field fracture the second-order Euler equation d − l 2 d = 0 in B0 with
an exponential-function-type profile of the phase field d as
We summarize the basic ingredients of a purely geometric visualized in Fig. 2b. We refer to Miehe et al. (2010b) for
approach to the phase-field modeling of fracture, which was more details.
outlined in detail in Miehe et al. (2014).

2.1 Regularization of sharp crack topology 2.2 Evolution of the regularized crack surface

Let B0 ⊂ Rδ be the reference configuration of a material Based on the global constitutive equation for the evolution of
body with dimension δ ∈ [2, 3] in space and ∂B0 ⊂ Rδ−1 the regularized crack surface functional postulated in Miehe
its surface as depicted in Fig. 1. We consider the crack phase et al. (2014), a local evolution of the phase field is given
field d : B0 × T → [0, 1] characterizing for d(X, t) = 0 the as

123
482 A. Raina, C. Miehe

d(x) lution of the regularized crack surface. This is realized by

1
expressing the local crack driving force H by the maximum

value of the associated crack driving state function D

H(X, t) = max D(state(X, s)) ≥ 0 , (4)

(a) s∈[0,t]

d(x)
1 obtained in the full process history s ∈ [0, t]. This func-
tion must be a monotonous function that depends on state
variables state of the mechanical bulk response, satisfying

(b) x
l l
unbroken
D| state

broken
= 0 and D| state = ∞ . (5)

Fig. 2 Sharp and diffusive crack modeling. a Sharp crack at x = 0. b

Diffusive crack at x = 0 modeled with the length scale l. Regularized
Note that the above ansatz is consistent with the local evolu-
curves tion equation
 obtained from minimization principle of diffusive crack topol-
ogy B0 γl d V → Min! with crack surface density function γl . Thick
line crack surface density function γl (1) with regularization profile
1 

− lδd γl ≥ 0
ex p[−|x|/l] satisfying d(0) = 1, dotted line higher order regulariza- ḋ = (1 − d) D (6)
tion approach as suggested in Borden et al. (2014) η

where
x := (x + |x|)/2 is the Macauley bracket.1 Hence,
ηḋ = (1 − d)H − [ d − l 2 d ] (2)
     a non-smooth evolution of the crack phase field takes place
evolution driving force geometric resistance when the driving force exceeds the geometric crack resis-
tance δd γl . For the rate-independent limit η → 0, the
in the domain B0 , along with the homogeneous Neumann associated local evolution equation is
condition ∇d · n0 = 0 on ∂B0 for the crack phase field. The
˙ = d/dt (·) refers to the time derivative of the
notation (·)
1 Generalized Ginzburg–Landau Equation. Equation (6) was
underlying field. Here, (1 − d)H is the local crack driving
defined in Miehe et al. (2010b) as a generalized Ginzburg–Landau equa-
force and is further assumed to be governed by the constitu- tion for the phase-field evolution, when it is represented
tive function

1
tot
ḋ = − δd ψ (7)

H = H(state(X, s), 0 ≤ s ≤ t)) , (3) η

in terms of the variational derivative of a total energy density function

where state stands for additional variables determined by tot by the phase field d. In contrast to classical phase-field equations
ψ
the model for the anisotropic bulk response under considera- as reviewed in Gurtin (1996), the above equation accounts for the irre-
tion, such as the local energetic state or the local stress state versibility due to the Macauley bracket. The total energy must contain
contributions due to the degrading bulk response and the regularized
in the solid. The left-hand side ηḋ of the local evolution in fracture surface energy. A simple example for finite elasticity at frac-
(2) represents the viscous crack resistance, where η ≥ 0 is a ture is
material parameter. The geometric crack resistance is related
to the variational derivative of the crack surface density func- tot (F, d) = (1 − d)2 ψ
ψ
(F) + gc γl (d, ∇d) , (8)
tion, i.e., d − l 2 d = lδd γl (d, ∇d). The Eq. (2) makes the
evolution of the crack phase field dependent on the difference where the parameter gc is a critical energy release rate. For this ansatz,
between the crack driving force and the geometric crack resis- (7) results in the evolution Eq. (6) with the crack driving state function
tance. Note carefully, that the above framework also covers

the rate-independent limit for η = 0, where the crack surface

= 2ψ (F)
D (9)
gc /l
is simply defined by an equilibrium for crack driving force
and geometric crack resistance.
(F) of the undamaged material. In
that contains the nominal energy ψ
this variationally consistent setting, the crack driving state function is
2.3 Irreversibility constraint on phase-field evolution derived from a total energy, see Miehe et al. (2014) for more advanced
definitions. In contrast, when starting from (2) and using (4), the defi-
nition of the crack driving is state function D
is not constrained to be
Within this work, we focus on an irreversibility of the crack related to a variational derivative of an energy. This allows the incorpo-
evolution, governed by the constraint Γ˙l (d) ≥ 0 on the evo- ration of the convenient stress-based function outlined in (19) below.

123
A phase-field model for fracture in biological tissues 483

ḋ ≥ 0, [(1 − d) D
− lδd γl ] ≤ 0, An incremental quadratic potential is introduced in Miehe
(10)

− lδd γl ] = 0
ḋ [(1 − d) D et al. (2014) as

Πdτ (d) = πdτ (cd )d V = 0 (11)
2.4 Time integration and operator split B0

The integration of (2) in the time interval [tn , tn+1 ] with incre- in terms of the quadratic potential density per unit of the
ment τ = tn+1 − tn > 0 gives an update for the regularized reference volume
crack surface. In what follows, all variables without a sub-
η
script are time discrete values at the current solution time πdτ (cd ) = (d − dn )2 + lγl (d, ∇d) − lγl (dn , ∇dn )
tn+1 . A robust algorithm for the update of the phase field is 2τ
obtained by an operator splitting, keeping the driving force 1
− [ (d − dn ) − (d 2 − dn2 ) ]H (12)
H = Hn constant in the time step [tn , t]. Note that this driving 2
force provides the impact from the bulk response to the crack
evolution. that depends on current phase field and its gradient cd :=
The above assumption covers for the one-path algorithm, a {d, ∇d}. The potential Πdτ (d) must be stationary, attaining
semi-implicit time integration of the phase field, for the multi- the minimum value zero for the current phase field d. Hence,
path algorithm, where H is corrected, an implicit update the phase field is determined by the minimization principle
scheme with Gauss–Seidel-type iterations between the phase
field d(X, t) and the state variables state(X, t) of the bulk d = Arg{inf Π τd (d)}. (13)
response. In both cases, H = Hn induces an algorith- d

mic decoupling of updates for the phase field and the bulk
response in the time interval under consideration, and is the The corresponding Euler equation of this minimization prin-
key ingredient of a modular implementation of phase-field ciple is the linear update equation for the crack phase field
fracture. The algorithm in Table 1 provides a crack update
module applicable to a wide spectrum of problems. The one-
η(d − dn )/τ = (1 − d)H − [ d − l 2 d ] (14)
path algorithm for H = Hn adopted here is in line with the      
treatment of Miehe et al. (2010a). crack update driving force geometric resistance

in the domain B0 , along with the Neumann condition ∇d ·

Table 1 Evolution and update of crack surface for given force H n0 = 0 on ∂B0 . Note that this provides in the one-path setting
with H = Hn a semi-implicit integration of the continuous
1. Incremental update of crack surface. Semi-implicit time Eq. (2). The link to the bulk response is exclusively governed
integration of (2) in [tn , t] with crack driving force
by the driving force H. Keeping it constant, fully decouples
H = const. gives linear update of crack surface
the update of the crack phase field from the update of the
constitutive response of the bulk.
η(d − dn )/τ = (1 − d)H − [ d − l 2 d ]
2.5 Finite element discretization of crack phase field

2. Finite element update of crack surface. Space discretization of The finite element implementation of the potential Πdτ in (11)
phase-field state cdh := {d h , ∇d h } = B d d d and optimization
of the incremental potential is straightforward, yielding the space-discrete formulation

Πdτ h (d d ) = Bh πdτ (B d d d )d V 
0
with the potential density function Πdτ h (d d ) = πdτ (B d d d )d V = 0 (15)
πdτ = η
(d − dn )2 + lγl (d, ∇d) − lγl (dn , ∇dn ) B0h
2τ
−[ (d − dn ) − 21 (d 2 − dn2 ) ]H
results in linear update of nodal degrees of the phase field in based on the discretization cdh (X, t) = B d (X)d d (t) in terms
[tn , t]
of the finite element interpolation matrix B d and the nodal
values d d of the phase field in a typical finite element mesh.
The nodal values follow from the optimization problem
d d = −[ Πdτ h ]−1 τh
,d d d d [ Πd ],d d

d d = Arg{ inf Π τd h (d d ) }, (16)

dd

123
484 A. Raina, C. Miehe

for the quadratic space-discrete potential, yielding the closed- driving force function. Note that this criterion is characterized
form update of the nodal values for the phase field by a threshold, that guarantees the existence of non-damaged
zones where the crack phase field is zero. To model the fail-
 −1
ure response in anisotropic biological tissues, the criterion
dd = − B dT [∂c2d cd πdτ (B d d d n )B d d V
Bh (19) needs to account for the inherent tissue morphology,
 giving rise to an orthotropic behavior, for accurate predic-
B dT [∂cd πdτ (B d d d n )] d V . (17) tion of the failure response as presented in the following
Bh
sections.
Reader is referred to Miehe et al. (2010b) for more details
on the finite element implementation of the phase-field frac-
3.2 Alignment of principal orthotropy axes in soft
ture. Note again that this solution is independent from the
tissues
finite element updates of bulk response due to the algorith-
mic decoupling within the time step under consideration. The
The histology of biological tissues makes them ideally an
only interface to the bulk response is provided by the driving
orthotropic material due to the presence of two identical fiber
force Hn that enters the potential πdτ in (12); see Table 1.
families with different directions a(1) and a(2) inclined at a
certain angle as discussed in detail in Sect. 4.2. To incorporate
an orthotropic failure criterion, local orthotropy axes g i are
3 Driving force for brittle failure and transition defined in addition to the global Cartesian coordinate axes ei
rules for i = 1, . . . , 3 as follows,
The impact from the bulk response on the crack propagation
is governed by the crack driving force H that is linked to a(1) + a(2) a(2) − a(1)
g 1 =  (1) , g 2 =  (2)  and

a constitutive crack driving state function D(state). Vice  a + a(2)   a − a(1) 

versa, the crack phase field enters the state functions of the g × g2
g3 =  1 , (20)
bulk material via constitutive transition rules 
f (state; d) f-  g1 × g2
rom the unbroken to the broken state. This section outlines
general structures of these functions, which are specified for
which are orthonormal with unit magnitude. See Fig. 3 for an
a model problem in Sect. 5.
illustration of the different coordinate systems. Based on the
local orthotropic triad {g i }i=1,3 , a second-order anisotropy
3.1 The maximum principal stress driving force
tensor A is introduced as
criterion

The key aspect is the modeling of the constitutive function H, A = a1 g 1 ⊗ g 1 + a2 g 2 ⊗ g 2 + a3 g 3 ⊗ g 3 (21)

introduced in (4), that governs the local nominal crack driving
force depending on the stress state. A brief summary of the
principal stress-based growth function developed in Miehe e2
et al. (2014), that depends on a spectral decomposition of the g2
g1
nominal Cauchy stress of the undamaged solid material is
a(1)
given below
θ

3
plane
of
ies a(2)
σ = σa na ⊗ na , (18) famil
fiber

a=1
e3 , g 3 e1
where {σa } are the principal stresses and {na } the eigenvec- 2γ
tors. The crack driving state function considered is
 3  

σa  2

iso = ζ
D −1 (19)
σcrit Fig. 3 Illustration of the global Cartesian coordinate system spanned
a=1
by the orthonormal unit vectors {ei }i=1,3 with a local orthotropic axes
spanned by the orthogonal unit vectors {g i }i=1,3 . The unit tangents a(1)
and models an isotropic failure surface in the principal stress and a(2) to the two fiber families are inclined at an angle of 2γ to each
state. σcrit > 0 is a critical fracture tensile stress and ζ > 0 other. The 1and2 directions span the plane of fiber families (plane of
a material parameter that enforces the slope of the quadratic paper) and 3 direction points into the plane of the paper

123
A phase-field model for fracture in biological tissues 485

where the coefficients {ai }i=1,3 are the scaling factors along 1
Φ= 2
(σ + : I A : σ + ) with
the orthotropy axes {g i }i=1,3 , respectively. The tensor A in 2σcrit
(21) automatically satisfies the conditions of invariance as 1
I A = (Aik A jl + Ail A jk )ei ⊗ e j ⊗ ek ⊗ el , (25)
imposed in (24). For the plane problems of interest here, the 2
components of the anisotropy tensor A for an arbitrary coun-
terclockwise rotation by angle θ of orthogonal axes {g i }i=1,2 where σcrit is the reference uniaxial critical fracture ten-
along e3 axis are given as sile stress and Ai j are the components of the second-order
anisotropy tensor A for i, j = 1, . . . , 3 introduced in (21).
⎡ ⎤
a1 cos2 θ + a2 sin2 θ (a1 − a2 ) cos θ sin θ 0 The symmetry of fourth-order tensor I A with coefficients
Ai j = ⎣ (a1 − a2 ) cos θ sin θ a1 sin2 θ + a2 cos2 θ 0⎦ . {ai }i=1,3 > 0 ensures its positive semi-definite form and
0 0 1 thereby a convex failure surface resulting from (25). Expan-
(22) sion of the first equation in (25), simplified to the case when
orthotropy axes {g i }i=1,3 and the principal axes {ni }i=1,3
It needs to be emphasized here that the above representation coincide with the global axes {ei }i=1,3 , reduces to
recovers the case of transverse isotropy and isotropy by the  2  2  2
following:
σ1 
σ2 
σ3 
Φ= + + . (26)
σcrit /a1 σcrit /a2 σcrit /a3
– Transverse isotropy: By setting a1 = a2 .
– Isotropy: By setting a1 = a2 = a3 . Substitution of (26) in (23) yields an expression for an
orthotropic crack driving state function, similar to (19), as
The coefficients {ai }i=1,3 enter as the material parameters in  3  
 
the model. However, only two independent coefficients need

ani = ζ
σi  2
D −1 . (27)
to be prescribed for a three-dimensional problem and only σcrit /ai
i=1
one for a two-dimensional problem. The remaining coeffi-
cient is set to unity for the given direction of the critical An illustration of the failure loci is shown in Fig. 4 where
fracture tensile stress σcrit . as the corresponding surface of crack driving state function
is shown in Fig. 5. Notice that for ai = 1 for i = 1, . . . , 3,
3.3 Anisotropic tensile stress-based driving force the criterion for isotropic failure surface (19) is retrieved.
criterion The simulations presented in Sect. 5 employ the general
orthotropic failure criterion (25).
In this section, we derive a general orthotropic failure crite-
rion based on the tensile part of the symmetric Cauchy stress 3.4 Transition rules from unbroken to fully broken
tensor σ . The crack driving state function in (4), alternative to response
(19), with a threshold is obtained from a scalar-valued tensor
function Φ(σ + ) as Depending on the crack phase field d ∈ [0, 1], constitutive
functions in the bulk undergo a phase transition from the
 

ani = ζ Φ(σ + ) − 1 with

D unbroken to the fully broken state. This includes conceptu-
ally a degradation, e.g., for the mechanical stresses, but also a

3
σ+ =
σa  na ⊗ na . (23) growth, e.g., for the non-mechanical fluxes at the crack faces.
a=1 A general structure of this phase transition for a generic con-
stitutive function f for the bulk response is provided by a
The argument σ + of function Φ is characterized by a symme- transition rule
try group RO(3) ⊂ SO(3), which is a subset of the special
orthogonal group SO(3), and satisfies the following invari- f (state; d) = g s (d)

 f s (state) + g c (d)

f c (state) (28)
ance condition
that depends on the phase-field and a non-specified set
Φ(σ + ) = Φ(R σ + R T ) ∀ R ∈ RO(3). (24) (state) of constitutive state variables for the bulk response.
Here,
f s (state) and

f c (state) are effective functions in the

The function Φ is restricted to quadratic terms in σ + by unbroken and the fully broken case. The dual weight func-
introducing a symmetric fourth-order anisotropy tensor I A tions may have the simple form
(Steinmann et al. 1994) and a specific invariant form of failure
surface as g s (d) = (1 − d)m and g c (d) = d m . (29)

123
486 A. Raina, C. Miehe

σ2 σ2

σcrit /a2

a1 = 2, a2 = 1 a1 = 1, a2 = 2

Dani > 0 σcrit /a2 Dani > 0

elastic range elastic range

Dani = 0 Dani = 0
σ1 σ1
σcrit /a1 σcrit /a1

(a) (b)

Fig. 4 Failure surfaces for orthotropic tensile stress failure criterion second-order anisotropy tensor (21) for the case when orthotropy axes
(27) in principal plane stress space for different degrees of anisotropies g i are coincident with the global axes ei for i = 1, . . . , 3
with coefficients a a1 = 2, a2 = 1 and b a1 = 1, a2 = 2 used in

γ̄ per unit volume. The overall response of the body under

fracture is described by the crack phase field d and the defor-
mation field ϕ written as

B0 × T → Rδ
ϕ: (30)
Dani (X, t) → x = ϕ(X, t).

ϕ maps at time t ∈ T the materials points X ∈ B0 of refer-

ence configuration B0 ∈ Rδ to spatial points x ∈ Bt in the
0
current configuration Bt ∈ Rδ . The deformation gradient is
σcrit /a2 F = ∇ϕ(X, t) ∈ GL(δ), which is imposed to constraint
σ1
σ2 J := det [F] > 0. Furthermore, let g ∈ Sym + (δ) be the
σcrit /a1
standard metric of current configuration Bt , the pull back of
Fig. 5 The surface of orthotropic crack driving state function D
ani which yields the right Cauchy-Green tensor C = F T g F. A
is shown in the principal plane stress space by plotting Eq. (27) with schematic illustration of the two primary unknowns with the
scaling coefficients a1 = 2, a2 = 1 corresponding boundary conditions is shown in Fig. 6. For
Here, g s (d) recovers for m = 1 the classical (1 − d)-theory the soft biological tissue under consideration, the body B0
of damage and for m = 2 the variational theory of brittle is assumed to consist of two fiber families denoted by unit
fracture in elastic solids outlined in Miehe et al. (2010b). vectors a(a) for a = 1, 2 with the corresponding structural
tensors given as
4 Finite elasticity coupled with phase-field fracture
M (1) = a(1) ⊗ a(1) and M (2) = a(2) ⊗ a(2) . (31)
This section focuses on the description of bulk response
In order to construct the constitutive equations for hypere-
undergoing fracture which is coupled to the phase field with
lastic response, an objective anisotropic free energy function
a constitutive setup of degrading energy. Section 4.1 presents
is introduced which depends on the deformation gradient F
primary field variables and coupled governing equations.
and the structural tensors M (1,2) as
Section 4.2 presents a brief overview of the tissue histol-
ogy concerning human iliac artery. The free energy function
(F, M (1,2) ; d, ∇d),
Ψ =Ψ (32)
characterizing the anisotropic behavior of the tissues is pre-
sented in Sect. 4.3.
and represents the energy stored per unit reference volume.
4.1 Primary fields and governing equations The dependence on the fracture phase field d and its gradient
∇d is considered as a geometric property and hence separated
Consider the motion of a body B0 given by a nonlinear defor- by a semicolon. The degradation of the stored energy due to
mation map ϕ undergoing finite strains with given body force evolving phase field d takes the form

123
A phase-field model for fracture in biological tissues 487

P ·n0 =T̄ ∇d ·n0 =0

Γ
X ∈ B0 n0 n0 X ∈ B0

ϕt d H

l
ϕ=ϕ̄
(a) (b) (c)
Fig. 6 A schematic of primary fields in body B0 undergoing fracture. b The fracture phase field is constrained by the possible Dirichlet-type
The deformation field ϕ, the fracture phase field d and the history boundary condition d = 1 on Γ and the Neumann condition ∇d ·n0 = 0
field H are defined on the solid domain B0 . a The deformation field on the full surface ∂ B0 . c The history field H defined in Eq. (4) is based
is constrained by the Dirichlet- and Neumann-type boundary condi- on the maximum of stress-based criterion in the deformation history. It
ϕ ϕ
tions ϕ = ϕ̄ on ∂ B0 and P · n0 = T̄ on ∂ B0T with ∂ B0 = ∂ B0 ∪ ∂ B0T . drives the evolution of the fracture phase field d via Eqs. (19) or (23)

 = [g s (d) + k]Ψ
Ψ
(F, M (1,2) ) (33) Table 2 Coupled governing equations of finite elasticity with fracture

following the argument in (28). Here, Ψ

is energy stored in 1. Balance of linear momentum
undamaged material with a specific form given in Sect. 4.3 Div[ P] + γ̄ = 0
and g s (d) is a monotonically decreasing function defined 2. Irreversible phase-field evolution
in (29). At a fully broken state d = 1, the full degra- ηḋ =
(1 − d)H − lδd γl  ≥ 0
dation of the energy Ψ  in (33) is circumvented by the

for small parameter
artificial elastic rest energy density k Ψ
k ≈ 0. which updates the current nodal values of the deformation
The final sets of governing coupled equations of finite elas- field at time tn+1 .
ticity with phase-field fracture are shown in Table 2 in their
strong form. The first Piola–Kirchhoff stress tensor is intro-
4.2 Morphology of human iliac artery
duced as P = [g s (d)+k]
P where
P for undamaged material
is defined in (38). Following the finite element discretization
One of the most interesting soft biological tissues to study is
of crack phase field in Sect. 2.5, the current deformation
blood-carrying arteries due to their significance in medical
field ϕ ≈ ϕ h is approximated from the minimization of dis-
examination of relevant diseases as well as their similarity to
cretized potential density πϕτ as
other tissues. Among the family of arteries, the elastic iliac
   artery is discussed here, although other arteries such as aorta
d ϕ = Arg inf πϕτ (B ϕ d ϕ ) d V with or carotid have similar structures. A detailed discussion of
dϕ B0h the microscopic structure and deformation mechanisms of
(F h , M (1,2) ; d) − γ̄ · ϕ h such tissues can be found in Rhodin (1980) and Silver et al.
πϕτ = Ψ (34)
(1989). For conciseness, only a brief summary is presented
here. An artery can be thought of as a composite structure
given in terms of the discretized deformation gradient cϕh :=
which mainly consists of three different layers: the adventi-
F h = B ϕ d ϕ . Here, B ϕ is a global strain displacement tia, the media and the intima. These layers are made up of
matrix consisting of shape functions and its derivatives, and collagen fibers of type I, surrounded by a groundmatrix con-
d ϕ ∈ Rδ is a global array of nodal position vector. A sisting of elastin, endothelial cells, fibroblast and fibrocytes.
Newton–Raphson iteration scheme is used for the solution A schematic illustration of the arterial structure is shown in
of the associated Euler equations of nonlinear elasticity at Fig. 7. The adventitia is the outermost layer whose thick-
finite strains as ness is not well defined and depends on its topographical
 −1 site and its physiological function. It is surrounded by loose
dϕ ⇐ dϕ − B ϕT [∂c2ϕ cϕ πϕτ ]B ϕ dV connective tissue with primary collagen fibers arranged in
Bh
 a double helix structure around the circumferential direc-
B ϕT [∂cϕ πϕτ ] d V , (35) tion and contributes significantly to the artery strength. The
Bh medial layer comes next to adventitial layer and is in the mid-

123
488 A. Raina, C. Miehe

where the volumetric and deviatoric parts of the isotropic

(2
) matrix are not necessarily decoupled. The volumetric part of

a
(1
a
the isotropic matrix assures the weak compressibility of the

)
γ

a
γ material in terms of low shear modulus μ ≥ 0 and parame-
ter β = 2ν/(1 − 2ν) given in terms of the Poisson ratio ν.
c i rc u m f .

The deviatoric part represents the energy stored in the soft

isotropic matrix. The polyconvexity of the free energy con-
axial
Intima Media Adventitia tribution corresponding to isotropic matrix in (36) is well
understood.
Fig. 7 Left illustration of the geometrical characteristics of an iliac
The anisotropic part Ψ
(a) in (36) is the energy stored
artery where three concentric layers of intima, media and adventitia aniso
form the overall tubular structure. Middle illustration of a dissected in wavy collagen fiber families and has been extensively
adventitial strip for tensile tests where fiber families a(1) and a(2) are proposed in the literature. Here, we follow the potential pro-
inclined at an angle γ to the circumferential direction. Right illustration posed in Billiar and Sacks (2000) for the anisotropic part to
of a medial strip dissected radially for peeling test with a single fiber
direction a
write

 
k1

(a) (F, M (a) ) =
Ψ (a)
exp {k2 (C : M − 1)} − 1
dle of the three aforementioned layers. In the media too, the aniso
k2
collagen fibers form a double helix structure, however, with
− k1 (C : M (a) − 1). (37)
a smaller pitch as compared to the adventitia, and also con-
tributes significantly to the arterial strength. The intima is the
innermost layer and mostly consists of endothelial cells. It is The scalar product C : M (a) of the right Cauchy-Green
observed that intima thickens and stiffens with age and these tensor C = F T F and structural tensor M (a) is the fourth
pathological changes are often associated with the arterial invariant which represents the square of the stretch along
disease called atherosclerosis. It involves deposition of fatty the fiber directions a(a) for a = 1, 2. The additional mate-
substances, calcium, collagen fibers, cellular waste products rial parameters k1 > 0 and k2 > 0 are introduced in
and fibrin, generally referred to as atherosclerotic plaque. (37) where the former represents fiber stiffness and latter
This changes geometrical and mechanical properties of the a dimensionless constant. The key difference of anisotropic
arterial wall and affects blood-carrying function by partial potential in (37) to the potentials introduced in Holzapfel
blockage of the lumen. A simulation of balloon angioplasty, et al. (2000) and Gasser et al. (2006) lies in the energy
which is a surgical intervention to treat atherosclerosis, is storage in compression of collagen fibers. Here, a very
performed in Sect. 5.3 with a real arterial cross section from small compressive stiffness is made available to collagen
Holzapfel et al. (2004) to model tissue rupture. fibers for better stability of numerical simulations at certain
deformations. The convexity of the anisotropic free energy
contribution in (37) can be easily established for any general
4.3 Free energy functions for soft biological tissues case by simply evaluating its Hessian to be positive semi-
definite.
The microstructure of arterial walls is composed of wavy The construction of the overall free energy function
collagen fibers embedded in a soft isotropic matrix. These (36) automatically fulfills the requirements of polyconvex-
fiber families are predominantly identified by two unit fiber ity, whose results of surface plots and convex contours are
directions a(1) and a(2) aligned at an angle γ with respect to shown in Figs. 8 and 9, respectively. A deformation gradient
circumferential direction as illustrated in Fig. 7. An undam- F = diag(λ1 , λ2 , 1) for 0.5 ≤ λ1,2 ≤ 1.5 is used with mate-
aged hyperelastic constitutive storage function for arterial rial parameters μ = 10 kPa, ν = 0.3, k1 = 20 kPa, k2 = 2
walls is postulated with the following split and γ = 40◦ . The derivate of the storage function (36) with
respect to the deformation gradient F gives the undamaged

μ μ first Piola–Kirchhoff stress

P as

(F, M (1,2) ) =
Ψ [(det F)−β − 1] + (F : F − 3)
β 2
 
isotropic matrix

P = −μ(det F)−β F −T + μF

2 
2
 
+
(a) (F, M (a) )
Ψ (36) + 2 k1 exp {k2 (C : M (a) −1)}−1 (Fa(a) ⊗ a(a) ).
aniso
a=1 a=1
 
anisotropic fibers (38)

123
A phase-field model for fracture in biological tissues 489

given in Table 3 which are in close agreement with the

similar parameters used in Gasser et al. (2006). The finite
element computation is performed with approximately 30000
displacement-based (Q1) four-noded quadrilateral elements
with the plane strain assumption and different values for vis-
cosity η = 0, 10−6 and 10−5 kNs/mm2 with length scale
Ψ

l = 0.06 mm. The problem setup is shown on the left

of Fig. 10 where a vertical displacement u is applied at
top and bottom surfaces in 104 solution steps up to com-
plete failure.The horizontal displacement at these surfaces is
1.5
1.25 1.5
not constrained. The resulting load-displacement curves up
1 1.25 to complete failure along axial and circumferential loading
1
0.75 0.75 directions are shown in Fig. 11a, b, respectively.
λ2 0.5 0.5
λ1 A quantitative comparison of load-displacement curves
with the experimental results validates the choice of mate-
Fig. 8 For the anisotropic free energy function (36) smoothly contin- rial parameters used in Table 3. The anisotropic response
uous surface plot with a softer compressive response is shown in the
stretch space {λ1 , λ2 }
with lower stiffness along the circumferential direction is
very well captured. The crack driving force obtained by the
1.5 anisotropic failure criterion (23) results in crack initiation
from the middle right finite element and propagation toward

1.25 Table 3 Material properties for uniaxial tensile tests in adventitial strips
No. Parameter Dimension Name Value

μ
λ2

1 1. kPa Shear modulus 35.71

2. ν – Poisson ratio 0.45
3. γ Degree Angle 50
0.75 4. k1 kPa Fiber stiffness 1200
5. k2 – Exponent 3.5
6. ζ – Failure slope 1.0
0.5 7. σcrit kPa Critical stress 350
0.5 0.75 1 1.25 1.5
8. a1 – Scaling coefficient 0.7
λ1

Fig. 9 The locally convex contours of the anisotropic free energy func- (a) (b)
tion (36) in the stretch space {λ1 , λ2 } is shown u u u

5 Representative numerical simulations 0.4 mm

This section presents the modeling capabilities of the frame-

work of phase-field fracture applied to the transversel-
)
a (2
a
(1

y isotropic biological soft tissues particularly the human

)
10 mm

4 mm

iliac artery wall through representative numerical examples

a

γ
shown in Sects. 5.1, 5.2 and 5.3. γ

5.1 Uniaxial tensile tests of iliac artery adventitial strips

Strips of artery wall from adventitia layer are dissected along

the circumferential and axial directions as shown in Fig. 7 3 mm 1.2 mm
and subjected to monotonic uniaxial loading tests in the
u
experimental investigations of Holzapfel et al. (2004). In
the numerical simulations presented in this section, strips Fig. 10 Problem setup: geometry and loading conditions for a uniaxial
of dimensions 10 × 3 mm2 are used and subjected to a nomi- tension test with fiber family a(1,2) orientation separated by angle 2γ
nal strain of up to 50 %. The material parameters used are and b peeling test with single fiber family a

123
490 A. Raina, C. Miehe

axial direction γ = 50◦ circumferential direction γ = 40◦

(a) (b)
1 1
0.7
η=0 0.6
η=0
0.9 0.65 η=10−6 0.9 0.58 η=10−6
η=10−5 0.56 η=10−5
0.8 0.6 0.8
exper. 0.54 exper.
0.7 0.55 0.7
0.52

0.6 0.5 0.6 0.5

Load [N]

Load [N]
15 15.5 16 16.5 17 37 37.5 38 38.5 39

0.5 0.5

0.4 0.4

0.3 0.3

0.2 0.2

0.1 0.1

0 0
0 5 10 15 20 25 0 10 20 30 40 50
Nominal Strain [%] Nominal Strain [%]

Fig. 11 Uniaxial tension test: plots of the load at top surface of the the respective mean experimental load–strain curves from Holzapfel
specimen against the applied strain for a axial direction and b circumfer- et al. (2004) is provided
ential direction, with different values of viscosity η. A comparison with

(a) (b)

0.0 d [-] 1.0 0.0 σ2 [kPa] 300

Fig. 12 Uniaxial tensile test: Simulation results of loading in the axial direction are shown with a contours of crack phase field d and b contours
of principal Cauchy stress σ2 , at different loading stages

left as soon as the stresses reach the failure surface. The snapshots are captured at different applied nominal strains to
abrupt loss of stiffness or load-bearing capacity is shown in show the evolution of crack phase field and stresses. Starting
the numerical load-displacement curves which is in accord as a small crack, the failure immediately propagates through
with the experimental observation. The influence of viscos- the tissue with a similar behavior observed in circumferential
ity parameter can be seen from the plots in Fig. 11 where direction.
the higher viscosity slightly delays crack propagation as
expected. The resulting failure contours of crack phase field d
and maximum principal Cauchy stress σ2 are shown in Fig. 12 5.2 Peeling test of iliac artery medial strips
for axial direction on the left and right, respectively, where
the elements with d > 0.9 are blanked for visualization. The Study of arterial dissection has enormous relevance of sur-
gical interest which involves tear of delicate internal lining

123
A phase-field model for fracture in biological tissues 491

of arterial wall. In this section, the experimental investiga- quadrilateral elements with the plane strain assumption and
tions of Sommer et al. (2008) are modeled where the strips viscosity η = 10−5 kNs/mm2 with length scale l = 0.06
of medial aorta are peeled in the circumferential direction to mm. The transversely isotropic form of anisotropic stress
study their dissection properties. See Fig. 7 for an illustra- driving force criterion (23) is employed for the evolution of
tion of tissue morphology. In the numerical simulations of crack phase field as the peeling predominantly takes place in
phase-field fracture presented in this section, a medial strip a single radial plane.
of dimensions 4 × 1.2 mm2 is used with an initial verti- The corresponding material parameters used for the peel-
cal notch of size 0.4 mm at the top surface as illustrated on ing simulation are given in Table 4 which agree closely with
the right of Fig. 10. The bottom surface is fixed, whereas the parameters in Gasser and Holzapfel (2006). The load per
a horizontal displacement u = 4.5 mm is applied at the unit width of one of the arms is plotted against the applied
two arms on the top surface in 103 solution steps. Contrary displacement and shown on the right of Fig. 14 with a quanti-
to the morphology in Sect. 5.1, the medial strip along the tative comparison with experimental data from Sommer et al.
radial direction of the artery in this section possesses a sin- (2008). Before the crack starts to propagate from the tip of
gle fiber family a aligned perpendicular to the direction of notch, almost linear response is observed till load per unit
applied displacement u. The spatial discretization consists width of approximately 23 mN/mm, which agrees well with
of approximately 7000 displacement-based (Q1) four-noded the experimental results of Sommer et al. (2008) and numer-
ical results from Gasser and Holzapfel (2006) and Ferrara
and Pandolfi (2010). At this nearly constant load, the crack
Table 4 Material properties for peeling test of medial strips continues to develop downwards vertically till the end of
No. Parameter Dimension Value simulation. The contours of the crack phase field and the
maximum principal Cauchy stress from peeling simulation
1. μ kPa 16.2
are shown in Fig. 13a, b, respectively, in the deformed con-
2. ν – 0.4
figuration at different stages of loading.
3. k1 kPa 45.1
4. k2 – 2.1
5. ζ – 10.0 5.3 Inflation test of atherosclerotic iliac artery
6. σcrit kPa 100
7. a1 – 1.0 The final example of phase-field modeling of fracture in soft
biological tissues is presented in this section where the pro-

0 d [-] 1.0

(a)

0 σ2 [kPa] 80

(b)
Fig. 13 Peeling test: contours of a crack phase field d and b maximum principal Cauchy stress σ2 , shown at different dissection lengths in the
deformed configuration

123
492 A. Raina, C. Miehe

30 and fibrotic intima and diseased fibrotic media are treated as

simulation single layer, i.e., fibrotic media m 3 , as followed in Balzani
experiment
25 et al. (2012). The overall solution process for the inflation test
is divided into two stages. First, a thermal boundary value
Load/width [mN/mm]

20
problem is solved in Sect. 5.3.1 to interpolate fiber orienta-
tions in fibrous layers m 1 , . . . , m 4 . Based on the computed
orientations, a mechanical boundary value problem of angio-
15
plasty is solved next in Sect. 5.3.2 to allow for the realistic
development of phase-field fracture.
10

5.3.1 Thermal boundary value problem

5

The complexity of the morphology of the atherosclerotic

0 artery requires some steps of preprocessing where the fiber
0 1 2 3 4 5
Displacement [mm]
orientations inside the material layers m 1 ,. . ., m 4 need to be
computed from given boundary conditions. In order to do so,
Fig. 14 Peeling test: plot of load per unit width on one side of the notch a heat-conduction-like problem is solved with the Dirichlet
at top surface against the applied displacement with a comparison to boundary conditions as shown in the following steps.
mean experimental value from Sommer et al. (2008)

1. The geometry in Fig. 15 is first discretized with 6086

m4 m3 four-noded quadrilateral (Q1) elements with materials
tags m 1 , . . . , m 6 corresponding to each layer.
2. For all nodes i = 1, . . . , n a lying on the closed
m5 boundaries ∂1 B, ∂2 B and ∂3 B, compute the x- and y-
components of tangent vectors with the line-weighted
m1 averaging algorithm, which is similar to area-weighted
averaging in Max (1999), as
m2
p i−1 X i + i X i−1
θ̄ Xi = and
i−1 + i
i−1 Yi + i Yi−1
θ̄Yi = . (39)
i−1 + i
m6
A schematic of the result is shown in Fig. 16a. In Eq. (39),
(·)k = (·)k − (·)k+1 are the x- and y-components of a
Fig. 15 Inflation test: cross-sectional view of human atherosclerotic vector between nodes k and k +1 normalized by its length
artery with six different material layers, namely adventitia (m 1 ), media k . For further calculations, unit tangents are adopted by
(m 2 ), fibrotic media (m 3 ), fibrous cap (m 4 ), lipid pool (m 5 ) and calci-
fication (m 6 ), as given in Holzapfel et al. (2004) normalization.
3. For the rest of n b nodes inside body B with material
tags m 1 , . . . , m 4 , compute the tangents by solving the
cedure of balloon angioplasty for an atherosclerotic artery corresponding thermal boundary value problem as
is simulated. As discussed in Sect. 1, angioplasty is a surgi-
cal intervention to mechanically widen the stenotic lumen by div(K ∇θ ) + Q = ρcθ̇ with θ = θ̄ on ∂B. (40)
inflating and implanting a stent inside the artery wall. In this
section, a two-dimensional cross section of a real atheroscle- Steady-state heat conduction (θ̇ = 0) is assumed with
rotic artery, obtained by high-resolution magnetic resonance zero body heat (Q = 0). Here, ρ is the mass density
imaging, from Holzapfel et al. (2004) is taken as the starting and c is the specific heat. For an isotropic constant ther-
point for numerical analysis, which differs from a simplified mal conductivity K = K 1, the above thermal boundary
symmetric tubular structure as shown on the left of Fig. 7. value problem (40) reduces to the homogeneous Lapla-
The actual cross section of a diseased artery consists of var- cian θ = 0 with Dirichlet boundary conditions θ = θ̄
ious layers, such as adventitia m 1 , media m 2 , fibrotic media on ∂B and is solved for x- and y-components separately.
m 3 , fibrous cap m 4 , lipid pool m 5 and calcifications m 6 , as 4. Project the nodal quantities θ A for A = 1, . . . , n a + n b
marked in Fig. 15. A layer of non-diseased intima is ignored, of Step 2 and Step 3 to the integration point I of those

123
A phase-field model for fracture in biological tissues 493

∂1 B

∂2 B

∂3 B

(a) (b) (c)

Fig. 16 Inflation test: a boundary ∂ B = ∂1 B ∪ ∂2 B ∪ ∂3 B of fibrous representing fiber directions, c zoom of expected zone of failure with
materials used in thermal problem where tangents (39) are applied, b finer mesh/fiber density from encircled part in (b)
solution of thermal problem (40) in the form of tangents at all nodes

finite element with material tags m 1 , . . . , m 4 using the Table 5 Material properties for inflation test of atherosclerotic artery
isoparametric projection as Layer μ (kPa) ν (−) k1 (kPa) k2 (−) σcrit (kPa)

a +n b m1 4.1 0.45 67 2.8 350

n
m2 10.6 0.45 112 2.8 202
θ X,Y
I
= N A θ X,Y
A
, (41)
A=1
m3 20.1 0.45 136 2.8 254
m4 22.71 0.45 175 2.8 155
where NA are the corresponding shape functions, same
as used to define the geometry. Each pair {θ XI , θYI } are
the components of a unit tangent a defining the fiber as harder isotropic neo-Hookean solids with shear modulus
orientation at each quadrature point of the discretized μ = 2 MPa and Poisson ratio ν = 0.3. Linear line ele-
body B. The resulting distribution of fiber orientations in ments are used on the boundary ∂2 B to apply a hydrostatic
fibrous layers m 1 , . . . , m 4 arranged circumferentially is pressure load of p = 15 kPa ≈ 112.5 mmHg in 104 solu-
shown in Fig. 16b. tion steps, which falls in the range of average physiological
pressure. The hydrostatic pressure load is non-conservative
The steps (1)–(4) complete the preprocessing step of com- which stays normal to the deformed surface, thereby, result-
puting fiber orientation distribution which is followed by the ing into an unsymmetrical tangent moduli. To account for
solution of mechanical boundary value problem as shown in the effect of surrounding tissues and fluids, the nodes on the
the next subsection. boundary ∂1 B of the arterial wall are constrained by linear
springs with Young’s modulus E = 0.1 kPa pointing out-
5.3.2 Mechanical boundary value problem wards.
In the two-dimensional numerical analysis presented in
We start with the same finite element discretization as used in this section, the effect of circumferential residual stresses
the previous subsection. The layers of adventitia m 1 , media (Vaishnav and Vossoughi 1987) is neglected due to its sig-
m 2 , fibrotic media m 3 and fibrous cap m 4 are treated as nificantly lower values (Balzani et al. 2012). During the
incompressible transversely isotropic fibrous layers where simulation, it is observed that peak stresses develop on the
the free energy (36) applies, with the corresponding material outer wall of the fibrous cap m 4 where lipid pool m 5 lies
parameters given in Table 5, where the critical stress val- in contact. The snapshots in Fig. 17 depict the contours of
ues are chosen from Holzapfel et al. (2004). Only the fiber crack phase field at different stages of internal pressures.
direction a(1) := a, aligned circumferentially, is applica- Depending upon the transversely isotropic form of driving
ble in (36), which is obtained from Sect. 5.3.1. The lipid force (23), one can clearly observe the initiation and propa-
pool m 5 is treated as a soft incompressible isotropic neo- gation of fracture from outer wall of fibrous cap toward the
Hookean material with shear modulus μ = 0.1 kPa and inner wall, at its thinnest cross section. The lack of availability
Poisson ratio ν = 0.4. The calcifications m 6 are treated of experimental data presents these simulation results merely

123
494 A. Raina, C. Miehe

(a) (b)

(c) (d)

Fig. 17 Inflation test: contours of crack phase field d are shown at different inflation pressures (zoom in the inset) with the maximum damage
obtained in fibrous cap at p = 15 kPa shown in (d). Elements with d > 0.9 are blanked out for visualization

as a qualitative measure of actual fracture during angioplasty, crack surface functional, governed by a crack phase field,
which can be compared with the simulation results of Fer- which converges to sharp crack topology for vanishing
rara and Pandolfi (2010) and Balzani et al. (2012). It is to length-scale parameter. A rate-dependent incremental vari-
be emphasized here that the analysis presented here is only ational framework for the evolution of crack phase field is
a mechanical simulation of the actual balloon angioplasty employed, which is driven by the maximum stress-based
procedure and ignores the associated biochemical stimuli, driving force from the deformation history. In addition, a
mass transport phenomena and three-dimensional geometri- general orthotropic failure criterion is developed and used
cal constraints. to compute the driving forces for phase-field growth in soft
biological tissues. The construction of constitutive relations
allows the degradation of full anisotropic energy storage
6 Conclusion function, which takes into account the contribution of col-
lagen fibers embedded in a soft surrounding matrix, with
We presented an application of recently developed ther- evolving crack phase field. Representative numerical simu-
modynamically consistent continuum phase-field models lations of fracture in soft biological tissues are run which
for fracture to soft biological tissues with specific focus provide quantitative comparisons with experimental data to
on human iliac arteries. One of the key characteristics of confirm the accuracy and applicability of phase-field models
the phase-field model is the derivation of the regularized of fracture.

123
A phase-field model for fracture in biological tissues 495

Acknowledgments Support for this research was provided by the Ferrara A, Pandolfi A (2010) A numerical study of arterial media dis-
German Research Foundation (DFG) for the Cluster of Excellence Exc section processes. Int J Fract 166:21–33
310 Simulation Technology at the University of Stuttgart. Flory PJ (1969) Statistical mechanics of chain molecules. Wiley,
Chichester-New York
Flory PJ, Rehner J (1943) Statistical mechanics of cross-linked polymer
networks: I. rubberlike elasticity. J Chem Phys 11:512–520
References Francfort GA, Marigo JJ (1998) Revisiting brittle fracture as an energy
minimization problem. J Mech Phys Solids 46:1319–1342
Alastrué V, Martínez VA, Doblaré M, Menzel A (2009) Anisotropic Frémond M (2002) Non-smooth thermomechanics. Springer, Berlin
micro-sphere-based finite elasticity applied to blood vessel mod- Frémond M, Nedjar B (1996) Damage, gradient of damage and principle
elling. J Mech Phys Solids 57:178–203 of virtual power. Int J Solids Struct 92:178–192
Ambrosio L, Tortorelli VM (1990) Approximation of functionals Fung YC, Fronek K, Patitucci P (1979) Pseudoelasticity of arteries
depending on jumps by elliptic functionals via γ -convergence. and the choice of its mathematical expression. Am J Physiol
Commun Pure Appl Math 43:999–1036 237:H620–H631
Arruda EM, Boyce MC (1993) A three-dimensional constitutive model Gasser TC, Holzapfel GA (2006) Modeling the propagation of arterial
for the large stretch behavior of rubber elastic materials. J Mech dissection. Eur J Mech A Solids 25:617–633
Phys Solids 41:389–412 Gasser TC, Holzapfel GA (2007) Finite element modeling of balloon
Balzani D, Neff P, Schröder J, Holzapfel GA (2006) A polyconvex angioplasty by considering overstretch of remnant non-diseased
framework for soft biological tissues. Adjustment to experimental tissues in lesions. Comput Mech 40:47–60
data. Int J Solids Struct 43:6052–6070 Gasser TC, Ogden RW, Holzapfel GA (2006) Hyperelastic modelling
Balzani D, Brinkhues S, Holzapfel GA (2012) Constitutive frame- of arterial layers with distributed collagen fibre orientations. J R
work for the modeling of damage in collagenous soft tissues with Soc Interface 3:15–35
application to arterial walls. Comput Methods Appl Mech Eng Gelsea K, Pöschl E, Aignera T (2003) Collagens-structure, function,
213–216:139–151 and biosynthesis. Adv Drug Deliv Rev 55:1531–1546
Belytschko T, Black T (1999) Elastic crack growth in finite elements Gertz SD, Roberts WC (1990) Hemodynamic shear force in rupture of
with minimal remeshing. Int J Numer Methods Eng 45:601–620 coronary arterial atherosclerotic plaques. Am J Physiol 66:1368–
Billiar KL, Sacks MS (2000) Biaxial mechanical properties of the native 1372
and glutaraldehyde-treated aortic valve cusp : part ii—a structural Govindjee S, Simo JC (1991) A micro-mechanically based contin-
constitutive model. J Biomech Eng 122:327–335 uum damage model for carbon black-filled rubbers incorporating
Bischoff JE, Arruda EM, Grosh K (2002) A microstructurally based mullins effect. J Mech Phys Solids 39:87–112
orthotropic hyperelastic constitutive law. J Appl Mech 69:570– Gundiah N, Ratcliffe MB, Pruitt LA (2007) Determination of strain
579 energy function for arterial elastin: experiments using histology
Borden MJ, Hughes TJR, Landis CM, Verhoosel CV (2014) A higher- and mechanical tests. J Biomech 70:586–594
order phase-field model for brittle fracture: formulation and Gürses E, Miehe C (2009) A computational framework of three
analysis within the isogeometric analysis framework. Comput dimensional configurational force driven brittle crack propagation.
Methods Appl Mech Eng 273:100–118 Comput Methods Appl Mech Eng 198:1413–1428
Bourdin B, Francfort GA, Marigo JJ (2000) Numerical experiments in Gurtin ME (1996) Generalized Ginzburg-Landau and Cahn-Hilliard
revisited brittle fracture. J Mech Phys Solids 48:797–826 equations based on a microforce balance. Phys D Nonlinear Phe-
Bourdin B, Francfort GA, Marigo JJ (2008) Special invited exposition: nom 92:178–192
the variational approach to fracture. J Elast 91:5–148 Hakim V, Karma A (2009) Laws of crack motion and phase-field models
Braides DP (1998) Approximation of free discontinuity problems. of fracture. J Mech Phys Solids 57:342–368
Springer, Berlin Holzapfel G, Sommer G, Regitnig P (2004) Anisotropic mechanical
Braides DP (2002) Γ -convergence for beginners. Oxford University properties of tissue components in human atherosclerotic plaques.
Press, New York J Biomech Eng 126:657–665
Buliga M (1999) Energy minimizing brittle crack propagation. J Elast Holzapfel GA, Gasser TC, Ogden RW (2000) A new constitutive
52:201–238 framework for arterial wall mechanics and a comparative study
Capriz G (1989) Continua with microstructure. Springer, Berlin of material models. J Elast 61:1–48
Castaneda-Zuniga WR, Formanek A, Tadavarthy M, Vlodaver Z, Humphrey JD (2002) Cardiovascular solid mechanics, cells, tissues,
Edwards J, Zollikofer C, Amplatz K (1980) The mechanism of and organs. Springer, New York
balloon angioplasty. Radiology 135:565–571 James H, Guth E (1943) Theory of elastic properties of rubber. J Chem
Chuong C, Fung Y (1984) Compressibility and constitutive equation of Phys 11:455–481
arterial wall in radial compression experiments. J Biomech 17:35– Kachanov LM (1986) Introduction to continuum damage mechanics,
40 1st edn. Springer, Netherlands
Dal Maso G (1993) An introduction to Γ -convergence. Birkhäuser, Kuhl E, Garikipati K, Arruda EM, Grosh K (2005) Remodeling of bio-
Boston logical tissue: mechanically induced reorientation of a transversely
Dal Maso G, Toader R (2002) A model for the quasistatic growth of isotropic chain network. J Mech Phys Solids 53:1552–1573
brittle fractures: existence and approximation results. Arch Ration Lanir Y (1979) A structural theory for the homogeneous biaxial stress–
Mech Anal 162:101–135 strain relationships in flat collagenous tissues. J Biomech 12:423–
Delfino A, Stergiopulos N, Moore J, Meister JJ (1997) Residual strain 436
effects on the stress field in a thick wall finite element model of Lendon CL, Davies MJ, Born GVR, Richardson PD (1991) Atheroscle-
the human carotid bifurcation. J Biomech 30:777–786 rotic plaque caps are locally weakened when macrophages density
Demiray H (1972) A note on the elasticity of soft biological tissues. J is increased. Atherosclerosis 87:87–90
Biomech 5(3):309–311 Linder C, Armero F (2007) Finite elements with embedded strong dis-
Ferrara A, Pandolfi A (2008) Numerical modelling of fracture in human continuities for the modeling of failure in solids. Int J Numer
arteries. Comput Methods Biomech Biomed Eng 11:553–567 Methods Eng 72:1391–1433

123
496 A. Raina, C. Miehe

Linder C, Raina A (2013) A strong discontinuity approach on multiple Mumford D, Shah J (1989) Optimal approximations by piecewise
levels to model solids at failure. Comput Methods Appl Mech Eng smooth functions and associated variational problems. Commun
253:558–583 Pure Appl Math 42:577–685
Loree HM, Kamm RD, Atkinson CM, Lee RT (1991) Turbulent pressure Oliver J (1996) Modelling strong discontinuities in solid mechanics via
fluctuations on surface of model vascular stenoses. Am J Physiol strain softening constitutive equations. Part1: fundamentals. Int J
261:H644–H650 Numer Methods Eng 39:3575–3600
Loree HM, Kamm RD, Stringfellow RG, Lee RT (1992) Effects of Ortiz M, Pandolfi A (1999) Finite-deformation irreversible cohesive
fibrous cap thickness on peak circumferencial stress in model elements for three-dimensional crack-propagation analysis. Int J
atherosclerotic vessels. Circ Res 71:850–858 Numer Methods Eng 44:1267–1282
Mariano PM (2001) Multifield theories in mechanics of solids. Arch Ortiz M, Quigley JJ (1991) Adaptive mesh refinement in strain local-
Appl Mech 38:1–93 ization problems. Comput Methods Appl Mech Eng 90:781–804
Max N (1999) Weights for computing vertex normals from facet nor- Raina A, Linder C (2014) A homogenization approach for nonwoven
mals. J Gr Tools 4(2):1–6 materials based on fiber undulations and reorientation. J Mech
Menzel A, Steinmann P (2001) A theoretical and computational frame- Phys Solids 65:12–34
work for anisotropic continuum damage mechanics at large strains. Rhodin JAG (1980) Architecture of the vessel wall. Compr Physiol
Int J Solids Struct 38:9505–9523 Handb Physiol Cardiovasc Syst Vasc Smooth Muscle 2:1–31
Miehe C (2011) A multi-field incremental variational framework for Schulze-bauer CAJ, Regitnig P, Holzapfel GA (2002) Mechanics of the
gradient-extended standard dissipative solids. J Mech Phys Solids human femoral adventitia including the high-pressure response.
59:898–923 Am J Physiol Heart Circ Physiol 282:2427–2440
Miehe C, Schänzel L (2014) Phase field modeling of fracture in rubbery Silver FH, Christiansen DL, Buntin CM (1989) Mechanical properties
polymers. Part I: finite elasticity coupled with brittle failure. J Mech of the aorta: a review. Crit Rev Biomed Eng 17:323–358
Phys Solids 65:93–113 Simo JC, Oliver J, Armero F (1993) An analysis of strong discontinuities
Miehe C, Göktepe S, Lulei F (2004) A micro–macro approach to rubber- induced by strain-softening in rate-independent inelastic solids.
like materials. Part I: the non-affine micro-sphere model of rubber Comput Mech 12:277–296
elasticity. J Mech Phys Solids 52:2617–2660 Sommer G, Gasser TC, Regitnig P, Auer M, Holzapfel G (2008) Dissec-
Miehe C, Hofacker M, Welschinger F (2010a) A phase field model tion properties of the human aortic media: an experimental study.
for rate-independent crack propagation: Robust algorithmic imple- J Biomech Eng 130:021007
mentation based on operator splits. Comput Methods Appl Mech Steinmann P, Miehe C, Stein E (1994) On the localization analysis
Eng 199:2765–2778 of orthotropic Hill type elastoplastic solids. J Mech Phys Solids
Miehe C, Welschinger F, Hofacker M (2010b) Thermodynamically 42:1969–1994
consistent phase-field models of fracture: variational principles Vaishnav RN, Vossoughi J (1987) Residual stress and strain in aortic
and multi-field fe implementations. Int J Numer Methods Eng segments. J Biomech 20(3):235–239
83:1273–1311 Vaishnav RN, Young JT, Patel DJ (1973) Distribution of stresses and of
Miehe C, Schänzel L, Ulmer H (2014) Phase field modeling of frac- strain-energy density through the wall thickness in a canine aortic
ture in multi-physics problems. Part I. Balance of crack surface segment. Circ Res 32:577–583
and failure criteria for brittle crack propagation in thermo-elastic Wells GN, Sluys LJ (2001) A new method for modelling cohesive cracks
solids. Comput Methods Appl Mech Eng. doi:10.1016/j.cma.2014. using finite elements. Int J Numer Methods Eng 50:2667–2682
11.016 Xu XP, Needleman A (1994) Numerical simulations of fast crack growth
Moës N, Dolbow J, Belytschko T (1999) A finite element method in brittle solids. J Mech Phys Solids 42:1397–1434
for crack growth without remeshing. Int J Numer Methods Eng
46:131–150

123