# Coupling of the phase field approach to the Armstrong-Frederick model for
the simulation of ductile damage under cyclic load

Serhat Aygün1, Tillmann Wiegold2 and Sandra Klinge1,

1Chair of Structural Mechanics and Analysis, TU Berlin,
Straße des 17. Juni 135, 10623 Berlin, Germany
2Institute of Mechanics, TU Dortmund University,
Leonhard-Euler-Straße 5, 44227 Dortmund, Germany

Correspondence<EMAIL_ADDRESS>
###### Abstract

The present contribution proposes a thermodynamically consistent model for the
simulation of the ductile damage. The model couples the phase field method of
fracture to the Armstrong-Frederick plasticity model with kinematic hardening.
The latter is particularly suitable for simulating the material behavior under
a cyclic load. The model relies on the minimum principle of the dissipation
potential. However, the application of this approach is challenging since
potentials of coupled methods are defined in different spaces: The dissipation
potential of the phase field model is expressed in terms of rates of internal
variables, whereas the Armstrong-Frederick model proposes a formulation
depending on thermodynamic forces. For this reason, a unique formulation
requires the Legendre transformation of one of the potentials. The present
work performs the transformation of the Armstrong-Frederick potential, such
that final formulation is only expressed in the space of rates of internal
variables. With the assumption for the free energy and the joint dissipation
potential at hand, the derivation of evolution equations is straightforward.
The application of the model is illustrated by selected numerical examples
studying the material response for different load constellations and sample
geometries. The paper provides a comparison with the experimental results as
well.

#### Keywords:

crack mechanics; cyclic loading; ductility; elastic-plastic material; fracture

##  1 Introduction

Fatigue failure significantly influences the safety measures during the
service life of engineering structures. It particularly holds for structures
subjected to alternate loading with high amplitudes, which is an often
situation when speaking of transportation vehicles. Dependent on the number of
cycles and load amplitude, the interplay of the brittle and ductile failure
takes place. However, the present contribution focuses on the low cycle
fatigue where the role of plasticity is essential. The high complexity of the
ductile damage under the cyclic load strongly motivates the virtual testing of
materials on the basis of the computer simulations. The present work envisages
a coupling of the phase field method of fracture to the Armstrong-Frederick
model of plasticity with the kinematic hardening to this end.

The base for the phase field model of fracture has been set in the work by
Francfort and Marigo (1998) proposing the formulation of brittle fracture
based solely on Griffith’s idea of competition between elastic and fracture
energy. This work initiated a large group of further contributions related to
the numerical implementation of the concept proposed (Bourdin and Chambolle,
2000; Negri and Paolini, 2001; Fraternali, 2007; Schmidt et al., 2009).
Amongst others, Ambrosio-Tortorelli regularizations have become ubiquitous
(Ambrosio and Tortorelli, 1990, 1992; Bourdin et al., 2000). These approaches
are nowadays known as phase-field models of fracture and share several common
features with the approaches resulting from Ginzburg-Landau models for phase
transition (Karma et al., 2001). These regularizations have been applied to a
wide variety of fracture problems including fracture of thermal and drying
cracks (Maurini et al., 2013; Bourdin et al., 2014), ferro-magnetic and piezo-
electric materials (Abdollahi and Arias, 2012; Wilson et al., 2013) and
hydraulic fracturing (Bourdin et al., 2012; Wheeler et al., 2014; Wilson and
Landis, 2016). These models also been enhanced to account for cohesive effects
(Crismale and Lazzaroni, 2016; Freddi and Iurlano, 2017), ductile behavior
(Alessi et al., 2014; Miehe et al., 2015; Ambati et al., 2015a), large
deformations (Ambati et al., 2015b; Miehe et al., 2016; Borden et al., 2016),
quasi-brittle damage (Narayan and Anand, 2019) and anisotropy (Li et al.,
2014).

Some of previous works have also considered coupling of plasticity to the
phase field method (Alessi et al., 2014; Ambati et al., 2015a, b; Borden et
al., 2016; Miehe et al., 2016, 2015). However, the current contribution will
focus more precisely on the plasticity under the cyclic load where the
Armstrong-Frederick model is a widely accepted concept (Bari and Hassan, 2001;
Khan and Jackson, 1999; Kobayashi and Ohno, 2002; Puzrin and Houlsby, 2001).
Some of the works in this field deal with so-called generalized
representations, where either the classical Armstrong-Frederick model is
extended to include several Armstrong-Frederick terms (Bari and Hassan, 2001),
or where the application of multiple yield surfaces is envisaged (Puzrin and
Houlsby, 2001). Another frequently used variant of the Armstrong-Frederick
model has been presented by Ohno and Wang (1993) and by Kobayashi and Ohno
(2002). The effects of the kinematic hardening are nowadays also investigated
within the framework of large deformations (Lührs et al., 1997; Svendsen et
al., 1998; Lion, 2000; Mollica et al., 2001).

The Armstrong-Frederick model coupled to the phase field method is a promising
approach inheriting the advantages of both incorporated techniques. However,
the numerical implementation of this coupled approach brings with it several
challenges such as the definition of a unique framework for both setups, the
derivation of coupled evolution equations, the distinction between the tension
and compression mode and certainly the development of the computationally
efficient algorithm. Some of these issues are discussed in the present
contribution which is structured as follows. An overview on the phase field
method of fracture is presented in Sect. 2, whereas Sect. 3 introduces the
Armstrong-Frederick model of plasticity. Furthermore, the Armstrong-Frederick
evolution equations are derived by using two versions of the minimum principle
of dissipation potential (Sects. 4 and 5). The concept for the coupling of two
methods is presented in Sect. 6 and complemented by explaining the scheme for
the calculation of the plasticity multiplier on the basis of the consistency
condition (Sect. 7). Details on the numerical implementation including the
global part and the material point part are presented in Sect. 8. Eventually,
the paper also includes numerical examples simulating the tests with the
uniformly increasing and the cyclic load. The contribution finishes with
conclusions and an outlook.

##  2 Diffusive crack topology

The derivation of the coupled method requires a short recapitulation of
separate approaches, which is done in two subsequent sections. Sect. 2
provides an overview of the phase field method, whereas Sect. 3 focuses on the
Armstrong-Frederick model.

The present work uses the phase field framework proposed by Miehe and
coworkers (Miehe et al., 2010a, b) as a basis. This approach supposes the
diffusive (regularized) crack topology instead of the sharp crack topology
leading to serious difficulties in the numerical implementation due to the
lack of the continuity and differentiability. It relies on the evaluation of
the damage variable $d$, which continuously changes in the range $[0,1]$,
where zero-value corresponds to the intact material and the unity-value to the
fully broken state. The approach starts by postulating an expression for the
crack surface density function depending on the damage variable $d$

$\gamma(d,\nabla d)=\frac{1}{2l}d^{2}+\frac{l}{2}\lvert\nabla d\rvert^{2},$
(1)

which is furthermore used to define the total crack surface $A$ for the entire
body $\mathcal{B}$

$A(d)=\int_{\mathcal{B}}\gamma(d,\nabla d)\,dV.$ (2)

An assumption that the crack propagation is a fully dissipative process yields
to the conclusion that the constitutive dissipation potential has to be
proportional to the rate of crack surface density

$\Phi=g_{\mathrm{c}}\dot{\gamma}(\dot{d},\nabla\dot{d},d,\nabla d),$ (3)

where parameter $g_{\mathrm{c}}$ is related to the critical Griffith-type
fracture energy and can be seen as a constitutive threshold value. However,
this potential still does not guarantee that the damage is an increasing
function. For this reason, the potential is extended by introducing the
penalty term

$P(\dot{d})=\frac{k_{\mathrm{p}}}{2}\langle\dot{d}\rangle^{2}_{-}.$ (4)

Here, the negative Macaulay brackets
$\langle\bullet\rangle_{-}=(\bullet-\lvert\bullet\rvert)/2$ assure that the
penalty term only activates for a negative damage evolution and the constant
$k_{\mathrm{p}}$ has to be chosen as high as possible in order to stipulate
the condition $\dot{d}\geq 0$. Along with the penalty term (4), the damage
potential turns into

$\Phi^{\mathrm{d}}=\Phi+P(\dot{d})=g_{\mathrm{c}}\dot{\gamma}(\dot{d},\nabla\dot{d},d,\nabla
d)+\frac{k_{\mathrm{p}}}{2}\langle\dot{d}\rangle^{2}_{-}.$ (5)

In a further step, the model focuses on the definition of the free energy and
its reduction on the basis of function $\omega=((1-d)^{2}+k_{\mathrm{d}})$.
Here, the positive constant $k_{\mathrm{d}}$ prevents the energy to become
identical to zero at a fully broken state. Moreover, the reduction function
only affects the tension part of the energy which corresponds to the concept
of the anisotropic degradation of energy

${\Psi}^{\mathrm{d}}(\boldsymbol{\epsilon},d)=\,\omega(d)\,{\Psi}_{+}(\boldsymbol{\epsilon})+{\Psi}_{-}(\boldsymbol{\epsilon}).$
(6)

The additive split of the energy into a tension part
${\Psi}_{+}(\boldsymbol{\epsilon})$ and a compression part
${\Psi}_{-}(\boldsymbol{\epsilon})$, is of the special interest for simulating
the cyclic behavior where the crack closure in the compression mode has to be
considered. The definitions of the positive and negative energy parts of the
free energy are based on the spectral decomposition of the strain tensor

$\boldsymbol{\epsilon}_{\pm}=\sum_{i}\langle\epsilon_{i}\rangle_{\pm}\boldsymbol{n}_{i}\otimes\boldsymbol{n}_{i},$
(7)

where $\epsilon_{i}$ are the principal strains, $\boldsymbol{n}_{i}$ are the
principal strain directions and $i$ is the summation index. Equation (7) uses
positive and negative Macaulay brackets. Positive ones are defined as
$\langle\bullet\rangle_{+}=(\bullet+\lvert\bullet\rvert)/2$, whereas the
negative ones have already been used in Eq. (4). The stresses corresponding to
(7) are then defined as

$\boldsymbol{\sigma}=\frac{\partial\Psi^{\mathrm{d}}}{\partial{\boldsymbol{\epsilon}}}=\omega(d)\,\boldsymbol{\sigma}_{+}+\boldsymbol{\sigma}_{-}$
(8)

such that it holds
$\boldsymbol{\sigma}_{\pm}=\frac{\partial\Psi_{\pm}}{\partial\boldsymbol{\epsilon}}$.
For a case of an isotropic material characterized by Lamé constants $\lambda$
and $\mu$, the energy spit turns into

$\Psi_{\pm}(\boldsymbol{\epsilon})=\frac{1}{2}\lambda(\mathrm{tr}\langle\boldsymbol{\epsilon}\rangle_{\pm})^{2}+\mu\mathrm{tr}(\langle\boldsymbol{\epsilon}\rangle^{2}_{\pm})$
(9)

with the corresponding equilibrium equation and stresses

$\displaystyle\mathrm{Div}[\,\omega(d)\,\boldsymbol{\sigma}_{+}+\boldsymbol{\sigma}_{-}]=\boldsymbol{0},$
(10)
$\displaystyle\boldsymbol{\sigma}_{\pm}=\sum_{i}{[\lambda\mathrm{tr}\langle\boldsymbol{\epsilon}\rangle_{\pm}+2\mu\langle\epsilon_{i}\rangle_{\pm}]\,\boldsymbol{n}_{i}\otimes\boldsymbol{n}_{i}}.$
(11)

Since the free energy and the dissipation potential are known, the minimum
principle of dissipation potential is used to derive the evolution equation
for the internal parameter, namely

$\displaystyle\underset{\dot{d}}{\text{min}}(\mathcal{L}^{\mathrm{MDP}}=\dot{\Psi}^{\mathrm{d}}(\dot{\boldsymbol{\epsilon}},\dot{d})+\Phi^{\mathrm{d}}(d,\dot{d}))\qquad\Rightarrow\qquad\frac{g_{c}}{l}[d-l^{2}\Delta
d]+[\omega^{\prime}\,\Psi_{+}+k_{\mathrm{p}}\langle\dot{d}\rangle_{-}]=0.$
(12)

The equilibrium equation (10) together with the evolution equation (12)b
defines the strong problem of brittle fracture.

##  3 The Armstrong-Frederick kinematic hardening model

The Armstrong-Frederick kinematic hardening model (Armstrong and Frederick,
1966) is widely accepted to simulate the characteristic phenomena of the
hardening behavior of metals, namely the Bauschinger and the ratcheting
effect. Within the framework of the small strain plasticity, the formulation
of the Armstrong-Frederick model starts with the typical assumption for the
additive decomposition of strain tensor $\boldsymbol{\epsilon}$ into an
elastic part $\boldsymbol{\epsilon}^{\mathrm{e}}$ and a plastic part
$\boldsymbol{\epsilon}^{\mathrm{p}}$, namely
$\boldsymbol{\epsilon}=\boldsymbol{\epsilon}^{\mathrm{e}}+\boldsymbol{\epsilon}^{\mathrm{p}}.$
However, the description of the hardening effects requires the introduction of
an additional strain-like internal variable $\boldsymbol{\xi}$. According to
Lion (Lion, 2000), $\xi$ relates to local viscoelastic deformations induced by
dislocations. In rheological models, it is simulated as a damping element
connected in series to an elastic spring. These are additionally coupled to a
friction element in parallel (Dettmer and Reese, 2004). This internal variable
is responsible for the change of material stiffness depending on the loading
history and is incorporated in the Helmholtz free energy, which now includes
two parts: the elastic energy
$\Psi^{\mathrm{el}}(\boldsymbol{\epsilon}^{\mathrm{e}})$ and free energy due
to the hardening $\Psi^{\mathrm{hard}}(\boldsymbol{\xi})$, both in the
quadratic form

$\Psi^{\mathrm{AF}}=\Psi^{\mathrm{el}}(\boldsymbol{\epsilon}^{\mathrm{e}})+\Psi^{\mathrm{hard}}(\boldsymbol{\xi})=\frac{1}{2}(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}}):\mathbb{C}:(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}})+\frac{1}{2}c\
\boldsymbol{\xi}:\boldsymbol{\xi}.$ (13)

Here, $\mathbb{C}$ is the fourth order elasticity tensor and $c$ represents
the kinematic hardening modulus. The Helmholtz free energy serves furthermore
as a basis for the formulation of dissipation which, in the case of an
isothermal process, only depends on the elastic power
$(\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}})$ and free energy rate
$(\dot{\Psi}^{\mathrm{AF}})$

$\mathcal{D}=\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}-\dot{\Psi}^{\mathrm{AF}}\geq
0.$ (14)

In the concrete case of the Armstrong-Frederick model, the dissipation turns
into:

$(\boldsymbol{\sigma}-\frac{\partial\Psi^{\mathrm{AF}}}{\partial\boldsymbol{\epsilon}}):\dot{\boldsymbol{\epsilon}}{-\frac{\partial\Psi^{\mathrm{AF}}}{\partial\boldsymbol{\epsilon}^{\mathrm{p}}}}:\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}{-\frac{\partial\Psi^{\mathrm{AF}}}{\partial\boldsymbol{\xi}}}:\dot{\boldsymbol{\xi}}\geq
0,$ (15)

which yields two groups of consequences. The first consequence is the
constitutive law
$\boldsymbol{\sigma}=\frac{\partial\Psi^{\mathrm{AF}}}{\partial\boldsymbol{\epsilon}}$,
whereas the second consequence defines the driving forces of internal
variables:

$\displaystyle\boldsymbol{q}_{\boldsymbol{\epsilon}^{\mathrm{p}}}=-\frac{\partial\Psi^{\mathrm{AF}}}{\partial\boldsymbol{\epsilon}^{\mathrm{p}}}=\mathbb{C}:[\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}}]=\boldsymbol{\sigma},$
(16)
$\displaystyle\boldsymbol{q}_{\boldsymbol{\xi}}=-\frac{\partial\Psi^{\mathrm{AF}}}{\partial\boldsymbol{\xi}}=-c\
\boldsymbol{\xi}=\boldsymbol{\chi}.$ (17)

The later definition corresponds to the back stress and gives insight into the
physical meaning of material parameter $c$. It is a proportionality constant
relating the back stress $\chi$ to the strain-like quantity $\xi$. In analogy
to the constitutive law of elastic materials, $c$ represents the kinematic
hardening modulus. A shorter notation for internal variables
$\boldsymbol{\nu}=\\{{\boldsymbol{\epsilon}^{\mathrm{p}}},{\boldsymbol{\xi}}\\}$
and for corresponding driving forces
$\boldsymbol{q}=\\{\boldsymbol{\sigma},\boldsymbol{\chi}\\}$, together with
the constitutive law for stresses, yields the result for the so-called reduced
dissipation

$\mathcal{D}^{\mathrm{red}}=\boldsymbol{q}:\dot{\boldsymbol{\nu}}=\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}+\boldsymbol{\chi}:\dot{\boldsymbol{\xi}}\geq
0.$ (18)

The Armstrong-Frederick model is eventually completed by introducing the yield
locus formula

$\Omega^{\mathrm{AF}}=\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert-\bar{\sigma}_{\mathrm{Y}}=0,\qquad\bar{\sigma}_{\mathrm{Y}}=\sqrt{2/3}\sigma_{\mathrm{Y}},$
(19)

where symbol $\rvert\rvert\bullet\lvert\lvert=\sqrt{\bullet:\bullet}$ denotes
the Frobenius norm, $\boldsymbol{\bar{\sigma}}$ and $\boldsymbol{\bar{\chi}}$
are deviatoric parts of stresses $\boldsymbol{\sigma}$ and back stresses
$\boldsymbol{\chi}$ and $\sigma_{\mathrm{Y}}$ is the yield limit. Condition
(19) defines the admissible set of driving forces yielding inelastic
deformations.

##  4 Derivation of the evolution equations based on the minimum principle for
the dissipation potential in terms of driving forces

The completion of the Armstrong-Frederick model previously described also
requires the derivation of evolution equations which can be accomplished by
using different concepts. One possibility is to follow the minimum principle
of dissipation potential (MDP) as proposed in works by Dettmer and Reese
(2004) and by Aygün and Klinge (2020)

$\displaystyle\underset{\boldsymbol{q}}{\text{min}}\left(\mathcal{L}^{\mathrm{MDP^{*}}}=-\boldsymbol{q}:\dot{\boldsymbol{\nu}}+\Phi^{\mathrm{AF^{*}}}\left(\boldsymbol{\nu},\boldsymbol{q}\right)\right),$
(20)
$\displaystyle\Phi^{\mathrm{AF^{*}}}(\boldsymbol{\nu},\boldsymbol{q})=a\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert+\frac{1}{2b}\lvert\lvert\boldsymbol{\bar{\chi}}\rvert\rvert^{2}.$
(21)

The dissipation potential (21) is inspired by the expression for the yield
locus (19) but additionally includes a term depending on the norm of driving
force $\boldsymbol{\chi}$. Symbols $a$ and $b$ denote material parameters.
However, the subsequent derivations will show that parameter $a$ does not
influence the evolution equations, whereas parameter $b$ plays an important
role and represents pseudo-viscoelasticty. Superscript $*$ indicates that a
formulation in terms of driving forces is chosen. The Lagrangian corresponding
to the minimization problem (20) and (21) has the form

$\displaystyle\mathcal{L}^{\mathrm{MDP^{*}}}=-\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}-\boldsymbol{\chi}:\dot{\boldsymbol{\xi}}+\Phi^{\mathrm{AF^{*}}}(\boldsymbol{\nu},\boldsymbol{q}),$
(22)

which yields the following stationary conditions

$\displaystyle\frac{\partial\mathcal{L}^{\mathrm{MDP^{*}}}}{\partial\boldsymbol{q}}=\left(\begin{array}[]{c}\frac{\partial\mathcal{L}^{\mathrm{MDP^{*}}}}{\partial\boldsymbol{\sigma}}\\\\[5.69046pt]
\frac{\partial\mathcal{L}^{\mathrm{MDP^{*}}}}{\partial\boldsymbol{\chi}}\\\
\end{array}\right)=\left(\begin{array}[]{c}-\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+\frac{\partial\Phi^{\mathrm{AF^{*}}}}{\partial\boldsymbol{\sigma}}\\\
-\dot{\boldsymbol{\xi}}+\frac{\partial\Phi^{\mathrm{AF^{*}}}}{\partial\boldsymbol{\chi}}\\\
\end{array}\right)=\boldsymbol{0}.$ (27)

The first equation in (27) determines the evolution of plastic strains and
shows that the rate $\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}$ is equal to
the derivative of dissipation potential with respect to corresponding driving
force:
$\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}=\frac{\partial\Phi^{\mathrm{AF^{*}}}}{\partial\boldsymbol{\sigma}}=\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert}.$
However, the received derivative is not uniquely defined and only contains the
information on the direction of the flow of the plastic strains, not on the
magnitude. For that reason, the right-hand side expression is scaled by the
plastic multiplier $\lambda$ which yields the standard solution

$\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}=\frac{\partial\Phi^{\mathrm{AF^{*}}}}{\partial\boldsymbol{\sigma}}=\lambda\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert}.$
(28)

The second equation in (27) determines the evolution of the internal variable
$\boldsymbol{\xi}$. Here, the same argumentation as in Eq. (28) is used for
the first term, whereas the second term is uniquely determined and does not
need to be scaled with the multiplier $\lambda$

$\dot{\boldsymbol{\xi}}=\frac{\partial\Phi^{\mathrm{AF^{*}}}}{\partial\boldsymbol{\chi}}=-\lambda\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert}+\frac{1}{b}\boldsymbol{\bar{\chi}}=-\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+\frac{1}{b}\boldsymbol{\bar{\chi}}\qquad\Rightarrow\qquad\boldsymbol{\bar{\chi}}=b(\dot{\boldsymbol{\xi}}+\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}).$
(29)

Equations (28) and (29)a show the deviatoric character of internal variables
${\boldsymbol{{\epsilon}}}^{\mathrm{p}}$ and ${\boldsymbol{\xi}}$. Moreover,
Eq. (29)b represents a constitutive law typical of a viscous material with the
viscosity $b$. Finally, Eq. (17) along with Eq. (29)a provides the evolution
equation for back stress $\chi$ which can also be identified as a deviatoric
quantity:

$\dot{\boldsymbol{\chi}}=-c\dot{\boldsymbol{\xi}}=c\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}-\frac{c}{b}\boldsymbol{\bar{\chi}}=\dot{\bar{\boldsymbol{\chi}}}.$
(30)

The MDP-approach differs from the common approach treating the yield locus
equation $\Omega$ as a subsidiary condition within the Lagrange formalism of
the constrained optimization. The main difference manifests itself in the
evolution equations taking the form
$\nobreak{\dot{\boldsymbol{\nu}}=\lambda\frac{\partial\Omega}{\partial\boldsymbol{q}}}$
with $\lambda$ as the Lagrange multiplier. An important advantage of the MDP-
approach compared to the constrained optimization method is that it allows a
standardized variable transformation, which is often a useful tool in coupling
strategies.

##  5 Derivation of the evolution equations based on the minimum principle for
the dissipation potential in terms of rates of internal variables

Both models summarized in Sects. 2 and 4 apply the minimum principle of
dissipation potential. However, the phase field model uses a formulation in
terms of the rate of internal variables, whereas the Armstrong-Frederick model
minimizes the dissipation potential in terms of driving forces. Naturally, the
coupling procedure requires a unique formulation in a single space. The
present contribution deals with the approach in terms of velocities, such that
the transformation of the potential
$\Phi^{\mathrm{AF^{*}}}(\boldsymbol{\nu},\boldsymbol{q})$ into the space of
rates is necessary. This type of exchange is conducted on the basis of
Legendre transformation (LT) defining the new potential as follows

$\displaystyle\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}})=\underset{{\boldsymbol{q}}}{\text{max}}\\{\mathcal{L}^{\mathrm{LT}}=\boldsymbol{q}:\dot{\boldsymbol{\nu}}-\Phi^{\mathrm{AF^{*}}}(\boldsymbol{\nu},{\boldsymbol{q}})\\},$
(31)
$\displaystyle\mathcal{L}^{\mathrm{LT}}=\boldsymbol{\sigma}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+\boldsymbol{\chi}:\dot{\boldsymbol{\xi}}-a||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||-\frac{1}{2b}||\boldsymbol{\bar{\chi}}||^{2}.$
(32)

The maximization procedure relies on two stationary conditions, the first of
which gives a relationship for rates of plastic deformations

$\frac{\partial\mathcal{L}^{\mathrm{LT}}}{\partial\boldsymbol{\sigma}}=\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}-\lambda\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert}=0.$
(33)

Bearing in mind that the multiplier $\lambda$ is a scalar, Eq. (33) shows that
rate $\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}$ and difference
$\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}$ are coaxial. In other
words, it holds that
$\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}=\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert}.$
Moreover, the plastic flow only occurs if the stress state fulfills the yield
locus equation (19) such that the norm
$||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||$ can be replaced by the
yield limit $\bar{\sigma}_{\mathrm{Y}}$. Equation (33) is then rewritten as
follows

$\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}\qquad\Rightarrow\qquad\boldsymbol{\bar{\sigma}}=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}+\boldsymbol{\bar{\chi}}.$
(34)

On the other hand, the second stationary condition

$\frac{\partial\mathcal{L}^{\mathrm{LT}}}{\partial\boldsymbol{\chi}}=\dot{\boldsymbol{\xi}}+\lambda\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||}-\frac{1}{b}\boldsymbol{\bar{\chi}}=0\qquad\Rightarrow\qquad\boldsymbol{\bar{\chi}}=b\dot{\boldsymbol{\xi}}+b\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}$
(35)

together with (34)b leads to the final expression for deviatoric stresses

$\boldsymbol{\bar{\sigma}}=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}+b\dot{\boldsymbol{\xi}}+b\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}.$
(36)

Driving forces (35)b and (36) are now inserted into (32), which reads the
desired dissipation potential in terms of the internal variables and their
rates

$\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}})=\bar{\sigma}_{\mathrm{Y}}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||+\frac{b}{2}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||^{2}+\frac{b}{2}||\dot{\boldsymbol{\xi}}||^{2}+b\,\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}\\!:\\!\dot{\boldsymbol{\xi}}.$
(37)

More details on derivation are presented in Appendix A. Result (37) now
enables the formulation of the new minimization problem

$\underset{\dot{\boldsymbol{\nu}}}{\text{min}}\left(\mathcal{L}^{\mathrm{MDP}}=\dot{\Psi}^{\mathrm{AF}}+\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}})\right),\qquad\mathcal{L}^{\mathrm{MDP}}=\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}-\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}-\boldsymbol{\chi}:\dot{\boldsymbol{\xi}}+\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}}).$
(38)

Bearing in mind the deviatoric character of quantities
$\dot{\boldsymbol{\epsilon}}$ and $\dot{\boldsymbol{\xi}}$, the problem can
also be written as:

$\underset{\dot{\boldsymbol{\nu}}}{\text{min}}\left(\mathcal{L}^{\mathrm{MDP}}=\dot{\Psi}^{\mathrm{AF}}+\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}})\right),\qquad\mathcal{L}^{\mathrm{MDP}}=\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}-\boldsymbol{\bar{\sigma}}:\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}-\boldsymbol{\bar{\chi}}:\dot{\boldsymbol{\xi}}+\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}}).$
(39)

Equation (39) is the counterpart of the MDP-problem (20) in terms of
velocities and provides expressions for driving forces according to
relationship
$\boldsymbol{q}=\frac{\partial\Phi^{\mathrm{AF}}}{\partial{\dot{\boldsymbol{\nu}}}}$:

$\displaystyle\boldsymbol{q}_{\boldsymbol{\epsilon}^{\mathrm{p}}}=\boldsymbol{\bar{\sigma}}=\frac{\partial\Phi^{\mathrm{AF}}}{\partial{\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}}}=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}+b{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}+b{\dot{\boldsymbol{{\xi}}}},$
(40)
$\displaystyle\boldsymbol{q}_{\boldsymbol{\xi}}=\boldsymbol{\bar{\chi}}=\frac{\partial\Phi^{\mathrm{AF}}}{\partial{\dot{\boldsymbol{\xi}}}}=b{\dot{\boldsymbol{{\xi}}}}+b{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}.$
(41)

A transformation of the system (40) and (41) yields the evolution equations
identical to the ones from Sect. 4 (Eqs. (28) and (29)):

${\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}=\lambda\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||},\qquad\dot{\boldsymbol{\xi}}=-{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}+\frac{1}{b}\boldsymbol{\bar{\chi}}.$
(42)

Note that the insertion of (41) into (40) yields the intermediate result

$\boldsymbol{\bar{\sigma}}=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}+\boldsymbol{\chi}\qquad\Rightarrow\qquad\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||},$
(43)

which is only valid for $||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||\neq
0$. By taking the norm of (43)b, it follows that
${||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||}=\bar{\sigma}_{\mathrm{Y}}$.
This consequence indicates that dissipation potential (37) intrinsically
includes the yield locus condition. Accordingly, it can be summed up that the
dissipation potential (37) describes the same problem as the dissipation
potential (21) along with the yield locus function (19).

Finally, the strong form corresponding to the Armstrong-Frederick problem can
be recapitulated as follows:

$\displaystyle\mathrm{Div}\,\boldsymbol{\sigma}=\boldsymbol{0},$ (44)
$\displaystyle\boldsymbol{\sigma}=\mathbb{C}:(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}}),$
(45)
$\displaystyle\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}=\lambda\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||},$
(46)
$\displaystyle\dot{\boldsymbol{\chi}}=c\left(\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}-\
\frac{1}{b}\boldsymbol{\bar{\chi}}\right),$ (47)
$\displaystyle\Omega^{\mathrm{AF}}=||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||-\bar{\sigma}_{\mathrm{Y}},$
(48) $\displaystyle\lambda\geq 0,\ \Omega^{\mathrm{AF}}\leq 0,\
\lambda\Omega^{\mathrm{AF}}=0.$ (49)

The previous system includes the equilibrium equation, the constitutive law,
two evolution equations and the Karush-Kuhn-Tucker conditions defining the
plastic domain. Equation (47) is derived by using the time derivative of
definition (17) and rate (42)b. The body forces in the equilibrium equation
are neglected. Evolution equation (46) and Karush-Kuhn-tucker conditions
depend on plastic multiplier $\lambda$ which is typically determined from the
consistency condition.

##  6 Coupling of the Armstrong-Frederick model to the phase-field approach

The coupling of the two methods starts by writing the Armstrong-Frederick
energy in a form splitting the tension and compression part of elastic energy
in order to introduce the damage influence as was done in Sect. 2

$\Psi^{\mathrm{AF}}(\boldsymbol{\epsilon}^{\mathrm{e}},\boldsymbol{\xi})=\Psi^{\mathrm{el}}(\boldsymbol{\epsilon}^{\mathrm{e}})+\Psi^{\mathrm{hard}}(\boldsymbol{\xi})=\Psi^{\mathrm{el}}_{+}(\boldsymbol{\epsilon}^{\mathrm{e}})+\Psi^{\mathrm{el}}_{-}(\boldsymbol{\epsilon}^{\mathrm{e}})+\Psi^{\mathrm{hard}}(\boldsymbol{\xi}).$
(50)

However, damage also influences the hardening energy in the case where the
tension mode is active. For this reason, the present model introduces a
function distinguishing the pure compression mode from the pure tension and
from the mixed modes

$\bar{\omega}=\begin{cases}\omega(d)&\text{if}\quad\text{max}\\{(\epsilon_{i})_{i=1,2,3},0\\}>0,\\\
1&\text{if}\quad\text{max}\\{(\epsilon_{i})_{i=1,2,3},0\\}\leq 0.\end{cases}$
(51)

By using this new notation, the coupled free energy is written as

$\Psi^{\mathrm{c}}(\boldsymbol{\epsilon}^{\mathrm{e}},d,\boldsymbol{\xi})=\omega\Psi^{\mathrm{el}}_{+}(\boldsymbol{\epsilon}^{\mathrm{e}})+\Psi^{\mathrm{el}}_{-}(\boldsymbol{\epsilon}^{\mathrm{e}})+\bar{\omega}\Psi^{\mathrm{hard}}(\boldsymbol{\xi})$
(52)

and corresponding constitutive laws take the form

$\displaystyle\boldsymbol{\sigma}=\frac{\partial\Psi^{\mathrm{c}}}{\partial\boldsymbol{\epsilon}}=\omega\frac{\partial\Psi^{\mathrm{el}}_{+}}{\partial\boldsymbol{\epsilon}}+\frac{\partial\Psi^{\mathrm{el}}_{-}}{\partial\boldsymbol{\epsilon}}=\omega\boldsymbol{\sigma}_{+}+\boldsymbol{\sigma}_{-},$
(53)
$\displaystyle\boldsymbol{\sigma}_{\pm}=\frac{\partial\Psi^{\mathrm{el}}_{\pm}}{\boldsymbol{\epsilon}}=\mathbb{C}:\boldsymbol{\epsilon}_{\pm}^{\mathrm{e}}=\mathbb{C}:\langle\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}}\rangle_{\pm},$
(54)
$\displaystyle\boldsymbol{\chi}=\frac{\partial\Psi^{\mathrm{c}}}{\partial{\boldsymbol{\xi}}}=-\bar{\omega}c\boldsymbol{\xi},$
(55)

where the spectral decomposition of elastic strains has the standard form

$\boldsymbol{\epsilon}_{\pm}^{\mathrm{e}}=(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}})_{\pm}=\sum_{i}\langle\epsilon^{\mathrm{e}}_{i}\rangle_{\pm}\boldsymbol{n}_{i}\otimes\boldsymbol{n}_{i}=\sum_{i}\langle(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}})_{i}\rangle_{\pm}\boldsymbol{n}_{i}\otimes\boldsymbol{n}_{i}.$
(56)

Bearing in mind stress definition (53), the equilibrium equation turns into

$\mathrm{Div}[\omega(d)\boldsymbol{\sigma}_{+}+\boldsymbol{\sigma}_{-}]=\boldsymbol{0}.$
(57)

The damage influence on the dissipation functional related to the plastic
deformations follows the same argumentation as in the case of the hardening
energy, such that the coupled dissipation potential $\Phi^{\mathrm{c}}$
consists of two terms where the second term is weighted by the function
$\bar{\omega}$

$\Phi^{\mathrm{c}}=\Phi^{\mathrm{d}}+\bar{\omega}\,\Phi^{\mathrm{AF}}.$ (58)

With the definitions (52) and (58), the minimization of the corresponding
Lagrange function

$\mathcal{L}^{\mathrm{MDP}}=\dot{\Psi}^{\mathrm{c}}+\Phi^{\mathrm{c}}=\dot{\Psi}^{\mathrm{c}}+\Phi^{\mathrm{d}}+\bar{\omega}\Phi^{\mathrm{AF}}$
(59)

yields the following system of equations which, together with equilibrium
equation (57), defines the strong form of the coupled problem

$\displaystyle\boldsymbol{\sigma}=\omega(d)\boldsymbol{\sigma}_{+}+\boldsymbol{\sigma}_{-},$
(60)
$\displaystyle\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}=\lambda\bar{\boldsymbol{N}},\qquad\bar{\boldsymbol{N}}=\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||},$
(61)
$\displaystyle\dot{\boldsymbol{\xi}}=\frac{1}{b\bar{\omega}}\boldsymbol{\bar{\chi}}-\dot{\boldsymbol{\epsilon}}^{\mathrm{p}},$
(62)
$\displaystyle\dot{\boldsymbol{\chi}}=-\dot{\bar{\omega}}c\boldsymbol{\xi}+c\bar{\omega}\left(\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}-\
\frac{1}{b\bar{\omega}}\ \boldsymbol{\bar{\chi}}\right),$ (63)
$\displaystyle\Omega^{\mathrm{c}}=||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||-\bar{\omega}\bar{\sigma}_{\mathrm{Y}},$
(64) $\displaystyle\lambda\geq 0,\ \Omega^{\mathrm{c}}\leq 0,\
\lambda\Omega^{\mathrm{c}}=0,$ (65)
$\displaystyle\frac{g_{c}}{l}[d-l^{2}\Delta
d]+[\omega^{\prime}\,\Psi^{\mathrm{el}}_{+}+\bar{\omega}^{\prime}\,\Psi^{\mathrm{hard}}+k_{\mathrm{p}}\langle\dot{d}\rangle_{-}]=0.$
(66)

Here, a shorter notation
$\bar{\boldsymbol{N}}=\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||}$
is introduced to denote the so-called direction tensor. Details on the
derivation of evolution equations and of the yield criterion are provided in
Appendix B.

##  7 Determination of the plastic multiplier

Evolution equation (61) and the Karush-Kuhn-Tucker conditions (65) depend on
plastic multiplier $\lambda$ which is commonly determined from the consistency
condition $\lambda\dot{\Omega}^{\mathrm{c}}=0$. An appropriate form of this
condition is obtained by taking the time derivative of the yield locus formula

$\lvert\lvert\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}\rvert\rvert^{2}=\bar{\omega}^{2}\bar{\sigma}_{\mathrm{Y}}^{2}\qquad\Rightarrow\qquad(\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}):(\dot{\bar{\boldsymbol{\sigma}}}-\dot{\bar{\boldsymbol{\chi}}})=0.$
(67)

Since difference $(\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}})$ is
deviatoric, the condition above also can be written as follows:

$(\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}):(\dot{\boldsymbol{\sigma}}-\dot{\boldsymbol{\chi}})=0.$
(68)

In a further step, the stress rate

$\dot{\boldsymbol{\sigma}}=\omega^{\prime}\,\boldsymbol{\sigma}_{+}+\omega(d)\,\mathbb{C}\\!:\\!(\dot{\boldsymbol{\epsilon}}-\dot{\boldsymbol{\epsilon}}^{\mathrm{p}})_{+}+\mathbb{C}\\!:\\!(\dot{\boldsymbol{\epsilon}}-\dot{\boldsymbol{\epsilon}}^{\mathrm{p}})_{-}$
(69)

along with Eqs. (61) and (63) is introduced in (68) which yields

$(\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}):\left[\omega^{\prime}\,\boldsymbol{\sigma}_{+}+\omega\,\mathbb{C}\\!:\\![\dot{\boldsymbol{\epsilon}}-\lambda\bar{\boldsymbol{N}}]_{+}+\mathbb{C}\\!:\\![\dot{\boldsymbol{\epsilon}}-\lambda\bar{\boldsymbol{N}}]_{-}+\dot{\bar{\omega}}\,c\boldsymbol{\xi}-c\bar{\omega}\lambda\bar{\boldsymbol{N}}+\frac{c}{b}\boldsymbol{\bar{\chi}}\right]=0.$
(70)

However, the split of the elastic strain rate into a positive and a negative
part cannot be performed as a superposition of positive/negative parts of the
total and plastic strains. For this reason, the previous equation can only be
solved numerically by using an iterative procedure. For this purpose, the
initial guess $\lambda_{0}$ is obtained by neglecting the terms including the
stress rates

$(\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}):\left[\omega^{\prime}\,\boldsymbol{\sigma}_{+}+\dot{\bar{\omega}}\,c\boldsymbol{\xi}-c\bar{\omega}\lambda_{0}\bar{\boldsymbol{N}}+\frac{c}{b}\boldsymbol{\bar{\chi}}\right]=0.$
(71)

Thereafter, all subsequent steps follow the iteration rule for the calculation
of the updated value $\lambda_{i}$

$(\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}):\left[\omega^{\prime}\,\boldsymbol{\sigma}_{+}+\omega\,\mathbb{C}\\!:\\![\dot{\boldsymbol{\epsilon}}-\lambda_{i-1}\bar{\boldsymbol{N}}]_{+}+\mathbb{C}\\!:\\![\dot{\boldsymbol{\epsilon}}-\lambda_{i-1}\bar{\boldsymbol{N}}]_{-}+\dot{\bar{\omega}}\,c\boldsymbol{\xi}-c\bar{\omega}\lambda_{i}\bar{\boldsymbol{N}}+\frac{c}{b}\boldsymbol{\bar{\chi}}\right]=0.$
(72)

The iterative process stops if the prescribed accuracy is achieved.

##  8 Numerical implementation

###  8.1 General approach

The minimum principle of dissipation potential along with the time-incremental
variational principles represents an important tool for the numerical solution
of boundary value problems. Within this concept, the Lagrange function
$\cal{L}^{\mathrm{MDP}}$ corresponding to the body $\Omega$ is integrated over
a single time increment $[t_{n},t_{n+1}]$. In its original form, this integral
depends on external variables ($\boldsymbol{u}$), internal variables
$(\boldsymbol{\nu})$ and on the rates of internal variables, so-called
velocities $(\dot{\boldsymbol{\nu}})$. However, the velocities can be replaced
by the Euler-forward approximation
$\dot{\boldsymbol{\nu}}=(\boldsymbol{\nu}_{n+1}-\boldsymbol{\nu}_{n})/\Delta
t$, where $\Delta t=t_{n+1}-t_{n}$ is the time increment. As a consequence,
the result only depends on discrete values of external and internal variables
such that the time integration is performed as follows

$\int_{t_{n}}^{t_{n+1}}\\!\\!\int_{\Omega}{\mathcal{L}^{\mathrm{MDP}}}dVdt\approx\int_{\Omega}\left\\{\Psi(\boldsymbol{\epsilon}^{\mathrm{e}}_{n+1},\boldsymbol{\nu}_{n+1})-\Psi(\boldsymbol{\epsilon}^{\mathrm{e}}_{n},\boldsymbol{\nu}_{n})+\Delta
t\,\Phi(\boldsymbol{\nu}_{n+1},(\boldsymbol{\nu}_{n+1}-\boldsymbol{\nu}_{n+1})/\Delta
t)\right\\}dV.$ (73)

Here, displacement $\boldsymbol{u}$ is the only external variable since an
isothermal process is considered. In a further step, the minimization with
respect to the displacements and internal variables in the current time step
$n+1$ yields the sought solution. Within this procedure, term
$\Psi(\boldsymbol{\epsilon}^{\mathrm{e}}_{n},\boldsymbol{\nu}_{n})$ can be
neglected since it only depends on values in the previous time step $n$. By
adding the potential of external forces $l(t_{n+1},\boldsymbol{u}_{n+1})$ in
(73), the new, combined Lagrangian is constructed

$\mathcal{L}^{\mathrm{comb}}=\int_{\Omega}\left\\{\Psi(\boldsymbol{\epsilon}^{\mathrm{e}}_{n+1},\boldsymbol{\nu}_{n+1})+\Delta
t\,\Phi(\boldsymbol{\nu}_{n+1},(\boldsymbol{\nu}_{n+1}-\boldsymbol{\nu}_{n+1})/\Delta
t)\right\\}dV+l(t_{n+1},\boldsymbol{u}_{n+1}).$ (74)

Within the present work, the minimization of (74) is performed by an approach
consisting of two parts: the global level solution part and the material point
solution part. The former calculates deformation and damage, whereas the
latter evaluates the internal variables $\boldsymbol{\epsilon}_{\mathrm{p}}$
and $\boldsymbol{\xi}$ by using the predictor-corrector scheme. In the
continuation, each solution part is explained separately. The indexes related
to step $n+1$ are omitted in order to achieve a concise representation.

###  8.2 Global solution part

The definition of the global solution part (gl) starts with the reduced
formulation of the Lagrangian by only including the reduced elastic energy and
damage potential

$\mathcal{L}^{\mathrm{gl}}=\int_{\Omega}\mathcal{L}dV=\int_{\Omega}\omega(d)\Psi_{+}^{\mathrm{el}}(\boldsymbol{\epsilon}^{\mathrm{e}})+\Psi_{-}^{\mathrm{el}}(\boldsymbol{\epsilon}^{\mathrm{e}})dV+\int_{\Omega}\Delta
t\,\Phi^{\mathrm{d}}(d_{n},\nabla d_{n},d,\nabla d)dV+l(t,\boldsymbol{u}).$
(75)

By using the expression for the crack surface density (1) and definition (5)
the increment $\Delta\Phi^{\mathrm{d}}$ can be reconstructed as shown in Miehe
et al. (2010a, b)

$\Delta t\,\Phi^{\mathrm{d}}(d_{n},\nabla d_{n},d,\nabla
d)=\frac{g_{c}}{2l}(d^{2}-d_{n}^{2})+\frac{g_{c}l}{2}(|\nabla d|^{2}-|\nabla
d_{n}|^{2})+\frac{k_{\mathrm{p}}}{2\Delta t}\langle d-d_{n}\rangle^{2}_{-}.$
(76)

The FE-implementation now requires the derivatives of the Lagrangian
$\mathcal{L}$. The first derivatives are required in order to form the
residual

$\displaystyle\partial_{\boldsymbol{\epsilon}}\mathcal{L}=\omega\,\mathbb{C}\\!:\\!\boldsymbol{\epsilon}^{\mathrm{e}}_{+}+\mathbb{C}\\!:\\!\boldsymbol{\epsilon}^{\mathrm{e}}_{-},$
(77)
$\displaystyle\partial_{d}\mathcal{L}=\omega^{\prime}\,\,\Psi^{\mathrm{el}}_{+}+\frac{g_{c}}{l}d+\frac{k_{\mathrm{p}}}{\Delta
t}\langle d-d_{n}\rangle_{-},\qquad\partial_{\nabla d}\mathcal{L}=g_{c}l\nabla
d,$ (78)

whereas the second derivatives are needed for the definition of the stiffness
matrix

$\displaystyle\partial^{2}_{\boldsymbol{\epsilon}\boldsymbol{\epsilon}}\mathcal{L}=\omega\,\mathbb{C}\\!:\\!\boldsymbol{I}_{+}^{\boldsymbol{\epsilon}^{\mathrm{e}}}+\mathbb{C}\\!:\\!\boldsymbol{I}_{-}^{\boldsymbol{\epsilon}^{\mathrm{e}}},\qquad\hskip
35.56593pt\partial^{2}_{\boldsymbol{\epsilon}d}\mathcal{L}=\omega^{\prime}\,\mathbb{C}\\!:\\!\boldsymbol{\epsilon}_{+}^{\mathrm{e}},$
(79)
$\displaystyle\partial^{2}_{dd}\mathcal{L}=\omega^{\prime\prime}\,\Psi_{+}^{\mathrm{el}}+\frac{g_{c}}{l}+\frac{k_{\mathrm{p}}}{\Delta
t}I^{d},\qquad\qquad\hskip 5.69054pt\partial^{2}_{\nabla d\nabla
d}\mathcal{L}=g_{c}l.$ (80)

These relationships use the following indicator functions

$\displaystyle\boldsymbol{I}_{+}^{\boldsymbol{\epsilon}^{\mathrm{e}}}=\frac{\partial\boldsymbol{\epsilon}^{\mathrm{e}}_{+}}{\partial\boldsymbol{\epsilon}^{\mathrm{e}}}=\sum_{i}I_{i+}^{\mathrm{e}}\boldsymbol{m}_{i}\otimes\boldsymbol{m}_{i},\qquad\boldsymbol{m}_{i}=\boldsymbol{n}_{i}\otimes\boldsymbol{n}_{i},$
(81)
$\displaystyle\boldsymbol{I}_{-}^{\boldsymbol{\epsilon}^{\mathrm{e}}}=\frac{\partial\boldsymbol{\epsilon}^{\mathrm{e}}_{-}}{\partial\boldsymbol{\epsilon}^{\mathrm{e}}}=\sum_{i}(1-I_{i+}^{\mathrm{e}})\boldsymbol{m}_{i}\otimes\boldsymbol{m}_{i},$
(82) $\displaystyle
I_{i+}^{\mathrm{e}}=\begin{cases}1&\text{if}\quad\epsilon_{i}>0,\\\
0&\text{if}\quad\epsilon_{i}\leq 0,\end{cases}\qquad
I^{d}=\begin{cases}1&\text{if}\quad d<d_{n},\\\ 0&\text{if}\quad d\geq
d_{n}.\end{cases}$ (83)

The remaining part of the numerical procedure at this level encompasses the
standard steps typical of an FE-model. To this end, the bilinear shape
functions corresponding to a quadrilateral element and staggered solution
scheme are applied in the present work. This approach staggers between the
displacement and the phase-field and is advantageous compared to the
monolithic solution type due to its higher robustness. However, both
approaches are well established nowadays as shown in works applying monolithic
schemes (Miehe et al., 2010b; Kuhn and Müller, 2010; Schlueter et al., 2014;
Msekh et al., 2015) or using its counterpart (Miehe et al., 2010a; Borden et
al., 2012; Hofacker and Miehe, 2013). The standard Newton-Raphson procedure is
applied for the solution of nonlinear systems of equations.

###  8.3 Material point solution part

The material point solution part evaluates internal variables
$\boldsymbol{\epsilon}^{\mathrm{p}}$ and $\boldsymbol{\xi}$ by using a
predictor-corrector scheme if damage and deformations $(d,\boldsymbol{u})$ as
well as rate $\dot{\bar{\omega}}=\omega^{\prime}\,\frac{(d-d_{n})}{\Delta t}$
are known from the global level solution. The predictor step calculates the
trial stress by assuming that plastic deformations do not change in comparison
with the previous time step $n$

$\boldsymbol{\sigma}^{\mathrm{tr}}=\omega(d)\
\mathbb{C}\\!:\\!(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}}_{n})_{+}+\mathbb{C}\\!:\\!(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}}_{n})_{-}.$
(84)

The trial stress is furthermore introduced in the yield locus function to
check whether the prediction is true

$\Phi^{\mathrm{c,tr}}=||\boldsymbol{\bar{\sigma}}^{\mathrm{tr}}-\boldsymbol{\bar{\chi}}_{n}||-\bar{\omega}\
\bar{\sigma}_{\mathrm{Y}}.$ (85)

If the control value $\Phi^{\mathrm{c,tr}}$ is negative, the prediction (84)
is correct and the plastic deformations do not evolve in the present step.
Otherwise, a corrector step is needed to update the plastic deformations. The
corrector step firstly calculates the plastic multiplier according to the
scheme shown in Sect. 7. The initial value $\lambda_{0}$ is determined by
solving the simplified problem (71)

$\displaystyle\lambda_{0}=\frac{(\boldsymbol{\bar{\sigma}}^{\mathrm{tr}}-\boldsymbol{\bar{\chi}}_{n}):\left[\omega^{\prime}\,\boldsymbol{\bar{\sigma}}_{+}^{\mathrm{tr}}+\dot{\bar{\omega}}\,c\boldsymbol{\xi}_{n}+\frac{c}{b}\boldsymbol{\bar{\chi}}_{n}\right]}{c\bar{\omega}\,(\boldsymbol{\bar{\sigma}}^{\mathrm{tr}}-\boldsymbol{\bar{\chi}}_{n}):\bar{\boldsymbol{N}}^{\mathrm{tr}}},$
(86)

whereas the later iterations follow the rule

$\displaystyle\lambda_{i}=\frac{(\boldsymbol{\bar{\sigma}}^{\mathrm{tr}}-\boldsymbol{\bar{\chi}}_{n}):\boldsymbol{R}}{c\bar{\omega}\,(\boldsymbol{\bar{\sigma}}^{\mathrm{tr}}-\boldsymbol{\bar{\chi}}_{n}):\bar{\boldsymbol{N}}^{\mathrm{tr}}},$
(87)
$\displaystyle\boldsymbol{R}=\left[\omega^{\prime}\,\boldsymbol{\sigma}_{+}^{\mathrm{tr}}+\omega\,\mathbb{C}\\!:\\![\frac{({\boldsymbol{\epsilon}}-{\boldsymbol{\epsilon}}_{n})}{\Delta
t}-\lambda_{i-1}\bar{\boldsymbol{N}}^{\mathrm{tr}}]_{+}\right.\left.+\,\mathbb{C}\\!:\\![\frac{({\boldsymbol{\epsilon}}-{\boldsymbol{\epsilon}}_{n})}{\Delta
t}-\lambda_{i-1}\bar{\boldsymbol{N}}^{\mathrm{tr}}]_{-}+\dot{\bar{\omega}}\,c\boldsymbol{\xi}_{n}+\frac{c}{b}\,\boldsymbol{\bar{\chi}}_{n}\right].$
(88)

The previous expressions use the approximation
$\dot{\boldsymbol{\epsilon}}=\frac{({\boldsymbol{\epsilon}}-{\boldsymbol{\epsilon}}_{n})}{\Delta
t}$ and notation
$\bar{\boldsymbol{N}}^{\mathrm{tr}}=\frac{(\boldsymbol{\bar{\sigma}}^{\mathrm{tr}}-\boldsymbol{\bar{\chi}}^{\mathrm{n}})}{\lvert\lvert\boldsymbol{\bar{\sigma}}^{\mathrm{tr}}-\boldsymbol{\bar{\chi}}^{\mathrm{n}}\rvert\rvert}$.
Note that index $i$ is related to the iterative solution for $\lambda_{i}$,
whereas index $n$ denotes the time step. Quantities without any index are
related to the step $n+1$ or are constant. Finally, the solution of the
iterative procedure (87)-(88) is denoted $\lambda=\lambda_{i}$ and used to
update plastic deformations $\boldsymbol{\epsilon}^{\mathrm{p}}$, driving
forces $\boldsymbol{\chi}$ and stress response $\boldsymbol{\sigma}$

$\displaystyle\dot{\boldsymbol{\epsilon}}_{\mathrm{p}}=\frac{\boldsymbol{\epsilon}^{\mathrm{p}}-\boldsymbol{\epsilon}^{\mathrm{p}}_{n}}{\Delta\mathrm{t}}=\lambda\bar{\boldsymbol{N}}^{\mathrm{tr}}\hskip
79.09866pt\qquad\Rightarrow\quad\boldsymbol{\epsilon}^{\mathrm{p}}=\Delta
t\lambda\bar{\boldsymbol{N}}^{\mathrm{tr}}+\boldsymbol{\epsilon}^{\mathrm{p}}_{n},$
(89)
$\displaystyle\dot{\boldsymbol{\xi}}=\frac{\boldsymbol{\xi}-\boldsymbol{\xi}_{n}}{\Delta
t}=\frac{1}{b\bar{\omega}}\boldsymbol{\bar{\chi}}_{n}-\dot{\boldsymbol{\epsilon}}_{\mathrm{p}}\hskip
62.59605pt\qquad\Rightarrow\quad\boldsymbol{\xi}=\frac{\Delta
t}{b\bar{\omega}}\boldsymbol{\bar{\chi}}_{n}-\Delta
t\,\lambda\bar{\boldsymbol{N}}^{\mathrm{tr}}+\boldsymbol{\xi}_{n},$ (90)
$\displaystyle\dot{\boldsymbol{\chi}}=\dot{\bar{\boldsymbol{\chi}}}=\frac{\boldsymbol{\bar{\chi}}-\boldsymbol{\bar{\chi}}_{n}}{\Delta\mathrm{t}}=-\dot{\bar{\omega}}\,c\boldsymbol{\xi}_{n}+c{\bar{\omega}}\dot{\boldsymbol{\epsilon}}_{\mathrm{p}}-\frac{c}{b}\,\boldsymbol{\chi}_{n}\quad\Rightarrow\quad\boldsymbol{\chi}=\boldsymbol{\bar{\chi}}=-\Delta
t\,\dot{\bar{\omega}}\,c\boldsymbol{\xi}_{n}+c\bar{\omega}\,\Delta
t\,\lambda\bar{\boldsymbol{N}}^{\mathrm{tr}}-\frac{c\Delta
t}{b}\boldsymbol{\bar{\chi}}_{n}+\boldsymbol{\bar{\chi}}_{n},$ (91)
$\displaystyle\boldsymbol{\sigma}_{\mathrm{}}=\omega(d)\
\mathbb{C}\\!:\\!(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}})_{+}+\mathbb{C}\\!:\\!(\boldsymbol{\epsilon}-\boldsymbol{\epsilon}^{\mathrm{p}})_{-}.$
(92)

##  9 Representative numerical examples

The model developed inherits the advantages of both methods that it
incorporates. That makes its application field large, as demonstrated in the
subsequent sections studying the behavior of the cold-formed carbon steel (CS)
and of the cold-formed stainless steel (SS). Among others, selected numerical
examples deal with the purely elasto-plastic material behavior (Sect. 9.1),
with the simulation of crack propagation on a notched sample (Sect. 9.2) and
with the study of the life time in the LCF-mode (Sect.9.3).

###  9.1 Elasto-plastic behavior of the carbon steel and of the stainless
steel

The first group of examples illustrating the application of the model studies
a purely elasto-plastic material behavior under a cyclic load by neglecting
the damage effects. To this end, simulations at a single material point are
performed for parameter sets shown in Tab. 1. The prescribed strains gradually
change from zero to 1 % and thereafter decrease to -1 %. The cycle again
closes at zero load level. Two cycles with identical load paths are carried
out. Each cycle takes 2 s to complete, whereas a single time step is $\Delta
t=0.01$ s.

Table 1: Elastic and plastic material parameters for CS and SS. Parameter | Value CS | Value SS | Unit | Name
---|---|---|---|---
$E$ | 212910 | 197960 | [MPa] | Young’s modulus
$\nu$ | 0.3 | 0.3 | [–] | Poisson’s ratio
$\sigma_{\mathrm{Y}}$ | 451 | 552 | [MPa] | yield stress
$b$ | 5 | 6 | [Ns/mm2] | pseudo-viscosity
$c$ | 30000 | 70000 | [N/mm] | kinematic hardening parameter

Simulations are carried out for two kinds of steel, namely for the cold-formed
carbon steel (CS) and the cold-formed stainless steel (SS). The results are
furthermore compared to the experimental findings by Nip et al. (2010). This
comparison is presented in Fig. 1 and approves an excellent agreement. The
stainless steel hardens stronger than the carbon steel, such that parameter
$c$ is higher in this case. More details on the parameter identification in
the context of the Armstrong-Frederick model are for example provided in Wolff
et al. (2010); Lubarda and Benson (2002).

Figure 1: Comparison of experimental and numerical results at a material
point. (a) Stress-strain hysteresis for carbon steel (CS). (b) Stress-strain
hysteresis for stainless steel (SS). Plotted are the 11-components of the
stress tensor and the strain tensor. Experimental results are taken from Nip
et al. (2010).

###  9.2 Crack propagation on a notched sample

The geometry of the sample for the second group of simulations is shown in
Fig. 2. The chosen square plate has dimensions 1 mm $\times$ 1 mm and a
thickness of 1 mm and represents a cutout of a 3D body with the large
thickness, which substantiates the assumption of plane strains (Fig. 2a). The
model does not include any material length scale, such that sample size and
prescribed displacements can be scaled in a straightforward manner. The left
edge of the sample is fixed in horizontal direction. The vertical displacement
at the bottom left corner is additionally suppressed. Horizontal displacements
$\bar{u}_{x}$ are prescribed at the right edge. The plate has a vertical
initial crack from the middle of the bottom edge to the midpoint of the
sample. It is discretized by a mesh with approximately 22500 quadrilateral
elements such that a fine discretization is performed in the areas where crack
propagation is expected (Fig. 2b). In this area, the effective element size
$h$ fulfills the condition $h\approx 0.001\text{ mm}<l/2$ and is significantly
less than the minimum size required to achieve reasonable accuracy in the
crack zone (Miehe et al., 2010b). Here, $l$ denotes the characteristic crack
width typical of the phase field method.

Figure 2: (a) Geometry of the notched plate with prescribed horizontal
displacements $\bar{u}_{x}$. (b) Discretization of the plate with
approximately 22500 elements.

The section simulates the crack propagation on the notched example and
investigates the influence of material parameters on this process. The carbon
steel is assumed as the representative material this time. Apart from elastic
and plastic material parameters (Tab. 1), simulations require damage
parameters. The values chosen to this end are listed in Tab. 2. Here, the
values proposed by Miehe et al. (2010b) are used for the critical energy
$g_{c}$ and characteristic crack width $l$. Constant $k_{\mathrm{p}}$ of the
penalty function (Eq. (4)) is chosen to be $1\times 10^{8}$, which is
sufficiently large to suppress the negative damage rate. Constant
$k_{\mathrm{d}}$ has the value of $1\times 10^{-8}$ and prevents the energy
from becoming identical to zero in the case of full material damage. In all
examples of this group, the load increment and damage parameters $g_{c}$ and
$l$ are assumed such that the crack formation and propagation can be simulated
within a fairly small number of time steps in order to reduce the
computational effort and to provide a qualitative analysis of relevant
phenomena according to the results from Miehe et al. (2010b).

Table 2: Material parameters of the damage model. Parameter | Value | Unit | Name
---|---|---|---
$g_{c}$ | 0.27 | [N/mm2] | critical energy release rate
$l$ | 0.0375 | [mm] | characteristic crack width
$k_{\mathrm{p}}$ | $1\times 10^{8}$ | [–] | penalty constant
$k_{\mathrm{d}}$ | $1\times 10^{-8}$ | [–] | numerical constant

The behavior of the notched sample (Fig. 2) is investigated for three load
constellations: uniformly increasing tension, cyclic load with a tension and
compression phase and cyclic load with an increasing amplitude in the tension
regime.

In the first case, the prescribed displacement linearly increases up to the
maximum value of $\bar{u}_{x}=0.005$ mm. The constant time increment is
$\Delta t~{}=~{}0.005$ s and the total loading time amounts to 1 s. The
results at a quadrature point of an element directly located at the end of the
initial crack are monitored for the illustration in Fig. 3. This figure shows
the change of damage and stress state during tension tests for different pairs
of parameters $b$ and $c$. Here, rapidly hardening materials with a high
hardening modulus $c$ lead to a faster increase of damage, whereas the higher
pseudo-viscosity slows down the damage evolution and postpones softening.
Damage increases from moderate values to the maximum, such that the material
loses its strength and stress falls to the nearly zero value. The stress drop
takes place in approximately 0.05 s. Simulations are performed for different
time increments in order to check the accuracy of the solution. For a time
increment $\Delta t~{}=~{}0.005$ s and less, identical results are achieved.

Figure 3: Results of tensile tests for different parameter pairs $b$ and $c$.
Prescribed displacements uniformly change up to the value of 0.005 mm. (a)
Applied load as a function of time and damage evolution over time. (b) The
11-component of the stress tensor over the time. Hardening modulus $c$ is
expressed in N/mm and pseudo-viscosity $b$ is expressed in Ns/mm2.

The same setup (Fig. 2) is used to investigate influences of a cyclic load as
shown in Fig. 4. Here, two load cycles are performed with the total duration
of 4 s. The time increment is $\Delta t=0.005$ s. The horizontal displacement
$\bar{u}_{x}$ changes in the range [0.0025 mm, -0.0025 mm], implying that a
cycle includes both: the tension and the compression mode. The load amplitudes
are constant in both cycles. Figure 4a shows the load path and the damage
evolution. The damage variable evolves during the tension phase, whereas the
unloading and the compression mode do not affect it. This behavior goes back
to the split of the energy into a tension and a compression part (Eq. (6)).
Even though the same displacement is applied in every cycle, an increase in
damage is observed during each loading cycle. Figure 4b shows the strong
influence of damage on the corresponding stress-strain hysteresis. Here, the
Baushinger effect is hardly noticeable, although a displacement controlled
test with constant amplitudes is simulated. Fig. 4b depicts the ratcheting
rather than the Bauschinger effect.

Figure 4: (a) Applied load as a function of time and the damage evolution for
two load cycles. (b) The 11-component of the stress tensor over the
11-component of the strain tensor.

The last load constellation applied to the setup from Fig. 2 deals with the
cyclic load in the tension regime with the increasing amplitude (Fig. 5). The
amplitude increment amounts to $2\times 10^{-4}$ mm per cycle. Time increment
$\Delta t=0.005$ s is assumed. Three damage contour plots are chosen
presenting states of the crack during propagation for parameters $c=15000$
N/mm (Fig. 6a) and $c=30000$ N/mm (Fig. 6b). The first snapshot indicates the
state at which damage variable $d$ reaches a value of one at the already
existing crack tip. The second plot is taken when the crack has propagated
roughly half the way through the plate, and the last plot is taken when the
plate is fully torn in half. For the first parameter value ($c=15000$ N/mm),
the crack opening starts in the 19th load cycle and is fully propagated
throughout the plate in the 23rd load cycle (Fig. 6a). For the second
parameter value ($c=30000$ N/mm), damage evolves earlier and faster: here, the
complete crack propagation takes place within three load cycles (Fig. 6b).

Figure 5: Increasing load amplitudes as a function of time.  Figure 6: (a)
Crack propagation takes place in five cycles ($c=15000$ N/mm, $b=2$ Ns/mm2).
(b) Crack propagation takes place in three cycles ($c=30000$ N/mm, $b=2$
Ns/mm2).

###  9.3 Life time of the cold steel and of the stainless steel in the LCF-
mode

Amongst others, the model proposed enables the estimation of the life time in
the LCF-mode. To this end, the number of cycles up to the total failure is
evaluated dependent on the amplitude of the cyclic load applied. The latter is
kept constant during the test. The material failure is caused by the material
defects such as production process induced pores at a small length scale. In
most cases, surface porosity is the critical factor for the fatigue
phenomenon. However, the present paper assumes that different methods of the
surface treatment can resolve this type of imperfection and that a defect in
the bulk of material causes the crack initiation.

This part of the analysis assumes the sample geometry (Fig. 7a) proposed by
Nip et al. (2010). The horizontal displacements are constrained at the left
boundary and the cyclic horizontal displacements are prescribed at the right
boundary. In addition, the vertical displacement at the lower left corner is
constrained to suppress the rigid body motion. For the chosen boundary
condition, the largest deformations are expected in the center of the sample,
which is discretized by a fine mesh. The effective element size in this area
is h $\approx$ 0.25 mm, such that the crack width length $l=0.6$ is chosen.
The critical energy is set to $g_{c}=100$ N/mm2 which corresponds to a ductile
material behavior, whereas the remaining material parameters are kept as in
Sects. 9.1 and 9.2. Simulations are additionally performed for higher values
of the hardening modulus $c$ in order to study the influence of this
parameter. A defect is created in an element in the center of the specimen,
which initiates crack growth in the otherwise homogeneous material. In this
element, the values of the Young’s modulus and the yield stress are set to 75%
of the actual material parameters.

The very first step of the analysis applies a strain of 3 %. The assumed time
increment is $\Delta t=4.28\times 10^{-4}$ s and 2000 time steps simulate a
single cycle. As expected, a vertical crack propagates through the middle of
the sample, which is shown in Fig. 7b. Here, an early stage of the crack is
already noticeable in the 47th cycle, whereas the total failure occurs in
cycle number 96. Figure 7b corresponds to the carbon steel with the hardening
modulus $c=30000$ N/mm.

Figure 7: (a) Geometry of the simulated sample and prescribed boundary
conditions. (b) Different stages of the crack propagation under the cyclic
load. Simulations are performed for the carbon steel with the hardening
modulus $c=30000$ N/mm.

Furthermore, the tests are repeated by increasing the load amplitude to the
strain of 4 %, 5 %, 6 % and 7 %. The same setup (Fig. 7a) is applied to this
end. The number of time steps within a cycle is kept constant (2000), however,
the time increment varies in the range from $\Delta t=4.28\times 10^{-4}$ s to
$\Delta t=1\times 10^{-3}$ s. Accordingly, the duration of one cycle takes the
values between 0.856 s and 2 s.

The results of simulations together with the experimental results by Nip et
al. (2010) are presented in Fig. 8. First, the behavior of carbon steel is
studied and compared to the Coffin-Manson curve relating the applied strains
to the number of cycles, both in logarithmic scales (Fig. 8a). Numerical
results show an excellent agreement with the experimental findings, in
particular for higher strains. The dependence between the applied strains and
the number of cycles is linear, however, the number of cycles is minimally
overestimated. The discrepancy between the experimental and numerical results
slightly increases with the decreasing strains. The same kind of simulations
is repeated for a higher value of the hardening parameter, which leads to a
significant decrease of the number of cycles up to the failure. This tendency
goes back to the fact that a higher hardening modulus is related to the faster
damage evolution and consequently causes a reduction of the number of cycles.
Simulations for a higher hardening modulus approves the linear dependency
between data, however, they indicate that a change of the strain amplitude has
a higher influence on the change of number of cycles in this case.

The observations mentioned previously have also been approved by repeating the
same kind of tests for the stainless steel as shown in Fig. 8b. Here, the
Manson-Coffin curve has a slightly higher slope, but numerical results
compared to the experimental ones show the same tendencies as for the carbon
steel.

Figure 8: Comparison of the experimental Coffin-Manson curve (Nip et al.,
2010) to the numerical results. (a) Results for carbon steel with two
different hardening parameters $c=30000$ N/mm and $c=300000$ N/mm. (b) Results
for the stainless steel with the hardening parameters $c=70000$ N/mm and
$c=700000$ N/mm.

##  10 Conclusions and outlook

The present work couples the phase field method of fracture to the Armstrong-
Frederick model of plasticity with the kinematic hardening. The chosen concept
inherits the advantages of both techniques and is aimed at the study of LCF
effects in ductile materials. However, the numerical implementation of this
promising approach faces several challenges, such as the definition of a
unique framework for both setups, the derivation of coupled evolution
equations, the distinction between tension and compression mode and the
development of a computationally efficient algorithm. The basis for the phase
field fracture model are the Griffith’s theory and the dissipation potential
relying on the assumption of a crack surface function. The derivation of
evolution equations uses the minimum principle of the dissipation potential,
which requires to express the dissipation potential of the classical
Armstrong-Frederick model in terms of the internal variable rates by using the
Legendre transformation. The model developed also takes into account that the
unloading and the compression mode have no effect on the damage evolution. The
approach is eventually implemented in the FE-program FEAP where the
displacements and damage are calculated at the global level by using a
staggered scheme, whereas inelastic internal variables are updated locally, at
each Gauss point.

The application of the model is illustrated by three groups of examples
related to the cold-formed carbon steel and the cold-formed stainless steel.
The first group of tests investigates the purely elasto-plastic material
behavior without the influence of damage and focuses on the numerical
reconstruction of the stress-strain hysteresis due to a cyclic load. The
second group of tests simulates the crack propagation on a notched sample and
particularly analyses the influence of plastic parameters on the damage
evolution. The simulations show that the rapidly hardening materials are prone
to a faster damage evolution. The ability of the model to constrain the damage
evolution in the unloading and compression stages of a loading cycle is
demonstrated as well. The last group of tests studies the life time of the
carbon steel in the LCF-mode and compares the numerical results to the
experimental findings by Nip et al. (2010). The results show an excellent
qualitative and quantitative agreement, and the linear dependency between the
data in a logarithmic scale is approved.

The model proposed is a promising tool with regard to the simulation of
fatigue effects, giving rise to many new issues. Among others, it can be
extended to capture “hidden” aspects of the fatigue process such as the
microcrack initiation based on the theory of persistent slip bands and the
microcrack propagation along the crystallographic planes. Moreover, the
additional effects typical of high and very high cycle fatigue, as well as the
threshold for the damage initiation have to be incorporated. A comprehensive
validation of the model with respect to the experimental results as well as
its application and calibration for further materials are also envisaged.

##  Appendix A Legendre transformation

The formulation of the MDP in terms of rates of internal variables is based on
Legendre transformation:

$\displaystyle\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}})=\underset{{\boldsymbol{q}}}{\text{max}}\\{\mathcal{L}^{\mathrm{LT}}=\boldsymbol{q}:\dot{\boldsymbol{\nu}}-\Phi^{\mathrm{AF^{*}}}(\boldsymbol{\nu},{\boldsymbol{q}})\\},$
(93)
$\displaystyle\mathcal{L}^{\mathrm{LT}}=\boldsymbol{\sigma}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+\boldsymbol{\chi}:\dot{\boldsymbol{\xi}}-a||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||-\frac{1}{2b}||\boldsymbol{\bar{\chi}}||^{2}.$
(94)

The corresponding stationary conditions have already been derived in Sect.5
(Eqs.(34)b and(35)b)

$\boldsymbol{\bar{\sigma}}=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}+b\dot{\boldsymbol{\xi}}+b\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}},\qquad\boldsymbol{\bar{\chi}}=b\dot{\boldsymbol{\xi}}+b\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}.$
(95)

Their implementation in (94) now yields

$\displaystyle\mathcal{L}^{\mathrm{LT}}=\boldsymbol{\sigma}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+\boldsymbol{\chi}:\dot{\boldsymbol{\xi}}-a||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||-\frac{1}{2b}||\boldsymbol{\bar{\chi}}||^{2}=$
$\displaystyle=\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+b\dot{\boldsymbol{\xi}}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+b\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+b\dot{\boldsymbol{\xi}}:\dot{\boldsymbol{\xi}}+b\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}:\dot{\boldsymbol{\xi}}-a||\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}||-\frac{1}{2b}\left[b^{2}||\dot{\boldsymbol{\xi}}||^{2}+2b^{2}\dot{\boldsymbol{\xi}}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+b^{2}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||^{2}\right]$
$\displaystyle=\bar{\sigma}_{\mathrm{Y}}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||+b\dot{\boldsymbol{\xi}}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}+b||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||^{2}+b||\dot{\boldsymbol{\xi}}||^{2}+b\dot{\boldsymbol{\xi}}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}-a\bar{\sigma}_{\mathrm{Y}}-\frac{b}{2}||\dot{\boldsymbol{\xi}}||^{2}-b\dot{\boldsymbol{\xi}}:\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}-\frac{b}{2}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||^{2}$
$\displaystyle=\bar{\sigma}_{\mathrm{Y}}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||+\frac{b}{2}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||^{2}+\frac{b}{2}||\dot{\boldsymbol{\xi}}||^{2}+b\,\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}\\!:\\!\dot{\boldsymbol{\xi}}-a\bar{\sigma}_{\mathrm{Y}}.$
(96)

The last term is a constant and thus can be neglected, since the evolution
equations are obtained by its minimization of the potential. Accordingly, the
sought potential turns into

$\Phi^{\mathrm{AF}}(\boldsymbol{\nu},\dot{\boldsymbol{\nu}})=\bar{\sigma}_{\mathrm{Y}}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||+\frac{b}{2}||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||^{2}+\frac{b}{2}||\dot{\boldsymbol{\xi}}||^{2}+b\,\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}\\!:\\!\dot{\boldsymbol{\xi}}.$
(97)

##  Appendix B Evolution equations of the coupled problem

The potential of the coupled problem has the form

$\mathcal{L}^{\mathrm{MDP}}=\dot{\Psi}^{\mathrm{c}}+\Phi^{\mathrm{c}}=\dot{\Psi}^{\mathrm{c}}+\Phi^{\mathrm{d}}+\bar{\omega}\Phi^{\mathrm{AF}},$
(98)

where $\dot{\Psi}^{\mathrm{c}}$ is defined in Eq. (52), $\Phi^{\mathrm{d}}$ is
defined by Eq. (5) and $\Phi^{\mathrm{AF}}$ by Eq. (37). Driving forces
corresponding to ${\boldsymbol{\epsilon}^{\mathrm{p}}}$ and
${\boldsymbol{\xi}}$ are then calculated as following derivatives

$\displaystyle\boldsymbol{q}_{\boldsymbol{\epsilon}^{\mathrm{p}}}=\boldsymbol{\bar{\sigma}}=\frac{\partial\Phi^{\mathrm{c}}}{\partial{\dot{\boldsymbol{\epsilon}}^{\mathrm{p}}}}=\bar{\omega}\left[\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}+b{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}+b{\dot{\boldsymbol{{\xi}}}}\right],$
(99)
$\displaystyle\boldsymbol{q}_{\boldsymbol{\xi}}=\boldsymbol{\bar{\chi}}=\frac{\partial\Phi^{\mathrm{c}}}{\partial{\dot{\boldsymbol{\xi}}}}=\bar{\omega}\left[b{\dot{\boldsymbol{{\xi}}}}+b{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}\right],$
(100)

such that a transformation of the system (99) and (100) yields the evolution
equations

${\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}=\lambda\frac{\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}}{||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||},\qquad\dot{\boldsymbol{\xi}}=-{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}+\frac{1}{\bar{\omega}b}\boldsymbol{\bar{\chi}},$
(101)

whereas the insertion of (100) into (99) provides the yield locus equation

$\boldsymbol{\bar{\sigma}}=\bar{\omega}\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}+\boldsymbol{\chi}\qquad\Rightarrow\qquad\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}=\bar{\omega}\bar{\sigma}_{\mathrm{Y}}\frac{\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}}{||\dot{\boldsymbol{{\epsilon}}}^{\mathrm{p}}||}\qquad\Rightarrow\qquad||\boldsymbol{\bar{\sigma}}-\boldsymbol{\bar{\chi}}||=\bar{\omega}\bar{\sigma}_{\mathrm{Y}}.$
(102)

Finally, constitutive law (55) is used to derive the evolution equation for
the back stress:

$\dot{\boldsymbol{\chi}}=-\dot{\bar{\omega}}c\boldsymbol{\xi}-\bar{\omega}c\dot{\boldsymbol{\xi}}=-\dot{\bar{\omega}}c\boldsymbol{\xi}+c\bar{\omega}\left(\dot{\boldsymbol{\epsilon}^{\mathrm{p}}}-\
\frac{1}{b\bar{\omega}}\ \boldsymbol{\bar{\chi}}\right).$ (103)

## Acknowledgments

The authors are thankful to Frank Walther and Mustafa Awd (TU Dortmund
University, Dortmund, Germany) for their helpful discussions on the fatigue
behavior of metal materials.

## References

  * Abdollahi and Arias (2012) Abdollahi, A., Arias, I., 2012\.  Phase-field modeling of crack propagation in piezoelectric and ferroelectric materials with different electromechanical crack conditions.  J. Mech. Phys. Solids 60, 2100–2126.
  * Alessi et al. (2014) Alessi, R., Marigo, J.J., Vidoli, S., 2014.  Gradient damage models coupled with plasticity and nucleation of cohesive cracks.  Arch. Rat. Mech. Anal. 214, 575–615.
  * Ambati et al. (2015a) Ambati, M., Gerasimov, T., De Lorenzis, L., 2015a.  Phase-field modeling of ductile fracture.  Comput. Mech. 55, 1017–1040.
  * Ambati et al. (2015b) Ambati, M., Kruse, R., De Lorenzis, L., 2015b.  A phase-field model for ductile fracture at finite strains and its experimental verification.  Comput. Mech. 57, 149–167.
  * Ambrosio and Tortorelli (1990) Ambrosio, L., Tortorelli, V.M., 1990\.  Approximation of functionals depending on jumps by elliptic functionals via t-convergence.  Commun- Pur. Appl. Math. 43, 999–1036.
  * Ambrosio and Tortorelli (1992) Ambrosio, L., Tortorelli, V.M., 1992\.  On the approximation of free discontinuity problems.  Boll. Un. Mat. Ital. B(7) 6, 105–123.
  * Armstrong and Frederick (1966) Armstrong, P.J., Frederick, C.O., 1966\.  A mathematical representation of the multiaxial Bauschinger effect.  Central Electricity Generating Board and Berkeley Nuclear Laboratories 731.
  * Aygün and Klinge (2020) Aygün, S., Klinge, S., 2020\.  Continuum mechanical modeling of strain-induced crystallization in polymers.  Int. J. Solids Struct. 196–197, 129–139.
  * Bari and Hassan (2001) Bari, S., Hassan, T., 2001.  Kinematic hardening rules in uncoupled modeling for multiaxial ratcheting simulation.  Int. J. Plast. 17, 885–905.
  * Borden et al. (2012) Borden, M., Verhoosel, C., Scott, M., Hughes, T., Landis, C., 2012.  A phase-field description of dynamic brittle fracture.  Comput. Methods Appl. Mech. Eng. 217–220, 77–95.
  * Borden et al. (2016) 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 derivation, plastic degradation, and stress triaxiality effects.  Comput. Methods Appl. Mech. Eng. 312, 130–166.
  * Bourdin and Chambolle (2000) Bourdin, B., Chambolle, A., 2000\.  Implementation of an adaptive finite-element approximation of the mumford-shah functional.  Numer. Math. 85, 609–646.
  * Bourdin et al. (2012) Bourdin, B., Chukwudozie, C., Yoshioka, K., 2012.  A variational approach to the numerical simulation of hydraulic fracturing, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. pp. 1–9.
  * Bourdin et al. (2000) Bourdin, B., Francfort, G.A., Marigo, J.J., 2000.  Numerical experiments in revisited brittle fracture.  J. Mech. Phys. Solids 48, 797–826.
  * Bourdin et al. (2014) Bourdin, B., Marigo, J.J., Maurini, C., Sicsic, P., 2014\.  Morphogenesis and propagation of complex cracks induced by thermal shocks.  Phys. Rev. Lett. 112, 014301\.
  * Crismale and Lazzaroni (2016) Crismale, V., Lazzaroni, G., 2016\.  Viscous approximation of quasistatic evolutions for a coupled elastoplastic-damage model.  Calc. Var. Partial Dif. 55, 1–54.
  * Dettmer and Reese (2004) Dettmer, W., Reese, S., 2004\.  On the theoretical and numerical modelling of Armstrong-Frederick kinematic hardening in the finite strain regime.  Comput. Methods Appl. Mech. Eng. 193, 87–116.
  * Francfort and Marigo (1998) Francfort, G.A., Marigo, J.J., 1998\.  Revisiting brittle fracture as an energy minimization problem.  J. Mech. Phys. Solids. 46, 1319–1342.
  * Fraternali (2007) Fraternali, F., 2007.  Free discontinuity finite element models in two-dimensions for in-plane crack problems.  Theor. Appl. Fract. Mec. 47, 274–282.
  * Freddi and Iurlano (2017) Freddi, F., Iurlano, F., 2017\.  Numerical insight of a variational smeared approach to cohesive fracture.  J. Mech. Phys. Solids 98, 156–171.
  * Hofacker and Miehe (2013) Hofacker, M., Miehe, C., 2013\.  A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns.  Int. J. Numer. Methods Eng. 93, 276–301.
  * Karma et al. (2001) Karma, A., Kessler, D.A., Levine, H., 2001.  Phase-field model of mode III dynamic fracture.  Phys. Rev. Lett. 87, 045501\.
  * Khan and Jackson (1999) Khan, A.S., Jackson, K.M., 1999\.  On the evolution of isotropic and kinematic hardening with finite plastic deformation, Part I: Compression/tension loading of OFHC copper cylinders.  Int. J. Plast. 15, 1265–1275.
  * Kobayashi and Ohno (2002) Kobayashi, M., Ohno, N., 2002\.  Implementation of cyclic plasticity models based on a general form of kinematic hardening.  Int. J. Numer. Meth. Eng. 53, 2217–2238.
  * Kuhn and Müller (2010) Kuhn, C., Müller, R., 2010\.  A continuum phase field model for fracture.  Eng. Fract. Mech. 77, 3625–3634.  Computational Mechanics in Fracture and Damage: A Special Issue in Honor of Prof. Gross.
  * Li et al. (2014) Li, B., Peco, C., Millan, D., Arias, I., Arroyo, M., 2014\.  Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy.  Int. J. Numer. Meth. Eng. 102, 711–727.
  * Lion (2000) Lion, A., 2000.  Constitutive modelling in finite thermoviscoplasticity: A physical approach based on nonlinear rheological models.  Int. J. Plast. 16, 469–494.
  * Lubarda and Benson (2002) Lubarda, V., Benson, D., 2002\.  On the numerical algorithm for isotropic–kinematic hardening with the Armstrong–Frederick evolution of the back stress.  Comput. Methods Appl. Mech. Eng. 191, 3583–3596.
  * Lührs et al. (1997) Lührs, G., Hartmann, S., Haupt, P., 1997.  On the numerical treatment of finite deformations in elastoviscoplasticity.  Comput. Methods Appl. Mech. Eng. 144, 1–21.
  * Maurini et al. (2013) Maurini, C., Bourdin, B., Gauthier, G., V., L., 2013\.  Crack patterns obtained by unidirectional drying of a colloidal suspension in a capillary tube: Experiments and numerical simulations using a two-dimensional variational approach.  Int. J. Fracture 184, 75–91.
  * Miehe et al. (2016) Miehe, C., Aldakheel, F., Raina, A., 2016.  Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory.  Int. J. Plast. 84, 1–32.
  * Miehe et al. (2015) Miehe, C., Hofacker, M., Schänzel, L.M., Aldakheel, F., 2015\.  Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids.  Comput. Methods Appl. Mech. Eng. 294, 486–522.
  * Miehe et al. (2010a) Miehe, C., Hofacker, M., Welschinger, F., 2010a.  A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits.  Comput. Methods Appl. Mech. Eng. 199, 2765–2778.
  * Miehe et al. (2010b) Miehe, C., Welschinger, F., Hofacker, M., 2010b.  Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations.  Int. J. Numer. Meth. Eng. 83, 1273–1311.
  * Mollica et al. (2001) Mollica, F., Rajagopal, K.R., Srinivasa, A.R., 2001.  The inelastic behavior of metals subject to loading reversal.  Int. J. Plast. 17, 1119–1146.
  * Msekh et al. (2015) Msekh, M., Sargado, J.M., Jamshidian, M., Areias, P., Rabczuk, T., 2015.  Abaqus implementation of phase-field model for brittle fracture.  Comput. Mater. Sci. 96, 472–484.
  * Narayan and Anand (2019) Narayan, S., Anand, L., 2019\.  A gradient damage theory of fracture of quasi-brittle materials.  J. Mech. Phys. Solids 129, 119–146.
  * Negri and Paolini (2001) Negri, M., Paolini, M., 2001\.  Numerical minimization of the mumford-shah functional.  Calcolo 38, 67–84.
  * Nip et al. (2010) Nip, K.H., Gardner, L., Davies, C.M., Elghazouli, A.Y., 2010\.  Extremely low cycle fatigue tests on structural carbon steel and stainless steel.  J. Constr. Steel Res. 66, 96–110.
  * Ohno and Wang (1993) Ohno, N., Wang, J.D., 1993.  Kinematic hardening rules with critical state of dynamic recovery: Part I: Formulation and basic features for ratcheting behavior, Part II: Application to experiments of ratcheting behavior.  Int. J. Plast. 9, 375–403.
  * Puzrin and Houlsby (2001) Puzrin, A.M., Houlsby, G.T., 2001\.  Fundamentals of kinematic hardening hyperplasticity.  Int. J. Solids Struct. 38, 3771–3794.
  * Schlueter et al. (2014) Schlueter, A., Willenbücher, A., Kuhn, C., Müller, R., 2014\.  Phase field approximation of dynamic brittle fracture.  Comput. Mech. 54, 1141–1161.
  * Schmidt et al. (2009) Schmidt, B., Fraternali, F., Ortiz, M., 2009.  Eigenfracture: An eigendeformation approach to variational fracture.  Multiscale Model. Sim. 7, 1237–1266.
  * Svendsen et al. (1998) Svendsen, B., Arndt, S., Klingbeil, D., Sievert, R., 1998\.  Hyperelastic models for elastoplasticity with nonlinear isotropic and kinematic hardening at large deformation.  Int. J. Solids Struct. 35, 3363–3389.
  * Wheeler et al. (2014) Wheeler, M.F., Wick, T., Wollner, W., 2014.  An augmented-lagrangian method for the phase-field approach for pressurized fractures.  Comput. Methods Appl. Mech. Eng. 271, 69–85.
  * Wilson et al. (2013) Wilson, Z.A., Borden, M.J., Landis, C.M., 2013.  A phase-field model for fracture in piezoelectric ceramics.  Int. J. Fracture 193, 135–153.
  * Wilson and Landis (2016) Wilson, Z.A., Landis, C.M., 2016\.  Phase-field modeling of hydraulic fracture.  J. Mech. Phys. Solids 96, 264–290.
  * Wolff et al. (2010) Wolff, M., Suhr, B., Simsir, C., 2010.  Parameter identification for an Armstrong–Frederick hardening law for supercooled austenite of SAE 52100 steel.  Comput. Mater. Sci. 50, 487–495.