A new perspective in linear Cauchy Elasticity: variational minimum principles for statics, dynamics, and heterogeneous materials

Abstract

A variational minimum principle for linear elastodynamics of a possibly heterogeneous material without a stored energy function is developed. It involves a change of variables to dual fields, and results in a degenerate elliptic Euler-Lagrange system, even when the primal elastodynamics is hyperbolic. Uniqueness assertions for the dual dynamic and static problems and implications of the degenerate ellipticity are sketched. Some implications pertaining to heterogeneous materials and ones with indefinite elastic moduli are discussed.

1 Introduction

This work describes an application to Linear Cauchy Elasticity of a technique [1, 3, 2] for associating a family of variational minimum (i.e., convex) principles to a system of equations with boundary and initial conditions when appropriate. The given ‘primal’ equations may not have a variational structure in terms of the primal fields or variables they are originally posed in, e.g., elasticity without a strain energy function. The procedure involves an adapted change of variables to a new set of ‘dual’ fields in terms of which a functional can be posed whose Euler-Lagrange equations become the primal set of equations, interpreted through the change of variables mapping involved. This approach has been further developed and demonstrated through computation in [8, 9, Ach6, 15, 10, 16], with rigorous results developed in [15, AGS24, AG25, 10, 17, ASZ24]. Here, we work out the details of the scheme in the context of the governing system being that of linear elastodynamics of a generally heterogeneous material possibly without a strain energy function. As is well understood, accounting for initial conditions of an initial value problem is problematic in a variational principle - including Hamilton’s celebrated principle - see, e.g., [5] for a discussion; Hamilton’s principle also has a problem with causality due to the need for an acausal final-time boundary condition on position (cf. [19], [4, 7]). The way out of this difficulty in linear problems has been to resort to working in the Laplace transform domain [5, 19], or at ‘fixed frequency’ [11], i.e., for elastodynamic fields satisfying a special ansatz for its time-dependent part, assumed separable from its space-dependent part; also see [7] in the context of Hamiltonian particle mechanics using an adaptation of the Schwinger-Keldysh formalism which involves a doubling of variables and a backward-in-time path. The approach adopted here works in real space-time domains without involving convolutions, assumptions of separability, or a doubling of variables in which the original system is posed, and works as well for nonlinear and dissipative time-dependent equations as shown in [9, 16, 17, 15].

The work [15] develops the dual approach for nonlinear elasticity, focusing primarily on elastostatics; the form of the dual functional is developed for quadratic nonlinearity in stress. The emphasis in this work is to explicitly work out the case of linear elasticity in full detail, which is also of interest to a much broader class of researchers interested in its theoretical and/or practical aspects. As examples, the specialization to linear elastodynamics greatly facilitates the explorations of the type provided in Sections 3 and 4. For comparison to other variational techniques for linear elastodynamics in the literature, to our knowledge this work provides the first convex variational principle for linear elastodynamics (with a strain energy function and without) in the time domain without approximation, recovering field equations, boundary conditions, and initial conditions without any loss of causality. In situations when a stored energy function exists, the recovery of initial conditions, a causal formulation, and the convexity of the variational principle for the whole range of positive-definite to indefinite stored energy function are some of the main contributions of this work.

An outline of the paper is as follows: in Sec. 2 the variational principle for elastodynamics is developed. In Sec. 3 a uniqueness result for the Euler-Lagrange system of the variational principle is studied; some discussion of its ellipticity properties (in space-time domains) is also provided. In Sec. 4 implications of the scheme for heterogeneous/composite materials are discussed. Some concluding remarks are recorded in Sec. 5.

2 A convex variational principle for linear elastodynamics

Given a 3-dimensional domain $\Omega$ with boundary $\partial\Omega=\overline{\partial\Omega_{s}\cup\partial\Omega_{u}}$ the closure of a union of disjoint sets, and an interval of time $[0,T]$, consider the system of linear elastodynamics whose given elastic moduli and density fields, $(x,t)\mapsto C(x),\rho(x)$ can vary in space. $C$ has minor symmetries and may or may not have major symmetries; we will use the notation $(C^{T})_{ijkl}=C_{klij}$. The fields $u$, $v$ will denote the displacement and velocity, and $E$ and $W$ the symmetric and antisymmetric parts of the displacement gradient. When using direct notation a “$\ \cdot\$” will denote the inner product of two vectors and a “$\ :\$” that of two second-order tensors. For fourth-order tensors, $A,B,C$, $(AB)_{ijkl}=A_{ijrm}B_{rmkl}$ and $(ABC)_{ijkl}=A_{ijrm}B_{rmps}C_{pskl}$. For this work, it will be convenient to work with a non-dimensional representation of the equations of linear elasticity. Using $E^{*}$ as a physical scale for stress, $L^{*}$ as some physical scale for length, and $T^{*}$ as a scale for time, the equations of linear elastodynamics can be written as

$$\begin{cases}&\partial_{\tilde{t}}\tilde{u}_{i}-\tilde{v}_{i}=0\\ &\tilde{\rho}\partial_{\tilde{t}}\tilde{v}_{i}-\partial_{\tilde{x}_{j}}(\tilde{C}_{ijkl}E_{kl})-\tilde{b}_{i}=0\\ &\partial_{\tilde{x}_{l}}\tilde{u}_{k}-E_{kl}-W_{kl}=0\\ \end{cases}\qquad\mbox{ on }\Omega\times(0,T)$$ (1)

$$\tilde{u}_{i}=\tilde{u}_{i}^{(b)}\mbox{ on }\partial\Omega_{u}\times(0,T);\qquad\qquad(\tilde{C}_{ijkl}E_{kl})n_{j}=\tilde{t}_{i}\mbox{ on }\partial\Omega_{s}\times(0,T)$$

$$\tilde{u}_{i}=\tilde{u}_{i}^{(0)}\mbox{ on }\Omega\times\{0\};\qquad\qquad\tilde{v}_{i}=\tilde{v}_{i}^{(0)}\mbox{ on }\Omega\times\{0\}.$$

where

$$\tilde{\rho}=\frac{\rho L^{*2}}{E^{*}T^{*2}},$$

$$\tilde{x}_{i}=\frac{x_{i}}{L^{*}},\quad\tilde{t}=\frac{t}{T^{*}},\quad\tilde{u}_{i}=\frac{u_{i}}{L^{*}},\quad\tilde{v}_{i}=\frac{v_{i}T^{*}}{L^{*}},\quad\tilde{C}_{ijkl}=\frac{C_{ijkl}}{E^{*}},$$

$$\tilde{u}^{(0)}_{i}=\frac{u_{i}^{(0)}}{L^{*}},\quad\tilde{v}^{(0)}_{i}=\frac{v_{i}^{(0)}T^{*}}{L^{*}},\quad\tilde{u}^{(b)}_{i}=\frac{u_{i}^{(b)}}{L^{*}},\quad\tilde{t}_{i}=\frac{t_{i}}{E^{*}},\quad\tilde{b}_{i}=\frac{b_{i}L^{*}}{E^{*}}.$$

Here, $b$ the body force density, $u^{(b)},t$ the displacement and traction boundary conditions, and $u^{(0)},v^{(0)}$ the displacement and velocity initial conditions. Henceforth, we work with the system (1) dropping all tildes for convenience:

$$\begin{cases}&\partial_{t}u_{i}-v_{i}=0\\ &\rho\partial_{t}v_{i}-\partial_{j}(C_{ijkl}E_{kl})-b_{i}=0\\ &\partial_{l}u_{k}-E_{kl}-W_{kl}=0\\ \end{cases}\qquad\mbox{ on }\Omega\times(0,T)$$ (2)

$$u_{i}=u_{i}^{(b)}\mbox{ on }\partial\Omega_{u}\times(0,T);\qquad\qquad(C_{ijkl}E_{kl})n_{j}=t_{i}\mbox{ on }\partial\Omega_{s}\times(0,T)$$

$$u_{i}=u_{i}^{(0)}\mbox{ on }\Omega\times\{0\};\qquad\qquad v_{i}=v_{i}^{(0)}\mbox{ on }\Omega\times\{0\}.$$

Using the notation $U:=(u,v,E,W)$, and $\bar{U}:=(\bar{u},\bar{v},\bar{E},\bar{W})$ where $(x,t)\mapsto\bar{U}(x,t)$ are arbitrarily specified functions referred to as base states (in space-time domains), we introduce an auxiliary potential

$$H(U,\bar{U}) =\frac{1}{2}\bigg((E_{ij}-\bar{E}_{ij})A_{ijkl}(E_{kl}-\bar{E}_{kl})\ +\ (W_{kl}-\bar{W}_{kl})(W_{kl}-\bar{W}_{kl})$$ (3)

$$+(u_{i}-\bar{u}_{i})(u_{i}-\bar{u}_{i})\ +\ \mathsf{r}(v_{i}-\bar{v}_{i})(v_{i}-\bar{v}_{i})\bigg)$$

for the system of linear elastodynamics, where $A,\mathsf{r}$ are nondimensional fields, with freedom to be designed as needed. Due to the quadratic form through which it enters, $A$ is assumed to have major symmetry without loss of generality.

Consider now a pre-dual functional

$$\widehat{S}[U,\gamma,\lambda,\mu] =\int_{0}^{T}\int_{\Omega}\Big(-u_{i}\partial_{t}\gamma_{i}-\gamma_{i}v_{i}\quad-\quad v_{i}\partial_{t}(\rho\lambda_{i})+E_{kl}C_{ijkl}\partial_{j}\lambda_{i}$$

$$\qquad\qquad\quad\quad-\quad u_{k}\partial_{l}\mu_{kl}-\mu_{kl}(E_{kl}+W_{kl})\quad-\quad H(U,\bar{U})\Big)\,dxdt$$

$$\quad-\int_{0}^{T}\int_{\Omega}\lambda_{i}b_{i}\,dxdt-\int_{\Omega}\Big(\gamma_{i}(x,0)u_{i}^{(0)}(x)+\rho(x)\lambda_{i}(x,0)v_{i}^{(0)}(x)\Big)\,dx$$ (4)

$$\quad-\int_{0}^{T}\int_{\partial\Omega_{s}}\lambda_{i}t_{i}\,dadt+\int_{0}^{T}\int_{\partial\Omega_{u}}\mu_{kl}u^{(b)}_{k}n_{l}\,dadt$$

with

$$\gamma_{i}=0\mbox{ on }\Omega\times\{T\}; \lambda_{i}=0\mbox{ on }\Omega\times\{T\};$$ (5)

$$\lambda_{i}=0\mbox{ on }\partial\Omega_{u}\times(0,T); \mu_{kl}n_{l}=0\mbox{ on }\partial\Omega_{s}\times(0,T).$$

The fields $D:=(\gamma,\lambda,\mu)$ are referred to as dual fields. The values (right hand-sides) of the dual (Dirichlet) space-time boundary conditions (5) can be arbitrarily specified, without loss of generality (note that in such cases no extra terms are, nevertheless, included in the pre-dual functional which is free to design.)

We refer to the first space-time bulk integrand in (4) as the Lagrangian

$$\mathcal{L}(U,\mathcal{D},\bar{U})\qquad\mbox{ where}\qquad\mathcal{D}:=(D,\nabla D,\partial_{t}D).$$ (6)

Defining the DtP (dual-to-primal) map $(\mathcal{D},\bar{U})\mapsto\hat{U}(\mathcal{D},\bar{U})$ as the solution of

$$\frac{\partial\mathcal{L}}{\partial U}\big(U,\mathcal{D},\bar{U}\big)=0$$

for $U$ in terms of $(\mathcal{D},\bar{U})$ so that

$$\frac{\partial\mathcal{L}}{\partial U}\big(\hat{U}(\mathcal{D},\bar{U}),\mathcal{D},\bar{U}\big)=0,$$ (7)

we have

$$\hat{u}_{m}-\bar{u}_{m} =\left(-\partial_{t}\gamma_{m}-\partial_{l}\mu_{ml}\right)$$ (8a)

$$\hat{v}_{m}-\bar{v}_{m} =\frac{1}{\mathsf{r}}\left(-\gamma_{m}-\rho\partial_{t}\lambda_{m}\right)$$ (8b)

$$\hat{E}_{rs}-\bar{E}_{rs} =A^{-1}_{rsmn}(C_{ijmn}\partial_{j}\lambda_{i}-\mu_{(mn)})$$ (8c)

$$\hat{W}_{mn}-\bar{W}_{mn} =-\mu_{[mn]},$$ (8d)

where square brackets around a pair of indices represents antisymmetrization e.g., $\left(A_{[ij]}=\frac{1}{2}(A_{ij}-A_{ji})\right)$ and ordinary brackets, symmetrization $\left(A_{(ij)}=\frac{1}{2}(A_{ij}+A_{ji})\right)$. Also, for any second-order tensor, say $D$, $D^{(s)},D^{(a)}$ will represent its symmetric and antisymmetric parts, respectively.

Now define the dual functional, $S$, by substituting $\hat{U}$ for $U$ in the pre-dual functional $\widehat{S}$:

$$S[D] =\int_{0}^{T}\int_{\Omega}\mathcal{L}(\hat{U}(\mathcal{D},\bar{U}),\mathcal{D},\bar{U})\,dxdt$$ (9)

$$-\int_{0}^{T}\int_{\Omega}\lambda_{i}b_{i}\,dxdt-\int_{\Omega}\gamma_{i}(x,0)u_{i}^{(0)}(x)+\rho(x)\lambda_{i}(x,0)v_{i}^{(0)}(x)\,dx$$

$$-\int_{0}^{T}\int_{\partial\Omega_{s}}\lambda_{i}t_{i}\,dadt+\int_{0}^{T}\int_{\partial\Omega_{u}}\mu_{kl}u^{(b)}_{k}n_{l}\,dadt$$

with the essential space-time boundary conditions (5). Then, noting (7) and that the Lagrangian (6) is affine in $\mathcal{D}$ by construction, it is easily seen that the Euler-Lagrange equations of the functional $S$ is the governing set of equations and side conditions of linear Cauchy elastodynamics (2) with the replacement $U\to\hat{U}$, regardless of the choices of the specified functions

$$(x,t)\quad\mapsto\quad A(x,t),\quad\mathsf{r}(x),\quad\bar{u}(x,t),\quad\bar{v}(x,t),\quad\bar{E}(x,t),\quad\bar{W}(x,t),$$

subject to $\mathsf{r}\neq 0$ and $A$ merely invertible on the space of symmetric second-orders tensors.

Let us now assume that $\mathsf{r}>0$ and $A$ is positive-definite on the space of symmetric second order tensors and note that

$$\widehat{S}[U,D]=\int_{0}^{T}\int_{\Omega}\mathcal{L}(U,\mathcal{D},\bar{U})\,dxdt+\mbox{ boundary terms involving dual and specified fields}.$$

Then, due to the lack of any differential constraints on $U$ in the definition of $\mathcal{L}$,

$$\int_{0}^{T}\int_{\Omega}\mathcal{L}(\hat{U}(\mathcal{D},\bar{U}),\mathcal{D},\bar{U})\,dxdt=\int_{0}^{T}\int_{\Omega}\sup_{U}\mathcal{L}(U,\mathcal{D},\bar{U})\,dxdt=\sup_{U}\int_{0}^{T}\int_{\Omega}\mathcal{L}(U,\mathcal{D},\bar{U})\,dxdt.$$

Therefore,

$$S[D]=\sup_{U}\widehat{S}[U,D],$$

and because $\widehat{S}$ is affine in $D$ by construction, and therefore convex in $D$, defining our problem statement as finding the infimum of $S$, i.e.,

$$\inf_{D}S[D]\ =\ \inf_{D}\sup_{U}\,\widehat{S}[U,D];\qquad D^{*}=\arg\!\min_{D}\ S[D]$$

corresponds to a minimization problem for a convex functional, solving formally the problem of linear Cauchy Elastodynamics (statics is a special case) through the use of the DtP map, including for heterogeneous materials.

To obtain the explicit form of the dual functional, the (dual) Lagrangian may be expressed, using the DtP map (8), as

$$\mathcal{L}(\hat{U},\mathcal{D},\bar{U})+H(\hat{U},\bar{U}) =(\hat{u}_{i}-\bar{u}_{i})\hat{u}_{i}+\mathsf{r}(\hat{v}_{i}-\bar{v}_{i})\hat{v}_{i}$$

$$\quad+\hat{E}_{mn}A_{mnrs}(\hat{E}_{rs}-\bar{E}_{rs})+(\hat{W}_{mn}-\bar{W}_{mn})\hat{W}_{mn}$$

$$=2H(\hat{U},\bar{U})$$

$$\quad\ +(\hat{u}_{i}-\bar{u}_{i})\bar{u}_{i}+\mathsf{r}(\hat{v}_{i}-\bar{v}_{i})\bar{v}_{i}$$

$$\quad\ +A_{mnrs}(\hat{E}_{rs}-\bar{E}_{rs})\bar{E}_{mn}+(\hat{W}_{mn}-\bar{W}_{mn})\bar{W}_{mn}$$

$$=(\partial_{t}\gamma_{m}+\partial_{l}\mu_{ml})(\partial_{t}\gamma_{m}+\partial_{l}\mu_{ml})+\frac{1}{\mathsf{r}}\left(\gamma_{m}+\rho\partial_{t}\lambda_{m}\right)\left(\gamma_{m}+\rho\partial_{t}\lambda_{m}\right)$$

$$\quad+\left(C_{rsmn}\partial_{s}\lambda_{r}-\mu_{(mn)}\right)A^{-1}_{mnij}\,A_{ijkl}\,A^{-1}_{klpq}\left(C_{rspq}\partial_{s}\lambda_{r}-\mu_{(pq)}\right)+\mu_{[mn]}\mu_{[mn]}$$

$$\quad-(\partial_{t}\gamma_{m}+\partial_{l}\mu_{ml})\bar{u}_{m}-\left(\gamma_{m}+\rho\partial_{t}\lambda_{m}\right)\bar{v}_{m}$$

$$\quad+(C_{mnij}\partial_{j}\lambda_{i}-\mu_{(mn)})\bar{E}_{mn}-\mu_{[mn]}\bar{W}_{mn}.$$

Thus,

$$S[D] =\frac{1}{2}\int_{0}^{T}\int_{\Omega}\Big(|\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu|^{2}+\frac{1}{\mathsf{r}}|\gamma+\rho\partial_{t}\lambda|^{2}$$

$$\qquad\qquad\qquad+\left(C^{T}\nabla\lambda-\mu^{(s)}\right):A^{-1}\left(C^{T}\nabla\lambda-\mu^{(s)}\right)\ +\ \mu^{(a)}:\mu^{(a)}\Big)\,dxdt$$

$$\quad-\int_{0}^{T}\int_{\Omega}\Big((\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu)\cdot\bar{u}+(\gamma+\mathsf{r}\partial_{t}\lambda)\cdot\bar{v}+(\mu-C\nabla\lambda):\bar{E}+\mu^{(a)}:\bar{W}\Big)\,dxdt$$

$$\quad-\int_{0}^{T}\int_{\Omega}\lambda\cdot b\,dxdt-\int_{\Omega}\Big(\gamma(x,0)\cdot u^{(0)}+\rho\lambda(x,0)\cdot v^{(0)}\Big)\,dx$$ (11)

$$\quad-\int_{0}^{T}\int_{\partial\Omega_{s}}\lambda\cdot t\,dadt+\int_{0}^{T}\int_{\partial\Omega_{u}}(\mu n)\cdot u^{(b)}\,dadt.$$

with its fields satisfying the following essential boundary conditions on the space-time domain $\Omega\times(0,T)$:

$$\gamma=0\mbox{ on }\Omega\times\{T\}; \lambda=0\mbox{ on }\Omega\times\{T\};$$ (12)

$$\lambda=0\mbox{ on }\partial\Omega_{u}\times(0,T); \mu n=0\mbox{ on }\partial\Omega_{s}\times(0,T).$$

For $A$ positive-definite and $\mathsf{r}>0$, the semi-definiteness of the second variation of the representation (11) of the dual functional furnishes another (direct) proof of its convexity, regardless of the symmetry or positivity of $C$ or $\rho$.

The DtP map for the problem is given by

$$\hat{u} =-(\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu)+\bar{u}$$ (13)

$$\hat{v} =-\frac{1}{\mathsf{r}}\left(\gamma+\rho\,\partial_{t}\lambda\right)+\bar{v}$$

$$\hat{E} =-A^{-1}\left(\mu^{(s)}-C^{T}\nabla\lambda\right)+\bar{E}$$

$$\hat{W} =-\mu^{(a)}+\bar{W}.$$

The system of linear elastodynamics along with its initial and boundary conditions, expressed in terms of the dual fields, is

$$\begin{cases}&\partial_{t}\big(\bar{u}-(\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu)\big)-\big(\bar{v}-\frac{1}{\mathsf{r}}\left(\gamma+\rho\,\partial_{t}\lambda\right)\big)=0\\ &\rho\partial_{t}\big(\bar{v}-\frac{1}{\mathsf{r}}\left(\gamma+\rho\,\partial_{t}\lambda\right)\big)-\mathop{\rm div}\nolimits\Big(C\big(\bar{E}-A^{-1}\left(\mu^{(s)}-C^{T}\nabla\lambda\right)\big)\Big)-b=0\\ &\nabla\big(\bar{u}-(\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu)\big)-\big(\bar{E}-A^{-1}\left(\mu^{(s)}-C^{T}\nabla\lambda\right)\big)-\big(\bar{W}-\mu^{(a)}\big)=0\end{cases}\qquad\mbox{ on }\Omega\times(0,T)$$

$$\bar{u}-(\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu)=u^{(b)}\qquad\mbox{ on }\partial\Omega_{u}\times(0,T)$$

$$\Big(C\Big(\bar{E}-A^{-1}\big(\mu^{(s)}-C^{T}:\nabla\lambda\big)\Big)\Big)n=t\qquad\mbox{ on }\partial\Omega_{s}\times(0,T)$$

$$\bar{u}-(\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu)=u^{(0)}\qquad\mbox{ on }\Omega\times\{0\}$$ (14)

$$\bar{v}-\frac{1}{\mathsf{r}}\left(\gamma+\rho\,\partial_{t}\lambda\right)=v^{(0)}\qquad\mbox{ on }\Omega\times\{0\}$$

$$\begin{cases}\gamma=0\qquad\mbox{ on }\Omega\times\{T\}\\ \lambda=0\qquad\mbox{ on }\Omega\times\{T\}\\ \lambda=0\qquad\mbox{ on }\partial\Omega_{u}\times(0,T)\\ \mu n=0\qquad\mbox{ on }\partial\Omega_{s}\times(0,T)\end{cases}\quad\mbox{(these space-time b.c.s can be arbitrarily assigned).}$$

The following observations are in order:

1. When the base state, $\bar{U}$, is a solution to the system (2), $D=0$ is a solution to the dual system (14) and an (absolute) minimizer of the dual functional, $S$ (11). This is an important property of the dual scheme for both theoretical and practical purposes [17], as it motivates the reasoning that if a base state $\bar{U}$ is ‘close’ to a solution of (2), then it may be possible to obtain a solution of (2) by using $D=0$ as a good guess to the corresponding dual problem (14).

2. The first variation of the dual functional $S$ (9) imposes the linear elastodynamics equations (2) in the weak form - thus, all jump conditions of the primal problem are preserved by any extremal of the dual variational problem.

3. When solutions of the primal system (2) are unique, any solution to the dual system (14), generated from, e.g., different choices of boundary conditions (5) or $(A,\mathsf{r},\bar{U})$, must generate the unique primal solution through the use of the DtP map. Examples of this fact are provided in [8] in the context of the heat equation. However, if the primal system is known to have non-unique solutions, the auxiliary potential, $H$ parametrized by $(A,\mathsf{r},\bar{U})$, can be used to obtain such solutions through the dual scheme in a ‘stable’ manner - e.g., even when $C$ is indefinite. In other words, the use of the auxiliary potential, $H$, acts as a selection criterion for solutions of the primal problem.

4. Given that the primal problem (2) is an initial value-problem a natural question is whether prescribing final-time Dirichlet boundary conditions on the dual fields interferes in any way with causality of the primal solutions obtained from the above dual approach. Such interference does not arise because the DtP map (8) involves time derivatives of the dual fields, and prescribing the values of dual fields at final-time $T$ allows these time derivatives (at time $T$) to adjust so as to recover the correct primal solution. This feature has been explained and demonstrated in [2, Sec. 7], [8, 9, 16], with an exact solution for the first-order initial value problem $\dot{u}=au,u(0)=u_{0}$ worked out by the dual scheme in [16, Sec. 3.1].

5. Due to the availability of the functional $S$, a gradient descent algorithm in a fake time-like variable, can be formulated to obtain solutions to the primal system (2), see [17, Sec. 3]. The fact that $S$ is convex also guarantees that the norm of its gradient is non-increasing along a gradient flow - this is of great importance in a solution scheme that employs more than one $S$ functional (based on choices of base states) while keeping the gradient continuous at such switches, with the final goal of reaching a vanishing norm of the gradient (a weak solution to the primal equations is then obtained). Such ideas can be useful for iterative schemes for linear elasticity of heterogeneous materials.

6. The point of view adopted here is that all physics is contained in the primal system and the variational scheme is simply a mathematical device to solve that system. Thus, the dual fields are unlikely to have much physical significance - e.g., they admit acausal boundary conditions. Also, in nonlinear problems the use of a whole sequence of dual functionals, parametrized by base states, is employed to obtain a single primal solution; physically important quantities are not usually found in such abundance.

7. From a practical standpoint, having to solve a boundary-value-problem on a large space-time domain as might be the case to probe long-time behavior can be prohibitive. Solution of the dual system (14) can be divided into space-time domain slices, with an arbitrary, finite, sub-division of the time interval $[0,T]$. In each such space-time slice the dual problem can be solved, the primal solution recovered at the final time of the slice and used to define the initial conditions for the next slice. This strategy for solving the dual problem has been demonstrated in the solution of Euler’s ODE system for motion of a rigid body about a fixed point (with and without damping), the heat equation, and the linear transport equation in [8], and in the context of the (in)viscid Burgers equation in [9].

3 Formal uniqueness for the dual system (14) and its degenerate ellipticity

For practical purposes related to the use of the dual system (14), it is reassuring to have an assertion of uniqueness of solutions for any specific set of dual Dirichlet boundary conditions (5) and base states $\bar{U}$, at least when the elastic modulus $C$ is positive definite and has major symmetry. We provide such a sufficient condition in this Section. In the static case, uniqueness holds simply with positive definiteness without an assumption of major symmetry and we sketch that proof as well.

Consider two solutions of system (14) for $(\gamma,\lambda,\mu)$ and let their difference be denoted as $(\gamma^{*},\lambda^{*},\mu^{*})$. Then

$$\begin{cases}&\partial_{t}(\partial_{t}\gamma^{*}+\mathop{\rm div}\nolimits\mu^{*})-\frac{1}{\mathsf{r}}\left(\gamma^{*}+\rho\,\partial_{t}\lambda^{*}\right)=0\\ &\rho\partial_{t}\big(\frac{1}{\mathsf{r}}\left(\gamma^{*}+\rho\,\partial_{t}\lambda^{*}\right)\big)-\mathop{\rm div}\nolimits\Big(C\big(A^{-1}\left(\mu^{*(s)}-C^{T}\nabla\lambda^{*}\right)\big)\Big)=0\\ &\nabla(\partial_{t}\gamma^{*}+\mathop{\rm div}\nolimits\mu^{*})-A^{-1}\left(\mu^{*(s)}-C^{T}\nabla\lambda^{*}\right)-\mu^{*(a)}=0\\ \end{cases}\qquad\mbox{ on }\Omega\times(0,T).$$

Also, the difference fields satisfy the following space-time boundary conditions:

$$\partial_{t}\gamma^{*}+\mathop{\rm div}\nolimits\mu^{*} =0\mbox{ on }\partial\Omega_{u}\times(0,T)$$ (16a)

$$\left(C\left(A^{-1}\big(\mu^{*(s)}-C^{T}\nabla\lambda^{*}\big)\right)\right)n =0\mbox{ on }\partial\Omega_{s}\times(0,T)$$ (16b)

$$(\partial_{t}\gamma^{*}+\mathop{\rm div}\nolimits\mu^{*}) =0\mbox{ on }\Omega\times\{0\}$$ (16c)

$$\frac{1}{\mathsf{r}}\left(\gamma^{*}+\rho\,\partial_{t}\lambda^{*}\right) =0\mbox{ on }\Omega\times\{0\}$$ (16d)

$$\gamma^{*} =0\mbox{ on }\Omega\times\{T\}$$ (16e)

$$\lambda^{*} =0\mbox{ on }\Omega\times\{T\}$$ (16f)

$$\lambda^{*} =0\mbox{ on }\partial\Omega_{u}\times(0,T)$$ (16g)

$$\mu^{*}n =0\mbox{ on }\partial\Omega_{s}\times(0,T).$$ (16h)

On forming products of the three sets of equations in (15) with $(\gamma^{*},\lambda^{*},\mu^{*})$, respectively, and using the space-time boundary conditions (16a)-(16d), one obtains,

$$0 =\int_{0}^{T}\int_{\Omega}\Big(|\partial_{t}\gamma^{*}+\mathop{\rm div}\nolimits\mu^{*}|^{2}+\frac{1}{\mathsf{r}}|\gamma^{*}+\rho\partial_{t}\lambda^{*}|^{2}$$

$$+\left(C^{T}\nabla\lambda^{*}-\mu^{*(s)}\right):A^{-1}\left(C^{T}\nabla\lambda^{*}-\mu^{*(s)}\right)\ +\ \mu^{*(a)}:\mu^{*(a)}\Big)\,dxdt,$$

and due to the positive-definiteness of $A$ on the space of symmetric tensors and $\mathsf{r}>0$,

$$\partial_{t}\gamma^{*}+\mathop{\rm div}\nolimits\mu^{*} =0$$ (17a)

$$\gamma^{*}+\rho\,\partial_{t}\lambda^{*} =0$$ (17b)

$$C^{T}\nabla\lambda^{*}-\mu^{*(s)} =0$$ (17c)

$$\mu^{*(a)} =0$$ (17d)

hold $a.e.$ on $\Omega\times(0,T)$. Then (17) (all four equations together) gives

$$\rho\partial_{tt}\lambda^{*}=\mathop{\rm div}\nolimits(C^{T}\nabla\lambda^{*}),$$ (18)

which further gives, after forming a scalar product with $\partial_{t}\lambda^{*}$ and integrating over $\Omega$, when $C$ has major symmetry,

$$\partial_{t}\int_{\Omega}\left(\frac{1}{2}\rho\,|\partial_{t}\lambda^{*}|^{2}\ +\ \frac{1}{2}\nabla\lambda^{*}:C^{T}\nabla\lambda^{*}\right)\,dx-\int_{\partial\Omega}\partial_{t}\lambda^{*}\cdot((C^{T}\nabla\lambda^{*})n)\,da=0.$$ (19)

Now, from (16g) one has that

$$\partial_{t}\lambda^{*}=0\mbox{ on }\partial\Omega_{u}\times(0,T).$$ (20)

Assuming that (using (17d) as motivation)

$$\mu^{*(a)}n=0\mbox{ on }\partial\Omega_{s}\times(0,T)$$ (21)

holds, and combining with (16h) gives

$$\mu^{*(s)}n=0\mbox{ on }\partial\Omega_{s}\times(0,T).$$ (22)

Assuming as well that (using (17c) as motivation),

$$\left(C^{T}\nabla\lambda^{*}-\mu^{*(s)}\right)\,n=0\mbox{ on }\partial\Omega_{s}\times(0,T),$$ (23)

(20)-(22)-(23) imply that the boundary term in (19) vanishes so that

$$\int_{\Omega\times\{T\}}\left(\frac{1}{2}\rho\,|\partial_{t}\lambda^{*}|^{2}\ +\ \frac{1}{2}\nabla\lambda^{*}:C^{T}\nabla\lambda^{*}\right)\,dx-\int_{\Omega\times\{t\}}\left(\frac{1}{2}\rho\,|\partial_{t}\lambda^{*}|^{2}\ +\ \frac{1}{2}\nabla\lambda^{*}:C^{T}\nabla\lambda^{*}\right)\,dx=0$$ (24)

for all $t\in(0,T)$. Assuming (using (17b) as motivation)

$$\gamma^{*}+\rho\,\partial_{t}\lambda^{*}=0\mbox{ on }\partial\Omega\times\{T\}$$ (25)

and noting (16e) and (16f), the first integral in (24) vanishes and when $\rho>0$ and $C$ is positive-definite (major symmetry has been assumed), we have that

$$\partial_{t}\lambda^{*} =0\quad\mbox{ on }\Omega\times(0,T)$$ (26a)

$$\quad\nabla\lambda^{*(s)} =0\quad\mbox{ on }\Omega\times(0,T).$$ (26b)

This implies, from (17c)-(17d) (and an assumption of continuous extension to the boundary) that

$$\mu^{*}=0\mbox{ on }\bar{\Omega}\times[0,T],$$

which combined with (17a) and (16e) gives

$$\gamma^{*}=0\mbox{ on }\bar{\Omega}\times[0,T].$$

Of course, (26a) along with (16f) implies that

$$\lambda^{*}=0\mbox{ on }\bar{\Omega}\times[0,T].$$

The (formal) argument above for dual elastodynamics may be summarized as follows:

Proposition 3.1

Consider a class of functions $(\gamma,\lambda,\mu)$, $\gamma,\lambda:\bar{\Omega}\times[0,T]\to\mathbb{R}^{3}$, $\mu:\bar{\Omega}\times[0,T]\to\mathbb{R}^{3\times 3}$, such that the difference $(\Delta\gamma,\Delta\lambda,\Delta\mu)$ of any two functions in the class satisfy the conditions (16), (21), (23), (25) with the replacement $(\gamma^{*},\lambda^{*},\mu^{*})\to(\Delta\gamma,\Delta\lambda,\Delta\mu)$. Furthermore, assume that $\Delta\mu$ is continuous on $\bar{\Omega}\times[0,T]$.

When $\rho>0$ and $C$ is positive-definite with major symmetry, there can be at most one solution in this class of functions to the system (14).

Uniqueness for dual Elastostatics: In the problem of dual elastostatics, the $\gamma,\gamma^{*}$ dual fields are redundant, and the difference solution $(\lambda^{*},\mu^{*})$ satisfies

$$\begin{cases}&-\mathop{\rm div}\nolimits\Big(C\big(A^{-1}\left(\mu^{*(s)}-C^{T}\nabla\lambda^{*}\right)\big)\Big)=0\\ &\nabla\mathop{\rm div}\nolimits\mu^{*}-A^{-1}\left(\mu^{*(s)}-C^{T}\nabla\lambda^{*}\right)-\mu^{*(a)}=0\\ \end{cases}\qquad\mbox{ on }\Omega,$$

with the boundary conditions

$$\left(C\left(A^{-1}\big(\mu^{*(s)}-C^{T}\nabla\lambda^{*}\big)\right)\right)n =0\mbox{ on }\partial\Omega_{s}$$ (28a)

$$\lambda^{*} =0\mbox{ on }\partial\Omega_{u}$$ (28b)

$$\mu^{*}n =0\mbox{ on }\partial\Omega_{s}.$$ (28c)

Then, following the previous arguments for the dynamic case, one has

$$\mathop{\rm div}\nolimits\mu^{*} =0\mbox{ on }\Omega$$ (29a)

$$C^{T}\nabla\lambda^{*}-\mu^{*(s)} =0\mbox{ on }\Omega$$ (29b)

$$\mu^{*(a)} =0\mbox{ on }\Omega,$$ (29c)

which further implies

$$0=\mathop{\rm div}\nolimits(C^{T}\nabla\lambda^{*})\qquad\mbox{ on }\Omega.$$ (30)

Forming the scalar product with $\lambda^{*}$, integrating by parts and using (28b) gives

$$-\int_{\Omega}\nabla\lambda^{*}:C^{T}\nabla\lambda^{*}\,dx+\int_{\partial\Omega_{s}}\lambda^{*}\cdot(C^{T}\nabla\lambda^{*})n\,da=0$$ (31)

and assuming (29b)-(29c) also hold on $\partial\Omega_{s}$, and combining with (28c), the boundary term in (31) vanishes.

Then, assuming $C$ to be positive-definite (but not necessarily symmetric, a difference from what was used in the proof of uniqueness for the dynamic case), one has $\nabla\lambda^{*(s)}=0$ on $\Omega$ and combining with (28b), uniqueness for $\lambda$ holds. But then, from (29b)-(29c), uniqueness for $\mu$ holds as well.

The (formal) argument above for dual elastostatics may be summarized as follows:

Proposition 3.2

Consider a class of functions $(\lambda,\mu)$, $\lambda:\bar{\Omega}\to\mathbb{R}^{3}$, $\mu:\bar{\Omega}\to\mathbb{R}^{3\times 3}$, such that the difference $(\Delta\lambda,\Delta\mu)$ of any two functions in the class satisfy the conditions i) (29b) and (29c) on $\partial\Omega_{s}$ and ii) (28), with the replacement $(\lambda^{*},\mu^{*})\to(\Delta\lambda,\Delta\mu)$.

When $C$ is positive-definite (major symmetry is not required), there can be at most one solution in this class of functions to the system

$$\begin{cases}&-\mathop{\rm div}\nolimits\Big(C\big(\bar{E}-A^{-1}\left(\mu^{(s)}-C^{T}\nabla\lambda\right)\big)\Big)-b=0\\ &\nabla\big(\bar{u}-\mathop{\rm div}\nolimits\mu\big)-\big(\bar{E}-A^{-1}\left(\mu^{(s)}-C^{T}\nabla\lambda\right)\big)-\big(\bar{W}-\mu^{(a)}\big)=0\end{cases}\qquad\mbox{ on }\Omega$$

$$\bar{u}-\mathop{\rm div}\nolimits\mu=u^{(b)}\qquad\mbox{ on }\partial\Omega_{u}$$

$$\Big(C\Big(\bar{E}-A^{-1}\big(\mu^{(s)}-C^{T}:\nabla\lambda\big)\Big)\Big)n=t\qquad\mbox{ on }\partial\Omega_{s}$$

$$\begin{cases}\lambda=0\qquad\mbox{ on }\partial\Omega_{u}\\ \mu n=0\qquad\mbox{ on }\partial\Omega_{s}\end{cases}\quad\mbox{(these b.c.s can be arbitrarily assigned).}$$

Returning to the full dual elastodynamic problem, it is to be noted that while uniqueness in the absence of strict convexity might seem surprising, degenerate elliptic equations are known to have uniqueness in many circumstances, see, e.g., [13] (here we deal with a system of equations).

Purely from the standpoint of the operator containing the highest order derivatives of the dual system (14) given by

$$-\begin{pmatrix}\partial_{tt}&0&\partial_{t}\partial_{m}\\ &&\\ 0&\qquad\frac{\rho_{o}^{2}}{\mathsf{r}}\partial_{tt}\qquad&\partial_{j}(C:A^{-1}:C^{T})_{ijmn}\partial_{n}\\ &&\\ \partial_{j}\partial_{t}&0&\partial_{j}\partial_{m}\end{pmatrix}\begin{pmatrix}\gamma_{i}\\ \\ \lambda_{m}\\ \\ \mu_{im}\end{pmatrix}$$

with the corresponding quadratic form

$$\frac{1}{2}\int_{0}^{T}\int_{\Omega}\Big(|\partial_{t}\gamma+\mathop{\rm div}\nolimits\mu|^{2}\ +\ \frac{\rho^{2}}{\mathsf{r}}|\partial_{t}\lambda|^{2}\ +\ \nabla\lambda:(C:A^{-1}:C^{T})\nabla\lambda\,dxdt,$$ (33)

the system (14) is at worst degenerate-elliptic since the quadratic form (33) is positive semi-definite, regardless of whether $C$ is positive-definite and $\rho$ is positive (assuming $A$ is positive definite and $\mathsf{r}>0$, which are free to choose). This property is true of the full quadratic form in (11) as well.

Thus, it is worthy of note that the degenerate ellipticity of the dual problem holds even when the primal elastodynamic problem is hyperbolic or ill-posed hyperbolic with no continuous dependence on initial data when $C$ is indefinite. Initial demonstrations of success in solving such problems in the context of the linear transport equation, the (in)viscid Burgers equation, and elastodynamics of a bar with non-convex strain energy have been provided in [8, 9, 15].

We end this Section by exploring the following question: even strictly convex linear elastodynamics admits stress-waves induced by boundary/initial conditions (or body forces), which travel at the linear elastic wave speed of the material. This naturally must involve a propagating strain discontinuity as well. How can a ‘close-to’ elliptic PDE formulation recover such intrinsically ‘hyperbolic’ behavior?

Consider the dual ‘action’ (11) (for $C_{ijkl}=C_{klij}$ and $C$ positive-definite) and one asks about solutions that have vanishing bulk ‘dual energy/action.’ This necessarily requires the system (17) to hold (where the ^∗ fields now represents the 0-dual-action solution under consideration), which further implies that (18) holds. But this directly implies, because of its hyperbolic nature, that $\nabla\lambda$ admits rank-one discontinuities in the usual way across the characteristic surfaces of the physical linear elastodynamic equation and, through the DtP map (8), in the physical strain field $\hat{E}$ produced by the dual scheme, regardless of the chosen $A$, $\mathsf{r}$, $\bar{U}$. On the other hand, in the static situation there can be no discontinuities in $\lambda$ (for $C$ +ve-definite) by the usual properties of elliptic equations. Thus, we conclude that the degenerate ellipticity of the dual system (14) is crucial for recovering physically mandated behavior as per classical linear elasticity, and ellipticity of the dual system for elastodynamics must fail. These are all manifestations of a convex minimum principle which is not strictly convex (cf., [3, Sec. 3]).

The above arguments lead to a natural conjecture: that the dual formulation, due to its degenerate ellipticity selects only ‘good’ solutions, i.e. ones with fairly mild discontinuities, for elastodynamics with indefinite $C$ tensor.

4 Heterogeneous materials

For application to heterogeneous materials, i.e., $C,\rho$ vary in space, it is useful to write the system (14) in the form (we focus only on the governing equations)

$$-\partial_{tt}\gamma-\partial_{t}\mathop{\rm div}\nolimits\mu+\gamma \ =\ -\partial_{t}\bar{u}+\bar{v}-\frac{1}{\mathsf{r}}\gamma\+-\frac{\rho}{\mathsf{r}}\partial_{t}\lambda+\gamma$$

$$-\frac{\rho^{2}}{\mathsf{r}}\partial_{tt}\lambda-\mathop{\rm div}\nolimits\Big(\Big({\color[rgb]{0,0,0}CA^{-1}C^{T}}\Big)\nabla\lambda\Big) \ =\ -\rho\partial_{t}\bar{v}+\mathop{\rm div}\nolimits(C\bar{E})+\frac{\rho}{\mathsf{r}}\partial_{t}\gamma-\mathop{\rm div}\nolimits\Big(\Big({\color[rgb]{0,0,0}CA^{-1}}\Big)\mu^{(s)}\Big)+b$$

$$-\nabla\partial_{t}\gamma-\nabla\mathop{\rm div}\nolimits\mu+\mu \ =\ -\nabla\bar{u}+\bar{E}+\bar{W}+\mu^{(s)}-A^{-1}\mu^{(s)}+({\color[rgb]{0,0,0}A^{-1}C^{T}})\nabla\lambda.$$

In the above, in the first equation we have added a term $\gamma$ to both sides of the equation; for the third equation we have done the same using $\mu^{(s)}$.

With the goal of converting the left-hand-side operator of (34) to a constant coefficient one (the ‘comparison medium’), let $S^{(0)}$ be a positive-definite, spatially constant, compliance tensor with major and minor symmetries which is otherwise arbitrary and $\rho_{0}>0$ a spatially constant density field. Define

$$C^{(0)}=\left(S^{(0)}\right)^{-1}$$

and assume $C$ (and $C^{T}$) are invertible on the space of second order tensors. One then makes the choice

$$A:={\color[rgb]{0,0,0}C^{T}S^{(0)}C}\qquad\mbox{ and }\qquad\mathsf{r}:=\frac{\rho^{2}}{\rho_{0}},$$ (35)

resulting in (34) taking the form

$$-\partial_{tt}\gamma-\partial_{t}\mathop{\rm div}\nolimits\mu+\gamma=-\partial_{t}\bar{u}+\bar{v}-\frac{\rho_{0}}{\rho^{2}}\gamma-\frac{\rho_{0}}{\rho}\partial_{t}\lambda+\gamma$$ (36a)

$$-\rho_{0}\partial_{tt}\lambda-\mathop{\rm div}\nolimits\Big(C^{(0)}\nabla\lambda\Big)=-\rho\partial_{t}\bar{v}+\mathop{\rm div}\nolimits(C\bar{E})+\frac{\rho_{0}}{\rho}\partial_{t}\gamma-\mathop{\rm div}\nolimits\Big(\Big({\color[rgb]{0,0,0}C^{(0)}({C^{T}})^{-1}}\Big)\mu^{(s)}\Big)+b$$ (36b)

$$-\nabla\partial_{t}\gamma-\nabla\mathop{\rm div}\nolimits\mu+\mu=-\nabla\bar{u}+\bar{E}+\bar{W}+\mu^{(s)}-\Big({\color[rgb]{0,0,0}C^{-1}C^{(0)}(C^{T})^{-1}}\Big)\mu^{(s)}+\Big({\color[rgb]{0,0,0}C^{-1}C^{(0)}}\Big)\nabla\lambda.$$ (36c)

It is worthy of note that obtaining a constant-coefficient highest order operator on the l.h.s. of (36) does not come at the expense of introducing a highest order operator on the r.h.s. of the equation, as is effectively the case in the conventional method used to solve the linear elastic problem for heterogeneous materials [6, 19, 12] (we will refer to this as the ‘established formulation’). There, the addition arises from the (negative) divergence of the stress polarization tensor in the static case:

$$\mathop{\rm div}\nolimits(C\nabla u)=\mathop{\rm div}\nolimits(C^{(0)}\nabla u)+\mathop{\rm div}\nolimits\tau\ ;\qquad\qquad\tau=\left(C-C^{(0)}\right)\nabla u.$$

To elaborate a bit further, consider the elastostatic case and ignore the base states $(\bar{U})$ for the moment in (36). On defining the polarization tensor here to be

$${\color[rgb]{0,0,0}\tau=\left(C^{(0)}\Big(C^{T}\Big)^{-1}\right)\mu^{(s)}},$$

(36b) appears to have a very similar structure to the established formulation for solving the problem for a given stress polarization field, as described in the above references. Here, the polarization depends on the tensor field $\mu$ which, very roughly speaking, from (36c) would seem to have one less spatial derivative than $\nabla\lambda$, the latter being the analogous situation in the case in the established formulation. Thus, it is conceivable that solving for $\nabla\lambda$ and $\mu$ here, using a natural adaptation of the Moulinec-Suquet [12] iterative scheme for solving such integral equations, and through these fields the solution for the physical strain field through the DtP map (8c), has better properties. A similar argument would also be applicable to the problem of heterogeneous elastodynamics.

We note the following features, potentially important for possible practical application:

1. A gradient flow algorithm [17] based on the convex dual functional $S$ (11) can be utilized as an alternate methodology to the Moulinec-Suquet integral equation based iterative scheme [12] for obtaining (approximate) solutions to system (2) for a heterogeneous body, employing the choices (35).

2. In either approach, changes of base states in an iterative scheme [17, Sec. 3], based on progressive iterates for the primal fields, can be invoked to enhance convergence to a solution.

5 Concluding remarks

A scheme for generating variational minimum principles for Cauchy elastodynamics has been described as an application of a broader program to ease the solution/approximation of difficult (linear and nonlinear) PDE problems [3]. This is achieved through a conversion of the question to a convex variational one which allows the definition of a weaker notion of solutions (Variational dual solutions [15, 17, ASZ24]) to PDE than weak solutions, thus allowing for the use of tools from the Calculus of Variations and Convex Analysis in PDE. Such solutions have the consistency that when they can be proven to be sufficiently regular, they define genuine weak solutions of the primal PDE problem. The development lends itself to natural computational implementation for practical purposes.

Several positive aspects of the development in the context of Linear Elasticity have been pointed out. An unpalatable feature is that the transformation to an at most first-order differential system, due to the requirement that the pre-dual Lagrangian have no derivatives on primal variables appearing in it, increases the number of field degrees of freedom. The number of degrees of freedom in classical linear elastodynamics is 6 (we count velocity/momentum as a separate field as that is the optimal setting, especially when it comes to questions of heterogeneous materials, see. e.g., [18]). In contrast, the present formulation needs a minimum of 12 (6 for $\mu$ when (2)₃ is posed in symmetrized form, and three each for displacement, and velocity). Thus, while providing conceptual benefits and exposing unexpected connections, the present approach is sub-optimal for standard problems of linear elasticity from a practical point of view. However, for non-standard or difficult problems, e.g., a convex variational formulation of ‘Odd Elasticity’ [14], problems with indefinite elastic moduli (when the dynamic Cauchy problem changes type and is ill-posed due to a lack of continuous dependence w.r.t. initial data), problems with non-unique solutions, or for heterogeneous materials with high contrast, the methodology can be potentially beneficial. ‘Metamaterials’ are included in this class of problems.

While the availability of a convex variational principle for a general PDE system is a definite advantage, the considerations related to degenerate ellipticity at the end of Sec. 3 make it clear that weak coercivity of the convex dual functional, i.e., roughly speaking, the growth of the functional to $+\infty$ as a relevant norm of its argument function tends to $\infty$, is not to be expected. Thus, one of the three criteria for use of the ‘big hammer’ of the Calculus of Variations (cf., [20] for the Generalized Weierstrass Theorem) to prove existence of minimizers is not likely to be satisfied in many situations, when working within the present dual formulation. The formulation of Vorotnikov [V22], extending the pioneering work of Brenier [CMP18], is to be noted in this regard, but in any case, the gradient flow scheme proposed in [17] is perhaps the most pragmatic strategy to work with in efforts to obtain weak solutions of the original primal PDE through the dual scheme. That strategy is also immediately transferable to a computational algorithm for approximate solutions.

Finally, a speculation from a non-expert in the theory of homogenization: it seems like the natural transformation of the elastostatic(dynamic) problem of an arbitrary linear elastic composite to a space-time boundary value problem with constant-coefficient highest-order operator (36), along with a (formal) uniqueness guarantee in the usually encountered situations (major symmetry and +ve definiteness for all phases), may be expected to provide some advantages for mathematical homogenization schemes for such materials.