RÃ©gularisation de Lipschitz pour les modÃ©les de matÃ©riaux adoucissants :
l’approche Lip-field N<EMAIL_ADDRESS>1Ecole Centrale de
Nantes, GeM Institute, UMR CNRS 6183, 1 rue de la NoÃ«, 44321 Nantes, France
N<EMAIL_ADDRESS>
# Lipschitz regularization for softening material models: the Lip-field
approach

Nicolas Moës  Ecole Centrale de Nantes, GeM Institute, UMR CNRS 6183, 1 rue de
la NoÃ«, 44321 Nantes, France  Institut Universitaire de France (IUF)  [
Nicolas Chevaugeon  [

###### Abstract

Softening material models are known to trigger spurious localizations. This
may be shown theoretically by the existence of solutions with zero dissipation
when localization occurs and numerically with spurious mesh dependency and
localization in a single layer of elements. We introduce in this paper a new
way to avoid spurious localization. The idea is to enforce a Lipschitz
regularity on the internal variables responsible for the material softening.
The regularity constraint introduces the needed length scale in the material
formulation. Moreover, we prove bounds on the domain affected by this
constraint. A first one-dimensional finite element implementation is proposed
for softening elasticity and softening plasticity.

###### keywords:

softening, localization, damage, plasticity, Lipschitz, Lip-field

Les modÃ¨les de matÃ©riaux adoucissants sont connus pour dÃ©clencher des
localisations parasites. Cela peut Ãªtre dÃ©montrÃ© thÃ©oriquement par
l’existence de solutions avec une dissipation nulle lors de la localisation et
numÃ©riquement avec une dÃ©pendance de maillage et une localisation dans une
seule couche d’Ã©lÃ©ments. Nous introduisons dans cet article une nouvelle
faÃ§on d’Ã©viter les localisations parasites. L’idÃ©e est d’imposer une
rÃ©gularitÃ© de Lipschitz sur les variables internes responsables de
l’adoucissement. La contrainte de rÃ©gularitÃ© introduit l’Ã©chelle de
longueur nÃ©cessaire dans la formulation du matÃ©riau. De plus, nous prouvons
des bornes sur le domaine affectÃ© par cette contrainte. Une premiÃ¨re mise en
Åuvre par Ã©lÃ©ments finis unidimensionnels est proposÃ©e pour l’adoucissement
Ã©lastique ou plastique. adoucissement, localisation, endommagement,
plasticitÃ©, Lipschitz, Lip-field

##  1 Introduction

With softening, the stress that a material can sustain is diminishing as the
strain increases. This phenomenon exists both for elasticity and plasticity.
For elasticity, the stiffness is decreasing as the strain increases whereas
for plasticity, the yield stress is diminishing as the plastic strain
accumulates. Dealing with these types of models in finite element analysis is
a challenge. From the mathematical point of view, these models loose the nice
convex properties of classical elastic models or hardening plasticity models.
Existence and uniqueness of solutions is no longer ensured. Limit points may
exist: the solution does not exist beyond some loads. Bifurcation points are
also possible: at some stage in time, several solutions start to exist (stable
or not). Among the many solutions that may exist, some are called spurious
localizations. In the one-dimensional setting, these localizations are
characterized by a softening occuring in a single point. For higher
dimensions, the localization occurs on zero-measure domains. As a consequence,
the dissipation involved in these zones is zero. In other words, the external
energy or load needed to diminish the bearing capacity of the structure is
highly under-estimated. The mathematical difficulties of softening material
models have a direct impact on their numerical treatment. Spurious mesh
dependency is observed with finite elements: the mesh orientation has a strong
impact on the results and, as the mesh is refined, only a single element or a
layer of elements, depending on the problem dimension, is affected by
localization.

This paper introduces a new way to eliminate these spurious localizations from
the model. It does not mean that the solution is now unique but, at least,
these unwanted solutions are removed. The design of approaches to remove
spurious localizations has been going on for about fourty years now.

Regarding quasi-static analysis of time-independent models, several remedies
have been studied in the literature. They all share in common the fact that a
regularizing length scale is injected in the model. For the so-called non-
local integral damage model, the damage evolution at a given point is governed
by a driving force which is the average of the local driving force over some
distance around that point [1, 2, 3, 4]. In higher order, kinematically based,
gradient models, the length scale is introduced through the inclusion of
higher order deformation gradient in the energy expression [5, 6, 7] or
through additional rotational degrees of freedom [8]. For higher order, damage
based, gradient models, the energy depends on the gradient of the damage thus
involving again a length scale [9, 10, 11, 12].

Brittle fracture was revisited as an energy minimization problem in [13, 14]
giving the so-called variational approach to fracture [15]. At about the same
time, the phase-field approach was emanating from the physics community [16,
17] and then developped for mechanics applications [18, 19, 20, 21, 22].

Yet another way to introduce a length scale is the Thick Level Set approach to
fracture. The approach was introduced for brittle damage in [23, 24]. The
damage evolution is tied to a distance field. The damage front on which the
damage starts is the level set zero and it grows to a value of 1 (fully damage
state) at some distance $l_{\mathrm{c}}$ in the wake of the front. That
amounts to imposing the norm of the damage gradient on zones with strictly
positive damage (equality constraint). The possibility of an inequality
constraint on the damage gradient was also considered in [25, 26] allowing to
combine diffuse and localizing damage fields. The advantage of this inequality
was further stressed in [27], where it was shown that the inequality
constraint is convex on the contrary to the equality constraint. This paper
also demonstrates that a level set field is not mandatory: a variational
approach with Lagrange multipliers enforcing the inequality constraint may be
used. Following the ideas of [27], an implementation is provided in [28].

The Lipschitz regularization introduced in this paper enforces a regularity on
the damage field. The obtained field is called a Lip-field. The Lipschitz
regularity does not require the existence of the damage gradient. Yet, the
Lip-field is differentiable almost everywhere. To present the Lip-field
concept, the appealing framework of incremental energetic variational
potential is considered. The interest of this framework was first demonstrated
for visco-plasticity [29, 30] and was later reused as a building block for the
variational approach to fracture [15]. When softening may occur, the
incremental potential is not convex but only separately convex with respect to
the displacement (and non-softening internal variable) on one-side and
softening variable on the other side. It is thus natural to proceed with an
alternate minimization of successive convex problems [31]. The Lip-field
convex constraint is added as an extra constraint in the alternate
minimization. We demonstrate in this paper upper and lower bounds for the
minimization over the damage field. These bounds reduce drastically the zone
over which the Lipschitz condition needs to be activated. The proof holds in
any spatial dimension.

A one-dimensional finite element implementation is provided for the Lip-field
approach. The potential needs to be minimized for the nodal values of the
displacement and for the internal variables located at the integration point
in each finite element, including the damage state variable. In other words,
with the Lip-field approach, the damage variable may be kept at the
integration points with the other internal variables and does not need to be
stored at the nodes, thus following the common practice for nonlinear finite
element analysis.

Finally, note that there is also a current interest in Lipschitz
regularization to improve the robustness to adversarial perturbations for
learning framework [32].

The paper is organized as follows. The next section describes the classical
mechanical formulation for non-softening material models. In Section 3, a
softening variable is introduced and the Lipschitz constraint is imposed.
Elastic and plastic softening models are presented in Section 4 for the one-
dimensional setting. Finite element analysis are carried out in Section 5.
Discussion and future works are provided in Section 6.

##  2 Generalized standard materials

We consider the deformation of a body initially occupying a domain ${\Omega}$
through a displacement field $u$. For simplicity, we assume small, quasi-
static deformations. The Cauchy stress is denoted $\sigma$ and the strain
$\epsilon$

$\mbox{\boldmath$\epsilon$}(\mbox{\boldmath$u$})=\frac{1}{2}(\nabla\mbox{\boldmath$u$}+(\nabla\mbox{\boldmath$u$})^{\mathrm{T}})$
(1)

where $\nabla$ indicates the gradient operator. Regarding the boundary
conditions, the displacement is controlled on a part of the boundary denoted
${\Gamma}_{u}$ assumed fixed in time. On the rest of boundary, zero traction
force are assumed (again for simplicity). To be kinematically admissible at
some instant $t$, the displacement field must belong to $U(t)$:

$U(t)=\\{\mbox{\boldmath$u$}\in
H^{1}({\Omega}):\mbox{\boldmath$u$}=\mbox{\boldmath$u$}_{d}(t)\text{\ on\
}{\Gamma}_{u}\\}$ (2)

The equilibrium condition reads

$\int_{{\Omega}}\mbox{\boldmath$\sigma$}:\mbox{\boldmath$\epsilon$}(\mbox{\boldmath$u$}^{*})\mbox{\,
$\mathrm{d}$}{\Omega}=0,\quad\forall\mbox{\boldmath$u$}^{*}\in U^{*}$ (3)

where

$U^{*}=\\{\mbox{\boldmath$u$}\in H^{1}({\Omega}):\mbox{\boldmath$u$}=0\text{\
on\ }{\Gamma}_{u}\\}$ (4)

Kinematics and equilibrium equations (2-3) must be complemented with the
constitutive model. We consider the formalism of generalized standard material
introduced in [33, 34]. The set of internal variables is denoted $\alpha$. It
is a generic notation that describes a set of scalar, vectorial or tensorial
variables. The model is characterized by a free energy potential
$\varphi(\mbox{\boldmath$\epsilon$},\alpha)$ and a dissipation potential
$\psi(\dot{\alpha},\alpha)$. We then introduce an implicit time-discretization
and use the energetic variational approach. Consider the displacement and
internal variables $(\mbox{\boldmath$u$}_{n},\alpha_{n})$ known at some
instant $t_{n}$. Finding the pair $(\mbox{\boldmath$u$}_{n+1},\alpha_{n+1})$
at the next instant $t_{n+1}=t_{n}+\Delta t$ amounts to a minimization problem

$(\mbox{\boldmath$u$}_{n+1},\alpha_{n+1})=\arg\min_{\begin{subarray}{c}\mbox{\boldmath$u$}^{\prime}\in
U_{n}\\\ \alpha^{\prime}\in
A_{n}\end{subarray}}F(\mbox{\boldmath$u$}^{\prime},\alpha^{\prime};\mbox{\boldmath$u$}_{n},\alpha_{n},\Delta
t)$ (5)

where $U_{n}$ is a short-hand notation for $U(t_{n+1})$ and $A_{n}$ indicates
the restriction on $\alpha_{n+1}$. The incremental potential $F$ involves the
energy and dissipation potentials $\varphi$ and $\psi$. Examples of
incremental potentials will be given in Section 4.

For simplicity, we shall consider time-independent material models. In this
case, the $F$ expression does not depend explicitely on
$\mbox{\boldmath$u$}_{n},\alpha_{n}$ and $\Delta t$. The extension to time-
dependent models does not introduce difficulties. Also, to lighten notation,
we drop the $n+1$ indices. The minimum problem is then

$(\mbox{\boldmath$u$},\alpha)=\arg\min_{\begin{subarray}{c}\mbox{\boldmath$u$}^{\prime}\in
U_{n}\\\ \alpha^{\prime}\in
A_{n}\end{subarray}}F(\mbox{\boldmath$u$}^{\prime},\alpha^{\prime})$ (6)

We assume that the domains $U_{n}$ and $A_{n}$ are convex and that $F$ is
strictly convex with respect to the pair $(\mbox{\boldmath$u$},\alpha)$ over
$U_{n}\times A_{n}$. These conditions ensures that the minimization yields an
existing unique solution at each time-step. This defines a non-softening
model.

The minimization above is traditionally solved by a repeated sequence of two
steps: the computation of the internal variables (and stress) for a given
displacement field, followed by the correction of the displacement field. At
iteration $m$, the two steps are given below

$\displaystyle\alpha^{m+1}$ $\displaystyle=\arg\min_{\alpha\in
A_{n}}F(\mbox{\boldmath$u$}^{m},\alpha)$ (7)
$\displaystyle\\{\mbox{\boldmath$u$}^{m+1}\\}$
$\displaystyle=\\{\mbox{\boldmath$u$}^{m}\\}+(K^{m})^{-1}R^{m},\quad\mbox{\boldmath$u$}^{m+1}\in
U_{n}$ (8)

The first step, (7), is purely local and may be carried out independently at
each material point. The second step, (8), involves a linear solve which
updates the degrees of freedom {$u$} associated to the field $u$. The matrix
$K$ depends on the current internal variables $\alpha^{m+1}$ and the residual
vector $R$ depends on the current stress $\mbox{\boldmath$\sigma$}^{m+1}$. The
matrix can be the algorithmic tangent operator [35] or some approximation of
it. The linear solver may possibly be followed by a line search to further
improve the solution. Expressions of $K$ and $R$ are given in the Appendix for
the models at stake in this paper.

##  3 Softening variable and Lipschitz regularization

We now consider that the model has an extra scalar variable $d$ responsible
for softening. The optimization problem becomes

$(\mbox{\boldmath$u$},\alpha,d)=\arg\min_{\begin{subarray}{c}\mbox{\boldmath$u$}^{\prime}\in
U_{n}\\\ \alpha^{\prime}\in A_{n}\\\ d^{\prime}\in
D_{n}\end{subarray}}F(\mbox{\boldmath$u$},\alpha^{\prime},d^{\prime})$ (9)

It is no longer convex, we can expect several local minima and also non-
uniqueness to the global minimum (i.e. several solutions leading to the same
global minimum). Even though $F$ is not convex with respect to the triple
$(\mbox{\boldmath$u$},\alpha,d)$, we ask the optimization to be convex with
respect to the couple $(\mbox{\boldmath$u$},\alpha)$ for all $d\in D_{n}$ and
to be convex with respect to $d$ for all $(\mbox{\boldmath$u$},\alpha)\in
U_{n}\times A_{n}$.

Given a $d$ field over the domain ${\Omega}$, the Lipschitz constant
associated to this field is the minimum $M$ value such that the following
holds

$\lvert d(\mbox{\boldmath$x$})-d(\mbox{\boldmath$y$})\rvert\leq
M\,\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$}),\quad\forall\mbox{\boldmath$x$},\mbox{\boldmath$y$}\in{\Omega}$
(10)

where $\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$})$ is the minimal
length of the path inside ${\Omega}$ joining $x$ and $y$ (the distance is
considered infinite if the two points cannot be connected inside ${\Omega}$).
The $\mathrm{dist}$ function is a metric since it satisfies for all
$\mbox{\boldmath$x$},\mbox{\boldmath$y$},\mbox{\boldmath$z$}\in{\Omega}$:

$\displaystyle\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$})=0\Leftrightarrow\mbox{\boldmath$x$}=\mbox{\boldmath$y$}$
(11)
$\displaystyle\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$})=\mathrm{dist}(\mbox{\boldmath$y$},\mbox{\boldmath$x$})$
(12)
$\displaystyle\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$})\leq\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$z$})+\mathrm{dist}(\mbox{\boldmath$z$},\mbox{\boldmath$y$})$
(13)

The value $M$ defined above is denoted $\mathrm{lip}(d)$. We define the
regularization space for the damage field

$\mathrm{L}=\\{d\in L^{\infty}({\Omega}):\mathrm{lip}(d)\leq 1/l\\}$ (14)

where $l$ is the regularizing length. The set $L$ is convex. We seek the
solution as "one" of the global minima

$(\mbox{\boldmath$u$},\alpha,d)=\arg\min_{\begin{subarray}{c}(\mbox{\boldmath$u$},\alpha)\in
U_{n}\times A_{n}\\\ d\in
D_{n}\cap\mathrm{L}\end{subarray}}F(\mbox{\boldmath$u$},\alpha,d)$ (15)

We proceed by alternate minimization.

$\displaystyle(\mbox{\boldmath$u$}^{k+1},\alpha^{k+1})$
$\displaystyle=\arg\min_{(\mbox{\boldmath$u$},\alpha)\in U_{n}\times
A_{n}}F(\mbox{\boldmath$u$},\alpha,d^{k})$ (16) $\displaystyle d^{k+1}$
$\displaystyle=\arg\min_{d\in
D_{n}\cap\mathrm{L}}F(\mbox{\boldmath$u$}^{k+1},\alpha^{k+1},d)$ (17)

For the first minimization, the damage variable is frozen and the problem is
thus identical to problem (6). It is a classical non-softening step. The
second minimization, (17), is less common. The objective function is convex
(and separable) as well as the constraint $D_{n}$. The Lip constraint is non-
local as it ties the damage variables between points. The optimization to find
$d$ is thus potentially time-consuming when turning to a numerical
implementation. The good news is that the quest for $d$ may be decomposed into
three steps reducing dramatically the cost of the optimization. The first step
is to create a trial $d$ field denoted $\overline{d}$ by ignoring the
Lipschitz constraint and performing a decoupled minimization at each point

$\overline{d}=\arg\min_{d\in
D_{n}}F(\mbox{\boldmath$u$}^{k+1},\alpha^{k+1},d)$ (18)

If the trial damage $\overline{d}$ satisfies the Lipschitz constraint, it is
the solution we are looking for as indicated in Figure 1 (left). If not, the
optimal damage field will be different from $\overline{d}$.

$d_{n}$
---
$d_{n}$
---
$\bar{d}=d^{k+1}$
---
$\bar{d}$
---
$L$
---
$L$
---
$\pi^{u}\bar{d}$
---
$d^{k+1}$
---
$\pi^{l}\bar{d}$
---
Figure 1: A sketch of the local update $\overline{d}$ from the previous known
damage field $d_{n}$. If the update satifies the Lip constaint (left), we have
directly the solution $d^{k+1}$. Othewise (right), the local update needs to
be projected back to the Lip constraint while minimizing the objective
function $F$. The upper and lower projection give bounds to the solution.

We can find an upper bound of the domain over which $d^{k+1}$ will differ from
$\overline{d}$. We define two projections onto $L$, a lower projection
$\pi^{\mathrm{l}}$ and an upper projection $\pi^{\mathrm{u}}$:

$\displaystyle\pi^{\mathrm{l}}\overline{d}(\mbox{\boldmath$x$})$
$\displaystyle=\min_{\mbox{\boldmath$y$}\in{\Omega}}\;(\overline{d}(\mbox{\boldmath$y$})+\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$}))$
(19) $\displaystyle\pi^{\mathrm{u}}\overline{d}(\mbox{\boldmath$x$})$
$\displaystyle=\max_{\mbox{\boldmath$y$}\in{\Omega}}\;(\overline{d}(\mbox{\boldmath$y$})-\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$}))$
(20)

We prove in the appendix that these projections satisfy the following
inequality:

$d_{n}\leq\pi^{\mathrm{l}}\overline{d}\leq\overline{d}\leq\pi^{\mathrm{u}}\overline{d}\leq
1$ (21)

and provide bounds for the optimal damage

$\pi^{\mathrm{l}}\overline{d}\leq d^{k+1}\leq\pi^{\mathrm{u}}\overline{d}$
(22)

As an important consequence, the trial and optimal solutions coincide wherever
the bounds are equal

$\pi^{\mathrm{l}}\overline{d}(\mbox{\boldmath$x$})=\pi^{\mathrm{u}}\overline{d}(\mbox{\boldmath$x$})\quad\Rightarrow\quad
d^{k+1}(\mbox{\boldmath$x$})=\overline{d}(\mbox{\boldmath$x$})$ (23)

A sketch of the projections is given in Figure 1 (right) and they are also
illustrated on the one-dimensional example, Figure 2.

Figure 2: A sketch of the bracketing capability of the upper and lower
projections. Outside the interval $[a,b]$, the projections are equal and thus
$d^{k+1}=\overline{d}$. Inside the interval $[a,b]$, the projections bracket
$\overline{d}$ and $d^{k+1}$.

##  4 Elastic and plastic softening models

Consider the one-dimensional model of a bar of length $L$ and unit section
attached at its left side and pulled with an imposed displacement
$u_{\mathrm{d}}(t)$ at it’s right end. The displacement at time $t$ must
belong to the set

$U(t)=\\{u\in C([0,L]):u(0)=0,\;u(L)=u_{\mathrm{d}}(t)\\}$

whereas the set of admissible displacement variations is given by

$U^{*}=\\{u\in C([0,L]):u(0)=0,\;u(L)=0\\}$

We now detail several incremental potentials to be used in the simulation.

###  4.1 Softening elasticity

The elastic potential reads

$F_{\mathrm{e}}(u)=\int_{0}^{L}f_{e}({\epsilon}(u))\mbox{\,
$\mathrm{d}$}x,\quad
f_{\mathrm{e}}({\epsilon}(u))=\frac{1}{2}E\epsilon(u)^{2}$ (24)

where $\epsilon(u)=\frac{\mathrm{d}u}{\mathrm{d}x}$ and $E$ is the Young
modulus. The elastic softening model, affects $E$ and adds a dissipation term.
For simplicity, we do not consider tension compression dissymmetry in the
model (the bar will always be in tension in the simulation). The softening
elasticity potential is

$F_{\mathrm{se}}(u,d)=\int_{0}^{L}f_{\mathrm{se}}({\epsilon}(u),d)\mbox{\,
$\mathrm{d}$}x,\quad
f_{\mathrm{se}}({\epsilon}(u),d)=(1-d)^{2}f_{\mathrm{e}}({\epsilon}(u))+Y_{\mathrm{c}}h(d)$
(25)

where $Y_{\mathrm{c}}$ is the critical energy release rate and the convex
function $h(d)$ defines the softening behavior. Damage can only grow and
cannot go beyond 1. The convex constraint set $D_{n}$ for damage is:

$D_{n}=\\{d\in L^{\infty}([0,L]):d_{n}\leq d\leq 1\\}$ (26)

The power 2 over the factor $(1-d)$ ensures convexity and a finite opening as
$d\rightarrow 1$ [36]. Note that the proper choice of the power is also an
issue with gradient damage models for which a value of 2 also ensures a finite
opening [37]. We consider two choices for $h(d)$

$\displaystyle h_{1}(d)$ $\displaystyle=2d+3d^{2}$ (27) $\displaystyle
h_{2}(d)$ $\displaystyle=\frac{2d-d^{2}}{(1-d+{\lambda}d^{2})^{2}}$ (28)

Both functions are convex (provided ${\lambda}\leq 1/2$). The condition
$h^{\prime}(0)=2$ ensures that when damage starts
$E\epsilon^{2}/2=Y_{\mathrm{c}}$. The second choice allows to mimic a linear
cohesive-zone model [38] where ${\lambda}$ is defined by

${\lambda}=2\frac{Y_{\mathrm{c}}l}{G_{\mathrm{c}}}$ (29)

with $G_{\mathrm{c}}$ denoting the toughness (energy per section area needed
to break the bar). We shall see that the choice of $h_{2}(d)$ allows to choose
independently $Y_{\mathrm{c}}$, $G_{\mathrm{c}}$ and $l$.

To get a better insight on the potential, we write the Karush-Kuhn-Tucker
(KKT) conditions associated to the optimization. We define the dual quantities
to the strain and damage by taking the derivative of $f_{\mathrm{se}}$. We get
the stress ${\sigma}$ and a variable associated to damage denoted $\mu$
(damage criterion):

$\displaystyle\sigma$ $\displaystyle=(1-d)^{2}E\epsilon$ (30)
$\displaystyle\mu$
$\displaystyle=-(1-d)E\epsilon^{2}+Y_{\mathrm{c}}h^{\prime}(d)$ (31)

where $h^{\prime}(d)$ denotes the derivative of $h$ with respect to $d$. The
KKT conditions read

$\displaystyle\mu-{\lambda}_{1}+{\lambda}_{2}=0$ (32)
$\displaystyle{\lambda}_{1}\geq 0,\quad d-d_{n}\geq
0,\quad{\lambda}_{1}(d-d_{n})=0$ (33) $\displaystyle{\lambda}_{2}\geq 0,\quad
1-d\geq 0,\quad{\lambda}_{2}(1-d)=0$ (34)

where ${\lambda}_{1}$ and ${\lambda}_{2}$ are Lagrange multipliers associated
to the constraints. We note that $\mu=0$ when damage is growing (and stays
below 1), and $\mu\geq 0$ when damage does not grow (and is different from 1).
Figure 3 shows the stress-strain relation for a strain loading-unloading
history.

Figure 3: Stress-strain relations for the softening elastic model
(${\sigma}_{\mathrm{c}}=\sqrt{2EY_{\mathrm{c}}},\epsilon_{\mathrm{c}}={\sigma}_{\mathrm{c}}/E$).
The strain loading/unloading follows the peak values sequence:
${\epsilon}/\epsilon_{\mathrm{c}}=(0,2,0,4)$.

The minimization of $F_{se}$ under the constraint $(u,d)\in U_{n}\times D_{n}$
is not a convex problem. But, the minimization of $F_{se}$ with respect to
$u\in U_{n}$, for $d\in D_{n}$ fixed is a convex problem. The minimization of
$F_{se}$ with respect to $d\in D_{n}$ for $u\in U_{n}$ fixed is also a convex
problem. Finally, note that in the presentation above, we did not care in
presenting separately the free energy and dissipation potentials. This is not
essential since we are not studying the temperature evolution. Moreover, for a
given stress-strain relation this choice is in general non-unique.

###  4.2 Softening elasticity with hardening plasticity

The introduction of softening in elasto-plastic models is a complex topic. The
goal here is not to find the appropriate model for a given situation but
rather to discuss how the Lip-field approach is behaving in the presence of
plastic internal variables. We consider a basic von Mises plasticity model
with elastic softening and, in the next section, a softening plasticity model
with preserved elasticity. The incremental potential for a von Mises isotropic
hardening plasticity model is

$F_{\mathrm{p}}(u,\epsilon_{p},p)=\int_{0}^{L}f_{\mathrm{p}}({\epsilon}(u),\epsilon_{p},p)\mbox{\,
$\mathrm{d}$}x,\quad
f_{\mathrm{p}}({\epsilon}(u),\epsilon_{p},p)=\frac{1}{2}E(\epsilon(u)-\epsilon_{p})^{2}+\sigma_{\mathrm{y}}(p+\frac{k}{2}p^{2})$
(35)

where $\sigma_{\mathrm{y}}$ is the yield stress and $k$ the isotropic
hardening parameter. Regarding the constraints, plasticity internal variables
must belong to the following convex set

$A_{n}=\\{(\epsilon_{p},p)\in(L^{\infty}([0,L]))^{2}:p-p_{n}\geq\lvert\epsilon_{p}-\epsilon_{pn}\rvert\\}$
(36)

The minimization of $F_{\mathrm{sp}}$ with respect to $(u,(\epsilon_{p},p))\in
U_{n}\times A_{n}$ is a convex problem. We introduce damage with a
multiplicative approach, also called effective stress approach:

$\displaystyle F_{\mathrm{sep}}(u,\epsilon_{p},p,d)$
$\displaystyle=\int_{0}^{L}f_{\mathrm{sep}}({\epsilon}(u),\epsilon_{p},p,d)\mbox{\,
$\mathrm{d}$}x,$ (37) $\displaystyle
f_{\mathrm{sep}}({\epsilon}(u),\epsilon_{p},p,d)$
$\displaystyle=(1-d)^{2}f_{\mathrm{p}}({\epsilon}(u),\epsilon_{p},p)+Y_{\mathrm{c}}h(d)$
(38)

The damage must belong to the set $D_{n}$ defined in (26).

Figure 4: Stress-strain relations for the softening elasticity model with
hardening plasticity with the choice
$\epsilon_{\mathrm{c}}/\epsilon_{\mathrm{y}}=2$ where
$\epsilon_{\mathrm{y}}=\sigma_{\mathrm{y}}/E$. The strain loading/unloading
follows the peak value sequence
${\epsilon}/\epsilon_{\mathrm{y}}=(0,2,e_{1},4,e_{2},5,e_{3})$ where $e_{1}$,
$e_{2}$, $e_{3}$ indicate the strain needed to reach full unloading
(${\sigma}=0$).

Taking the derivative of $f_{\mathrm{sep}}$ with respect to the plastic
strain, cumulative plasticity and damage, gives the stress, current yield
stress and the damage criterion, respectively

$\displaystyle\sigma$ $\displaystyle=(1-d)^{2}E(\epsilon-\epsilon_{p})$ (39)
$\displaystyle R$ $\displaystyle=\sigma_{\mathrm{y}}(1-d)^{2}(1+kp)$ (40)
$\displaystyle\mu$
$\displaystyle=-(1-d)E(\epsilon(u)-\epsilon_{p})^{2}-2(1-d)\sigma_{\mathrm{y}}(p+kp^{2}/2)+Y_{\mathrm{c}}h^{\prime}(d)$
(41)

They are involved in the KKT conditions

$\displaystyle{\sigma}+{\lambda}_{\mathrm{p}}s(\epsilon_{\mathrm{p}}-\epsilon_{\mathrm{p}n})=0$
(42) $\displaystyle R-{\lambda}_{\mathrm{p}}=0$ (43)
$\displaystyle{\lambda}_{\mathrm{p}}\geq 0,\quad
p-p_{n}-\lvert\epsilon_{\mathrm{p}}-\epsilon_{\mathrm{p}n}\rvert\geq
0,\quad{\lambda}_{\mathrm{p}}(p-p_{n}-\lvert\epsilon_{\mathrm{p}}-\epsilon_{\mathrm{p}n}\rvert)=0$
(44)

where $s$ is the multi-valued signed function ($s(x)=-1$ if $x<0$, $s(x)=+1$
if $x>0$ and $s(0)\in[-1,1]$). We observe that the effective stress
${\sigma}/(1-d)^{2}$ and variables $\epsilon_{\mathrm{p}}$ and $p$ may be
computed from the strain independently of the damage variable. Stress-strain
relations for a strain loading-unloading history are shown in Figure 4. The
effect of damage may be observed in the unloading phase.

###  4.3 Softening plasticity

Finally, we consider an elasto-plastic model in which the softening affects
only the yield stress leaving elasticity unchanged:

$\displaystyle
F_{\mathrm{sp}}(u,\epsilon_{p},p,d)=\int_{0}^{L}f_{\mathrm{sp}}({\epsilon}(u),\epsilon_{p},p,d)\mbox{\,
$\mathrm{d}$}x,$ $\displaystyle
f_{\mathrm{sp}}({\epsilon}(u),\epsilon_{p},p,d)=\frac{1}{2}E(\epsilon(u)-\epsilon_{p})^{2}+(1-d)^{2}\sigma_{\mathrm{y}}(p+kp^{2}/2)+\sigma_{\mathrm{y}}g(d)$

where $g(d)$ describes softening and is chosen as $g(d)=d^{2}$.

Figure 5: Stress-strain relations for the softening plastic model
($\epsilon_{\mathrm{y}}=\sigma_{\mathrm{y}}/E)$. The strain loading/unloading
follows the peak value sequence ${\epsilon}/\epsilon_{\mathrm{y}}=(0,1.5,e,5)$
where $e$ indicates the value needed to reach ${\sigma}=0$.

Taking the derivative of $f_{\mathrm{sp}}$, we get

$\displaystyle\sigma$ $\displaystyle=E(\epsilon-\epsilon_{p})$ (45)
$\displaystyle R$ $\displaystyle=\sigma_{\mathrm{y}}(1-d)^{2}(1+kp)$ (46)
$\displaystyle\mu$
$\displaystyle=-2\sigma_{\mathrm{y}}(1-d)(p+kp^{2}/2)+\sigma_{\mathrm{y}}g^{\prime}(d)$
(47)

and the associated KKT conditions

$\displaystyle{\sigma}+{\lambda}_{\mathrm{p}}s(\epsilon_{\mathrm{p}}-\epsilon_{\mathrm{p}n})=0$
(48) $\displaystyle R-{\lambda}_{\mathrm{p}}=0$ (49)
$\displaystyle{\lambda}_{\mathrm{p}}\geq 0,\quad
p-p_{n}-\lvert\epsilon_{\mathrm{p}}-\epsilon_{\mathrm{p}n}\rvert\geq
0,\quad{\lambda}_{\mathrm{p}}(p-p_{n}-\lvert\epsilon_{\mathrm{p}}-\epsilon_{\mathrm{p}n}\rvert)=0$
(50)

We have omitted above the equations associated to $d$. They are identical to
the ones presented in (32-34). The $d$ variable may be interpreted as a void
volume fraction quantity used in ductile failure models [39, 40]. Its growth
is linked to the cumulative plastic strain and no longer to the elastic strain
(as indicated by the difference between (41) and (47)). Stress-strain
relations for a strain loading-unloading history are shown in Figure 5. The
elasticity is not affected as can be observed from the unloading phase.

$\epsilon_{p},p,d$
---
$\epsilon_{p},p,d$
---
$\epsilon_{p},p,d$
---
$\epsilon_{p},p,d$
---
$u=0$
---
$u$
---
$u$
---
$u$
---
$u$
---
$u=u_{d}(t)$
---
Figure 6: A one-dimensionbal bar discretized with equal-sized finite elements.
Displacements are stored at the node whereas internal variables are stored at
the element integration point (depicted with a vertical bar).

##  5 Lipschtiz regularization over a bar

The bar is discretized with $N$ finite elements of equal size $h=L/N$. The
displacement is linear over each element between nodal values. The internal
variables are stored at each element centroid (element integration point) as
indicated in Figure 6. The $i$ index is used either for the element numbering
( $d_{i},i=1,\ldots,N$ denote the damage at the centroid of each element) or
for the node numbering ($u_{i},i=1,\ldots,N+1$ are the nodal displacements).
The Lipschitz constraint implies the following inequalities defining the set
$L$

$\displaystyle d_{i}-d_{i+1}-h/l$ $\displaystyle\leq 0,\quad i=1,\ldots,N-1$
$\displaystyle d_{i}-d_{i-1}-h/l$ $\displaystyle\leq 0,\quad i=2,\ldots,N$

The optimization problem reads for the softening elasticity problem

$\min_{\begin{subarray}{c}\mbox{\boldmath$u$}\in U_{n}\\\ d\in D_{n}\cap
L\end{subarray}}\sum_{i=1}^{N}hf_{\mathrm{se}}(\frac{u_{i+1}-u_{i}}{h},d_{i})$
(51)

and

$\min_{\begin{subarray}{c}(\mbox{\boldmath$u$},(\epsilon_{\mathrm{p}},p))\in
U_{n}\times A_{n}\\\ d\in D_{n}\cap
L\end{subarray}}\sum_{i=1}^{N}hf_{\mathrm{x}}(\frac{u_{i+1}-u_{i}}{h},\epsilon_{\mathrm{p}i},p_{i},d_{i})$
(52)

for the plasticity models ($\mathrm{x}$ stands for $\mathrm{sep}$ or
$\mathrm{sp}$). The set $U_{n}$ enforces that $u_{0}=0$ and
$u_{N+1}=u_{\mathrm{d}}(t_{n+1})$. The set $A_{n}$ is related to (36) and the
set $D_{n}$ enforces the fact that $d_{i}$ must be above its previous time-
step value and below 1. For all examples treated, the softening elasticity
function is chosen as $h_{2}$.

|
---|---
|
Figure 7: Response for the bar with the elastic softening behavior and imposed
displacement: non-regularized model (left) and Lip-field model (right). The
top row gives the stress-strain average response and the bottom one the damage
profile at the color dots given in the middle row. The mesh used for the
middle row is the most refined of the top row. Parameters are: $L$=1, $l$ =
0.5, $E$ = 1, $Y_{\mathrm{c}}$ = 1, ${\lambda}$ = 0.3, $G_{\mathrm{c}}$ =
10/3.

The optimization problems above are solved by alternate minimization as
explained in Section 3. The minimization with respect to the damage variable
is performed using the scipy python library. More precisely, the
scipy.optimize.minimize function is used. It is described in [41]. It offers
an interface to the SLSQP (Sequential Least Square Quadratic Programming)
routine created by Dieter Kraft in [42].

Because of its homogeneity, the optimation problem at stake has no unique
global minimizer and many local minimizers. Aware of this situation, we
trigger localization by introducing a slightly higher damage on top of $d_{n}$
in the middle of the bar at each step. This is just an initial guess for the
alternate minimization. This strategy does not require to alter the stiffness
or initial damage on the bar.

For the elastic softening model, Figure 7 gives the bar response for the non-
regularized and Lip-field models. The regularizing effect of the latter is
clear. Note that the case of two nodes for the discretization (a single
element) gives the homogeneous bar response. No snap-back appears in this
example because the regularization length is pretty large. A smaller length
and toughness are used for Figure 8 and a snap-back may be observed. To handle
the snap-back in the simulation, the displacement at the end of the bar is no
longer imposed but controlled by limiting the strain increment over the bar at
each time step. More evolved control could be used [43, 44]. Figure 8 shows
for the non-regularized case (left column) a convergence towards a solution
without any dissipation (zero area under the converged curve). On the
contrary, the converged solution, (right column) indicates a non-zero
dissipation (whose dimensional value is close to the toughness
$G_{\mathrm{c}}$).

|
---|---
|
Figure 8: Response for the bar with the elastic softening behavior and snap-
back control: non-regularized model (left) and Lip-field model (right),
convergence of the stress-average strain response (top) and damage evolution
(bottom) corresponding the color dots in the middle figures. The middle
figures reproduce the curve of the top figures for the most refined grid.
Parameters are: $L$=1, $l$ = 0.1, $E$ = 1, $Y_{\mathrm{c}}$ = 1, ${\lambda}$ =
0.3, $G_{\mathrm{c}}$ = 2/3.

Still for the elastic softening model, we show, Figure 9, that the stress-
strain response does not vary much with the regularization length
(demonstrating the nice capability of the $h_{2}$ softening function). The
response are obtained for a non-varying toughness, $G_{\mathrm{c}}$, and
strength $Y_{\mathrm{c}}$. Only the regularization length is varied (and the
mesh size is adjusted to always have the same number of elements over $l$). We
observe a dimensional area under the curve close to $G_{\mathrm{c}}$.

Figure 9: Lip-field model response with the elastic softening behavior (choice
$h_{2}$) for several regularization length values. Stress-average strain
response (above) and damage profile (below). Parameters are: $L$ = 1., $l$ =
(0.4, 0.2, 0.1), $E$ = 1., $Y_{\mathrm{c}}$ = 1., ${\lambda}$ = (0.4,0.2,0.1),
$G_{\mathrm{c}}$ = 2. The number of elements is 51, 101 and 201 so that we
have $h/l\sim 20$.

Regarding the softening elastic model with hardening plasticity, the response
is given in Figure 10 (left column). We observe that plasticity proceeds first
in an homogeneous fashion, without any damage. During this homogeneous phase,
the stress-average strain curve is rising. As this curve reaches its peak,
localization of the damage and plasticity start to occur. For the softening
plasticity model, the response is given in Figure 10 (right column). We
observe that plasticity and damage proceed first in an homogeneous manner due
to the initial hardening effect of plasticity, then concentrates. The
departure from the homogeneous response also occurs when the stress-average
strain reaches its limit point.

|
---|---
Figure 10: Lip-field model response for the softening elasticity model with
plasticity hardening (left column) and softening plasticity (right column):
stress-strain average curves (top row), damage field (middle row) and
cumulative plasticity (bottom row). For both models, $L=1$, $l=0.5$ and the
mesh size is $h=L/64$. Model specific parameters are: $E=2$,
$\sigma_{\mathrm{y}}$ = 1, $Y_{\mathrm{c}}$ = 1, k=1, ${\lambda}$ = 1/3 for
the hardening plasticity (corresponding to the model the orange curve in
Figure 4) and E = 1., $\sigma_{\mathrm{y}}$ = 1/16, k = 4 for the softening
plasticity model (corresponding to the model with the dashed green curve in
Figure 5).

We provide a last example with a non-uniform loading caused by a volumic force
$f(x/L)=0.1*\sin(8*\pi*x/L)$. This loading is added to the imposed
displacement. The model parameters are the same as for Figure 10 (right
column) except that $l=0.25$ and the number of nodes is 256. The damage
develops in a complex pattern. The right-column gives the results for the non-
regularized case (local approach $l=0$). We observe the Lip-field results
(right column) are the same as the local approach as long as the stress-
average strain curve did not reach the peak point. This demonstrates the fact
that the Lip-field approach presevers the local solution when it is not
localizing. This is quite different from the other regularizing approaches
(non-local, damage gradient, phase-field) that will modify the local solution
even prior to localization. It is also clear from the damage profile that it
does not satisfy any specific boundary conditions.

|
---|---
Figure 11: Response for a softening plasticity model with a volumic
oscillatory force. No regularization, $l=0$ (left) and Lip regularization
(right). Stress-average strain results (top), damage (middle) and cumulative
plasticity profile (bottom). Parameters are: $L$ = 1., $l$ = 0.25, $h$ =
1/255, $E$ = 1., $\sigma_{\mathrm{y}}$ = 1/16, $k$ = 4.

##  6 Conclusion and future works

A new regularization approach has been introduced to alleviate spurious
localization with softening material models. Both elastic and plastic
softening models have been considered. The regularization enforces a Lipschitz
condition on the field responsible for the softening. In doing so, a length is
introduced in the model.

Compared to the gradient damage or phase-field approaches, the Lip-field
approach introduces an extra constraint but does not affect the expression of
the incremental potential (the objective function in the minimization is not
modified). In other words, the energy does not depend on the damage gradient
but solely on the damage. Its expression is the one of the local model. As a
consequence, the Lip-field does not introduce an extra partial differential
equation with its questionable boundary conditions.

Compared to the non-local integral approach, no widening of the fully damage
zone is observed and the implementation only requires an element to be related
to its neighbors. Finally, compared to the Thick Level Set approach, the level
set know-how is no longer needed.

Future works will be dedicated to the extension of the results to two- and
three-dimensional problems. For this extension, the bounds demonstrated in
this paper (valid for any dimensions) will allow one to predict a priori the
zones over which damage is affected by the Lipschitz constraint. This will
reduce the cost of the damage optimization.

Finally, the paper did not discuss the possibility of allowing displacement
jumps during the simulation. This issue is important when dealing with
fragmentation or crack under large deformation and needs to be studied in the
future for the Lip-field approach.

## Acknowledgements

The authors wish to thank J. Dolbow, B. Le, B. Masseron, G. Rastiello, J.
Rethore, A. Salzman, I. Stefanou, A. Stershic and C. Stolz, for their advice
on ways to improve the manuscript.

##  Appendix A Proof of the bounds

Consider a field $d$ defined over ${\Omega}$. To this field, we associate two
other fields denoted, $\pi^{\mathrm{u}}d$ and $\pi^{\mathrm{l}}d$, called
upper and lower projections, respectively.

$\displaystyle\pi^{\mathrm{u}}:L^{\infty}({\Omega})\mapsto\mathrm{L}:d(\mbox{\boldmath$x$})\mapsto\pi^{\mathrm{u}}d(\mbox{\boldmath$x$})=\max_{\mbox{\boldmath$y$}\in{\Omega}}\;(d(\mbox{\boldmath$y$})-\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$})),$
(53)
$\displaystyle\pi^{\mathrm{l}}:L^{\infty}({\Omega})\mapsto\mathrm{L}:d(\mbox{\boldmath$x$})\mapsto\pi^{\mathrm{l}}d(\mbox{\boldmath$x$})=\min_{\mbox{\boldmath$y$}\in{\Omega}}\;(d(\mbox{\boldmath$y$})+\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$}))$
(54)

They satisfy the following properties:

  1. (a)

$\pi^{\mathrm{l}}d\in\mathrm{L}$, $\pi^{\mathrm{u}}d\in\mathrm{L}$

  2. (b)

$d\in L\Rightarrow\pi^{\mathrm{l}}d=\pi^{\mathrm{u}}d=d$

  3. (c)

$\pi^{\mathrm{l}}d\leq d\leq\pi^{\mathrm{u}}d$

  4. (d)

$d_{1}\leq
d_{2}\Rightarrow\pi^{\mathrm{l}}d_{1}\leq\pi^{\mathrm{l}}d_{2},\;\pi^{\mathrm{u}}d_{1}\leq\pi^{\mathrm{u}}d_{2}$
(monotonicity)

To prove that the mapping falls indeed in $\mathrm{L}$ (property (a)), we take
the difference

$\displaystyle\pi^{\mathrm{u}}d(\mbox{\boldmath$x$}^{1})-\pi^{\mathrm{u}}d(\mbox{\boldmath$x$}^{2})=\max_{\mbox{\boldmath$y$}\in{\Omega}}\;(d(\mbox{\boldmath$y$})-\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$}^{1},\mbox{\boldmath$y$}))-\max_{\mbox{\boldmath$y$}\in{\Omega}}\;(d(\mbox{\boldmath$y$})-\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$}^{2},\mbox{\boldmath$y$}))$
(55)

Using the triangular inequality written as

$\mathrm{dist}(\mbox{\boldmath$x$}^{1},\mbox{\boldmath$y$})\leq\mathrm{dist}(\mbox{\boldmath$x$}^{1},\mbox{\boldmath$x$}^{2})+\mathrm{dist}(\mbox{\boldmath$x$}^{2},\mbox{\boldmath$y$})$

we get

$\pi^{\mathrm{u}}d(\mbox{\boldmath$x$}^{1})-\pi^{\mathrm{u}}d(\mbox{\boldmath$x$}^{2})\geq-\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$}^{1},\mbox{\boldmath$x$}^{2})$

and using the triangular inequality as

$\mathrm{dist}(\mbox{\boldmath$x$}^{2},\mbox{\boldmath$y$})\leq\mathrm{dist}(\mbox{\boldmath$x$}^{2},\mbox{\boldmath$x$}^{1})+\mathrm{dist}(\mbox{\boldmath$x$}^{1},\mbox{\boldmath$y$})$

we get

$\pi^{\mathrm{u}}d(\mbox{\boldmath$x$}^{1})-\pi^{\mathrm{u}}d(\mbox{\boldmath$x$}^{2})\leq\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$}^{1},\mbox{\boldmath$x$}^{2})$

The proof is similar for $\pi^{\mathrm{l}}d$. Property (b) is obtained from
the definition of $d\in\mathrm{L}$:

$d(\mbox{\boldmath$y$})-\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$})\leq
d(\mbox{\boldmath$x$})\leq
d(\mbox{\boldmath$y$})+\frac{1}{l}\mathrm{dist}(\mbox{\boldmath$x$},\mbox{\boldmath$y$})$
(56)

Computing the minimum and maximum over $y$ for the upper and lower bounds,
respectively yields the property. Property (c) is obtained directly from the
definition of the projections by testing $y$ as $x$ in the max and min.
Similarly, property (d) stems from the definition of the projections.

We now discuss the consequences of the properties. Properties (a) and (b)
indicate that the projections are idempotent (applied twice yield the same
result as applied once). Property (c) indicates that at any point where the
projections are equal, they are equal to $d$. Combining properties (c) and (d)
and using the fact that $d_{n}\in\mathrm{L}$ (as well as noting that any
uniform function is Lipschitz) yield the result (21).

We are now ready to prove the bounds result (22). Consider a damage field
$d^{*}\in D_{n}\cap\mathrm{L}$ lying ouside the bounds, we associate to this
field another field by clipping it to the bounds

$d^{**}(\mbox{\boldmath$x$})=\max(\pi^{\mathrm{l}}\overline{d}(\mbox{\boldmath$x$}),\min(d^{*}(\mbox{\boldmath$x$}),\pi^{\mathrm{u}}\overline{d}(\mbox{\boldmath$x$}))$
(57)

Since the mininum of Lipschitz functions (with the same constant) is also
Lipschitz (and similarly for maximum), the field $d^{**}$ belongs also to
$D_{n}\cap\mathrm{L}$. The field $d^{**}$ is a better solution than $d^{*}$
because it lowers the objective function:

$F(d^{**})\leq F(d^{*})$ (58)

Indeed, the objective function $F$ is an integral over ${\Omega}$ of a convex
local function with respect to $d$ denoted $f$ and we show that the inequality
also holds at the local level

$f(d^{**})\leq f(d^{*})$ (59)

We have two possible orderings

$d^{*}<d^{**}=\pi^{\mathrm{l}}\overline{d}\leq\overline{d}\quad\text{\ or\
}\quad\overline{d}\leq\pi^{\mathrm{u}}\overline{d}=d^{**}<d^{*}$ (60)

So, there exist ${\lambda}\in]0,1[$ such that

$d^{**}={\lambda}\overline{d}+(1-{\lambda})d^{*}$ (61)

leading to, by strict convexity of $f$

$f(d^{**})<{\lambda}f(\overline{d})+(1-{\lambda})f(d^{*})=f(d^{*})+{\lambda}(f(\overline{d})-f(d^{*}))<f(d^{*})$
(62)

where the last inequality is obtained from the fact that $\overline{d}$ at $x$
is the minimum of $f$ at $x$. Thus, to any field satisfying the constraints
but lying outside the bounds, we can associate a better one inside the bounds,
it proves that the optimal solution is inside the bounds.

##  Appendix B Nonlinear solver

The nonlinear solver used for (16) alternates between a local and global step.
The $m$ exponent is used to denote these iterations. The local step computes
at each integration point the stress ${\sigma}^{m}$ and internal variables
$\epsilon_{\mathrm{p}}^{m},p^{m}$, as well as the tangent operator $T^{m}$
knowing $u^{m},\epsilon_{\mathrm{p}}^{n},p^{n}$ and $d^{k}$. The $n$ exponent
indicates the previous known time-step values whereas the $k$ exponent
indicates the latest computed value of damage from (17) in the alternate
minimization.

The global step involves a linear solve which finds $u^{m+1}\in U_{n}$ such
that

$\int_{0}^{L}T^{m}\epsilon(u^{m+1}-u^{m})\epsilon(u^{*})\mbox{\,
$\mathrm{d}$}x=-\int_{0}^{L}{\sigma}^{m}\epsilon(u^{*})\mbox{\,
$\mathrm{d}$}x,\;\forall u^{*}\in U_{0}$

When the correction $u^{m+1}-u^{m}$ is below a threshold, we set
$u^{k+1},\epsilon_{\mathrm{p}}^{k+1},p^{k+1}$ as the current $m$ values. The
expressions of ${\sigma}^{m}$ and $T^{m}$ are given below for each model.

###  B.1 Softening elasticity

We have

${\sigma}^{m}=(1-d^{k})^{2}E{\epsilon}^{m},\quad{\epsilon}^{m}=\epsilon(u^{m}),\quad
T^{m}=(1-d^{k})^{2}E$ (63)

Since the model is linear, a single "$m$" iteration is needed to converge. If
dissymmetric tension-compression behavior was considered for damage, more
iterations would be needed to converge.

###  B.2 Softening elasticity with hardening plasticity

The trial stress is computed as

${\sigma}^{t}=E(1-d^{k})^{2}(\epsilon^{m}-\epsilon_{\mathrm{p}}^{n})$

followed by the plasticity criterion

$f^{t}=\lvert{\sigma}^{t}\rvert-\sigma_{\mathrm{y}}(1-d^{k})^{2}(1+kp^{n})$

giving two cases

if $\displaystyle
f^{t}<=0\quad\epsilon_{\mathrm{p}}^{m}=\epsilon_{\mathrm{p}}^{n},\;p^{m}=p^{n},\;{\sigma}^{m}={\sigma}^{t},T^{m}=E(1-d^{k})^{2}$
if $\displaystyle f^{t}>0\quad
p^{m}=\frac{\lvert{\sigma}^{t}\rvert-\sigma_{\mathrm{y}}+Ep^{n}}{E+\sigma_{\mathrm{y}}k},\quad
T^{m}=\frac{E\sigma_{\mathrm{y}}k(1-d^{k})^{2}}{E+\sigma_{\mathrm{y}}k}$
$\displaystyle\quad\epsilon_{\mathrm{p}}^{m}=\epsilon_{\mathrm{p}}^{n}+(p^{m}-p^{n})\frac{{\sigma}^{t}}{\lvert{\sigma}^{t}\rvert},\quad{\sigma}^{m}=E(1-d^{k})^{2}(\epsilon^{m}-\epsilon_{\mathrm{p}}^{m})$

###  B.3 Softening plasticity

The trial stress is computed as

${\sigma}^{t}=E(\epsilon^{m}-\epsilon_{\mathrm{p}}^{n}).$

followed by the plasticity criterion

$f^{t}=\lvert{\sigma}^{t}\rvert-\sigma_{\mathrm{y}}(1-d^{k})^{2}(1+kp^{n})$

and two cases

if $\displaystyle
f^{t}<=0:\quad\epsilon_{\mathrm{p}}^{m}=\epsilon_{\mathrm{p}}^{n},\;p^{m}=p^{n},\;{\sigma}^{m}={\sigma}^{t},\;T^{m}=E$
if $\displaystyle f^{t}>0:\quad
p^{m}=\frac{\lvert{\sigma}^{t}\rvert-\sigma_{\mathrm{y}}(1-d^{k})^{2}+Ep^{n}}{E+\sigma_{\mathrm{y}}k(1-d^{k})^{2}},\quad
T^{m}=\frac{E\sigma_{\mathrm{y}}k(1-d^{k})^{2}}{E+\sigma_{\mathrm{y}}k(1-d^{k})^{2}}$
$\displaystyle\quad\quad\quad\quad\epsilon_{\mathrm{p}}^{m}=\epsilon_{\mathrm{p}}^{n}+(p^{m}-p^{n})\frac{{\sigma}^{t}}{\lvert{\sigma}^{t}\rvert},\quad{\sigma}^{m}=E(\epsilon^{m}-\epsilon_{\mathrm{p}}^{m})$

## References

  * [1] Z. Bazant, T. Belytschko, T. Chang, “Continuum theory fo strain-softening”, _Journal of Engineering Mechanics_ 110 (1984), p. 1666-1692.
  * [2] G. Pijaudier-Cabot, Z. Bazant, “Non-local damage theory”, _Journal of Engineering Mechanics_ 113 (1987), p. 1512-1533.
  * [3] Z. Bazant, M. Jirasek, “Nonlocal Integral Formulations of Plasticity and Damage : Survey of Progress”, _Journal of Engineering Mechanics_ 128 (2002), no. November, p. 1119-1149.
  * [4] E. Lorentz, S. Andrieux, “Analysis of non-local models through energetic formulations”, in _International Journal of Solids and Structures_ , vol. 40, 2003, p. 2905-2936.
  * [5] E. C. Aifantis, “On the Microstructural Origin of Certain Inelastic Models”, _Journal of Engineering Materials and Technology_ 106 (1984), no. 4, p. 326-330, https://doi.org/10.1115/1.3225725, https://arxiv.org/abs/https://asmedigitalcollection.asme.org/materialstechnology/article-pdf/106/4/326/5788761/326_1.pdf.
  * [6] A. E. Triantafyllidis N, “A gradient approach to localization of deformation: I. Hyperelastic model”, _Journal of Elasticity_ 16 (1986), p. 225-237.
  * [7] C. Z. Schreyer H, “One-dimensional softening with localization”, _Journal of Applied Mechanics_ 53 (1986), p. 891-979.
  * [8] H. B. Mühlhaus, L. Vardoulakis, “The thickness of shear bands in granular materials”, _Geotechnique_ 37 (1987), p. 271-283, http://cat.inist.fr/?aModele=afficheN&cpsidt=7575585.
  * [9] M. Frémond, B. Nedjar, “Damage, gradient of damage and principle of virtual power”, _International Journal of Solids and Structures_ 33 (1996), no. 8, p. 1083-1103, http://linkinghub.elsevier.com/retrieve/pii/0020768395000747.
  * [10] G. Pijaudier-Cabot, N. Burlion, “Damage and localisation in elastic materials with voids”, _Cohesive-frictional Materials_ 144 (1996), p. 129-144, http://www3.interscience.wiley.com/journal/24869/abstract.
  * [11] R. Peerlings, M. Geers, R. De Borst, W. Brekelmans, “A critical comparison of nonlocal and gradient-enhanced softening continua”, _International Journal of Solids and Structures_ 38 (2001), no. 44-45, p. 7723-7746.
  * [12] Q.-S. Nguyen, S. Andrieux, “The non-local generalized standard approach: a consistent gradient theory”, _Comptes rendus de l’Académie des sciences: Mécanique, physique, chimie, astronomie_ 333 (2005), p. 139-145, http://linkinghub.elsevier.com/retrieve/pii/S1631072104002244.
  * [13] G. A. Francfort, J.-J. Marigo, “Revisiting brittle fracture as an energy minimization problem”, _Journal of the Mechanics and Physics of Solids_ 46 (1998), p. 1319-1412.
  * [14] B. Bourdin, G. A. Francfort, J.-J. Marigo, “Numerical experiments in revisited brittle fracture”, _Journal of the Mechanics and Physics of Solids_ 48 (2000), no. 4, p. 797-826, http://linkinghub.elsevier.com/retrieve/pii/S0022509699000289.
  * [15] B. Bourdin, G. a. Francfort, J.-J. Marigo, “The Variational Approach to Fracture”, _Journal of Elasticity_ 91 (2008), no. 1-3, p. 5-148, http://link.springer.com/10.1007/s10659-007-9107-3.
  * [16] A. Karma, D. Kessler, H. Levine, “Phase-Field Model of Mode III Dynamic Fracture”, _Physical Review Letters_ 87 (2001), no. 4, p. 045501, http://link.aps.org/doi/10.1103/PhysRevLett.87.045501.
  * [17] V. Hakim, A. Karma, “Laws of crack motion and phase-field models of fracture”, _Journal of the Mechanics and Physics of Solids_ 57 (2009), no. 2, p. 342-368, http://linkinghub.elsevier.com/retrieve/pii/S0022509608001774.
  * [18] H. Amor, J.-J. Marigo, C. Maurini, “Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments”, _Journal of the Mechanics and Physics of Solids_ 57 (2009), no. 8, p. 1209-1229, http://linkinghub.elsevier.com/retrieve/pii/S0022509609000659.
  * [19] C. Miehe, F. Welschinger, M. Hofacker, “Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations”, _International Journal For Numerical Methods in Engineering_ 83 (2010), no. 10, p. 1273-1311.
  * [20] C. Miehe, M. Hofacker, F. Welschinger, “A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits”, _Computer Methods in Applied Mechanics and Engineering_ 199 (2010), no. 45-48, p. 2765-2778, http://linkinghub.elsevier.com/retrieve/pii/S0045782510001283.
  * [21] C. Kuhn, R. Müller, “A continuum phase field model for fracture”, _Engineering Fracture Mechanics_ 77 (2010), no. 18, p. 3625-3634, http://linkinghub.elsevier.com/retrieve/pii/S0013794410003668.
  * [22] M. Ambati, T. Gerasimov, L. De Lorenzis, “Phase-field modeling of ductile fracture”, _Computational Mechanics_ 55 (2015), no. 5, p. 1017-1040, http://dx.doi.org/10.1007/s00466-015-1151-4.
  * [23] N. Moës, C. Stolz, P.-E. Bernard, N. Chevaugeon, “A level set based model for damage growth : the thick level set approach”, _International Journal For Numerical Methods in Engineering_ 86 (2011), p. 358-380.
  * [24] C. Stolz, N. Moës, “A new model of damage : a moving thick layer approach.”, _International Journal of Fracture_ 174 (2012), no. 1, p. 49-60, http://link.springer.com/10.1007/s10704-012-9693-3.
  * [25] N. Moës, C. Stolz, N. Chevaugeon, “Coupling local and non-local damage evolution with The Thick Level Set model”, _Advanced Modelling and Simulation in Engineering Sciences_ 2 (2014), no. 16, p. 21.
  * [26] A. J. Stershic, J. E. Dolbow, N. Moës, “The Thick Level-Set Model for Dynamic Fragmentation”, _Engineering Fracture Mechanics_ 172 (2017), p. 39-60, http://dx.doi.org/10.1016/j.engfracmech.2016.12.012.
  * [27] M. Frémond, C. Stolz, “On Alternative Approaches for Graded Damage Modelling”, in _Models, Simulation, and Experimental Issues in Structural Mechanics_ (Cham) (M. Frémond, F. Maceri, G. Vairo, eds.), Springer Series in Solid and Structural Mechanics, vol 8, Springer International Publishing, 2017, p. 87-104, https://doi.org/10.1007/978-3-319-48884-4_4.
  * [28] N. Valoroso, C. Stolz, “Progressive Damage in Quasi-brittle Solids”, in _Proceedings of XXIV AIMETA Conference 2019_ (A. Carcaterra, A. Paolone, G. Graziani, eds.), AIMETA 2019, springer, 2020, https://doi.org/10.1007/978-3-030-41057-5_34.
  * [29] P. Mialon, “ElÃ©ments d’analyse et de rÃ©solution numÃ©rique des relations de l’Ã©lasto-plasticitÃ©”, _Bulletin de la Direction des Ã©tudes et recherches - ElectricitÃ© de France. SÃ©rie C, mathÃ©matiques, informatique_ (1986) (French).
  * [30] M. Ortiz, L. Stainier, “The variational formulation of viscoplastic constitutive updates”, _Computer Methods in Applied Mechanics and Engineering_ 171 (1999), no. 3, p. 419 - 444, http://www.sciencedirect.com/science/article/pii/S0045782598002199.
  * [31] L. Boyd, Stephen P.; Vandenberghe, _Convex Optimization_.
  * [32] V. Krishnan, A. A. A. Makdah, F. Pasqualetti, “Lipschitz Bounds and Provably Robust Training by Laplacian Smoothing”, 2020, https://arxiv.org/abs/2006.03712.
  * [33] B. Halphen, Q.-S. Nguyen, “Sur les matériaux standards généralisés”, _Journal de Mécanique_ 14 (1975), no. 1, p. 39-63.
  * [34] P. Germain, P. Suquet, Q. S. Nguyen, “Continuum thermodynamics”, _ASME Journal of Applied Mechanics_ 50 (1983), p. 1010-1020.
  * [35] J. Simo, T. Hughes, _Computational inelasticity_ , Springer, 1997, 392 pages.
  * [36] J. Zghal, K. Moreau, N. Moës, D. Leguillon, C. Stolz, “Analysis of the failure at notches and cavities in quasi-brittle media using the Thick Level Set damage model and comparison with the coupled criterion”, _International Journal of Fracture_ 211 (2018), no. 1-2, p. 253-280, https://doi.org/10.1007/s10704-018-0287-6.
  * [37] E. Lorentz, V. Godard, “Gradient damage models: Toward full-scale computations”, _Computer Methods in Applied Mechanics and Engineering_ 200 (2011), no. 21-22, p. 1927-1944, http://linkinghub.elsevier.com/retrieve/pii/S0045782510001921.
  * [38] A. Parrilla Gómez, N. Moës, C. Stolz, “Comparison between thick level set (TLS) and cohesive zone models”, _Advanced Modeling and Simulation in Engineering Sciences_ 2 (2015), no. 1, p. 18, http://www.amses-journal.com/content/2/1/18.
  * [39] V. Tvergaard, A. Needleman, “Analysis of the cup-cone fracture in a round tensile bar”, _Acta Metallurgica et Materialia_ 32 (1984), p. 157-169.
  * [40] R. Becker, A. Needleman, O. Richmond, V. Tvergaard, “Void growth and failure in notched bars”, _Journal of the Mechanics and Physics of Solids_ 36 (1988), no. 3, p. 317-351.
  * [41] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”, _Nature Methods_ 17 (2020), p. 261-272.
  * [42] D. Kraft _et al._ , “A software package for sequential quadratic programming”, (1988).
  * [43] M. A. Crisfield, “A Fast Incremental /Iterative solution procedure that handles "snap-through"”, _Computers & Structures_ 13 (1981), p. 55-62.
  * [44] R. De Borst, “Computation of post-bifurcation and post-failure behavior of strain softening solids”, _Computers & Structures_ 2 (1987), no. 2, p. 211-224.