Computational Mechanics (2022) 69:151–175

[Link]

ORIGINAL PAPER

Crack phase-field model equipped with plastic driving force and

degrading fracture toughness for ductile fracture simulation
Jike Han1 · Seishiro Matsubara2 · Shuji Moriguchi3 · Michael Kaliske4 · Kenjiro Terada3

Received: 20 April 2021 / Accepted: 21 August 2021 / Published online: 23 September 2021
© The Author(s), under exclusive licence to Springer-Verlag GmbH Germany, part of Springer Nature 2021

Abstract
This study presents a novel phase-field model for ductile fracture by the introduction of both the plastic driving force
and the degrading fracture toughness into crack phase-field computations based on the phenomenological justification for
ductile fracture in elastoplastic materials. Assuming that the constitutive work density consists of elastic, pseudo-plastic and
crack components, we derive the governing equations from local and global optimization problems within the continuum
thermodynamics framework. In addition to the elastic strain energy, the plastic strain energy also works as a driving force to
sustain damage evolution. Additionally, we introduce a degrading fracture toughness to reflect the evolution of micro-defects
and their coalescences into each other that are caused by accumulated plastic deformation. Equipped with these ingredients,
the proposed model realizes the reduction of both stiffness and fracture toughness to simulate the failure phenomena of
elastoplastic materials. Several numerical examples are presented to demonstrate the capability of the proposed model in
reproducing some typical ductile fracture behaviors. The findings and perspectives are subsequently summarized.

Keywords Crack phase-field · Ductile fracture · Plastic driving force · Degrading fracture toughness · Gradient plasticity

1 Introduction tal and numerical studies have been focused on understanding

failure phenomena. The underlying principle, in theory, is
Fracture is one of the most critical failure mechanisms in var- that brittle fracture evolution requires a portion of stored
ious engineering fields. Since Griffith [1] laid the foundation strain energy to generate a new pair of crack surfaces when
of classical fracture mechanics, a vast number of experimen- the crack evolving condition is satisfied; that is, when an
instantaneous energy release rate G∗ is equal or larger than
B Kenjiro Terada the critical value (G∗ ≥ Gc ), the crack begins to propagate. In
tei@[Link] case of ductile materials, like steels, since plastic deformation
Jike Han accompanied by hardening occurs before crack initiation, a
[Link].s7@[Link] portion of energy is dissipated by the sizable inelastic defor-
Seishiro Matsubara mation before fracture.
[Link]@[Link] Within the framework of the finite element method (FEM),
Shuji Moriguchi approaches to represent crack propagations can be classi-
s_mori@[Link] fied into two major categories. One is a group of methods
Michael Kaliske that incorporate discontinuous displacement fields into FE
[Link]@[Link] approximations, and the other is a group of continuum
damage models that represent a crack by strain-softening
1 Department of Civil and Environmental Engineering, Tohoku behavior associated with fracture by establishing appropri-
University, Sendai, Japan
ate constitutive laws. Early developments for the former
2 Department of Mechanical Systems Engineering, Nagoya include the node-releasing FEM [2,3], the embedded FEM
University, Nagoya, Japan
[4], the X-FEM [5,6], the partition-of-unity-based FEM [7]
3 International Research Institute of Disaster Science, Tohoku and the finite cover method (FCM) [8]. While these meth-
University, Sendai, Japan
ods represent discrete cracks with high accuracy, their ability
4 Institute for Structural Analysis, Technische Universität to trace complex crack propagations in three-dimensional
Dresden, Dresden, Germany

123
152 Computational Mechanics (2022) 69:151–175

cases is limited. The continuum damage models, in which into the crack phase-field modeling to establish the crack
discrete cracks are replaced by the material behavior, have kinematics.
also been studied for a long time because of their amenabil- Implementation of phase-field crack models covers a
ity to FEM. Some of these models have long histories: broad range of discretization methods. In the higher-order
the Gurson-Tverg-aard-Needleman (GTN) model [9,10], the phase-field formulation presented by Borden et al. [31], more
Lemaitre model [11] and the microplane damage model regular and faster converging solutions of the variational
[12] were first reported in the 1980s. However, these dam- problem are obtained by adopting the isogeometric analysis
age models are known to suffer from numerical instabilities (IGA) framework [32] as well. Meanwhile, Kakouris et al.
caused by ill-posed partial differential equations [13]. This [33] incorporated the crack phase-field theory into the mate-
deficit inevitably yields pathological mesh-dependent solu- rial point method (MPM), and showed valuable comparisons
tions, and even results in diverged solutions immediately between FEM and MPM. Also, Roy et al. [34] introduced the
after the onset of damage evolution. Therefore, some reg- crack phase-field theory into the peridynamics framework.
ularization techniques are required to overcome these issues, While most of the developments target quasi-static states,
which have been reported to originate from non-local integral some approaches to dynamics have also been proposed [35–
algorithms [14,15] and gradient enhanced approximations 38]. Additionally, a study on the degradation function was
[16,17], respectively. conducted by Sargado et al. [39], and an intensive inves-
During the last decade or two, one of the non-local tigation study for the length scale parameter was provided
methods, the “crack phase-field”, has received considerable by Zhang et al. [40]. Tanné et al. [41] conducted numerical
attention due to its ability to predict arbitrary crack propaga- investigations focussing on the crack nucleation of several
tion and its compatibility with classical fracture mechanics. phase-field models.
The pioneering work was done by Francfort and Marigo [18], Phase-field approaches to ductile fracture were indepen-
who take a variational approach to associate brittle fracture dently pioneered by Ambati et al. [42,43] and Miehe et al.
with the minimization problem of potential energy based on [44,45]. On the understanding that functionals relevant to
the aforementioned Griffith theory. The geometry represen- plastic deformation are not available, Ambati et al. [42] devel-
tation of a crack is subsequently regularized by adopting a oped a coupling strategy between a phase-field brittle model
phase-field approximation so that the variational problem is and an elastoplastic constitutive law by making the degra-
solved with little difficulty [19–21]. As an alternative, Miehe dation function dependent on plastic deformation to realize
et al. [22,23] proposed a formulation based on thermody- ductile fracture. Here, the fundamental idea is that the amount
namics within the continuum mechanics framework and also of degradation depends on both the phase-field parameter d
introduced a history variable to represent the maximum elas- and the accumulated plastic strain α. This recasting of the
tic strain energy in tension ever experienced to ensure the degradation function makes it possible to control the damage
irreversibility of crack propagations. The introduction of this evolution according to the increase of plastic deformation. In
history variable prevents the imposition of Karush-Kuhn- contrast, Miehe et al. [44] introduced a pseudo-plastic energy
Tucker condition on damage evolution and results in a simple based on variational gradient plasticity [46,47] so that plastic
and robust staggered solution scheme for the coupled two- deformations also contribute to a driving force for damage
field problem. In the sequel, various achievements have been evolution. In another attempt to redefine the driving force,
accomplished for phase-field modeling of brittle fracture. A Borden et al. [48] introduced the effect of stress triaxiality
comprehensive review was published by Ambati et al [24]. to plastic strain energy as a parameter to control the damage
More recently, Wu et al. [25,26] proposed peculiar phase- evolution. A review of several major phase-field models of
field models based on cohesive zone model. ductile fracture was done by Alessi et al. [49]. More recently,
Particular attention is paid to the definition of fracture driv- Dittmann et al. [50] and Yin and Kaliske [51] made the frac-
ing forces among the intense studies on the crack phase-field ture toughness into a decreasing function of accumulated
modeling. Henry and Levine [27] identified the crack states plastic strain α as G c (α) to consider the effect of plastic
by decomposing a strain tensor into tensile and compres- deformation.
sive parts, while Amor et al. [28] proposed an energy split In the present work, we propose a new crack phase-field
based on the volumetric-deviatoric decomposition. Miehe model for ductile fracture, into which both the plastic driving
et al. [22] defined the driving force based on the spectral force and degrading fracture toughness are introduced along
split principle. Along with a concise review for some popu- with variational gradient enhanced plasticity. On the assump-
lar strategies, Steinke and Kaliske [29] proposed a split based tion that the constitutive work density consists of elastic and
on the directional stress decomposition and Storm et al. [30] pseudo-plastic components (for convenience, it is referred to
has introduced a homogenization-type approach by incorpo- simply as “plastic” hereafter) as well as crack component,
rating the concept of a representative crack element (RCE) the governing equations are derived from local and global
optimization problems within the continuum thermodynam-

123
Computational Mechanics (2022) 69:151–175 153

ics framework. In addition to the elastic strain energy, the

plastic strain energy also works as a driving force to sustain
damage evolution. Also, we originally introduce a degrading
fracture toughness to reflect the evolution of micro-defects
and their coalescences into each other due to accumulated
plastic deformation. Thanks to these ingredients, the pro-
posed model is capable of representing the reduction of both
stiffness and fracture toughness within a damaged region to
simulate the failure phenomena of elastoplastic materials.
The structure of this paper is as follows: Sect. 2 pro- (a) (b)
vides a review of an issue concerning the straightforward
incorporation of elastoplastic material behavior into con-
ventional crack phase-field modeling. Then, a strategy to
settle the issue, which is based on the phenomenological
justification for ductile fracture, is also explained. Section
3 is devoted to the formulation, starting with the definition
of global state variables. The governing equations are then
derived by solving the corresponding local and global opti- (c)
mization problems. Also, the thermodynamic consistency is Fig. 1 Schematic diagram of stress–strain relations for three deforma-
demonstrated in this section. In Sect. 4, a spatial discretiza- tion cases
tion is provided within the IGA framework. In Sect. 5, several
numerical examples are presented to assess the performance
yields
of the proposed model and demonstrate its capability in
reproducing some typical ductile fracture behaviors. Finally,
in Sect. 6, the findings are summarized, and the contribution Ψ̂0e+
d= , (2)
of this study is considered. G c /2lf + Ψ̂0e+

which clearly implies that the increase of stored elastic strain

energy Ψ̂0e+ directly determines the damage evolution. In
2 Issue on crack phase-field modeling for other words, as the area of the gray-colored region in Fig. 1a
ductile fracture increases, the damage evolution becomes more significant.
When plastic deformation is considered, the straightfor-
As well known, an issue is raised when applying conventional ward incorporation of elastoplastic material behavior into
crack phase-field modeling to ductile fracture in elastoplastic conventional phase-field models cannot provide adequate
materials. In this section, we address it in a one-dimensional damage evolution in the above setting for brittle fracture.
setting and devise our strategy that can be justified from the For example, if the material exhibits a linear plastic hard-
phenomenological viewpoint. ening response as depicted in Fig. 1b, the increase rate of
We begin with providing the following governing equa- elastic strain energy density ΔΨ̂0e+ is considerably less than
tions of phase-field modeling for brittle fracture, which are that of the purely elastic case. Therefore, to obtain the same
consistent with those in Miehe et al. [22,23]: amount of damage evolution as that computed in the elastic
case, an extremely large amount of plastic deformation εp
is required. Additionally, in the elastic-perfectly plastic case
⎧
⎨∇ · σ = 0, shown in Fig. 1c, no damage growth occurs during the plas-
G   (1) tic deformation because ΔΨ̂0e+ = 0. Adequate modifications
⎩ c d − lf2 d = 2 (1 − d) Ψ̂0e+ . are therefore required to realize a realistic damage evolution
lf
for ductile fracture.
Typically, an elastoplastic material experiences the tran-
Here, ∇ and  represent the first-order and second-order sitional process from the intact state to absolute failure as
spatial gradients, respectively, and σ is the stress in a general depicted in Fig. 2. This can be itemized as follows:
sense. Also, G c , lf , d and Ψ̂0e+ are the fracture toughness,
crack length scale parameter, phase-field variable and elastic A Void nucleation: When the accumulation of dislocations
strain energy in tension, respectively. In case of unidirectional reaches its limit under tensile loading, small defects such
uniform deformation, i.e., d = 0, the phase-field variable as cracks and voids are generated at the micro-scale and

123
154 Computational Mechanics (2022) 69:151–175

Fig. 2 Micro-mechanisms of
ductile fracture and its
numerical approximation with
crack phase-field modeling

evolve as the amount of macroscopic plastic deformation evolution. The situation is schematically shown in Fig. 2.
increases; As can be seen, within the region of localized plastic defor-
B Void coalescence: Micro-defeccts selectively coalesce mation, the evolution of small defects causes the decrease
into larger ones, resulting in the macroscopic reduction of in effective surfaces at the micro-scale, which results in the
material stiffness because of the decrease of the effective corresponding reduction of load bearing capacity. This kind
area that can sustain the stress; of phenomena is commonly represented by the reduction
C Crack propagation: Before long, some of the coales- of material properties according to the effective stress con-
cent defects rapidly grow into visible cracks around the cept in conventional continuum damage mechanics (CDM)
stiffness-reduced domain at the macro-scale; [52,53]. In analogy with this, the degradation function can be
D Separation: Complete material separations are observed a multiplier of the energy functions as formulated in crack
and the material totally loses the resistance to external phase-field modeling [44]. It is thus reasonable to assume
forces. that the fracture toughness also gradually decreases as the
damage evolves, and eventually becomes zero at the fully
Thus, from the phenomenological viewpoint, the accumu- broken state from the phenomenological viewpoint.
lation of plastic deformation induces crack growths at the The noteworthy aspect, again, is that a sufficient amount of
macro-scale, and at the same time serves as a driving force plastic deformation induces the generation of micro-cracks
for damage evolution from the computational viewpoint. To and voids, which causes the deterioration of material prop-
this end, strategies for incorporating plastic deformation into erties. Therefore, we assume that the fracture toughness is
the driving force have been proposed by Miehe et al. [44,45] degraded by the accumulated plastic strain α in this study.
and Borden et al. [48]. Although, to the best of our knowledge, no such discussion
However, the crack growth mechanism corresponding to has been published in the literature, Dittmann et al. [50]
the progressive deterioration of the material seems to be over- and Yin and Kaliske [51] present different function forms of
looked in the previous crack phase-field models. Specifically, fracture toughness which decrease with plastic deformation.
the fracture toughness G c has been treated as being con- Hence, the present study is the first to introduce the multiple
stant throughout the computations. If the transitional process effects of plastic deformation, i.e., the plastic driving force
from intact to fully broken regions is infinitesimal as in brit- and degrading fracture toughness, into the modeling.
tle fracture, such a treatment is feasible and leads to lots of
archeivements [24]. However, when large plastic deforma-
tions are involved, the material behavior should adequately
reflect the so-called process zone, in which the material
properties are gradually changed due to void nucleation and

123
Computational Mechanics (2022) 69:151–175 155

3 Crack phase-field modeling for ductile

fracture with variational gradient
enhanced plasticity

The phase-field modeling for brittle fracture is based on

the variational principle that minimizes the total energy, as
demonstrated by Francfort and Marigo [18], and has subse-
quently been regularized by Bourdin et al. [21]. In the present
work, we propose phase-field modeling for ductile fracture
coupled with variational gradient enhanced plasticity accom-
modating a driving force for both the damage evolution and Fig. 3 Diffusive crack topology in crack phase-field modeling
the degradation of fracture toughness. After prescribing state
variables, we define the energy functions for the formulation
based on thermodynamics. Then, the thermodynamic con- on the boundary of the body in the reference configuration,
sistency is demonstrated, and the governing equations are ∂B0 = ∂Bα0 ∪ ∂B∇α 0 . In the weak formulation presented in
derived as stationary conditions of these energy functions. the next section, only the Neumann condition is adopted for
convenience.
3.1 State variables Meanwhile, following the crack phase-field theory, we
introduce the phase-field variable d (X, t) ∈ [0, 1] as a
Let B0 and Bt be the reference and current configurations of global variable to describe the damage evolution that approx-
a solid body, respectively, and let initial material positions imates the crack topology. The phase-field approximation is
X ∈ B0 be mapped onto current positions x ∈ Bt at time to regularize the geometries of discrete discontinuities, and
t by motions x = ϕ (X, t), so that those at time t0 can be d evolves according to relevant driving forces, as will be
identified with X = ϕ (X, t0 ). The deformation gradient is described later, along with the irreversible and initial condi-
defined as F := ∇ϕ (X, t) = ∂ x/∂ X with its determinant tions, ḋ ≥ 0 and d (X, t0 ) = 0, respectively. As shown in
J := det [F] > 0. The Dirichlet and Neumann conditions are Fig. 3, the crack phase-field represents a regularized crack
imposed on the boundary surface of the body in the reference surface such that
ϕ
configuration, ∂B0 = ∂B0 ∪ ∂B0T as
d 2 + lf2 ||∇d||2
Γlf (d) = γlf dV with γlf = , (6)
ϕ B0 2lf
ϕ = ϕ̄ (X, t) on ∂B0 and P · N = T (X, t) on ∂B0T ,
(3)
where γlf is a second-order functional of crack surface den-
sity and lf is the crack length scale parameter to represent
where ϕ̄ and T are a prescribed motion and an external
the width of a regularized crack domain. As pointed out
force vector, respectively. Also, P and N are the first Piola-
by Borden et al. [31], the actual crack surface domain Γ
Kirchhoff stress tensor and the outward unit normal vector
is approximately obtained by the volumetric domain Γlf as
on the boundary in the reference configuration.
Γ ≈ lim Γlf . Possible Dirichlet and Neumann conditions
Under consideration of incompressible plastic deforma- lf →0
tion, the deformation gradient F is multiplicatively decom- are given as
posed as
d = 1 on Γ ⊂ B0 and ∇d · N = 0 on ∂B∇d
0 , (7)
 
F=F ·F e p
with det F p
= 1, (4)
the latter of which is solely used in this study.
where Fe and Fp are the elastic and plastic deformation gra- In summary, the multi-field setting reads the following set
dients, respectively. Within the framework of the variational of global field variables:
gradient enhanced plasticity theory, we employ the accumu-
lated plastic strain α (X, t) as a global field variable in this U := {ϕ, α, d} (8)
study, which satisfies irreversibility α̇ ≥ 0 along with the
initial condition α (X, t0 ) = 0. Also, the Dirichlet and Neu- Also, the set of state variables for the subsequent constitutive
mann conditions are employed as in the literature [22] as equations yields

α = 0 on ∂Bα0 and ∇α · N = 0 on ∂B∇α

0 (5) C := F, F p , α, ∇α, d, ∇d . (9)

123
156 Computational Mechanics (2022) 69:151–175

Nevertheless, for easier numerical implementation, we employ function forms for degradation are subjective as far as they
the idea of micromorphic regularization outlined in [54–56] meet the requirement in Eq. (13).
to only name a few, so that the modified set of state variables To describe the elastic response, we employ the standard
is set to be neo-Hookean constitutive law that consists of volumetric and
deviatoric parts as follows:
C := F, F p , ᾱ, α, ∇α, d, ∇d , (10)
Ψ̂0e = Ψ̂0,vol
e
+ Ψ̂0,dev
e
for the phase-field modeling of ductile fracture, where ᾱ is the  
local accumulated plastic strain corresponding to the global κ J e2 − 1 μ  (15)
= − lnJ e + I e −3 ,
counterpart α. 2 2 2 b̄

Remark 1 Recently, various damage models within the finite

where μ and κ are the elastic shear and bulk moduli, respec-
deformation framework, including the gradient-enhanced e
tively, and I b̄e is the first invariant of b̄ = J e−2/3 be with be
approach or micromorphic regularization, have been pro-
denoting the elastic left Cauchy-Green tensor. Then, the sec-
posed to realize mesh-independent solutions; see References
ond Piola-Kirchhoff stress tensor S0 for the intact material
[57–60] to only name a few.
is obtained straightforwardly from the derivative of Ψ̂0e with
respect to the right Cauchy-Green tensor C as
3.2 Energy functions

First, following some previous works [44,45,48], we define ∂ Ψ̂0,e dev

∂ Ψ̂0,e vol
S0 = 2 +2 = S0, vol + S0, dev
the constitutive work density functional as ∂C ∂C  
κ  e2  I e
  = J − 1 C −1 + μJ e−2/3 C p−1 − b C −1 ,
Ψ̂ (C) = Ψ̂ e F, F p , d + Ψ̂ p (ᾱ, α, ∇α, d) + Ψ̂ f (α, d, ∇d), 2 3
   
stored energy accumulated dissipation (16)
(11)
where C p is the plastic part of right Cauchy-Green tensor. By
which consists of elastic Ψ̂ e , plastic Ψ̂ p and crack Ψ̂ f
compo- the standard push-forward process, the corresponding Kirch-
nents. The elastic contribution follows a tensile-compressive hoff stress tensor yields
split as
κ  e2   e
  τ 0 = F · S0 · F T = J − 1 1 + μdev b̄ . (17)
Ψ̂ e
F, F p , d = g (d) Ψ̂0e+ + Ψ̂0e− , (12) 2

where Ψ̂0e+ and Ψ̂0e− are the undamaged tensile and com- Following the assumption above, the degraded elastic strain
pressive components, respectively, and g (d) is a degradation energy and Kirchhoff stress are obtained as, respectively,
function that represents the reduction of material stiffness.
Following the general description of crack phase-field, the Ψ̂ e = g (d) Ψ̂0e+ + Ψ̂0e−
⎧  
degradation has been assumed only for the tensile compo- ⎪
⎪ κ J e2 − 1 μ 
⎪
⎪g (d) − lnJ + g (d)
e I e − 3 for J e ≥ 1,
nent of the elastic strain energy Ψ̂0e+ . For the degradation ⎨ 2 2 2 b̄
=  
function, the following conditions are necessarily imposed: ⎪ κ J e2 − 1  
⎪
⎪ − e + g (d) μ I e − 3 for J e < 1,
⎪
⎩ lnJ
  2 2 2 b̄
1 for d = 0 ∂ g (d)  (18)
g (d) = and ≤ 0,
0 for d = 1 ∂d 0≤d≤1
(13) and

which ensure the monotonic decrease of the material stiff-

τ = g (d) τ +
0 + τ0
−
ness. The most widely-used choice of the degradation func- ⎧ κ  e2   e
tion is the quadratic form as ⎨g (d) J − 1 1 + g (d) μdev b̄ for J e ≥ 1,
= κ 2 e2   e
⎩ J − 1 1 + g (d) μdev b̄ for J e < 1.
g (d) = (1 − d)2 , (14) 2
(19)
which was originally introduced by Bourdin et al. [21]. Also,
several different degradation functions have been reported Next, following Miehe et al. [55,56], the gradient-type
in the literature [29,61,62]. To the best of our knowledge, second-order functional associated with plastic deformation

123
Computational Mechanics (2022) 69:151–175 157

is introduced as elastic and plastic parts as

e p
Ψ̂ p (ᾱ, α, ∇α, d) l = Ḟ · F e−1 + F e · Ḟ · F
  ·F
p−1 e−1
= l e + l p, (23)
 
ᾱ   lp2 pp :=L p
= g (d) ŷ ᾱ˜ d ᾱ˜ + y0 ||∇α||2 + (ᾱ − α)2 ,
0 2 2 when the body experiences elastoplastic deformations. It
(20) should be noted that l p is the push-forward of the plas-
tic velocity gradient L p from the intermediate configuration
where lp ≥ 0 is the plastic length scale parameter in gra- defined by the multiplicative decomposition Eq. (4). Also,
dient plasticity, y0 is the initial yield stress, and pp plays a based on this decomposition, the elastic left Cauchy-Green
role of a penalty parameter to link the local variable ᾱ to tensor be yields
the global micromorphic field variable α. Also, ŷ (ᾱ) is the
empirically defined hardening function. A point to remember be = F e · F eT = F · C p−1 · F T
(24)
here is that, as illustrated in Fig. 2, since the discrete crack with C p = F pT · F p and C p (X, t0 ) = 1.
neither stores the energy nor sustains stress, the elastic and
plastic components of the constitutive work density in fully On the other hand, the evolution of the stored energy is
broken regions must be zero, i.e., Ψ̂ e + Ψ̂ p = 0. In order represented by the material time derivative as
to approximate this situation in our crack phase-field mod-
eling, we have multiplied the degradation function g (d) by d e ∂ Ψ̂ e e ∂ Ψ̂ e
p
the pure plastic contribution Ψ̂0 as in this expression. This Ψ̂ = e : ḃ + ḋ. (25)
dt ∂b ∂d
setup corresponds to the postulate made in Sect. 2.
Finally, we define the crack contribution to the constitutive Here, the material time derivative of be is expanded as
work density as
e  
ḃ = l · be + be · l T + L be
  (26)
d 2 + lf2 ||∇d||2 = l · be + be · l T + F · Ċ
p−1
· FT,
Ψ̂ (α, d, ∇d) = G c (α)
f
, (21)
2lf
 
where L be denotes the Lie derivative of be . Then, the first
where G c (α) is the fracture toughness that is supposed to term of Eq. (25) reads
degrade depending on the amount of plastic deformation,      
as mentioned in Sect. 2. For the sake of convenience, we ∂ Ψ̂ e e ∂ Ψ̂ e e l + lT ∂ Ψ̂ e e
: ḃ = 2 e · b : − 2 e · b : dp
assume that G c (α) is a decreasing function of a scalar vari- ∂ be ∂b 2 ∂b
able α ∗ calculated from the accumulated plastic strain. The  
function form of the fracture toughness is specified in Sect. = τ : d − dp ,
5. Although similar ideas of degrading fracture toughness are (27)
presented in the literature [50,51], the distinctive feature of
the proposed model is obvious since the amount of plastic where d p is the symmetric Eulerian plastic rate of deforma-
deformation is employed as a driving force for both the dam- tion tensor defined as
age evolution and the degradation of fracture toughness as   1 p−1
advocated in Sect. 2. d p := F e · sym L p · F e−1 = − F · Ċ · F T · be−1
2
1   1
· C p · F −1 .
p−1
= − L be · be−1 = − F · Ċ
3.3 Thermodynamic consistency 2 2
(28)
Following some previous studies [45,50], the thermodynamic
consistency of the proposed model is confirmed based on the In consideration of Eq. (25) along with Eq. (27), the dissipa-
following energy dissipation rate functional tion rate functional Eq. (22) yields

d e Dpf = τ : d p + τf ḋ, (29)

Dpf = τ : d − Ψ̂ (C) , (22)
dt
where we have defined τf as
where τ is the Kirchhoff stress tensor conjugate to the sym-
metric Eulerain rate of the deformation tensor d =: sym [l]. ∂ Ψ̂ e
τf := − ≥ 0. (30)
Here, the spatial velocity gradient l is decomposed into the ∂d

123
158 Computational Mechanics (2022) 69:151–175

In this study, we assume a threshold function for plasticity where λp is the viscoplastic multiplier defined as
as
 3
  2 p λp := Φ̂ p . (37)
Φ̂ τ dev , r p := ||τ dev || −
p
r (31) 2ηp
3

in which the thermodynamic force r p associated with hard- Also, from Eq. (37), we have the loading-unloading condition
ening (thermodynamic stress that is energy-conjugated to the as
internal variable α) is defined as
2
λp ≥ 0, Φ̃ p = Φ̂ p − ηp λp ≤ 0, λp Φ̃ p = 0, (38)
∂ Ψ̂ p 3
r :=
p
= g (d) ŷ (ᾱ) + pp (ᾱ − α) , (32)
∂ ᾱ
where, in particular, Φ̃ p = 0 corresponds to the plastic yield
Similarly, we assume a threshold function for fracture as condition; see, e.g., Reference [64] for a similar setting. On
  the other hand, the phase-field variable satisfies irreversibility
Φ̂ f τf − r f := τf − r f (33) ḋ ≥ 0 thanks to the linearity of Φ̂ f with respect to τf − r f .
Lastly, by substituting Eq. (36) into Eq. (29), we can ensure
where r f is the damage resistance force (thermodynamic the thermodynamic consistency as
stress that is energy-conjugated to the internal variable d)
and is defined as follows: ∂ Φ̂ p Φ̂ f ∂ Φ̂ f
Dpf = τ : d p + τf ḋ = τ : λp + τf  
∂τ ηf ∂ τf − r f
∂ Ψ̂ p ∂ Ψ̂ f d ∂ Ψ̂ f
r f := + − · . (34) Φ̂ f
∂d ∂d d X ∂∇d = λp ||τ dev || + τf ≥ 0,
ηf
To prescribe the sufficient conditions for the Clausius- (39)
Planck inequality in Eq. (29), we follow the principle of
maximum dissipation as Remark 2 Thermodynamic consistency is often confirmed
  using Ψ̂ instead of Ψ̂ e in Eq. (22); see References [53,63]
V̂ Ċ, τ , r p , τf − r f to only name a few. Specifically, from the viewpoints of
   phenomenological (continuum mechanical) modeling, the
= sup τ : d p − r p ᾱ˙ + τf − r f ḋ dissipation rate can also be defined as
τ ,r p ,τf −r f (35)
    
−
3
Φ̂ p τ dev , r p 2 −
1
Φ̂ f τf − r f 2 . d  
4ηp 2ηf Dpf = τ : d − Ψ̂ (C) = τ : d p − r p ᾱ˙ + τf − r f ḋ,
dt
(40)
Here, the 4th and 5th terms in the brackets are the penalization
functions, which control the plastic deformation and damage
evolution with the numerical viscosities ηp > 0 and ηf > 0, from which the first three terms in the right-hand side of Eq.
(35) are straightforwardly obtained by some tensor algebraic
respectively. Also, the convex functions Φ̂ p and Φ̂ f are placed
manipulations.
in the Macaulay bracket • ± = (• ± | • |) /2. It should be
noted here that positive values of Φ̂ p and Φ̂ f are acceptable
and converge to zero for ηp → 0 and ηf → 0 [63]. Then, 3.4 Strong forms
the stationary conditions for the unconstrained optimization
problem in Eq. (35) with respect to the state variables, τ , r p To derive the strong forms of the governing equations, let
and τf − r f , yield the following evolution equations for the us consider the rate potential density per unit volume for the
plastic deformation and crack phase-field evolution: constitutive relations as
⎧    
⎪
⎪ ∂ Φ̂ p d
π Ċ, τ , r p , τf − r f = Ψ̂ (C) + V̂ Ċ, τ , r p , τf − r f
⎪
⎪ d p = λp ,
⎪
⎪ ∂τ dt
⎨ ∂ Φ̂ p (41)
ᾱ˙ = −λp p , (36)
⎪
⎪ ∂r
⎪
⎪ Φ̂ f ∂ Φ̂ f
⎪
⎪  , in terms of the basic constitutive functions Ψ̂ and V̂ defined
⎩ḋ =
ηf ∂ τf − r f in Eq. (11) and Eq. (35), respectively. The spatial integration

123
Computational Mechanics (2022) 69:151–175 159

of Eq. (41) yields the global rate potential such that [L] Accumulated strain:
   
∂ Φ̂ p
Π Ċ, τ , r p , τf − r f = π Ċ, τ , r p , τf − r f dV −ᾱ˙ − λp = 0, (50)
B0 ∂r p
− B · ϕ̇dV − T · ϕ̇dA, [L] Accumulated damage:
B0 ∂B0T
(42)
Φ̂ f ∂ Φ̂ f
ḋ −   = 0 with Φ̂ f = τf − r f . (51)
ηf ∂ τf − r f
where B and T are the body and traction forces in the ref-
erence configuration, respectively. Then, the involved field
Following Miehe et al. [45], we choose Eqs. (44)∼(46) as
variables for the continuum body are determined from the
the global governing equations for spatial descretization and
optimization problem of the above global rate potential as
the others are considered local ones. Then, a local return
  mapping algorithm is conducted to update internal vari-
˙ α̇, ḋ, Ḟ p , τ , r p , τf − r f
ϕ̇, ᾱ, ables by using Eqs. (47)∼(50), whose details are provided in
  “Appendix A”.
(43)
= arg inf  sup Π Ċ, τ , r p , τf − r f .
˙ α̇,ḋ, Ḟ
ϕ̇,ᾱ,
p
(τ ,r p ,τf −r f )

4 Numerical implementation within IGA

The stationary conditions for this optimization problem yield
framework
the strong forms of the governing equations with respect to
the eight fields as This section is devoted to the numerical implementation for
[G] Mechanical field:
the global governing equations in Eqs. (44)∼(46) within the
IGA framework.
d
· P + B = 0 with P = g (d) P + −
0 + P0 , (44)
dX
4.1 Weak forms
[G] Plastic field:
First, from the spatial integrations of Eqs. (3) and (44)
∂ Ψ̂ p d ∂ Ψ̂ p ∂ Ψ̂ f multiplied by the virtual displacement δu, the weak form
− · + = 0, (45) concerning the linear momentum is obtained as
∂α d X ∂∇α ∂α

[G] Phase-field: P : ∇δudV − B · δudV − T · δudA = 0.

B0 B0 ∂B0T

∂ Ψ̂ e ∂ Ψ̂ p ∂ Ψ̂ f d ∂ Ψ̂ f   (52)
+ + − · + τf − r f = 0, (46)
∂d ∂d ∂d d X ∂∇d
Also, the substitution of Eq. (20) into the spatial integration
[L] Hardening force: of Eq. (45) multiplied by the virtual quantity δα yields the
weak form with respect to micromorphic hardening
∂ Ψ̂ p    
− r p = 0, (47) y0 lp2 ∂G c α ∗
∂ ᾱ (α − ᾱ) δα + ∇α · ∇δα + γlf δα dV = 0,
B0 pp ∂α
[L] Plastic force: (53)

∂ Ψ̂ e e in which the Neumann condition in Eq. (5) has been reflected.

τ −2 · b = 0, (48) Similarly, by substituting Eqs. (21) and (51) into the spatial
∂ be
integration of Eq. (46) multiplied by the virtual quantity δd,
[L] Plastic strain: we have the weak form concerning the phase-field as

G c (α)  
∂ Φ̂ p ηf ḋδd + Ĥf dδd + lf2 ∇d · ∇δd
d −λ
p p
=0 B0 lf
∂τ 
 (49) (54)
3 2 p ∂ Ψ̂ e+ + Ψ̂ p
with λ = p
Φ̂ p and Φ̂ p = ||τ dev || − r , + δd dV = 0,
2ηp 3 ∂d

123
160 Computational Mechanics (2022) 69:151–175

Table 1 Material parameters for uni-directional uniform response

Parameter Value Unit

Young’s modulus E 200000 MPa

Poisson’s ratio ν 0.3 –
Elastic threshold Ψ̂cre 5.0 N mm
Initial yield stress y0 1000 MPa
Linear hardening parameter h 10000 MPa
Plastic viscosity ηp 10−7 MPa s
Plastic length scale parameter lp 2.0 mm
p
Plastic threshold Ψ̂cr 20/40/60 N mm
Fig. 4 Geometry and boundary conditions for uniform response
Penalty parameter pp 1000 MPa
Initial fracture toughness G c0 1000 kJ/m2
in consideration of the Neumann condition in Eq. (7). Here, Critical fracture toughness G c∞ 0.1 kJ/m2
Ĥf is the Heaviside step function, which has the value 1 if Saturation parameter βg 20/30/40/50 –
Φ̂ f has a positive value. Damage viscosity ηf 10−7 MPa s
Finally, the spatial counterparts of Eqs. (52)∼(54) are Crack length scale parameter lf 2.0 mm
obtained by the push-forward operations as, respectively,

σ : ∇ x δudv − b · δudv − t · δuda = 0, (55)

Bt Bt ∂Btt

1 y0 lp2
(α − ᾱ) δα + (∇ x α · F) · (∇ x δα · F)
Bt J pp
!
∂G c (α)
+ γlf δα dv = 0 (56)
∂α

and


1 G c (α)  
ηf ḋδd + Ĥf dδd + lf2 (∇ x d · F) · (∇ x δd · F)
Bt J lf

Fig. 5 Issue concerning the incorporation of elastoplastic material
∂ Ψ̂ e+ + Ψ̂ p
+ δd dv = 0, behavior into conventional crack phase-field modeling
∂d
(57)
functions R I as
where σ is the Cauchy stress. Here, ∇ denotes the spatial n#
node n#
node
gradient with respect to the reference configuration, whereas u= R I u I , δu = R I δu I ,
∇ x represents that for the current configuration. Also, b and I =1 I =1
t are the body and traction forces in the current configuration n#node n#node
defined as, respectively, α= R I α I , δα = R I δα I , (59)
I =1 I =1
n#node n#
node
B T d= R I d I , δd = R I δd I ,
b= and t = " . (58)
J J N · C −1 · N I =1 I =1

where the index n node represents the number of control points

in an element. Also, we approximate the rate quantity ḋ
4.2 Discretization within a finite time increment Δtn+1 = tn+1 − tn as

The global variables and their variations are approximated dn+1 − dn

ḋ = . (60)
in a general manner by using second-order NURBS basis Δtn+1

123
Computational Mechanics (2022) 69:151–175 161

The element residual vectors corresponding to the weak evaluated at the previous time step tn , so that Eq. (62) reads
forms (52), (53), (54) for the displacement, gradient enhanced
plasticity and crack phase-field are given as, respectively, ∂G c (αn+1 ) ∂G c (αn )
≈ = 0. (65)
  ∂α ∂α
∂R I
RϕI i = − σi j dv, (61)
Bt ∂x j The tangent matrices are simply obtained by differenti-
⎧⎛ ⎞ ating the above residual vectors with respect to the global
1 ⎨⎝ #
n node
RαI = − J ⎠ R I − ᾱ
R J αn+1 I variables. In this work, we introduce the staggered scheme
n+1 R
Bt J ⎩ J =1 for stable damage computations, as has been detailed in the
⎛ ⎞ ⎫ previous study [24]. Then, two tangent matrices
y0 lp2 n#
node
∂R J ∂R I ∂G (α) ⎬
+ ⎝ J ⎠F
αn+1 Fbk +
c
γlf R I dv
pp ∂ xa
ak
∂ xb ∂α ⎭ , 
J =1 Kϕϕ Kϕα  
Kep = and Kpf = Kdd (66)
(62) Kαϕ Kαα

and are considered for the displacement/plastic strain and crack

 n  phase-field, respectively. The components in Kep are written
1 ηf # node n#
node
in index notation as
RdI =− R J J
dn+1 − R J dnJ RI
Bt J Δtn+1  
J =1 J =1
n  ∂R I ∂R J 1 ∂τik
KϕI Ji ϕ j =
= G c (α) # node
Bt ∂ xk
Cik jl
∂ xl
dv with Cik jl =
J ∂ F ja
Fla − σil δ jk ,
+ R J J
dn+1 R I
lf (67)
J =1  
n  ∂ ᾱn+1

y0 lp2 ∂R I ∂R J
#node
∂R J J ∂R I IJ =
Kαα
1
1− RI R J + Fak Fbk dv,
+G c (α)lf dn+1 Fak Fbk Bt J ∂α pp ∂ x a ∂ xb
∂ xa ∂ xb (68)
J =1
+n node J J    
∂g J =1 R dn+1  e p  I ∂R I ∂σi j
+ Hn+1 + Hn+1 R dv. KϕI Ji α = R J dv, (69)
∂ x j ∂α
∂d 
Bt

J
(63) IJ =
Kαϕ
1
−R I ∂ ᾱn+1 ∂R F (70)
ki dv,
j
Bt J ∂ F ji ∂ xk

Here, by following the idea of [23], we have utilized the

characteristic history variables for the driving forces in Eq. where the details of Eqs. (67)∼(70) are provided in “Appendix
(63) as B”. On the other hand, the component in Kpf is written as

 1 ηf Gc I J
Ψ̂0,e+n+1 − Ψ̂cre if Ψ̂0,e+n+1 − Ψ̂cre > Hen Kdd
IJ
= RI R J + R R
Hen+1 = Bt J Δtn+1 lf
Hen otherwises
 ∂R I ∂R J
p
Ψ̂0, n+1 − Ψ̂cr
p
if
p
Ψ̂0, n+1 − Ψ̂cr > Hn
p p +G c lf Fak Fbk (71)
p
Hn+1 = ∂ xa ∂ xb
p
Hn otherwises !
∂2g  e p  I J
(64) + 2 Hn+1 + Hn+1 R R dv.
∂d

to ensure the irreversibility of damage evolution. Here, Ψ̂cre Note that a solver for unsymmetrical systems of equations is
p
and Ψ̂cr are material constants to control the activations of required to solve the set of fully coupled discretized equa-
elastic and plastic driving forces, respectively. tions, because the whole tangent matrix is unsymmetric.
A point to notice is that, different from standard gradient It is pertinent to note that the convergence of a staggered
enhanced plasticity, the derivative of the degrading fracture iterative process must be confirmed with an adequate crite-
toughness with respect to α appears in Eq. (62). Dittmann et rion. With reference to the convergence study conducted by
al. [50] employ the fracture toughness at the current time step [24], we employ the convergence criterion as
tn+1 in damage computations by ignoring this term. In this

regard, Yin and Kaliske [51] made the assumption such that Go to next loading step if log10 χk+1 ≤ Tolerance
∂G c /∂α = 0 by mentioning that the degrading tendency of ,
One more cycle if log10 χk+1 > Tolerance
the fracture toughness has a special non-linearity. Similarly,
we assume that the degrading fracture toughness could be (72)

123
162 Computational Mechanics (2022) 69:151–175

where χk+1 is the error at k + 1-th iteration defined as promote the damage evolution to some extent, according to
,   Eq. (2), such a deliberate setting is physically unacceptable
 Πk+1 − Πk  and useless for elastic-perfectly plastic materials.
χk+1 
= min   , |Πk+1 − Πk | . (73)
Πk+1  For this reason, we have introduced a plastic driving force
p
along with threshold Ψ̂cr for its activation. To demonstrate its
If the error becomes smaller than Tolerance, we go to the effect, neglecting the crack viscous parameter, we explicitly
next loading step. Otherwise, a new staggered iteration step derive the damage variable under uniform tensile deforma-
is restarted. Within each staggered iteration step, the Newton- tion as
Raphson iteration is carried out to make the residuals (52),
- . - .
(53), (54) small enough. p p
Ψ̂0e+ − Ψ̂cre + Ψ̂0 − Ψ̂cr
d= - . - .. (75)
p p
G c /2lf + Ψ̂0e+ − Ψ̂cre + Ψ̂0 − Ψ̂cr
5 Representative numerical examples
Here, the plastic driving force enables the phase-field variable
Three representative numerical examples are presented to to evolve depending on the amount of plastic deformation.
p
demonstrate the capability and performance of the pro- Also, by adjusting the critical value Ψ̂cr working as a thresh-
posed model. We utilize the second-order NURBS basis old, we can control the activation of the plastic driving force.
functions in our in-house code that ensures C 1 -continuity Such expected tendencies are exemplified in Fig. 6, which
at a minimum for the discretized weak forms. In order to reveals reduction of stiffness for both linear hardening and
save computation time, the convergence process of the stag- perfectly plastic cases with different threshold values. Nev-
gered iterative scheme is conducted only when dmax > 0.2. ertheless, since this modification realizes a gentle decrease
Also, to maintain computational stability, the evolution of only in all the cases, another strategy is needed to promote
plastic deformation is terminated at integration points with the failure.
sufficiently large damage values of d. Further details are Based on the conceptual background in Sect. 2 and the
explained in “Appendix C”. In the post-processes, we project numerical results presented above, we make the fracture
solutions from each control point to the corresponding FEM toughness be a degrading function. To be specific, the func-
node just for visualization purpose. tion form
   
5.1 Tensile failure of a homogeneous domain cos −βg α ∗ + 1
subjected to unidirectional uniform deformation G c (α) = (G c0 − G c∞ ) + G c∞ (76)
2

The first numerical example aims to explore the issue

is introduced, where G c0 and G c∞ are the initial and criti-
addressed in Sect. 2 through numerical analyses for unidi-
cal fracture toughness, respectively, and βg is the saturation
rectional uniform deformation by using the square element
parameter. Also, the independent variable is defined as α ∗ =
shown in Fig. 4. In this particular example, we employ the
α − αcr , where the accumulated plastic strain α is chosen in
linear hardening rule as
this study with αcr being a constant to control-the initiation
. of
p p
ŷ (ᾱ) = y0 + h ᾱ. (74) degradation. Here, αcr is the value of α when Ψ̂0 − Ψ̂cr first
has a non-zero value. Then, the damage variable is written as
The material parameters are given in Table 1. The two length - . - .
p p
scale parameters are set within acceptable ranges but have no Ψ̂0e+ − Ψ̂cre + Ψ̂0 − Ψ̂cr
influence on the numerical results since uniform deformation d= - . - .. (77)
p p
ensuring α = 0 and d = 0 is considered in this numerical G c (α) /2lf + Ψ̂0e+ − Ψ̂cre + Ψ̂0 − Ψ̂cr
examination.
Figure 5 shows the reaction force- Figure 7 shows the reaction force versus
displacement relations for the cases with h = 10000 MPa displacement for the cases with four different values of sat-
and h = 0 MPa, the latter of which corresponds to perfectly uration parameter βg . As can be seen, owing to the nature
plastic behavior. Here, the solid lines indicate the responses of the cosine function in Eq. (76), the fracture toughness
without damage for comparison. As can be seen from the monotonously decreases with the accumulation of plastic
figure, only slight damage evolution is observed for linear deformation, and as a result, a more rapid drop is obtained
plastic hardening case ΔΨ̂ e > 0, and no damage evolution with the increase in βg , so that the ductility can be controlled.
is observed for the elastic-perfectly plastic case ΔΨ̂ e = 0. Thus, the proposed model enables us to mimic the failure
Although a smaller value of fracture toughness can be used to phenomena of elastoplastic materials.

123
Computational Mechanics (2022) 69:151–175 163

Fig. 6 Reaction
force-displacement relations for
different critical plastic stored
p
energies, Ψ̂cr [N mm]

(a) (b)

Fig. 7 Reaction
force-displacement relations for
different saturation parameters,
βg [-]

(a) (b)

Fig. 8 Geometry, boundary

conditions and IGA
discretization for
asymmetrically notched
specimen

(a) (b)

5.2 Mixed-mode failure for asymmetrically notched of the specimen is fully fixed, and a forced displacement is
specimen imposed on the other edge. The material parameters are pro-
vided in Table 2. In this example, we adopt the standard Voce
The second numerical example is motivated to examine mesh hardening law
sensitivity and conduct parametric studies for the proposed
model. We target the asymmetrically notched specimen
shown in Fig. 8 with its boundary conditions. The bottom   
ŷ (ᾱ) = y0 + h ᾱ + (y∞ − y0 ) 1 − exp −βg ᾱ , (78)

123
164 Computational Mechanics (2022) 69:151–175

(a) (b)

(c)

Fig. 9 Crack evolution for asymmetrically notched specimen: mesh sensitivity study

(a) (b) (c)

Fig. 10 Damage distribution around the round notch

Table 2 Material parameters for asymmetrically notched specimen: and the same degrading fracture toughness Eq. (76) provided
mesh sensitivity study in the previous example.
Parameter Value Unit First, we investigate the relationship between crack length
scale parameter lf and mesh size h e . Fig. 9 shows the crack
Young’s modulus E 220,000 MPa
evolutions with three different mesh sizes of h e ≈ 0.80,
Poisson’s ratio ν 0.3 –
0.40, 0.20 mm with the crack length scale parameter being
Elastic threshold Ψ̂cre 0 N mm fixed at lf = 1.5 mm. For each case, the same failure pro-
Initial yield stress y0 810 MPa cess is observed. To be more specific, cracks initiate at both
Linear hardening parameter h 1500 MPa notches and, propagate towards the center of the specimen
Critical yield stress y∞ 1000 MPa and coalesce with each other to be a single crack path. How-
Saturation parameter βy 100 – ever, as can be seen in Fig. 10, which shows magnifications
Plastic viscosity ηp 10−6 MPa s of the cracked regions, the result of the coarsest mesh size
Plastic length scale parameter lp 2.0 mm h e ≈ 0.80 mm exhibits only a rough crack path. This is due to
Plastic threshold Ψ̂cr
p
0 N mm the inappropriateness of the magnitude relationship between
Penalty parameter pp 10000 MPa lf and h e in the crack phase-field modeling. In fact, according
Initial fracture toughness G c0 2000 kJ/m2 to the previous study [22], the inequality h e ≤ lf /2 is used
Critical fracture toughness G c∞ 0.00001 kJ/m2
as a rough indication to obtain appropriate crack topologies.
Saturation parameter βg 20 –
Nevertheless, the coarsest case does not satisfy the require-
ment; lf /2 = 0.75 mm < h e ≈ 0.80 mm. This defect is
Damage viscosity ηf 10−6 MPa s
also noticeable in the reaction force-displacement relation
Crack length scale parameter lf 1.5 mm
shown in Fig. 11a. The result is depicted with the solid line,
where the load fails to reach zero, since the poor connection

123
Computational Mechanics (2022) 69:151–175 165

Fig. 11 Reaction
force–displacement relations for
asymmetrically notched
specimen

(a) (b)

(c) (d)

(a) (b)
Fig. 12 Evolutions of degrading fracture toughness and constitutive work for the finest mesh size, h e ≈ 0.20 [mm]

between the damaged elements does not allow the formation a finer mesh size h e ≈ 0.20 mm that provides the appropri-
of a continuous crack path. ate crack path with completely connected damaged elements
To circumvent the unsatisfactory situation illustrated and the proper reaction force, as shown in Figs. 10c and 11a,
above, a finer mesh must be used. In this context, let us exam- respectively. It is, therefore, pertinent to point out that fur-
ine the result with the medium mesh size h e ≈ 0.40 mm ther attention should be paid for the relationship between lf
that fulfills the condition h e ≤ lf /2. However, as shown in and h e in terms of element topologies (triangle, rectangle
Fig.9b, the reaction force does not reach zero even when the etc.), the order of polynomials (linear, quadratic etc.), inter-
two cracks coalesce with each other. This is probably because polation scheme (FEM or IGA) and other factors. However,
the element alignment is not consistent with the crack path. this is beyond the scope of this study and will therefore be
As far as we know, the direction of elements is the same as the addressed in future work.
direction of crack propagation reported by [22] in their anal- Instead, we further investigate the result with the finest
ysis of the aforementioned condition. Although such a setup mesh size h e ≈ 0.20 mm in terms of the evolutions of
seems to be amenable to the condition, a more careful choice the degrading fracture toughness G c (α ∗ ) and the combined
is needed. This can probably be explored for the case with stored energy Ψ̂ e + Ψ̂ p shown in Fig. 12. It can be seen

123
166 Computational Mechanics (2022) 69:151–175

from Fig. 12a along with Fig. 9c that the simulated failure Table 3 Material parameters for asymmetrically notched specimen:
process seems close to reality. More specifically, the frac- parametric study
ture toughness is degraded due to the accumulation of plastic Parameter Value Unit
deformation even prior to damage evolution. In the damaged
Young’s modulus E 220000 MPa
regions representing crack paths, the fracture toughness tends
Poisson’s ratio ν 0.3 –
to become zero in accordance with the actual crack opening
area where no material is present. Similarly, as can be seen Elastic threshold Ψ̂cre 0 N mm
from Fig. 12b, the combined stored energy degrades along Initial yield stress y0 810 MPa
with damage evolution, and no energy exists in the fully dam- Linear hardening parameter h 1500 MPa
aged region. This is reasonable because no material rests on Critical yield stress y∞ 1000 MPa
the region between two discrete crack surfaces in a reality Saturation parameter βy 100 –
where stored energy cannot be defined. It should be recalled Plastic viscosity ηp 10−6 MPa s
that, with the intention of representing these phenomena, the Plastic length scale parameter lp 2.0/2.5/3.0 mm
degradation function is multiplied by not only the elastic Plastic threshold Ψ̂cr
p
0/10/20 N mm
stored energy, but also the plastic one in the proposed model; Penalty parameter pp 10000 MPa
see Eqs. (19) and (20) that were formulated based on the Initial fracture toughness G c0 2000 kJ/m2
discussion in Sect. 2. Critical fracture toughness G c∞ 0.00001 kJ/m2
Second, we conduct a parametric study to validate the Saturation parameter βg 10/15/20 –
capability and performance of the proposed model. The mate-
Damage viscosity ηf 10−6 MPa s
rial parameters for this study are provided in Table 3. Figure
Crack length scale parameter lf 1.5 mm
11b shows the reaction force versus displacement for three
different values of saturation parameter βg = 10, 15, 20.
As mentioned in the previous numerical example, this sat-
uration parameter controls the ductility. The lowest case of attention is paid to the relation between lf and lp although a
βg = 10 shows more ductile behavior and sustains larger variable investigation is provided by Miehe et al. [45].
enforced displacement than the other cases. Also, Fig. 13a–c Figures 11c and 13g–i show the reaction force
show the evolution of damage distributions. As can be seen -displacement relations and damage evolutions for three dif-
p
from these figures, the result of the case with βg = 10 tends ferent plastic threshold Ψ̂cr = 0, 10, 20 N mm, respectively.
to be diffusive, while the increase in βg makes the damage With the introduction of the plastic threshold, it is possi-
region less diffusive. This is because the smaller the value of ble to control the noticeable stiffness degradation. On the
βg is, the larger the plastic deformation is, and therefore its other hand, no remarkable difference in the damage evolu-
distribution is uniform in most cases. Owing to this tendency, tion is observed in these figures. This is probably due to the
fracture toughness is degraded more diffusively. attached notches that induce the crack paths. Thus, the point
Figure 11c shows the reaction force-displacement rela- that activates the plastic driving force does not affect the
tion for three different plastic length scale parameters lp = crack geometry. This is probably because the geometry of
2.0, 2.5, 3.0 mm and Fig. 13d–f show the corresponding the specimen is more dominant to determine the crack path
evolution of damage distributions. As can be seen from Fig. than the parameters for hardening behavior and damage evo-
11c, it is inferred that the increase in the plastic length scale lution. With respect to this point, we will make an additional
parameter lp delays the damage evolution. This is due to the note in the next subsection.
fact that the increase in the plastic length scale parameter Finally, from the results of these investigations, it can be
lp promotes the diffusive behavior of the global variable α concluded that the proposed model is promising to capture
that determines the decrease of the fracture toughness. This mixed-mode failure typical for elastoplastic materials. In par-
characteristic is reflected in Fig. 13d–f. To be specific, the ticular, the ductility and the onset of intense damage evolution
crack becomes wider as lp increases. Note that the origins of can be numerically calibrated by properly adjusting relevant
the gradient enhanced plasticity and the crack phase-field are parameters.
completely different despite the similarities in their govern-
ing equations. The crack phase-field stems from the diffusive 5.3 Tensile failure for flat I-shaped specimen
crack topology that provides a physical interpretation to crack
length scale parameter lf . In contrast, a conventional gradi- The last example is devoted to demonstrating the reproduc-
ent enhanced plasticity is a numerical approach based on the ing capability of crack patterns in a flat I-shaped specimen
basic idea of Taylor’s expansion for a local variable. Also, subjected to tensile loading. The target problem is borrowed
to our knowledge, the plastic length scale parameter itself from references [65,66], in which a typical crack pattern in
does not involve any physical background. Therefore, little shear-band for DP980 steel is reported. The geometry of the

123
Computational Mechanics (2022) 69:151–175 167

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)

p
Fig. 13 Crack evolution for asymmetrically notched specimen: parametric study with different values of βg , lp and Ψ̂cr

specimen is shown in Fig. 14 along with the boundary condi- characteristics, the Swift-Voce hardening law is adopted in
tions, in which the width and thickness are denoted by x and this particular numerical example as
y, respectively. The bottom of the specimen is fully fixed, and
a forced displacement is given on the top surface. In addition ŷ (ᾱ) = (1 − ζ ) ya (αb + ᾱ)βc
  (79)
to the original thin specimen (x = 12.5 mm, y = 1.2 mm) + ζ y∞ − (y∞ − y0 ) exp −β y ᾱ .
in [66], we also consider a thick specimen (x = 5.0 mm,
y = 3.0 mm) to examine the effect of the thickness on the The material parameters are provided in Table 4, in which the
crack patterns as reported by [67]. Here, the same cross- hardening parameters are taken from [65] and only the rate
sectional area A = 12.5 × 1.2 = 5 × 3 = 15 mm2 is set coefficient ζ is readjusted. Also, we assume the degrading
for both the specimens. To simulate the hardening behavior fracture toughness to be
relevant to the subsequent necking, localization and failure

123
168 Computational Mechanics (2022) 69:151–175

Fig. 15 Reaction force–displacement relations for flat I-shaped speci-

men (Thin (Exp.) [65])

Fig. 14 Geometry and boundary conditions for flat I-shaped specimen 

  
cos −βg1 α ∗ + 1
G c (α) = (G c0 − G c∞ )
2
Table 4 Material parameters for flat I-shaped specimen     (80)
cos −βg2 α − αcr2 + 1
Parameter Value Unit + G c∞ ,
2
Young’s modulus E 225800 MPa
Poisson’s ratio ν 0.39 – where βg2 is a saturation parameter. Here, the second bracket
Elastic threshold Ψ̂cre 0 N mm is appended to Eq. (76) to promote the evolution of degrada-
Initial yield stress y0 608.9 MPa tion of the fracture toughness with the critical accumulated
Critical yield stress y∞ 1071.8 MPa plastic strain αcr2 working as a threshold. Note that the func-
Saturation parameter βy 76.4 – tion above has been introduced for convenience, particularly
Strength coefficient ya 1573.3 MPa for the purpose of realizing fine calibrations in this numerical
Pre-strain parameter αb 0.000049 – example, although it should be determined by reference to
Hardening parameter βc 0.12936 – actually measured data. Additional experiment-based inves-
Rate coefficient ζ 0.4 – tigations are planned for future work.
Plastic viscosity ηp 10−3 MPa s
Figure 15 shows the reaction-force versus displacement
relations obtained from different thicknesses in comparison
Plastic length scale parameter lp 1.0 mm
p to the experimentally measured one (thin specimen only).
Plastic threshold Ψ̂cr 180 N mm
As can be seen in Fig. 15, both cases show almost the same
Penalty parameter pp 1000 MPa
response before necking; that is, the range of displacement is
Initial fracture toughness G c0 2000 kJ/m2 from u = 0 mm to around u ≈ 5 mm. This is because the two
Critical fracture toughness G c∞ 0.001 kJ/m2 specimens have the same cross-sectional area as mentioned
Saturation parameter βg1 5 – above. However, after they exhibit necking behavior, differ-
Saturation parameter βg2 200 – ent failure modes are observed, as shown in Figs. 16 and 19.
Plastic threshold αcr2 0.4 – These show detailed evolutions of accumulated plastic strain,
Damage viscosity ηf : 10−3 MPa s phase-field variable, degrading fracture toughness and con-
Crack length scale parameter lf 0.5 mm stitutive work density developed in the two specimens.
The numerical result for the thin specimen shown in
Fig. 15 matches the experimental one even after the neck-
ing and subsequent failed states. To be more specific, the
result shows a stable decrease of the reaction force within
u ≈ 5 ∼ 6 mm and an unstable one after u ≈ 6 mm. As can
be seen from Fig. 16, the specimen starts to show the necking
behavior around u ≈ 5.0 mm, and then forms a shear band

123
Computational Mechanics (2022) 69:151–175 169

Fig. 16 Detailed evolution of accumulated plastic strain α, phase-field variable d, fracture toughness G c and constitutive work Ψ̂ e + Ψ̂ p for thin
specimen with x = 12.5 [mm], y = 1.2 [mm]

Fig. 17 Experimental result for sively stretched region and finally makes a slanted crack path.
thin specimen subjected to
tensile loading: a fragment after Here, we remove the domain whose phase-field variable is
fracture reported in [66] greater than 0.95 in the post-processing for a better under-
standing of crack propagation. The simulated crack pattern is
in close agreement with both the experimental one reported
by [66] as shown in Fig. 17 and the other numerical result sim-
ulated by [65]. Also, the fracture toughness is degraded prior
to the damage evolution as in the previous numerical exam-
ple for the asymmetrically notched specimen. Similarly, the
combined stored energy degrades along with damage evo-
lution, and no energy exists in the fully damaged region.
This realizes the actual situation because the fully broken
region is incapable of storing energy. It is also worth men-
tioning that although the slanted crack has been obtained in
the present simulation, the horizontal crack reported typi-
Fig. 18 Various crack patterns for the thin specimen obtained during cally in the previous studies [43,51] can also be simulated by
calibration process adjusting function forms of the hardening and degradation
functions as well as their parameters. This is because, unlike
the previous example, the crack path is less affected by the
around u ≈ 5.5 mm. After that, strain localization appears geometry of the flat I-shaped specimen than the aforemen-
around the center of the specimen, and the plastic deforma- tioned various factors that determine the evolution of plastic
tion further localized selectively to one of the shear bands. deformation and damage. In fact, Fig. 18 shows several crack
Subsequently, an asymmetric evolution of accumulated plas- patterns obtained during our calibration process for Figs. 15
tic strain is obtained. Following this strain localization, the and 19. Thus, these investigations using DP980 steel are a
crack initiates from the center, propagates along the inten-

123
170 Computational Mechanics (2022) 69:151–175

Fig. 19 Detailed evolution of accumulated plastic strain α, phase-field variable d, fracture toughness G c and constitutive work Ψ̂ e + Ψ̂ p for thick
specimen x = 5 [mm], y = 3 [mm]

Fig. 20 Experiment result for ter of the specimen and spreads out toward the four edges on
thick specimen subjected to
tensile loading [67] the horizontal plane, in which the plastic strain is accumu-
lated. Again, the crack pattern conforms to the experimental
one provided in Fig. 20, although neither the dimensions nor
the material properties are the same as what we used here.
This kind of difference in crack pattern was referred to as the
thickness effect [67]. Specifically, with the x/y rate increas-
ing to 1, the effect of shear stress disappears and results in
the strain localization occurring in the horizontal direction
instead of the shear band, which causes the crack surface
clear demonstration of the capability in adequately simulat- perpendicular to the loading direction. Thus, the proposed
ing typical failure processes in metallic materials. model has successfully demonstrated the thickness effect.
On the other hand, the reaction force-displacement rela- From these results, it is reasonable to conclude that the
tion for the thick specimen shown in Fig. 15 is similar to that proposed model has the capabilities to realize various failure
of the thin one and also exhibits typical elastic-plastic mate- patterns that are often experimentally measured for elasto-
rial behavior, in which a stable reduction of material stiffness plastic materials. This can be confirmed by further studies to
is observed for u ≈ 5 ∼ 5.8 mm until the unstable failure gain similar perspectives in experimental investigations.
stage. Nevertheless, the crack pattern in Fig. 19 is different
from that of the thin specimen; while necking is observed,
it does not form a shear band. Instead, strain localization
develops horizontally. This tendency continues in the crack
evolution; that is, the critical damage initiates from the cen-

123
Computational Mechanics (2022) 69:151–175 171

6 Conclusion and flow tensor n are considered as internal variables to be

determined. Also, to ensure incompressibility of the plastic
Crack phase-field modeling is a promising approach to deformation, an exponential mapping is utilized.
capture failure characteristics. In this study, a novel phase- The residuals of two local differential equations yield
field model is proposed for ductile fracture by introducing   
both the plastic driving force and the degrading fracture 2 p 2ηp p
rkλ
p
= − ||τ 0dev,k || − r − λ (A.1)
toughness into crack phase-field computations based on the 3 0,k 3g (d) k
phenomenological justification for ductile fracture.  
τ 0dev,k
The governing equations in the proposed model are dis- rnk =− − nk , (A.2)
cretized within the IGA framework, and three representative ||τ 0dev,k ||
numerical examples are studied to demonstrate its capability
where the index k represents the iteration of the local Newton-
and performance. First, the tensile failure simulations with a
Raphson scheme. By numerical linearization, its tangent
single element reveals the underlying issue of the straightfor-
matrix yields
ward combined application of the elastoplastic constitutive
law and conventional crack phase-field modeling. Then, by , 
kλp λp kλp n
introducing the plastic driving force and the degrading frac- k= (A.3)
knλp knn
ture toughness, we have shown that the model is capable of
representing the failure phenomena in typical elastoplastic with diagonal components
materials. The objective of the second numerical example is
  p  
to conduct studies on the sensitivities with respect to the ∂||τ 0dev,k || ∂ exp −2Δtn+1 λk nk · be,tr
n+1
mesh size and the parameters relevant to the failure pro- kλp λp = :
∂ be ∂λp
cess in an asymmetrically notched specimen whose material  p
exhibits elastic-plastic material behavior. From the mesh sen- 2 ∂r0,k 2ηp
− −
sitivity study, we realized the necessity of paying attention 3 ∂λ p 3g (d)
to the relationship between lf and h e in addtion to the gen- ∂ exp (−2Δtn+1 λp n) · be,tr
n+1
eral requirement h e ≤ lf /2. In addition, a parametric study with
∂λp
is made to confirm the ability to control the ductility and   e,tr
= D exp −2Δtn+1 λp n : (−2Δtn+1 n) · bn+1
the onset of intense damage evolution, which is of particu-
lar concern in fine numerical calibrations. The last example (A.4)
is devoted to demonstrating the capable in reproducing the
and
crack patterns often observed in actual experiments. The
typical failure along the shear band is simulated for a thin 1 ∂τ 0dev,k τ 0dev,k ∂||τ 0dev,k ||
knn = − ⊗ − 1sym.
specimen due to strain localization, whereas a thick bar has a ||τ 0dev,k || ∂n ||τ 0dev,k ||2 ∂n
horizontal crack path typical of the so-called thickness effect.  
1 ∂τ 0dev,k τ 0dev,k ∂||τ 0dev,k ||
Additionally, we have confirmed that the proposed model = − ⊗
||τ 0dev,k || ∂ be ||τ 0dev,k ||2 ∂ be
enables us to obtain different crack patterns by changing the  
p
activation conditions of the plastic driving force and degrad- ∂ exp −2Δtn+1 λk nk · be,tr n+1
ing fracture toughness. : − 1sym.
∂n
Nevertheless, further mesh sensitivity and parametric  
∂ exp −2Δtn+1 λp n · be,tr n+1
studies need to be made for engineering applications involv- with
 ∂n
ing more complex geometries and material behavior. Besides,    
= D exp −2Δtn+1 λp n : −2Δtn+1 λp 1sym. ∗ be,tr
this work is formulated only for the quasi-static case, and  
n+1

needs to be extended to a dynamic setup. Further experi- iakl aj
mental investigations into the variation of fracture toughness (A.5)
during plastic deformation will also be conducted in future
studies. and coupling components
 p  e,tr
∂||τ 0dev,k || ∂ exp −2Δtn+1 λk nk · bn+1
kλp n = : (A.6)
A Return mapping algorithm for internal ∂ be ∂n
variables and
We show a pseudo-code for the local return mapping algo- 1 ∂τ 0dev,k τ 0dev,k ∂||τ 0dev,k ||
knλp = −
rithm in Fig. 21. During this process, the plastic multiplier λp ||τ 0dev,k || ∂λp ||τ 0dev,k ||2 ∂λp

123
172 Computational Mechanics (2022) 69:151–175

Fig. 21 Return mapping algorithm

 
1 ∂τ 0dev,k τ 0dev,k ∂||τ 0dev,k || with
= − ⊗
||τ 0dev,k || ∂ be ||τ 0dev,k ||2 ∂ be
⎧   ⎫
⎨ ∂ exp −2Δtn+1 λpk nk ⎬  e

: · e,tr ∂τ 0 ∂τ 0vol ∂ Je ∂τ 0dev ∂ b̄ ∂ Fe
⎩ ∂λp
bn+1 ⎭ . (A.7) = ⊗ + e : ∂ Fe :
∂F ∂ Je ∂ Fe ∂ b̄ ∂F

= κ J e2 1 ⊗ F e−T + μJ e−2/3 (B.2)

, !
2 ∂ Fe
Here, D exp denotes the derivative of the second order expo- dev 1 ⊗ F e + F e ⊗ 1 − be ⊗ F e−T : .
3 ∂F
nential tensor and 1sym. is the fourth order symmetric identity
tensor. Also, to show formulations explicitly in tensor nota-
tions, we define an operator ∗ that denotes an inner product Using the total differential, the last term in Eq. (B.2) yields
of the second basis of a fourth order tensor X and the first
basis of a second order tensor Y , which is shown in index
 
notation as X ia jk Ya j . ∂ Fe ∂ F e  ∂ F e  ∂λp
= + ⊗
∂F ∂ F λp ,n=const. ∂λp n,F=const. ∂ F
 (B.3)
∂ F e  ∂n
+  :
∂n p ∂F
F,λ =const.

B Components accounting for global

tangent matrix with

The component in Eq. (67) is written in an explicit form as


∂ F e    p−T
= exp −Δtn+1 λp n ⊗ F n , (B.4)
∂ F λp ,n=const.

⎧ ∂ F e   
= D exp −Δtn+1 λp n : (−Δtn+1 n) · F e,tr , (B.5)
⎪ ∂τ ∂τ
⎨g (d) 0vol + g (d) 0dev for J e ≥ 1 ∂λp n,F=const.
∂τ ∂ ∂F
    
= F ∂ F e  
∂F ⎪ ∂τ 0vol ∂τ 0dev  = D exp −Δtn+1 λp n : −Δtn+1 λp 1sym. ∗ F e,tr  .
⎩ + g (d) for J e < 1 ∂n F,λp =const.  
∂F ∂F iakl
aj
(B.1) (B.6)

123
Computational Mechanics (2022) 69:151–175 173

∂λp ∂n with
Then, and are obtained from the total differential /  0
∂F ∂F 2ηp p
of two local residuals λp 2 p
dα r = dα −||τ 0dev,n+1 || + r + λ =0
3 0,n+1 3g (d) n+1
/  0
∂rλ ∂rλ ∂λp ∂rλ
p p p
2 p 2ηp p ∂n
d F rλ = d F −||τ 0dev,n+1 || +
p
r + λ =0 → + + :
3 0,n+1 3g (d) n+1 ∂α ∂λp ∂α ∂n ∂α
 p
∂rλ
p
∂rλ ∂λp
p
∂rλ ∂n
p 2 ∂r0 ∂λ p ∂n
→ + + : = − kλp λp − kλp n : =0
∂F ∂λp ∂ F ∂n ∂ F 3 ∂α ∂α ∂α
(B.12)
∂||τ 0dev || ∂λp ∂n
=− − kλp λp − kλp n : =0 (B.7)
∂F ∂F ∂F
and
, 
τ 0dev,n+1
and dα rn = dα − +n =0
||τ 0dev,n+1 ||
, 
(B.13)
τ 0dev,n+1 ∂rn ∂λp ∂rn ∂n ∂λp ∂n
n
dF r = dF − +n =0 → + : = −knλp − knn : = 0,
||τ 0dev,n+1 || ∂λp ∂α ∂n ∂α ∂α ∂α
∂rn ∂rn ∂λp ∂rn ∂n Taking into account of the symmetry of flow tensor, an equation in Vogit
→ + p ⊗ + : (B.8)
∂F ∂λ ∂F ∂n ∂ F notation for Eqs. (B.12) and (B.13) is obtained
∂ τ 0dev ∂λp ∂n
=− − knλp ⊗ − knn : = 0. ⎡ ⎤
∂ F ||τ 0dev || ∂F ∂F ∂λp
⎡ ⎤ ⎢ ∂α ⎥
kλp λp kλp n kλp n kλp n 2kλp n 2kλp n 2kλp n ⎢ ∂n11 ⎥
11 22 33 23 13 12 ⎢ ⎥
⎢ k ⎥ ⎢ ⎥
⎢ n λp kn11 n11 kn11 n22 kn11 n33 2kn11 n23 2kn11 n13 2kn11 n12
⎥ ⎢ ∂α ⎥
Taking into account the symmetry of flow tensor, an equation in Vogit ⎢ 11 ⎥ ⎢ ∂n22 ⎥
⎢ k kn22 n11 kn22 n22 kn22 n33 2kn22 n23 ⎥
2kn22 n13 2kn22 n12 ⎢ ⎥
notation for Eqs. (B.7) and (B.8) is obtained ⎢ n22 λp
⎢
⎥
⎥
⎢
⎢ ∂α ⎥
⎥
⎢ k kn33 n11 kn33 n22 kn33 n33 2kn33 n23 ⎥
2kn33 n13 2kn33 n12 ⎢ ∂n33 ⎥
⎢ n33 λp ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ∂α ⎥
⎡ ⎤ ⎢ 2kn λp 2kn23 n11 2kn23 n22 2kn23 n33 4kn23 n23 4kn23 n13 4kn23 n12 ⎥ ⎢ ∂n23 ⎥
∂λp ⎢ 23 ⎥ ⎢ ⎥
⎢ 2k 2kn13 n11 ⎥
2kn13 n22 2kn13 n33 4kn13 n23 4kn13 n13 4kn13 n12 ⎦ ⎢ ∂α ⎥
⎡ ⎤ ⎢ ⎥ ⎣ n13 λp ⎢ ∂n13 ⎥
kλp λp kλp n kλp n kλp n 2kλp n 2kλp n 2kλp n ⎢ ∂F ⎥ ⎢ ⎥
⎢ ∂n11 ⎥ 2kn λp 2kn12 n11 2kn12 n22 2kn12 n33 4kn12 n23 4kn12 n13 4kn12 n12 ⎢ ∂α ⎥
⎢ k
11 22 33 23 13
⎥
12
⎢ ⎥ ⎣ ⎦
⎢ n λp kn11 n11 kn11 n22 kn11 n33 2kn11 n23 2kn11 n13 2kn11 n12
⎥ ⎢ ∂F ⎥
12
  ∂n12
⎢ 11 ⎥ ⎢ ∂n22 ⎥ ∂α
⎢ k kn22 n11 kn22 n22 kn22 n33 2kn22 n23 ⎥
2kn22 n13 2kn22 n12 ⎢ ⎥ k
⎢ n22 λp ⎥ ⎢ ∂F ⎥
⎢ ⎥ ⎢ ∂n33 ⎥ ⎡ p ⎤
⎢ k kn33 n11 kn33 n22 kn33 n33 2kn33 n23 ⎥
2kn33 n13 2kn33 n12 ⎢ ⎥ 2 ∂r0
⎢ n33 λp ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ∂F ⎥ ⎢ ⎥
⎢ 2kn λp 2kn23 n11 2kn23 n22 2kn23 n33 4kn23 n23 4kn23 n13 4kn23 n12 ⎥ ⎢ ∂n23 ⎥ ⎢ 3 ∂α ⎥
⎢ 23 ⎥ ⎢ ⎥ ⎢ ⎥
⎢ 2k 2kn13 n11 ⎥
2kn13 n22 2kn13 n33 4kn13 n23 4kn13 n13 4kn13 n12 ⎦ ⎢ ∂F ⎥ ⎢ 0 ⎥
⎣ n13 λp ⎢ ∂n13 ⎥ ⎢ ⎥
⎢ ⎥ ⎢ 0 ⎥
2kn λp 2kn12 n11 2kn12 n22 2kn12 n33 4kn12 n23 4kn12 n13 4kn12 n12 ⎢ ⎥ =⎢ ⎥
⎣ ∂F ⎦ ⎢ 0 ⎥
12
  ∂n12 ⎢ ⎥
⎢ 0 ⎥
∂F ⎢ ⎥
k ⎣ 0 ⎦
⎡ ∂||τ 0dev || ⎤
0
⎢ ∂F ⎥
⎢ ∂ τ0dev,11
⎢
⎥
⎥
(B.14)
⎢ ∂ F ||τ ⎥
⎢ 0dev || ⎥
⎢ ∂ τ0dev,22 ⎥
⎢ ⎥ with
⎢ ∂ F ||τ ⎥
⎢ 0dev || ⎥
⎢ ∂ τ0dev,33 ⎥
⎢ ⎥
= − ⎢ ∂ F ||τ
⎢ 0dev ||
⎥
⎥
(B.9) p
∂r0 ∂
⎢ ∂ τ0dev,23
⎢2
⎥
⎥ = ŷ (ᾱ) + pp (ᾱ − α) = − pp . (B.15)
⎢ ∂ F ||τ ⎥ ∂α ∂α
⎢
⎢ ∂ τ0dev,130dev || ⎥
⎥
⎢2 ⎥
⎢ ∂ F ||τ ⎥
⎢
⎣ ∂ τ0dev,120dev || ⎥
⎦
∂ ᾱ ∂ ᾱ
In addition, and are obtained as
2
∂ F ||τ 0dev || ∂α ∂F
   
∂ ᾱ ∂ 2 p 2 ∂λp
with = ᾱn + λ Δtn+1 = Δtn+1 (B.16)
∂α ∂α 3 3 ∂α
∂||τ 0dev || ∂||τ 0dev ||   p−T
= : exp −Δtn+1 λp n ⊗ F n , and
∂ F n+1 ∂ F en+1
     
∂ τ 0dev 1 ∂τ 0dev τ 0dev ∂||τ 0dev || ∂ ᾱ ∂ 2 p 2 ∂λp
= − ⊗ = ᾱn + λ Δtn+1 = Δtn+1 . (B.17)
∂ F n+1 ||τ 0dev || ||τ 0dev || ∂ F en+1 ||τ 0dev ||2 ∂ F en+1 ∂F ∂F 3 3 ∂F
  p−T
: exp −Δtn+1 λp n ⊗ F n .
(B.10)
C An issue about crack phase-field modeling
∂ Fe
In a same manner, yields Generally, once the material experiences plastic deformation, it does
∂α
not show elastic behavior again except unloadings. Nevertheless, due
  to the setup of damage modeling incorporated into the elastoplastic
∂ Fe ∂ F e  ∂λp ∂ F e  ∂n
= + : (B.11) constitutive model such as [42,50,51], the material shows elastic behav-
∂α ∂λp n=const. ∂α ∂n λp =const. ∂α ior again when the damage accumulates significantly. We explain this

123
174 Computational Mechanics (2022) 69:151–175

elastic-plastic-elastic transition using the same setup as in Sect. 2. We 5. Moës N, Dolbow J, Belytschko T (1999) A finite element method
assume the material is under the plastic state and damaged by external for crack growth without remeshing. Int J Numer Methods Eng
loadings. Then, the damage variable without the plastic driving force at 46(1):131–150
the end of the loading step n yields 6. Moës N, Belytschko T (2002) Extended finite element method for
cohesive crack growth. Eng Fract Mech 69(7):813–833
7. Wells GN, Sluys L (2001) A new method for modelling cohe-
Ψ̂ne
dn = sive cracks using finite elements. Int J Numer Methods Eng
G c /2lf + Ψ̂ne (C.1) 50(12):2667–2682
1 p un 8. Terada K, Asai M, Yamagishi M (2003) Finite cover method for
with Ψ̂ne = Eεne 2 , εne = εn − εn and εn = , linear and non-linear analyses of heterogeneous solids. Int J Numer
2 L + un
Methods Eng 58(9):1321–1346
where L and u n denote the structure’s length and the total displacement 9. Gurson AL (1977) Continuum theory of ductile rupture by void
nucleation and growth: part I-yield criteria and flow rules for porous
by the loading step n. When the next loading is introduced at the loading
ductile media. J Eng Mater Technol 99(1):2–15
step n + 1, we start from checking the deformation state by using the
p 10. Tvergaard V, Needleman A (1984) Analysis of the cup-cone fracture
yield function Φ̂n+1,k=1 at the first global Newton-Raphson iteration
in a round tensile bar. Acta Metall 32(1):157–169
k = 1. By using the variables obtained in the previous loading step n, 11. Lemaitre J (1985) A continuous damage mechanics model for duc-
the yield function is written as tile fracture. J Eng Mater Technol 107:83–89
 12. Bažant ZP, Prat PC (1988) Microplane model for brittle-plastic
p p
Φ̂n+1,k=1 = (1 − dn )2 σn+1,k=1 − σ0 + hεn material: I. Theory. J Eng Mech 114(10):1672–1688
 p  13. Bazant ZP, Belytschko TB, Chang TP et al (1984) Continuum theory
with σn+1,k=1 = E εn+1 − εn and εn+1 (C.2) for strain-softening. J Eng Mech 110(12):1666–1692
u n + Δu n+1 14. Pijaudier-Cabot G, Bažant ZP (1987) Nonlocal damage theory. J
= ,
L + u n + Δu n+1 Eng Mech 113(10):1512–1533
15. Bažant ZP, Pijaudier-Cabot G (1988) Nonlocal continuum damage,
where Δu n+1 is the new displacement increment. Here, Eq. (C.2) may localization instability and convergence. J Appl Mech 55(2):287–
be lower than zero due to the updated damage variable dn , which repre- 293
sents that the material shows elastic behavior despite new loading. By 16. Peerlings RH, de Borst R, Brekelmans WM, De Vree J (1996a)
conducting some algebra, we obtain the explicit displacement increment Gradient enhanced damage for quasi-brittle materials. Int J Numer
Methods Eng 39(19):3391–3403
 p 17. Peerlings Rd, Borst Rd, Brekelmans Wd, Vree Jd, Spee I (1996b)
Cn+1 (L + u n ) − u n σ0 + hεn p Some observations on localization in non-local and gradient dam-
Δu n+1 ≥ with Cn+1 = + εn ,
Cn+1 − 1 (1 − dn )2 E age models. Eur J Mech A Solids 15(6):937–953
(C.3) 18. Francfort GA, Marigo JJ (1998) Revisiting brittle fracture as an
energy minimization problem. J Mech Phys Solids 46(8):1319–
which leads the material into the elastic state again. Thus, we may infer 1342
that this setup is not absolutely appropriate to explain ductile fracture 19. Mumford DB, Shah J (1989) Optimal approximations by piecewise
in elastoplastic materials. On the other hand, we also realize that this smooth functions and associated variational problems. Commun
unrealistic transition of deformation state contributes to computational Pure Appl Math
20. Ambrosio L, Tortorelli VM (1990) Approximation of functional
stability. As widely known, the violent plastic deformation due to the
depending on jumps by elliptic functional via t-convergence. Com-
severely damaged mesh’s softening usually causes unstable tendencies,
mun Pure Appl Math 43(8):999–1036
especially in the local return mapping process. This problem may be
21. Bourdin B, Francfort GA, Marigo JJ (2000) Numerical experiments
avoided by the unreal state transition since the severely damaged mesh in revisited brittle fracture. J Mech Phys Solids 48(4):797–826
does not show the plastic behavior. 22. Miehe C, Welschinger F, Hofacker M (2010a) Thermodynamically
However, the proposed model introduces the plastic deformation, consistent phase-field models of fracture: variational principles
which results in the inevitable of the violent plastic deformation in and multi-field FE implementations. Int J Numer Methods Eng
the severely damaged region. Therefore, we set a threshold around 83(10):1273–1311
d ≈ 0.9 ∼ 0.95 to terminate the plastic deformation to maintain com- 23. Miehe C, Hofacker M, Welschinger F (2010b) A phase field model
putational stability. for rate-independent crack propagation: robust algorithmic imple-
mentation based on operator splits. Comput Methods Appl Mech
Eng 199(45–48):2765–2778
24. Ambati M, Gerasimov T, De Lorenzis L (2015) A review on phase-
References field models of brittle fracture and a new fast hybrid formulation.
Comput Mech 55(2):383–405
1. Griffith AA (1921) Vi. The phenomena of rupture and flow in solids. 25. Wu JY (2017) A unified phase-field theory for the mechanics of
Philos Trans R Soc Lond Ser A 221(582–593):163–198 damage and quasi-brittle failure. J Mech Phys Solids 103:72–99
2. Ngo D, Scordelis AC (1967) Finite element analysis of reinforced 26. Wu JY, Nguyen VP (2018) A length scale insensitive phase-field
concrete beams. J Proc 64:152–163 damage model for brittle fracture. J Mech Phys Solids 119:20–42
3. Ingraffea AR, Saouma V (1985) Numerical modeling of dis- 27. Henry H, Levine H (2004) Dynamic instabilities of fracture
crete crack propagation in reinforced and plain concrete. Fracture under biaxial strain using a phase field model. Phys Rev Lett
mechanics of concrete: structural application and numerical calcu- 93(10):105504
lation. Springer, Netherlands, pp 171–225 28. Amor H, Marigo JJ, Maurini C (2009) Regularized formulation
4. Jirásek M (2000) Comparative study on finite elements with embed- of the variational brittle fracture with unilateral contact: numerical
ded discontinuities. Comput Methods Appl Mech Eng 188(1– experiments. J Mech Phys Solids 57(8):1209–1229
3):307–330

123
Computational Mechanics (2022) 69:151–175 175

29. Steinke C, Kaliske M (2019) A phase-field crack model based on 50. Dittmann M, Aldakheel F, Schulte J, Wriggers P, Hesch C (2018)
directional stress decomposition. Comput Mech 63(5):1019–1046 Variational phase-field formulation of non-linear ductile fracture.
30. Storm J, Supriatna D, Kaliske M (2020) The concept of represen- Comput Methods Appl Mech Eng 342:71–94
tative crack elements for phase-field fracture: anisotropic elasticity 51. Yin B, Kaliske M (2020) A ductile phase-field model based on
and thermo-elasticity. Int J Numer Methods Eng 121(5):779–805 degrading the fracture toughness: theory and implementation at
31. Borden MJ, Hughes TJ, Landis CM, Verhoosel CV (2014) A small strain. Comput Methods Appl Mech Eng 366:113068
higher-order phase-field model for brittle fracture: Formulation 52. Kachanov L (1986) Introduction to continuum damage mechanics,
and analysis within the isogeometric analysis framework. Comput vol 10. Springer, Berlin
Methods Appl Mech Eng 273:100–118 53. Murakami S (2012) Continuum damage mechanics: a continuum
32. Cottrell JA, Hughes TJ, Bazilevs Y (2009) Isogeometric analysis: mechanics approach to the analysis of damage and fracture, vol
toward integration of CAD and FEA. Wiley, Hoboken 185. Springer, Berlin
33. Kakouris E, Triantafyllou SP (2017) Phase-field material 54. Forest S (2009) Micromorphic approach for gradient elasticity, vis-
point method for brittle fracture. Int J Numer Methods Eng coplasticity, and damage. J Eng Mech 135(3):117–131
112(12):1750–1776 55. Miehe C, Teichtmeister S, Aldakheel F (2016) Phase-field mod-
34. Roy P, Pathrikar A, Deepu S, Roy D (2017) Peridynamics damage elling of ductile fracture: a variational gradient-extended plasticity-
model through phase field theory. Int J Mech Sci 128:181–193 damage theory and its micromorphic regularization. Philos Trans
35. Larsen CJ, Ortner C, Süli E (2010) Existence of solutions to a R Soc A Math Phys Eng Sci 374(2066):20150170
regularized model of dynamic fracture. Math Models Methods Appl 56. Miehe C, Aldakheel F, Teichtmeister S (2017) Phase-field model-
Sci 20(07):1021–1048 ing of ductile fracture at finite strains: a robust variational-based
36. Bourdin B, Larsen CJ, Richardson CL (2011) A time-discrete model numerical implementation of a gradient-extended theory by micro-
for dynamic fracture based on crack regularization. Int J Fract morphic regularization. Int J Numer Methods Eng 111(9):816–863
168(2):133–143 57. Saanouni K, Hamed M (2013) Micromorphic approach for finite
37. Borden MJ, Verhoosel CV, Scott MA, Hughes TJ, Landis CM gradient-elastoplasticity fully coupled with ductile damage: for-
(2012) A phase-field description of dynamic brittle fracture. Com- mulation and computational aspects. Int J Solids Struct 50(14–
put Methods Appl Mech Eng 217:77–95 15):2289–2309
38. Hofacker M, Miehe C (2012) Continuum phase field modeling of 58. Nguyen VD, Lani F, Pardoen T, Morelle X, Noels L (2016) A large
dynamic fracture: variational principles and staggered FE imple- strain hyperelastic viscoelastic–viscoplastic-damage constitutive
mentation. Int J Fract 178(1–2):113–129 model based on a multi-mechanism non-local damage continuum
39. Sargado JM, Keilegavlen E, Berre I, Nordbotten JM (2018) High- for amorphous glassy polymers. Int J Solids Struct 96:192–216
accuracy phase-field models for brittle fracture based on a new 59. Diamantopoulou E, Liu W, Labergere C, Badreddine H, Saanouni
family of degradation functions. J Mech Phys Solids 111:458–489 K, Hu P (2017) Micromorphic constitutive equations with damage
40. Zhang X, Vignes C, Sloan SW, Sheng D (2017) Numerical evalu- applied to metal forming. Int J Damage Mech 26(2):314–339
ation of the phase-field model for brittle fracture with emphasis on 60. Brepols T, Wulfinghoff S, Reese S (2020) A gradient-extended two-
the length scale. Comput Mech 59(5):737–752 surface damage-plasticity model for large deformations. Int J Plast
41. Tanné E, Li T, Bourdin B, Marigo JJ, Maurini C (2018) Crack 129:102635
nucleation in variational phase-field models of brittle fracture. J 61. Kuhn C, Schlüter A, Müller R (2015) On degradation functions in
Mech Phys Solids 110:80–99 phase field fracture models. Comput Mater Sci 108:374–384
42. Ambati M, Gerasimov T, De Lorenzis L (2015) Phase-field model- 62. Yin B, Steinke C, Kaliske M (2020) Formulation and implemen-
ing of ductile fracture. Comput Mech 55(5):1017–1040 tation of strain rate-dependent fracture toughness in context of the
43. Ambati M, Kruse R, De Lorenzis L (2016) A phase-field model for phase-field method. Int J Numer Methods Eng 121(2):233–255
ductile fracture at finite strains and its experimental verification. 63. Simo JC, Hughes TJ (2006) Computational inelasticity, vol 7.
Comput Mech 57(1):149–167 Springer, New York
44. Miehe C, Hofacker M, Schänzel LM, Aldakheel F (2015) Phase 64. Matsubara S, Terada K (2021) A variationally consistent formu-
field modeling of fracture in multi-physics problems. Part ii. Cou- lation of the thermo-mechanically coupled problem with non-
pled brittle-to-ductile failure criteria and crack propagation in associative viscoplasticity for glassy amorphous polymers. Int J
thermo-elastic–plastic solids. Comput Methods Appl Mech Eng Solids Struct 212:152–168
294:486–522 65. Djouabi M, Ati A, Manach PY (2019) Identification strategy
45. Miehe C, Aldakheel F, Raina A (2016) Phase field modeling of influence of elastoplastic behavior law parameters on Gurson–
ductile fracture at finite strains: a variational gradient-extended Tvergaard–Needleman damage parameters: Application to DP980
plasticity-damage theory. Int J Plast 84:1–32 steel. Int J Damage Mech 28(3):427–454
46. Miehe C (2014) Variational gradient plasticity at finite strains. 66. Chung K, Lee C, Kim H (2014) Forming limit criterion for ductile
Part i: mixed potentials for the evolution and update problems of anisotropic sheets as a material property and its deformation path
gradient-extended dissipative solids. Comput Methods Appl Mech insensitivity, part ii: boundary value problems. Int J Plast 58:35–65
Eng 268:677–703 67. Ikeda K, Okazawa S, Terada K, Noguchi H, Usami T (2001)
47. Miehe C, Welschinger F, Aldakheel F (2014) Variational gradient Recursive bifurcation of tensile steel specimens. Int J Eng Sci
plasticity at finite strains. Part ii: Local-global updates and mixed 39(17):1913–1934
finite elements for additive plasticity in the logarithmic strain space.
Comput Methods Appl Mech Eng 268:704–734
48. Borden MJ, Hughes TJ, Landis CM, Anvari A, Lee IJ (2016)
Publisher’s Note Springer Nature remains neutral with regard to juris-
A phase-field formulation for fracture in ductile materials: finite
dictional claims in published maps and institutional affiliations.
deformation balance law derivation, plastic degradation, and stress
triaxiality effects. Comput Methods Appl Mech Eng 312:130–166
49. Alessi R, Ambati M, Gerasimov T, Vidoli S, De Lorenzis L (2018)
Comparison of phase-field models of fracture coupled with plastic-
ity. In: Advances in computational plasticity. Springer, pp 1–21

123