# Analytical and numerical relaxation results for models in soil mechanics

Florian Behr  Faculty of Mathematics; University of Regensburg, 93053
Regensburg, Germany Georg Dolzmann  Faculty of Mathematics; University of
Regensburg, 93053 Regensburg, Germany Klaus Hackl  Institute of Mechanics of
Materials; Ruhr University Bochum, 44780 Bochum, Germany  Ghina Jezdan
Institute of Mechanics of Materials; Ruhr University Bochum, 44780 Bochum,
Germany

###### Abstract

A variational model of pressure-dependent plasticity employing a time-
incremental setting is introduced. A novel formulation of the dissipation
potential allows one to construct the condensed energy in a variationally
consistent manner. For a one-dimensional model problem, an explicit expression
for the quasiconvex envelope can be found which turns out to be essentially
independent of the original pressure-dependent yield surface. The model
problem can be extended to higher dimensions in an empirical manner. Numerical
simulation exhibit well-posed behavior showing mesh-independent results.

##  1 Introduction

Classical models for soils and granular materials are formulated in the
framework of linearized plasticity based on pressure-dependent flow rules.
While this approach works very well in engineering applications, it cannot
explain complex structures seen in experiments conducted in granular media
like clay, silt or sand, see in particular [39, Figure 12]. The aim of this
paper is both the investigation and the modeling of these complex structures
via relaxation methods, a well-established approach based on variational
analysis. In a first step, the analysis will be based on a time discretization
scheme which involves a variational formulation for the problems at given
finite time-increments. The corresponding variational problems will be studied
concerning their stability properties and the main features of minimizing
microstructures, which are closely related to the complexity of the
experimental data. Moreover, the dependence of the main features of the
minimizing structures on the given yield functions will be investigated. The
goal is to employ the variational approach in order to describe the essential
model features with as few parameters as possible, thus establishing this type
of models for plasticity in soil mechanics.

Starting with [1, 2, 11] it was realized that experimentally observed
microstructures may be explained based on energy minimization. This
observation made it possible to use variational methods in order to understand
complex material behavior. For example, the shape memory effect is related to
an accommodation of macroscopic deformations on a microscopic scale through
complex patterns. These patterns are determined through mathematical
compatibility conditions which require in the nonlinear theory that two
adjacent elastic strains $\mathbf{A}$ and $\mathbf{B}$ be rank-one connected
in the sense that $\mathbf{A}-\mathbf{B}=\mathbf{a}\otimes\mathbf{n}$. The
vector $\mathbf{n}$ provides geometric information since it represents the
orientation of planes that can connect the two states $\mathbf{A}$ and
$\mathbf{B}$. Minimization of the original energy with respect to such
compatible microstructures leads to the relaxed energy. This approach allows
one on the one hand to formulate well-posed boundary-value problems that lend
themselves to numerical solutions. The minimizing microstructures, on the
other hand, will provide information on the expected deformation patterns.
Applications of this approach specifically to elastoplastic materials can be
found in [25, 21].

Since the explicit construction of microstructures requires to determine rank-
one connections, which is algebraically a challenging problem, it is not
surprising that there are only a few examples with applications to elastic
materials for which a closed form of $W^{qc}$ has been found [10, 17, 37, 38,
12, 23, 24]. In the linearized setting, the compatibility condition concerns
symmetric stains $\mathbf{A}$ and $\mathbf{B}$ whose difference is symmetrized
rank-one convex, i.e., $\mathbf{A}-\mathbf{B}=\mathbf{a}\odot\mathbf{b}$ with
$\mathbf{a}\odot\mathbf{b}=\frac{1}{2}(\mathbf{a}\otimes\mathbf{b}+\mathbf{b}\otimes\mathbf{a})$.
Examples for relaxation results in this context can be found in [27, 29, 26,
36]. Numerical relaxation schemes are employed in [6]. A combination of an
analytical approximation of the relaxed energy density combined with a finite
element simulation can be found in [8, 30, 31]. A novel algorithmic approach
for the computation of the relaxed energy density with applications to
problems in phase transformations and models in plasticity was recently
proposed in [13, 14, 15], see also [4, 5].

The paper is organized as follows: In Sec. 2 standard models of pressure
dependent plasticity are reviewed. In Sec. 3 we bring these models into a
variational setting suitable for treatment within relaxation theory. In Sec. 4
this setting is specialized in order to obtain a one-dimensional model problem
for which analytical results can be obtained. Sec. 5 contains the main results
of this paper from the mathematical point of view. An explicit expression for
the relaxed energy is derived in the one-dimensional case, i.e., the
quasiconvex envelope of the energy density, which, in the given situation with
one spatial variable, is equal to the convex envelope. In Sec. 6 the one-
dimensional relaxed energy is, in a heuristic manner, extended to three
dimensions and in Sec. 7 numerical results are presented for both models,
first for the relaxed energy derived for the one-dimensional problem, then for
an extension of this energy to three-dimensions, for which no rigorous
mathematical results are available yet. More extensive results will be given
in [22]. We will close with conclusion and outlook in Sec. 8.

##  2 Pressure-dependent plasticity

The most important class of materials exhibiting pressure-dependent plastic
behavior are soils and granular media. They are usually heterogeneous mixtures
of fluids (usually air and water) and particles (usually clay, silt, sand, and
gravel) that have little to no cementation [33]. The shear strength is
provided by friction and interlocking of the particles.

The understanding of deformations of solid materials is of key importance in
applications. In contrast to elastic materials, which return to their original
configuration upon unloading, elasto-plastic materials may undergo a permanent
change of their configuration if large shear strains occur in the material; in
contrast, it is frequently assumed that they may support arbitrary hydrostatic
pressures. Usually this material behaviour is modeled in a linearized setting
via a convex set $\mathbb{K}\subset{\mathbb{R}}^{n\times n}_{\mathrm{sym}}$
which describes the set of admissible stresses in the material which lead to
elastic deformations. Only for stresses on $\partial\mathbb{K}$, permanent
plastic deformation occurs and this deformation is formulated via a flow rule
which determines the evolution of the plastic part of the total elastoplastic
deformation. Consequently, the set $\mathbb{K}$ is an infinite cylinder

$\displaystyle\mathbb{K}=\bigcup\limits_{\pi\in{\mathbb{R}}}\pi\mathbf{Id}\oplus
K=\bigcup\limits_{\pi=\pi\in{\mathbb{R}}}^{\pi_{\mathrm{max}}}\\{\pi\mathbf{Id}+\mathbf{E},\mathbf{E}\in
K\\}\,,$ (2.1)

where $K$ is a compact subset in ${\mathbb{R}}^{n\times n}_{\mathrm{D}}$ and

$\displaystyle{\mathbb{R}}^{n\times
n}_{\mathrm{D}}=\\{\mathbf{E}\in{\mathbb{R}}^{n\times
n}\colon\mathbf{E}=\mathbf{E}^{T},\,\operatorname{tr}(\mathbf{E})=0\\}$ (2.2)

is the set of all symmetric and trace free matrices. Here and in the following
we will use boldface symbols for vector-valued or matrix-valued objects.

Thus the basic problem in the theory of linear elastoplasticity can be
formulated as follows. Given a material body described in its reference
configuration by a domain $\Omega\subset{\mathbb{R}}^{n}$, find a deformation
$\mathbf{y}:\Omega\to{\mathbb{R}}^{n}$ and a decomposition of the symmetrized
displacement gradient $\boldsymbol{\epsilon}$ defined by

$\displaystyle\mathbf{y}(\mathbf{x})=\mathbf{id}(\mathbf{x})+\mathbf{u}(\mathbf{x})\,,\quad\boldsymbol{\epsilon}(\mathbf{x})=\frac{1}{2}(D\mathbf{u}(\mathbf{x})+D\mathbf{u}(\mathbf{x})^{T})$
(2.3)

into an elastic and a plastic part,
$\boldsymbol{\epsilon}=\boldsymbol{\epsilon}_{e}+\boldsymbol{\epsilon}_{p}$,
$\boldsymbol{\epsilon}_{p}\in{\mathbb{R}}^{n\times n}_{\mathrm{D}}$, such that
the following system of equations holds:

  * (i)

constitutive equation
$\boldsymbol{\sigma}=\mathbb{C}\boldsymbol{\epsilon}_{e}$,

  * (ii)

balance of momentum
$\mathbf{u}_{tt}-\operatorname{div}\boldsymbol{\sigma}=\mathbf{f}$ or the
equilibrium condition $-\operatorname{div}\boldsymbol{\sigma}=\mathbf{f}$,

  * (iii)

stress constraint $\boldsymbol{\sigma}\in\mathbb{K}$,

  * (iv)

associated flow rule
$(\boldsymbol{\xi}-\boldsymbol{\sigma}):\dot{\boldsymbol{\epsilon}}_{p}\leq 0$
for every $\boldsymbol{\xi}\in\mathbb{K}$.

It is important to note that the assumption (iv) is not realistic for
materials in which frictional effects dominate as in clay, sands, gravel or
other granular materials. Here pressure-dependence needs to be taken into
account and this leads to a flow rule that cannot be modeled by an infinite
cylinder; instead one is led to consider cones in the space of principle
stresses, classical models being the Drucker-Prager and Mohr-Coulomb criteria
[28]. To be specific, let $\pi=-1/3\,\operatorname{tr}\boldsymbol{\sigma}$
denote the pressure, then the convex elastic domain is given as

$\displaystyle\mathbb{K}=\bigcup\limits_{\pi=\pi_{\mathrm{min}}}^{\pi_{\mathrm{max}}}\pi\mathbf{Id}\oplus
K_{\pi}=\bigcup\limits_{\pi=\pi_{\mathrm{min}}}^{\pi_{\mathrm{max}}}\\{\pi\mathbf{Id}+\mathbf{E},\mathbf{E}\in
K_{\pi}\\}\,,$ (2.4)

where $K_{\pi}$ is a compact subset in ${\mathbb{R}}^{n\times n}_{\mathrm{D}}$
and $-\infty\leq\pi_{\mathrm{min}}<\pi_{\mathrm{max}}\leq\infty$. In order to
comply with experimental observations, plastic deformation should still be
volume preserving. For this purpose, the flow rule (iv) is then replaced by
the non-associated flow rule and the full system of equations in its
quasistatic form is given by

  * (i)

constitutive equation
$\boldsymbol{\sigma}=\mathbb{C}\boldsymbol{\epsilon}_{e}$,

  * (ii)

the equilibrium condition $-\operatorname{div}\boldsymbol{\sigma}=\mathbf{f}$,

  * (iii)

stress constraint $\boldsymbol{\sigma}\in\mathbb{K}$,

  * (iv’)

non-associated flow rule
$(\boldsymbol{\xi}-\boldsymbol{\sigma}):\dot{\boldsymbol{\epsilon}}_{p}\leq 0$
for every $\xi\in\mathbb{K}_{\pi}$.

Note that $\dot{\boldsymbol{\epsilon}}_{p}$ is now normal to the boundary of
$\mathbb{K}_{\pi}$ as subset of ${\mathbb{R}}^{n\times n}_{\mathrm{D}}$ but
not normal to $\mathbb{K}$ as subset of ${\mathbb{R}}^{n\times
n}_{\mathrm{sym}}$. This observation lies at the heart of the problems
concerning a variational formulation to be discussed in the following.

##  3 Variational model

In order to be able to apply relaxation theory, we have to reformulate the
setting given in the previous section. For this purpose, we will employ a
variational approach for the description of inelastic processes, see [34, 9,
32, 35, 19, 20]. Let us consider a physical system described by (sets of)
external, i.e., controllable, state variables, in our case given by the total
strain $\boldsymbol{\varepsilon}$, and internal state variables $\mathbf{z}$.
Here $\mathbf{z}=\\{\boldsymbol{\varepsilon}_{\mathrm{p}},p\\}$, where
$\boldsymbol{\varepsilon}_{\mathrm{p}}$ denotes plastic strain and $p\geq 0$ a
hardening variable. We assume that the system behavior may be defined using
only two scalar potentials: free energy
$\psi(\boldsymbol{\varepsilon},\mathbf{z})$ and dissipation potential
$\Delta(\boldsymbol{\varepsilon},\mathbf{z},\dot{\mathbf{z}})$.

The evolution of the internal variables is described by

${\mathbf{0}}\in\frac{\partial\psi}{\partial{\mathbf{z}}}+\frac{\partial\Delta}{\partial\dot{\mathbf{z}}}.$
(3.1)

Note, that Eq. (3.1) may be written as a stationarity condition of the
minimization problem

$\dot{\mathbf{z}}=\arg\inf_{\dot{\mathbf{z}}}\left\\{{\dot{\psi}+\Delta}\right\\}.$
(3.2)

We assume a free energy of the form

$\psi=\frac{K}{2}\,{\operatorname{tr}\boldsymbol{\varepsilon}}^{2}+\mu\left\|\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p}}\right\|^{2}+\rho(\operatorname{tr}\boldsymbol{\varepsilon})\,p+\frac{\beta}{2}\left\|\boldsymbol{\varepsilon}_{\mathrm{p}}\right\|^{2}.$
(3.3)

Here, $p$ is an auxiliary internal variable, $K$ the bulk-modulus, $\mu$ the
shear-modulus, and $\beta$ a hardening-modulus. The function
$\rho(\operatorname{tr}\boldsymbol{\varepsilon})$ governs pressure-dependence.
For constant $\rho(\operatorname{tr}\boldsymbol{\varepsilon})$, Eq. (3.3)
corresponds to the free energy of an isotropic elastoplastic material with
kinematic hardening. In this work, we will assume $\rho$ to be a concave
function inside an interval $[\xi_{\mathrm{min}},\xi_{\mathrm{max}}]$ and zero
outside of it. The function $\rho(\operatorname{tr}\boldsymbol{\varepsilon})$
for a typical model, e.g. capped Drucker-Prager, see [33], can be seen in Fig.
1. The stress is given as

Figure 1: typical pressure-dependent yield surface.

$\boldsymbol{\sigma}=\frac{\partial\psi}{\partial\boldsymbol{\varepsilon}}=(K\,\operatorname{tr}\boldsymbol{\varepsilon}+\rho^{\prime}(\operatorname{tr}\boldsymbol{\varepsilon})\,p)\mathbf{I}+2\mu\left(\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p}}\right).$
(3.4)

Note, that the pressure is given by the relation
$\pi=-K\operatorname{tr}\boldsymbol{\varepsilon}-\rho^{\prime}(\operatorname{tr}\boldsymbol{\varepsilon})p$.
Here, the second term may be interpreted as a consolidation pressure due to an
inelastic volume change governed by the parameter $p$. A variationally more
consistent approach would be to replace Eq. (3.3) by an energy of the form

$\psi_{\mathrm{consol}}=\frac{K}{2}\,(\operatorname{tr}\boldsymbol{\varepsilon}-\theta_{\mathrm{p}})^{2}+\mu\left\|\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p}}\right\|^{2}+\rho(\theta_{\mathrm{p}})\,p+\frac{\beta}{2}\left\|\boldsymbol{\varepsilon}_{\mathrm{p}}\right\|^{2},$
(3.5)

where $\theta_{\mathrm{p}}$ now represents the inelastic volume change
explicitly. However, we would like to keep our formulation simple within the
present work and plan to discuss this formulation in a subsequent paper.

The dissipation potential will be assumed as

$\Delta(\dot{\boldsymbol{\varepsilon}}_{\mathrm{p}},\dot{p})=\begin{cases}0&\text{for
}\|\dot{\boldsymbol{\varepsilon}}_{\mathrm{p}}\|\leq\dot{p}\\\
\infty&\text{otherwise}\end{cases}.$ (3.6)

This procedure differs from the standard one, as for example implemented in
[9]. In the standard approach, Eq. (3.3) is replaced by

$\psi_{\mathrm{alt}}=\frac{K}{2}\,{\operatorname{tr}\boldsymbol{\varepsilon}}^{2}+\mu\left\|\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p}}\right\|^{2}+p+\frac{\beta}{2}\left\|\boldsymbol{\varepsilon}_{\mathrm{p}}\right\|^{2},$
(3.7)

and Eq. (3.6) by

$\Delta_{\mathrm{alt}}(\operatorname{tr}\boldsymbol{\varepsilon},\dot{\boldsymbol{\varepsilon}}_{\mathrm{p}},\dot{p})=\rho(\operatorname{tr}\boldsymbol{\varepsilon})\begin{cases}0&\text{for
}\|\dot{\boldsymbol{\varepsilon}}_{\mathrm{p}}\|\leq\dot{p}\\\
\infty&\text{otherwise}\end{cases}.$ (3.8)

However, the formulation in Eq. (3.8) would be problematic with respect to the
variational formulation which does not allow the dissipation potential to
depend on an external variable. Moreover, it would lead to complications when
it comes to integrating it up in order to obtain a dissipation distance, as
will be done later in the text. Minimizing trajectories connecting different
states will be non-monotone leading to behavior that is questionable from a
physical point of view.

Minimization in Eq. (3.2) gives, taking into account the expressions in Eqs.
(3.3) and (3.6),

$\mathrm{dev}\boldsymbol{\sigma}-\beta\boldsymbol{\varepsilon}_{\mathrm{p}}\in\rho(\operatorname{tr}\boldsymbol{\varepsilon})\operatorname{sign}\dot{\boldsymbol{\varepsilon}}_{\mathrm{p}},$
(3.9)

and

$\dot{p}=\|\dot{\boldsymbol{\varepsilon}}_{\mathrm{p}}\|.$ (3.10)

Note, that Eq. (3.9) is equivalent to the yield condition

$\Phi:=\|\mathrm{dev}\boldsymbol{\sigma}-\beta\boldsymbol{\varepsilon}_{\mathrm{p}}\|-\rho(\operatorname{tr}\boldsymbol{\varepsilon})\leq
0,$ (3.11)

and the flow rule

$\dot{\boldsymbol{\varepsilon}}_{\mathrm{p}}=\kappa\,(\mathrm{dev}\boldsymbol{\sigma}-\beta\boldsymbol{\varepsilon}_{\mathrm{p}}),$
(3.12)

where $\kappa\geq 0$ is a consistency parameter satisfying the Kuhn-Tucker
conditions

$\kappa\geq 0,\qquad\Phi\leq 0,\qquad\kappa\,\Phi=0.$ (3.13)

Note, that Eq. (3.10) defines $p$ as equivalent plastic strain.

Having established a variational scheme for pressure-dependent plasticity now,
we will employ a time-incremental approach as proposed in [9, 35] in order to
investigate the evolution of material microstructures. For this purpose, we
consider a time-increment $\Delta t=[t_{n},t_{n+1}]$. Given the external
forces and boundary conditions at time $t_{n+1}$ and the values of the
internal variables $\mathbf{z}_{n}$ at time $t_{n}$, we seek the values
$\boldsymbol{\varepsilon}_{n+1}$ and $\mathbf{z}_{n+1}$ at time $t_{n+1}$.

A time-incremental version of the dissipation-potential is given by the
dissipation-distance The extremal principle Eq. (3.2) can now be replaced by

$\\{\boldsymbol{\varepsilon}_{\mathrm{p},n+1},p_{n+1}\\}=\arg\inf_{\boldsymbol{\varepsilon}_{\mathrm{p}},p}\left\\{{\psi(\boldsymbol{\varepsilon},\boldsymbol{\varepsilon}_{\mathrm{p}},p)+D(\boldsymbol{\varepsilon}_{\mathrm{p},n},p_{n},\boldsymbol{\varepsilon}_{\mathrm{p}}},p)\right\\}.$
(3.14)

Minimization with respect to $p$ gives immediately

$p-p_{n}=\|\boldsymbol{\varepsilon}_{\mathrm{p}}-\boldsymbol{\varepsilon}_{\mathrm{p},n}\|.$
(3.15)

Then, substitution of Eq. (3.15) into Eq. (3.14) and minimization with respect
to $\boldsymbol{\varepsilon}_{\mathrm{p}}$ yields

$2\mu(\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p}})-\beta\boldsymbol{\varepsilon}_{\mathrm{p}}\in\rho(\operatorname{tr}\boldsymbol{\varepsilon})\operatorname{sign}(\boldsymbol{\varepsilon}_{\mathrm{p}}-\boldsymbol{\varepsilon}_{\mathrm{p},n}).$
(3.16)

Solving for $\boldsymbol{\varepsilon}_{\mathrm{p}}$ in Eq. (3.16)yields

$\boldsymbol{\varepsilon}_{\mathrm{p}}=\boldsymbol{\varepsilon}_{\mathrm{p},n}+\frac{1}{2\mu+\beta}\,\left[2\mu\|\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p},n}\|-\rho(\operatorname{tr}\boldsymbol{\varepsilon})\right]_{+}\operatorname{sign}(\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p},n}),$
(3.17)

where $\left[x\right]_{+}=\max\\{0,x\\}$. This allows to introduce the
condensed energy, see [9, 35], via

$\psi_{\mathrm{cond}}(\boldsymbol{\varepsilon},\boldsymbol{\varepsilon}_{\mathrm{p},n},p_{n})=\inf_{\boldsymbol{\varepsilon}_{\mathrm{p}},p}\left\\{{\psi(\boldsymbol{\varepsilon},\boldsymbol{\varepsilon}_{\mathrm{p}},p)+D(\boldsymbol{\varepsilon}_{\mathrm{p},n},p_{n},\boldsymbol{\varepsilon}_{\mathrm{p}}},p)\right\\}.$
(3.18)

Substitution of Eqs. (3.15) and (3.17) into Eq. (3.18) gives, up to an
irrelevant constant term of the form
$\frac{\beta}{2}\left\|\boldsymbol{\varepsilon}_{\mathrm{p},n}\right\|^{2}$,

$\psi_{\mathrm{cond}}(\boldsymbol{\varepsilon},\boldsymbol{\varepsilon}_{\mathrm{p},n},p_{n})=\frac{K}{2}\,{\operatorname{tr}\boldsymbol{\varepsilon}}^{2}+\mu\left\|\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p},n}\right\|^{2}+\rho(\operatorname{tr}\boldsymbol{\varepsilon})p_{n}\\\
-\frac{1}{2}\,\frac{1}{2\mu+\beta}\,\left[2\mu\|\mathrm{dev}\boldsymbol{\varepsilon}-\boldsymbol{\varepsilon}_{\mathrm{p},n}\|-\rho(\operatorname{tr}\boldsymbol{\varepsilon})\right]_{+}^{2}.$
(3.19)

The condensed energy can now be considered as formally elastic and dependent
on the parameters $\boldsymbol{\varepsilon}_{\mathrm{p},n},p_{n}$, governing
the behavior during the time-increment $\Delta t=[t_{n},t_{n+1}]$. At the end
of the time-step the internal variables can be updated using Eqs. (3.15) and
(3.17). This way, the inelastic behavior can be approximated by solving a
sequence of elastic problems.

##  4 A one-dimensional model problem

It is well-known that quasiconvex envelopes can be calculated analytically
only in rare special cases. In order to circumvent this difficulty, we
restrict ourselves to a one-dimensional problem for which the quasiconvex
envelope coincides with the convex envelope. We consider a state of specific
deformation corresponding to a thin sheet of material under uniform tension or
compression and shear. We assume the sheet to be supported by stiff fibers in
a way that prevents bending deformations. The non-vanishing displacements and
strains are in Cartesian coordinates

$u_{x}=u(x),\qquad u_{y}=v(x)+y\,u^{\prime}(x),\qquad u_{z}=z\,u^{\prime}(x)$
(4.1)
$\varepsilon_{xx}=\varepsilon_{yy}=\varepsilon_{zz}=u^{\prime},\qquad\varepsilon_{yx}=\varepsilon_{xy}=\frac{1}{2}\,v^{\prime},$
(4.2)

where we neglected terms of order $u^{\prime\prime}$. For the plastic strains,
we assume that only the components
$\varepsilon_{\mathrm{p}}:=\varepsilon_{\mathrm{p}yx}=\varepsilon_{\mathrm{p}xy}$
are different from zero. The free energy assumes the form

$\psi=\frac{K}{2}\,\varepsilon_{xx}^{2}+\mu\left(\left(\varepsilon_{yx}-\varepsilon_{\mathrm{p}yx}\right)^{2}+\left(\varepsilon_{xy}-\varepsilon_{\mathrm{p}xy}\right)^{2}\right)+\rho(\varepsilon_{xx})\,p+\beta\,\varepsilon_{\mathrm{p}xy}^{2}.$
(4.3)

Substitution of Eq. (4.2) into Eq. (4.3) gives

$\psi=\frac{K}{2}\left(u^{\prime}\right)^{2}+2\,\mu\left(\frac{1}{2}v^{\prime}-\varepsilon_{\mathrm{p}}\right)^{2}+\rho\left(u^{\prime}\right)\,p+\beta\,\varepsilon_{\mathrm{p}}^{2}.$
(4.4)

As before, minimization with respect to $\varepsilon_{\mathrm{p}}$ and $p$
once again leads to the condensed energy

$\psi_{\mathrm{cond}}(u^{\prime},v^{\prime},\varepsilon_{\mathrm{p},n},p_{n})=\frac{K}{2}\,\left(u^{\prime}\right)^{2}+2\,\mu\,\left(\frac{1}{2}\,v^{\prime}-\varepsilon_{\mathrm{p},n}\right)^{2}+\rho\left(u^{\prime}\right)p_{n}\\\
-\frac{1}{2}\,\frac{1}{2\mu+\beta}\,\left[2\mu\left|\frac{1}{2}\,v^{\prime}-\varepsilon_{\mathrm{p},n}\right|-\rho\left(u^{\prime}\right)\right]_{+}^{2}.$
(4.5)

For the remainder of the paper, we will focus on the case $p_{n}=0$. This will
simplify the mathematical treatment significantly and will allow to keep the
derivations stringent and instructive. In mechanical terms this means, that we
will describe the initiation of microstructures rather than their evolution.
However, the general case will remain subject of future investigations.

It is convenient for the mathematical treatment following in the subsequent
section to introduce a dimensionless formulation. For this purpose, we define
modified variables and quantities by

$y_{1}=\sqrt{\frac{K}{2\mu}}\,u^{\prime},\quad
y_{2}=\frac{1}{2}\,v^{\prime},\quad
r(y_{1})=\frac{1}{2\mu}\,\rho\left(u^{\prime}\right),\quad
b=\frac{\beta}{2\mu},$ (4.6)

giving, using $p_{n}=0$,

$\psi_{\mathrm{cond}}(y_{1},y_{2},\varepsilon_{\mathrm{p},n},p_{n})=2\mu\,\left[\frac{1}{2}\,y_{1}^{2}+\frac{1}{2}\left(y_{2}-\varepsilon_{\mathrm{p},n}\right)^{2}-\frac{1}{2}\,\frac{1}{b+1}\,\left[\,|y_{2}-\varepsilon_{\mathrm{p},n}|-r(y_{1})\right]_{+}^{2}\right].$
(4.7)

Finally we may drop the multiplicative constant $2\mu$ and write
$z_{n}=\varepsilon_{\mathrm{p},n}$. Since $\psi_{\mathrm{cond}}$ is the free
energy density in the variational functional that needs to be minimized in the
first time step, it is crucial to understand the convexity properties of
$\psi_{\mathrm{cond}}$ in order to employ the direct method in the calculus of
variations. It turns out that the concavity of $r$ leads to a lack of
convexity of $\psi_{\mathrm{cond}}$ and the approach via relaxation [16]
requires one to characterize the quasiconvex envelope of
$\psi_{\mathrm{cond}}$, which, in the scalar case, reduces to the convex
envelope. This task is accomplished in the next section.

##  5 Characterization of the relaxed energy for the one-dimensional model

For any function $f:{\mathbb{R}}^{2}\to{\mathbb{R}}$ we denote the convex
envelope of $f$ with $f_{\mathrm{c}}$ and, if $f$ is differentiable at
$y\in{\mathbb{R}}^{2}$, the linear Taylor polynomial of $f$ about the point
$y$ with $T_{y}f$. Fix $b>0$ and
$y_{\mathrm{min}},y_{\mathrm{max}}\in{\mathbb{R}}$ with
$y_{\mathrm{min}}<y_{\mathrm{max}}$. Define the set of admissible dissipation
functions by

$\displaystyle\mathcal{R}=\\{r\in
C({\mathbb{R}},[0,\infty)),\,\operatorname{supp}(r)=[y_{\mathrm{min}},y_{\mathrm{max}}],r|_{[y_{\mathrm{min}},y_{\mathrm{max}}]}\text{
concave}\\}\,.$ (5.1)

Fix $y_{\mathrm{mid}}:=(y_{\mathrm{min}}+y_{\mathrm{max}})/2$ and define the
piecewise affine function $r_{0}\in\mathcal{R}$ by

$\displaystyle
r_{0}(y_{1})=\sqrt{b}\cdot\left(\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}-|y_{1}-y_{\mathrm{mid}}|\right)\cdot\chi_{[y_{\mathrm{min}},y_{\mathrm{max}}]}(y_{1})\,.$
(5.2)

The condensed energy for some $r\in\mathcal{R}$ is then given by

$\displaystyle
f^{(r)}_{z_{n}}(y)=\frac{1}{2}(y_{1}^{2}+(y_{2}-z_{n})^{2})-\frac{(|y_{2}-z_{n}|-r(y_{1}))_{+}^{2}}{2(b+1)}\,.$
(5.3)

Since $y_{2}$ appears only shifted by $z_{n}$ it is sufficient to calculate
the convex envelope $f_{\mathrm{c}}^{(r)}$ of $f^{(r)}:=f^{(r)}_{0}$ and
obtain the convex envelope of $f^{(r)}_{z_{n}}(y_{1},y_{2})$ by
$f_{\mathrm{c}}^{(r)}(y_{1},y_{2}-z_{n})$.
If the dissipation function is apparent from the context, we simply write $f$
instead of $f^{(r)}$. Define the disjoint sets

$\displaystyle\begin{aligned} \widetilde{Y}&=\\{y\in{\mathbb{R}}^{2}\colon
y_{1}\leq y_{\mathrm{min}}\text{ or }y_{1}\geq y_{\mathrm{max}}\\}\,,\\\
Y_{0}&=\\{y\in{\mathbb{R}}^{2}\colon
y_{1}\in(y_{\mathrm{min}},y_{\mathrm{max}}),\,|y_{2}|<r(y_{1})\\}\,,\\\
Y_{\pm}&=\\{y\in{\mathbb{R}}^{2}\colon
y_{1}\in(y_{\mathrm{min}},y_{\mathrm{max}}),\,\pm y_{2}\geq
r(y_{1})\\}\,,\end{aligned}$ (5.4)

whose union is ${\mathbb{R}}^{2}$, see Figure 2.

Figure 2: Definition of the condensed energy for the model problem. The set
$Y_{0}$ in Eq. (5.4) is enclosed by the blue curves ($\pm r$), the sets
$Y_{\pm}$ correspond to the sets above and below $Y_{0}$, respectively, and
the set $\widetilde{Y}$ is the complement of these three sets in
${\mathbb{R}}^{2}$.

###### Lemma 1

Consider two dissipation functions $r,\bar{r}\in\mathcal{R}$. If
$y_{1}\in{\mathbb{R}}$ and $r(y_{1})\leq\bar{r}(y_{1})$, then for all
$y_{2}\in{\mathbb{R}}$ the inequality $f^{(r)}(y_{1},y_{2})\leq
f^{(\bar{r})}(y_{1},y_{2})$ holds. In particular, if $r\leq\bar{r}$, then
$f^{(r)}\leq f^{(\bar{r})}$ and $g:{\mathbb{R}}^{2}\to{\mathbb{R}},\
g(y)=\frac{1}{2}y_{1}^{2}+\frac{1}{2}\frac{b}{b+1}y_{2}^{2}$ is a convex lower
bound of $f^{(r)}$, which coincides on $\widetilde{Y}$ with $f^{(r)}$. □

###### Proof

For any $y\in\mathbb{R}^{2}$ with $r(y_{1})\leq\bar{r}(y_{1})$ it follows that
$|y_{2}|\geq|y_{2}|-r(y_{1})\geq|y_{2}|-\overline{r}(y_{1})$. Therefore we
have
$y_{2}^{2}\geq(|y_{2}|-r(y_{1}))_{+}^{2}\geq(|y_{2}|-\overline{r}(y_{1}))_{+}^{2}$
and the inequality

$\displaystyle g(y)\leq f^{(r)}(y_{1},y_{2})\leq
f^{(\overline{r})}(y_{1},y_{2})$ (5.5)

is obtained from Eq. (5.3). ■

Since Lemma 1 establishes a lower bound on $f_{\mathrm{c}}$, we may use the
results in [18, Theorem 2.1]. For a lower semicontinuous function
$f:{\mathbb{R}}^{d}\to{\mathbb{R}}$, which satisfies the growth condition
$\lim_{|\delta|\to\infty}f(\delta)/|\delta|=\infty$, the convex envelope of
$f$ is a real-valued convex function and for any $y\in{\mathbb{R}}^{d}$ the
value $f_{\mathrm{c}}(y)$ is a convex combination of at most $d+1$ function
values $f(y_{i})$, $1\leq i\leq q\leq d+1$. Moreover, the subgradient
$\partial f_{\mathrm{c}}(y)$ is nonempty and for any $v\in\partial f(y)$

$\displaystyle f(y_{i})-\langle v,y_{i}\rangle=f(y)-\langle
v,y\rangle\,,\qquad i\in\\{1,...,q\\}$ (5.6)

and

$\displaystyle\forall\bar{y}\in{\mathbb{R}}^{d}:f(\bar{y})-\langle
v,\bar{y}\rangle\geq f(y_{i})-\langle v,y_{i}\rangle\,,\qquad
i\in\\{1,...,q\\}\,.$ (5.7)

Since Eq. (5.6) implies $f(y_{i})-\langle v,y_{i}\rangle=f(y_{j})-\langle
v,y_{j}\rangle$ for all $i,j\in\\{1,...,q\\}$ and Eq. (5.7) is equivalent to
$v\in\partial f(y_{i})$, a natural approach to calculate the convex envelope
of a continuous function $f:{\mathbb{R}}^{d}\to{\mathbb{R}}$ is therefore to
find pairwise distinct points $y_{1},...,y_{q}\in{\mathbb{R}}^{d}$ ($q\leq
d+1$), which satisfy the subgradient condition

$\displaystyle\exists v\in\bigcap\limits_{i=1}^{q}\partial f(y_{i}):\forall
i,j\in\\{1,...,q\\}:f(y_{i})-\langle v,y_{i}\rangle=f(y_{j})-\langle
v,y_{j}\rangle$ (SG)

and define for all $y\in\text{Conv}(\\{y_{1},...,y_{q}\\})$, the convex hull
of the points $y_{1},\ldots,y_{q}$, the convex envelope of $f$ by

$\displaystyle y=\sum\limits_{i=1}^{q}\lambda_{i}y_{i}\,,\quad
f_{\mathrm{c}}(y)=\sum\limits_{i=1}^{q}\lambda_{i}f(y_{i})\,.$ (5.8)

Finally, if for some direction $w\in{\mathbb{R}}^{d}$ and some
$y\in{\mathbb{R}}^{d}$ the directional derivative $D_{w}f(y)$ exists, then
$v\in\partial f(y)$ with $|v|=1$ implies $\langle v,w\rangle=D_{w}f(y)$.
Furthermore, if $f$ is differentiable in $y$, then $v\in\partial f(y)$ implies
$v=\nabla f(y)$. This leads to the observation that, if
$y_{1},...,y_{q}\in{\mathbb{R}}^{d}$ is a solution of Eq. (SG) and if $f$ is
differentiable in some $y_{i}$, then for all
$y\in\text{Conv}(\\{y_{1},...,y_{q}\\})$ the convex envelope of $f$ is given
by

$\displaystyle f_{\mathrm{c}}(y)=(T_{y_{i}}f)(y)=f(y_{i})+\langle\nabla
f(y_{i}),y-y_{i}\rangle\,.$ (5.9)

The calculation of the subgradient involves an evaluation of the values of $f$
on the whole space, and the solution of Eq. (SG) is in general challenging.

Based on these observations, for the construction of the convex envelope of
the condensed energy $f^{(r)}$ we first derive pairs and triples of points
$y_{1},...,y_{q}\in{\mathbb{R}}^{2}$ ($q=2,3$), which have equal partial
derivatives in the $y_{2}$ direction in which $f^{(r)}$ is differentiable.
Then we derive a candidate for the direction $v\in{\mathbb{R}}^{d}$ satisfying
the equation in Eq. (SG). However, at this point we do not know whether this
$v$ in Eq. (SG) is indeed a subgradient of $f$ in all the points
$y_{1},...,y_{q}$. For those pairs and triples for which we expect this to be
true, we do not show explicitly that the first part in Eq. (SG) is satisfied,
rather we apply the construction described in Eqs. (5.8) and (5.9) and prove
in Theorem 1 that this construction leads to the convex envelope. It is a
remarkable result of the subsequent analysis that the explicit formulas for
the relaxed energy show a very different behavior with respect to the
dependence on the dissipation function $r$. In fact, for small values of $b$
relative to $r$, i.e., for

$\displaystyle\sqrt{b}\cdot\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}\leq
r\Bigl{(}\frac{y_{\mathrm{min}}+y_{\mathrm{max}}}{2}\Bigr{)}\,,$ (5.10)

which is the case of relevance in soil mechanics, the relaxed energy is in
fact independent of $r$. Therefore we restrict the analysis in this paper to
the case of small $b$, the general case will be treated in [7]. In fact, this
assumption implies $r_{0}\leq r$ by concavity of $r$. In the following we
consider $r$ as fixed, write $f=f^{(r)}$ and assume without loss of generality
that $r_{0}<r$ on $(y_{\mathrm{min}},y_{\mathrm{max}})$.

In order to implement the strategy just outlined, we construct a candidate
$f_{\mathrm{a}}$ for $f_{\mathrm{c}}$. In a first step, we search for points
with common supporting planes in $Y_{0}$ and $\widetilde{Y}$, see Figure 2.
Assume thus that $(y,\widetilde{y})\in Y_{0}\times\widetilde{Y}$ satisfy Eq.
(SG). Since $f$ is differentiable in $y$, $v=\nabla f(y)$ is the only possible
choice for $v$ in Eq. (SG). If
$\widetilde{y}\in\operatorname{int}(\widetilde{Y})$, then $f$ is
differentiable in $\widetilde{y}$ and $\nabla f(y)=v\in\partial
f(\widetilde{y})\subset\\{\nabla f(\widetilde{y})\\}$ implies
$\widetilde{y}_{1}=\partial_{1}f(\widetilde{y})=\partial_{1}f(y)=y_{1}$, a
contradiction. Consequently we assume that
$\widetilde{y}\in\partial\widetilde{Y}$, that is,
$\widetilde{y}_{1}\in\\{y_{\mathrm{min}},y_{\mathrm{max}}\\}$. The existence
of $\partial_{2}f(\widetilde{y})$ implies
$\partial_{2}f(\widetilde{y})=v_{2}=\partial_{2}f(y)$ and equivalently

$\displaystyle\frac{b}{b+1}\widetilde{y}_{2}=y_{2}\quad\Leftrightarrow\quad\widetilde{y}_{2}-y_{2}=\frac{1}{b}y_{2}\,.$
(5.11)

The equation in Eq. (SG) is equivalent to
$f(\widetilde{y})=(T_{y}f)(\widetilde{y})$. By definition of $f$ and since
$f|_{Y_{0}}$ is a quadratic function, these conditions lead to

$\displaystyle\frac{1}{2}\widetilde{y}_{1}^{2}+\frac{1}{2}\frac{b}{b+1}\widetilde{y}_{2}^{2}=f(\widetilde{y})=(T_{y}f)(\widetilde{y})=\frac{1}{2}\widetilde{y}_{1}^{2}+\frac{1}{2}\widetilde{y}_{2}^{2}-\frac{1}{2}\nabla^{2}f(y)[\widetilde{y}-y,\widetilde{y}-y]\,.$
(5.12)

In view of $\nabla^{2}f(y)=\operatorname{Id}$, the algebraic condition can be
simplified to

$\displaystyle(\widetilde{y}_{1}-y_{1})^{2}+(\widetilde{y}_{2}-y_{2})^{2}=\frac{1}{b+1}\widetilde{y}_{2}^{2}\,.$
(5.13)

One defines $s=|\widetilde{y}_{1}-y_{1}|$ and substitutes the relation between
$\widetilde{y}_{2}$ and $y_{2}$ to obtain

$\displaystyle y_{2}=\pm\sqrt{b}s\quad\Leftrightarrow\quad
y_{2}=\pm\sqrt{b}s\,.$ (5.14)

Since $\widetilde{y}_{1}\in\\{y_{\mathrm{min}},y_{\mathrm{max}}\\}$, one
obtains four one-parameter families of pairs of points,

$\displaystyle\alpha_{\mathrm{min}}^{\pm}(s)=\left(\begin{array}[]{cc}y_{\mathrm{min}}+s\\\
\pm\sqrt{b}\cdot s\end{array}\right)$
$\displaystyle\,,\quad\beta_{\mathrm{min}}^{\pm}(s)=\left(\begin{array}[]{cc}y_{\mathrm{min}}\\\
\pm\frac{b+1}{\sqrt{b}}s\end{array}\right)\,,$ (5.19)
$\displaystyle\alpha_{\mathrm{max}}^{\pm}(s)=\left(\begin{array}[]{cc}y_{\mathrm{max}}-s\\\
\pm\sqrt{b}\cdot s\end{array}\right)$
$\displaystyle\,,\quad\beta_{\mathrm{max}}^{\pm}(s)=\left(\begin{array}[]{cc}y_{\mathrm{max}}\\\
\pm\frac{b+1}{\sqrt{b}}s\end{array}\right)\,.$ (5.24)

where $s$ can take all values in $(y_{\mathrm{min}},y_{\mathrm{max}})$ with
$\alpha_{\mathrm{min}}^{\pm}(s)\in Y_{0}$ (or
$\alpha_{\mathrm{max}}^{\pm}(s)\in Y_{0}$). For the construction of the convex
envelope only the values $0<s<(y_{\mathrm{max}}-y_{\mathrm{min}})/2=:s_{\ast}$
will be relevant, since the curves $\alpha_{\mathrm{min}}^{\pm}$ and
$\alpha_{\mathrm{max}}^{\pm}$ intersect in the point

$\displaystyle
y_{\ast}^{\pm}:=\alpha_{\mathrm{min}}^{\pm}(s_{\ast})=\alpha_{\mathrm{max}}^{\pm}(s_{\ast})=(y_{\mathrm{mid}},\pm\sqrt{b}(y_{\mathrm{max}}-y_{\mathrm{min}})/2)\in
Y_{0}\,.$ (5.25)

For each $y\in{\mathbb{R}}^{2}$, which lies on a connecting straight line of
two corresponding points of the curves $\alpha_{\mathrm{min}}^{\pm}$ and
$\beta_{\mathrm{min}}^{\pm}$, i.e.,
$y=t\cdot\alpha_{\mathrm{min}}^{\pm}(s)+(1-t)\cdot\beta_{\mathrm{min}}^{\pm}(s)$
with $s\in(0,s_{\ast})$ and $t\in[0,1]$, we define $f_{\mathrm{a}}(y)=t\cdot
f(\alpha_{\mathrm{min}}^{\pm}(s))+(1-t)\cdot f(\beta_{\mathrm{min}}^{\pm}(s))$
according to Eq. (5.8) and we use the same construction for
$\alpha_{\mathrm{max}}^{\pm}$ and $\beta_{\mathrm{max}}^{\pm}$. Note, that
$f_{\mathrm{a}}$ is affine along these connecting lines and the area, which is
covered by those lines consists of four triangular shaped sets denoted by
$Y_{2}$ in Figure 3. The fact that $\alpha_{\mathrm{min}}^{\pm}$ and
$\alpha_{\mathrm{max}}^{\pm}$ intersect in $y_{\ast}^{\pm}$ means that the
tangential plane of $f$ in $y_{\ast}^{\pm}$ touches the graph in two other
points, namely

$\displaystyle y_{\ast\mathrm{min}}^{\pm}$
$\displaystyle:=\beta_{\mathrm{min}}^{\pm}(s_{\ast})=(y_{\mathrm{min}},\pm\frac{b+1}{\sqrt{b}}(y_{\mathrm{max}}-y_{\mathrm{min}})/2)\,,$
(5.26) $\displaystyle y_{\ast\mathrm{max}}^{\pm}$
$\displaystyle:=\beta_{\mathrm{max}}^{\pm}(s_{\ast})=(y_{\mathrm{max}},\pm\frac{b+1}{\sqrt{b}}(y_{\mathrm{max}}-y_{\mathrm{min}})/2)\,.$
(5.27)

Since the pairs
$(\alpha_{\mathrm{min}}^{\pm}(s),\beta_{\mathrm{min}}^{\pm}(s))$ and
$(\alpha_{\mathrm{max}}^{\pm}(s),\beta_{\mathrm{max}}^{\pm}(s))$ satisfy the
equality in Eq. (SG) for all $s\in(0,s_{\ast})$, the triples
$(y_{\ast}^{\pm},y_{\ast\mathrm{min}}^{\pm},y_{\ast\mathrm{max}}^{\pm})$ also
satisfy Eq. (SG). Therefore we define $f_{\mathrm{a}}$ for $y$ in one of the
triangles
$\text{Conv}(\\{y_{\ast}^{+},y_{\ast\mathrm{min}}^{+},y_{\ast\mathrm{max}}^{+}\\})\cup\text{Conv}(\\{y_{\ast}^{-},y_{\ast\mathrm{min}}^{-},y_{\ast\mathrm{max}}^{-}\\})$
($Y_{3}$ in Figure 3) according to Eq. (5.9) as the unique affine function
which coincides with $f$ in the vertices of each of these triangles,

$\displaystyle
f_{\mathrm{a}}(y)=\left\\{\begin{array}[]{cl}(T_{y_{\ast}^{+}}f)(y)&\text{ on
}\text{Conv}(\\{y_{\ast}^{+},y_{\ast\mathrm{min}}^{+},y_{\ast\mathrm{max}}^{+}\\})\,,\\\
(T_{y_{\ast}^{-}}f)(y)&\text{ on
}\text{Conv}(\\{y_{\ast}^{-},y_{\ast\mathrm{min}}^{-},y_{\ast\mathrm{max}}^{-}\\})\,.\end{array}\right.$
(5.30)

Finally, consider a pair of two different points
$(\widetilde{y}^{(1)},\widetilde{y}^{(2)})\in\partial\widetilde{Y}$ satisfying
Eq. (SG). For the existence of a supporting hyperplane of $f$, which touches
the graph of $f$ in $\widetilde{y}^{(1)}$ and $\widetilde{y}^{(2)}$, we
require equal partial derivatives in the second component, i.e.,

$\displaystyle\widetilde{y}^{(1)}_{2}=\frac{b+1}{b}\partial_{2}f(\widetilde{y}^{(1)})=\frac{b+1}{b}\partial_{2}f(\widetilde{y}^{(2)})=\widetilde{y}^{(2)}_{2}\,.$
(5.31)

With
$\widetilde{y}^{(1)}_{1},\widetilde{y}^{(2)}_{1}\in\\{y_{\mathrm{min}},y_{\mathrm{max}}\\}$
and $\widetilde{y}^{(1)}\neq\widetilde{y}^{(2)}$ we get
$\widetilde{y}^{(1)}=y_{\mathrm{min}}$ and
$\widetilde{y}^{(2)}_{1}=y_{\mathrm{max}}$ and obtain for $s>0$ the pairs of
curves

$\displaystyle\alpha_{\infty}^{\pm}(s):=(y_{\mathrm{min}},\pm
s)\,,\quad\beta_{\infty}^{\pm}(s):=(y_{\mathrm{max}},\pm s)\,.$ (5.32)

Since we already constructed our candidate $f_{\mathrm{a}}$ for the convex
envelope on the set
$(y_{\mathrm{min}},y_{\mathrm{max}})\times[-\frac{b+1}{\sqrt{b}}s_{\ast},\frac{b+1}{\sqrt{b}}s_{\ast}]$,
we restrict $s$ to have values $s>\frac{b+1}{\sqrt{b}}s_{\ast}$. For each
$y\in{\mathbb{R}}^{2}$, which lies on a connecting straight line of two
corresponding points of the curves $\alpha_{\infty}^{\pm}$ and
$\beta_{\infty}^{\pm}$, i.e.,
$y=t\cdot\alpha_{\infty}^{\pm}(s)+(1-t)\cdot\beta_{\infty}^{\pm}(s)$ with
$s\in(\frac{b+1}{\sqrt{b}}s_{\ast},\infty)$ and $t\in[0,1]$ (represented by
the two unbounded sets forming $Y_{4}$ in Figure 3), we define
$f_{\mathrm{a}}(y)=t\cdot f(\alpha_{\infty}^{\pm}(s))+(1-t)\cdot
f(\beta_{\infty}^{\pm}(s))$ according to Eq. (5.8). We summarize the result of
the construction in the following theorem; the explicit formulas in the
different regions will be verified in the proof.

###### Theorem 1

Define

$\displaystyle Y_{1}$ $\displaystyle=\\{y\in{\mathbb{R}}^{2},\
|y_{2}|<r_{0}(y_{1})\\},$ (5.33) $\displaystyle Y_{2}$
$\displaystyle=\\{y\in(y_{\mathrm{min}},y_{\mathrm{max}})\times{\mathbb{R}},\
r_{0}(y_{1})\leq|y_{2}|<\frac{b+1}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}-\frac{r_{0}(y_{1})}{b}\\},$
(5.34) $\displaystyle Y_{3}$
$\displaystyle=\\{y\in(y_{\mathrm{min}},y_{\mathrm{max}})\times{\mathbb{R}},\
\frac{b+1}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}-\frac{r_{0}(y_{1})}{b}\leq|y_{2}|\leq\frac{b+1}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}\\},$
(5.35) $\displaystyle Y_{4}$
$\displaystyle=\\{y\in(y_{\mathrm{min}},y_{\mathrm{max}})\times{\mathbb{R}},\
\frac{b+1}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}<|y_{2}|\\},$
(5.36)

as shown in Figure 3. Then the convex envelope $f_{\mathrm{c}}$ of $f$ is
given by

$\displaystyle
f_{\mathrm{c}}(y)=\left\\{\begin{array}[]{cl}\displaystyle\frac{1}{2}y_{1}^{2}+\frac{1}{2}\frac{b}{b+1}y_{2}^{2}&\text{
on }\widetilde{Y}\,,\\\\[7.22743pt]
\displaystyle\frac{1}{2}y_{1}^{2}+\frac{1}{2}y_{2}^{2}&\text{ on
}Y_{1}\,,\\\\[7.22743pt]
\displaystyle\frac{1}{2}y_{1}^{2}+\frac{1}{2}y_{2}^{2}-\frac{(|y_{2}|-r_{0}(y_{1}))^{2}}{2(b+1)}&\text{
on }Y_{2}\,,\\\\[7.22743pt] \displaystyle\begin{aligned}
&\frac{1}{2}y_{1}^{2}+\frac{1}{2}y_{2}^{2}-\frac{(|y_{2}|-r_{0}(y_{1}))^{2}}{2(b+1)}\\\
\displaystyle&\qquad-\frac{1}{2}\frac{b}{b+1}\left(|y_{2}|-\frac{b+1}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}+\frac{r_{0}(y_{1})}{b}\right)^{2}\end{aligned}&\text{
on }Y_{3}\,,\\\\[7.22743pt]
\displaystyle\frac{1}{2}y_{1}^{2}+\frac{1}{2}\frac{b}{b+1}y_{2}^{2}+\frac{1}{2}(y_{1}-y_{\mathrm{min}})\cdot(y_{\mathrm{max}}-y_{1})&\text{
on }Y_{4}\,.\end{array}\right.$ (5.42)

□

###### Proof

Denote by $f_{\mathrm{a}}$ the formula on the right-hand side in the assertion
of the theorem. We are going to show in Step 1 that $f_{\mathrm{a}}$ is
differentiable on
${\mathbb{R}}^{2}\setminus(\\{y_{\mathrm{min}},y_{\mathrm{max}}\\}\times{\mathbb{R}})$
and convex. In Step 2 we verify that $f_{\mathrm{a}}$ is indeed the expression
we obtain from the previously described construction and $f_{\mathrm{c}}\leq
f_{\mathrm{a}}$. Finally in Step 3 we prove $f_{\mathrm{a}}\leq f$ and
conclude, that $f_{\mathrm{a}}$ coincides with $f_{\mathrm{c}}$ by
$f_{\mathrm{c}}\leq f_{\mathrm{a}}\leq f$ and $f_{\mathrm{a}}$ being convex.

Step 1: Differentiability and convexity of $f_{\mathrm{a}}$. The
differentiability of $f_{\mathrm{a}}$ on
$\operatorname{int}(\widetilde{Y})\cup(y_{\mathrm{min}},y_{\mathrm{max}})\times{\mathbb{R}}$
follows from

$\displaystyle\nabla
f_{\mathrm{a}}(y)=\left\\{\begin{array}[]{cl}\displaystyle(y_{1},\frac{b}{b+1}y_{2})&\text{
on }\mathring{\widetilde{Y}}\,,\\\\[7.22743pt]
\displaystyle(y_{1},y_{2})&\text{ on }Y_{1}\,,\\\\[7.22743pt]
\displaystyle(y_{1}+\frac{\operatorname{sign}(y_{\mathrm{mid}}-y_{1})\sqrt{b}}{b+1}(|y_{2}|-r_{0}(y_{1})),\frac{b}{b+1}y_{2}+\frac{\operatorname{sign}(y_{2})r_{0}(y_{1})}{b+1})&\text{
on }Y_{2}\,,\\\\[7.22743pt]
\displaystyle(y_{\mathrm{mid}},\operatorname{sign}(y_{2})\sqrt{b}\cdot\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2})&\text{
on }Y_{3}\,,\\\\[7.22743pt]
\displaystyle(y_{\mathrm{mid}},\frac{b}{b+1}y_{2})&\text{ on
}Y_{4}\,,\end{array}\right.$ (5.48)

and verifying the continuity of $\nabla f_{\mathrm{a}}$ along the boundaries
of the different regions. With the notation $e_{1}\odot
e_{2}:=\frac{1}{2}(e_{1}\otimes e_{2}+e_{2}\otimes e_{1})$, the second
derivatives are given by

$\displaystyle\nabla^{2}f_{\mathrm{a}}(y)=\left\\{\begin{array}[]{cl}\displaystyle
e_{1}\otimes e_{1}+\frac{b}{b+1}e_{2}\otimes e_{2}&\text{ on
}\operatorname{int}(\widetilde{Y})\,,\\\\[7.22743pt] \displaystyle
e_{1}\otimes e_{1}+e_{2}\otimes e_{2}&\text{ on
}\operatorname{int}(Y_{1})\,,\\\\[7.22743pt]
\displaystyle\frac{1}{b+1}\cdot\left(e_{1}\otimes
e_{1}+\frac{2\sqrt{b}\operatorname{sign}(y)}{\operatorname{sign}(y_{\mathrm{mid}}-y_{1})}e_{1}\odot
e_{2}+be_{2}\otimes e_{2}\right)&\text{ on
}\operatorname{int}(Y_{2})\,,\\\\[7.22743pt] \displaystyle 0&\text{ on
}\operatorname{int}(Y_{3})\,,\\\\[7.22743pt]
\displaystyle\frac{b}{b+1}e_{2}\otimes e_{2}&\text{ on
}\operatorname{int}(Y_{4})\,.\end{array}\right.$ (5.54)

In fact, $\nabla^{2}f_{\mathrm{a}}$ is positively semidefinite on the interior
of each domain $\widetilde{Y},Y_{1},Y_{2},Y_{3},Y_{4}$, consequently
$f_{\mathrm{a}}$ is locally convex on those domains. Since $f_{\mathrm{a}}$ is
differentiable on $(y_{\mathrm{min}},y_{\mathrm{max}})\times{\mathbb{R}}$, an
iterative application of Lemma 3 gives us the local convexity (and hence the
convexity by Lemma 6) of $f_{\mathrm{a}}$ on
$(y_{\mathrm{min}},y_{\mathrm{max}})\times{\mathbb{R}}$. A short calculation
shows, that $g:{\mathbb{R}}^{2}\to{\mathbb{R}},\
g(y)=\frac{1}{2}y_{1}^{2}+\frac{1}{2}\frac{b}{b+1}y_{2}^{2}$ is a convex lower
bound of $f_{\mathrm{a}}$, and since $f_{\mathrm{a}}$ and $g$ coincide on
$\widetilde{Y}$ and $f_{\mathrm{a}}$ is convex on
${\mathbb{R}}\setminus\widetilde{Y}=(y_{\mathrm{min}},y_{\mathrm{max}})\times{\mathbb{R}}$,
the convexity of $f_{\mathrm{a}}$ follows by Lemma 7.

Step 2: $f_{\mathrm{c}}\leq f_{\mathrm{a}}$. Since the calculation of the
explicit formula of the candidate $f_{\mathrm{a}}$ derived from the previously
described construction is extensive, we just prove, that the given formula is
the one, which is obtained by this construction. Refer to Figure 3 for a
sketch of the various domains.

On $\widetilde{Y}\cup Y_{1}$: Since by assumption $r_{0}\leq r$ and hence
$Y_{1}\subset Y_{0}$, $f_{\mathrm{a}}$ coincides with $f$ in $Y_{1}$ and in
$\widetilde{Y}$ and $f_{\mathrm{c}}\leq f=f_{\mathrm{a}}$, as asserted.

On $Y_{2}$: The set $Y_{2}$ is exactly the set of points, which are covered by
the straight line segments connecting the corresponding pairs of points on the
four pairs of curves
$(\alpha_{\mathrm{min}}^{\pm},\beta_{\mathrm{min}}^{\pm})$ and
$(\alpha_{\mathrm{max}}^{\pm},\beta_{\mathrm{max}}^{\pm})$. Assume $y\in
Y_{2}$ with $y_{1}<y_{\mathrm{mid}}$ and $y_{2}>0$ (the other cases are
analogous). There exists $s\in(0,(y_{\mathrm{max}}-y_{\mathrm{min}})/2)$ and
$t\in(0,1)$ with
$y=t\cdot\alpha_{\mathrm{min}}^{+}(s)+(1-t)\cdot\beta_{\mathrm{min}}^{+}(s)$.
To show that the given expression for $f_{\mathrm{a}}(y)$ results from the
construction described above, we notice that
$f_{\mathrm{a}}(\alpha_{\mathrm{min}}^{+}(s))=f(\alpha_{\mathrm{min}}^{+}(s))$
and
$f_{\mathrm{a}}(\beta_{\mathrm{min}}^{+}(s))=f(\beta_{\mathrm{min}}^{+}(s))$
and that $f_{\mathrm{a}}$ is affine along the connecting line between
$\alpha_{\mathrm{min}}^{+}(s)$ and $\beta_{\mathrm{min}}^{+}(s)$. This can be
seen by
$\beta_{\mathrm{min}}^{+}(s)-\alpha_{\mathrm{min}}^{+}(s)=(s-y_{\mathrm{min}})\cdot(-1,\sqrt{b}^{-1})^{T}$
and using $y_{2}$, $y_{\mathrm{mid}}-y_{1}>0$ to obtain

$\displaystyle(-1,\sqrt{b}^{-1})\nabla^{2}f_{\mathrm{a}}(y)\begin{pmatrix}-1\\\
\sqrt{b}^{-1}\end{pmatrix}=\frac{1}{b+1}(-1,\sqrt{b}^{-1})\begin{pmatrix}1&\sqrt{b}\\\
\sqrt{b}&b\end{pmatrix}\begin{pmatrix}-1\\\ \sqrt{b}^{-1}\end{pmatrix}=0\,.$
(5.55)

Since $f_{\mathrm{c}}$ must lie below the linear interpolation, we obtain
$f_{\mathrm{c}}\leq f_{\mathrm{a}}$.

On $Y_{3}$: By construction, $f_{\mathrm{a}}$ and $f$ have the same values in
$y_{\ast}^{\pm},y_{\ast\mathrm{min}}^{\pm},y_{\ast\mathrm{max}}^{\pm}$ and on
$\text{Conv}(y_{\ast}^{+},y_{\ast\mathrm{min}}^{+},y_{\ast\mathrm{max}}^{+})$
and
$\text{Conv}(y_{\ast}^{-},y_{\ast\mathrm{min}}^{-},y_{\ast\mathrm{max}}^{-})$
the function $f_{\mathrm{a}}$ is affine due to $\nabla^{2}f_{\mathrm{a}}=0$ on
$\operatorname{int}(Y_{3})$. By $y_{\ast}^{\pm}\in Y_{0}$, $f$ is
differentiable in $y_{\ast}^{\pm}$ with $\nabla f(y_{\ast}^{\pm})=\nabla
f_{\mathrm{a}}(y_{\ast}^{\pm})$. Now it follows in view of Eq. (5.30) that
$f_{\mathrm{a}}$ is the above constructed function.

On $Y_{4}$: Similar to the argument on $Y_{2}$, we recognize that for any
$y\in Y_{4}$ with $\pm y_{2}>0$ there is an
$s\in(0,\frac{b+1}{\sqrt{b}}(y_{\mathrm{max}}-y_{\mathrm{min}})/2)$ and
$t\in(0,1)$ with $y=t\cdot\alpha^{\pm}(s)+(1-t)\cdot\beta^{\pm}(s)$. We have
$f_{\mathrm{a}}(\alpha_{\infty}^{\pm}(s))=f(\alpha_{\infty}^{\pm}(s))$ and
$f_{\mathrm{a}}(\beta_{\infty}^{\pm}(s))=f(\beta_{\infty}^{\pm}(s))$, while
$f_{\mathrm{a}}$ is affine along the connecting line between
$\alpha_{\infty}^{\pm}(s)$ and $\beta_{\infty}^{\pm}(s)$ since
$\beta_{\infty}^{\pm}(s)-\alpha_{\infty}^{\pm}(s)=(y_{\mathrm{max}}-y_{\mathrm{min}},0)$
and in the formula of $f_{\mathrm{a}}$ on $Y_{4}$ the quadratic term in
$y_{1}$ cancels out. Finally, we get $f_{\mathrm{c}}\leq f_{\mathrm{a}}$ since
any value of $f_{\mathrm{a}}$ is a convex combination of two or three function
values of $f$.

Step 3: The inequality $f_{\mathrm{a}}\leq f$. We show $f_{\mathrm{a}}\leq
f^{(r_{0})}\leq f$, where the second inequality is an immediate consequence of
Lemma 1. For the first inequality, recognize, that
$f_{\mathrm{a}}=f^{(r_{0})}$ holds on $\widetilde{Y}\cup Y_{1}\cup Y_{2}$ and
the inequality only has to be shown on $Y_{3}$ and $Y_{4}$.
On $Y_{3}$ the estimate follows by
$f_{\mathrm{a}}(y)=f^{(r_{0})}(y)-\frac{1}{2}\frac{b}{b+1}\left(|y_{2}|-\frac{b+1}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}+\frac{r_{0}(y_{1})}{b}\right)^{2}$.
On $Y_{4}$ we can calculate

$\displaystyle
f^{(r_{0})}(y)-\left(\frac{1}{2}y_{1}^{2}+\frac{1}{2}\frac{b}{b+1}y_{2}^{2}\right)$
$\displaystyle=\frac{y_{2}^{2}}{2(b+1)}-\frac{(|y_{2}|-r_{0}(y_{1}))_{+}^{2}}{2(b+1)}\geq-\frac{r_{0}(y_{1})^{2}}{2(b+1)}+\frac{r_{0}(y_{1})}{b+1}|y_{2}|$
(5.56)
$\displaystyle\geq-\frac{r_{0}(y_{1})^{2}}{2b}+\frac{r_{0}(y_{1})}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}\,.$
(5.57)

For $y_{1}\in(y_{\mathrm{min}},y_{\mathrm{mid}}]$ we have

$\displaystyle-\frac{r_{0}(y_{1})^{2}}{2b}+\frac{r_{0}(y_{1})}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}$
$\displaystyle=-\frac{b(y_{1}-y_{\mathrm{min}})^{2}}{2b}+\frac{\sqrt{b}(y_{1}-y_{\mathrm{min}})}{\sqrt{b}}\cdot\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}$
(5.58)
$\displaystyle=\frac{1}{2}(y_{1}-y_{\mathrm{min}})(y_{\mathrm{max}}-y_{1})$
(5.59)

and for $y_{1}\in(y_{\mathrm{mid}},y_{\mathrm{max}})$ we have

$\displaystyle-\frac{r_{0}(y_{1})^{2}}{2b}+\frac{r_{0}(y_{1})}{\sqrt{b}}\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}$
$\displaystyle=-\frac{b(y_{\mathrm{max}}-y_{1})^{2}}{2b}+\frac{\sqrt{b}(y_{\mathrm{max}}-y_{1})}{\sqrt{b}}\cdot\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}$
(5.60)
$\displaystyle=\frac{1}{2}(y_{1}-y_{\mathrm{min}})(y_{\mathrm{max}}-y_{1})\,,$
(5.61)

which gives us the desired estimate on $Y_{4}$. ■

Figure 3: Construction of the relaxed energy

###### Corollary 1

For all $r\in\mathcal{R}$ with
$\sqrt{b}\cdot\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{2}\leq
r\Bigl{(}\frac{y_{\mathrm{min}}+y_{\mathrm{max}}}{2}\Bigr{)}$ the equality

$\displaystyle f_{\mathrm{c}}^{(r)}=f_{\mathrm{c}}^{(r_{0})}$ (5.62)

holds. □

###### Proof

In the proof of the previous theorem we already showed the inequality

$\displaystyle f_{\mathrm{c}}^{(r)}\leq f^{(r_{0})}\leq f^{(r)}\,,$ (5.63)

which implies the claimed equality by the monotonicity in taking convex
envelopes. ■

###### Remark 1

More generally, one can even show that without any restrictions on $r$ the
convex envelope of $f^{(r)}$ does not change while replacing $r$ by
$\min\\{r,r_{0}\\}$. Corollary 1 is a special case of this statement, since
Eq. (5.10) implies $\min\\{r,r_{0}\\}=r_{0}$ by concavity of $r$ on
$(y_{\mathrm{min}},y_{\mathrm{max}})$. □

##  6 Extension of the results to three dimensions

The convex envelope in the one-dimensional case for the density
$f(y_{1},y_{2},\varepsilon_{\mathrm{p},n},p_{n})$ given in Theorem 1 may be
generalized in a natural way to three dimensions by employing a substitution
inspired by the transformation of variables defined in Eq. (4.6). In this
sense, we define a formally relaxed energy as

$\psi_{\mathrm{rel}}(\boldsymbol{\varepsilon},\boldsymbol{\varepsilon}_{\mathrm{p},n},p_{n})=f\left(\sqrt{\frac{K}{2\mu}}\operatorname{tr}\boldsymbol{\varepsilon},\|\mathrm{dev}\boldsymbol{\varepsilon}\|,\boldsymbol{\varepsilon}_{\mathrm{p},n},p_{n}\right).$
(6.1)

Observe that, when applying the same substitution to the one-dimensional
condensed energy given in Eq. (4.7), we recover the original three-dimensional
condensed energy given in Eq. (3.19).

The relations of the relaxed energies according to Fig. 3, and their
derivatives with respect to the strain tensor, i.e., the corresponding relaxed
stress tensor, are given below.

Domain $\mathrm{Y_{1}}$:

$\psi_{\mathrm{1,rel}}=\frac{K}{2}(\operatorname{tr}\boldsymbol{\varepsilon})^{2}\,+\,\mu\,\|\,\mathrm{dev}\boldsymbol{\varepsilon}\|^{2},\quad\sigma_{\mathrm{1,rel}}=K(\operatorname{tr}\boldsymbol{\varepsilon})\,\mathrm{I}\,+\,2\mu\,\mathrm{dev}\boldsymbol{\varepsilon}\,.\hfill$
(6.2)

Domains $\mathrm{Y_{2}}$:

$\displaystyle\psi_{\mathrm{2,rel}}$
$\displaystyle=\frac{K}{2}(\operatorname{tr}\boldsymbol{\varepsilon})^{2}\,+\,\mu\,\|\mathrm{dev}\boldsymbol{\varepsilon}\|^{2}-\frac{4\mu^{2}}{2(2\mu+\beta)}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|-\frac{1}{2\mu}\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})\Big{)}^{2}\,,$
(6.3) $\displaystyle\begin{aligned}
\sigma_{\mathrm{2,rel}}=K(\operatorname{tr}\boldsymbol{\varepsilon})\,\mathrm{I}\,+\,2\mu\,\mathrm{dev}\boldsymbol{\varepsilon}-\,\frac{4\mu^{2}}{2\mu+\beta}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|-\frac{1}{2\mu}\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})\Big{)}\frac{\mathrm{dev}\boldsymbol{\varepsilon}}{\|\mathrm{dev}\boldsymbol{\varepsilon}\|}\\\
\qquad-\,\frac{2\mu\,\sqrt{\beta\,K}}{2\mu+\beta}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|-\frac{1}{2\mu}\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})\Big{)}\,\mathrm{I}\,.\end{aligned}$
(6.4)

Domains $\mathrm{Y_{3}}$:

$\psi_{\mathrm{3,rel}}=\frac{K}{2}(\operatorname{tr}\boldsymbol{\varepsilon})^{2}\,+\,\mu\,\|\mathrm{dev}\boldsymbol{\varepsilon}\|^{2}-\frac{4\mu^{2}}{2(2\mu+\beta)}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|-\frac{1}{2\mu}\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})\Big{)}^{2}\\\
-\frac{2\mu\beta}{2(2\mu+\beta)}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|+\frac{1}{\beta}\,\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})-\frac{2\mu+\beta}{2\mu}\,\sqrt{\frac{K}{\beta}}\,\frac{\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{max}}-\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{min}}}{2}\Big{)}^{2}\,,$
(6.5)
$\sigma_{\mathrm{3,rel}}=K(\operatorname{tr}\boldsymbol{\varepsilon})\,\mathrm{I}\,+\,2\mu\,\mathrm{dev}\boldsymbol{\varepsilon}-\,\frac{4\mu^{2}}{2\mu+\beta}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|-\frac{1}{2\mu}\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})\Big{)}\frac{\mathrm{dev}\boldsymbol{\varepsilon}}{\|\mathrm{dev}\boldsymbol{\varepsilon}\|}\\\
-\,\frac{2\mu\,\sqrt{\beta\,K}}{2\mu+\beta}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|-\frac{1}{2\mu}\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})\Big{)}\,\mathrm{I}\\\
-\frac{2\mu\sqrt{\beta\,K}}{2\mu+\beta}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|+\frac{1}{\beta}\,\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})-\frac{2\mu+\beta}{2\mu}\,\sqrt{\frac{K}{\beta}}\,\frac{\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{max}}-\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{min}}}{2}\Big{)}\,\mathrm{I}\\\
-\,\frac{2\mu\beta}{2\mu+\beta}\Big{(}\|\mathrm{dev}\boldsymbol{\varepsilon}\|+\frac{1}{\beta}\,\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})-\frac{2\mu+\beta}{2\mu}\,\sqrt{\frac{K}{\beta}}\,\frac{\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{max}}-\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{min}}}{2}\Big{)}\frac{\mathrm{dev}\boldsymbol{\varepsilon}}{\|\mathrm{dev}\boldsymbol{\varepsilon}\|}\,.$
(6.6)

Domains $\mathrm{Y_{4}}$:

$\displaystyle\psi_{\mathrm{4,rel}}=\frac{K}{2}(\operatorname{tr}\boldsymbol{\varepsilon})^{2}+\frac{2\mu\beta}{2(2\mu+\beta)}\|\mathrm{dev}\boldsymbol{\varepsilon}\|^{2}+\frac{K}{2}\Big{(}\operatorname{tr}\boldsymbol{\varepsilon}-\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{min}}\Big{)}\Big{(}\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{max}}-\operatorname{tr}\boldsymbol{\varepsilon}\Big{)}\,,\hfill$
(6.7)

$\sigma_{\mathrm{4,rel}}=K(\operatorname{tr}\boldsymbol{\varepsilon})\,\mathrm{I}+\frac{2\mu\beta}{2\mu+\beta}\mathrm{dev}\boldsymbol{\varepsilon}+K\Big{(}\frac{\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{max}}+\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{min}}}{2}-\operatorname{tr}\boldsymbol{\varepsilon}\Big{)}\,\mathrm{I}\,.\hfill$
(6.8)

for

$\rho_{0}(\operatorname{tr}\boldsymbol{\varepsilon})=\sqrt{\beta
K}\Big{(}\frac{\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{max}}-\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{min}}}{2}-|\operatorname{tr}(\boldsymbol{\varepsilon})-\frac{\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{max}}+\operatorname{tr}\boldsymbol{\varepsilon}_{\mathrm{min}}}{2}|\Big{)}\hfill$
(6.9)

##  7 Numerical results

###  7.1 The one-dimensional model problem

We make use of the dimensionless formulation introduced in Eq. (4.6). The
function $r(y_{1})$ is selected to fit the yield surface of a modified
Drucker-Prager model employing two quadratic functions over two intervals as
given below.

$\displaystyle
r(y_{1})=r_{\mathrm{max}}\left\\{\begin{array}[]{lll}\displaystyle
1-\frac{\big{(}y_{1}-y_{0}\big{)}^{2}}{\big{(}y_{0}-y_{\mathrm{min}}\big{)}^{2}}&\text{for}&y_{\mathrm{min}}\leq
y_{1}\leq y_{0}\\\ \\\ \displaystyle
1-\frac{\big{(}y_{1}-y_{0}\big{)}^{2}}{\big{(}y_{\mathrm{max}}-y_{0}\big{)}^{2}}&\text{for}&y_{0}\leq
y_{1}\leq y_{\mathrm{max}}\\\ \\\ 0&\text{otherwise}\end{array}\right.,$ (7.6)

where $y_{\mathrm{min}}=-0.058$, $y_{\mathrm{max}}=0.00107$, $y_{0}=-0.0385$
and $r_{\mathrm{max}}=0.016$, see Fig. 4. We choose a hardening parameter
$b=0.095$ to insure that we are within the regime of small values of $b$
treated in Sec. 5.

Figure 4: Elastic region given by the function $r(y_{1})$.

We model a one-dimensional specimen of length $L=1$, represented by a domain
$\Omega=(0,1)$, employing a discretization using $n=80$ standard linear finite
elements. The specimen is fixed at the left-hand side, $u(0)=v(0)=0$, and
subjected to given displacements $u(L)=u_{\mathrm{ext}}$ and
$v(L)=v_{\mathrm{ext}}$ in $x$\- and $y$-direction, respectively, at the
right-hand side. The affine solution is given by
$x\mapsto(u(x),v(x))=x(u_{\mathrm{ext}},v_{\mathrm{ext}})$ with constant
derivative $(u_{\mathrm{ext}},v_{\mathrm{ext}})$. We fix $u_{\mathrm{ext}}$
throughout our calculations and vary $v_{\mathrm{ext}}$ exploring four
different regions in the phase diagramme Fig. 3 which correspond to the four
panels in Fig. 5: $(u_{\mathrm{ext}},v_{\mathrm{ext}})$ in (a) $Y_{1}$, (b)
$Y_{2}\cap\\{y_{2}\leq r(y_{1})\\}$, (c) $Y_{2}\setminus\\{y_{2}\leq
r(y_{1})\\}$ and (d) $Y_{3}$, respectively.

The total energies are then functions of the vectors of nodal displacements
$\mathbf{u}=\\{u_{1},\ldots,u_{n+1}\\}$,
$\mathbf{v}=\\{v_{1},\ldots,v_{n+1}\\}$ with $u_{1}=v_{1}=0$ and
$u_{n+1}=u_{\mathrm{ext}}$, $v_{n+1}=v_{\mathrm{ext}}$, given as

$W_{\mathrm{cond,rel}}(\mathbf{u},\mathbf{v})=\sum_{j=1}^{n}\psi_{\mathrm{cond,rel}}\Big{(}\frac{u_{j+1}-u_{j}}{L/n},\frac{v_{j+1}-v_{j}}{L/n}\Big{)}\,\frac{L}{n},$
(7.7)

for the condensed and the relaxed energy, respectively. The minimization is
performed using the Mathematica function FindMinimum, which employs a
gradient-based local search. Note, that for the relaxed energy, the minimum is
obtained for the affine displacements $u_{i}=\frac{i-1}{n}u_{\mathrm{ext}}$,
$v_{i}=\frac{i-1}{n}v_{\mathrm{ext}}$.

We add a random perturbation of amplitude $\alpha$ to the exact solution given
above as an initial guess for the minimization procedure. Hence, let
$\mathrm{rand}(i)$ be a random function possessing a uniform distribution
within the interval $[-1,1]$. Then the initial guess is

$u^{\mathrm{ini}}_{i}=\frac{i-1}{n}u_{\mathrm{ext}}+\alpha\frac{L}{n}\mathrm{rand}(i),\quad
v^{\mathrm{ini}}_{i}=\frac{i-1}{n}v_{\mathrm{ext}}+\alpha\frac{L}{n}\mathrm{rand}(n+i),i=2,\ldots,n.$
(7.8)

Results are depicted in Fig. 5. The numerical simulations reflect the
predictions based on the phase diagram in Figure 3. In case (a) with
$v_{\mathrm{ext}}\in[0,0.0036]$, both the relaxed and the condensed energy are
convex with $f=f_{\mathrm{c}}$ and the corresponding minimizers are affine
with gradient $(u_{\mathrm{ext}},v_{\mathrm{ext}})$. For larger values of
$v_{\mathrm{ext}}\in[0.0036,0.014]$ the function $f$ is locally convex but
different from its relaxation $f_{\mathrm{c}}$ which is obtained by a mixture
of two states and which is affine on lines parallel to the boundary of the
region $Y^{4}$, see the red lines in the left panel in Figure 3. Thus the
calculation with the condensed energy is expected to discover the
microstructure supported on two points while the calculation with
$f_{\mathrm{c}}$ can use all the gradients on the affine part. In fact, the
second panel in Figure 5 shows exactly this behavior. The case (c) corresponds
to $v_{\mathrm{ext}}\in[0.014,0.063]$ and shows the same behavior of case (b),
the only difference being that $f$ is not locally convex at the point
$(u_{\mathrm{ext}},v_{\mathrm{ext}})$. Finally, for
$v_{\mathrm{ext}}\in[0.0163,0.15]$ the condensed energy $f$ is not locally
convex and the relaxed energy $f_{\mathrm{c}}$ is affine on the triangle
$Y_{3}$. Therefore the range of the gradient of the numerical solution for the
calculation with the condensed energy is expected to be located in the three
corners of the triangle while the gradient of the solution with the relaxed
energy can explore the full triangle. The last panel in Fig. 5 confirms this
prediction.

The procedure employing the condensed energy does not recover the exact
minimum, but produces values which are between 10% and 20% higher. However,
the values are lower than the energy of the condensed energy evaluated at the
affine solution. The approximate minimizer attempts to mimic the construction
of the relaxed energy by concentrating points at the corresponding positions
within the $y_{1}$-$y_{2}$-plane. The procedure employing the relaxed energy
returns the exact minimum with working precision in all cases.

Figure 5: Minimization with exact initial guess plus random perturbation for
increasing boundary values $u_{\mathrm{ext}}$, $v_{\mathrm{ext}}$. Left
column: $y_{2}=v^{\prime}(x)$ versus position $x$ of the beam, right column:
points of the approximate minimizers $(y_{1},y_{2})$ within relaxation domains
as depicted in Fig. 3. Comparison from the condensed and relaxed energies for
a mesh with 80 elements.

###  7.2 Three-dimensional model

The generalization to a three-dimensional model derived in Sec. 6 is tested
through a two-dimensional boundary value problem, namely the square plate with
a circular hole. The plate is fixed at the left-hand side and subject to
displacements in compression in $x$-direction along the right-hand side
increasing in time. Two unstructured meshes are computed, a coarse one with
1793 elements and a finer mesh with 4917 elements. The meshes are shown in
Fig. 6.

Finite element computations are performed using the finite element analysis
program (FEAP). Standard hexahedral (8-Nodes) elements are used, they are
restricted in the $3^{rd}$ direction employing plane strain conditions, which
apply throughout all examples. Stresses are calculated as analytical
derivatives of the energies with respect to the strains. Tangent operators are
calculated as second derivatives of the energies with respect to the strains,
numerically via perturbation. Equilibrium states are obtained by standard
Newton-Raphson procedure.

Figure 6: A plate with a hole is fixed on the left and subject to a given
displacement in $x$-direction on the right-hand side. The figure to the left
shows the coarse mesh with 1793 elements and that to the right the fine mesh
with 4917 elements. The elements employed for the stress plots are marked in
red.

For the computations, we use the same material parameters as in the one-
dimensional case taken from a modified Mohr-Coulomb material model with bulk
modulus $K=3.9\,\text{GPa}$, shear modulus $\mu=2.8\,\text{GPa}$, internal
friction angle $\phi=32^{\circ}$, cohesion $c=25\,\text{MPa}$, Poisson’s ratio
$\nu=0.25$, and tensile strength $\sigma_{\mathrm{t}}=5\,\text{MPa}$, see [3].
Under increasing load, the elements in the neighborhood of the hole start to
plastify first. Therefore, we compare the results for the model employing the
relaxed energy at the marked elements adjacent to the hole for both meshes.
The values are taken for a specific integration point at these elements. We
would like to stress, that a comparison with the model employing the condensed
energy is not possible, as the lack of convexity causes severe numerical
instabilities. Therefore, only the results for the relaxed model can be
inferred from Fig. 7.

Figure 7: Square plate with circular hole: from top left to bottom right:
stresses $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$ and displacement $u_{x}$
at the marked elements as a function of the external boundary displacement
$u_{x,\mathrm{ext}}$, path of average stress $\sigma_{yy}$ versus average
strain $\varepsilon_{yy}$, path of average stress $\sigma_{xy}$ versus average
strain $\varepsilon_{xy}$

The final external displacement is reached employing 100 load steps. The
corresponding strains switch from domain ${Y_{1}}$ into domain ${Y_{2}}$,
inducing an initiation of microstructure. It can be seen that the behavior of
the displacement field and the stresses are almost identical for both the
coarse and fine mesh, even after plastification. This demonstrates, that the
relaxed model is mesh independent and capable of capturing the material
behavior in different zones of relaxation.

The distribution of the stresses $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$
can be seen in Fig. 8. The relaxed model shows similar stress distribution for
the coarse and fine mesh. Microstructure is initiating at the same positions
and load step despite the different spatial discretization.

Figure 8: Contours from the plate with a hole. Distribution of the stresses
$\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$ is shown. Coarse mesh (left) and
fine mesh (right).

##  8 Conclusion and outlook

We established a variational model for pressure-dependent plasticity as common
in models of soil mechanics. Our approach employed a novel aspect by using a
characteristic function as dissipation potential. That way, it was possible to
establish a dissipation distance as well as a condensed energy connected to a
time-incremental setting in a consistent manner. This allowed us to
investigate the initiation of microstructure via relaxation theory by
calculating the quasiconvex envelope of the condensed energy. For a one-
dimensional model problem, for which the quasiconvex envelope coincides with
the convex envelope, we succeeded in arriving at a closed-form expression.
Interestingly, the quasiconvex envelope turned out to be largely independent
of the yield surface of the original model for small values of the hardening
parameter $b$. Numerical simulations confirmed the predictions of the
analytical results and illustrated the qualitatively different behaviour of
the variational models using the condensed and the relaxed energy,
respectively. In particular the simulations employing the condensed energy
succeeded in finding the necessary oscillations. It was possible to extend the
one-dimensional formulation to higher dimensions in an empirical manner.
Numerical simulations for a paradigmatic model problem provided strong
evidence for the superior predictive power of the relaxed model.

However, problems left for future investigations are manifold. The
characterization of the relaxed energy for large values of $b$ is a
challenging task since the affine curves in Figure 3 bounding the domain with
$f=f_{\mathrm{c}}$ intersect the graphs of the functions $\pm r$ and become
nonlinear. It is expected that the qualitative behaviour for the relaxation
involving at most three points in the phase diagram remains unchanged for an
intermediate range of values for $b$ and that for sufficiently large values of
$b$ two points will suffice.

From the mathematical side, the treatment of the higher-dimensional case is
completely open. In fact, our heuristic model provides a guaranteed lower
bound for the energy in the system but may, in general, fail to be the largest
lower bound and thus may not be optimal. Concerning variational modeling,
there are still problems to be solved connected to the non-associated
character of pressure-dependent plasticity. Finally, the treatment of several
successive time-increments is still open as well as establishing fully
evolutionary models.

Acknowledgment: The authors gratefully acknowledge the funding by the German
Research Foundation (DFG) within the Priority Program 2256 “Variational
Methods for Predicting Complex Phenomena in Engineering Structures and
Materials” within project 441211072/441468770 “Variational modeling of
pressure-dependent plasticity - a paradigm for model reduction via
relaxation”.

##  Appendix A Convex analysis

The construction of $f_{\mathrm{c}}$ in Theorem 1 was based on the
construction of a candidate $f_{\mathrm{a}}$ and a decomposition of
${\mathbb{R}}^{2}$ into disjoint regions with the property that
$f_{\mathrm{a}}$ was globally continuous and convex on the interior of the
different regions. This appendix collects the results that are needed in order
to prove that $f_{\mathrm{a}}$ is convex on ${\mathbb{R}}^{2}$. The proofs
follow from classical results in convex analysis, see [7] for more
information.

###### Definition 1

Suppose that $\Omega\subset{\mathbb{R}}^{d}$ is open and that
$f\colon\Omega\to{\mathbb{R}}$. Then $f$ is said to be locally convex, if for
all $x\in\Omega$ there exists an $r_{x}>0$ with $B(x,r_{x})\subset\Omega$ and
$f$ convex on $B(x,r_{x})$. A vector $v\in{\mathbb{R}}^{d}$ is called a local
subgradient for $f$ at $x\in\Omega$ if there exists an $r_{x}>0$ with
$B(x,r_{x})\subset\Omega$ and if for all $y\in B(x,r_{x})$ the inequality
$f(y)\geq f(x)+\langle v,y-x\rangle$ holds. The set of all local subgradients
is denoted by $\partial_{\mathrm{loc}}f(x)$ and referred to as the local
subdifferential. □

The subdifferential in the sense of convex analysis is denoted by $\partial
f(x)$. Then, by definition, $\partial f(x)\subset\partial_{\mathrm{loc}}f(x)$
and if $f$ is differentiable, then
$\partial_{\mathrm{loc}}f(x)\subset\\{f^{\prime}(x)\\}$; however, the local
subdifferential may be empty if $f$ the graph of $f$ does not lie above the
tangent plane in a sufficiently small neighbourhood of $x$.

The construction of the convex envelope of the condensed energy relies on a
local construction and the verification, that the local construction leads to
a convex function. The proof uses the following results.

###### Lemma 2

Fix $a$, $b\in{\mathbb{R}}$, $a<0<b$, and $f_{a}\in
C^{1}([a,0];{\mathbb{R}})$, $f_{b}\in C^{1}([0,b];{\mathbb{R}})$ convex with
$f_{a}(0)=f_{b}(0)$, $f_{a}^{\prime}(0)=f_{b}^{\prime}(0)$. Then
$f=f_{a}\chi_{[a,0)}+f_{b}\chi_{[0,b]}\in C^{1}([a,b];{\mathbb{R}})$ is
convex. □

###### Lemma 3

Suppose that $f\in C({\mathbb{R}}^{d})$, $B$, $C\subset{\mathbb{R}}^{d}$ open
and disjoint, $C$ convex and $f$ locally convex on $B\cup C$, $f\in
C^{1}(\operatorname{int}(B\cup\overline{C}))$. Then $f$ is locally convex on
$\operatorname{int}(B\cup\overline{C})$. □

###### Lemma 4

Suppose that $d\geq 2$, $\Omega\subset{\mathbb{R}}^{d}$ is open and
$f:\Omega\to{\mathbb{R}}$ is continuous. The following statements are
equivalent:

  * (i)

$f$ is locally convex in $x\in\Omega$;

  * (ii)

for $x\in\Omega$ there exists an $r_{x}>0$ with $B(x,r_{x})\subset\Omega$ such
that $\partial_{\mathrm{loc}}f(y)\neq\emptyset$ for all $y\in B(x,r_{x})$.

□

###### Lemma 5

Suppose that $a_{i}$, $b_{i}\in{\mathbb{R}}$, $i=1,2$ with
$a_{1}<b_{1}<a_{2}<b_{2}$ and that $f\in C^{0}([a_{1},b_{2}])$ with
$f|_{[a_{1},a_{2}]}$, $f|_{[b_{1},b_{2}]}$ convex. Then $f$ is convex. □

###### Lemma 6

Suppose that $\Omega\subset{\mathbb{R}}^{n}$ is open and that $f\in C(\Omega)$
is locally convex. Then for any pair $a$, $b\in\Omega$ with $\\{\lambda
a+(1-\lambda)b\colon\lambda\in[0,1]\\}\subset\Omega$, the convexity inequality
holds on $[a,b]$. In particular, if $\Omega$ is convex, then $f$ is convex on
$\Omega$. □

###### Lemma 7

Suppose $\Omega\subset{\mathbb{R}}^{d}$ is convex and $f\in
C({\mathbb{R}}^{d})$. If $f\big{|}_{\Omega}$ is convex and if there is a
convex function $g:{\mathbb{R}}^{d}\to{\mathbb{R}}$ with $g\leq f$ on
${\mathbb{R}}^{d}$ and $g=f$ on ${\mathbb{R}}^{d}\setminus\Omega$, then $f$ is
convex. □

## References

  * [1] Ball, J. M., and James, R. D. Fine phase mixtures as minimizers of energy.  Arch. Rational Mech. Anal. 100 (1987), 13–52.
  * [2] Ball, J. M., and James, R. D. Proposed experimental tests of a theory of fine microstructure and the two-well problem.  Phil. Trans. R. Soc. Lond. A 338 (1992), 389–450.
  * [3] Baotang Shen, J. S., and Barton, N. An approximate nonlinear modified mohr-coulomb shear strength criterion with critical state for intact rocks.  Journal of Rock Mechanics and Geotechnical Engineering 10 (2018), 645–652.
  * [4] Bartels, S. Linear convergence in the approximation of rank-one convex envelopes.  M2AN Math. Model. Numer. Anal. 38 (2004), 811–820.
  * [5] Bartels, S. Reliable and efficient approximation of polyconvex envelopes.  SIAM J. Numer. Anal. 43 (2005), 363–385.
  * [6] Bartels, S., Carstensen, C., Hackl, K., and Hoppe, U. Effective relaxation for microstructure simulations: algorithms and applications.  Computer Methods in Applied Mechanics and Engineering 193, 48 (2004), 5143–5175.  Advances in Computational Plasticity.
  * [7] Behr, F. Variational models in pressure dependent plasticity.  PhD thesis, Universität Regensburg, 2023 (in preparation).
  * [8] Carstensen, C., Conti, S., and Orlando, A. Mixed analytical-numerical relaxation in finite single-slip crystal plasticity.  Cont. Mech. Thermodyn. 20 (2008), 275–301.
  * [9] Carstensen, C., Hackl, K., and Mielke, A. Non-convex potentials and microstructures in finite-strain plasticity.  Proc. Roy. Soc. A 458 (2002), 299–317.
  * [10] Carstensen, C., and Plecháč, P. Numerical solution of the scalar double-well problem allowing microstructure.  Math. Comp. 66 (1997), 997–1026.
  * [11] Chipot, M., and Kinderlehrer, D. Equilibrium configurations of crystals.  Arch. Rational Mech. Anal. 103 (1988), 237–277.
  * [12] Conti, S., and Dolzmann, G. Relaxation of a model energy for the cubic to tetragonal phase transformation in two dimensions.  Math. Models. Methods App. Sci. 24 (2014), 2929–2942.
  * [13] Conti, S., and Dolzmann, G. An adaptive relaxation algorithm for multiscale problems and application to nematic elastomers.  J. Mech. Phys. Solids 113 (2018), 126–143.
  * [14] Conti, S., and Dolzmann, G. Numerical study of microstructures in single-slip finite elastoplasticity.  J. Optim. Theory Appl. 184, 1 (2020), 43–60.
  * [15] Conti, S., and Dolzmann, G. Numerical study of microstructures in multiwell problems in linear elasticity.  In Variational views in mechanics, Adv. Mech. Math. Birkhäuser/Springer, Cham, [2021] ©2021, pp. 1–29.
  * [16] Dacorogna, B. Direct methods in the calculus of variations, vol. 78 of Applied Mathematical Sciences.  Springer-Verlag, Berlin, 1989.
  * [17] DeSimone, A., and Dolzmann, G. Macroscopic response of nematic elastomers via relaxation of a class of $\rm SO(3)$-invariant energies.  Arch. Ration. Mech. Anal. 161 (2002), 181–204.
  * [18] Griewank, A., and Rabier, P. J. On the smoothness of convex envelopes.  Trans. Am. Math. Society 322, 2 (1990), 691–709.
  * [19] Hackl, K. Generalized standard media and variational principles in classical and finite strain elastoplasticity.  J. Mech. Phys. Solids 45 (1997), 667–688.
  * [20] Hackl, K., and Fischer, F. D. On the relation between the principle of maximum dissipation and inelastic evolution given by dissipation potentials.  Proceeding of the Royal Society of London A 464 (2008), 117–132.
  * [21] Hackl, K., Heinz, S., and Mielke, A. A model for the evolution of laminates in finite-strain elastoplasticity.  ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 92, 11-12 (2012), 888–909.
  * [22] Jezdan, G. Relaxation-based modeling of inelastic materials.  PhD thesis, Ruhr-Universität Bochum, 2023 (in preparation).
  * [23] Khan, M. S., and Hackl, K. Modeling of Microstructures in a Cosserat Continuum Using Relaxed Energies.  Springer International Publishing, Cham, 2018, pp. 103–125.
  * [24] Khan, M. S., and Hackl, K. Modeling of Microstructures in a Cosserat Continuum Using Relaxed Energies: Analytical and Numerical Aspects.  Springer International Publishing, Cham, 2021, pp. 57–87.
  * [25] Kochmann, D., and Hackl, K. The evolution of laminates in finite crystal plasticity: a variational approach.  Continuum Mech. Therm. 23, 1 (2011), 63–85.
  * [26] Kohn, R. V. The relaxation of a double-well energy.  Contin. Mech. Thermodyn. 3 (1991), 193–236.
  * [27] Kohn, R. V., and Strang, G. Optimal design and relaxation of variational problems. II.  Comm. Pure Appl. Math. 39 (1986), 139–182.
  * [28] Lubliner, J. Plasticity Theory.  Dover books on engineering. Dover Publications, Mineola, NY, 2008.
  * [29] Lurie, K. A., and Cherkaev, A. V. On a certain variational problem of phase equilibrium.  In Material instabilities in continuum mechanics (Edinburgh, 1985–1986), Oxford Sci. Publ. Oxford Univ. Press, New York, 1988, pp. 257–268.
  * [30] Miehe, C., and Lambrecht, M. Analysis of microstructure development in shearbands by energy relaxation of incremental stress potentials: large-strain theory for standard dissipative solids.  Internat. J. Numer. Methods Engrg. 58 (2003), 1–41.
  * [31] Miehe, C., Lambrecht, M., and Gürses, E. Analysis of material instabilities in inelastic solids by incremental energy minimization and relaxation methods: evolving deformation microstructures in finite plasticity.  J. Mech. Phys. Solids 52 (2004), 2725–2769.
  * [32] Miehe, C., Schotte, J., and Lambrecht, M. Homogenization of inelastic solid materials at finite strains based on incremental minimization principles. application to the texture analysis of polycrystals.  Journal of the Mechanics and Physics of Solids 50, 10 (2002), 2123–2167.
  * [33] Mitchell, J. K., Soga, K., et al. Fundamentals of soil behavior, vol. 3.  John Wiley & Sons, Hoboken, NJ, 2005.
  * [34] Moreau, J. La notion de sur-potentiel et les liaisons unilatérales en elastostatiques.  C. R. Acad. Sci. Paris, Série A 267 (1968), 954–957.
  * [35] Ortiz, M., and Repetto, E. Nonconvex energy minimization and dislocation structures in ductile single crystals.  Journal of the Mechanics and Physics of Solids 47, 2 (1999), 397–462.
  * [36] Pipkin, A. C. Elastic materials with two preferred states.  Quart. J. Mech. Appl. Math. 44 (1991), 1–15.
  * [37] Šilhavý, M. Rank-1 convex hulls of isotropic functions in dimension 2 by 2.  In Proceedings of Partial Differential Equations and Applications (Olomouc, 1999) (2001), vol. 126, pp. 521–529.
  * [38] Šilhavý, M. Ideally soft nematic elastomers.  Netw. Heterog. Media 2 (2007), 279–311.
  * [39] Wolf, H., König, D., and Triantafyllidis, T. Experimental investigation of shear band patterns in granular material.  J. Struct. Geol. 25 (2003), 1229 – 1240.