A Hidden Convexity of Nonlinear Elasticity

A technique for developing convex dual variational principles for the governing PDE of non-linear 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 solu-tion 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.


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 and the equations of equilibrium for the primal model of elasticity invoked.Even in the case of hyperelasticity, if the primal energy functional is not lower-semicontinuous then classical variational techniques cannot be used to prove existence of minimizers.In this case, the classical approach is to consider the relaxation of the energy functional.However, minimizers of the relaxed energy functional typically do not satisfy the original PDE system.
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 [Zub70,Koi73,dV72,BM79,KTW03]; 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 [Ach22a, Sec.6.1]- [Ach23,Ach24], and connects to the ideas of 'hidden convexity in nonlinear PDE' advanced by Y. Brenier [Bre20,Bre18].
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.

Dual nonlinear Elastostatics
Consider the problem of nonlinear elastostatics with boundary conditions of place and dead loads (for simplicity): P ij n j = p i on ∂Ω p ; i on ∂Ω y , ∂Ω = ∂Ω p ∪ Ω 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 Ω (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 F T F for S being any arbitrary symmetric tensor valued function of a symmetric tensor; this allows frameindifference 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 λ and µ, respectively, integrating by parts, and imposing the prescribed b.c.s: (2)

Define
U := (y, F ); D := (λ, µ); D := (∇λ, λ, ∇µ, µ); and require the choice of the potential H to be such that it facilitates the existence of a function When such a 'change of variables' mapping, U H = y (H) , F (H) , exists, defining the dual functional as i µ ij n j da with λ specified on ∂Ω\∂Ω p and µ specified on ∂Ω\∂Ω y , and noting (4), the first variation of S H (about a state D in the direction δD, the latter constrained by δλ = 0 on ∂Ω\∂Ω p and δµ = 0 on ∂Ω\∂Ω y ), is given by Noting, now, that L H is necessarily affine in 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 (D, x).
It is this simple idea that we exploit to develop variational principles for (Cauchy) Elasticity.

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 → ȳ(x); x → ∇ȳ(x) =: F (x), we assume that the nominal stress fields, x → P (x), for all candidate pairs (y, F ) that belong to the domain of the functional S are well-represented by a quadratic expansion around (ȳ, F ): i.e.
where 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 P , A, B may be interpreted as those of a truncated Taylor expansion of the nominal stress response function about 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: where c y , c F ≫ 1 are real-valued scalars.Then, the implementation of (4) in this context means algebraically solving for y (Q) , F (Q) in terms of (x, ∇λ, ∇µ, µ) from the dual-to-primal (DtP) mapping equations which can be rearranged as the following relations (suppressing the explicit x dependence of the base state for brevity) where Then 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): and using (6) again, a dual functional for nonlinear (Cauchy) elastostatics is given by

Formal second variation and stability of dual solutions
From the DtP equations (6), a primal pair y (Q) , F (Q) corresponding to the dual fields is given by (x → ȳ(x), x → F (x)).
We are now interested in obtaining the second variation of the dual functional S Q about the dual state D, along sufficiently smooth perturbations (x → δD(x) , x → dD(x)) for c F , c y ≫ 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 the first variation of S Q , at D and in the direction δD, is given by with the evaluation Our interest is in examining the second variation of S Q , at the dual state D and in arbitrarily fixed directions δ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 D belongs to this class), and then find that where the coefficients of c −2 F , c −3 F are finite (the negative powers on c F arise from repeated appearances of ∂K ∂∇λ in the integrands comprising these terms), so that for c F ≫ 1 the second variation is controlled by the first term on the right hand side of (9).Furthermore, for any dual state for which ∇λ is bounded on the domain, c F can be chosen large enough so that K and K −1 are positive-definite.
Then, for c F , c y ≫ 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 ' ∂L ∂U = 0.' The second variation at the state D is given, regardless of the magnitude of c F , by 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 (ȳ, F ) is chosen as a solution to the primal problem -even in a situation where the physical moduli field x → A| 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 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.

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 where [0, T ] is an interval of time, P satisfies the same conditions as in Sec. 2 and ρ 0 is a given timeindependent, positive, scalar-valued function on Ω representing the mass density in the reference configuration.
Following the steps of Sec. 2, we define the pre-dual functional 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 Ω × (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 B := (∂Ω × (0, T )) ∪ (Ω × {0, T }).Defining now and requiring the choice of the potential H to be such that it facilitates the existence of a function it is observed that defining the dual functional as with (arbitrarily) specified λ on B\( (∂Ω p × (0, T )) ∪ (Ω × {0}) ) and and noting (12), the first variation of S H (about a state D in the direction δD, with δλ and δγ constrained to vanish on parts of the boundary of B on which λ and γ are specified, respectively), is given by 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 B) -of the dual functional (13) are exactly the system (10) with U substituted by U H (D, x, t).
As argued in [Ach22b, Sec.7] and demonstrated in [KA23], 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: We make an arbitrary choice of a base state and adopt the following shifted quadratic for the potential H: Next, we form the Lagrangian corresponding to the system (10) and the potential (15): The DtP mapping corresponding to ( 16)-( 14)-( 15) is then given by the relations where K is defined in (6c).Utilizing these ingredients, the dual action, corresponding to 'quadratically' nonlinear elastodynamics (10) and ( 15) is given by (again dropping the superscript (Q) for brevity)

Local degenerate ellipticity of dual nonlinear Elasticity -statics and dynamics
This treatment is a special case of that in [Ach24, Sec.3] and is included here for a self-contained account.
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, • • • , N }, indexing the components of the N × 1 array of primal variables, U , with, possibly, a conversion to first-order form as necessary.Now consider the system of primal PDE where upper-case Greek indices index the number of equations involved, after conversion to firstorder form when needed.For example, consider the three conservation laws for nonlinear elastodynamics (10a).The conversion to first-order form (so that the Lagrangian 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 where A, B are constant matrices (with 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 → ∂ 2 F Γ ∂U P ∂U R (U ) and U → ∂ 2 G Γ ∂U P ∂U R (U ) are bounded functions on their domains.
Let D be the N * × 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 Then the corresponding DtP mapping, obtained by 'solving ∂L ∂U = 0 for U in terms of (∇D, D, Ū ),' is given by the implicit equation By the fundamental property of the dual scheme, the dual E-L equation is then given by (where we have dropped the superscript (Q) for notational convenience), whose ellipticity is governed by the term .
From (20) we have and so which is positive semi-definite on the space of N * × 3 (or N * × 4) matrices.This establishes the degenerate ellipticity of the dual system at the state x → D(x) = 0.
To examine the ellipticity-related properties of the system in a bounded neighborhood, say N , of (D = 0, ∇D = 0) ∈ R N * × R N * × ᾱ, ᾱ = 3, 4, we define and note that where M −1 exists and is positive definite by the boundedness of N and the second derivatives of the functions 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 N is now defined as the positive semi-definiteness of the matrix A on the space R N * × ᾱ of matrices, and this in turn is governed by the matrix By the positive definiteness of the matrix [M P S ] in the neighborhood N , it follows that 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 A on the space of tensor products from R N * ⊗ R ᾱ, and not directly comparable to the strong-ellipticity condition, since it is weaker than the latter when restricted to the space R N * ⊗ R ᾱ but simultaneously requiring semi-definiteness on the larger space of R N * × ᾱ.Also of note is that degenerate ellipticity does not preclude the failure of ellipticity characterized by the condition det[A Γ αΠµ n α n µ ] ̸ = 0 for all unit direction n ∈ R ᾱ, ᾱ = 3 or 4, thus allowing for weak (gradient) discontinuities of weak solutions x → Ḋ(x) of the linearization of ( 21).If a solution of the primal system is close to the base state Ū , 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 Ū , for instance by an iterative scheme starting from a guess (D = 0, U = Ū ).
Our experience (Sec.6 of this work, [KA23, KA24, Aro23]) shows that this observation is of great practical relevance in using the dual scheme, and we consistently exploit it in all our computational approximations.

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 for fixed values of (D, x).If for all D the maximizer is unique then we call the mapping (D(x), x) → U H (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 = (λ, µ) as D → (x → U H (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) but the definition of SH 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 SH that we adopt in this section.We note that the functional ŜH is affine in the functions λ, ∇λ, div(µ) and µ.In particular, if formulated in appropriate spaces it will be continuous with respect to weak convergence of these quantities.This implies that SH as the sup over these weakly continuous functions is weakly lower semicontinuous.This renders the dual problem SH 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 SH 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 SH ) 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: 2. Prove up to a non-relabeled subsequence (weak) compactness for the minimizing sequence (λ k , µ k ) ⇀ (λ, µ).
Note that by the argumentation above the dual functional SH is (at least formally) lower semicontinuous 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 ŜH is affine in λ and µ the function SH 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 SH .
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 (λ, µ) is a variational dual solution to the equation (1) if it minimizes the dual functional SH .Moreover, we say that a pair (λ, µ) is a dual solution to (1) if the dual functional SH is differentiable at (λ, µ), the pair (λ, µ) satisfies the weak form of the corresponding Euler-Lagrange equation for SH and there exists a smooth dual-to-primal mapping U H corresponding to SH .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 SH allows for nontrivial minimizers.For example, for H = 0 it is easy to check that for the examples discussed below it follows SH [D] = +∞ except for D = (0, 0), where SH [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 SH 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 +∞ 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 Ω y .Much research in this direction was inspired by the seminal works on functionals with quasiconvex [Mor52] or polyconvex energy densities [Bal76,MQY94].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 [Bal02] and the references therein (in the case of a convex integrand see, for example, [KK22,CKPdN15] 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 SH 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 SH 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 SH can also be made rigorous and hence existence of minimizers of SH 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 Ω ⊆ R d is open, connected and bounded with Lipschitz boundary.

The Saint Venant-Kirchhoff model
We consider the Saint Venant-Kirchhoff energy where C = F T F is the right Cauchy-Green tensor, E = 1 2 (C − Id) is the Green-Lagrange strain tensor, the shear modulus G > 0 and G + d 2 L > 0. Correspondingly, we set For y (b) ∈ W 1,4 (Ω; R d ) and p (b) ∈ L 4/3 (∂Ω; R d ) we look for weak solutions to the equations (1) which take the form Precisely, we look for y ∈ W 1,4 (Ω; R d ) and F ∈ L 4 (Ω; R d×d ) such that for all λ ∈ W 1,4 (Ω; R d ) with λ = 0 on ∂Ω y and µ ∈ L 4/3 (Ω; R d×d ) with div(µ) ∈ L 4/3 (Ω; R d ) and µn = 0 on ∂Ω p that Ω div(µ) • y + µ : Note that the boundary terms are well-defined but have to be understood in the sense of traces i.e., for p ∈ [1, ∞] there are continuous trace operators [Leo17,Theorem 18.40] and [DL99, 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 [LDR95]).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 and consider the functional ), µn = 0 on ∂Ω p and λ = 0 on ∂Ω y } → R defined as We will first show that I agrees with the dual functional SH .
Clearly, it holds that sup We will now prove the reverse inequality.Let O = {O ⊆ Ω : O open} and consider the countable family of functions q n χ En and y E n is a cuboid with rational side lengths and corners in Moreover, it holds by definition for each (y, F ) ∈ A and O ∈ O that Then it follows from [Bra02, Lemma 15.2] that This proves the reverse inequality.
As explained above, in order to prove existence of minimizers for SH it suffices by the direct method of the Calculus of Variations to prove that SH 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, Proof.By optimizing in y it follows that g(A, a, B) where we used that where c > 0 is a universal constant.
Consequently, we find for a constant c > 0 (that depends on G and L and might change from line to line) that Remark 5.1.Similarly, one can show that it also holds which shows that the dual functional I = SH is well-defined and its minimization is non-trivial.
Let us now show the coercivity of the dual functional I = SH .

Incompressible Neo-Hookean model in two dimensions
Let us now consider the stress in an incompressible Neo-Hookean model where G > 0 is the shear modulus, det(F ) = 1 and p is the pressure which is to be determined.For simplicity, let G = 1.Moreover, set P (F, p) = F − p cof(F ).
Note that in two dimensions P depends linearly on F and for F with det(F ) = 1 it holds P (F, p) = P (F, p).Then for y (b) ∈ W 1,4 (Ω; R d ) and p (b) ∈ L 4/3 (∂Ω; R d ) we look for weak solutions to Namely, we seek to find y ∈ L 2 (Ω; R 2 ) and 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 det(F ) = 1 is treated through an extra dual variable ϑ, whereas the pressure p is an additional primal variable.Following the same approach, we define the function g : where H(y, p, F ) = 1 2 y 2 + 1 2 |F | 2 + 1 4 p 4 .Correspondingly, we define the functional By the same argument as for the Saint Venant -Kirchhoff material in Proposition 5.1 it holds for all (µ, λ, ϑ) = sup The right hand side of the above equation is the analogue of SH including the extra variable ϑ.
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.
Proof.First note that by optimizing in y we find g(A, a, B, s) = sup We will now only work with the function h.We observe that it holds for all F ∈ R 2×2 that F : cof(F ) = 2 det(F ).In particular, | det(F )| ≤ 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 The case s < −1 can be treated similarly.This shows 1. Next, let us consider |s| ≤ 1.First we assume that 1 4 √ 8 |B| 2 ≤ |A + B|.Then set F = (A + B)/2 and p = 0 to obtain Combining ( 32), ( 33) and (34) yields which is 2. It remains to show 3. We estimate h(A, B, s) ≤ sup 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 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 ỹ ∈ L 2 (Ω; R 2 ) and F ∈ L 2 (Ω; R 2×2 ) and note for H(F, y, p; x) = H(F − F (x), y − ỹ(x), p) that S H (λ, µ, ϑ) This functional has the same coercivity and lower-semicontinuity properties as SH .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 ϑ = ±1 the function g is not differentiable.Making good choices for F and ỹ could enable us to guarantee that a minimizer (λ, µ, ϑ) of S H only takes values such that g is actually differentiable in the points (µ − F , div(µ) − ỹ, ∇λ, ϑ).
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.
It is clear from (37) and ( 38) that û = ū and ê = ē is a value of the DtP mapping for D = 0.
We substitute the DtP mapping (38) in the functional S to define the dual functional The requirement (37) ensures that the dual E-L system of (39) is simply the primal system (36), with the replacements u → û, 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.

Finite element algorithm
Denoting the dual fields as 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: A Dirichlet boundary condition is applied only on µ.A Neumann boundary condition is applied on λ, 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 (•) , where (•) 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 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 (D) given by and its variation in the direction dD, the Jacobian that are used in a Newton Raphson scheme.Here, R A i , J AB ij are given by the expressions and The Newton-Raphson scheme involves solving the following matrix equation for the k th correction to an initial guess and this is continued until convergence of |R(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 (û, ê) with physical meaning, and infer initial guesses for the dual fields consistent with it as a solution of the DtP mapping.Choosing initial guesses û(i) and ê(i) the initial guess for the dual fields, D B i = 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: Algorithm Initialization: Set D B(0) j = 0. Choose a value for tol.
1. Newton's Method is used to evaluate the results.Set k = 1 < tol then go to step 2 else do k = k + 1 and go to step 1i 2. Evaluate the values of û and ê at the nodes using the L 2 projection.These nodal values establish the solution.
• tol represents the tolerance for the iterative solve of the Newton Raphson iterations; residuals with |R(D)| < tol are accepted as having approximately satisfied the corresponding equations.We use the definition |R(D)| := max A,i R A i (D), where (A, i) ranges over all dual nodal degrees of freedom on which a Dirichlet boundary conditions are not specified.
• 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: • In all our calculations we seek weak solutions of the E-L system so that jump conditions (like traction continuity) are automatically imposed.

Results
The algorithm of Sec.6.2 is utilized to solve key problems demonstrating the theory.It is wellunderstood 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 [Bha03]).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 = β are solutions provided β = α = α * .

Stress-free solution
We consider the primal energy functional and compute solutions to its Euler-Lagrange equation ( 36) with α = α * = 1.A stress-free equilibrium solution of the primal problem is given by 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. 1a and Fig. 1b.The base state provided for an ē profile is the derivative of the corresponding ū profile.The ū 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. 1a.These values are of the same order for the other sub-case (Fig. 1b) as well.Table 2: Error reduction in L 1 norms of displacement and strain with mesh refinement.These values are for the case corresponding to Fig. 1a.

Stressed solutions with spatially homogeneous and inhomogeneous strain configurations
We now consider the primal functional with e (t) = 0.5, u (t) = 0.5x, as the exact target solution.We choose a base state near this solution as shown in Fig. 2a and start from an initial guess coincident with the base state.The base state provided for an ē profile is the derivative of the corresponding ū profile.The ū profiles satisfy the displacement b.c.s of the target solution.The dual scheme is able to recover the target solution, presented in Fig. 2b.The L 1 displacement and strain errors relative to the target solution are summarized in  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 Clearly, a homogeneous strain solution is not possible in this case due to the mismatch α ̸ = α * 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. 4b.The difference between the base state chosen for this case and the solution can be observed in Fig. 4a.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 û(ne) corresponds to the discrete displacement solution obtained using n e number of elements and similarly for the strain, ê.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 α = 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 ) = αx 0 and u(0) = α * = 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, σ, 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 ≤ x ≤ x 0 would consist of a typical 'staircase' profile of slopes 0, 2 along the backbone piece of slope α, as shown in Fig. 3.However, due to the presence of the forcing term (u − α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 ̸ = 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 1 2 backbone up to x 0 = 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 = 2 3 , with constant slopes of 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.

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: 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, (ϕ), its derivative, the stress (ϕ ′ ), and its second derivative, the stiffness (ϕ ′′ ) in Fig. 5.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, ē = 1 − a and ē = 1 + a.These two profiles are equally distributed in the different parts of the domain, as shown in Fig. 6a and Fig. 7a.Displacement boundary conditions are as in (43).When 0 ≤ a ≤ 0.3, the former (uniform e profile) solution is obtained; when 0.9 ≤ a ≤ 5, the latter (hat function profile) solution is obtained.Our scheme did not converge for 0.3 ≤ a ≤ 0.9.Convergence with respect to mesh refinement for the two cases is shown in Tables.5 and 6.

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 ⋆ in Fig. 5) resulting in uniform stress over the domain are chosen in distinct parts of the domain, as shown in Fig. 9b 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.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 σ(e) x − ρ 0 v t = 0 in Ω × (0, T ) (44a) The base state chosen for e is a shifted cubic and for v a shifted quadratic: 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 scheme1 and the dual scheme.It is observed in Fig. 10b 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. 10c); it stays stable at the initial equilibrium configuration even after a significant duration of time.

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.

)
Note that A is countable.Moreover, it holds for all (a, A, B) ∈ R d × R d×d × R d×d and x ∈ Ω sup y,F ∈D a • y(x) + A : F (x) + P (F (x)) : B − H(y(x), F (x)) = g(A, a, B).For O ∈ O let us now define F(O) = sup y∈L 4 ,F ∈L 4 O y • div(µ) + µ : F + P (F ) : ∇λ − H(y, F ) dx.It can be shown that for O, O ′ ∈ O such that dist(O, O ′ ) > 0 it follows that

Figure 1 :
Figure 1: (a) and (b) show the initial guesses, ē, for two cases and (c) shows the solution for the primal fields.Fig. (c) contains both û and ê 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.

Figure 2 :
Figure 2: (a) Base state ē and (b) Solution for the primal fields.

Figure 3 :Figure 4 :
Figure 3: The displacement profile for the case x 0 = 2 3 .The staircase profile (green) rides on a line of slope 1 2 (blue) until x 0 , followed by an orange line of slope 2.

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

Figure 8 :Figure 9 :
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.

Figure 11 :
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.

Table 4 :
Mesh convergence for computations shown in Fig.4.

Table 6 :
Mesh convergence for computations shown in Fig.7.