# A Hidden Convexity of Nonlinear Elasticity

Siddharth Singh Janusz Ginster Amit Acharya Department of Civil &
Environmental Engineering, Carnegie Mellon University, Pittsburgh, PA 15213,
email<EMAIL_ADDRESS>für Mathematik, Humboldt-Universität zu
Berlin, Unter den Linden 6, 10117 Berlin, Germany, email: ginsterj@hu-
berlin.de.Department of Civil & Environmental Engineering, and Center for
Nonlinear Analysis, Carnegie Mellon University, Pittsburgh, PA 15213, email:
<EMAIL_ADDRESS>
###### Abstract

A technique for developing convex dual variational principles for the
governing PDE of nonlinear elastostatics and elastodynamics is presented. This
allows the definition of notions of a variational dual solution and a dual
solution corresponding to the PDEs of nonlinear elasticity, even when the
latter arise as formal Euler-Lagrange equations corresponding to non-
quasiconvex elastic energy functionals whose energy minimizers do not exist.
This is demonstrated rigorously in the case of elastostatics for the Saint-
Venant Kirchhoff material (in all dimensions), where the existence of
variational dual solutions is also proven. The existence of a variational dual
solution for the incompressible neo-Hookean material in 2-d is also shown.
Stressed and unstressed elastostatic and elastodynamic solutions in 1 space
dimension corresponding to a non-convex, double-well energy are computed using
the dual methodology. In particular, we show the stability of a dual
elastodynamic equilibrium solution for which there are regions of non-
vanishing length with negative elastic stiffness, i.e. non-hyperbolic regions,
for which the corresponding primal problem is ill-posed and demonstrates an
explosive ‘Hadamard instability;’ this appears to have implications for the
modeling of physically observed softening behavior in macroscopic mechanical
response.

##  1 Introduction

The goal of this paper is to demonstrate a technique for obtaining solutions
of the governing PDE of nonlinear elasticity. The elasticity may be of Cauchy
type, with stress not necessarily arising as a derivative of an energy density
function of the deformation gradient. The methodology is of a ‘dual’
variational nature, transforming the given ‘primal’ PDE problem written in
first-order form, in terms of the physical deformation, deformation gradient,
and velocity fields, to a dual variational and PDE problem in terms of dual
Lagrange multiplier fields. The fundamental insight is to treat the primal PDE
under consideration as constraints and invoke a more-or-less arbitrarily
designable strictly convex, auxiliary potential to be optimized. Then, a dual
variational principle for the Lagrange multiplier (dual) fields can be
designed involving a dual-to-primal (DtP) mapping (i.e. an adapted change of
variables), with the special property that its Euler-Lagrange equations are
exactly the primal PDE system, interpreted as equations for the dual fields
using the DtP mapping, and this, even though the primal system may not have
the required symmetries necessary to be the Euler-Lagrange equations of any
objective functional of the primal fields alone. As we show, the dual
variational principle is convex, with its formal Euler-Lagrange system being
locally degenerate elliptic, regardless of the monotonicity properties of the
primal PDE system. Thus, to the extent that coercivity of the dual functional
can be established, existence of minimizers of the dual functional, our
variational dual solutions that are formally consistent with weak solutions of
the primal PDE system, can be established. Since these ideas do not even
depend on the existence of an energy density function for the primal problem,
the issue of existence of minimizers of a primal energy functional is not
relevant to our technique; instead it works with the structure of the stress
response function for the primal model of elasticity invoked.

The question of dual variational principles for nonlinear elasticity - in the
form of the possibility of an extremum principle of a complementary energy -
has been explored in [25, 19, 15, 9, 20]; the question of invertibility of the
first Piola-Kirchhoff stress as a function of the deformation gradient (or the
question of expressing the right Cauchy-Green tensor in terms of the first
Piola Kirchhoff stress) for physically realistic elastic response functions
has been a fundamental obstacle in these earlier studies. Our concerns are
different from these earlier works and not focused on a complementary energy
principle - as one difference, we require our dual solution to recover as its
Euler-Lagrange equations the primal system consisting of compatibility and
equilibrium (in the case of statics), not by imposing the latter as a feature
of the class of dual fields under consideration. Also, there is no analog in
the earlier works of a ‘free’ choice of an auxiliary potential that we
crucially exploit. Our work builds on ideas from [1, Sec. 6.1]-[4, 3], and
connects to the ideas of ‘hidden convexity in nonlinear PDE’ advanced by Y.
Brenier [12, 11].

An outline of the paper is as follows: in Secs. 2-4 we lay out the formal
calculations motivating and explaining the details of our approach. In Sec. 5
these ideas are made rigorous in the context of nonlinear elastostatics. Sec.
6 contains computational results demonstrating various features of the
developed methodology. Sec. 7 contains some concluding remarks.

##  2 Dual nonlinear Elastostatics

Consider the problem of nonlinear elastostatics with boundary conditions of
place and dead loads (for simplicity):

$\displaystyle\partial_{j}P_{ij}(F)+b_{i}$ $\displaystyle=0\ \mbox{in}\
\Omega$ (1a) $\displaystyle\partial_{j}y_{i}-F_{ij}$ $\displaystyle=0\
\mbox{in}\ \Omega$ (1b) $\displaystyle P_{ij}n_{j}=p_{i}\ \mbox{on}\
\partial\Omega_{p}\ $ $\displaystyle;\qquad y_{i}=y^{(b)}_{i}\ \mbox{on}\
\partial\Omega_{y},\quad\partial\Omega=\partial\Omega_{p}\cup\Omega_{y}.$

In the above, $P$ is a given tensor valued function of invertible tensors $F$
(delivering the First Piola-Kirchhoff stress tensor for a prescribed
deformation gradient), and all derivatives are w.r.t. rectangular Cartesian
coordinates on a fixed reference configuration $\Omega$ ($P$ is not assumed to
necessarily be a gradient of a scalar valued function on the space of
invertible tensors). The only restrictions on $P$ we require are that it be
sufficiently smooth in its argument and that it be of the form

$P(F)=F\,S\\!\left(F^{T}F\right)$

for $S$ being any arbitrary symmetric tensor valued function of a symmetric
tensor; this allows frame-indifference to be satisfied. Furthermore, $b$ is
the specified body force density per unit volume of the reference
configuration field and it is not required that the latter arise from a
potential.

Define the pre-dual functional by forming the scalar products of (1a) and (1b)
with the dual fields $\lambda$ and $\mu$, respectively, integrating by parts,
and imposing the prescribed b.c.s:

$\displaystyle\widehat{S}_{H}[y,F,\lambda,\mu]$
$\displaystyle=\int_{\Omega}-P_{ij}\big{|}_{F}\partial_{j}\lambda_{i}+\lambda_{i}b_{i}\,dx+\int_{\Omega}(-y_{i}\partial_{j}\mu_{ij}-F_{ij}\mu_{ij})\,dx+\int_{\Omega}H(y,F;x)\,dx$
(2) $\displaystyle\quad+\int_{\partial\Omega_{p}}p_{i}\lambda_{i}\,da\
+\int_{\partial\Omega_{y}}y^{(b)}_{i}\mu_{ij}n_{j}\,da.$

Define

$\displaystyle U:=(y,F);\qquad
D:=(\lambda,\mu);\qquad\mathcal{D}:=(\nabla\lambda,\lambda,\nabla\mu,\mu);$
(3)
$\displaystyle\mathcal{L}_{H}(U,\mathcal{D},x):=-P_{ij}\big{|}_{F}\partial_{j}\lambda_{i}+\lambda_{i}b_{i}-y_{i}\partial_{j}\mu_{ij}-F_{ij}\mu_{ij}+H(y,F;x),$

and require the choice of the potential $H$ to be such that it facilitates the
existence of a function

$U=U_{H}(\mathcal{D},x)$

which satisfies

$\frac{\partial\mathcal{L}_{H}}{\partial
U}\left(U_{H}(\mathcal{D},x),\mathcal{D},x\right)=0\qquad\forall\
(\mathcal{D},x).$ (4)

When such a ‘change of variables’ mapping,
$U_{H}=\left(y^{(H)},F^{(H)}\right)$, exists, defining the dual functional as

$\displaystyle
S_{H}[D]:=\widehat{S}_{H}\left[y^{(H)},F^{(H)},\lambda,\mu\right]=\int_{\Omega}\mathcal{L}_{H}\left(U_{H}(\mathcal{D},x),\mathcal{D},x\right)\,dx+\int_{\partial\Omega_{p}}p_{i}\lambda_{i}\,da\
+\int_{\partial\Omega_{y}}y^{(b)}_{i}\mu_{ij}n_{j}\,da$
$\displaystyle\mbox{with $\lambda$ specified on
$\partial\Omega\backslash\partial\Omega_{p}$ and $\mu$ specified on
$\partial\Omega\backslash\partial\Omega_{y}$},$

and noting (4), the first variation of $S_{H}$ (about a state $D$ in the
direction $\delta D$, the latter constrained by $\delta\lambda=0$ on
$\partial\Omega\backslash\partial\Omega_{p}$ and $\delta\mu=0$ on
$\partial\Omega\backslash\partial\Omega_{y}$), is given by

$\delta
S_{H}=\bigintsss_{\Omega}\frac{\partial\mathcal{L}_{H}}{\partial\mathcal{D}}\left(U_{H}(\mathcal{D},x),\mathcal{D},x\right)\cdot\delta\mathcal{D}\,dx+\int_{\partial\Omega_{p}}p_{i}\delta\lambda_{i}\,da\
+\int_{\partial\Omega_{y}}y^{(b)}_{i}\delta\mu_{ij}n_{j}\,da.$

Noting, now, that $\mathcal{L}_{H}$ is necessarily affine in $\mathcal{D}$,
its second argument, it can be checked that the Euler-Lagrange equations and
natural boundary conditions of the dual functional $S_{H}$ are exactly the
system (1) with $U$ substituted by $U_{H}(\mathcal{D},x)$.

It is this simple idea that we exploit to develop variational principles for
(Cauchy) Elasticity.

###  2.1 Quadratic nonlinearity about arbitrarily fixed base states

In this section we implement the generalities above in a context that allows
to demonstrate an explicitly written out dual functional. This is only for
conceptual simplicity in conveying ideas; in fact, the examples we
successfully compute in Sec. 6 will be more general than the assumptions made
here.

For solutions of (1) ‘close’ to an arbitrarily chosen (and fixed thereafter)
pair of fields we refer to as a base state given by

$x\mapsto\bar{y}(x);\qquad x\mapsto\nabla\bar{y}(x)=:\bar{F}(x),$

we assume that the nominal stress fields, $x\mapsto P(x)$, for all candidate
pairs $(y,F)$ that belong to the domain of the functional $\widehat{S}$ are
well-represented by a quadratic expansion around $(\bar{y},\bar{F})$: i.e.

$\displaystyle P_{ij}(F(x))$
$\displaystyle=\bar{P}_{ij}\left(\bar{F}(x)\right)+A_{ijkl}|_{\bar{F}(x)}\left(F_{kl}(x)-\bar{F}_{kl}(x)\right)$
(5)
$\displaystyle\quad+\frac{1}{2}B_{ijklmn}|_{\bar{F}(x)}\left(F_{kl}(x)-\bar{F}_{kl}(x)\right)\left(F_{mn}(x)-\bar{F}_{mn}(x)\right),$

where $\bar{P},A,B$ are given functions from the space of second order tensors
to the space of second, fourth, and sixth order tensors, respectively. We
assume the symmetries $B_{ijklmn}=B_{ijmnkl}$, without loss of generality. The
coefficients $\bar{P},A,B$ may be interpreted as those of a truncated Taylor
expansion of the nominal stress response function about $\bar{F}$, but not
necessarily. In the following, we will assume these functions to be bounded.

We note that the assumption on the form of the stress allows for nominal
stress fields that access ‘monotone increasing’ and ‘monotone decreasing’
parts of the ‘stress-strain’ response of the elastic material at different
spatial locations of the domain.

We now choose a ‘shifted quadratic’ form for the potential $H$:

$Q(y,F;x)=\frac{1}{2}\left(c_{y}\left|y-\bar{y}(x)\right|^{2}+c_{F}\left|F-\bar{F}(x)\right|^{2}\right),$

where

$c_{y},c_{F}\gg 1$

are real-valued scalars. Then, the implementation of (4) in this context means
algebraically solving for $\left(y^{(Q)},F^{(Q)}\right)$ in terms of
$(x,\nabla\lambda,\nabla\mu,\mu)$ from the dual-to-primal (DtP) mapping
equations

$\displaystyle\frac{\partial\mathcal{L}_{Q}}{\partial y_{i}}:$
$\displaystyle\qquad-\partial_{j}\mu_{ij}+c_{y}\left(y_{i}^{(Q)}-\bar{y}_{i}(x)\right)=0$
$\displaystyle\frac{\partial\mathcal{L}_{Q}}{\partial F_{rs}}:$
$\displaystyle\qquad-\partial_{j}\lambda_{i}\left(A_{ijrs}|_{\bar{F}(x)}+B_{ijklrs}|_{\bar{F}(x)}\left(F_{kl}^{(Q)}-\bar{F}_{kl}(x)\right)\right)-\mu_{rs}+c_{F}\left(F_{rs}^{(Q)}-\bar{F}_{rs}(x)\right)=0.$

which can be rearranged as the following relations (suppressing the explicit
$x$ dependence of the base state for brevity)

$\displaystyle c_{y}\left(y^{(Q)}_{i}-\bar{y}_{i}\right)$
$\displaystyle=\partial_{j}\mu_{ij}$ (6a) $\displaystyle
c_{F}\,\mathbb{K}_{rskl}|_{\bar{F}}\left(F_{kl}^{(Q)}-\bar{F}_{kl}\right)$
$\displaystyle=\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}|_{\bar{F}}$ (6b)
$\displaystyle\mbox{where}\qquad\mathbb{K}_{rskl}\big{|}_{(\bar{F},\nabla\lambda)}$
$\displaystyle:=\delta_{kr}\delta_{ls}-\frac{1}{c_{F}}\partial_{j}\lambda_{i}B_{ijklrs}|_{\bar{F}}.\quad\mbox{is
affine in}\ \nabla\lambda.$ (6c)

Then

$\displaystyle\widehat{S}_{Q}\left[y^{(Q)},\lambda,F^{(Q)},\mu\right]=$
$\displaystyle\quad\int_{\Omega}-\partial_{j}\lambda_{i}\biggl{(}\bar{P}_{ij}\left(\bar{F}\right)+A_{ijkl}|_{\bar{F}}\left(F_{kl}^{(Q)}-\bar{F}_{kl}\right)+\frac{1}{2}B_{ijklmn}|_{\bar{F}}\left(F_{kl}^{(Q)}-\bar{F}_{kl}\right)\left(F_{mn}^{(Q)}-\bar{F}_{mn}\right)\biggr{)}\,dx$
$\displaystyle+\int_{\Omega}\left(-\mu_{ij}\left(F_{ij}^{(Q)}-\bar{F}_{ij}\right)-\mu_{ij}\bar{F}_{ij}\right)\,dx\quad+\int_{\Omega}\left(-\left(y_{i}^{(Q)}-\bar{y}_{i}\right)\partial_{j}\mu_{ij}-\bar{y}_{i}\partial_{j}\mu_{ij}\right)\,dx$
$\displaystyle+\int_{\Omega}\frac{1}{2}\left(c_{y}\left|y^{(Q)}-\bar{y}\right|^{2}+c_{F}\left|F^{(Q)}-\bar{F}\right|^{2}\right)\,dx\quad+\int_{\partial\Omega_{p}}p_{i}\lambda_{i}\,da\quad+\int_{\partial\Omega_{y}}y^{(b)}_{i}\mu_{ij}n_{j}\,da+\int_{\Omega}\lambda_{i}b_{i}\,dx.$

For notational simplicity, in the following computations, we will omit mention
of the superscript (Q).

We obtain the dual functional through the following computation utilizing (6):

$\displaystyle S_{Q}[\lambda,\mu]$
$\displaystyle=\int_{\Omega}\left(-\bar{P}_{ij}\big{|}_{\bar{F}}\partial_{j}\lambda_{i}-\bar{F}_{ij}\mu_{ij}-\bar{y}_{i}\partial_{j}\mu_{ij}+\lambda_{i}b_{i}\right)\,dx\quad-\int_{\Omega}\Big{(}\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\Bigr{)}\left(F_{rs}-\bar{F}_{rs}\right)\,dx$
(7)
$\displaystyle\quad+\frac{1}{2}\int_{\Omega}\Bigl{(}c_{F}\delta_{kr}\delta_{ls}-\partial_{j}\lambda_{i}B_{ijklrs}\big{|}_{\bar{F}}\Bigr{)}\left(F_{kl}-\bar{F}_{kl}\right)\left(F_{rs}-\bar{F}_{rs}\right)\,dx$
$\displaystyle\quad-\int_{\Omega}\partial_{j}\mu_{ij}(y_{i}-\bar{y}_{i})\,dx+\frac{1}{2}\int_{\Omega}c_{y}(y_{i}-\bar{y}_{i})(y_{i}-\bar{y}_{i})\,dx\quad+\int_{\partial\Omega_{p}}p_{i}\lambda_{i}\,da\quad+\int_{\partial\Omega_{y}}y^{(b)}_{i}\mu_{ij}n_{j}\,da$
$\displaystyle=\int_{\Omega}\left(-\bar{P}_{ij}\big{|}_{\bar{F}}\partial_{j}\lambda_{i}-\bar{F}_{ij}\mu_{ij}-\bar{y}_{i}\partial_{j}\mu_{ij}+\lambda_{i}b_{i}\right)\,dx\quad+\int_{\partial\Omega_{p}}p_{i}\lambda_{i}\,da\quad+\int_{\partial\Omega_{y}}y^{(b)}_{i}\mu_{ij}n_{j}\,da$
$\displaystyle\quad-\frac{1}{2}\int_{\Omega}\Big{(}\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\Bigr{)}\left(F_{rs}-\bar{F}_{rs}\right)\,dx\quad-\frac{1}{2}\int_{\Omega}\partial_{j}\mu_{ij}(y_{i}-\bar{y}_{i})\,dx,$

and using (6) again, a dual functional for nonlinear (Cauchy) elastostatics is
given by

$\displaystyle S_{Q}[\lambda,\mu]=$ (8)
$\displaystyle\qquad-\frac{1}{2}\int_{\Omega}\biggl{[}\Big{(}\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\Bigr{)}\frac{1}{c_{F}}\left(\mathbb{K}\big{|}_{(\bar{F},\nabla\lambda)}^{-1}\right)_{rskl}\Big{(}\mu_{kl}+\partial_{j}\lambda_{i}A_{ijkl}\big{|}_{\bar{F}}\Bigr{)}\quad+\frac{1}{c_{y}}\partial_{j}\mu_{ij}\partial_{k}\mu_{ik}\biggr{]}\,dx$
$\displaystyle\qquad+\int_{\Omega}\left(-\bar{P}_{ij}\big{|}_{\bar{F}}\partial_{j}\lambda_{i}-\bar{F}_{ij}\mu_{ij}-\bar{y}_{i}\partial_{j}\mu_{ij}+\lambda_{i}b_{i}\right)\,dx\quad+\int_{\partial\Omega_{p}}p_{i}\lambda_{i}\,da\quad+\int_{\partial\Omega_{y}}y^{(b)}_{i}\mu_{ij}n_{j}\,da.$

###  2.2 Formal second variation and stability of dual solutions

From the DtP equations (6), a primal pair $\left(y^{(Q)},F^{(Q)}\right)$
corresponding to the dual fields

$\bar{D}:=(x\mapsto\lambda(x)=0,x\mapsto\mu(x)=0)$

is given by $(x\mapsto\bar{y}(x),x\mapsto\bar{F}(x))$.

We are now interested in obtaining the second variation of the dual functional
$S_{Q}$ about the dual state $\bar{D}$, along sufficiently smooth
perturbations $(x\mapsto\delta D(x)\,,\,x\mapsto dD(x))$ for $c_{F},c_{y}\gg
1$, as required.

We make the (physically realistic) assumption that the nominal stress response
of the material is such that $A,B$ are bounded functions on the space of
second order tensors.

Noting

$\displaystyle\mathbb{K}_{rskl}\big{|}_{(\bar{F},\nabla\lambda)}=\delta_{kr}\delta_{ls}-\frac{1}{c_{F}}\partial_{j}\lambda_{i}B_{ijrskl}\big{|}_{\bar{F}}$
$\displaystyle;\qquad\mathbb{K}_{rskl}\big{|}_{(\bar{F},\nabla\lambda=0)}=\delta_{kr}\delta_{ls}=:\mathbb{I}_{rskl}$
$\displaystyle\mathbb{K}\mathbb{K}^{-1}=\mathbb{I}$
$\displaystyle\Longrightarrow\frac{\partial\mathbb{K}^{-1}}{\partial\nabla\lambda}=-\mathbb{K}^{-1}\frac{\partial\mathbb{K}}{\partial\nabla\lambda}\mathbb{K}^{-1}$
$\displaystyle\left.\frac{\partial\mathbb{K}}{\partial\nabla\lambda}\right|_{(\bar{F},\nabla\lambda)}$
$\displaystyle=-\frac{1}{c_{F}}B\big{|}_{\bar{F}}$
$\displaystyle\mathbb{K}_{rskl}\big{|}_{(\bar{F},\nabla\lambda)}$
$\displaystyle=\mathbb{K}_{klrs}\big{|}_{(\bar{F},\nabla\lambda)},$

the first variation of $S_{Q}$, at $D$ and in the direction $\delta D$, is
given by

$\displaystyle\delta S_{Q}\big{|}_{\delta D}[D]=$
$\displaystyle\qquad-\int_{\Omega}\biggl{(}\frac{1}{c_{F}}\left(\mathbb{K}\big{|}_{(\bar{F},\nabla\lambda)}^{-1}\right)_{rskl}\left(\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\right)\left(\delta\mu_{kl}+\partial_{j}\delta\lambda_{i}A_{ijkl}\big{|}_{\bar{F}}\right)\
+\ \frac{1}{c_{y}}\partial_{j}\mu_{ij}\partial_{k}\delta\mu_{ik}\biggr{)}\,dx$
$\displaystyle\qquad+\int_{\Omega}\biggl{(}\frac{1}{2}\left(\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\right)\left(\mu_{kl}+\partial_{j}\lambda_{i}A_{ijkl}\big{|}_{\bar{F}}\right)\frac{1}{c_{F}}\biggl{(}\mathbb{K}^{-1}_{rsmn}\frac{\partial\mathbb{K}_{mnij}}{(\partial\nabla\lambda)_{pq}}\mathbb{K}^{-1}_{ijkl}\biggr{)}\bigg{|}_{(\bar{F},\nabla\lambda)}\partial_{q}\delta\lambda_{p}\biggr{)}\,dx$
$\displaystyle\qquad+\int_{\Omega}\left(-\bar{P}_{ij}\big{|}_{\bar{F}}\partial_{j}\delta\lambda_{i}-\bar{F}_{ij}\delta\mu_{ij}-\bar{y}_{i}\partial_{j}\delta\mu_{ij}+b_{i}\delta\lambda_{i}\right)\,dx\quad+\int_{\partial\Omega_{p}}p_{i}\delta\lambda_{i}\,da\quad+\int_{\partial\Omega_{y}}y^{(b)}_{i}\delta\mu_{ij}n_{j}\,da$

with the evaluation

$\delta S_{Q}\big{|}_{\delta
D}[\bar{D}]=\int_{\Omega}\left(-\bar{P}_{ij}\big{|}_{\bar{F}}\partial_{j}\delta\lambda_{i}-\bar{F}_{ij}\delta\mu_{ij}-\bar{y}_{i}\partial_{j}\delta\mu_{ij}+b_{i}\delta\lambda_{i}\right)\,dx\quad+\int_{\partial\Omega_{p}}p_{i}\delta\lambda_{i}\,da\quad+\int_{\partial\Omega_{y}}y^{(b)}_{i}\delta\mu_{ij}n_{j}\,da.$

Our interest is in examining the second variation of $S_{Q}$, at the dual
state $\bar{D}$ and in arbitrarily fixed directions $\delta D$ and $dD$ (which
we will assume to be sufficiently regular so that integrals containing mutual
products of them and their derivatives further multiplied by bounded
functions, are finite). To do so, it suffices to consider states $D$ with its
derivatives to have similar regularity (surely $\bar{D}$ belongs to this
class), and then find that

$\displaystyle d\delta S_{Q}\big{|}_{(\delta D,dD)}[D]=$ (9)
$\displaystyle\qquad-\int_{\Omega}\biggl{(}\frac{1}{c_{F}}\left(d\mu_{rs}+\partial_{j}d\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\right)\left(\mathbb{K}\big{|}_{(\bar{F},\nabla\lambda)}^{-1}\right)_{rskl}\left(\delta\mu_{kl}+\partial_{j}\delta\lambda_{i}A_{ijkl}\big{|}_{\bar{F}}\right)\
+\
\frac{1}{c_{y}}\partial_{j}d\mu_{ij}\partial_{k}\delta\mu_{ik}\biggr{)}\,dx$
$\displaystyle\qquad+\frac{1}{c_{F}^{2}}\bigl{(}\cdots\bigr{)}\quad+\quad\frac{1}{c_{F}^{3}}\bigl{(}\cdots\bigr{)},$

where the coefficients of $c_{F}^{-2},c_{F}^{-3}$ are finite (the negative
powers on $c_{F}$ arise from repeated appearances of
$\frac{\partial\mathbb{K}}{\partial\nabla\lambda}$ in the integrands
comprising these terms), so that for $c_{F}\gg 1$ the second variation is
controlled by the first term on the right hand side of (9). Furthermore, for
any dual state for which $\nabla\lambda$ is bounded on the domain, $c_{F}$ can
be chosen large enough so that $\mathbb{K}$ and $\mathbb{K}^{-1}$ are
positive-definite.

Then, for $c_{F},c_{y}\gg 1$, the second variation at such states is non-
positive and it is this feature of the dual nonlinear elastostatics problem
that we term as its hidden convexity (with a trivial change of sign).

This result is also obtained rigorously and without approximation in Sec. 5,
with a slight reinterpretation of the DtP mapping-determining condition
‘$\frac{\partial\mathcal{L}}{\partial U}=0$.’

The second variation at the state $\bar{D}$ is given, regardless of the
magnitude of $c_{F}$, by

$d\delta S_{Q}\big{|}_{(\delta
D,dD)}[\bar{D}]=-\int_{\Omega}\biggl{(}\frac{1}{c_{F}}\left(d\mu_{rs}+\partial_{j}d\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\right)\left(\delta\mu_{rs}+\partial_{j}\delta\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\right)\
+\
\frac{1}{c_{y}}\partial_{j}d\mu_{ij}\partial_{k}\delta\mu_{ik}\biggr{)}\,dx$

and is manifestly semi-definite, regardless of the definiteness properties of
the physical tensor of elastic moduli $A$.

An important corollary of the above observation is that if the primal base
state $(\bar{y},\bar{F})$ is chosen as a solution to the primal problem - even
in a situation where the physical moduli field $x\mapsto A|_{\bar{F}(x)}$
takes positive-definite and negative definite values in distinct parts of the
domain so that the primal problem (1) is unstable at this state \- a
corresponding dual solution, given explicitly by the dual state $\bar{D}$, is
neutrally stable to the class of variations considered and, in the sense of
the DtP correspondence adopted here, confers a certain stability to the primal
solution, when viewed as obtained from the dual variational principle.

##  3 Dual nonlinear Elastodynamics

An appealing feature of the dual methodology is that it seamlessly extends to
time-dependent primal problems, leading to well-defined boundary-value-
problems in space-time domains. We consider the system of elastodynamic
equations

$\displaystyle\partial_{j}P_{ij}(F)+b_{i}-\rho_{0}\,\partial_{t}v_{i}$
$\displaystyle=0\ \mbox{in}\ \Omega\times(0,T)$ (10a)
$\displaystyle\partial_{j}y_{i}-F_{ij}$ $\displaystyle=0\ \mbox{in}\
\Omega\times(0,T)$ (10b) $\displaystyle\partial_{t}y_{i}-v_{i}$
$\displaystyle=0\ \mbox{in}\ \Omega\times(0,T)$ (10c) $\displaystyle
P_{ij}n_{j}=p_{i}\ \mbox{on}\ \partial\Omega_{p}\times(0,T)\qquad$
$\displaystyle;\qquad y_{i}=y^{(b)}_{i}\ \mbox{on}\
\partial\Omega_{y}\times(0,T)$ (10d)
$\displaystyle\left.y_{i}(x,0)=y_{i}^{(0)}\right|_{x}\quad$
$\displaystyle;\quad\left.v_{i}(x,0)=v_{i}^{(0)}\right|_{x}\ \mbox{on}\
\Omega,$ (10e)

where $[0,T]$ is an interval of time, $P$ satisfies the same conditions as in
Sec. 2 and $\rho_{0}$ is a given time-independent, positive, scalar-valued
function on $\Omega$ representing the mass density in the reference
configuration.

Following the steps of Sec. 2, we define the pre-dual functional

$\displaystyle\widehat{S}_{H}[v,y,F,\lambda,\gamma,\mu]=$
$\displaystyle\qquad\int_{\Omega}\int^{T}_{0}\left(-P_{ij}\big{|}_{F}\partial_{j}\lambda_{i}+\lambda_{i}b_{i}+\rho_{0}\,v_{i}\partial_{t}\lambda_{i}\right)\,dxdt\quad+\quad\int_{\Omega}\int^{T}_{0}(-y_{i}\partial_{j}\mu_{ij}-F_{ij}\mu_{ij})\,dxdt$
$\displaystyle\quad+\int_{\Omega}\int^{T}_{0}\left(-y_{i}\partial_{t}\gamma_{i}-\gamma_{i}v_{i}\right)\,dxdt\quad+\quad\int_{\Omega}\int^{T}_{0}H(v,y,F;x,t)\,dxdt$
$\displaystyle\quad+\int_{\partial\Omega_{p}}\int^{T}_{0}p_{i}\lambda_{i}\,da\
+\int_{\partial\Omega_{y}}\int^{T}_{0}y^{(b)}_{i}\mu_{ij}n_{j}\,da+\int_{\Omega}v^{(0)}_{i}(x)\lambda_{i}(x,0)\,dx-\int_{\Omega}y^{(0)}_{i}(x)\gamma_{i}(x,0)\,dx.$

In the above, for any boundary term obtained on integration by parts, all
primal boundary conditions and initial conditions are imposed. Any integral on
a part of the boundary of the space-time region $\Omega\times(0,T)$ where no
primal ‘boundary’ data is available is simply ignored at this stage.

We define the boundary of the space-time domain as ${\cal
B}:=(\partial\Omega\times(0,T))\cup(\Omega\times\\{0,T\\})$.

Defining now

$\displaystyle U:=(v,y,F);\qquad
D:=(\lambda,\gamma,\mu);\qquad\mathcal{D}:=(\nabla\lambda,\lambda,\partial_{t}\lambda,\nabla\mu,\mu,\partial_{t}\gamma,\gamma);$
(11)
$\displaystyle\mathcal{L}_{H}(U,\mathcal{D},x,t):=-P_{ij}\big{|}_{F}\,\partial_{j}\lambda_{i}+b_{i}\lambda_{i}+\rho_{0}v_{i}\partial_{t}\lambda_{i}-y_{i}\partial_{j}\mu_{ij}-\mu_{ij}F_{ij}-y_{i}\partial_{t}\gamma_{i}-v_{i}\gamma_{i}+H(v,y,F;x,t),$

and requiring the choice of the potential $H$ to be such that it facilitates
the existence of a function

$U=U_{H}(\mathcal{D},x,t)$

which satisfies

$\frac{\partial\mathcal{L}_{H}}{\partial
U}\left(U_{H}(\mathcal{D},x,t),\mathcal{D},x,t\right)=0\qquad\forall\
(\mathcal{D},x,t),$ (12)

it is observed that defining the dual functional as

$\displaystyle S_{H}[D]$
$\displaystyle:=\widehat{S}_{H}\left[v^{(H)},y^{(H)},F^{(H)},\lambda,\mu,\gamma\right]$
(13)
$\displaystyle=\int_{\Omega}\int^{T}_{0}\mathcal{L}_{H}\left(U_{H}(\mathcal{D},x),\mathcal{D}\right)\,dxdt+\int_{\partial\Omega_{p}}\int^{T}_{0}p_{i}\lambda_{i}\,dadt$
$\displaystyle\qquad+\int_{\partial\Omega_{y}}\int^{T}_{0}y^{(b)}_{i}\mu_{ij}n_{j}\,dadt+\int_{\Omega}v^{(0)}_{i}(x)\lambda_{i}(x,0)\,dx-\int_{\Omega}y^{(0)}_{i}(x)\gamma_{i}(x,0)\,dx$
with (arbitrarily) specified $\lambda$ on ${\cal B}\backslash(\
(\partial\Omega_{p}\times(0,T))\cup(\Omega\times\\{0\\})\ )$ and
$\displaystyle\qquad\gamma\ \mbox{on}\ {\cal B}\backslash(\
(\partial\Omega_{y}\times(0,T))\cup(\Omega\times\\{0\\})\ ),$

and noting (12), the first variation of $S_{H}$ (about a state $D$ in the
direction $\delta D$, with $\delta\lambda$ and $\delta\gamma$ constrained to
vanish on parts of the boundary of ${\cal B}$ on which $\lambda$ and $\gamma$
are specified, respectively), is given by

$\displaystyle\delta S_{H}$
$\displaystyle=\bigintsss_{\Omega}\bigintsss_{0}^{T}\frac{\partial\mathcal{L}_{H}}{\partial\mathcal{D}}\left(U_{H}(\mathcal{D},x),\mathcal{D},x,t\right)\cdot\delta\mathcal{D}\,dxdt+\int_{\partial\Omega_{p}}\int^{T}_{0}p_{i}\delta\lambda_{i}\,dadt$
$\displaystyle\quad\
+\int_{\partial\Omega_{y}}\int^{T}_{0}y^{(b)}_{i}\delta\mu_{ij}n_{j}\,dadt+\int_{\Omega}v^{(0)}_{i}(x)\delta\lambda_{i}(x,0)\,dx-\int_{\Omega}y^{(0)}_{i}(x)\delta\gamma_{i}(x,0)\,dx.$

For the same reasons as in Sec. 2, it can again be checked that the Euler-
Lagrange equations and natural boundary conditions (on the space-time domain
${\cal B}$) - of the dual functional (13) are exactly the system (10) with $U$
substituted by $U_{H}(\mathcal{D},x,t)$.

As argued in [2, Sec. 7] and demonstrated in [17], the imposition of final-
time conditions on the dual fields does not in any way constrain the dual
problem from generating solutions to the primal initial-boundary-value-problem
(10) through the DtP mapping.

We now develop the explicit form of the dual action for nonlinear
elastodynamics for the quadratically nonlinear stress:

$\displaystyle P_{ij}(F(x,t))$
$\displaystyle=\bar{P}_{ij}\left(\bar{F}(x,t)\right)+A_{ijkl}|_{\bar{F}(x,t)}\left(F_{kl}(x,t)-\bar{F}_{kl}(x,t)\right)$
(14)
$\displaystyle\quad+\frac{1}{2}B_{ijklmn}|_{\bar{F}(x,t)}\left(F_{kl}(x,t)-\bar{F}_{kl}(x,t)\right)\left(F_{mn}(x,t)-\bar{F}_{mn}(x,t)\right),$

We make an arbitrary choice of a base state

$(\ (x,t)\mapsto\bar{v}(x,t)=\partial_{t}\bar{y}(x,t),\
(x,t)\mapsto\bar{y}(x,t),\ (x,t)\mapsto\bar{F}(x,t)=\nabla\bar{y}(x,t)\ )$

and adopt the following shifted quadratic for the potential $H$:

$Q(v,y,F;x,t)=\frac{1}{2}\left(c_{v}|v-\bar{v}(x,t)|^{2}+c_{y}\left|y-\bar{y}(x,t)\right|^{2}+c_{F}\left|F-\bar{F}(x,t)\right|^{2}\right).$
(15)

Next, we form the Lagrangian corresponding to the system (10) and the
potential (15):

$\displaystyle\mathcal{L}_{Q}(v,y,F,\nabla\lambda,\partial_{t}\lambda,\nabla\mu,\mu,\partial_{t}\gamma,\gamma)$
$\displaystyle=-P_{ij}\big{|}_{F}\,\partial_{j}\lambda_{i}+b_{i}\lambda_{i}+\rho_{0}v_{i}\partial_{t}\lambda_{i}-y_{i}\partial_{j}\mu_{ij}-\mu_{ij}F_{ij}-y_{i}\partial_{t}\gamma_{i}-v_{i}\gamma_{i}$
(16)
$\displaystyle\quad+\frac{1}{2}\left(c_{v}|v-\bar{v}(x,t)|^{2}+c_{y}\left|y-\bar{y}(x,t)\right|^{2}+c_{F}\left|F-\bar{F}(x,t)\right|^{2}\right).$

The DtP mapping corresponding to (16)- (14)-(15) is then given by the
relations

$\displaystyle c_{v}\left(v_{i}^{(Q)}-\bar{v}_{i}\right)$
$\displaystyle=\gamma_{i}-\rho_{0}\,\partial_{t}\lambda_{i}$ $\displaystyle
c_{y}\left(y_{i}^{(Q)}-\bar{y}_{i}\right)$
$\displaystyle=\partial_{j}\mu_{ij}+\partial_{t}\gamma_{i}$ $\displaystyle
c_{F}\,\mathbb{K}_{rskl}\big{|}_{(\bar{F},\nabla\lambda)}\left(F_{kl}^{(Q)}-\bar{F}_{kl}\right)$
$\displaystyle=\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}|_{\bar{F}},$

where $\mathbb{K}$ is defined in (6c). Utilizing these ingredients, the dual
action,

$S_{Q}[D]:=\widehat{S}_{Q}\left[v^{(Q)},y^{(Q)},F^{(Q)},\lambda,\mu,\gamma\right]$

corresponding to ‘quadratically’ nonlinear elastodynamics (10) and (15) is
given by (again dropping the superscript (Q) for brevity)

$\displaystyle S_{Q}[D]$
$\displaystyle=-\int_{\Omega}\int^{T}_{0}\bar{P}_{ij}\big{|}_{\bar{F}}\partial_{j}\lambda_{i}+\bar{v}_{i}\left(\gamma_{i}-\rho_{0}\,\partial_{t}\lambda_{i}\right)+\bar{y}_{i}\left(\partial_{j}\mu_{ij}+\partial_{t}\gamma_{i}\right)+b_{i}\lambda_{i}\,dxdt$
$\displaystyle\quad-\frac{1}{2}\int_{\Omega}\int^{T}_{0}\frac{1}{c_{v}}\left(\gamma_{i}-\rho_{0}\,\partial_{t}\lambda_{i}\right)\left(\gamma_{i}-\rho_{0}\,\partial_{t}\lambda_{i}\right)\
+\
\frac{1}{c_{y}}\left(\partial_{j}\mu_{ij}+\partial_{t}\gamma_{i}\right)\left(\partial_{j}\mu_{ij}+\partial_{t}\gamma_{i}\right)\,dxdt$
$\displaystyle\quad-\frac{1}{2}\int_{\Omega}\int^{T}_{0}\Big{(}\mu_{rs}+\partial_{j}\lambda_{i}A_{ijrs}\big{|}_{\bar{F}}\Big{)}\frac{1}{c_{F}}\left(\mathbb{K}\big{|}_{(\bar{F},\nabla\lambda)}^{-1}\right)_{rskl}\Big{(}\mu_{kl}+\partial_{j}\lambda_{i}A_{ijkl}\big{|}_{\bar{F}}\Big{)}\,dxdt$
$\displaystyle\quad+\int_{\partial\Omega_{p}}\int^{T}_{0}p_{i}\lambda_{i}\,dadt\
+\int_{\partial\Omega_{y}}\int^{T}_{0}y^{(b)}_{i}\mu_{ij}n_{j}\,dadt$
$\displaystyle\quad+\int_{\Omega}v^{(0)}_{i}(x)\lambda_{i}(x,0)\,dx-\int_{\Omega}y^{(0)}_{i}(x)\gamma_{i}(x,0),\,dx.$

##  4 Local degenerate ellipticity of dual nonlinear Elasticity - statics and
dynamics

It is easier to see the following argument for a general setup, which is what
we employ. For this section, let Greek lower-case indices belong to the set
$\\{0,1,2,3\\}$ representing Rectangular Cartesian space-time coordinates; $0$
represents the time coordinate when the PDE is time-dependent. Let upper-case
Latin indices belong to the set $\\{1,2,3,\cdots,N\\}$, indexing the
components of the $N\times 1$ array of primal variables, $U$, with, possibly,
a conversion to first-order form as necessary. Now consider the system of
primal PDE

$\partial_{\alpha}({\cal
F}_{\Gamma\alpha}(U))+G_{\Gamma}(U)=0,\qquad\Gamma=1,\ldots,N^{*}$ (17)

where upper-case Greek indices index the number of equations involved, after
conversion to first-order form when needed. For example, consider the three
conservation laws for nonlinear elastodynamics (10a). The conversion to first-
order form (so that the Lagrangian $\mathcal{L}$ (19) below contains no
derivatives in the primal variables) requires the addition of nine more primal
variables $F$ and three more as $v$ (10b), (10c). These additional twelve
relations can be written in the form

$\mathcal{A}_{\Gamma I\alpha}\partial_{\alpha}U_{I}-\mathcal{B}_{\Gamma
I}U_{I}=0,\qquad\Gamma=4,\ldots 15,$ (18)

where $\mathcal{A},\mathcal{B}$ are constant matrices (with $\mathcal{B}$
diagonal in many cases) that define the augmentation of the primal list from
$(y)$ to $(y,F,v)$, and define the augmenting primal variables as, in general,
linear combinations of the partial derivatives of components of $U$. The
equation set (18) can be expressed in the form (17).

We assume that the functions $U\mapsto\frac{\partial^{2}{\cal
F}_{\Gamma}}{\partial U_{P}\partial U_{R}}(U)$ and
$U\mapsto\frac{\partial^{2}G_{\Gamma}}{\partial U_{P}\partial U_{R}}(U)$ are
bounded functions on their domains.

Let $D$ be the $N^{*}\times 1$ array of dual fields and, as earlier, let us
consider a shifted quadratic for the potential $H$, characterized by a
diagonal matrix $[a_{kj}]$ with constant positive diagonal entries so that the
Lagrangian takes the form

$\mathcal{L}(U,D,\nabla D,\bar{U}):=-{\cal
F}_{\Gamma\alpha}(U)\partial_{\alpha}D_{\Gamma}+D_{\Gamma}G_{\Gamma}(U)+\frac{1}{2}(U_{k}-\bar{U}_{k})a_{kj}(U_{j}-\bar{U}_{j}).$
(19)

Then the corresponding DtP mapping, obtained by ‘solving
$\frac{\partial\mathcal{L}}{\partial U}=0$ for $U$ in terms of $(\nabla
D,D,\bar{U})$,’ is given by the implicit equation

$U_{J}^{(Q)}(\nabla
D,D,\bar{U})=\bar{U}_{J}+(a^{-1})_{JK}\left(\frac{\partial{\cal
F}_{\Gamma\alpha}}{\partial U_{K}}\bigg{|}_{U^{(Q)}(\nabla
D,D,\bar{U})}\partial_{\alpha}D_{\Gamma}-D_{\Gamma}\frac{\partial
G_{\Gamma}}{\partial U_{K}}\bigg{|}_{U^{(Q)}(\nabla D,D,\bar{U})}\right).$
(20)

By the fundamental property of the dual scheme, the dual E-L equation is then
given by

$\partial_{\alpha}\Big{(}{\cal F}_{\Gamma\alpha}\big{(}U(\nabla
D,D,\bar{U}))\Big{)}+G_{\Gamma}(U(\nabla D,D,\bar{U}))=0$ (21)

(where we have dropped the superscript (Q) for notational convenience), whose
ellipticity is governed by the term

$\mathbb{A}_{\Gamma\alpha\Pi\mu}(\nabla D,D,\bar{U}):=\frac{\partial{\cal
F}_{\Gamma\alpha}}{\partial U_{P}}\bigg{|}_{U(\nabla
D,D,\bar{U})}\,\frac{\partial U_{P}}{\partial(\nabla
D)_{\Pi\mu}}\bigg{|}_{U(\nabla D,D,\bar{U})}.$

From (20) we have

$\displaystyle
a^{-1}_{PR}\left(\delta_{\Gamma\Pi}\delta_{\mu\alpha}\frac{\partial{\cal
F}_{\Gamma\alpha}}{\partial
U_{R}}\,+\,\partial_{\alpha}D_{\Gamma}\frac{\partial^{2}{\cal
F}_{\Gamma\alpha}}{\partial U_{R}\partial U_{S}}\frac{\partial
U_{S}}{\partial(\nabla
D)_{\Pi\mu}}-D_{\Gamma}\frac{\partial^{2}G_{\Gamma}}{\partial U_{R}\partial
U_{S}}\frac{\partial U_{S}}{\partial(\nabla D)_{\Pi\mu}}\right)$
$\displaystyle=\frac{\partial U_{P}}{\partial(\nabla D)_{\Pi\mu}}$
$\displaystyle\Longrightarrow\left(\delta_{PS}-a^{-1}_{PR}\partial_{\alpha}D_{\Gamma}\frac{\partial^{2}{\cal
F}_{\Gamma\alpha}}{\partial U_{R}\partial
U_{S}}+a^{-1}_{PR}D_{\Gamma}\frac{\partial^{2}G_{\Gamma}}{\partial
U_{R}\partial U_{S}}\right)\frac{\partial U_{S}}{\partial(\nabla D)_{\Pi\mu}}$
$\displaystyle=a^{-1}_{PR}\frac{\partial{\cal F}_{\Pi\mu}}{\partial U_{R}},$

and so

$\mathbb{A}_{\Gamma\alpha\Pi\mu}(0,0,\bar{U})=\frac{\partial{\cal
F}_{\Gamma\alpha}}{\partial
U_{P}}\bigg{|}_{\bar{U}}\,a^{-1}_{PR}\,\frac{\partial{\cal
F}_{\Pi\mu}}{\partial U_{R}}\bigg{|}_{\bar{U}},$

which is positive semi-definite on the space of $N^{*}\times 3$ (or
$N^{*}\times 4$) matrices. This establishes the degenerate ellipticity of the
dual system at the state $x\mapsto D(x)=0$.

To examine the ellipticity-related properties of the system in a bounded
neighborhood, say $\mathcal{N}$, of $(D=0,\nabla
D=0)\in\mathbb{R}^{N^{*}}\times\mathbb{R}^{N^{*}\times\bar{\alpha}},\bar{\alpha}=3,4$,
we define

$M_{PS}:=\delta_{PS}-a^{-1}_{PR}\left(\partial_{\alpha}D_{\Gamma}\frac{\partial^{2}{\cal
F}_{\Gamma\alpha}}{\partial U_{R}\partial
U_{S}}+D_{\Gamma}\frac{\partial^{2}G_{\Gamma}}{\partial U_{R}\partial
U_{S}}\right),$

and note that

$\frac{\partial U_{P}}{\partial(\nabla
D)_{\Pi\mu}}=M^{-1}_{PQ}a^{-1}_{QR}\frac{\partial{\cal F}_{\Pi\mu}}{\partial
U_{R}},$

where $M^{-1}$ exists and is positive definite by the boundedness of
$\mathcal{N}$ and the second derivatives of the functions ${\cal F}$ and $G$,
along with an appropriately large choice of the elements of the diagonal
matrix $[a_{ij}]$.

The degenerate ellipticity or ‘convexity’ of the system (17) in the
neighborhood $\mathcal{N}$ is now defined as the positive semi-definiteness of
the matrix $\mathbb{A}$ on the space $\mathbb{R}^{N^{*}\times\bar{\alpha}}$ of
matrices, and this in turn is governed by the matrix

$\mathbb{A}^{(sym)}_{\Gamma\alpha\Pi\mu}\bigg{|}_{(\nabla
D,D,\bar{U})}=\frac{\partial{\cal F}_{\Gamma\alpha}}{\partial
U_{P}}\bigg{|}_{U(\nabla
D,D,\bar{U})}\frac{1}{2}\left(M^{-1}_{PQ}\bigg{|}_{U(\nabla
D,D,\bar{U})}a^{-1}_{QR}+M^{-1}_{RQ}\bigg{|}_{U(\nabla
D,D,\bar{U})}a^{-1}_{QP}\right)\frac{\partial{\cal F}_{\Pi\mu}}{\partial
U_{R}}\bigg{|}_{U(\nabla D,D,\bar{U})}.$

By the positive definiteness of the matrix $[M_{PS}]$ in the neighborhood
$\mathcal{N}$, it follows that

$\xi^{\Gamma\alpha}\ \mathbb{A}^{(sym)}_{\Gamma\alpha\Pi\mu}\bigg{|}_{(\nabla
D,D,\bar{U})}\xi^{\Pi\mu}\geq 0\qquad\forall\qquad(D,\nabla
D)\in\mathcal{N},\quad\xi\in\mathbb{R}^{N^{*}\times\bar{\alpha}}$

which establishes a ‘local’ degenerate ellipticity of the system (17). We note
that degenerate ellipticity is stronger than the Legendre-Hadamard condition
given by the requirement of positive semi-definiteness of $\mathbb{A}$ on the
space of tensor products from
$\mathbb{R}^{N^{*}}\otimes\mathbb{R}^{\bar{\alpha}}$, and not directly
comparable to the strong-ellipticity condition, since it is weaker than the
latter when restricted to the space
$\mathbb{R}^{N^{*}}\otimes\mathbb{R}^{\bar{\alpha}}$ but simultaneously
requiring semi-definiteness on the larger space of
$\mathbb{R}^{N^{*}\times\bar{\alpha}}$. Also of note is that degenerate
ellipticity does not preclude the failure of ellipticity characterized by the
condition $\det[\mathbb{A}_{\Gamma\alpha\Pi\mu}n_{\alpha}n_{\mu}]\neq 0$ for
all unit direction $n\in\mathbb{R}^{\bar{\alpha}}$, $\bar{\alpha}=3$ or $4$,
thus allowing for weak (gradient) discontinuities of weak solutions
$x\mapsto\dot{D}(x)$ of the linearization of (21).

If a solution of the primal system is close to the base state $\bar{U}$, then
it seems natural to expect, due to this local degenerate ellipticity, that
such a solution can be obtained in a ‘stable’ manner by the dual formulation
designed by the choice of the auxiliary potential $H$ as a shifted quadratic
about the base state $\bar{U}$, for instance by an iterative scheme starting
from a guess $(D=0,U=\bar{U})$.

Our experience (Sec. 6 of this work, [17, 16, 5]) shows that this observation
is of great practical relevance in using the dual scheme, and we consistently
exploit it in all our computational approximations.

##  5 Variational Analysis of dual nonlinear Elastostatics

In the remaining sections of the paper we will assume that the body force
$b=0$.

In this section, we will follow the approach outlined in Section 2 and study
the corresponding dual functional $S_{H}$. First, note that under sufficient
growth and regularity assumptions on $H$ the relation (4) is satisfied if
$U_{H}(\mathcal{D},x)$ maximizes $-\mathcal{L}_{H}(\cdot,\mathcal{D},x)$ for
fixed values of $(\mathcal{D},x)$. If for all $D$ the maximizer is unique then
we call the mapping $(\mathcal{D}(x),x)\mapsto U_{H}(\mathcal{D}(x),x)$ the
dual-to-primal mapping. With a little abuse of notation we also use the same
notation $U_{H}$ for the mapping between functions defined for a pair of
functions $D=(\lambda,\mu)$ as $D\to(x\to U_{H}(\mathcal{D}(x),x))$. In this
setting the dual functional can equivalently be written as (c.f. Proposition
5.1, (29) and (31) for the precise statement for the examples considered
below)

$\tilde{S}_{H}[D]:=\sup_{y,F}-\widehat{S}_{H}[y,F,\lambda,\mu]=-S_{H}[D],$
(22)

but the definition of $\tilde{S}_{H}$ makes sense more generally even in the
absence of a well-defined mapping $U_{H}$, which suffices for the notion of a
variational dual solution defined below. It is this definition of
$\tilde{S}_{H}$ that we adopt in this section. We note that the functional
$\hat{S}_{H}$ is affine in the functions $\lambda$, $\nabla\lambda$,
$\operatorname{div}(\mu)$ and $\mu$. In particular, if formulated in
appropriate spaces it will be continuous with respect to weak convergence of
these quantities. This implies that $\tilde{S}_{H}$ as the $\sup$ over these
weakly continuous functions is weakly lower semicontinuous. This renders the
dual problem $\tilde{S}_{H}$ to be rather well-suited for applying the direct
method of the Calculus of Variations to prove the existence of a minimizer. If
the dual functional $\tilde{S}_{H}$ is regular enough, such a minimizer would
then satisfy the Euler-Lagrange equations of $S_{H}$ in a weak form (which are
the same as for $\tilde{S}_{H}$) which correspond to the equations (1) if
there is a well-defined mapping from the dual states to the primal states as
explained above. The direct method of the Calculus of Variations essentially
corresponds to performing the following three steps:

  1. 1.

Show that $\inf_{\lambda,\mu}\tilde{S}_{H}>-\infty$ and take a minimizing
sequence $(\lambda_{k},\mu_{k})$ such that
$\tilde{S}_{H}(\lambda_{k},\mu_{k})\to\inf_{\lambda,\mu}\tilde{S}$.

  2. 2.

Prove up to a non-relabeled subsequence (weak) compactness for the minimizing
sequence $(\lambda_{k},\mu_{k})\rightharpoonup(\lambda,\mu)$.

  3. 3.

Show that $S_{H}$ is sequentially lower-semicontinuous with respect to the
convergence established in 2. to conclude that
$\tilde{S}_{H}(\lambda,\mu)\leq\liminf_{k\to\infty}\tilde{S}_{H}(\lambda_{k},\mu_{k})=\inf_{\lambda,\mu}\tilde{S}_{H}$
i.e., $(\lambda,\mu)$ minimizes $\tilde{S}_{H}$.

Note that by the argumentation above the dual functional $\tilde{S}_{H}$ is
(at least formally) lower semi-continuous with respect to weak convergence.
Consequently, the main work in this setting is to establish step 2.

Let us also note that as the function $\hat{S}_{H}$ is affine in $\lambda$ and
$\mu$ the function $\tilde{S}_{H}$ as the supremum over affine functions is
automatically convex (however, most likely not strictly convex). This means
that there is a large toolbox for the computational task to find a minimizer
of the dual functional $\tilde{S}_{H}$.

In light of the discussion above, we define the notion of a variational dual
solution to the equation (1) and a dual solution.

######  Definition 5.1.

Given a Borel-measurable function $H$ we say that a pair $(\lambda,\mu)$ is a
variational dual solution to the equation (1) if it minimizes the dual
functional $\tilde{S}_{H}$. Moreover, we say that a pair $(\lambda,\mu)$ is a
dual solution to (1) if the dual functional $\tilde{S}_{H}$ is differentiable
at $(\lambda,\mu)$, the pair $(\lambda,\mu)$ satisfies the weak form of the
corresponding Euler-Lagrange equation for $\tilde{S}_{H}$ and there exists a
smooth dual-to-primal mapping $U_{H}$ corresponding to $\tilde{S}_{H}$.

By the considerations above, a dual solution gives rise to a weak solution of
(1). Note that for the notion of a variational dual solution in principle we
do not require $H$ to satisfy any additional assumption. However, the notion
of a dual variational solution is only meaningful if the function $H$ is
chosen in a way such that the dual functional $\tilde{S}_{H}$ allows for
nontrivial minimizers. For example, for $H=0$ it is easy to check that for the
examples discussed below it follows $\tilde{S}_{H}[D]=+\infty$ except for
$D=(0,0)$, where $\tilde{S}_{H}[D]=0$. In particular, the dual functional is
minimized at $(0,0)$ but every pair of primal variables is a maximizer for the
middle term of (22). This means that not every variational dual solution
corresponds to a dual solution in the sense of Definition 5.1 or gives rise to
a weak solution of (1). Consequently, ideally we would like to find a function
$H$ such that the dual functional $\tilde{S}_{H}$ is regular and the dual-to-
primal mapping is well-defined and regular. In order to obtain these two
properties, it appears to be reasonable (but in most cases not sufficient) to
choose $H$ such that $H$ has a strictly positive definite Hessian at all
points (in particular $H$ is convex) and the growth of $H$ at $+\infty$ is not
smaller than the growth of the nonlinearity $P(F)$ so that a maximizer for the
problem in the middle of (22) exists and might be unique (at least for a large
number of dual configurations).

In the case of a hyperelastic materials the typical variational strategy to
obtain a weak solution to (1) would be to prove existence of a minimizer for
its free energy functional subject to $y=y^{(b)}$ on $\Omega_{y}$. Much
research in this direction was inspired by the seminal works on functionals
with quasiconvex [23] or polyconvex energy densities [6, 24]. However, these
notions of convexity do not hold for numerous classical models, e.g. the Saint
Venant - Kirchhoff model discussed in the section below. In addition, even if
existence of a minimizer for the free energy can be established, it is not
clear that the minimizer satisfies the weak formulation of the corresponding
Euler-Lagrange equation, see the discussion in [7] and the references therein
(in the case of a convex integrand see, for example, [18, 13] for a positive
result).

Hence, we propose the concept of a variational dual solution as another
variational solution concept for nonlinear elastostatics given by (1).

In the following we will show that the concept of variational dual solutions
is applicable in a meaningful way to some standard models, namely, we will
prove, for auxiliary potentials $H$ that render the dual functional
$\tilde{S}_{H}$ non-trivial, existence of variational dual solutions to (1)
for the Saint Venant-Kirchhoff model (in all dimensions) and the
incompressible Neo-Hookean model (in two dimensions). We will focus on proving
coercivity for the dual functionals $\tilde{S}_{H}$ which guarantees that step
2. in the direct method can be performed. In these examples the formal
argument for the lower semi-continuity of the dual functional $\tilde{S}_{H}$
can also be made rigorous and hence existence of minimizers of $\tilde{S}_{H}$
can be established rigorously.

Following the arguments in the beginning of Section 3 one can set up
analogously a dual variational problem in the elastodynamic setting. The dual
variational problem is again (formally) lower semi-continuous and convex.
Understanding the coercivity of the dual variational problem in the dynamical
setting, which would in turn give rise to the existence of the analogue of
dual variational solutions, is currently work in progress.

In the following we will assume that $\Omega\subseteq\mathbb{R}^{d}$ is open,
connected and bounded with Lipschitz boundary.

###  5.1 The Saint Venant-Kirchhoff model

We consider the Saint Venant-Kirchhoff energy

$W(F)=\mathsf{G}|E|^{2}+\frac{\mathsf{L}}{2}\operatorname{tr}(E)^{2},$

where $C=F^{T}F$ is the right Cauchy-Green tensor, $E=\frac{1}{2}(C-Id)$ is
the Green-Lagrange strain tensor, the shear modulus $\mathsf{G}>0$ and
$\mathsf{G}+\frac{d}{2}\mathsf{L}>0$. Correspondingly, we set

$P(F)=DW(F)=F\left(\mathsf{G}\,F^{T}F+\frac{\mathsf{L}}{2}|F|^{2}Id-(\mathsf{G}+\frac{\mathsf{L}d}{2})Id\right).$

For $y^{(b)}\in W^{1,4}(\Omega;\mathbb{R}^{d})$ and $p^{(b)}\in
L^{4/3}(\partial\Omega;\mathbb{R}^{d})$ we look for weak solutions to the
equations (1) which take the form

$\begin{cases}&F=\nabla y\text{ in }\Omega,\\\ &\operatorname{div}P(F)=0\text{
in }\Omega,\\\ &P(F)\cdot n=p^{(b)}\text{ on }\partial\Omega_{p},\\\
&y=y^{(b)}\text{ on }\partial\Omega_{y}.\end{cases}$ (23)

Precisely, we look for $y\in W^{1,4}(\Omega;\mathbb{R}^{d})$ and $F\in
L^{4}(\Omega;\mathbb{R}^{d\times d})$ such that for all $\lambda\in
W^{1,4}(\Omega;\mathbb{R}^{d})$ with $\lambda=0$ on $\partial\Omega_{y}$ and
$\mu\in L^{4/3}(\Omega;\mathbb{R}^{d\times d})$ with
$\operatorname{div}(\mu)\in L^{4/3}(\Omega;\mathbb{R}^{d})$ and $\mu n=0$ on
$\partial\Omega_{p}$ that

$\int_{\Omega}\operatorname{div}(\mu)\cdot
y+\mu:F+P(F):\nabla\lambda\,dx=\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{d-1}+\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{d-1}.$

Note that the boundary terms are well-defined but have to be understood in the
sense of traces i.e., for $p\in[1,\infty]$ there are continuous trace
operators $W^{1,p}(\Omega;\mathbb{R}^{d})\to
W^{1-1/p,p}(\partial\Omega;\mathbb{R}^{d})\subseteq L^{p}(\partial\Omega)$ and
$\\{\mu\in L^{p}(\Omega;\mathbb{R}^{d}):\operatorname{div}(\mu)\in
L^{p}\\}\to\left(W^{1-1/p}(\partial\Omega;\mathbb{R}^{d})\right)^{\prime}$,
$\mu\to\mu n$, see [22, Theorem 18.40] and [14, Theorem 1, page 204] for the
case $p=2$ (the general case follows with the same proof). In particular, the
boundary terms are continuous with respect to weak convergence.

Next, let us mention that the stored energy density $W$ is not rank-$1$ convex
and hence neither polyconvex nor quasiconvex (see [21]). In particular, the
corresponding stored elastic energy is not weakly lower semi-continuous and
existence for (23) via the minimization of the corresponding energy cannot be
obtained by the direct method of the Calculus of Variations, c.f. the
discussion in the preamble of this section.

Following the approach explained above we define for

$H(y,F)=\frac{1}{4}|y|^{4}+\frac{\mathsf{G}+\mathsf{L}/2}{2}|F^{T}F|^{2}$

the function $g:\mathbb{R}^{d\times
d}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\to\mathbb{R}$ given by

$g(A,a,B)=\sup_{F\in\mathbb{R}^{d\times d},y\in\mathbb{R}^{d}}a\cdot
y+A:F+B:P(F)-H(y,F)$

and consider the functional $I:\\{(\mu,\lambda)\in
L^{4/3}(\Omega;\mathbb{R}^{d\times d})\times
W^{1,4}(\Omega;\mathbb{R}^{d}):\operatorname{div}(\mu)\in
L^{4/3}(\Omega;\mathbb{R}^{d}),\mu n=0\text{ on }\partial\Omega_{p}\text{ and
}\lambda=0\text{ on }\partial\Omega_{y}\\}\to\mathbb{R}$ defined as

$I(\mu,\lambda)=\int_{\Omega}g(\mu,\operatorname{div}(\mu),\nabla\lambda)\,dx-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{d-1}-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{d-1}.$ (24)

We will first show that $I$ agrees with the dual functional $\tilde{S}_{H}$.

######  Proposition 5.1.

It holds for all $\mu\in L^{4/3}(\Omega;\mathbb{R}^{d\times d})$, $\lambda\in
W^{1,4}(\Omega;\mathbb{R}^{d\times d})$ that

$\tilde{S}_{H}[\lambda,\mu]=\int_{\Omega}g(\mu,\operatorname{div}(\mu),\nabla\lambda)\,dx-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{d-1}-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{d-1}.$

###### Proof.

It suffices to show that for $\mu\in L^{4/3}(\Omega;\mathbb{R}^{d\times d})$,
$\lambda\in W^{1,4}(\Omega;\mathbb{R}^{d\times d})$ that it holds

$\sup_{y\in L^{4}(\Omega;\mathbb{R}^{d}),F\in L^{4}(\Omega;\mathbb{R}^{d\times
d})}\int_{\Omega}\operatorname{div}(\mu)\cdot
y+\mu:F+P(F):\nabla\lambda-H(y,F)\,dx=\int_{\Omega}g(\mu,\operatorname{div}(\mu),\nabla\lambda)\,dx.$
(25)

Clearly, it holds that

$\sup_{y\in L^{4}(\Omega;\mathbb{R}^{d}),F\in L^{4}(\Omega;\mathbb{R}^{d\times
d})}\int_{\Omega}\operatorname{div}(\mu)\cdot
y+\mu:F+P(F):\nabla\lambda-H(y,F)\,dx\leq\int_{\Omega}g(\mu,\operatorname{div}(\mu),\nabla\lambda)\,dx.$

We will now prove the reverse inequality. Let
$\mathcal{O}=\\{O\subseteq\Omega:O\text{ open}\\}$ and consider the countable
family of functions

$\displaystyle\mathcal{A}=\\{(y,F):$ $\displaystyle
F=\sum_{n=1}^{N}q_{n}\chi_{E_{n}}\text{ and
}y=\sum_{n=1}^{N}q_{n}^{\prime}\chi_{E_{n}},\text{ where
}q_{n}\in\mathbb{Q}^{d\times d},\,q_{n}^{\prime}\in\mathbb{Q}^{d}\text{ and }$
(26) $\displaystyle E_{n}\text{ is a cuboid with rational side lengths and
corners in }\mathbb{Q}^{d}\\}.$ (27)

Note that $\mathcal{A}$ is countable. Moreover, it holds for all
$(a,A,B)\in\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}^{d\times
d}$ and $x\in\Omega$

$\sup_{y,F\in\mathcal{D}}a\cdot y(x)+A:F(x)+P(F(x)):B-H(y(x),F(x))=g(A,a,B).$

For $O\in\mathcal{O}$ let us now define $\mathcal{F}(O)=\sup_{y\in L^{4},F\in
L^{4}}\int_{O}y\cdot\operatorname{div}(\mu)+\mu:F+P(F):\nabla\lambda-H(y,F)\,dx$.
It can be shown that for $O,O^{\prime}\in\mathcal{O}$ such that
$\operatorname{dist}(O,O^{\prime})>0$ it follows that

$\mathcal{F}(O\cup O^{\prime})=\mathcal{F}(O)+\mathcal{F}(O^{\prime}).$

Moreover, it holds by definition for each $(y,F)\in\mathcal{A}$ and
$O\in\mathcal{O}$ that

$\mathcal{F}(O)\geq\int_{O}y\cdot\operatorname{div}(\mu)+\mu:F+P(F):\nabla\lambda-H(y,F)\,dx.$

Then it follows from [10, Lemma 15.2] that

$\mathcal{F}(\Omega)\geq\int_{\Omega}\sup_{(y,F)\in\mathcal{A}}\left[y\cdot\operatorname{div}(\mu)+\mu:F+P(F):\nabla\lambda-H(y,F)\right]\,dx=\int_{\Omega}g(\mu,\operatorname{div}(\mu),\nabla\lambda)\,dx.$

This proves the reverse inequality. ∎

As explained above, in order to prove existence of minimizers for
$\tilde{S}_{H}$ it suffices by the direct method of the Calculus of Variations
to prove that $\tilde{S}_{H}$ is coercive. Hence, we will now estimate the
function $g$ from below.

######  Proposition 5.2.

There exists $c>0$ such that it holds for all
$(a,A,B)\in\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\times\mathbb{R}^{d\times
d}$

$g(A,a,B)\geq c\left(|a|^{4/3}+|A|^{4/3}+|B|^{4}\right)-\frac{1}{c}.$

###### Proof.

By optimizing in $y$ it follows that

$\displaystyle g(A,a,B)$ $\displaystyle=$
$\displaystyle\frac{3}{4}|a|^{4/3}+\sup_{F}\left[\left(A-(\mathsf{G}+\frac{\mathsf{L}d}{2})B\right):F+B:\left(\mathsf{G}\,FF^{T}F+\frac{\mathsf{L}}{2}|F|^{2}F\right)-\frac{\mathsf{G}+\mathsf{L}/2}{2}|F^{T}F|^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{3}{4}|a|^{4/3}+2(\mathsf{G}+\mathsf{L}/2)\sup_{F}\left[\tilde{A}:F+\frac{1}{2(\mathsf{G}+\mathsf{L}/2)}B:\left(\mathsf{G}\,FF^{T}F+\frac{\mathsf{L}}{2}|F|^{2}F\right)-\frac{1}{4}|F^{T}F|^{2}\right]$
$\displaystyle=:$
$\displaystyle\frac{3}{4}|a|^{4/3}+2(\mathsf{G}+\mathsf{L}/2)h(\tilde{A},B).$

where
$\tilde{A}=\frac{1}{2\mathsf{G}+\mathsf{L}}\left(A-(\mathsf{G}+\frac{\mathsf{L}d}{2})B\right)$.
Next, we estimate $h$ from below.

  1. 1.

If $|B|^{3}\leq\frac{27}{8^{3}}|\tilde{A}|$ we find for
$F=|\tilde{A}|^{-2/3}\tilde{A}$ that

$\displaystyle h(\tilde{A},B)$
$\displaystyle\geq|\tilde{A}|^{4/3}-|\tilde{A}||B|-\frac{1}{4}|F|^{4}$
$\displaystyle=\frac{3}{4}|\tilde{A}|^{4/3}-|\tilde{A}||B|\geq\frac{3}{8}|\tilde{A}|^{4/3}\geq\frac{3}{16}(|\tilde{A}|^{4/3}+|B|^{4}),$

where we used that $|F^{T}F|^{2}\leq|F|^{4}$.

  2. 2.

If $\frac{27}{8^{3}}|\tilde{A}|<|B|^{3}\leq\frac{8d}{9}|\tilde{A}|$ we find
for $F=\frac{1}{\sqrt{3}}|B|^{-1/2}|\tilde{A}|^{-1/2}\tilde{A}$ that

$\displaystyle h(\tilde{A},B)$
$\displaystyle\geq\frac{1}{\sqrt{3}}|\tilde{A}|^{3/2}|B|^{-1/2}-\frac{1}{3\sqrt{3}}|\tilde{A}|^{3/2}|B|^{-1/2}-\frac{1}{36}|\tilde{A}|^{2}|B|^{-2}$
$\displaystyle\geq\left(\frac{2}{3\sqrt{3}}-\frac{1}{36}\sqrt{\frac{8^{3}}{27}}\right)|\tilde{A}|^{3/2}|B|^{-1/2}$
$\displaystyle\geq\left(\frac{2}{3\sqrt{3}}-\frac{1}{36}\sqrt{\frac{8^{3}}{27}}\right)\left(\frac{1}{2}\left(\frac{9}{8d}\right)^{1/6}|\tilde{A}|^{4/3}+\frac{1}{2}\left(\frac{9}{8d}\right)^{3/2}|B|^{4}\right)$
$\displaystyle\geq c\left(|\tilde{A}|^{4/3}+|B|^{4}\right),$

where $c>0$ is a universal constant.

  3. 3.

If $\frac{8d}{9}|\tilde{A}|<|B|^{3}$ we find for $F=3B$ using
$\frac{1}{d}|B|^{4}\leq|B^{T}B|^{2}\leq|B|^{4}$ that

$\displaystyle h(\tilde{A},B)$
$\displaystyle\geq-3|\tilde{A}||B|+27|B^{T}B|^{2}-\frac{81}{4}|B^{T}B|^{2}$
$\displaystyle\geq\frac{27}{4d}|B|^{4}-3|\tilde{A}||B|$
$\displaystyle\geq\frac{27}{8d}|B|^{4}\geq\frac{27}{16d}|B|^{4}+\frac{27}{16d}\left(\frac{8d}{9}\right)^{4/3}|\tilde{A}|^{4/3}\geq
c\left(|\tilde{A}|^{4/3}+|B|^{4}\right),$

where $c>0$ is a universal constant.

Consequently, we find for a constant $c>0$ (that depends on $\mathsf{G}$ and
$\mathsf{L}$ and might change from line to line) that

$\displaystyle g(A,a,B)$ $\displaystyle\geq
c\left(|a|^{4/3}+|\tilde{A}|^{4/3}+|B|^{4}\right)$ $\displaystyle\geq
c\left(|a|^{4/3}+|A|^{4/3}-|B|^{4/3}+|B|^{4}\right)$ $\displaystyle\geq
c\left(|a|^{4/3}+|A|^{4/3}+|B|^{4}\right)-\frac{1}{c}.$ (28)

∎

######  Remark 5.1.

Similarly, one can show that it also holds

$g(A,a,B)\leq C\left(|a|^{4/3}+|A|^{4/3}+|B|^{4}+1\right),$

which shows that the dual functional $I=\tilde{S}_{H}$ is well-defined and its
minimization is non-trivial.

Let us now show the coercivity of the dual functional $I=\tilde{S}_{H}$.

######  Proposition 5.3.

The functional $I$ as defined in (24) is weakly coercive i.e., for every
sequence $(\mu_{k},\lambda_{k})_{k}\subseteq\\{(\mu,\lambda)\in
L^{4/3}(\Omega;\mathbb{R}^{d\times d})\times
W^{1,4}(\Omega;\mathbb{R}^{d}):\operatorname{div}(\mu)\in
L^{4/3}(\Omega;\mathbb{R}^{d}),\mu n=0\text{ on }\partial\Omega_{p}\text{ and
}\lambda=0\text{ on }\partial\Omega_{y}\\}$ such that
$\sup_{k}I(\mu_{k},\lambda_{k})<\infty$ there exists a (not relabeled)
subsequence and $(\mu,\lambda)\in\\{(\mu,\lambda)\in
L^{4/3}(\Omega;\mathbb{R}^{d\times d})\times
W^{1,4}(\Omega;\mathbb{R}^{d}):\operatorname{div}(\mu)\in
L^{4/3}(\Omega;\mathbb{R}^{d}),\mu n=0\text{ on }\partial\Omega_{p}\text{ and
}\lambda=0\text{ on }\partial\Omega_{y}\\}$ such that

$\mu_{k}\rightharpoonup\mu\text{ in }L^{4/3}(\Omega;\mathbb{R}^{d\times
d}),\;\operatorname{div}(\mu_{k})\rightharpoonup\operatorname{div}(\mu)\text{
in }L^{4/3}(\Omega;\mathbb{R}^{d})\text{ and
}\lambda_{k}\rightharpoonup\lambda\text{ in }W^{1,4}(\Omega;\mathbb{R}^{d}).$

###### Proof.

Let us first recall Poincaré’s inequality $\|\lambda\|_{W^{1,4}}\leq
C\|\nabla\lambda\|_{L^{4}}$ for all $\lambda\in
W^{1,4}(\Omega;\mathbb{R}^{d})$ such that $\lambda=0$ on $\partial\Omega_{y}$.
Hence, we obtain from Proposition 5.2 that

$\displaystyle I(\mu_{k},\lambda_{k})\geq$ $\displaystyle
c\left(\|\operatorname{div}(\mu_{k})\|_{L^{4/3}}^{4/3}+\|\mu_{k}\|_{L^{4/3}}^{4/3}+\|\nabla\lambda_{k}\|_{L^{4}}^{4}\right)-\frac{1}{c}\mathcal{L}^{d}(\Omega)$
$\displaystyle-\|\mu_{k}n\|_{(W^{1-1/4,4}(\partial\Omega))^{\prime}}\|y^{(b)}\|_{W^{1-1/4,4}(\partial\Omega)}-\|\lambda_{k}\|_{L^{4}(\partial\Omega)}\|p^{(b)}\|_{L^{4/3}(\partial\Omega)}$
$\displaystyle\geq$ $\displaystyle
c\left(\|\operatorname{div}(\mu_{k})\|_{L^{4/3}}^{4/3}+\|\mu_{k}\|_{L^{4/3}}^{4/3}+\|\lambda_{k}\|_{W^{1,4}}^{4}\right)-\frac{1}{c}\mathcal{L}^{d}(\Omega)$
$\displaystyle-C\left(\|\operatorname{div}(\mu_{k})\|_{L^{4/3}}+\|\mu_{k}\|_{L^{4/3}}+\|\lambda_{k}\|_{W^{1,4}}\right).$

Since $\sup_{k}I(\mu_{k},\lambda_{k})<\infty$, it follows immediately that
$\sup_{k}\|\operatorname{div}(\mu_{k})\|_{L^{4/3}}+\|\mu_{k}\|_{L^{4/3}}+\|\lambda_{k}\|_{W^{1,4}}<\infty$.
This implies the existence of $\mu\in L^{4/3}(\Omega;\mathbb{R}^{d\times d})$,
$\nu\in L^{4/3}(\Omega;\mathbb{R}^{d})$ and $\lambda\in
W^{1,4}(\Omega;\mathbb{R}^{d})$ such that (up to a not relabeled subsequence)
it holds $\mu_{k}\rightharpoonup\mu$ in $L^{4/3}$,
$\operatorname{div}(\mu_{k})\rightharpoonup\nu$ in
$L^{4/3}(\Omega;\mathbb{R}^{d})$ and $\lambda_{k}\rightharpoonup\lambda$ in
$W^{1,4}(\Omega;\mathbb{R}^{d})$. Then it is straightforward to check that
$\operatorname{div}(\mu)=\nu\in L^{4/3}(\Omega;\mathbb{R}^{d})$. Moreover, as
the trace operators are weakly continuous it follows that $\mu n=0$ on
$\partial\Omega_{p}$ and $\lambda=0$ on $\partial\Omega_{y}$ in the sense of
traces as discussed at the beginning of this section. ∎

###  5.2 Incompressible Neo-Hookean model in two dimensions

Let us now consider the stress in an incompressible Neo-Hookean model

$P(F,p)=\operatorname{det}(F)\left(\mathsf{G}\,F-pF^{-T}\right),$

where $\mathsf{G}>0$ is the shear modulus, $\operatorname{det}(F)=1$ and $p$
is the pressure which is to be determined. For simplicity, let $\mathsf{G}=1$.
Moreover, set

$\tilde{P}(F,p)=F-p\operatorname{cof}(F).$

Note that in two dimensions $\tilde{P}$ depends linearly on $F$ and for $F$
with $\operatorname{det}(F)=1$ it holds $\tilde{P}(F,p)=P(F,p)$. Then for
$y^{(b)}\in W^{1,4}(\Omega;\mathbb{R}^{d})$ and $p^{(b)}\in
L^{4/3}(\partial\Omega;\mathbb{R}^{d})$ we look for weak solutions to

$\begin{cases}F=\nabla y\text{ in }\Omega,\\\ \operatorname{det}(F)=1\text{ in
}\Omega,\\\ \operatorname{div}\tilde{P}(F,p)=0\text{ in }\Omega,\\\
\tilde{P}(F,p)\cdot n=p^{(b)}\text{ on }\partial\Omega_{p},\\\ y=y^{(b)}\text{
on }\partial\Omega_{y}.\end{cases}$

Namely, we seek to find $y\in L^{2}(\Omega;\mathbb{R}^{2})$ and $F\in
L^{2}(\Omega;\mathbb{R}^{2\times 2})$ such that it holds for all $\mu\in
L^{2}(\Omega;\mathbb{R}^{2\times 2})$ with $\operatorname{div}(\mu)\in
L^{2}(\Omega;\mathbb{R}^{2})$, $\lambda\in W^{1,4}(\Omega;\mathbb{R}^{2})$,
and $\vartheta\in L^{\infty}(\Omega)$ that

$\int_{\Omega}\operatorname{div}(\mu)\cdot
y+\mu:F+\vartheta(\operatorname{det}(F)-1)+\nabla\lambda:\tilde{P}(F)\,dx=\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{1}+\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{1}.$

As explained for the Saint Venant-Kirchhoff model, the boundary integrals have
to be understood in the sense of traces.

Note that compared to the general discussion in Section 2 the additional
constraint $\operatorname{det}(F)=1$ is treated through an extra dual variable
$\vartheta$, whereas the pressure $p$ is an additional primal variable.
Following the same approach, we define the function $g:\mathbb{R}^{d\times
d}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times
d}\times\mathbb{R}\to\mathbb{R}\cup\\{+\infty\\}$

$g(A,a,B,s)=\sup_{F\in\mathbb{R}^{d\times
d},y\in\mathbb{R}^{d},p\in\mathbb{R}}a\cdot
y+s(\operatorname{det}(F)-1)+A:F+B:P(F,p)-H(y,p,F),$

where $H(y,p,F)=\frac{1}{2}y^{2}+\frac{1}{2}|F|^{2}+\frac{1}{4}p^{4}$.
Correspondingly, we define the functional

$\displaystyle I:\big{\\{}(\mu,\lambda,\vartheta)\in
L^{4/3}(\Omega;\mathbb{R}^{d\times d})$ $\displaystyle\times
W^{1,4}(\Omega;\mathbb{R}^{d})\times L^{\infty}(\Omega):$
$\displaystyle\operatorname{div}(\mu)\in L^{4/3}(\Omega;\mathbb{R}^{d}),\mu
n=0\text{ on }\partial\Omega_{p}\text{ and }\lambda=0\text{ on
}\partial\Omega_{y}\big{\\}}\to\mathbb{R}$

by

$I(\mu,\lambda,\vartheta)=\int_{\Omega}g(\mu,\operatorname{div}(\mu),\nabla\lambda,\vartheta)\,dx-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{d-1}-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,\mathcal{H}^{d-1}.$ (29)

By the same argument as for the Saint Venant - Kirchhoff material in
Proposition 5.1 it holds for all
$(\mu,\lambda,\vartheta)\in\\{(\mu,\lambda,\vartheta)\in
L^{4/3}(\Omega;\mathbb{R}^{d\times d})\times
W^{1,4}(\Omega;\mathbb{R}^{d})\times
L^{\infty}(\Omega):\operatorname{div}(\mu)\in
L^{4/3}(\Omega;\mathbb{R}^{d}),\mu n=0\text{ on }\partial\Omega_{p}\text{ and
}\lambda=0\text{ on }\partial\Omega_{y}\\}$ that

$\displaystyle I(\mu,\lambda,\vartheta)$ (30) $\displaystyle=$
$\displaystyle\sup_{y\in L^{2}(\Omega;\mathbb{R}^{2}),F\in
L^{2}(\Omega;\mathbb{R}^{2\times 2}),p\in
L^{4}(\Omega)}\int_{\Omega}\operatorname{div}(\mu)\cdot
y+\vartheta(\operatorname{det}(F)-1)+\mu:F+\nabla\lambda:P(F,p)-H(y,p,F)\,dx$
$\displaystyle-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{1}-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{1}.$ (31)

The right hand side of the above equation is the analogue of $\tilde{S}_{H}$
including the extra variable $\vartheta$.

As in the section above, in order to show existence of minimizers of $I$ it
suffices to show that $I$ is coercive. Again, we first prove lower bounds for
the function $g$.

######  Proposition 5.4.

There exist constants $C>c>0$ such that the following is true.

  1. 1.

If $|s|>1$ then $g(A,a,B,s)=+\infty$ for all $(A,a,B)\in\mathbb{R}^{d\times
d}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}$.

  2. 2.

If $|s|\leq 1$ then it holds for $(A,a,B)\in\mathbb{R}^{d\times
d}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}$

$g(A,a,B,s)\geq c(|a|^{2}+|A+B|^{2}+|B|^{4})-1$

  3. 3.

If $|s|\leq 1$ then it holds for $(A,a,B)\in\mathbb{R}^{d\times
d}\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}$ that

$g(A,a,B,s)\leq
C\left(|a|^{2}+\frac{1}{1-|s|}|A+B|^{2}+\frac{1}{(1-|s|)^{2}}|B|^{4}+1\right).$

###### Proof.

First note that by optimizing in $y$ we find

$\displaystyle g(A,a,B,s)=$ $\displaystyle\sup_{F\in\mathbb{R}^{d\times
d},y\in\mathbb{R}^{d},p\in\mathbb{R}}a\cdot
y+s(\operatorname{det}(F)-1)+A:F+B:P(F,p)-H(y,p,F)$ $\displaystyle=$
$\displaystyle\frac{1}{2}|a|^{2}-s+\sup_{F\in\mathbb{R}^{d\times
d},p\in\mathbb{R}}s\operatorname{det}(F)+A:F+B:P(F,p)-\frac{1}{4}p^{4}-\frac{1}{2}|F|^{2}$
$\displaystyle=:$ $\displaystyle\frac{1}{2}|a|^{2}-s+h(A,B,s).$ (32)

We will now only work with the function $h$. We observe that it holds for all
$F\in\mathbb{R}^{2\times 2}$ that
$F:\operatorname{cof}(F)=2\operatorname{det}(F)$. In particular,
$|\operatorname{\det}(F)|\leq\frac{1}{2}|F|^{2}$. Let us now assume that
$s>1$. Then consider for $t>0$ the matrix $F=tId$. In addition, set $p=0$. It
follows that

$h(A,B,s)\geq st^{2}-2|A+B|t-t^{2}\stackrel{{\scriptstyle
t\to\infty}}{{\longrightarrow}}+\infty.$

The case $s<-1$ can be treated similarly. This shows 1. Next, let us consider
$|s|\leq 1$. First we assume that $\frac{1}{4\sqrt{8}}|B|^{2}\leq|A+B|$. Then
set $F=(A+B)/2$ and $p=0$ to obtain

$\displaystyle
h(A,B,s)\geq(A+B):F-\frac{1}{2}|F|^{2}+s\operatorname{det}(F)\geq\frac{|A+B|^{2}}{2}-|F|^{2}$
$\displaystyle=\frac{|A+B|^{2}}{4}$ (33)
$\displaystyle\geq\frac{1}{512}\left(|A+B|^{2}+|B|^{4}\right).$

For $|A+B|\leq\frac{1}{4\sqrt{8}}|B|^{2}$ consider $F\in\mathbb{R}^{2\times
2}$ such that $\operatorname{cof}(F)=-\frac{1}{\sqrt{8}}|B|B$ and
$p=\frac{1}{\sqrt{2}}|B|$. Then

$\displaystyle h(A,B,s)$
$\displaystyle\geq(A+B):F-B:(p\operatorname{cof}(F))+s\operatorname{det}(F)-\frac{1}{4}p^{4}-\frac{1}{2}|F|^{2}$
$\displaystyle\geq(A+B):F-B:(p\operatorname{cof}(F))-\frac{1}{4}p^{4}-|F|^{2}$
$\displaystyle\geq-\frac{1}{\sqrt{8}}|A+B|\,|B|^{2}+\frac{1}{4}|B|^{4}-\frac{1}{16}|B|^{4}-\frac{1}{8}|B|^{4}$
$\displaystyle\geq-\frac{1}{32}|B|^{4}+\frac{1}{16}|B|^{4}=\frac{1}{32}|B|^{4}\geq\frac{1}{64}\left(|A+B|^{2}+|B|^{4}\right)$
(34)

Combining (32), (33) and (34) yields

$g(A,a,B,s)\geq\frac{1}{2}|a|^{2}-|s|+h(A,B,s)\geq
c(|a|^{2}+|A+B|^{2}+|B|^{4})-1,$

which is 2.

It remains to show 3. We estimate

$\displaystyle h(A,B,s)$ $\displaystyle\leq\sup_{F\in\mathbb{R}^{2\times
2},p\in\mathbb{R}}F:(A+B-p\operatorname{cof}(B))-\frac{1}{4}p^{4}-\frac{1-|s|}{2}|F|^{2}$
$\displaystyle=\sup_{p\in\mathbb{R}}\frac{1}{2(1-|s|)}|A+B-p\operatorname{cof}(B)|^{2}-\frac{1}{4}p^{4}$
$\displaystyle\leq\frac{1}{(1-|s|)}|A+B|^{2}+\sup_{p\in\mathbb{R}}\frac{p^{2}}{(1-|s|)}|\operatorname{cof}(B)|^{2}-\frac{1}{4}p^{4}$
$\displaystyle\leq\frac{1}{(1-|s|)}|A+B|^{2}+\frac{1}{(1-|s|)^{2}}|B|^{4}.$

Combining this with (32) yields 3. ∎

Eventually, we state the coercivity statement for the dual functional $I$.

######  Proposition 5.5.

The functional $I$ as defined in (29) is weakly coercive i.e., for every
sequence
$(\mu_{k},\lambda_{k},\vartheta_{k})_{k}\subseteq\\{(\mu,\lambda,\vartheta)\in
L^{4/3}(\Omega;\mathbb{R}^{d\times d})\times
W^{1,4}(\Omega;\mathbb{R}^{d})\times
L^{\infty}(\Omega):\operatorname{div}(\mu)\in
L^{4/3}(\Omega;\mathbb{R}^{d}),\mu n=0\text{ on }\partial\Omega_{p}\text{ and
}\lambda=0\text{ on }\partial\Omega_{y}\\}$ such that
$\sup_{k}I(\mu_{k},\lambda_{k},\vartheta_{k})<\infty$ there exists a (not
relabeled) subsequence and
$(\mu,\lambda,\vartheta)\in\\{(\mu,\lambda,\vartheta)\in
L^{4/3}(\Omega;\mathbb{R}^{d\times d})\times
W^{1,4}(\Omega;\mathbb{R}^{d})\times
L^{\infty}(\Omega):\operatorname{div}(\mu)\in
L^{4/3}(\Omega;\mathbb{R}^{d}),\mu n=0\text{ on }\partial\Omega_{p}\text{ and
}\lambda=0\text{ on }\partial\Omega_{y}\\}$ such that

$\displaystyle\mu_{k}\rightharpoonup\mu\text{ in
}L^{2}(\Omega;\mathbb{R}^{d\times
d}),\;\operatorname{div}(\mu_{k})\rightharpoonup\operatorname{div}(\mu)\text{
in }L^{2}(\Omega;\mathbb{R}^{d}),\;\lambda_{k}\rightharpoonup\lambda\text{ in
}W^{1,4}(\Omega;\mathbb{R}^{d})$ $\displaystyle\text{ and
}\vartheta_{k}\stackrel{{\scriptstyle*}}{{\rightharpoonup}}\vartheta\text{ in
}L^{\infty}(\Omega).$

###### Proof.

The proof is similar to the proof of Proposition 5.3 using Proposition 5.4 1.
and 2. ∎

We end this section with a brief comment on changing the base states.

######  Remark 5.2.

Let us fix $\tilde{y}\in L^{2}(\Omega;\mathbb{R}^{2})$ and $\tilde{F}\in
L^{2}(\Omega;\mathbb{R}^{2\times 2})$ and note for
$\tilde{H}(F,y,p;x)=H(F-\tilde{F}(x),y-\tilde{y}(x),p)$ that

$\displaystyle S_{\tilde{H}}(\lambda,\mu,\vartheta)$ $\displaystyle:=$
$\displaystyle\sup_{y\in L^{2},F\in L^{2},p\in
L^{4}}\int_{\Omega}\operatorname{div}(\mu)\cdot
y+\vartheta(\operatorname{det}(F)-1)+\mu:F+\nabla\lambda:P(F,p)-\tilde{H}(F,y,p;x)\,dx$
$\displaystyle-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{1}-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{1}$ $\displaystyle=$ $\displaystyle\sup_{y\in L^{2},F\in
L^{2},p\in L^{4}}\int_{\Omega}\operatorname{div}(\mu)\cdot
y+\vartheta(\operatorname{det}(F)-1)+\mu:F+\nabla\lambda:P(F,p)-\tilde{H}(F,y,p;x)\,dx$
$\displaystyle-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{1}-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{1}$ $\displaystyle=$ $\displaystyle\sup_{y\in L^{2},F\in
L^{2},p\in L^{4}}\int_{\Omega}(\operatorname{div}(\mu)-\tilde{y})\cdot
y+\vartheta(\operatorname{det}(F)-1)+(\mu+\nabla\lambda-p\operatorname{cof}(\nabla\lambda)-\tilde{F}):F-H(F,y,p)\,dx$
$\displaystyle+\int_{\Omega}\frac{1}{2}|\tilde{F}|^{2}+\frac{1}{2}|\tilde{y}|^{2}\,dx-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,d\mathcal{H}^{1}-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{1}$ $\displaystyle=$
$\displaystyle\int_{\Omega}g(\mu-\tilde{F},\operatorname{div}(\mu)-\tilde{y},\nabla\lambda,\vartheta)\,dx-\int_{\partial\Omega_{y}}y^{(b)}\cdot\mu
n\,d\mathcal{H}^{d-1}-\int_{\partial\Omega_{p}}\lambda\cdot
p^{(b)}\,\mathcal{H}^{d-1}+\int_{\Omega}\frac{1}{2}|\tilde{F}|^{2}+\frac{1}{2}|\tilde{y}|^{2}\,dx.$

This functional has the same coercivity and lower-semicontinuity properties as
$\tilde{S}_{H}$. Hence, existence of minimizers for this functional is also
guaranteed. The question of whether such a minimzer gives rise to a weak
solution of the primal problem hinges on the regularity of the function $g$ on
the values of a minimizer. However, Proposition 5.4 shows that at least at
$\vartheta=\pm 1$ the function $g$ is not differentiable. Making good choices
for $\tilde{F}$ and $\tilde{y}$ could enable us to guarantee that a minimizer
$(\lambda,\mu,\vartheta)$ of $\tilde{S}_{\tilde{H}}$ only takes values such
that $g$ is actually differentiable in the points
$(\mu-\tilde{F},\operatorname{div}(\mu)-\tilde{y},\nabla\lambda,\vartheta)$.

##  6 A computational case study: non-convex elasticity in 1-d

In this section, we present computational approximations of our dual scheme
applied to nonconvex elastostatics and elastodynamics in 1 space dimension and
time; this provides intuition for its capabilities and evidence for its
practical feasibility.

The specific variational and PDE problems we study in the elastostatic case
are described in Sec. 6.1. The finite element method based algorithm used to
solve the problem is developed in Sec. 6.2. Results of the computational
studies are presented and discussed in Sec. 6.3. In particular, in Sec. 6.3.4
is presented and discussed the stability of some dual solutions to
elastodynamic problems that are ‘Hadamard unstable’ in the primal, classical
elasticity formulation, in the sense of lacking continuous dependence on
initial data because some regions of the domain have initial conditions where
hyperbolicity is lost. While it is not our goal to discuss the (de)merits of
conventional elasticity theory to model such situations, we give a physical
interpretation of the problem that is studied. This section also contains a
brief description of the extension of the computational scheme to
elastodynamics.

###  6.1 Formulation

We computationally approximate solutions to a double-well elastostatics
problem in 1 space dimension without any higher-order (strain-gradient)
regularization. The objective is to compute critical points of the following
energy functional:

$I=\int_{0}^{1}\left(\left(u_{x}-1\right)^{2}-1\right)^{2}+\left(u-\alpha
x\right)^{2}\,dx;\qquad u(0)=0;\quad u(1)=\alpha^{*}.$ (35)

We use the notation $\partial_{i}a=a_{i},i=x,t$, and $\alpha,\alpha^{*}$ are
specified constants. The function $u$ corresponds to the displacement. The
Euler-Lagrange equations of (35) are

$\displaystyle u_{x}-e$ $\displaystyle=0$ (36a) $\displaystyle
2\partial_{x}\left((e-1)^{3}-(e-1)\right)-\left(u-\alpha x\right)$
$\displaystyle=0$ (36b) $\displaystyle u(0)=0$ $\displaystyle;\qquad
u(1)=\alpha^{*},$ (36c)

where $e$ corresponds to the strain. The stress, $\sigma$, for this model is
given by

$\sigma=4(e-1)\left((e-1)^{2}-1)\right).$

Using the dual (Lagrange multiplier) fields $\lambda$ and $\mu$ for the two
equations, we define the pre-dual functional

$\widehat{S}\,\left[u,\lambda,e,\mu\right]\,=\,\int_{0}^{1}\,\,(-u\lambda_{x}-\lambda
e-2((e-1)^{3}-(e-1))\mu_{x}-\mu(u-\alpha x)+H(u,e))\,dx+\alpha^{*}\lambda(1),$

where the requirement on the function $H$ is that, defining

$\mathcal{D}:=\left(\lambda,\lambda_{x},\mu,\mu_{x}\right),$

the following equations are solvable for $u$, $e$ in terms of $\mathcal{D}:$

$\displaystyle-\lambda_{x}-\mu+\partial_{u}H$ $\displaystyle=0$ (37)
$\displaystyle-2(3(e-1)^{2}-1)\mu_{x}-\lambda+\partial_{e}H$
$\displaystyle=0.$

Choosing a ‘shifted quadratic’ and ‘shifted cubic’ form for the potential $H$:

$H(u,e;x)=\frac{1}{2}\left(c_{u}(u-\bar{u}(x))^{2}\right)+\frac{1}{2}\left(c_{e}(e-\bar{e}(x))^{2}\right)+\frac{1}{3}\left(c_{e}(e-\bar{e}(x))^{3}\right),$

we obtain the following dual-to-primal (DtP) mapping equations

$\displaystyle
u^{(H)}(\lambda_{x},\mu,x)=\frac{\lambda_{x}+\mu}{c_{u}}+\bar{u}(x)$ (38)
$\displaystyle-6\left(e^{(H)}(\mu_{x},\lambda,x)-1\right)^{2}\mu_{x}+2\mu_{x}-\lambda$
$\displaystyle+c_{e}\left(e^{(H)}(\mu_{x},\lambda,x)-\bar{e}(x)\right)\left|e^{(H)}(\mu_{x},\lambda,x)-\bar{e}(x)\right|+c_{e}\left(e^{(H)}(\mu_{x},\lambda,x)-\bar{e}(x)\right)=0.$

Henceforth, the function $u^{(H)}$ will be referred to as $\hat{u}$ with the
understanding that the symbol carries argument $\mathcal{D}$ and $\bar{u}$.
This also applies to $e^{(H)}$:

$\hat{u}(x):=u^{(H)}(\lambda_{x}(x),\mu(x),x);\qquad\hat{e}(x):=e^{(H)}(\mu_{x}(x),\lambda(x),x).$

It is clear from (37) and (38) that $\hat{u}=\bar{u}$ and $\hat{e}=\bar{e}$ is
a value of the DtP mapping for $\mathcal{D}=0$.

We substitute the DtP mapping (38) in the functional $\widehat{S}$ to define
the dual functional

$S\,\left[\lambda,\mu\right]\,=\,\int_{0}^{1}\,\,\left(-\hat{u}\lambda_{x}-\lambda\hat{e}-2\left((\hat{e}-1)^{3}-(\hat{e}-1)\right)\mu_{x}-\mu\left(\hat{u}-\alpha
x\right)+H\left(\hat{u},\hat{e}\right)\right)\,dx+\alpha^{*}\lambda(1).$ (39)

The requirement (37) ensures that the dual E-L system of (39) is simply the
primal system (36), with the replacements $u\to\hat{u},e\to\hat{e}$. Thus, the
critical points for the dual functional $S$ (39) is also a critical point of
the primal functional $I$ (35), interpreted through the DtP mapping. The
algorithm employed to obtain such critical points is discussed next.

###  6.2 Finite element algorithm

Denoting the dual fields as

$D=(D_{1},D_{2}):=(\lambda,\mu),$

the following weak form suffices to compute solutions to the dual problem
(40), which also constitutes the first variation of the dual functional (39)
set to 0:

$\displaystyle R[D;\delta D]$
$\displaystyle=-\int_{0}^{1}\,\hat{u}\delta\lambda_{x}\,dx\,-\int_{0}^{1}\,\hat{e}\delta\lambda\,dx$
(40)
$\displaystyle\quad-\int_{0}^{1}\,2\left(\left(\hat{e}-1\right)^{3}-(\hat{e}-1)\right)\delta\mu_{x}\,dx\,-\int_{0}^{1}\,\left(\hat{u}-\alpha
x\right)\delta\mu\,dx+\alpha^{*}\delta\lambda(1)=0;$
$\displaystyle\qquad\mu(0)=0,\qquad\mu(1)=0.$

A Dirichlet boundary condition is applied only on $\mu$. A Neumann boundary
condition is applied on $\lambda$, corresponding to the Dirichlet boundary
condition on $u$ in the primal problem.

We use the Finite Element (FE) method to discretize the problem. A linear span
of globally continuous, piecewise smooth finite element shape functions
corresponding to a FE mesh for the domain is used to achieve this
discretization. These shape functions are represented by
$N^{\left(\cdot\right)}$, where ($\cdot$) denotes the index of the node under
consideration. We use $C^{0}$ FE shape functions to approximate the dual
fields. Hence, the direct evaluation of primal fields through the DtP mapping
exhibit discontinuities in the domain in general, even when approximating
continuous primal solutions. This is because the DtP mappings involve
derivatives of the dual fields. To deal with this feature, we use an $L^{2}$
projection of the DtP mapping generated primal $u$ field from the dual
solution to define the corresponding continuous primal field.

The discretized dual fields and test functions are expressed as

$D_{i}(x)=D^{A}_{i}N^{A}(x);\qquad\delta D_{i}(x)=\delta
D^{A}_{i}N^{A}(x);\qquad dD_{i}(x)=dD^{A}_{i}N^{A}(x),\qquad i=1,2,$

where $D^{A}_{i}$ denote the finite element nodal degrees of freedom (which
are the coefficients of a linear combination of a finite dimensional
representation of the dual fields). The discretized version of (40) generates
a discrete residual $R^{A}_{i}(\mathcal{D})$ given by

$R[D;\delta D]=\delta D^{A}_{i}R^{A}_{i}(\mathcal{D})$

and its variation in the direction $dD$, the Jacobian

$J[D;\delta D,dD]=\delta D^{A}_{i}\,J^{AB}_{ij}(\mathcal{D})\,dD^{B}_{j}$

given by

$\displaystyle J[D;\delta D,dD]\,=$
$\displaystyle\,-\int_{0}^{1}\,\delta\lambda_{x}\left[\frac{\partial\hat{u}}{\partial\lambda_{x}}\,d\lambda_{x}+\frac{\partial\hat{u}}{\partial\mu}\,d\mu\right]\,dx\,-\int_{0}^{1}\,\delta\lambda\left[\frac{\partial\hat{e}}{\partial\mu_{x}}\,d\mu_{x}+\frac{\partial\hat{e}}{\partial\lambda}\,d\lambda\right]\,dx$
(41)
$\displaystyle\,-\int_{0}^{1}\,2\left(3\left(\hat{e}-1\right)^{2}-1\right)\delta\mu_{x}\left[\frac{\partial\hat{e}}{\partial\mu_{x}}\,d\mu_{x}+\frac{\partial\hat{e}}{\partial\lambda}\,d\lambda\right]\,dx$
$\displaystyle\,-$
$\displaystyle\int_{0}^{1}\,\delta\mu\left[\frac{\partial\hat{u}}{\partial\lambda_{x}}\,d\lambda_{x}+\frac{\partial\hat{u}}{\partial\mu}\,d\mu\right]\,dx$

that are used in a Newton Raphson scheme. Here, $R^{A}_{i},J^{AB}_{ij}$ are
given by the expressions

$\displaystyle
R^{A}_{1}=\,-\int_{0}^{1}\,\hat{u}N^{A}_{x}\,dx\,-\int_{0}^{1}\,\hat{e}N^{A}\,dx$
$\displaystyle
R^{A}_{2}=-\int_{0}^{1}\,2\left(\left(\hat{e}-1\right)^{3}-(\hat{e}-1)\right)N^{A}_{x}\,dx\,-\int_{0}^{1}\,\left(\hat{u}-\alpha
x\right)N^{A}\,dx$

and

$\displaystyle
J^{AB}_{11}=\,-\int_{0}^{1}\,\frac{\partial\hat{u}}{\partial\lambda_{x}}N^{A}_{x}N^{B}_{x}\,dx\,-\int_{0}^{1}\,\frac{\partial\hat{e}}{\partial\lambda}N^{A}N^{B}\,dx$
$\displaystyle
J^{AB}_{12}=\,-\int_{0}^{1}\,\frac{\partial\hat{e}}{\partial\mu_{x}}N^{A}N^{B}_{x}\,dx\,-\int_{0}^{1}\,\frac{\partial\hat{u}}{\partial\mu}N^{A}_{x}N^{B}\,dx$
$\displaystyle
J^{AB}_{21}=\,-\int_{0}^{1}\,\frac{\partial\hat{u}}{\partial\lambda_{x}}N^{A}N^{B}_{x}\,dx\,-\int_{0}^{1}\,2(3(\hat{e}-1)^{2}-1)\,\frac{\partial\hat{e}}{\partial\lambda}N^{A}_{x}N^{B}\,dx$
$\displaystyle
J^{AB}_{22}=\,-\int_{0}^{1}\,2(3(\hat{e}-1)^{2}-1)\,\frac{\partial\hat{e}}{\partial\mu_{x}}N^{A}_{x}N^{B}_{x}\,dx\,-\int_{0}^{1}\,\frac{\partial\hat{u}}{\partial\mu}N^{A}N^{B}\,dx$

The Newton-Raphson scheme involves solving the following matrix equation for
the $k^{th}$ correction to an initial guess $D^{B(0)}_{j}$

$\displaystyle-R^{A}_{i}\left(\mathcal{D}^{(k-1)}\right)$
$\displaystyle=J^{AB}_{ij}\left(\mathcal{D}^{(k-1)}\right)\,dD^{B}_{j}$
$\displaystyle D^{B(k)}_{j}$ $\displaystyle=D^{B(k-1)}_{j}+dD^{B}_{j},$

and this is continued until convergence of $|R(\mathcal{D})|$ to $0$, up to
the convergence threshold $tol$.

An initial guess for the dual fields $D$ is required by the Newton Raphson
(NR) algorithm. It is practical, and more feasible, to make an initial guess
for a primal fields $(\hat{u},\hat{e})$ with physical meaning, and infer
initial guesses for the dual fields consistent with it as a solution of the
DtP mapping. Choosing initial guesses $\hat{u}_{(i)}$ and $\hat{e}_{(i)}$

$\hat{u}_{(i)}=\bar{u},\qquad\hat{e}_{(i)}=\bar{e},$

the initial guess for the dual fields,

$D_{i}^{B}=0,$

satisfies the DtP mapping equations and this is considered by us as a good
initial guess to initiate the Newton-Raphson iterations. The algorithm
capturing these ideas has been outlined in Table. 1.

Some notational details related to our computational algorithm and results are
as follows:

  * •

$tol$ represents the tolerance for the iterative solve of the Newton Raphson
iterations; residuals with $|D|<tol$ are accepted as having approximately
satisfied the corresponding equations.

  * •

$u^{(t)},e^{(t)}$ will represent exact ‘target’ solutions for the context
under discussion.

  * •

We report the accuracy of our approximations in the $L^{1}$ norm:

$\|f\|_{1}=\int_{0}^{1}|f(x)|\,dx.$

  * •

In all our calculations we seek weak solutions of the E-L system so that jump
conditions (like traction continuity) are automatically imposed.

Algorithm

Initialization: Set $D^{B(0)}_{j}=0$. Choose a value for $tol$.

  1. 1.

Newton’s Method is used to evaluate the results. Set $k=1$

  2. For $k\geq 1$:

    1. i

Evaluate $R^{A}_{i}\left(\mathcal{D}^{(k-1)}\right)$ and
$J^{AB}_{ij}\left(\mathcal{D}^{(k-1)}\right)$.

    2. ii

Solve
$-R^{A}_{i}\left(\mathcal{D}^{(k-1)}\right)=J^{AB}_{ij}\left(\mathcal{D}^{(k-1)}\right)\,dD^{B}_{j}$,
with appropriate boundary conditions.

    3. iii

Set $D^{B(k)}_{j}=D^{B(k-1)}_{j}+dD^{B}_{j}$.

    4. iv

Evaluate $d^{(k)}=\max_{(A,i)}\,|R^{A}_{i}\left(\mathcal{D}^{(k)}\right)|$.

    5. if $d^{(k)}<tol$ then go to step 2

    6. else do $k=k+1$ and go to step 1i

  3. 2.

Evaluate the values of $\hat{u}$ and $\hat{e}$ at the nodes using the $L^{2}$
projection. These nodal values establish the solution.

Table 1: Algorithm to solve non-convex elastostatics in 1-d.

###  6.3 Results

The algorithm of Sec. 6.2 is utilized to solve key problems demonstrating the
theory. It is well-understood that due to the lack of lower semi-continuity of
the primal energy functional (35), limits of energy minimizing sequences do
not, in general, minimize the energy, i.e. the minimization problem does not
have any solution (for an elementary treatment, see [8]). The dual scheme, on
the other hand, defines the notion of a variational dual solution to the E-L
equations of the energy functional, as in Definition 1, Sec. 5 . Such
solutions need not have a connection to being minimizers of the primal energy
functional, and we compute such results in the following examples. Before
proceeding to the problems that follow, we make the observation that our goal
is to seek solutions to the system (36). In particular, homogeneous strain
profiles $u_{x}=\beta$ are solutions provided $\beta=\alpha=\alpha^{*}$.

####  6.3.1 Stress-free solution

We consider the primal energy functional

$I=\int_{0}^{1}\left(\left(u_{x}-1\right)^{2}-1\right)^{2}+(u-x)^{2}\,dx;\qquad
u(0)=0;\quad u(1)=1\ $

and compute solutions to its Euler-Lagrange equation (36) with
$\alpha=\alpha^{*}=1$. A stress-free equilibrium solution of the primal
problem is given by

$u^{(t)}=x,\qquad e^{(t)}=1,$

considered a target solution here, and we seek to compute it by the dual
scheme, probing its ‘basin of attraction’ by choosing initial guesses
sufficiently away from it. The boundary conditions chosen (profile shown in
Fig. 1) accommodate the stress-free configuration (which happens to be a limit
of an energy minimizing sequence). As mentioned in Sec. 6.2, we choose the
initial guess to coincide with the base state; we demonstrate the solution of
two dual problems, each with a different choice of base state as shown in Fig.
1(a) and Fig. 1(b). The base state provided for an $\bar{e}$ profile is the
derivative of the corresponding $\bar{u}$ profile. The $\bar{u}$ profiles
satisfy the displacement b.c.s of the target solution. The dual scheme is able
to recover the target solution shown in Fig. 1 without difficulty, even for an
initial guess amplitude that is $30\%$ away from the target. The $L^{1}$
displacement and strain errors relative to the target solution is summarized
in Table. 2 for the case corresponding to Fig. 1(a). These values are of the
same order for the other sub-case (Fig. 1(b)) as well.

(a)

(b)

(c)

Figure 1:  (a) and (b) show the initial guesses, $\bar{e}$, for two cases and (c) shows the solution for the primal fields. Fig. (c) contains both $\hat{u}$ and $\hat{e}$ fields with different y-axis labels. The color purple corresponds to information related to the displacement, while the color green corresponds to to that for strain. Number of Elements | $\|\hat{u}-u^{(t)}\|_{1}$ | $\|\hat{e}-e^{(t)}\|_{1}$
---|---|---
100 | $1\times 10^{-4}$ | $1\times 10^{-5}$
1600 | $4\times 10^{-7}$ | $2\times 10^{-7}$
8000 | $1\times 10^{-8}$ | $8\times 10^{-9}$
Table 2: Error reduction in $L^{1}$ norms of displacement and strain with mesh
refinement. These values are for the case corresponding to Fig. 1(a).

####  6.3.2 Stressed solutions with spatially homogeneous and inhomogeneous
strain configurations

We now consider the primal functional

$I=\int_{0}^{1}\left(\left(u_{x}-1\right)^{2}-1\right)^{2}+(u-0.5x)^{2}\,dx;\qquad
u(0)=0;\quad u(1)=0.5,$ (42)

with

$e^{(t)}=0.5,\qquad u^{(t)}=0.5x,$

as the exact target solution. We choose a base state near this solution as
shown in Fig. 2(a) and start from an initial guess coincident with the base
state. The base state provided for an $\bar{e}$ profile is the derivative of
the corresponding $\bar{u}$ profile. The $\bar{u}$ profiles satisfy the
displacement b.c.s of the target solution. The dual scheme is able to recover
the target solution, presented in Fig. 2(b). The $L^{1}$ displacement and
strain errors relative to the target solution are summarized in Table. 3.

(a)

(b)

Figure 2: (a) Base state $\bar{e}$ and (b) Solution for the primal fields. Number of Elements | $\|\hat{u}-u^{(t)}\|_{1}$ | $\|\hat{e}-e^{(t)}\|_{1}$
---|---|---
100 | $2\times 10^{-4}$ | $1\times 10^{-4}$
1600 | $1\times 10^{-6}$ | $4\times 10^{-7}$
8000 | $4\times 10^{-8}$ | $1\times 10^{-8}$
Table 3: Error reduction in $L^{1}$ norms of displacement and strain with mesh
refinement. These values are for the case corresponding to Fig. 2.

Next, we consider the energy functional

$I=\int_{0}^{1}\left(\left(u_{x}-1\right)^{2}-1\right)^{2}+(u-0.5x)^{2}\,dx;\qquad
u(0)=0;\quad u(1)=1.$

Clearly, a homogeneous strain solution is not possible in this case due to the
mismatch $\alpha\neq\alpha^{*}$ and hence a solution is not easy to guess, but
the expectation of the existence of a stressed, smooth solution of (36) seems
eminently reasonable (even though, as usual, energy minimizers for this
functional do not exist). However, a choice of a base state close to such an
expected solution is not easy to guess. The equilibrium solution obtained in
this case is shown in Fig. 4(b). The difference between the base state chosen
for this case and the solution can be observed in Fig. 4(a). Convergence is
obtained even with a $20\%$ difference in the average value of the initial
guess for strain and the obtained solution ($e$). The problem is solved with
varying mesh densities in the domain, starting with 8000 elements and
comparing to solutions obtained with 4000, 2000, and 100 elements. The $L^{1}$
norm of the difference between solutions obtained with different mesh
refinements is calculated and the results are presented in Table. 4, where
$\hat{u}^{(n_{e})}$ corresponds to the discrete displacement solution obtained
using $n_{e}$ number of elements and similarly for the strain, $\hat{e}$. As
the mesh is made finer the solutions appear to approach a fixed profile, this
being interpreted as an indication of convergence.

We have considered the question of finding weak solutions to the equations of
elastostatics (36) close to a guess for a limit of an energy minimizing
sequence for this problem. The displacement profile for such a limit is
expected to consist of a straight line of slope $\alpha=0.5$ up to some point
$x_{0}$ close to $x=1$, followed by an energy minimizing boundary layer
profile satisfying the boundary conditions $u(x_{0})=\alpha x_{0}$ and
$u(0)=\alpha^{*}=1$ (which can be obtained by solving the elastostatics
equations with these b.c.s for fixed $x_{0}$, and then optimizing w.r.t
$x_{0}$). However, since no considerations of force balance have entered, it
is very unlikely that the traction jump condition is satisfied (the jump
condition requires the discontinuity in stress, $\sigma$, in this 1-d case, to
vanish at the displacement kink at $x_{0}$ joining the constant strain segment
and boundary layer profile). Furthermore, the displacement profile for an
element of such a sequence for $0\leq x\leq x_{0}$ would consist of a typical
‘staircase’ profile of slopes $0,2$ along the backbone piece of slope
$\alpha$, as shown in Fig. 3. However, due to the presence of the forcing term
$(u-\alpha x)^{2}$ in the energy functional, constant strain segments do not
satisfy the E-L equations (36) (while the jump condition would be satisfied at
the interfaces of these ‘laminates’). Consequently, neither the elements of
such an energy minimizing sequence, nor the limit can be expected to strictly
satisfy the weak form of the E-L equations.

On setting up the base state and initial guess as the limit displacement
profile discussed above with $x_{0}\neq\frac{2}{3}$, we have not been able to
find solutions to the corresponding dual problem. On noting the fact that any
continuous displacement profile comprising constant strain segments of slope
$\\{0,2\\}$ does satisfy the traction jump condition, a staircase profile with
these slopes riding on a slope $\frac{1}{2}$ backbone up to
$x_{0}=\frac{2}{3}$ allows spanning the boundary layer with a constant strain
of $2$. Thus, this is a profile for which traction jumps vanish, and with
profiles close to this as a base state, the dual problem is able to find
approximations to the weak form of the E-L equations.

In this manner, we have succeeded in finding many equilibrium solutions by the
dual scheme close to a limit of a minimizing sequence, the limit in question
being a continuous displacement profile with a single kink at
$x_{0}=\frac{2}{3}$, with constant slopes of $\frac{1}{2}$ to the left and $2$
to the right of it. However, we have not succeeded in finding a dual (or a
primal) weak solution to (36) using this limit profile as a base state and/or
initial guess.

Figure 3: The displacement profile for the case $x_{0}=\frac{2}{3}$. The
staircase profile (green) rides on a line of slope $\frac{1}{2}$ (blue) until
$x_{0}$, followed by an orange line of slope 2.

(a)

(b)

(c)

Figure 4: (a) Base state $\bar{e}$, (b) Solution for the primal fields and (c) Difference between the linear profile for $u$ and the constant profile for $e$ is shown for clarity. Number of Elements ($n_{e}$) | $\|\hat{u}^{(n_{e})}-\hat{u}^{(8000)}\|_{1}$ | $\|\hat{e}^{(n_{e})}-\hat{e}^{(8000)}\|_{1}$
---|---|---
100 | $1\times 10^{-5}$ | $2\times 10^{-4}$
2000 | $4\times 10^{-8}$ | $5\times 10^{-7}$
4000 | $8\times 10^{-9}$ | $1\times 10^{-7}$
Table 4: Mesh convergence for computations shown in Fig. 4.

####  6.3.3 Tracking non-unique solutions through the choice of base states

For the next two examples, we ignore the bulk loading term in the primal
functional:

$I=\int_{0}^{1}\left(\left(u_{x}-1\right)^{2}-1\right)^{2}\,dx;\qquad
u(0)=0;\quad u(1)=1\ $ (43)

The goal is to show that we can recover distinct solutions in a seamless
manner by using the choice of base states which acts to select particular
solutions. As discussed in the preceding Sections 2.2 and 4, an important
factor that determines the stability of a primal solution is the
convexity/concavity of the functional at the obtained solution. We have
plotted the energy density, ($\phi$), its derivative, the stress
($\phi^{\prime}$), and its second derivative, the stiffness
($\phi^{\prime\prime}$) in Fig. 5.

(a)

(b)

(c)

Figure 5: Profile showing the (a) energy density, (b) its derivative (stress),
and (c) second derivative (stiffness). The region between the purple dashed
lines, $0.423\leq e\leq 1.577$, shows the strain range where stiffness is
negative. The black dashed lines mark the zero on the y-axis. The $\star$s
mark points that have the same stress, and are discussed in Sec. 6.3.4.

In this Section, we have obtained weak solutions to the primal E-L system (36)
for which all points in the domain are either in the negative stiffness region
or positive stiffness region using the dual scheme. Different solutions are
obtained depending on the base state used. In one instance, the chosen base
state and the obtained solution are shown in Fig. 6. Another solution for the
given boundary condition is the ‘hat function’ shown in Fig. 7. The base state
provided in this example is of the form, $\bar{e}=1-a$ and $\bar{e}=1+a$.
These two profiles are equally distributed in the different parts of the
domain, as shown in Fig. 6(a) and Fig. 7(a). Displacement boundary conditions
are as in (43). When $0\leq a\leq 0.3$, the former (uniform $e$ profile)
solution is obtained; when $0.9\leq a\leq 5$, the latter (hat function
profile) solution is obtained. Our scheme did not converge for $0.3\leq a\leq
0.9$. Convergence with respect to mesh refinement for the two cases is shown
in Tables. 5 and 6.

(a)

(b)

Figure 6: (a) Base state ($\bar{e}$) and (b) Solution for the primal fields. Number of Elements ($n_{e}$) | $\|\hat{u}^{(n_{e})}-\hat{u}^{(8000)}\|_{1}$ | $\|\hat{e}^{(n_{e})}-\hat{e}^{(8000)}\|_{1}$
---|---|---
100 | $3\times 10^{-4}$ | $3\times 10^{-8}$
2000 | $6\times 10^{-7}$ | $9\times 10^{-11}$
4000 | $1\times 10^{-7}$ | $1\times 10^{-11}$
Table 5: Mesh convergence for computations shown in Fig. 6.

(a)

(b)

Figure 7: (a) Base state $\bar{u}$ and $\bar{e}$ and (b) Solution for the primal fields. This $\hat{u}$ field profile is referred to as hat function because of its shape. Number of Elements ($n_{e}$) | $\|\hat{u}^{(n_{e})}-\hat{u}^{(8000)}\|_{1}$ | $\|\hat{e}^{(n_{e})}-\hat{e}^{(8000)}\|_{1}$
---|---|---
100 | $9\times 10^{-5}$ | $4\times 10^{-7}$
2000 | $2\times 10^{-7}$ | $2\times 10^{-8}$
4000 | $4\times 10^{-8}$ | $1\times 10^{-8}$
Table 6: Mesh convergence for computations shown in Fig. 7.

####  6.3.4 Stability of an equilibrium configuration with negative elastic
stiffness grain boundaries without regularization: elastostatics and
elastodynamics

We consider a 1-d idealization of an elastic body in the shape of a bar as
shown in Fig. 8, comprising three grains separated by grain boundary phases. A
high-angle grain boundary in an ordered medium is a strongly disordered region
of the overall nominally crystalline state. We model the ordered grains as
regions with positive elastic stiffness and the disordered grain boundary
regions as ones with negative elastic stiffness. The material is described
using the energy density profile shown in Fig. 5. As uniform stress over the
domain is one of the solutions to the elastostatic problem, we look for a
configuration with a constant stress field. Three strain ($e$) values (marked
by $\star$ in Fig. 5) resulting in uniform stress over the domain are chosen
in distinct parts of the domain, as shown in Fig. 9(b) and in Table 7. Strain
values in the stable region of the plots (positive stiffness), two of the
three marked, are used to form the grains. The remaining value is used to form
the unstable grain boundaries. The left boundary $(x=0)$ is held fixed and the
displacement b.c. at $x=1$ is determined from integrating the assigned strain
profile along $x$.

Figure 8: Schematic of a bar highlighting green grains separated by red grain boundaries. The left end is fixed and appropriate b.c.s are applied to the right end. Domain | $e^{(t)}$ | $u^{(t)}$
---|---|---
$0\leq x<0.3225$ | $0.115$ | $0.115x$
$0.3225\leq x<0.3325$ | $0.8$ | $0.8(x-0.3225)+0.037$
$0.3325\leq x<0.8275$ | $2.085$ | $0.045+2.085(x-0.3325)$
$0.8275\leq x<0.8875$ | $0.8$ | $1.077+0.8(x-0.8275)$
$0.8875\leq x\leq 1$ | $0.115$ | $1.125+0.115(x-0.8875)$
Table 7: Target solution for the case presented in Fig. 9.

Elastostatics: The base state chosen is as shown in Fig. 9, for the $\bar{e}$
(and consistent $\bar{u}$) and, following standard protocol for solving the
dual problem, we start from an initial guess coincident with it. The applied
boundary condition is $\hat{u}(0)=0$ and $\hat{u}(1)=1.138$. The target
solution is listed in Table 7. The result, Fig. 9, shows that the solution
obtained by the dual scheme has a significant ‘radius’ of its basin of
attraction, as gauged by the deviation of the initial guess (and base state)
from the obtained solution. A convergence study w.r.t mesh refinement is also
performed on the dual solution and presented in Table. 8.

(a)

(b)

Figure 9: (a) Base states $\bar{e}$ field and (b) shows solution for the primal fields. Grain boundaries are marked with red $\star$. Number of Elements | $\|\hat{u}-\hat{u}^{(t)}\|_{1}$ | $\|\hat{e}-\hat{e}^{(t)}\|_{1}$
---|---|---
100 | $2\times 10^{-5}$ | $8\times 10^{-7}$
1600 | $9\times 10^{-8}$ | $6\times 10^{-8}$
8000 | $1\times 10^{-8}$ | $3\times 10^{-9}$
Table 8: Error Reduction in $L^{1}$ norms of displacement and strain with mesh
refinement. These values are for the case corresponding to Fig. 9.

Elastodynamics: We probe the dynamic stability of the equilibrium solution
obtained in the elastostatic case above. As is well-understood, regions
corresponding to negative stiffness, as evaluated from the initial condition
for the dynamic problem, are locally elliptic and hence the primal problem is
ill-posed as an initial value problem and suffers from the ‘Hadamard
instability’ (or lack of continuous dependence w.r.t initial data). Solved by
the dual scheme, however, this problem, which has a reasonable physical
interpretation as mentioned above, is approachable as a stable one, i.e.
thinking of the dual elastodynamic scheme as a mathematical model for a body
with grain boundaries, represented here as containing regions with negative
stiffness, an initial condition corresponding to an elastostatic solution to
(44) remains stable under numerical evolution without further loading. In
contrast, the primal formulation of elastodynamics suffers from severe
numerical instability, resulting in (inevitably present) round-off errors from
computer arithmetic ‘blowing-up’ in finite (very short) time (as well-
understood, and expected).

The theory discussed in Sec. 3 is employed in order to solve the elastodynamic
problem. In 1 space dimension compatibility is not a constraint, and we
consider only two primal variables $e$ and $v$, corresponding to $F$ and
$v_{i}$, respectively. The system of equations solved, therefore, becomes

$\displaystyle\sigma(e)_{x}-\rho_{0}v_{t}$
$\displaystyle=0\qquad\mbox{in}\qquad\Omega\times(0,T)$ (44a) $\displaystyle
e_{t}-v_{x}$ $\displaystyle=0\qquad\mbox{in}\qquad\Omega\times(0,T)$ (44b)
$\displaystyle v(x,0)=v^{(0)}(x);\qquad e(x,0)=e^{(0)}(x)$
$\displaystyle\qquad v(0,t)=v^{(l)}(t);\qquad v(1,t)=v^{(r)}(t).$ (44c)

The base state chosen for $e$ is a shifted cubic and for $v$ a shifted
quadratic:

$H(\hat{v},\hat{e};x,t)=\frac{1}{2}\left(c_{v}(\hat{v}-\bar{v}(x,t))^{2}\right)+\frac{1}{2}\left(c_{e}(\hat{e}-\bar{e}(x,t))^{2}\right)+\frac{1}{3}\left(c_{e}(\hat{e}-\bar{e}(x,t))^{3}\right).$

The DtP mapping equations and the corresponding dual functional, using the
dual (Lagrange multiplier) fields $L$ and $P$, are given by

$\displaystyle\hat{v}$
$\displaystyle=\frac{\rho_{0}L_{t}-P_{x}}{c_{v}}+\bar{v}$ (45) $\displaystyle
c_{e}\frac{\hat{e}-\bar{e}}{|\hat{e}-\bar{e}|}$
$\displaystyle\left(\hat{e}-\bar{e}\right)^{2}+c_{e}\left(\hat{e}-\bar{e}\right)-P_{t}+L_{x}\left(\frac{d\sigma}{de}\Big{|}_{\hat{e}}\right)=0$
$\displaystyle S\,\left[L,P\right]$
$\displaystyle=\,\int_{0}^{T}\int_{0}^{1}\,\,dx\,dt\left[\left(-\rho_{0}L_{t}\hat{v}+\sigma
L_{x}\right)+\left(P_{x}\hat{v}-\hat{e}P_{t}\right)+H\left(\hat{v},\hat{e}\right)\right]$
(46)
$\displaystyle+\int_{0}^{T}\left(Pv^{(l)}\right)dt-\int_{0}^{T}\left(Pv^{(r)}\right)dt-\int_{0}^{1}\left(Pe^{(0)}\right)dx-\int_{0}^{1}\left(\rho_{0}Lv^{(0)}\right)dt.$

The numerical framework for solving the elastostatic problem discussed in Sec.
6.2 extends seamlessly to solve the dual elastodynamics problem. The latter is
a 2-d boundary value problem (in space-time), having $x$ and $t$ as the two
dimensions. All primal boundary and initial conditions are imposed. At any
boundary where boundary/initial conditions on a primal variable is not
available, we apply homogeneous zero Dirichlet boundary condition on the
corresponding dual variable. The residual for this case can again be obtained
by taking the first variation of the dual functional and the second variation
yields the Jacobian. The Newton-Raphson scheme is used to solve this system as
well.

To analyze the stability of the elastostatic solution obtained above in this
section, we provide that solution (see, Table 7) as an initial condition. Zero
velocity boundary conditions are imposed, i.e. $v^{(l)}(t)=v^{(r)}(t)=0$. This
is equivalent to evolving the system without any external influence. We evolve
this case, shown in Fig. 10, using both a primal scheme111The primal scheme
refers to the conventional second-order form of (44) solved by a standard
Galerkin discretization in space and an explicit centered difference scheme in
time with very conservative time stepping. and the dual scheme. It is observed
in Fig. 10(b) that the maximum strain magnitude ($|e|$) becomes unbounded
(‘NaN’ in computer output) within small times when evolved using the primal
scheme. This unbounded behaviour is not observed when the problem is evolved
using the dual scheme (Fig. 10(c)); it stays stable at the initial equilibrium
configuration even after a significant duration of time.

(a)

(b)

(c)

Figure 10: (a) Initial condition ($e$) for both primal and dual scheme. Over
time, the solution to the primal problem diverges (red profile) as shown in
(b), becoming unbounded (‘NaN’) after the next few time steps. Evolution
solved using the dual scheme, shown in (c), remains stable (green profile).
Figures (b) and (c) are plotted for different time intervals.

We next probe the evolution of a small perturbation (i.e. within the typical
elastic limit) of this equilibrium solution in the dual formulation induced by
a perturbation of the initial condition on strain and a corresponding
perturbation in the displacement profile. The initial condition for the
displacement provided w.r.t the equilibrium solution (Fig. 9(b)) is shown in
Fig. 11(a). The perturbations are provided in the stable phases of the domain.
As expected, the initial profile breaks up into two waves of approximately
half the amplitude of the initial displacement perturbation moving in opposite
directions, following the trend of a D’Alembert-like solution from linear
elastodynamics. The wave speed is observed to be higher in the region of
higher stiffness.

(a)

(b)

Figure 11: Evolution of displacement profiles relative to equilibrium: (a)
Initial perturbations, (b) Displacement behavior over time. $u^{(0)}$,
$u^{(eq)}$ represents initial and equilibrium profile of the displacement. $t$
corresponds to time stamp. Green and red color in (a) represents part of the
domain in stable and unstable phase, respectively.

##  7 Conclusion

We have demonstrated a methodology for developing a (family of) dual
variational principle(s) for nonlinear elastostatics and elastodynamics with
the property that its Euler-Lagrange equations are the governing PDE of
nonlinear elasticity with the primal fields represented in terms of a mapping
of the dual fields and their space-time gradient, pointwise. The dual
variational principles are convex by design, regardless of the properties of
the primal problem. Existence of minimizers of such a dual variational
principle is rigorously shown, and the notions of variational dual solutions
and dual solution to the PDE of nonlinear elastostatics are defined in terms
of the dual functional. Computed examples of illustrative cases corresponding
to noncovex elastostatics and elastodynamics, including a stable dual
calculation of an ill-posed, primal elastodynamic problem is also
demonstrated.

The computed examples, especially in Sec. 6.3.4, raise a philosophical
question about the optimal, minimalist, mathematical model that should be
adopted for modeling physically observed softening behavior in macroscopic
mechanical response. This work seems to suggest that existence of dual
solutions may not be a show-stopper, with the means of attaining uniqueness in
the model transferred to the choice of special types of auxiliary functions
$H$ employing physically meaningful base states which, nevertheless, need not
be solutions to singularly-perturbed higher order PDE when physical
justification for such is not available, but include experimentally observed
input on physically relevant length scales. A crucial test in this regard
would be to study the details of wave ‘propagation’ through the non-hyperbolic
regions in our model, and check how well such behavior squares with
experiment.

## Acknowledgment

We thank Robert V. Kohn for helpful discussion and suggestions. This work was
partially done while AA was on sabbatical leave in 2023 at a) the Max Planck
Institute for Mathematics in the Sciences in Leipzig, and b) the Hausdorff
Institute for Mathematics at the University of Bonn funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC-2047/1 – 390685813, as part of the Trimester Program
on Mathematics for Complex Materials. The support and hospitality of both
institutions are acknowledged. The collaboration between JG and AA on this
paper was initiated at the trimester program in Bonn. It is also a pleasure to
acknowledge discussions at these locations on the subject matter of this paper
with John M. Ball, Lukas Koch, and Angkana Rüland. AA was supported by a
Simons Pivot Fellowship grant # 983171. SS was supported by the NSF OIA-DMR
grant # 2021019 and the Center for Extreme Events in Structurally Evolving
Materials, ARL contract # W911NF2320073.

## References

  * [1] Amit Acharya.  An action for nonlinear dislocation dynamics.  Journal of the Mechanics and Physics of Solids, 161:104811, 2022\.
  * [2] Amit Acharya.  Variational principles for nonlinear PDE systems via duality.  Quarterly of Applied Mathematics, LXXXI:127–140, 2023, Article electronically published on September 26, 2022.
  * [3] Amit Acharya.  A hidden convexity in continuum mechanics, with application to classical, continuous-time, rate-(in) dependent plasticity.  arXiv preprint arXiv:2310.03201, 2023.
  * [4] Amit Acharya.  A dual variational principle for nonlinear dislocation dynamics.  Journal of Elasticity, https://doi.org/10.1007/s10659-023-09998-5, 2023.
  * Aro [23] Abhishek Arora.  A study of nonlinear deformations and defects in the actuation of soft membranes, rupture dynamics, and mesoscale plasticity.  Phd thesis, Carnegie Mellon University, September 2023.  Available at https://www.proquest.com/docview/2869411693?accountid=9902.
  * Bal [76] J. M Ball.  Convexity conditions and existence theorems in nonlinear elasticity.  Archive for rational mechanics and Analysis, 63:337–403, 1976.
  * Bal [02] J. M. Ball.  Some Open Problems in Elasticity, pages 3–59.  Springer New York, New York, NY, 2002.
  * Bha [03] Kaushik Bhattacharya.  Microstructure of martensite: why it forms and how it gives rise to the shape-memory effect, volume 2.  Oxford University Press, 2003.
  * BM [79] Victor L. Berdichevskii and V. A. Misiura.  On a dual variational principle in geometrically nonlinear elasticity theory.  Journal of Applied Mathematics and Mechanics, PMM, 43(2):343–352, 1979.
  * Bra [02] A. Braides.  $\Gamma$-convergence for beginners, volume 22.  Clarendon Press, 2002.
  * Bre [18] Y. Brenier.  The initial value problem for the Euler equations of incompressible fluids viewed as a concave maximization problem.  Communications in Mathematical Physics, 364(2):579–605, 2018.
  * Bre [20] Y. Brenier.  Examples of hidden convexity in nonlinear PDEs.  https://hal.science/hal-02928398/document, 2020.
  * CKPdN [15] M. Carozza, J. Kristensen, and A. Passarelli di Napoli.  On the validity of the Euler-Lagrange system.  Communications on Pure and Applied Analysis, 14(1):51–62, 2015\.
  * DL [99] R. Dautray and J.-L. Lions.  Mathematical analysis and numerical methods for science and technology: volume 3 Spectral Theory and Applications, volume 3.  Springer Science & Business Media, 1999.
  * dV [72] B. F. de Veubeke.  A new variational principle for finite elastic displacements.  International Journal of Engineering Science, 10(9):745–763, 1972\.
  * [16] Uditnarayan Kouskiya and Amit Acharya.  Inviscid Burgers as a degenerate elliptic problem.  In Preparation, 2023.
  * [17] Uditnarayan Kouskiya and Amit Acharya.  Hidden convexity in the heat, linear transport, and Euler’s rigid body equations: A computational approach.  to appear in Quarterly of Applied Mathematics, https://arxiv.org/abs/2304.09418, 2023.
  * KK [22] Lukas Koch and J. Kristensen.  On the validity of the Euler-Lagrange system without growth assumptions, 2022.  arXiv:2203.00333.
  * Koi [73] W. T. Koiter.  On the principle of stationary complementary energy in the nonlinear theory of elasticity.  SIAM Journal on Applied Mathematics, 25(3):424–434, 1973.
  * KTW [03] Robin J. Knops, C. Trimarco, and H. T. Williams.  Uniqueness and complementary energy in nonlinear elastostatics.  Meccanica, 38:519–534, 2003.
  * LDR [95] H. Le Dret and A. Raoult.  The quasiconvex envelope of the Saint Venant-Kirchhoff stored energy functional, pages 171–180.  De Gruyter, Berlin, New York, 1995.
  * Leo [17] G. Leoni.  A first course in Sobolev spaces.  American Mathematical Soc., 2017.
  * Mor [52] C. B. Morrey.  Quasi-convexity and the lower semicontinuity of multiple integrals.  Pacific Journal of Mathematics, 2:25–53, 1952.
  * MQY [94] S. Müller, Tang Qi, and B. S. Yan.  On a new class of elastic deformations not allowing for cavitation.  Annales de l’I.H.P. Analyse non linéaire, 11(2):217–243, 1994\.
  * Zub [70] L. M. Zubov.  The stationary principle of complementary work in nonlinear theory of elasticity.  Journal of Applied Mathematics and Mechanics, PMM, 34(2):228–232, 1970.