A variational minimization formulation for hydraulically

induced fracturing in elastic-plastic solids

∗
Daniel Kienle , Marc-André Keip

Institute of Applied Mechanics

Department of Civil and Environmental Engineering
University of Stuttgart, Stuttgart, Germany
arXiv:2103.14970v1 [math.NA] 27 Mar 2021

Abstract. A variational modeling framework for hydraulically induced fracturing of

elastic-plastic solids is developed in the present work. The developed variational structure
provides a global minimization problem. While fracture propagation is modeled by means
of a phase-field approach to fracture, plastic effects are taken into account by using a
Drucker–Prager-type yield-criterion function. This yield-criterion function governs the
plastic evolution of the fluid-solid mixture. Fluid storage and transport are described by a
Darcy–Biot-type formulation. Thereby the fluid storage is decomposed into a contribution
due to the elastic deformations and one due to the plastic deformations. A local return
mapping scheme is used for the update of the plastic quantities. The global minimization
structure demands a H(div)-conforming finite-element formulation. Furthermore this
is combined with an enhanced-assumed-strain formulation in order to overcome locking
phenomena arising from the plastic deformations. The robustness and capabilities of the
presented framework will be shown in a sequence of numerical examples.
Keywords: Variational principles, porous media, ductile fracture, hydraulic fracture,
phase-field modeling.

1 Introduction

One of the main applications for a model of hydraulically induced fractures can be found
in the recently utilized oil-production technique called fracking. During this process a
highly pressurized fluid is injected into a perforated bore hole. The goal is to induce frac-
tures in layers of the earth’s crust that store large amounts of oil and gas. The fractures
are created to increase the fluid permeability of the present soil or rock, therewith leading
to a higher flow of gas and oil into the bore hole where it is collected. The main criticisms
of this technique are related to the inducement of fractures that could lead to environ-
mental issues such as seismic activities or contamination of drinking water. Therefore,

∗
Corresponding author. E-mail address: kienle@mechbau.uni-stuttgart.de
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 2

in recent years various attempts have been made to better capture and understand the
underlying physical processes, for example by the development of models that can be used
in numerical simulations. The idea is that an accurate and robust model can help in fore-
casting both risks and potential of such a technique. To be successful, such a model has to
capture three distinct mechanisms: First, it has to describe the mechanical deformation
of the present porous medium (which could be soil or rock); second, it has to describe
the fluid transport throughout the intact and fractured porous medium; third, the model
needs to be capable of describing fracture initiation and propagation, for example driven
by fluid pressure.
To embed the present work into the literature, we comment on recent development
in the field. For the modeling of hydraulic fracturing, it is of substantial importance
to describe the mechanical deformation of the underlying porous media. Powerful tech-
niques in this direction are given on the one hand by model formulations within the
context of the so-called Theory of Porous Media, see de Boer [2000], Ehlers [2002] and
Bluhm and De Boer [1997]. An alternative approach that is suitable for the description
of fully saturated porous media with one fluid phase is formulated in the seminal works of
Terzaghi [1925] and Biot [1941]. Associated models reside in the area of the Biot theory
of consolidation, see Detournay and Cheng [1993] and Bear [1972] for a general overview.
When it comes to the modeling of fluid flow within a porous medium, it is often made
use of the law of Darcy [1856]. Darcy’s law provides a phenomenological approach to
the modeling of fluid flux that is driven by gradients of fluid pressure or, more precisely,
by gradients of the chemical potential. Darcy’s law is often referred to as a macroscopic
approach since it describes the fluid flow through the porous medium’s pore scale in an
averaged or homogenized sense. Fluid flow within fractures, i.e. regions without solid
content, could directly be modeled based on the Navier–Stokes-equation. Associated
simplifications can be based on the lubrication theory and lead to so-called Poiseuille-type
flow models. While the first approach is often used in combination with models related
to the Theory of Porous Media, the latter one can nicely be embedded in formulations
that relate to Biot’s theory of consolidation.
Since the present work aims at the description of fracturing processes in porous media,
we briefly refer to some related literature. Fundamental concepts for the description of
fracture were proposed by Griffith [1921] and Irwin [1958]. These works provide energy-
based fracture criteria for brittle materials and build the theoretical basis for the recently
developed phase-field models of fracturing. With regard to the latter we highlight the con-
tributions of Francfort and Marigo [1998], Bourdin et al. [2000, 2008], Amor et al. [2009],
Kuhn and Müller [2010], Miehe et al. [2010a] and Pham et al. [2011]. The phase-field
approach to fracture is extremely versatile and has been extended to the modeling of duc-
 D. Kienle & M.-A. Keip 3

tile fracture (Miehe et al. [2015a, 2016a,b,c], Ambati et al. [2015], Borden et al. [2016],
Alessi et al. [2018a], Steinke et al. [2020], Yin and Kaliske [2020] ), anisotropic fractur-
ing (Li et al. [2015], Teichtmeister et al. [2017], Storm et al. [2020] ) and fatigue fracture
(Alessi et al. [2018b], Lo et al. [2019], Schreiber et al. [2020], Carrara et al. [2020]), just
to name a few.
Next to the above mentioned extensions, the phase-field approach to fracture has
also seen pronounced activity in the field of porous media. Here we highlight the con-
tributions of Mikelic et al. [2015a,c,b], Miehe et al. [2015b], Wilson and Landis [2016],
Wu and De Lorenzis [2016], Mauthe and Miehe [2017] and Cajuhi et al. [2018] that are
embedded into Biot’s theory of consolidation. In that connection we also mention a recent
work on the numerical treatment of poro-elastic problems, in which an emphasis is put
on suitable finite-element formulations Teichtmeister et al. [2019]. Combinations of the
Theory of Porous Media with the phase-field approach to fracture have been provided by
Ehlers and Luo [2017, 2018a,b]; Heider and Markert [2017].
We note that the above mentioned models all relate to elastic deformations of the
material. However, there is experimental evidence that elastic material response alone
cannot capture all relevant effects Johnson et al. [1991]. This serves as a motivation to
the develop models that incorporate elastic-plastic deformations. Here we would like to
refer to the recent developments of Pise et al. [2019] and Aldakheel et al. [2020]. In the
given works, the plastic response is incorporated by introducing a Drucker–Prager-type
yield criterion function Drucker and Prager [1952]. We note that the works Pise et al.
[2019], Aldakheel et al. [2020] are not variational and consider a yield-criterion function
in terms of the effective stresses. The present work provides an alternative formulation
of elastic-plastic hydraulic fracturing in a rigorously variational setting motivated by the
ideas of Armero [1999]. In the latter work, an additive split of the fluid content into
an elastic and a plastic part is suggested, which leads to a constitutive fluid pressure as
a function of only elastic quantities. This concept is perfectly suitable for a variational
framework that combines the elastic-plastic deformation of the porous medium with the
fluid transport on the one hand and the fracture initiation and evolution on the other.
The present work is structured as follows. In Section 2 the unknown fields and associ-
ated kinematic relations are introduced. Thereafter, in 3, we provide the variational frame-
work and the constitutive functions. Here, we make use of a Darcy–Biot-type formulation
for the fluid transport in the porous medium and a Drucker–Prager-type yield-criterion
function for plastic flow. This then leads to the set of Euler equations that describe
the behavior of a fracturing porous-elastic-plastic solid. The numerical treatment of the
problem is based on a time-discrete incremental variational formulation combined with a
local return-mapping scheme and discussed in Section 4. We showcase the capabilities of
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 4

the presented numerical framework for hydraulically induced fracturing of porous-elastic-

plastic solids by means of some numerical examples to be discussed in Section 5. In detail
we consider a rigid-footing test as well as two fluid-injection tests. A summary and an
outlook will be provided in Section 6.

2 Independent and primary fields

In the present section we introduce the independent and primary fields of the porous-
elastic-plastic fracture model. They account for the elastic-plastic deformation of a body
B, for the fluid flux and storage in B as well as for the fracture initiation and evolution
in B. The surface of the body will in the following be denotes as ∂B.

Small-strain kinematics. The model presented in this work is formulated in the con-
text of the infinitesimal-strain theory. Hence, we consider the displacement field u(x, t)
at the material point x ∈ B and time t
(
B × T → R3
u: (1)
(x, t) 7→ u(x, t)

as independent field. Based on that, we introduce the infinitesimal strain tensor as a

primary variable. It cam be computed from the displacement gradient by

ε := 12 (∇u + ∇uT ). (2)

Plastic deformations. Since we are interested in modeling ductile response of the ma-
terial, we additively decompose the strain tensor (3) into elastic and plastic contributions

ε = εe + εp , (3)

wherein the plastic strain εp will be treated as a local internal variable. Further, to
describe local isotropic hardening, we introduce a local hardening field α formally by
(
B×T →R
α: (4)
(x, t) 7→ α(x, t).

Fluid mass and fluid flux. The initial density (mass per unit volume) of the fluid-
solid mixture is denoted by m0 . It can be computed from the densities of the fluid ρf and
 D. Kienle & M.-A. Keip 5

the solid ρs through the given porosity of the mixture ϕ as

Vpore
m0 = ρf ϕ + ρs (1 − ϕ) with ϕ = Vpore +Vsolid
. (5)

In the above equation, the porosity ϕ has been obtained from the given pore and solid
volume Vpore and Vsolid, respectively. The sum of solid and fluid volume is often denoted
as bulk volume Vbulk .
Following Biot’s approach to thermodynamically open systems, the mass balance in
its global and local representation reads (Biot [1984])
Z Z
d
m0 + m dV = − h · n dA ⇔ ṁ = − div h in B. (6)
dt B ∂B

Here, m denotes the change in the bulk’s density, which is caused by fluid flux through the
body’s surface. Thus, the vector h given on the right-hand side of (6) denotes a fluid-flux
vector. Formally, the latter two quantities can be introduced as
( (
B×T →R B × T → R3
m: and h: (7)
(x, t) 7→ m(x, t) (x, t) 7→ h(x, t).

In what follows we denote m as the change of fluid content. Analogous to the strain
tensor, we decompose m into an elastic and a plastic part

m = me + mp . (8)

Here, the plastic contribution mp describes a change in bulk density that is caused by
plastic deformations. It is associated with fluid that is remanently squeezed out or soaked
in due to plastic deformations.

Fracture phase field. As mentioned above, we will model cracks and their evolution
based on the phase-field approach to fracture. The fracture phase field d is thus formally
introduced as (
B × T → [0, 1]
d: (9)
(x, t) 7→ d(x, t).
It denotes with d = 0 an intact state and with d = 1 a broken state of the material. The
phase field is used to approximate a sharp crack interface Γ in a diffuse manner. This
results in the definition of a regularized crack surface Γl in terms of a crack-surface density
γ and a corresponding length-scale parameter l given by
Z
1 2
Γ ≈ Γl (d) := γ(d, ∇d) dV with γ(d, ∇d) := 2l
d + 2l |∇d|2, (10)
B
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 6
PSfrag replacements
σn = t̄

n
x∈B x∈B x∈B

u εp α
u = ū
displacement field plastic strain hardening

h · n = h̄ c∇d · n = 0

n n
x∈B x∈B x∈B x∈B

m mp h d
µ = µ̄ d = d¯
fluid content plastic fluid content fluid flux fracture phase-field

Figure 1: The unknown fields for porous-elastic-plastic solids at fracture. The

boundary ∂B is decomposed into Dirichlet and Neumann parts for the displacement
∂Bu ∪∂Bt , the fluid flux ∂Bh ∪∂Bµ and the fracture phase-field ∂Bd ∪∂Bk . For the frac-
ture phase-field, zero Neuman boundary conditions are assumed. Here c is a constant
depending on the model formulation.

Note that the sharp crack surface is recovered for vanishing length-scale parameter (l →
0 ⇒ Γl → Γ). Here we follow the notation of Miehe et al. [2010b].

3 Variational formulation of fracturing porous-elastic-

plastic solids

Based on the previous section, we are able to introduce the primary fields for the descrip-
tion of fracturing porous-elastic-plastic solids as

U := {u, h, d, εp , α, mp }. (11)

In order to formulate a rate-type variational principle we define the rate of the primary
fields as
U̇ := {u̇, h, d,
˙ ε̇p , α̇, ṁp }. (12)

Furthermore, we identify the constitutive state of the model and its evolution as

C := {ε, m, h, d, ∇d, εp , α, mp } and Ċ := {ε̇, h, div h, d,

˙ ∇d,
˙ ε̇p , α̇, ṁp }, (13)

respectively.
 D. Kienle & M.-A. Keip 7

3.1 Formulation of the rate-type potential

The general form of the rate-type potential is given by

d
Π(U̇; U) := E(Ċ; C) + D(Ċ) − Pext (U̇), (14)
dt
d
where dt E(Ċ; C) is the rate of the energy, D(Ċ) is the dissipation potential and Pext (U̇)
is the potential of the external loading. The rate of the energy is described in terms of
the energy density ψ(C) Z
d d
E(Ċ; C) = ψ(C) dV, (15)
dt dt B
which, by application of the chain rule yields
Z
d
E(Ċ; C) = (∂ε ψ : ε̇ − ∂m ψ div h + ∂d ψ d˙ + ∂εp ψ : ε̇p + ∂α ψ α̇ + ∂mp ψ ṁp ) dV. (16)
dt B

In the latter equation, we made use of the fluid mass balance (6)2 . Similarly to the rate
of energy the dissipation potential can be expressed in terms of a dissipation potential
density φ(Ċ) Z
D(Ċ) = φ(Ċ) dV. (17)
B

Combing the right-hand sides of equations (16) and (17) yields the internal rate-potential
density per unit volume

π(Ċ; C) = ∂ε ψ : ε̇ − ∂m ψ div h + ∂d ψ d˙ + ∂εp ψ : ε̇p + ∂mp ψ ṁp + ∂α ψ α̇ + φ(Ċ). (18)

so that Z
Π(U̇; U) := π(Ċ; C) dV − Pext (U̇). (19)
B

The particular forms of the energy density ψ(C), the dissipation-potential density φ(Ċ)
and the potential of the external loading Pext (U̇) will be discussed in the following sections.

3.1.1 Constitutive energy density

The energy density has two contributions, one from the solid ψsolid and one from the fluid
ψfluid
ψ(C) = ψsolid (ε − εp , α, d) + ψfluid (ε − εp , m − mp ). (20)
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 8

Energy density of the solid phase. The energy density of the solid phase ψsolid can
be decomposed into an effective elastic and a plastic part

ψbulk (εe , α, d) = ψeff (εe , d) + ψplast (α, d). (21)

Both parts depend on the fracture phase field d by means of a degradation function g(d).
The degraded effective elastic energy reads

ψeff (εe , d) = [g(d) + k]ψeff

0+ e 0− e
(ε ) + ψeff (ε ), (22)

where the superscript “0” indicates the energy density of the undamaged solid matrix.
In the latter equation, we have decomposed the undamaged energy into tensile and com-
pressive parts (indicated by the superscripts “+” and “-”, respectively), from which only
the tensile part is assumed to contribute to fracture propagation. The parameter k ≪ 1
ensures the well posedness of the problem. In what follows, we assume that g(d) = (1−d)2 .
0±
The undamaged effective energies ψeff represent the behavior of the elastic matrix and
take the simple quadratic forms

0± e

2 P
ψeff (ε ) = λ
2
tr (εe ) ±
+ G εe± : εe± with εe± = e
i=1,3 hλi i± ni ⊗ ni (23)

in terms of the elastic strain εe and the Lamé constants λ ≥ − 32 G and G ≥ 0. The tensile
and compressive strains εe± are given in terms of the eigenvalues of the strain tensor λei
and the ramp function hxi± = (x ± |x|)/2.
The plastic energy density considers isotropic saturation-type hardening and takes the
form

0 0
 
ψplast (α, d) = [g(d) + k]ψplast (α) with ψplast (α) = h2 α2 + σy α + ω1 exp(−ωα) − ω1 . (24)

Here, h ≥ 0 is the hardening modulus, σy is the saturated yield shift and ω is a sat-
uration parameter, see Kienle et al. [2019]. The resulting hardening function β(α, d) =
∂α ψplast (α, d) is shown in Figure 3.

Energy density of the fluid phase. Based on Biot’s theory of consolidation (Biot
[1941]) the fluid energy density is chosen as
 
me 2
ψfluid (εe , me ) = M
2
b tr(εe ) − ρf
, (25)
 D. Kienle & M.-A. Keip 9

where M is Biot’s modulus, b is Biot’s coefficient and ρf is the fluid density. It satisfies
the following conditions in terms of the fluid pressure p

p := − 1b ∂tr ε ψfluid = ρf ∂m ψfluid . (26)

We refer to Miehe et al. [2015b] for a more detailed discussion on the construction of
ψfluid . Note that according to (25) and (26) the fluid pressure depends only on the elastic
quantities εe and me .

3.1.2 Dissipation potential density

Similar to the energy density the dissipation potential density can be additively decom-
posed into individual contributions, here associated with dissipative effects arising from
fluid flow, fracture evolution and plastic deformations. The dissipation potential density
is formally given by

φ(Ċ) = φfluid (h) + φfrac(d,

˙ ∇d)
˙ + φplast (ε̇p , α̇, ṁp ). (27)

Dissipation-potential density for fluid flow. We follow Miehe et al. [2015b] and
employ a dissipation-potential density in terms of the fluid flux h given by

1
φfluid (h) = K −1 : (h ⊗ h) (28)
2

The permeability tensor K at a given state {ε, d, ∇d} is defined as

K(ε, d, ∇d) = [1 − f (d)]K 0 + f (d)K frac(ε, ∇d), (29)

where K 0 is the permeability tensor of the undamaged bulk and K frac is the permeability
tensor within a crack. Clearly, the function f (d) acts as an interpolation function between
intact and fully damaged states of the material. In what follows, we select f (d) = dǫ ,
where ǫ is an interpolation parameter.
While the permeability tensor of the undamaged bulk can be formulated in an isotropic
manner based on the spatial permeability K as K 0 = ρ2f K1 , the permeability tensor
within a crack can be expressed in terms of the fracture-opening function w(ε, ∇d)

2

K frac (ε, ∇d) = ρ2f w 12η
(ε,∇d)
f
1 −n⊗n with n := ∇d
|∇d|
, (30)

which is anisotropic in nature. In the above definition, ηf is the fluid’s dynamic viscosity.
Note that the representation of the permeability in (30) is derived based on the lubrica-
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 10

PSfrag replacements
w
L⊥

a) b)

Figure 2: Schematic representation of a) fluid flow in porous medium according to

Dracy’s law and b) within developing fractures according to Poiseuille–type law with
fracture opening w.

tion theory and relates to Poiseuille-type flow within the fractures. In this context, the
fracture-opening function is given as

w(ε, ∇d) = (n · ε · n)L⊥ , (31)

where L⊥ is the length of a line element that is perpendicular to the crack. In a finite-
element representation this can be identified as the element size he .

Dissipation-potential density for fracture evolution. The dissipation-potential

density associated with fracturing accounts for the dissipation of an evolving crack surface.
It preserves thermodynamical consistency by ensuring a local irreversibility condition of
the fracture phase field d˙ ≥ 0, which can be achieved by introducing the indicator function
˙
I(d).
For general Griffith-type fracturing the dissipation-potential density reads

dn o
φgfrac
c
(d,˙ ∇d)
˙ := ˙
gc γ(d, ∇d) + I(d), (32)
dt
˙ is given as
where gc is Griffith’s critical energy-release rate. The indicator function I(d)
(
˙ := 0 for d˙ ≥ 0,
I(d) (33)
∞ otherwise.

Note that the fracture evolution described by (32) does not include any threshold value
for the fracture evolution. In that case, damage occurs already at very small load levels.
Hence, to clearly separate the material response during loading in elastic and plastic
behavior as well as subsequent fracturing, the following dissipation-potential density will
be used n o
φfrac(d, ˙ := d [1 − g(d)]ψc + 2ψc lγ(d, ∇d) + I(d),
˙ ∇d) ˙ (34)
dt
where ψc is a threshold value.
 D. Kienle & M.-A. Keip 11

Dissipation-potential density for plastic response. The dissipation-potential den-

sity accounting for plastic behavior can be derived based on the principle of maximum
dissipation. For that, the thermodynamical driving forces for the plastic strain, the hard-
ening and the change of plastic fluid content need to be specified. Using the second law
of thermodynamics (Clausius–Planck inequality) yields the driving forces as thermody-
namical duals of εp , α and mp as

σ := −∂εp ψ = −∂εp ψeff − ∂εp ψfluid = σ eff − bρf µ1

β := −∂α ψ = −∂α ψplast (35)
µ := −∂mp ψ = −∂mp ψfluid .

We denote the set of diving forces as F = {σ, β, µ}, wherein σ is the Cauchy stress,
β is the hardening function and µ is the fluid potential. Based on that, we construct
the dissipation-potential density φplast , which is formulated in terms of the constrained
optimization problem

φplast (ε̇p , α̇, ṁp ) := sup[σ : ε˙p + β α̇ + µṁp ] . (36)

F∈E

In the above definition, the plastic driving forces are constrained to lie within the elastic
domain E := {F|f p (F) ≤ 0}, which is characterized by the yield function f p . The
latter will be specified at a later stage. By introducing the Lagrange multiplier λp the
constrained optimization in (36) can be rewritten as
 
φplast (ε̇p , α̇, ṁp ) := sup inf σ : ε˙p + β α̇ + µṁp − λp f p (F) , (37)
F λp ≥0

The latter representations are related to the rate-independent elastic-plastic material re-
sponse, which leads to non-smooth evolution of plastic deformations. The non-smoothness
can be relaxed by introducing viscous regularization, which yields the modified dissipation-
potential density
 
φplast (ε̇p , α̇, ṁp ) := sup σ : ε˙p + β α̇ + µṁp − 1
2ηp
hf p (F)i2+ . (38)
F | {z }
Dplast

From a mathematical point of view this density is obtained from an optimization proce-
dure with side condition, the latter of which is enforced by a quadratic penalty term. The
penalty parameter is given by the plastic viscosity ηp .
p
In the present work it is assumed that the yield function fsolid describes the plastic
response of the drained solid matrix. It is thus expressed in terms of the effective stresses
σ eff and the hardening function β. It should be equal to the yield function f p of the
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 12
q
2
β(α) β(α) 3 | dev[σ eff ]|
h
σy Mφ
1
elastic domain
smax 1
tr σeff
ω 3
s∗max

1 fp = 0
α∗ α
a) b)
Figure 3: Hardening function β(α; d = 0) with linear hardening h 6= 0 (solid line)
and without linear hardening h = 0 (dashed line), where α∗ = ( σhy + ω)−1 in a). In b)
the yield function in two-dimensional hydrostatic-deviatoric
q plane with and without
2
regularization (dashed/solid lines) where s∗max = smax − 3 q1 . In green with hardening
and no hardening in purple.

undrained bulk, which can be expressed in terms of the total stress σ, the hardening
function β and the fluid potential µ

p
fsolid (σ eff , β) = f p (σ, β, µ). (39)

A similar assumption has been made in Armero [1999]. The starting point for the deriva-
tion of an appropriate yield function for the presented model is given by the yield function
for frictional materials presented in Kienle et al. [2019], see also Vermeer and de Borst
[1984] and Lambrecht and Miehe [2001]. Making use of equations (35)1 and (39) yields
q q
p
f (σ, β, µ) = 32 | dev[σ]|2 + Mφ2 q12 − Mφ (smax − 13 tr σ − bρf µ) + βMh (σ, µ), (40)

where Mφ is related to the friction angle, q1 is a parameter related to the regularization

of the tip of the yield surface and smax is related to the cohesion of the material. The
hardening response is limited to friction hardening by choosing the material function
Mh (σ, µ) as q  
Mh (σ, µ) = 1 − 32 q1 exp 13 tr σ + bρf µ − smax . (41)
p
By inserting equation (35)1 in equation (40) the yield function fsolid (σ eff , β) can be recov-
ered. The yield function in terms of the effective stress σ eff is visualized in Figure 3.
 D. Kienle & M.-A. Keip 13

3.1.3 Potential of external loading

The external loading in form of mechanical tractions and fluid potential is formulated as
Z Z
Pext (U̇) = t̄ · u̇ dA − µ̄h · n dA, (42)
∂Bt ∂Bµ

where t̄ is the mechanical traction vector applied on the traction boundary ∂Bt of the
domain B. The fluid contribution is due to the fluid transport over the boundary ∂Bµ of
the domain B, where the fluid potential µ̄ is applied.

3.2 Minimization principle and mixed variational principle

Based on the above introduced functions we can introduce a rate-type minimization prin-
ciple that governs the boundary-value problems of porous-elastic-plastic solids at fracture


U̇∗ = arg inf Π(U̇; U) (43)
U̇∈W

where W = {Wu̇ , Wh , Wd˙, Wε̇p , Wα̇ , Wṁp } is the set of admissible spaces corresponding
to the set of the rate of unknowns U̇. The admissible spaces are given as

˙ on ∂Bu }, Wh := {h ∈ H(div, B)|h · n = h̄ on ∂Bh }

Wu̇ := {u̇ ∈ H 1 (B)|u̇ = ū
Wd˙ := {d˙ ∈ H 1 (B)}, Wε̇p := {εp ∈ L2 } (44)
Wα̇ := {α̇ ∈ L2 }, Wṁp := {ṁp ∈ L2 }.

Combining the global minimization principle (43) with the local maximization princi-
ple in (38) yields the mixed variational principle

 Z 
∗ ∗ ⋆
{U̇ , F } = arg inf sup π (Ċ, F; C) dV − Pext (U̇) (45)
U̇∈W F∈L2 B

where π ⋆ (Ċ, F; C) is the mixed potential density. It reads

π ⋆ (Ċ, F; C) = d
dt
ψ(C) + Dplast (ε̇p , α̇, ṁp , F) + φfluid (h) + φfrac(d,
˙ ∇d).
˙ (46)

Performing the variation of (46) at a fixed state C, we obtain the Euler equations of the
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 14

mixed variational principle (45) for the global unknowns as

div[∂ε ψ] = 0 in B
∇[∂m ψ] + ∂h φ = 0 in B
∂d ψ + ∂d˙φ − div[∂∇d
˙ φ] ∋ 0 in B
(47)
−∂m ψ + µ̄ = 0 on ∂Bµ
∂ε ψ · n − t̄ = 0 on ∂Bt
∂∇d
˙ φ·n = 0 on ∂Bk

and the local unknowns with the corresponding thermodynamic duals as

∂εp ψ + σ = 0 in B
∂α ψ + β = 0 in B
∂mp ψ + µ = 0 in B
(48)
ε̇p − λpv ∂σ f p = 0 in B
α̇ − λpv ∂β f p = 0 in B
ṁp − λpv ∂µ f p = 0 in B

The equations in (47) represent the global balance laws and the corresponding Neumann
boundary conditions. In (48) the definition of the thermodynamic duals of the local fields
as well as their evolution equations are given. In the latter, we introduced the visco-plastic
1
multiplier λpv := ηp
hf p i+ .

Condensation of local variables. The set of the rate of the primary fields can be
split into a local and global part. The former is given as U̇l = {ε̇p , α̇, ṁp } and the latter
arises as U̇g = U̇ \ U̇l = {u̇, d,˙ h}. The set of the rate of the local fields U̇l is governed
by the mixed variational principle


{U̇∗l , F∗ } = arg inf sup π ⋆ (Ċ, F; C) (49)
U̇l F

In order to obtain a solution for U̇g the reduced potential density is introduced as

⋆
πred (Ċred ; Cred ) = inf sup π ⋆ (Ċ, F; C), (50)
U̇l F
 D. Kienle & M.-A. Keip 15
PSfrag replacements
σ 0eff σ 0eff σ 0eff

0
ψplast =0 0
ψplast wplast

0
ε 0
ε 0
ε
ψeff ψeff ψeff
a) b) c)

Figure 4: Visualization of the different energy contributions to the crack driving

history field H. For the representation of H in (53)2 only the elastic energy and the
energy arising from the hardening, as show in a) and b), will drive the crack. For ideal
0
plasticity the plastic energy vanishes (ψplast = 0) leading to pure brittle fracture, see
a). By using (56) the crack is driven by the elastic energy and the full plastic work,
see c).

where the reduced constitutive state Cred = {ε, m, d, ∇d} and its evolution Ċred =
{ε̇, h, div h, d,
˙ ∇d}
˙ are introduced. The set of the rate of the global fields is given by
the following minimization principle

 Z 
U̇∗g = {u̇ , d , h } = arg
∗ ˙∗ ∗
inf inf inf ⋆
πred (Ċred ; Cred ) dV − Pext (U̇) (51)
˙
u̇∈Wu̇ d∈W
d˙
h∈Wh B

3.3 Modification of the fracture driving force

In this section a closer look at equation (47)3 is taken. Inserting the definitions of the
energy density and the dissipation potential density yields

0 0
− 2(1 − d)[ψeff + ψplast − ψc ] + 2ψc (d − l2 ∆d) + ∂d˙I ∋ 0. (52)

By introducing the crack driving history field H this is modified to

− 2(1 − d)H + 2ψc (d − l2 ∆d) = 0 with H := max hψeff

0 0
+ ψplast − ψc i+ . (53)
s∈[0,t]

This follows the notation for brittle fracture in Miehe et al. [2010a] and for ductile fracture
in Miehe et al. [2015a, 2016a].
For plasticity models with a yield limit that is independent of the stress state it is
possible to formulate a plastic energy density ψeplast
0
(α), which contains not only the work
of the hardening for the solid matrix but also the work of the ideal plastic deformation
of the solid matrix. One such example is given by von-Mises plasticity, for which we can
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 16

write Z
von Mises (49)
ψeplast
0
(α) = wplast with wplast = σ 0eff : ε̇p dt (54)

Note that the undamaged effective stress σ 0eff acting on the solid matrix is used here. The
plastic work wplast can alternatively be expressed in terms of the undamaged total stress
σ 0 and the fluid potential µ as
Z Z
wplast = σ 0eff p
: ε̇ dt = σ 0 : ε̇p + µṁp dt. (55)

Since the construction of an energy density ψeplast

0
like (54)1 is not possible for more com-
plicated plasticity models such as the Drucker–Prager model, the crack driving history
field in (53)2 is modified to

0
H = max hψeff + wplast − ψc i+ . (56)
s∈[0,t]

With the representation of the crack driving history field in (56) it is possible to model
ductile fracture evolution that is driven by the elastic and the ideal plastic deformation as
well as the hardening. With the representation in (53)2 and the definition of the plastic
energy in (24)2 it is only possible to model a ductile fracture evolution which is driven by
the elastic deformation and the hardening.

3.3.1 Relation between plastic strain and change of fluid content

Based on the construction of the yield function, see (39), the evolution of the plastic strain
(48)4 can be reformulated as

p
ε̇p = λpv ∂σ f p = λpv ∂σeff fsolid . (57)

Furthermore, reformulation of the evolution of the change of plastic fluid content (48)6
yields
p
ṁp = λpv ∂µ f p = λpv ρf b tr(∂σeff fsolid ). (58)

Combining the above two equations gives a relation between the evolution of the plastic
strain and the evolution of the change of the plastic fluid content

ṁp = ρf b tr ε̇p . (59)

We observe that the change of the plastic fluid content depends exclusively on volumetric
plastic deformation.
 D. Kienle & M.-A. Keip 17

4 Numerical Treatment

4.1 Incremental variational formulation

The incremental version of the rate-type potential introduced in Section 3 is obtained by

algorithmic time integration over a given time step τ = [tn , tn+1 ). For a pure Dirichlet
problem (Pext = 0) we arrive at
nZ tn+1 o Z
τ
Π (U, F) = Algo Π(U̇, F; C) dt = π ⋆τ (C, F; Cn ) dV, (60)
tn B

where π ⋆τ (C, F; Cn ) is the incremental potential density. It is given in terms of the energy
density ψ, the incremental fluid and fracture dissipation-potential density φτfluid and φτfrac ,
respectively, as well as the incremental dissipation density related to the visco-plastic
behavior Dτplast as

π ⋆τ (C, F; Cn ) = ψ(C) + Dτplast (εp , α, mp , F; εpn , αn , mpn ) + φτfluid (h) + φτfrac(d, ∇d, dn ). (61)

The individual dissipative contributions read

φτfluid = τ φfluid (h)

φτfrac = [1 − g(d)]ψc + 2ψc lγ(d, ∇d) + I τ (d, dn ) (62)
Dτplast = σ : (εp − εpn ) + β(α − αn ) + µ(mp − mpn ) + τ
2ηp
hf p (F)i2+ .

Furthermore the fluid mass balance (6)2 is satisfied by the implicit update

m = mn − τ div h. (63)

Condensation of local variables. Similar as in the continuous problem, the set of

primary fields can be decomposed into a local and global part. Again, the local fields are
identified as Ul = {εp , α, mp } and the global fields as Ug = U \ Ul = {u, d, h}. The local
fields are governed by the mixed variational principle

U∗l = arg{ inf sup π ⋆τ (C, F; Cn )} (64)

Ul F
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 18

Using the representation (38) the mixed variational principle (64) leads to the following
condition  
∂ ψ+σ
εp
 
 ∂α ψ + β 
 
 ∂mp ψ + µ 
 
∂Ul ,F π ⋆ =  p  = 0. (65)
 ε − εpn − γvp ∂σ f p 
 
 α − α − γp ∂ f p 
 n v β 
m − mn − γvp ∂µ f p
p p

τ
Here, γvp := τ λpv = ηp
hf p i+ is the incremental visco-plastic multiplier. The local system
of equations in (65) is solved via a general return mapping scheme summarized in Box 1.
 D. Kienle & M.-A. Keip 19

0. Get trial values [εe,tr ; me,tr ; αtr ]T = [(ε − εpn ); (m − mpn ); αn ]T

1. Set initial values [εe ; me ; α]T = [εe,tr ; me,tr ; αtr ] and γvp = 0.
2. Compute derivatives of energy density and yield function

   
ψ,εe ψ,εe εe ψ,εe me ψ,εe α
   
ψ(εe , me , α; d) s =  ψ,me  E =  ψ,me εe ψ,me me ψ,me α 
ψ,α ψ,αεe ψ,αme ψ,αα
   
f p ,σ f p ,σσ f p ,σµ f p ,σβ
 p   p 
f p (σ, µ, β) n =  f ,µ  F =  f ,µσ f p ,µµ f p ,µβ 
f p ,β f p ,βσ f p ,βµ f p ,ββ

3. Check for yielding. If yielding do a local Newton iteration

if f p < 0 then // elastic step

set [σ; ψ,m ]T = [ψ,εe ; ψ,me ]T and Eep = E
return
else // plastic step
compute residual vector
r := [(εpn − εp ); (mpn − mp ); (αn − α)]T + γvp n
check if local Newton is converged
q
if [ r T r + [ f p − ητp γvp ]2 < tol ] go to 4.

compute incremental plastic parameter

ηp
∆γvp = C1 [ f p − nT X r ] with C := nT X n + τ
and X := [ E−1 + γvp F ]−1
compute incremental strains, plastic fluid content and hardening variable
[∆εp ; ∆mp ; ∆α]T = −E−1 X [ r + ∆γvp n ]
update plastic quantities
[εp ; mp ; α; γvp ]T ⇐ [εp ; mp ; α; γvp ]T + [∆εp ; ∆mp ; ∆α; ∆γvp ]T
update elastic quantities
[εe ; me ]T = [(ε − εp ); (m − mp )]T
go to 2.

4. For plastic step: Obtain stresses and consistent moduli

1
[σ; ψ,m ]T = [ψ,εe ; ψ,me ]T and Eep = X − [X· n]⊗ [n·X]
C

Box 1: Return mapping and tangent moduli for poro-elasto-plasticity. It is based on the
algorithm for elasto-plasticity in Miehe [1998].
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 20

Reduced global problem. In order to obtain a solution for Ug the reduced potential
density is introduced
⋆τ
πred (Cred ) = inf sup π ⋆τ (C, F; Cn ), (66)
Ul F

where the reduced constitutive state Cred = {ε, h, div h, d, ∇d} is introduced. The global
fields are then given by the minimization principle

 Z 
U∗g = {u , d , h } = arg
∗ ∗ ∗
inf inf inf ⋆τ
πred (Cred ) dV , (67)
u∈Wu d∈Wd h∈Wh B

with the admissible spaces

Wu := {u ∈ H 1 (B)|u = ū on ∂Bu },
Wh := {h ∈ H(div, B)|h · n = h̄ on ∂Bh }, (68)
Wd := {d ∈ H 1 (B)}.

4.2 Space-time-discrete finite-element formulation

Considering a finite-element discretization T h (B), the discrete state vector d containing

the discrete values of {u, h, d} and the interpolation of the constitutive state Chred = Bd
the global minimization principle (67) can be written as
Z
∗ τh τh ⋆τ h
d = arg{inf Π (d)} with Π (d) = πred (Bd) dV. (69)
d B

Here, we employ the shape functions from Raviart and Thomas [1977] for the interpo-
lation of the fluid flux, see also Teichtmeister et al. [2019]. The nodal displacement
is interpolated by the shape functions of the enhanced-assumed-strain formulation, see
Simo and Rifai [1990]. The interpolation of the nodal phase-field values is done by the
standard Q1 -type shape functions.
The global algebraic minimization principle (69) leads to the following condition
Z
R := Πτ,dh = BT S dV = 0 (70)
B
 D. Kienle & M.-A. Keip 21

where the generalized array S is introduced. This array is defined as follows

 
ψ,εh
 
 −τ ψ,mh 
 
 τ φ,hh 
S := 


 (71)
 
 −2(1 − d)H + φτ ′ 
 ,dh 
τ
φ,∇dh

Note that this array is not obtained by straightforward differentiation of (69)2 , i.e. we
⋆τ h
have in general S 6= ∂Chred πred . This is due to the consideration of the history field H in
(53). Above, we have introduced the φτ ′ = φτ − I τ (d, dn ), where I τ (d, dn ) is the time
discrete version of the indicator function (33). The system of equations (70) is solved by
a Newton–Raphson-type iteration yielding
Z
−1
d ← d − K R with K := Πτ,dd
h
= BT CB dV. (72)
B

The generalized tangent array C is given as

 
 Eepεε −τ Eepεm 0 · · 
 
 −τ Emε τ Emm
ep 2 ep
0 · · 
 
 
C := ∂Chred S =  0 0 τ φ,hh hh 0 0 . (73)
 
 
 · · 0 2H + φτ,d′h dh 0 
 
 · · 0 0 τ
φ,∇dh∇dh 

Here ”0” indicates that the corresponding derivative does not exists. The derivatives
at the slots labeled with ”·” indeed exist but are not needed due to the modification
of the fracture driving force (53) as consequence of an operator split. The latter leads
to a decoupling of the related fields so that the two boxed sub-blocks in (73) can be
treated in separate solution steps. The so called one-pass solution strategy is utilized here
(Miehe et al. [2010a, 2015a]). This means that the displacement and flux is updated first
and then the fracture phase-field is updated. This might underestimate the speed of the
fracture evolution but can be controlled by the choice of an appropriated time step size τ
(Kienle et al. [2019]).
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 22

Table 1: Material parameters for rigid footing test.

Lamé parameter λ = 180.0 GN/m2 slope yield function Mφ = 0.6 –

2
Lamé parameter G = 31.0 GN/m position of peak smax = 4.0 MN/m2
hardening modulus h = 0.035 MN/m2 Biot’s modulus M = 25.0 GN/m2
saturated yield shift σy = 0.1 MN/m2 Biot’s coefficient b = 0.5 –
saturation parameter ω = 2.0 – fluid dyn. viscosity ηf = 1.0 · 10 −3
Ns/m2
plastic viscosity ηp = 5 · 10−6 s permeability K = 9.8 · 10−12 m3 s/kg
2
perturbation parameter q1 = 0.04 MN/m fluid density ρf = 1000.0 kg/m3

5 Numerical Examples

In the following we present a sequence of numerical examples that demonstrate the capa-
bilities of the model formulation. The examples start with a test of a porous-elastic-plastic
medium that is surcharged with a rigid footing leading to the creation of shear bands.
Furthermore, we analyze the effect of different driving forces on porous-elastic-plastic frac-
ture evolution. The latter analysis is extended to the comparison of porous-elastic and
porous-elastic-plastic materials response leading to the evolution of Hydraulically induced
fractures.

5.1 Rigid-footing test on porous-elastic-plastic medium

In the first example we consider a rigid footing test without fracture evolution. Our goal is
to analyze the effect of plasticity as well as fluid flux and storage on the system’s response.
We thus take into account three different material types:

i) drained elastic-plastic material,

ii) undrained porous-elastic material,
iii) undrained porous-elastic-plastic material with different permeabilities.

The drained elastic-plastic material is recovered by setting Biot’s modulus and coefficient
to zero (M = 0, b = 0). In order to model the undrained porous-elastic material, the
yield limit is increased to a very high value by setting smax = 1 · 104 MN/m2 . For
the undrained porous-elastic-plastic material all contributions of the model are active,
hence no artificial choice of any material parameter is necessary. In order to analyze the
effect of the permeability on the overall model response, we consider different magnitudes
of permeabilities given by an original value K as well as a reduced and an increased
permeability (K/5 and 5K, respectively). The chosen material parameters are listed in
Table 1.
The geometry and boundary conditions are shown in Figure 5. Due to the symmetry
of loading and geometry only one half of the specimen is discretized by 2, 376 quadri-
 D. Kienle & M.-A. Keip 23

ū ū
Sfrag replacements

a a/2
H H

W W/2

Figure 5: Rigid footing test: Geometry and boundary conditions. Bottom is mechan-
ically fixed, left and right edge is mechanically fixed in horizontal direction. Bottom,
left and right is impermeable.

9.4 · 10−3
a) α
PSfrag replacements
0

b)
Figure 6: Rigid footing test: Distribution of the hardening variable α in drained
elastic-plastic material in a) and undrained porous-elastic-plastic material in b).

lateral Raviart–Thomas-type enhanced-assumed-strain elements. The dimensions are

H = 4.758 m, W = 23.088 m and a = 4.587 m. The loading increment is ∆ū = 5·10−6 m.
The loading is linearly increased until a total displacement of ū = 0.0023 m is reached.
In Figure 6 the distribution of the hardening variable α for the drained elastic-plastic
material and undrained porous-elastic-plastic material is shown. It can be seen that the
plastic deformation is more pronounced in the case of the drained material. This leads to
the conclusion that the fluid within the material leads to an additional hardening effect.
The load-displacement curve in Figure 7 also shows this behaviour. Furthermore, a lower
permeability leads to more pronounced hardening. This can be explained by the fact that
the transport of the fluid is hindered and thus requires more work, see Figure 7.
Next, the distribution of the change of elastic fluid content me in the domain for
the undrained porous-elastic material is compared with the one of the undrained porous-
elastic-plastic material, see Figure 8.
Note that for the elastic material the fluid is squeezed out right underneath the area
where the loading is applied (see the negative change of the elastic fluid content at the
boundary of the applied footing depicted in Figure 8a). Opposed to that, in case of the
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 24

80

70

60

50

|F |/(MN/m)
40

30

20 drained elastic-plastic
undrained poro-elastic
10 undrained poro-elastic-plastic K
undrained poro-elastic-plastic 51 K
undrained poro-elastic-plastic 5K
0
0 0.0005 0.001 0.0015 0.002
|u|/m

Figure 7: Rigid footing: Load displacement curve for the different tested material
types. The load and displacement is taken from the area were the displacement is
prescribed, see Figure 5.

1.2 · 10−1
a) me /(kg/m3 )
PSfrag replacements
−5.6 · 10−1

b)
Figure 8: Rigid footing test. Change of the elastic fluid content me for undrained
porous-elastic material in a) and in b) for undrained porous-elastic-plastic material.

elastic-plastic material, the highest (negative) change of elastic fluid content is occurring
in a more diffuse region that also extends to the bulk (see Figure 8b). This is precisely the
area where most of the plastic deformation is happening, as can be observed in Figure 6
b). This phenomenon can be explained by the fact that the plastic deformation leads to
a positive change of the plastic fluid content. Due to the additive decomposition of the
change of the fluid content and the fact that no fluid is injected, the change of the elastic
fluid content becomes negative in the plastifying areas.
We depict the distribution of the change of the elastic me , the plastic mp and total
fluid content m for the undrained porous-elastic-plastic material in Figure 9. One can
observe that the change of the total fluid content m is strongly dominated by the change
of the plastic fluid content mp .
 D. Kienle & M.-A. Keip 25

1.2 · 10−1
me /(kg/m3 )

−5.6 · 10−1
a)
PSfrag replacements
2.8
mp /(kg/m3 )

0
b)

2.5
m/(kg/m3 )

−1.4 · 10−1
c)
Figure 9: Rigid footing test: Change of fluid content for undrained porous-elastic-
plastic material. Change of elastic fluid content me in a), change of plastic fluid content
mp in b) and change of total fluid content m in c).

Table 2: Material parameters for hydraulically induced fracture.

Lamé parameter λ = 180.0 GN/m2 Biot’s modulus M = 25.0 GN/m2

Lamé parameter G = 31.0 GN/m2 Biot’s coefficient b = 0.5 –
2
hardening modulus h = 5.0 MN/m fluid dyn. viscosity ηf = 1.0 · 10 −3
Ns/m2
saturated yield shift σy = 0.1 MN/m2 permeability K = 9.8 · 10−12 m3 s/kg
saturation ω = 2.0 – fluid density ρf = 1000.0 kg/m3
plastic viscosity ηp = 5 · 10−6 s crit. fracture energy ψc = 5.0 · 10−8 MN/m2
perturbation param. q1 = 2 · 10−5 MN/m2 length scale l = 0.5 m
slope yield function Mφ = 1.8 – residual stiffness k = 1 · 10 −5
–
2
position of peak smax = 2 · 10 −3
MN/m interpolation param. ǫ = 50 –

5.2 Comparison of different fracture driving forces for porous-

elastic-plastic fracturing

In the present example, we investigate the influence of the presented fracture driving
forces defined in (53) and (56). For that purpose, a squared domain with the dimensions
of 80 m × 80 m and a notch of the length a = 8 m in its center is considered, see Figure 10.
The fracture evolution is triggered by fluid injection into the notch.
Due to the symmetry of loading and geometry only one half of the domain is discretized
with 12, 060 quadrilateral Raviart–Thomas-type enhanced-assumed-strain elements. The
elements in the area surrounding the anticipated crack are refined yielding an element
size of he = 0.25 m in that region ([0 m, 40 m] × [31.875 m, 48.125 m] around the notch).
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 26

PSfrag replacements

˙
m̄ h̄
L L
a a/2

y
x

L L/2
Figure 10: Comparison of different fracture driving forces: Geometry and boundary
conditions. Due to the symmetry of both only one half of the specimen is discretized.
˙ = 0.01 kg/s. All edges are mechanically fixed and
Here the prescribed flux is h̄ = m̄
permeable.

3
The fluid injection is modeled by a prescribed fluid flux of h̄ = 0.01 kgs . The time step
is set to τ = 1 · 10−3 s and the material parameters are listed in Table 2.
The test was performed for the following two undrained settings with different choices
of fracture driving forces:
0 0
i) porous-elastic-plastic material with H = f (ψeff , ψplast ) according to (53)
0
ii) porous-elastic-plastic material with H = f (ψeff , wplast ) according to (56)
We now compare the hydraulically induced fracture lengths for the two porous-elastic-
plastic settings. As can be observed in Figure 11, both driving forces lead to the evolution
of cracks. The fracture evolution in consideration of the fracture driving force H =
0 0
f (ψeff , ψplast ) is however less prominent. Note that in that setting, only the elastic and
hardening energies contribute to the fracture driving force. Thus, we would not obtain
ductile fracture evolution in case of ideal plasticity with (h = 0, σy = 0; not investigated
here). In particular the latter observation justifies the presented modification of the
fracture driving force in (56).

5.3 Detailed analysis of hydraulically induced porous-elastic-

plastic fracture

Finally, we investigate the ductile fracture evolution driven by an injected fluid in detail.
The setup of the geometry and boundary conditions as well as the material parameters
are taken from the previous example (please refer to Figure 10 and Table 2).
The test was performed for two kinds of undrained materials:
i) porous-elastic material
0
ii) porous-elastic-plastic material with H = f (ψeff , wplast ).
 D. Kienle & M.-A. Keip 27

0 0
H(ψeff , ψplast )

0
PSfrag replacements

0
H(ψeff , wplast )

a) b)
Figure 11: Comparison of different fracture driving forces: Distribution of the frac-
ture phase-field for a porous-elastic-plastic with H = f (ψeff 0 , ψ0
plast ) and a porous-
0
elastic-plastic material with H = f (ψeff , wplast ) at two time steps: a) t = 0 s and b)
t = 90 s.

1
PSfrag replacements
d

a) b)
Figure 12: Hydraulically induced ductile fracture: Distribution of the fracture phase-
field in a) porous-elastic material in b) porous-elastic-plastic material at t = 90 s.
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 28

As can be seen in Figure 12, the length of the finally induced crack for the porous-elastic
material is much more pronounced that in case of the porous-elastic-plastic material. This
goes along with the observation of a higher fluid pressure inside the crack, see Figure 13.
We conclude that in case of an elastic-plastic material more fluid needs to be injected into
the crack to drive fracturing. For the elastic material we can observe a characteristic drop
of the pressure within the fracture at the onset of fracture propagation (injected fluid
volume V ≈ 0.00004 m3 ). In the elastic-plastic material this drop cannot be observed,
see again Figure 13.

0.14

0.12

0.1

0.08
p/Mpa

0.06

0.04

0.02
poro-elastic
poro-elastic-plastic
0
0 0.00015 0.0003 0.00045 0.0006 0.00075 0.0009
V /m3

Figure 13: Hydraulically induced ductile fracture: Fluid pressure within the fracture
at x = 1.0 m, y = 39.875 m over inject volume for porous-elastic and porous-elastic-
plastic material.

In Figure 14 the distribution of the change of the elastic fluid content is shown for the
final equilibrium state. The individual lengths of the cracks are clearly visible. Due to the
short crack length in the elastic-plastic case, the injected fluid is distributed over a smaller
region. This then gives rise to a higher change of the elastic fluid content, in particular
close to the fracture center. Note carefully that the change of the elastic fluid content in
front of the fracture tips is negative for the elastic-plastic material. This phenomenon is
investigated in a more detailed way in Figure 15, where the contributions of the change
of the fluid content in the elastic-plastic material are shown.
By taking a look at Figure 15c, it can be seen that most of the volumetric plastic
deformation occurs at the fracture tips (ṁp = ρf b tr ε̇p ). This leads to a positive change
of the plastic fluid content. Since the permeability in the bulk is comparably low, very
little amount of fluid diffuses from the fracture into the bulk. In other words, the change
of the fluid content at the fracture tips is almost zero. Due to the definition m = me + mp
the positive change of the plastic fluid content leads to a negative change of elastic fluid
content. Hence the fluid in the fully saturated medium, which is initially stored elastically,
is now stored plastically due to the plastic deformation of the solid matrix.
 D. Kienle & M.-A. Keip 29

4.2 · 10−2
PSfrag replacements
me /(kg/m3 )

−3.2 · 10−4

a) b)
Figure 14: Hydraulically induced ductile fracture: Change of elastic fluid content me
for porous-elastic material in a) and porous-elastic-plastic material in b).

PSfrag replacements
4.2 · 10−2 1.3 · 10−3

m, me /(kg/m3 ) mp /(kg/m3 )

−3.2 · 10−4 0.0

a) b) c)

Figure 15: Hydraulically induced ductile fracture: Change of total fluid content m
in a), change of elastic fluid content me in b) and change of plastic fluid content mp
in c) for porous-elastic-plastic material.

0.05 0.14 20
porous-elastic-plastic m porous-elastic-plastic poro-elastic-plastic
porous-elastic-plastic mp porous-elastic 18 poro-elastic
porous-elastic m 0.12
0.04 16
0.1
14
w/(1 · 10−6 m)

0.03 0.08
m/(kg/m3)

12
p/MPa

0.06 10
0.02 8
0.04
6
0.01 0.02
4
0
2
0
−0.02 0
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
x/m x/m x/m

Figure 16: Hydraulically induced ductile fracture: Change of fluid content m and
mp , fluid pressure p and fracture opening w for porous-elastic and porous-elastic-plastic
material over x at y = 40 m.

In Figure 16 we show the distribution of the change of fluid content as well as the
pressure and the fracture-opening width along x at y = 40 m. These distributions refer
to the final equilibrium state for both the porous-elastic and the porous-elastic-plastic
material. The difference in the crack length for the two different materials and the positive
change of the plastic fluid content as well as the negative change of the elastic fluid content
in front of the fracture tip in the elastic-plastic material can again be observed.
Finally, we depict a sequence of three snapshots in the course of the fracture evolution
of the porous-elastic-plastic material in Figure 17. To be specific, we show the fracture
phase field across the whole domain together with the change of fluid content, the pressure
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 30

0.05 0.05 0.05

m m m
mp mp mp
0.04 0.04 0.04

0.03 0.03 0.03

m/(kg/m3)

m/(kg/m3)
m/(kg/m3)
0.02 0.02 0.02

0.01 0.01 0.01

0 0 0
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
x/m x/m x/m

0.14 0.14 0.14

0.12 0.12 0.12

0.1 0.1 0.1

0.08 0.08 0.08

p/MPa

p/MPa
p/MPa

0.06 0.06 0.06

0.04 0.04 0.04

0.02 0.02 0.02

0 0 0

−0.02 −0.02 −0.02

0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
x/m x/m x/m
20 20 20

15 15 15
w/(1 · 10−6 m)

w/(1 · 10−6 m)

w/(1 · 10−6 m)
10 10 10

5 5 5

0 0 0
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
x/m x/m x/m

PSfrag replacements
1

a) b) c)
Figure 17: Hydraulically induced ductile fracture: Distribution of the fracture phase
field and change of fluid content m and mp , fluid pressure p and fracture opening w
over x at y = 40 m for three different time steps: a) t = 0 s, b) t = 45 s and c) t = 90 s.

and the fracture-opening width at three different time steps along x at y = 40 m. The first
time step is at t = 0 s, the second time step is at t = 45 s and the third time step is at the
final state (t = 90 s). As can be seen, all considered quantities are mainly concentrated
in the center of the fracture. Such a concentration is less prominent in case of an elastic
material, see Figure 16.
 D. Kienle & M.-A. Keip 31

6 Conclusion

A model for hydraulically induced fracturing in porous-elastic-plastic solids was developed

in the present work. It incorporates a phase-field approach to fracture that is combined
with a Drucker–Prager-type plasticity formulation and a Darcy–Biot-type fluid model.
The model exploits a variational structure leading to a global minimization formulation.
For this variational formulation it is crucial to introduce a plastic fluid content as an
additional unknown yielding a constitutive fluid pressure in terms of only the elastic
quantities. The global minimization structure demands the use of an H(div)-conforming
finite-element formulation, which has been implemented by means of Raviart–Thomas-
type shape functions. The locking phenomenon of the plasticity formulation is overcome
by using an enhanced-assumed-strain formulation for the deformation.
In the first numerical example a comparison of an undrained porous-elastic, an undrained
porous-elastic-plastic and a drained elastic-plastic formulation was performed. Here, the
different physical effects were investigated. It could be shown that the permeability of
porous media can be considered as an additional hardening parameter. The second exam-
ple shows the effect of the proposed modification of the fracture driving force. In the third
example, a hydraulically induced crack in a porous-elastic and a porous-elastic-plastic
medium was investigated. There it could be shown that neglecting the plastic effects
underestimates the pressure inside the fracture and overestimates the fracture length.
Acknowledgments. This work was funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Project Number 327154368 – SFB 1313. This
funding is gratefully acknowledged.

References
Aldakheel, F., Noii, N., Wick, T., Wriggers, P., [2020]. A global-local approach for hydraulic
phase-field fracture in poroelastic media. arXiv preprint arXiv:2001.06055.
Alessi, R., Marigo, J.-J., Maurini, C., Vidoli, S., [2018a]. Coupling damage and plasticity
for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional
examples. International Journal of Mechanical Sciences 149, 559–576.
Alessi, R., Vidoli, S., De Lorenzis, L., [2018b]. A phenomenological approach to fatigue
with a variational phase-field model: The one-dimensional case. Engineering Fracture
Mechanics 190, 53–73.
Ambati, M., Gerasimov, T., De Lorenzis, L., [2015]. Phase-field modeling of ductile frac-
ture. Computational Mechanics 55, 1017–1040.
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 32

Amor, H., Marigo, J., Maurini, C., [2009]. Regularized formulation of the variational brittle
fracture with unilateral contact: Numerical experiments. Journal of the Mechanics and
Physics of Solids 57, 1209–1229.
Armero, F., [1999]. Formulation and finite element implementation of a multiplicative
model of coupled poro-plasticity at finite strains under fully saturated conditions. Com-
puter Methods in Applied Mechanics and Engineering 171 (34), 205–241.
Bear, J., [1972]. Dynamics of Fluids in Porous Media. Dover Publications, New York.
Biot, M., [1941]. General theory of three-dimensional consolidation. Journal of applied
physics 12, 155–164.
Biot, M. A., [1984]. New variational-lagrangian irreversible thermodynamics with appli-
cation to viscous flow, reaction-diffusion, and solid mechanics. Advances in applied
mechanics. 24, 1–91.
Bluhm, J., De Boer, R., [1997]. The volume fraction concept in the porous media theory.
ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte
Mathematik und Mechanik 77 (8), 563–577.
Borden, M. J., Hughes, T. J. R., Landis, C. M., Anvari, A., Lee, I. J., [2016]. A phase-field
formulation for fracture in ductile materials: Finite deformation balance law deriva-
tion, plastic degradation, and stress triaxiality effects. Computer Methods in Applied
Mechanics and Engineering 312, 130–166.
Bourdin, B., Francfort, G. A., Marigo, J.-J., [2000]. Numerical experiments in revisited
brittle fracture. Journal of the Mechanics and Physics of Solids 48, 797–826.
Bourdin, B., Francfort, G. A., Marigo, J.-J., [2008]. The Variational Approach to Fracture.
Springer.
Cajuhi, T., Sanavia, L., De Lorenzis, L., [2018]. Phase-field modeling of fracture in variably
saturated porous media. Computational Mechanics 61 (3), 299–318.
Carrara, P., Ambati, M., Alessi, R., De Lorenzis, L., [2020]. A framework to model the fa-
tigue behavior of brittle materials based on a variational phase-field approach. Computer
Methods in Applied Mechanics and Engineering 361, 112731.
Darcy, H., [1856]. Les fontaines publiques de la ville de Dijon. Dalmont, Paris.
de Boer, R., [2000]. Theory of Porous Media. Springer, Berlin.
Detournay, E., Cheng, A. H.-D., [1993]. Fundamentals of poroelasticity. In: Fairhurst,
C. (Ed.), Comprehensive Rock Engineering: Principles, Practice and Projects, Vol. II,
Analysis and Design Method. Pergamon Press, Ch. 5, pp. 113–171.
Drucker, D. C., Prager, W., [1952]. Soil mechanics and plastic analysis or limit design.
 D. Kienle & M.-A. Keip 33

Quarterly of applied mathematics 10 (2), 157–165.

Ehlers, W., [2002]. Foundations of multiphasic and porous materials. In: Ehlers, W.,
Bluhm, J. (Eds.), Porous Media: Theory, Experiments and Numerical Applications.
Springer-Verlag, Berlin, pp. 3–86.
Ehlers, W., Luo, C., [2017]. A phase-field approach embedded in the theory of porous media
for the description of dynamic hydraulic fracturing. Computer Methods in Applied
Mechanics and Engineering 315, 348–368.
Ehlers, W., Luo, C., [2018a]. A phase-field approach embedded in the theory of porous
media for the description of dynamic hydraulic fracturing, part ii: The crack-opening
indicator. Computer Methods in Applied Mechanics and Engineering 341, 429–442.
Ehlers, W., Luo, C., [2018b]. A phase-field approach embedded in the theory of porous
media for the description of dynamic hydraulic fracturing, part ii: The crack-opening
indicator. Computer Methods in Applied Mechanics and Engineering 341, 429–442.
Francfort, G. A., Marigo, J.-J., [1998]. Revisiting brittle fracture as an energy minimiza-
tion problem. Journal of the Mechanics and Physics of Solids 46, 1319–1342.
Griffith, A. A., [1921]. The phenomena of rupture and flow in solids. Philosophical trans-
actions of the royal society of london. Series A, containing papers of a mathematical or
physical character, 163–198.
Heider, Y., Markert, B., [2017]. A phase-field modeling approach of hydraulic fracture in
saturated porous media. Mechanics Research Communications 80, 38–46.
Irwin, G. R., [1958]. Fracture. In: Flügge, S. (Ed.), Encyclopedia of Physics. Vol. 6,
Elasticity and Plasticity. Springer, pp. 551–590.
Johnson, E., Cleary, M. P., et al., [1991]. Implications of recent laboratory experimental
results for hydraulic fractures. In: Low permeability reservoirs symposium. Society of
Petroleum Engineers.
Kienle, D., Aldakheel, F., Keip, M.-A., [2019]. A finite-strain phase-field approach to
ductile failure of frictional materials. International Journal of Solids and Structures
172, 147–162.
Kuhn, C., Müller, R., [2010]. A continuum phase field model for fracture. Engineering
Fracture Machanics 77, 3625–3634.
Lambrecht, M., Miehe, C., [2001]. A note on formulas for localized failure of frictional
materials in compression and biaxial loading modes. International Journal for Numerical
and Analytical Methods in Geomechanics 25, 955–971.
Li, B., Peco, C., Millán, D., Arias, I., Arroyo, M., [2015]. Phase-field modeling and simu-
lation of fracture in brittle materials with strongly anisotropic surface energy. Interna-
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 34

tional Journal for Numerical Methods in Engineering 102 (3-4), 711–727.

Lo, Y.-S., Borden, M. J., Ravi-Chandar, K., Landis, C. M., [2019]. A phase-field model
for fatigue crack growth. Journal of the Mechanics and Physics of Solids 132, 103684.
Mauthe, S., Miehe, C., [2017]. Hydraulic fracture in poro-hydro-elastic media. Mechanics
Research Communications 80, 69–83.
Miehe, C., [1998]. A formulation of finite elastoplasticity based on dual co- and contra-
variant eigenvector triads normalized with respect to a plastic metric. Computer Meth-
ods in Applied Mechanics and Engineering 159, 223–260.
Miehe, C., Aldakheel, F., Raina, A., [2016a]. Phase field modeling of ductile fracture at
finite strains. a variational gradient-extended plasticity-damage theory. International
Journal of Plasticity 84, 1–32.
Miehe, C., Hofacker, M., Schänzel, L.-M., Aldakheel, F., [2015a]. Phase field modeling
of fracture in multi-physics problems. Part II. brittle-to-ductile failure mode transition
and crack propagation in thermo-elastic-plastic solids. Computer Methods in Applied
Mechanics and Engineering 294, 486–522.
Miehe, C., Hofacker, M., Welschinger, F., [2010a]. A phase field model for rate-independent
crack propagation: Robust algorithmic implementation based on operator splits. Com-
puter Methods in Applied Mechanics and Engineering 199, 2765–2778.
Miehe, C., Kienle, D., Aldakheel, F., Teichtmeister, S., [2016b]. Phase field modeling of
fracture in porous plasticity: A variational gradient-extended eulerian framework for
the macroscopic analysis of ductile failure. Computer Methods in Applied Mechanics
and Engineering 312, 3–50.
Miehe, C., Mauthe, S., Teichtmeister, S., [2015b]. Minimization principles for the coupled
problem of darcy-biot-type fluid transport in porous media linked to phase field modeling
of fracture. Journal of the Mechanics and Physics of Solids 82, 186–217.
Miehe, C., Teichtmeister, S., Aldakheel, F., [2016c]. Phase-field modeling of ductile frac-
ture: A variational gradient-extended plasticity-damage theory and its micromorphic
regularization. Philisophical Transactions of the Royal Society A 374.
Miehe, C., Welschinger, F., Hofacker, M., [2010b]. Thermodynamically consistent phase-
field models of fracture: Variational principles and multi-field fe implementations. In-
ternational Journal for Numerical Methods in Engineering 83, 1273–1311.
Mikelic, A., Wheeler, M. F., Wick, T., [2015a]. A phase-field method for propagating
fluid-filled fractures coupled to a surrounding porous medium. Multiscale Modeling &
Simulation 13 (1), 367–398.
Mikelic, A., Wheeler, M. F., Wick, T., [2015b]. A phase-field method for propagating
 D. Kienle & M.-A. Keip 35

fluid-filled fractures coupled to a surrounding porous medium. Multiscale Modeling &

Simulation 13 (1), 367–398.
Mikelic, A., Wheeler, M. F., Wick, T., [2015c]. Phase-field modeling of a fluid-driven frac-
ture in a poroelastic medium. Computational Geosciences, Springer Verlag (Germany).
Pham, K., Amor, H., Marigo, J.-J., Maurini, C., [2011]. Gradient damage models and
their use to approximate brittle fracture. International Journal of Damage Mechanics
20 (4), 618–652.
Pise, M., Bluhm, J., Schröder, J., [2019]. Elasto-plastic phase-field model of hydraulic
fracture in saturated binary porous media. International Journal for Multiscale Compu-
tational Engineering 17 (2).
Raviart, P. A., Thomas, J. M., [1977]. Primal hybrid finite element methods for 2nd order
elliptic equations. Mathematics of computation 31 (138), 391–413.
Schreiber, C., Kuhn, C., Müller, R., Zohdi, T., [2020]. A phase field modeling approach of
cyclic fatigue crack growth. International Journal of Fracture, 1–12.
Simo, J. C., Rifai, S., [1990]. A class of mixed assumed strain methods and the method of
incompatible modes. International Journal for Numerical Methods in Engineering 29,
1595–1638.
Steinke, C., Zreid, I., Kaliske, M., [2020]. Modelling of ductile fracture of strain-hardening
cement-based composites-novel approaches based on microplane and phase-field method.
In: Plasticity, Damage and Fracture in Advanced Materials. Springer, pp. 175–199.
Storm, J., Supriatna, D., Kaliske, M., [2020]. The concept of representative crack elements
for phase-field fracture: Anisotropic elasticity and thermo-elasticity. International Jour-
nal for Numerical Methods in Engineering 121 (5), 779–805.
Teichtmeister, S., Kienle, D., Aldakheel, F., Keip, M.-A., [2017]. Phase field modeling of
fracture in anisotropic brittle solids. International Journal of Non-Linear Mechanics 97,
1–21.
Teichtmeister, S., Mauthe, S., Miehe, C., [2019]. Aspects of finite element formulations
for the coupled problem of poroelasticity based on a canonical minimization principle.
Computational Mechanics 64 (3), 685–716.
Terzaghi, K., [1925]. Erdbaumechanik auf bodenphysikalischer Grundlage. F. Deuticke.
Vermeer, P., de Borst, R., [1984]. Non-associated plasticity for soils, concrete and rock.
Heron 29 (3), 1–64.
Wilson, Z. A., Landis, C. M., [2016]. Phase-field modeling of hydraulic fracture. ICES
Report 16-10.
 A variational minimization formulation for hydraulically induced fracturing in elastic-plastic solids 36

Wu, T., De Lorenzis, L., [2016]. A phase-field approach to fracture coupled with diffusion.
Computer Methods in Applied Mechanics and Engineering 312, 196–223.
Yin, B., Kaliske, M., [2020]. A ductile phase-field model based on degrading the fracture
toughness: Theory and implementation at small strain. Computer Methods in Applied
Mechanics and Engineering 366, 113068.