Energy minimising configurations of pre-strained multilayers

TL;DR

This paper studies the optimal configurations of pre-strained thin multilayer structures, revealing a phase transition from cylindrical to spherical shapes depending on the pre-strain strength, through theoretical analysis and numerical experiments.

Contribution

It introduces a family of von Kármán functionals interpolating between linearised regimes and rigorously analyzes the phase transition in optimal configurations.

Findings

Identification of a critical pre-strain level causing shape transition.

Rigorous convergence results for minimizers in asymptotic regimes.

Numerical evidence of a sharp transition at a specific parameter value.

Abstract

We investigate energetically optimal configurations of thin structures with a pre-strain. Depending on the strength of the pre-strain we consider a whole hierarchy of effective plate theories with a spontaneous curvature term, ranging from linearised Kirchhoff to von K\'arm\'an to linearised von K\'arm\'an theories. While explicit formulae are available in the linearised regimes, the von K\'arm\'an theory turns out to be critical and a phase transition from cylindrical (as in linearised Kirchhoff) to spherical (as in von linearised K\'arm\'an) configurations is observed there. We analyse this behavior with the help of a whole family $(\mathcal{I}^{\theta}_{\rm vK})_{\theta \in (0,\infty)}$ of effective von K\'arm\'an functionals which interpolates between the two linearised regimes. We rigorously show convergence to the respective explicit minimisers in the asymptotic regimes $\theta \to 0$ and $\theta \to \infty$. Numerical experiments are performed for general $\theta \in (0,\infty)$ which indicate a stark transition at a critical value of $\theta$.

Figures

Equations

$$\displaystyle\begin{split}u_{i}^{h}(x_{1},x_{2})&\coloneqq\frac{1}{(\sqrt{\theta}h)^{\gamma}}\int_{-1/2}^{1/2}\big{(}y^{h}_{i}(x_{1},x_{2},x_{3})-x_{i}\big{)}\,\mathrm{d}x_{3},\quad i=1,2,\\ v^{h}(x_{1},x_{2})&\coloneqq\frac{1}{(\sqrt{\theta}h)^{\alpha-2}}\int_{-1/2}^{1/2}y^{h}_{3}(x_{1},x_{2},x_{3})\,\mathrm{d}x_{3},\end{split}$$

$$\gamma=\left\{\begin{array}[]{rll}2(\alpha-2)&\text{if}&\alpha\in(2,3],\\ \alpha-1&\text{if}&\alpha\geq 3.\end{array}\right.$$

$$(x_{1},x_{2})\mapsto\big{(}(\sqrt{\theta}h)^{\gamma}u(x_{1},x_{2}),(\sqrt{\theta}h)^{\alpha-2}v(x_{1},x_{2})\big{)}.$$

$${\Omega }_h:=\omega \times (-h/2,h/2)\subset {\mathbb{R}}^3,$$

$$E_{\alpha }^h(y)={\int }_{{\Omega }_1}{W}_{\alpha }^h(x_3,{\partial }_{1}y,{\partial }_{2}y,h^{-1}{\partial }_{3}y),$$

$$W_{\alpha }^h(x_3,F)=W_0(x_3,F(I+h^{\alpha -1}{B}^h(x_3))),F\in {\mathbb{R}}^{3\times 3}.$$

$$W_{\alpha=3}^{h}(x_{3},F)=W_{0}\big{(}x_{3},F\big{(}I+h^{2}\sqrt{\theta}B^{h}(x_{3})\big{)}\big{)},\quad F\in\mathbb{R}^{3\times 3}.$$

$$Q_3(t,F):=D^2{W}_0(t,I)[F,F]=\frac{{\partial }^2{W}_0(t,I)}{\partial F_{ij}\partial F_{ij}}{F}_{ij}{F}_{ij},$$

$$Q_2(t,G):=\underset{c\in {\mathbb{R}}^3}{\min }{Q}_3(t,\hat{G}+c\otimes e_3),$$

$$Q_2(t,G)\le C|G{|}^2\forall G\in {\mathbb{R}}^{2\times 2}\text{ and }{Q}_2(t,G)\ge c|G{|}^2\forall G\in {\mathbb{R}}_{\mathrm{sym}}^{2\times 2}$$

$$\check{B}\in L^{\infty}\big{(}(-1/2,1/2),\mathbb{R}_{\operatorname{sym}}^{2\times 2}\big{)}.$$

$${\overline{Q}}_2[E,F]:={\int }_{-1/2}^{1/2}{Q}_2(t,E+tF+\check{B}(t))\mathrm{d}t,$$

$${\overline{Q}}_2^{\star }(F):=\underset{E\in {\mathbb{R}}_{\mathrm{sym}}^{2\times 2}}{\min }{\int }_{-1/2}^{1/2}{Q}_2(t,E+tF+\check{B}(t))\mathrm{d}t.$$

$$\mathcal{I}_{\rm lKi}(v)\coloneqq\left\{\begin{array}[]{rl}\frac{1}{2}\int_{\omega}\overline{Q}_{2}^{\star}(-\nabla^{2}v)&\text{ if }v\in W^{2,2}_{\rm sh}(\omega),\\ \infty&\text{ otherwise}.\end{array}\right.$$

$$\mathcal{I}^{\theta}_{\rm vK}(u,v)\coloneqq\left\{\begin{array}[]{l}\frac{1}{2}\int_{\omega}\overline{Q}_{2}[\theta^{1/2}(\nabla_{s}u+\tfrac{1}{2}\nabla v\otimes\nabla v),-\nabla^{2}v]\\ \text{{\hskip 50.00008pt}if }(u,v)\in W^{1,2}(\omega;\mathbb{R}^{2})\times W^{2,2}(\omega;\mathbb{R}),\\ \infty,\text{ otherwise}.\end{array}\right.$$

$$\mathcal{I}_{\rm lvK}(u,v)\coloneqq\left\{\begin{array}[]{l}\frac{1}{2}\int_{\omega}\overline{Q}_{2}[\nabla_{s}u,-\nabla^{2}v],\\ \text{{\hskip 50.00008pt}if }(u,v)\in W^{1,2}(\omega;\mathbb{R}^{2})\times W^{2,2}(\omega;\mathbb{R})\\ \infty,\text{ otherwise}.\end{array}\right.$$

$${\mathcal{I}}_{\mathrm{lKi}}(v)$$

$${\mathcal{I}}_{\mathrm{vK}}^{\theta }(u,v)$$

$${\mathcal{I}}_{\mathrm{lvK}}(u,v)$$

$$E\mapsto e:=(E_{11},E_{22},E_{12}),$$

$$Q_2(t,A)=a^{\top }M(t)a.$$

$$M_0:={\int }_{-1/2}^{1/2}M(t)\mathrm{d}t,M_1:={\int }_{-1/2}^{1/2}tM(t)\mathrm{d}t,M_2:={\int }_{-1/2}^{1/2}{t}^{2}M(t)\mathrm{d}t.$$

$$M^{\ast }:=M_2-M_1{M}_0^{-1}{M}_1$$

$$\displaystyle\int_{-1/2}^{1/2}\big{|}\big{(}tM^{1/2}(t)-M^{1/2}(t)\Lambda\big{)}x\big{|}^{2}\,\mathrm{d}t>0$$

$$0$$

$$\displaystyle=x^{\top}\big{(}M_{2}-\Lambda^{\top}M_{1}-M_{1}\Lambda+\Lambda^{\top}M_{0}\Lambda\big{)}\,x$$

$$0\char 60\relax x^{\top}\big{(}M_{2}-M_{1}M_{0}^{-1}M_{1}\big{)}x.$$

$${\overline{Q}}_2[E,F]\text{leavevmode}=\text{leavevmode}{\int }_{-1/2}^{1/2}{Q}_2(t,E+tF+\check{B}(t))\mathrm{d}t$$

$$\text{leavevmode}\text{leavevmode}=\text{leavevmode}{e}^{\top }{M}_{0}e+f^{\top }{M}_{2}f+{\beta }_0+2e^{\top }{M}_{1}f+2e^{\top }{b}_1+2f^{\top }{b}_2$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ =\leavevmode\nobreak\ \big{(}e+M_{0}^{-1}(M_{1}f+b_{1})\big{)}^{\top}M_{0}\big{(}M_{0}e+M_{0}^{-1}(M_{1}f+b_{1})\big{)}$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \qquad+\big{(}f+(M^{\ast})^{-1}(b_{2}-M_{1}M_{0}^{-1}b_{1})\big{)}^{\top}M^{\ast}\big{(}f+(M^{\ast})^{-1}(b_{2}-M_{1}M_{0}^{-1}b_{1})\big{)}$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \qquad-\big{(}M_{1}M_{0}^{-1}b_{1}\big{)}^{\top}(M^{\ast})^{-1}\big{(}M_{1}M_{0}^{-1}b_{1}\big{)}-b_{1}^{\top}M_{0}^{-1}b_{1}+\beta_{0}$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ =\leavevmode\nobreak\ \gamma+\big{(}e+M_{0}^{-1}(M_{1}f+b_{1})\big{)}^{\top}M_{0}\big{(}M_{0}e+M_{0}^{-1}(M_{1}f+b_{1})\big{)}$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \qquad+\big{(}f+(M^{\ast})^{-1}(b_{2}-M_{1}M_{0}^{-1}b_{1})\big{)}^{\top}M^{\ast}\big{(}f+(M^{\ast})^{-1}(b_{2}-M_{1}M_{0}^{-1}b_{1})\big{)},$$

$$\displaystyle\gamma:=-\big{(}M_{1}M_{0}^{-1}b_{1}\big{)}^{\top}(M^{\ast})^{-1}\big{(}M_{1}M_{0}^{-1}b_{1}\big{)}-b_{1}^{\top}M_{0}^{-1}b_{1}+\beta_{0}.$$

Full text

1. 1 Introduction

2. 2 Effective plate theories

3. 2.1 Dimension reduction for pre-strained multilayers

4. 2.2 Effective moduli and minimising strains

5. 3 Optimal configurations in the linearised and the asymptotic critical regimes

6. 4 Structure of minimisers for $\mathcal{I}^{\theta}_{\rm vK}$ for small $\theta$

7. 4.1 A branch of solutions for $\theta\ll 1$

8. 4.2 Uniqueness and globality of minimisers

9. 5 Discretisation of the interpolating theory

10. 5.1 Discretisation

11. 5.2 $\Gamma$-convergence of the discrete energies

12. 5.3 Discrete gradient flow

13. 5.4 Experimental results

1 Introduction

The topic of this paper is motivated by experimental observations on optimal energy configurations in thin (heterogeneous) structures with a pre-strain. The simplest example of such a structure is the classical bimetallic strip which consists of two strips of different materials with different thermal expansion coefficients joined together throughout their length. If heated or cooled, due to the misfit of equilibria, internal stresses develop. The flat reference configuration is no longer optimal and the strip bends in order to reduce elastic energy. This behavior can effectively be modelled with a 1d energy functional comprising a temperature dependent spontaneous curvature term.

In this paper we will investigate thin layers whose two lateral dimensions are much larger than their very small height and whose flat reference configuration is subject to internal stresses (one speaks of pre-strained or pre-stressed bodies). Examples of such structures are heated materials (with inhomogeneous expansion coefficients as in the bimetallic strip referred to above or homogeneous materials with a temperature gradient), crystallisations on top of a substrate as in epitaxially grown layers, or biological materials whose internal misfit is caused by swelling and growing tissue. Our main focus will be on multilayered heterogeneous plates, for which the effective plate theories have been provided in [11]. Our findings, however, apply equally to different situations as long as they are described by the same effective functionals, cf. Remark 1 below.

As a matter of fact, the situation is much more complicated and interesting for two dimensional plates than for one dimensional strips. It has been found that the assumed shape depends on the strength of the pre-strain and the aspect ratio of the specimen: Large pre-strains in very thin layers tend to cause cylindrical shapes whereas smaller pre-strains in thicker layers lead to spherical caps, [25, 31, 13, 14, 21, 12]. To explain this observation one argues that locally the energy is best released if a spherical shape is assumed. If, however, the aspect ratio is very small, i.e., the lateral dimensions are very large compared to the thickness, then this leads to geometric incompatibilities: non-zero Gauß curvature introduces a change of the metric which by far has too high elastic energy. In contrast, cylindrical shapes do not lead to such incompatibilities.

A thorough theoretical understanding of this mechanism through which ‘misfit’ of equilibria is converted into mechanical displacement is not only interesting from a mathematical point of view. In view of applications it has proved to constitute a convenient and feasible method to access and manipulate objects even at the nanoscale. By way of example we mention experiments on the self-organised fabrication of nano-scrolls, as reported in [34, 18, 27].

The aim of this paper is to shed light on the geometry of energetically optimal configurations of pre-strained heterostructures with the help of two-dimensional plate theories. More precisely, we consider effective plate theories for multilayers with reference configuration $\Omega_{h}=\omega\times(-h/2,h/2)$, $0\char 60\relax h\ll 1$, whose (small) misfit pre-strain is described by a matrix $h^{\alpha-1}B^{h}$, scaling with $h$.

The particular case $\alpha=2$ with a misfit of the order $h$ of the aspect ratio has been investigated in [32, 33, 7]. The appropriate plate theory is the nonlinear Kirchhoff theory (in the finite bending regime) and energy minimizers turned out to be (portions of) cylinders whose possible winding directions and radii are determined explicitly. Therefore, in order to be able to encounter different behavior one has to consider weaker scalings of the misfit.

In [11] – based on the homogeneous case explored in [15] – we have found a whole hierarchy of effective plate theories for the scalings $\alpha>2$. Suitably rescaled, one obtains only three different limiting plate theories: the linearised Kirchhoff theory for $\alpha\in(2,3)$, the von Kármán theory for $\alpha=3$ and the linearised von Kármán theory for $\alpha>3$. With a view to our present investigation, we have moreover derived a fine scale $\theta$ in the critical von Kármán scale which interpolates continuously between the the two linearised theories.

For such small misfits one is lead to describe a deformation $y^{h}:\Omega_{h}\to\mathbb{R}^{3}$ in terms of the scaled and averaged in-plane, respectively, out-of-plane displacements

[TABLE]

where $\theta\equiv 1$ unless $\alpha=3$ and

[TABLE]

A limiting plate theory in terms of the limiting quantities $(u,v)$ is then derived as the $\Gamma$-limit of the 3d nonlinearly elastic energy, rescaled by $h^{1-2\alpha}$, cf. [11]. For a minimizer $(u,v)$ of the limiting theory one obtains the shape of an optimal configuration at finite $0\char 60\relax h\ll 1$: After descaling, its $x_{3}$-averaged displacement is given approximately by

[TABLE]

Since $\gamma>\alpha-2$, the in-plane components are indeed much smaller than the out-of-plane component. In his sense, the shape is to leading order described by $v:\omega\to\mathbb{R}$ only.

In the linearised regimes our results give the following picture: If $\alpha\char 60\relax 3$, degenerate parabolas (infinitesimal parts of cylinders) are seen to be optimal, whereas for $\alpha>3$, non-degenerate parabolas (infinitesimal parts of an elliptical cap) are energy minimizers. Only in the latter case, however, the minimizer is unique (up to affine terms). Yet, even in case $\alpha\char 60\relax 3$ it turns out the geometric shape is uniquely determined as an infinitesimal part of a cylinder while the winding direction and radius may have several optimal values. In both cases we explicitly determine these minimizers. A basic observation shows that for $\alpha=3$ these configurations are still asymptotically optimal in the ‘almost linearised’ regimes $\theta\gg 1$ and $\theta\ll 1$, respectively.

The von Kármán regime is much more subtle. We focus on a prototypical functional in order to understand better the material response if the misfit (and hence $\theta$) is increased from [math] to a finite value. We show that for finite, although small, values of $\theta$ there is a unique branch of global minimizers emanating from a spherical cap. For a further study for general values of $\theta\in(0,\infty)$ we then rely on computer experiments. To this end, we develop a penalised, nonconforming finite element discretisation using $P^{1}$ elements and employ projected gradient descent to solve the ensuing nonlinear problems while ensuring constraints are met. We first show $\Gamma$-convergence of the discrete problems to the continuous one, then investigate the minimizers in their dependence on $\theta$. Interestingly, our results seem to indicate a stark change of material response at a critical value of $\theta$, showing a symmetry breaking ‘phase transition’ from a nearly spherical cap to an approximate cylinder.

Outline

We begin by recalling our main results from [11] in order to provide the appropriate plate theories in Section 2. There we also identify the effective elastic moduli and spontaneous curvature terms explicitly so as to transform the problem into a more amenable form to identify minimizers. We then discuss the linearised regimes $\alpha\in(2,3)$ and $\alpha>3$ as well as the asymptotic von Kármán regimes $\theta\to 0$ and $\theta\to\infty$ in Section 3. The structure of minimisers for small $\theta$ is investigated in Section 4. Finally, Section 5 contains our numerical findings.

2 Effective plate theories

We first recall the main results of our contribution [11] on a hierarchy of plate theories for pre-strained multilayers derived from non-linear three dimensional elasticity by $\Gamma$-convergence. We then determine the effective (homogenised) elastic moduli and corresponding quadratic energy desnities of the plates in terms of the moments of the pointwise elastic constants of the layers.

2.1 Dimension reduction for pre-strained multilayers

Working exactly in the setting of [11] we consider a thin domain

[TABLE]

where $\omega\subset\mathbb{R}^{2}$ is bounded with Lipschitz boundary, $0\char 60\relax h\ll 1$, subject to a deformation $w:\Omega_{h}\to\mathbb{R}^{3}$. Changing variables form $x_{3}$ to $x_{3}/h$ we obtain a deformation mapping $y(x)=w(x_{1},x_{2},hx_{3})$ and the energy per unit volume

[TABLE]

where the elastic energy density $W_{\alpha}^{h}$ depends on a scaling parameter $\alpha\in(2,\infty)$ and is given by

[TABLE]

for $\alpha\neq 3$, $B^{h}:\left(-1/2,1/2\right)\rightarrow\mathbb{R}^{3\times 3}$ describing the internal misfit and $W_{0}$ the stored energy density of the reference configuration. For $\alpha=3$ we include an additional parameter $\theta>0$ controlling further the amount of misfit in the model:

[TABLE]

We take $W_{0}$ fulfilling the usual assumptions of smoothness around $SO(3)$, frame invariance, boundedness and quadratic growth which are detailed in [11]. After linearising around the identity, one obtains the Hessian

[TABLE]

for $t\in\left(-1/2,1/2\right),F\in\mathbb{R}^{3\times 3}$ and defines $Q_{2}$ by minimising away the effect of transversal strain on $Q_{3}$:

[TABLE]

for $t\in\left(-1/2,1/2\right),G\in\mathbb{R}^{2\times 2}$, $e_{3}=(0,0,1)\in\mathbb{R}^{3}$, and $\hat{G}\in\mathbb{R}^{3\times 3}$ has $G$ as its upper left $2\times 2$ submatrix and zeros in the third column and the third row. The functions $Q_{2}(t,\cdot)$, $t\in(-1/2,1/2)$, are quadratic forms on $\mathbb{R}^{2\times 2}$ which are positive definite on $\mathbb{R}_{\operatorname{sym}}^{2\times 2}$ and vanish on antisymmetric matrices. Moreover, they satisfy the bounds

[TABLE]

for constants $c,C>0$ and a.e. $t\in(-1/2,1/2)$. We also denote by $\check{B}(t)$ the $2\times 2$ matrix which arises from $B(t)\in\mathbb{R}^{3\times 3}$ by deleting its last row and last column. Then

[TABLE]

From $Q_{2}(t,\cdot)$ and $\check{B}(t)$ we define the effective form:

[TABLE]

and its relaxation

[TABLE]

In [11] it is shown that $h^{2-2\alpha}E^{h}_{\alpha}$ $\Gamma$-converges for the convergence of the averaged in-plane and out-of-plane displacements $(u^{h},v^{h})\rightharpoonup(u,v)$ in $W^{1,2}(\omega;\mathbb{R}^{3})$ modulo a global rigid motion, cf. (1), to the following effective limiting functionals:

For the scaling $\alpha\in(2,3)$ as defined in [11] and convex $\omega$, the linearised Kirchhoff energy is given by

[TABLE]

For $\alpha=3$ we have the von Kármán type energy333As in [11] we slightly overload the notation in what would be a double definition of $\mathcal{I}_{3}^{h}$, using the letter in the subindex to dispel the ambiguity.

[TABLE]

Finally, in the regime $\alpha>3$ we have the linearised von Kármán energy

[TABLE]

Remark 1

The precise assumptions on $W^{h}_{\alpha}$ from [11] are not essential for the results of the present contribution. In what follows we will only need that the $Q_{2}(t,\cdot)$, $t\in(-1/2,1/2)$, are quadratic forms on $\mathbb{R}^{2\times 2}$ that vanish on antisymmetric matrices and satisfy (2) and that $\check{B}$ satisfies (3).

The existence of minimizers of (5) (6) and (7) follows by a standard application of the direct method or, in the setting of [11], as a direct consequence of $\Gamma$-convergence and compactness.

Example. For a homogeneous material $Q_{2}(t,A)=Q_{2}(A)$ with linear internal misfit $B(t)=tI$ one has

[TABLE]

for $v\in W^{2,2}_{\rm sh}(\omega)$, respectively, $(u,v)\in W^{1,2}(\omega;\mathbb{R}^{2})\times W^{2,2}(\omega;\mathbb{R})$. These functionals, where the elastic coefficients do not depend on the out-of-plane component, can model for instance a single-layer material under thermal stress. In Section 5, we will study the energy (8) as a function of $\theta$.

2.2 Effective moduli and minimising strains

This subsection serves to give explicit formulae relating the homogenised effective elastic moduli found above to the zeroth, first and second moment in $t$ of the individual $Q_{2}(t,\cdot)$. We also identify their pointwise minimiser so as to rewrite the effective quadratic forms in their most convenient form. The computations are completely elementary, we indicate the main steps.

Because $Q_{2}$ vanishes on antisymmetric matrices we may restrict our attention to $F\in\mathbb{R}^{2\times 2}_{\operatorname{sym}}$. From now on, we identify matrices $E=(E_{ij})_{i,j=1}^{2}\in\mathbb{R}^{2\times 2}_{\operatorname{sym}}$ with vectors in $\mathbb{R}^{3}$ via

[TABLE]

and analogously $F\mapsto f$, $\check{B}\mapsto b$, $A\mapsto a$. Then, for each $t\in\left(-1/2,1/2\right)$ there exists some symmetric, positive definite matrix $M(t)$ such that for all $A\in\mathbb{R}^{2\times 2}_{\operatorname{sym}}$:

[TABLE]

We define the moments of $M$ as

[TABLE]

It is easy to see that (2) implies that $M_{0}$ and $M_{2}$ are positive definite. We claim that also

[TABLE]

is positve definite. To see this, fix $\Lambda\in\mathbb{R}^{2\times 2}$ and note that for all $x\in\mathbb{R}^{2}\setminus\{0\}$

[TABLE]

since \big{(}tM^{1/2}(t)-M^{1/2}(t)\Lambda\big{)}x=0 for a.e. $t$ would imply that $(tI-\Lambda)x=0$ in contradiction to $\Lambda$ having at most two eigenvalues. Expanding the square we get

[TABLE]

and, choosing $\Lambda=M_{0}^{-1}M_{1}$,

[TABLE]

Let $\overline{Q}_{2}$ be given as in (2.1). Elementary calculations show that

[TABLE]

where

[TABLE]

We define the linear mappings $\mathcal{L}_{i},\mathcal{L}_{\ast}:\mathbb{R}^{2\times 2}_{\operatorname{sym}}\to\mathbb{R}^{2\times 2}_{\operatorname{sym}}$, $i=1,2,3$, by

[TABLE]

and the positive definite quadratic forms $Q_{2}^{0}$ and $Q_{2}^{\ast}$ on $\mathbb{R}^{2\times 2}_{\operatorname{sym}}$ by

[TABLE]

In terms of these quantities our computation reads

[TABLE]

with

[TABLE]

Minimizing out $E$ yields

[TABLE]

3 Optimal configurations in the linearised and the asymptotic critical regimes

In this section we develop a characterisation of minimisers for the lower range $\alpha\in(2,3)$ and for the upper range $\alpha>3$ of scalings. Recall from the discussion in Section 1 that we are primarily intested in the shape of the out-of-plane component $v$. The results indicate that the characteristic shapes in the limit $h\to 0$ are (infinitesimal) cylinders and paraboloids respectively. Invoking the $\Gamma$-convergence results with respect to the interpolation parameter $\theta$ from [11, Section 6] this will also shed light on the optimal shapes in the asymptotic regimes $\theta\to 0$ and $\theta\to\infty$ for the von Kármán scaling $\alpha=3$. We collect our results in the following three theorems, where indeed Theorem 1 is indeed rather an elementary observation based on our preparations form the previous section and Theorem 3 is a direct consequence of [11, Section 6]. We allow for a general bounded Lipschitz domain $\omega$ in these theorems.

Theorem 1

The minimisers of $\mathcal{I}_{\rm lvK}$, eq. (7), are of the form

[TABLE]

with $E_{0},F_{0}\in\mathbb{R}_{\operatorname{sym}}^{2\times 2}$ the constants from (13). $u$ is unique up to an infinitesimal rigid motion and $v$ up to the addition of an affine transformation.

Theorem 2

Up to the addition of an affine transformation, the minimisers of $\mathcal{I}_{\rm lKi}$, eq. (5), are of the form

[TABLE]

where $Q_{2}^{\ast},F_{0}$ are given in (11) and (13), respectively.

Remark 2

Describing symmetric $2\times 2$ matrices $A$ by vectors $a\in\mathbb{R}^{3}$ as in Section 2.2, the set $\mathcal{N}$ is the set of touching points of the two quadrics $\{a\in\mathbb{R}^{3}:a_{1}a_{2}-a_{3}^{2}=0\}$ (a cone) and $\{a\in\mathbb{R}^{3}:a^{\top}M^{\ast}a=c_{m}\}$ (an ellipsoid), where $c_{m}=Q_{2}^{\ast}(F-F_{0})$ with $F\in\mathcal{N}$. If $\#\mathcal{N}\geq 3$, intersecting with an affine plane $P$ containing three distinct points of $\mathcal{N}$ shows that $\mathcal{N}\cap P$ is an ellipse and then even $\mathcal{N}\subset P$. This shows that either $\#\mathcal{N}=1$ and there is a unique minimizer, or $\#\mathcal{N}=2$ and there are precisely two minimizers, or $\mathcal{N}$ is an affine ellipse and to each ‘winding direction’ $\mathbb{R}e$, $e\in S^{1}$, there is a unique curvature $\lambda=\lambda(e)$ such that $\nabla^{2}v\equiv\lambda e\otimes e$.

Theorem 3

Suppose that $(u^{\theta},v^{\theta})$ are minimisers of $\mathcal{I}^{\theta}_{\rm vK}$, eq. (6).

1. a)

As $\theta\to 0$, up to infinitesimal rigid motions in the in-plane component and up to the addition of affine transformations in the out-of-plane compenent, $(\theta^{1/2}u^{\theta},v^{\theta})\rightharpoonup(u,v)$ in $W^{1,2}(\omega,\mathbb{R}^{2})\times W^{2,2}(\omega;\mathbb{R})$ with $(u,v)$ as in (15). 2. b)

As $\theta\to\infty$, up to the addition of affine transformations in the out-of-plane component and up to passing to a subsequence, $v^{\theta}\rightharpoonup v$ in $W^{2,2}(\omega;\mathbb{R})$ with $v$ as in (16).

**Proof of Theorem 1 ** By (7) and (12)

[TABLE]

with $u\in W^{1,2}(\omega;\mathbb{R}^{2})$ and $v\in W^{2,2}(\omega;\mathbb{R})$ is minimal (with value $\gamma|\omega|/2$) if and only if $\nabla^{2}v=-F_{0}$ and $\nabla_{s}u=\mathcal{L}_{0}^{-1}\mathcal{L}_{1}F_{0}+E_{0}$ a.e. $\Box$

**Proof of Theorem 3 ** a) is immediate from [11, Theorems 7,10,11]. b) directly follows from [11, Theorems 7,8,9] if $\omega$ is convex. For general $\omega$ first note that the compactness result in [11, Theorem 7] does not use convexity, so that $v^{\theta}\rightharpoonup v$ in $W^{2,2}(\omega;\mathbb{R})$ for some $v\in W^{2,2}_{\rm sh}$. Now fix $F=(f_{ij})_{1\leqslant i,j\leqslant 2}\in\mathcal{N}$ and $\bar{v}(x)=\frac{1}{2}x^{\top}Fx$. Since $\det F=0$, the function $u^{\prime}(x)=-\frac{1}{3}f_{11}x_{1}^{3}(f_{11},f_{12})-f_{12}x_{1}^{2}x_{2}(f_{11},f_{12})-f_{12}x_{1}x_{2}^{2}(f_{12},f_{22})-\frac{1}{3}f_{22}x_{2}^{3}(f_{12},f_{22})$ satisfies $\nabla_{s}u^{\prime}+\frac{1}{2}\nabla\bar{v}\otimes\nabla\bar{v}=0$. Also choose $u^{\prime\prime}(x)=Ex$ with $E=\mathcal{L}_{0}^{-1}\mathcal{L}_{1}F+E_{0}$, cf. (12) and (13). Then for $\bar{u}=u^{\prime}+\theta^{-1/2}u^{\prime\prime}$ we have by (14)

[TABLE]

With the help of the Vitali covering theorem we can exhaust $\omega$ up to a set of negligible measure with disjoint convex subdomains $\omega_{1},\omega_{2},\ldots$. Denoting the accordingly restricted functionals by $\mathcal{I}^{\theta}_{\rm vK}(\ \cdot\ ;\omega_{n})$, $\mathcal{I}_{\rm lKi}(\ \cdot\ ;\omega_{n})$ we have

[TABLE]

where we have made use of the lower bound in the $\Gamma$-convegence of $\mathcal{I}^{\theta}_{\rm vK}(\cdot;\omega_{n})$ to $\mathcal{I}_{\rm lKi}(\cdot;\omega_{n})$, see [11, Theorem 8], in the third step and of Theorem 2 in the fourth step. So we must have $\mathcal{I}_{\rm lKi}(v;\omega_{n})=\mathcal{I}_{\rm lKi}(\bar{v};\omega_{n})$ for all $n$ and hence $\nabla^{2}v\in\mathcal{N}$ a.e. on $\omega$ and so the claim follows from Theorem 2. $\Box$

As for Theorem 2, it is straightforward to see that $v$ as defined in the theorem is a minimisers of $\mathcal{I}_{\rm lKi}$. However, the proof that every minimiser of $\mathcal{I}_{\rm lKi}$ is necessarily of this form needs some work. The difficulty lies in excluding the possibility of constructing a minimiser by piecing together functions whose Hessian belongs to the set $\mathcal{N}$, all with minimal energy but lacking a nice global structure. Yet it is possible to obtain a global representation of the Hessian which shows that it must be constant over $\omega$ so that minimisers are (up to an affine transformation) indeed cylindrical. In order to do this we require (cf. [28]):

Definition 1

Let $\omega^{\prime}\subset\mathbb{R}^{2}$ a convex bounded domain and $y\in W^{1,2}(\omega^{\prime},\mathbb{R}^{3})$ be an isometry. A connected maximal subdomain of $\omega^{\prime}$ where $\nabla y$ is constant and $y$ is affine whose boundary contains more than two segments inside $\omega^{\prime}$ is called a body. A leading curve is a curve orthogonal to the preimages of $\nabla y$ on the open regions where $\nabla y$ is not constant, parametrised by arc-length. We define an arm to be a maximal subdomain $\omega(\gamma)$ which is covered (parametrised) by some leading curve $\gamma$ as follows:

[TABLE]

where $\nu(t)=\gamma^{\prime}(t)^{\perp}$. We also speak of a covered domain.

The existence of covered domains for isometric immersions $y\in W^{1,2}$ is shown in [28, Corollary 1.2].

Proposition 1

Let $v\in W^{2,2}_{\rm sh}(\omega)$ and $x_{0}\in\omega$. There exists a neighbourhood $U$ of $x_{0}$ such that, if $\nabla^{2}v\neq 0$ a.e. in $U$, then for a suitable $\varepsilon>0$ there exist maps $\gamma\in W^{2,2}((-\varepsilon,\varepsilon);\mathbb{R}^{2})$ and $\lambda\in L^{2}((-\varepsilon,\varepsilon))$ such that $U\subset\{\gamma(t)+s\nu(t):s\in\mathbb{R},t\in(-\varepsilon,\varepsilon)\}$ and

[TABLE]

if $\gamma(t)+s\nu(t)\in U$.

Proof.

We may without loss of generality assume that $\omega$ is convex. Using [15, Theorem 10] take $v_{k}\in W^{2,2}\cap W^{1,\infty},S_{k}\subset\omega$ such that $x_{0}\in\operatorname{int}S_{k}$, $v_{k}=v$ on $S_{k}$ and $\|v_{k}\|_{1,\infty}\leqslant C$. By scaling $v_{k}$ with $\eta>0$ we can extend $\eta v_{k}$ to an isometry $y$ ([15, Theorem 7]) with $\eta v_{k}=y_{3}$. Then, because $y$ is an isometry:

[TABLE]

where $n=y_{,1}\wedge y_{,2}$ is the normal and $\operatorname{II}_{(y)}=(\nabla y)^{\top}\nabla n$ the second fundamental form of the surface $y(\omega)$. Since $\nabla^{2}y\neq 0$ a.e. near $x_{0}$, there is a neighbourhood $U$ of $x_{0}$ covered by some leading curve $\gamma$, that is: $U\subset\{\gamma(t)+s\nu(t):s\in\mathbb{R},t\in(-\varepsilon,\varepsilon)\}$ and, by [33, p. 111], on $U$ we have

[TABLE]

with $\tilde{\lambda}\in L^{2}$. Now, [19, Proposition 1, eq. (12)] shows that $\nabla y(\gamma(t)+s\nu(t))$ is independent of $s$, hence $n_{3}=(y_{,1}\wedge y_{,2})_{3}$ is also independent of $s$ and we can subsume it into the function $\tilde{\lambda}$. Setting $\lambda(t)=-n_{3}(t)\tilde{\lambda}(t)/\eta$ we obtain the representation (17). ∎

Finally, we come to:

**Proof of Theorem 2 ** To recapitulate, according to (5) and (14) the linearised Kirchhoff energy is given by

[TABLE]

for $v\in W^{2,2}_{\rm sh}$ (and $\infty$ otherwise).

We observe first that the set $\mathcal{N}=\operatorname{argmin}\{Q_{2}^{\ast}(F-F_{0}):F\in\mathbb{R}^{2\times 2}_{\operatorname{sym}},\det F=0\}$ is not empty because $F\mapsto Q_{2}^{\ast}(F-F_{0})$ is non-negative and strictly convex, but it also need not consist of just one point. Note next that $v$ is a minimiser of (18) iff $\nabla^{2}v(x)\in\mathcal{N}$ for almost every $x\in\omega$: On the one hand, every minimiser has finite energy and thus $\nabla^{2}v$ must be pointwise a.e. in the set $\{F\in\mathbb{R}^{2\times 2}_{\operatorname{sym}}:\det F=0\}$. On the other, any function $F:\omega\rightarrow\mathbb{R}^{2\times 2}_{\operatorname{sym}}$ with $F(x)\in\mathcal{N}$ a.e. minimises the integrand in (18) pointwise and thus the energy.

Next we show that any two elements $F,G$ of $\mathcal{N}$ are linearly independent. Indeed, by strict convexity we have for all $\lambda\in(0,1)$:

[TABLE]

Hence $\lambda F+(1-\lambda)G\not\in\mathcal{N}$ or else $F,G$ would not be minimisers. Because $Q_{2}^{\ast}$ attains a lower value here we must have $\det(\lambda F+(1-\lambda)G)\neq 0$. But then it cannot be that $G=\rho F$ for any scalar $\rho\in\mathbb{R}$ or else it would hold that $\det(\lambda F+(1-\lambda)G)=\det(\lambda F+(1-\lambda)\rho F)=C\det F=0$, a contradiction. Consequently, we have in particular $0\not\in\mathcal{N}$ unless $\mathcal{N}=\{0\}$. But in that case $\nabla^{2}v\equiv 0$ and the proof would be concluded.

Let now $v\in W^{2,2}_{\rm sh}$ be a minimiser for $\mathcal{I}_{\rm lKi}$. Note first that $\nabla v$ cannot be constant over open sets: indeed we just saw that w.l.o.g. $0\not\in\mathcal{N}$ and consequently the condition $\nabla^{2}v=0$ is excluded for a minimiser on any set of positive measure. Consider then some point $x_{0}\in\omega$ with a neighbourhood $U$ where $\nabla v$ is not constant and use the representation (17). We have that, pointwise a.e. and over $U$:

[TABLE]

If $\kappa(t)\neq 0$, by varying $s$ we obtain distinct, linearly dependent matrices $\nabla^{2}v(t,s)$. Because $\nabla^{2}v\in\mathcal{N}$ a.e., this shows that $\kappa(t)=0$ for a.e. $t$. As a consequence, $\gamma^{\prime}$ must be constant. But then $\lambda$ is also constant or again we would have points at which $\nabla^{2}v$ is linearly dependent. Since this holds locally around every $x=\gamma(t)+s\gamma^{\prime}(t)$, we deduce that $\nabla^{2}v$ is constant on $U$ and because we can cover $\omega$ in this manner, there exists $F\in\mathcal{N}$ such that $\nabla^{2}v\equiv F$ a.e. over $\omega$. $\Box$

4 Structure of minimisers for

$\mathcal{I}^{\theta}_{\rm vK}$ for small $\theta$

The second main contribution of this work is a first study of the properties of minimisers in the interpolating regime, “close” to the linearised von Kármán model. The results in Section 3 show that the transition from spherical to cylindrical shapes occurs in the interpolated von Kárma\a’an as the strength $\theta$ of the misfit increases. We will see that for small $\theta>0$ indeed there exists a unique stable branch of solutions emanating from a perfect spherical cap at $\theta=0$.

For the sake of clarity we restrict to the prototypical model from (8):

[TABLE]

Natural subsequent steps along this line of work, which we do not take here, are to consider the regime of large values of $\theta$ and to investigate the existence of the conjectured critical value $\theta_{c}$, as well as to consider the full model derived in (6).444In Section 5 we conduct numerical experiments supporting the conjecture that this critical value exists.

We recall that the existence of minimizers is guaranteed, cf. Remark 1. Without loss of generality wee assume that the barycenter of $\omega$ is [math]. So with $(f)_{\omega}:=\frac{1}{|\omega|}\int_{\omega}f(x)\,\mathrm{d}x$ for a function $f$ we in particular have $(x)_{\omega}=0$. In order to avoid ambiguities (and to apply Korn’s and Poincaré’s inequalities) we restrict the functions $w=(u,v)$ to lie in the Banach space

[TABLE]

with $X_{u},X_{v}$ as in

[TABLE]

and norm $\|(u,v)\|_{X}=(\|u\|_{1,2}^{2}+\|v\|_{2,2}^{2})^{1/2}$. By the arguments in [11, Remark 2] working with these spaces does not lead to a loss of generality either: For an affine function $g$, $\nabla(v+g)\otimes\nabla(v+g)-\nabla v\otimes\nabla v$ is a symmetrised gradient.

For small values of the parameter $\theta$ we have the following structural result on the set of minimizers showing the existence of a smooth branch of unique global minimisers. Let $v_{0}(x)=\tfrac{1}{2}|x|^{2}-c_{0}$ with $c_{0}=\frac{1}{2}(|x|^{2})_{\omega}$.

Theorem 4

There exists an $\varepsilon>0$, a unique point $u_{0}\in X_{u}$ and a uniquely determined $C^{1}$ map $\phi:[0,\varepsilon)\rightarrow X$ such that $\phi(0)=(u_{0},v_{0})$ and for each $\theta\in[0,\varepsilon)$:

[TABLE]

The proof is a direct consequence of Theorems 5 and 6 that are proved in the following two subsections. The main difficulty in obtaining a local branch of minimizers for $\theta\ll 1$ lies in the fact that minimisers at $\theta=0$ are not unique. Indeed,

[TABLE]

as can be readily checked. This is addressed in Subsection 4.1. The proof that in fact these minimisers are global is achieved by an application of a Taylor expansion for a carefully perturbed functional in Subsection 4.2.

4.1 A branch of solutions for

$\theta\ll 1$

**Notation **

In this section, the parameter $\theta$ will be explicitly included in the arguments of the functional and differentiation is understood to be with respect to the variables $w=(u,v)$, unless otherwise stated, i.e.

[TABLE]

We are interested in the existence and uniqueness of solutions $w=(u,v)$ to the equation

[TABLE]

as a function of $\theta\in[0,\varepsilon)$ with $\mathcal{I}^{\theta}_{\rm vK}$ given by (8). We will in fact prove the existence of a point $(u_{0},v_{0})\in X$ such that there exists a (locally) unique function $\phi(\theta)$, starting for $\theta=0$ at $(u_{0},v_{0})$, such that every $\phi(\theta)\in X$ is a critical point for $\mathcal{I}^{\theta}_{\rm vK}$. However, lack of uniqueness of minimisers at $\theta=0$, (19) will thwart what would be a natural application of the implicit function theorem. The problem manifests itself as a lack of injectivity of the first derivative at $(u,v)\in X$

[TABLE]

which for $\theta=0$ is

[TABLE]

and this vanishes at every $u\in X_{u}$ and the unique $v(x)=\tfrac{1}{2}|x|^{2}+a\cdot x+b$, $a\in\mathbb{R}^{2},b\in\mathbb{R}$, such that $(v)_{\omega}=0$ and $(\nabla v)_{\omega}=0$, i.e., $v=v_{0}$. Because of this the equation

[TABLE]

cannot be uniquely solvable for $(u,v)\in X$ as a function of $\theta$, even locally. Nevertheless, after some computations one can see that the problem is the presence of a leading factor $\theta$ which we can dispense with, because we may apply the implicit function theorem to the set of equivalent equations

[TABLE]

These equations are equivalent to $D\mathcal{I}^{\theta}_{\rm vK}(u,v;\theta)=0$ for any $\theta>0$ and by an application of the implicit function theorem around a specific point $(u_{0},v_{0};0)$ we determine the existence of a solution function $\phi:\Theta\rightarrow U\times V$ with $[0,\varepsilon)\subset\Theta,\varepsilon>0,U\times V\subset X$ open, $\phi(0)=(u_{0},v_{0})$ and $\left(\tfrac{1}{\theta}\partial_{u},\partial_{v}\right)\mathcal{I}^{\theta}_{\rm vK}(\phi(\theta);\theta)=0$. Then we have $D\mathcal{I}^{\theta}_{\rm vK}(\phi(\theta);\theta)=0$ for $\theta>0$ because of the equivalence mentioned and $D\mathcal{I}^{\theta}_{\rm vK}(\phi(0);0)=0$ by the choice of $(u_{0},v_{0})$.

Theorem 5

There exists an open set $W$ in $X$, an $\varepsilon>0$, a point $u_{0}\in X_{u}$ such that $w_{0}=(u_{0},v_{0})\in W$ and a uniquely determined $C^{1}$ map $\phi:\Theta\rightarrow W$ such that $\phi(0)=w_{0}$ and

[TABLE]

for all $w\in W$ and $\theta\in[0,\varepsilon)$.

Proof.

We first define a new set of equations to solve, then show that the second derivative of $\mathcal{I}^{\theta}_{\rm vK}$ is one to one and then the conclusion is exactly that of the implicit function theorem. For brevity we write

[TABLE]

These define a scalar product and a norm in $L^{2}(\omega;\mathbb{R}^{2\times 2}_{\operatorname{sym}})$ since $Q_{2}$ is by construction bilinear and symmetric and it is positive definite on this space. Even though $Q_{2}$ vanishes on antisymmetric matrices, during the proof we keep track of symmetrised arguments to these functions for the sake of clarity.

Step 1: Equivalent equations.

From the computations leading to (20) we have:

[TABLE]

and

[TABLE]

for all $(\varphi,\psi)\in X$. We observe first that, because $\left(\frac{1}{\theta}\partial_{u}\right)\mathcal{I}^{\theta}_{\rm vK}$ is independent of $\theta$ the right hand side makes sense even if $\theta=0$. Now, on the one hand, for any fixed value of $\theta\geqslant 0$ solving the system

[TABLE]

implies solving:

[TABLE]

where $f:X\times\mathbb{R}\rightarrow\mathcal{L}(X,\mathbb{R})$ is given by

[TABLE]

On the other hand, solving $f(u,v;\theta)=0$ for $\theta>0$ is equivalent to solving the original problem $D\mathcal{I}^{\theta}_{\rm vK}(u,v;\theta)=0$ as we desired.

Step 2: A zero and the derivative of $f$.

Since we are interested in the behaviour around $\theta=0$, we evaluate here and obtain

[TABLE]

We can compute a zero of $f(\cdot,\cdot;0)$ by first considering the last term, which vanishes for all $\psi\in X_{v}$ if and only if $v=v_{0}$. We next observe that the first term encodes the orthogonality of $\nabla_{s}u+\tfrac{1}{2}\nabla v_{0}\otimes\nabla v_{0}$ to the space of symmetrised gradients $\operatorname{SG}_{u}\coloneqq\left\{\nabla_{s}\varphi:\varphi\in X_{u}\right\}$ with respect to the scalar product induced by $Q_{2}$. The $u\in X_{u}$ realizing this is attained by projecting onto $\operatorname{SG}_{u}$, i.e.

[TABLE]

where $\pi:L^{2}(\omega;\mathbb{R}^{2\times 2}_{\operatorname{sym}})\rightarrow L^{2}(\omega;\mathbb{R}^{2\times 2}_{\operatorname{sym}})$ is the orthogonal projection onto $\operatorname{SG}_{u}$ given by

[TABLE]

By the Korn-Poincaré inequality this determines $u_{0}\in X_{u}$ uniquely. We have then a point $w_{0}=(u_{0},v_{0})$ such that

[TABLE]

Finally, we compute $\frac{\,\mathrm{d}}{\,\mathrm{d}\varepsilon}|_{\varepsilon=0}f(u_{0}+\varepsilon\varphi_{2},v_{0}+\varepsilon\psi_{2};0)[\varphi_{1},\psi_{1}]$ to have the derivative of $f$:

[TABLE]

Step 3: The map $F:X\rightarrow\mathcal{L}(X,\mathbb{R})$ is an isomorphism.

Note first that the map

[TABLE]

defines a scalar product in $X$, with positive-definiteness following from Korn-Poincaré’s and Poincaré’s inequality. Then we can write $F$ as

[TABLE]

where we defined $\tilde{\pi}\coloneqq\nabla_{s}^{-1}\circ\pi$, a continuous map from $L^{2}(\omega;\mathbb{R}^{2\times 2}_{\operatorname{sym}})$ to $X_{u}$. The Riesz representation for $F(\varphi_{2},\psi_{2})$ in $\mathcal{L}(X,\mathbb{R})$ is then $(\varphi_{2}+\tilde{\pi}((\nabla v_{0}\otimes\nabla\psi_{2})_{\operatorname{sym}}),\tfrac{1}{12}\psi_{2})$ and the map

[TABLE]

is clearly an isomorphism in $X$, with continuity for $\psi_{2}\mapsto\tilde{\pi}((\nabla v_{0}\otimes\nabla\psi_{2})_{\operatorname{sym}})$ following from the continuity of $\tilde{\pi}$ and the Sobolev embedding $W^{1,2}\hookrightarrow L^{4}$. ∎

4.2 Uniqueness and globality of minimisers

In addition to the previous local result, we can prove that the critical points found in the previous subsection are the unique global minimizers for small non zero values of the parameter $\theta$. We do this in two steps: close to the origin $(u_{0},v_{0})$ of the branch of solutions, we would like to perform a Taylor expansion and use that the second differential at $(u_{0},v_{0})$ is “almost” positive definite.

The key idea is to slightly modify the energy by a shift and a rescaling in order to obtain derivatives as those appearing in the equivalent equations (22) of Theorem 5, thus obtaining a positive definite second derivative. We set

[TABLE]

and then $(\tilde{u}_{\theta},\tilde{v}_{\theta})$ is a minimiser of $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$ if and only if $(u_{0}+\tilde{u}_{\theta}/\theta,\tilde{v}_{\theta})$ is a minimiser of $\mathcal{I}^{\theta}_{\rm vK}$. In other words, if $(u_{\theta},v_{\theta})$ is a minimiser of $\mathcal{I}^{\theta}_{\rm vK}$, then $\tilde{u}_{\theta}=\theta(u_{\theta}-u_{0})$ and $\tilde{v}_{\theta}=v_{\theta}$ minimise $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$.

We name $\tilde{w}_{0}$ the point around which we investigate the modified functional:

[TABLE]

Theorem 6

There exists $\theta_{c}>0$ and a neighborhood $\tilde{W}\subset X$ with $\tilde{w}_{0}\in\tilde{W}$ such that for every $\theta\in(0,\theta_{c})$, every critical point of $D\tilde{\mathcal{I}}^{\theta}_{\rm vK}$ is the unique global minimiser of $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$.

Proof.

We proceed in three steps. First we prove that there is some $\theta_{c}>0$ such that $D^{2}\tilde{\mathcal{I}}^{\theta}_{\rm vK}(\tilde{w})$ is positive definite for all $\theta\in(0,\theta_{c})$ if $\|\tilde{w}-\tilde{w}_{0}\|\char 60\relax\eta$ for some suitable $\eta>0$ and $\tilde{w}_{0}=(0,v_{0})$ as defined in (23). Then we use this to determine a neighbourhood of $\tilde{w}_{0}$ where (local) minimisers of $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$ will be global by first considering points close to one such minimiser and finally those far away. We will need the first two derivatives of $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$.

For the first differential we apply the chain rule to obtain $D_{u}\tilde{\mathcal{I}}^{\theta}_{\rm vK}(\tilde{u},\tilde{v})=\frac{1}{\theta}D_{u}\mathcal{I}^{\theta}_{\rm vK}\left(u_{0}+\frac{\tilde{u}}{\theta},\tilde{v}\right)$ and substitute:

[TABLE]

For the second differential we can compute another directional derivative:

[TABLE]

Step 1: Local positive definiteness.

We show there exist $\eta>0$ and $\theta_{c}>0$ s.t. $D^{2}\tilde{\mathcal{I}}^{\theta}_{\rm vK}(\tilde{w})$ is positive definite for all $\theta\char 60\relax\theta_{c}$ and all $\|\tilde{w}-\tilde{w}_{0}\|_{X}\char 60\relax\eta$. More precisely, we even show that there exists some $\bar{c}>0$ such that

[TABLE]

for all $\theta\char 60\relax\theta_{c}$, $\|\tilde{w}-\tilde{w}_{0}\|_{X}\leq\eta$ and $(\varphi,\psi)\in X$.

Let then $\eta>0$ be fixed and to be determined later and let $\tilde{w}=(\tilde{u},\tilde{v})\in X$ with $\|\tilde{w}-\tilde{w}_{0}\|_{X}\char 60\relax\eta$. We start by  bringing terms together in (24):

[TABLE]

Given $f,g\in W^{1,2}(\omega;\mathbb{R}^{2})$ we have, by the bounds (2) for $Q_{2}$ and Hölder (with the Sobolev embedding $W^{1,2}(\omega;\mathbb{R}^{2})\hookrightarrow L^{4}(\omega;\mathbb{R}^{2})$):

[TABLE]

Using this, the first and last term above can be estimated using Korn-Poincaré and Poincaré’s inequality:

[TABLE]

for constants $c_{1},C_{1},\tilde{C}_{1}>0$, where in the last step we used the assumption $\|\tilde{v}-v_{0}\|_{2,2}\char 60\relax\eta$ to bound $\|\tilde{v}\|_{2,2}^{2}$ by some constant independent of $\eta\leqslant 1$. For the second term, use Cauchy-Schwarz for $Q_{2}$, and the same ideas as above:

[TABLE]

Again, we used that by assumption $\|\tilde{u}\|_{1,2}\char 60\relax\eta$ and $\|\tilde{v}-v_{0}\|_{2,2}\char 60\relax\eta$.

Finally we estimate the third term in $D^{2}\tilde{\mathcal{I}}^{\theta}_{\rm vK}$ with analogous arguments and obtain $(c)\geqslant c_{2}\|\psi\|_{2,2}^{2}$, for a $c_{2}>0$. Bringing the previous computations together, with a $C_{2}>0$ we have:

[TABLE]

from which (25) follows if $\theta_{c}$ and $\eta$ are chosen sufficiently small.

From now on, we let $\tilde{w}_{\theta}=(\tilde{u}_{\theta},\tilde{v}_{\theta})$ be a critical point of $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$ with

[TABLE]

and we prove that it is in fact the unique global minimizer.

Step 2: Estimates close to $\tilde{w}_{\theta}$.

Consider first some $\tilde{w}\in X$ which is close to $\tilde{w}_{\theta}$:

[TABLE]

With a Taylor expansion and (25) we see:

[TABLE]

where $z\in\{\alpha\tilde{w}+(1-\alpha)\tilde{w}_{\theta}:\alpha\in[0,1]\}\subset B_{\eta}(\tilde{w}_{0})$ by (26) and (27). So

[TABLE]

Step 3: Estimates far away from $\tilde{w}_{\theta}$.

Consider now any $\tilde{w}\in X$ with

[TABLE]

which by (26) implies that $\|\tilde{w}-\tilde{w}_{0}\|_{X}>\eta/3$. We consider two cases:

Case 1: $\|\tilde{v}-v_{0}\|_{2,2}\geqslant\eta/6$: We discard the first term in the energy, recall that $v_{0}(x)=|x|^{2}/2-c_{0}$ and use the lower bound for $Q_{2}$ in (2) and Poincaré’s inequality:

[TABLE]

for a $c_{1}>0$. To compare this with the energy at $\tilde{w}_{0}$ we add and subtract $\tilde{\mathcal{I}}^{\theta}_{\rm vK}(\tilde{w}_{0})=\frac{\theta}{2}\langle\nabla_{s}u_{0}+\tfrac{1}{2}\nabla v_{0}\otimes\nabla v_{0}\rangle$:

[TABLE]

where the last line is due to the fact that $\tilde{w}_{\theta}$ minimises $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$ over the ball $B_{\frac{2}{3}\eta}(\tilde{w}_{\theta})$.

Case 2: $\|\tilde{v}-v_{0}\|_{2,2}\char 60\relax\eta/6$: In this case we also have $\|\tilde{u}\|_{1,2}\geq\eta/6$ by (26) and (28). We can estimate the energy for $\tilde{w}$ as follows:

[TABLE]

where we used the Cauchy-Schwarz inequality with $\varepsilon\coloneqq\frac{1}{4\theta}$. Both terms may be estimated once again by a combination of the bounds (2) for $Q_{2}$, Sobolev’s embedding $W^{1,2}(\omega)\hookrightarrow L^{4}(\omega)$ and Poincaré’s inequality:

[TABLE]

and

[TABLE]

since $\|\tilde{v}-v_{0}\|_{2,2}^{2}\char 60\relax\eta/6$. Now plug this back into the previous estimate and insert

[TABLE]

to obtain

[TABLE]

As above, the last line holds because $\tilde{w}_{\theta}$ minimises $\tilde{\mathcal{I}}^{\theta}_{\rm vK}$ in a $\frac{2}{3}\eta$-neighbourhood of itself. ∎

5 Discretisation of the interpolating theory

Our goal in this section is to study the qualitative behaviour of minimisers in the interpolating regime $\alpha=3$. To this end, we develop a simple numerical method to approximate minimisers and prove $\Gamma$-convergence to the continuous problem. Numerical computations are then conducted for the prototypical example from (8). We experimentally evaluate the conjectured existence of a critical value $\theta_{c}>0$ for which the symmetry of minimisers is “strongly” broken. We will not provide a full theoretical analysis, but instead adduce some empirical evidence to support the claim.

As can only be expected from a topic originating in structural mechanics, numerical methods for plate models are a vast field with a long history and as such a comprehensive review falls well beyond the scope of this contribution. However, it can be said that a significant portion of finite element approaches focus on the Euler-Lagrange equations. For von Kármán-like theories like our interpolating regime, these are transformed into an equivalent form in terms of the Airy stress function [20, §2.6.2]. The resulting system of equations is of fourth order and can be solved with conforming $C^{1}$ elements like Argyris or specifically taylored ones. To avoid the higher number of degrees of freedom, non-conforming methods can be used instead,555See [23, 24] for particular instances of a conforming and a non-conforming method respectively, as well as reviews of recent literature. but a poor choice of the discretisation can suffer from locking, as briefly described in Remark 4. Some successful classical methods employ $C^{0}$ Discrete Kirchhoff triangles (DKT), but it is also possible to employ standard Lagrange elements with penalty methods [8], as we will do.

A recent line of work, upon which we heavily build in this section, is that of [4, 6], where the author develops discrete gradient flows for the direct computation of (local) minimisers of non-linear Kirchhoff and von Kármán models. $\Gamma$-convergence and compactness results are also proved showing the convergence of the discrete energies to the continuous ones, as well as their respective minimisers.666For a concise introduction to $\Gamma$-convergence for Galerkin discretisations and quadrature approximations of energy functionals, see [26]. Crucially, these papers use DKTs for the discretisation of the out-of-plane displacements, allowing for a representation of derivatives at nodes in the mesh which is decoupled from function values. This enables e.g. the imposition of an isometry constraint for the non-linear Kirchhoff model, but also the computation of a discrete gradient $\nabla_{\varepsilon}$ projecting the true gradient $\nabla v_{\varepsilon}$ of a discrete function $v_{\varepsilon}$ into a standard piecewise $P_{2}$ space. The operator $\nabla_{\varepsilon}$ has good interpolation properties circumventing the lack of $C^{1}$ smoothness of DKTs which would otherwise make them unsuitable to approximate solutions in $H^{2}$. We refer to the book [5] for a systematic and mostly self-contained introduction to these methods.

5.1 Discretisation

We wish to investigate minimal energy configurations of the following functional:

[TABLE]

where $(u,v)\in W^{1,2}(\omega;\mathbb{R}^{2})\times W^{2,2}(\omega;\mathbb{R}^{2})$, cf. (6). We recall the representation of $\overline{Q}_{2}$ derived in (12), which in particular shows that $\overline{Q}_{2}$ is a strictly convex polynomial of degree $2$ on $\mathbb{R}^{2\times 2}_{\operatorname{sym}}\times\mathbb{R}^{2\times 2}_{\operatorname{sym}}$. It is extended to a convex quadratic function on $\mathbb{R}^{2\times 2}\times\mathbb{R}^{2\times 2}$ by our setting

[TABLE]

for $F,G\in\mathbb{R}^{2\times 2}$. We assume that $\omega\subset\mathbb{R}$ is a bounded simply connected domain with Lipschitz boundary and barycenter [math]. We implement (projected) gradient descent in a non-conforming method using $C^{0}$ linear Lagrange elements. The first step is to transform the problem into one of constrained minimisation reducing the order of the elements required.

Problem 1

Find minimisers of

[TABLE]

with $u,z\in W^{1,2}(\omega;\mathbb{R}^{2})$ and

[TABLE]

If $z\not\in Z$, then we set $J^{\theta}(u,z)=+\infty$.

Note that our assumptions on $\omega$ guarantee that $Z=\{\nabla v:v\in W^{2,2}(\omega)\}$. We can now use $H^{1}$-conforming elements but, for simplicity of implementation, instead of adding the constraint into the discrete spaces to obtain a truly conforming discretisation, we add a penalty term $\mu_{\varepsilon}\|\operatorname{curl}z_{\varepsilon}\|^{2}$ to ensure that the solutions $z_{\varepsilon}$ are close to gradients.

Assume from now on that $\omega$ is a polygonal domain. For fixed $\varepsilon>0$, introduce a quasi-uniform triangulation $\mathcal{T}_{\varepsilon}$ of $\omega$ with triangles $T$ of uniformly bounded diameter  $c^{-1}\varepsilon\leqslant\varepsilon_{T}\leqslant c\varepsilon$ for some $c>0$ and all $\varepsilon>0$ and $T\in\mathcal{T}_{\varepsilon}$.777 Note that this does not allow for arbitrary local refinements or grading (a different scaling of simplices along different directions as $\varepsilon\rightarrow 0$), but the fact that this is not optimal is not of concern here. Such a mesh is in particular said to be, in virtue of the uniform upper bound, shape-regular. We denote by $\mathcal{N}_{\varepsilon}$ the set of all nodes of the triangulation. Define $V_{\varepsilon}$ to be the standard piecewise affine, globally continuous Lagrange $P_{1}$ finite element space $\mathcal{S}^{1}(\mathcal{T}_{\varepsilon})$ in two dimensions:

[TABLE]

Quadrature rules will be chosen to be exact for this polynomial degree and the first integrand in the energy interpolated for this to apply by means of the interpolated quadratic function

[TABLE]

This is defined (with a slight abuse of notation) component-wise using the element-wise nodal interpolant $\hat{I}_{\varepsilon}$, defined for functions $v\in L^{\infty}(\omega)$ such that $v_{|T}\in C(\overline{T})$ for all $T\in\mathcal{T}_{\varepsilon}$ as

[TABLE]

where $\varphi_{z|T}$ is the truncation by zero outside $T$ of the global basis function $\varphi_{z}\in\mathcal{S}^{1}$. Because this is a linear combination of truncated global basis functions, the range of $\hat{I}_{\varepsilon}$ is the space $\hat{\mathcal{S}}^{1}(\mathcal{T}_{\varepsilon})$ of discontinuous, piecewise affine Lagrange elements.

In cases where the function to be interpolated is continuous, the element-wise nodal interpolant coincides with the standard nodal interpolant into the space $\mathcal{S}^{1}$ of globally continuous, piecewise affine functions, which is defined as

[TABLE]

Notice that the shape functions $\varphi_{z}$ are not truncated. In order to control the error incurred by the interpolation. When working with discontinuous functions in $\hat{\mathcal{S}}^{1}$, we will  use the following local result. This follows from standard nodal interpolation estimates (see e.g. [17, Theorem 4.28] or [9, (4.4.4)])

[TABLE]

or can be shown directly, e.g. in [5, Proposition 3.1].

Lemma 1** **(Local interpolation estimate)

Let $T\in\mathcal{T}_{\varepsilon}$ and $v\in C^{1}(\overline{T})$. If $\hat{I}_{\varepsilon}$ is the element-wise nodal interpolant (30), then

[TABLE]

The goal is to solve:

Problem 2

Let $\mu_{\varepsilon}>0$. Compute minimisers of the discrete energy

[TABLE]

for $(u_{\varepsilon},z_{\varepsilon})\in V_{\varepsilon}^{2}$. (As usual, if $(u_{\varepsilon},z_{\varepsilon})\in W^{1,2}(\omega;\mathbb{R}^{2})^{2}\backslash V^{2}_{\varepsilon}$, we set $J^{\theta}_{\varepsilon}(u_{\varepsilon},z_{\varepsilon})=+\infty$.)

Remark 3** **(Scaling of the constants)

The penalty $\mu_{\varepsilon}=\mu(\varepsilon)$ needs to explode as $\varepsilon\rightarrow 0$ in order for the functionals to $\Gamma$-converge (Theorem 7). However, large penalties negatively affect the condition number of the system, so that an adequate choice for $\mu_{\varepsilon}$, dependent on the mesh size $\varepsilon$, is required [17, p.416]. We have not explictly investigated how this requirement interacts with the $\Gamma$-convergence of the functionals, but in our proof we require only that $\mu_{\varepsilon}\rightarrow\infty$ not faster than $\varepsilon^{-2}$. In the implementation we use $\mu_{\varepsilon}=\varepsilon^{-1/2}$. Analogously, large values of the Lamé constants have a similar effect and therefore hinder convergence, so one needs to scale them to the order of the problem.

Remark 4** **(Common issues with FEM for plates)

Discretisations for lower dimensional theories can face complications due to the infamous locking phenomena. In a nutshell, these mean that as the thickness of the plate tends to zero, discrete solutions “lock” to stiff states of lower, or even vanishing, bending or shearing than the analytic ones.888We refer to [3] for a first rigorous definition of locking, to [29, Chapters 5 and 6] for detailed computations highlighting the issues with linear elements in the context of Timoshenko beams and to the thesis [30] for a thorough and detailed analysis of locking in shell models. Another instance of unexpected behaviour is known as the Babuška paradox [2], again a failure to converge as expected, which can happen in e.g. the Kirchhoff model when both vertical and tangential displacements are fixed at the boundaries of a polygonal domain: these so-called “hard” support constraints are not enforced in the same manner as in the continuous model because of the approximated domain.

There are two potential sources of locking in our setting: the penalty term $\mu_{\epsilon}$, which is akin to the shear strain in Timoshenko beams, and $\theta$. We have not obtained any a priori bounds on the error in this work, but a rigorous treatment of the problem would require estimates which are uniform in these parameters as the mesh diameter goes to zero. For the regimes studied and the geometries considered we have found the issue to be of moderate practical relevance, but it does manifest itself e.g. with more complicated domains or higher values of $\theta$.

Finally, our simulations will not suffer from Babuška’s paradox because we do not prescribe boundary conditions.

5.2 $\Gamma$-convergence of the discrete energies

The first step in the proof that $J^{\theta}_{\varepsilon}\overset{\Gamma}{\rightarrow}J^{\theta}$ is dispensing with the interpolation operators for numerical integration: due to the good properties of $\hat{I}_{\varepsilon}$, we can assume that we work with the true integrals $\int\overline{Q}_{2}$ instead of $\int\overline{Q}_{2}^{\varepsilon}$:

Lemma 2** **(Numerical integration)

Let $u_{\varepsilon},z_{\varepsilon}\in W^{1,2}(\omega;\mathbb{R}^{2})$ be uniformly bounded in $W^{1,2}$ and let $Q_{2}^{\varepsilon}=\hat{I}_{\varepsilon}\circ Q_{2}$ as above. Let A_{\varepsilon}\coloneqq\big{(}\theta^{1/2}(\nabla_{s}u_{\varepsilon}+\tfrac{1}{2}z_{\varepsilon}\otimes z_{\varepsilon}),-\nabla z_{\varepsilon}\big{)}. Then, as $\varepsilon\rightarrow 0$:

[TABLE]

Proof.

By the local interpolation estimate Lemma 1:

[TABLE]

Now, the first term is simply $\|1+|A_{\varepsilon}|\|_{0,2,\omega}\leq|\omega|^{1/2}+\|A_{\varepsilon}\|_{0,2,\omega}$ which is uniformly bounded since $\|z_{\varepsilon}\otimes z_{\varepsilon}\|_{0,2}=\|z_{\varepsilon}\|_{0,4}^{2}\lesssim\|z_{\varepsilon}\|_{1,2}^{2}$, and for the second we use that both $\nabla_{s}u_{\varepsilon}$ and $\nabla z_{\varepsilon}$ are piecewise constant so that for $i=1,2$,

[TABLE]

and

[TABLE]

A standard inverse estimate (see e.g. [9, Theorem 4.5.11]) provides the bound

[TABLE]

We plug this into the preceding computation to obtain

[TABLE]

The last two norms being uniformly bounded, we conclude:

[TABLE]

∎

The second step is, as usual, to ensure that we can focus on smooth functions for simplicity in the construction of the upper bound:

Lemma 3

The set $C^{\infty}(\overline{\omega},\mathbb{R}^{2})\cap Z$ is $W^{1,2}$-dense in $Z$.

Proof.

This follows from $Z=\{\nabla v:v\in W^{2,2}(\omega)\}$ and the density of $C^{\infty}(\overline{\omega})$ in $W^{2,2}(\omega)$. ∎

Theorem 7

Let $J^{\theta},J^{\theta}_{\varepsilon}$ be given by (29) and (32) respectively. Assume that $\mu_{\varepsilon}\to\infty$ such that $\mu_{\varepsilon}=o(\varepsilon^{-2})$ as $\varepsilon\rightarrow 0$. Then $J^{\theta}_{\varepsilon}\overset{\Gamma}{\rightarrow}J^{\theta}$ as $\varepsilon\rightarrow 0$ with respect to weak convergence in $W^{1,2}$.

Proof.

Because of Lemma 2 we can substitute $\overline{Q}_{2}$ for $\overline{Q}_{2}^{\varepsilon}$ in $J^{\theta}_{\varepsilon}$. Also, by Lemma 3 it is enough to consider smooth functions for the upper bound. Set

[TABLE]

Step 1: Upper bound.

Let $(u,z)\in W^{1,2}(\omega;\mathbb{R}^{2})\times Z$ be $C^{\infty}$ up to the boundary and define $u_{\varepsilon}\coloneqq I_{\varepsilon}(u),z_{\varepsilon}\coloneqq I_{\varepsilon}(z)$, where $I_{\varepsilon}$ is the nodal interpolant of (31). Note that because $u$ and $z$ are smooth, we can apply standard interpolation estimates to show strong convergence in $W^{1,2}$ of these sequences towards $u$ and $z$. By the compact Sobolev embedding $W^{1,2}\hookrightarrow L^{4}$ we have $z_{\varepsilon}\rightarrow z$ in $L^{4}$, and $z_{\varepsilon}\otimes z_{\varepsilon}\rightarrow z\otimes z$ in $L^{2}$, so we have that $A_{\varepsilon}\rightarrow A$ in $L^{2}$. Since $\overline{Q}_{2}$ is a polynomial of degree 2, this implies

[TABLE]

as $\varepsilon\rightarrow 0$. By the same interpolation estimate above and the assumption on $\mu_{\varepsilon}$ we have that $\mu_{\varepsilon}\|\operatorname{curl}(\hat{I}_{\varepsilon}(z)-z)\|^{2}_{0,2}=o(1)$ as $\varepsilon\rightarrow 0$, and consequently

[TABLE]

Step 2: Lower bound.

Let $u_{\varepsilon},z_{\varepsilon}\in V_{\varepsilon}\subset W^{1,2}$ with $u_{\varepsilon}\rightharpoonup u$, and $z_{\varepsilon}\rightharpoonup z$ weakly in $W^{1,2}$ to $u\in W^{1,2}(\omega;\mathbb{R}^{2}),z\in Z$. Because $z_{\varepsilon}\otimes z_{\varepsilon}\rightarrow z\otimes z$ in $L^{2}$, we have that $A_{\varepsilon}\rightharpoonup A$ in $L^{2}$. Moreover, $\operatorname{curl}z_{\varepsilon}\rightharpoonup\operatorname{curl}z$. If $\underset{\varepsilon\rightarrow 0}{\operatorname{linf}}J^{\theta}_{\varepsilon}=\infty$, the assertion is trivial. If not, then $\mu_{\varepsilon_{k}}\int_{\omega}|\operatorname{curl}z_{\varepsilon_{k}}|^{2}\,\mathrm{d}x\leqslant C$ and $\|\operatorname{curl}z_{\varepsilon_{k}}\|_{0,2}\rightarrow 0$ for a subsequence $\varepsilon_{k}\rightarrow 0$. But then $\operatorname{curl}z=0$. Dropping the (non-negative) $\operatorname{curl}$ term in $J^{\theta}_{\varepsilon}$ and by the weak sequential lower semicontinuity of all integrands involved ($\overline{Q}_{2}$ being a convex quadratic function), we then get

[TABLE]

∎

The final ingredient of this subsection is a proof that sequences with bounded energy are (weakly) precompact. The fundamental theorem of $\Gamma$-convergence then shows convergence of global minimisers. In order for this to work, we need to assume conditions in the space which provide Korn and Poincaré inequalities. We can do this using functions with zero mean, zero mean of the gradient or zero mean of the antisymmetric gradient as we do above, but including these conditions in the discrete spaces is not entirely trivial. Because the energies are invariant under the transformations which are factored out by taking quotient spaces as described in the sections mentioned, it is enough for our purposes to claim compactness modulo these transformations and to exclude them in the implementation via projected gradient descent.

Theorem 8** **(Compactness)

Let $(u_{\varepsilon},z_{\varepsilon})_{\varepsilon>0}$ be a sequence in $(V_{\varepsilon}\cap X_{u})^{2}$ with bounded energy. Then there exist $u\in W^{1,2},z\in Z$ such that $u_{\varepsilon}\rightharpoonup u$ and $z_{\varepsilon}\rightharpoonup z$. in $W^{1,2}$.

Proof.

As above, let A_{\varepsilon}\coloneqq\big{(}\theta^{1/2}(\nabla_{s}u_{\varepsilon}+\frac{1}{2}z_{\varepsilon}\otimes z_{\varepsilon}),-\nabla z_{\varepsilon}\big{)}. Note that we cannot use Lemma 2 to substitute $Q_{2}$ for $Q_{2}^{\varepsilon}$ since we do not have uniform bounds in $W^{1,2}$ by assumption, so we work directly with $J^{\theta}_{\varepsilon}$.

We begin by observing that, as $\overline{Q}_{2}:\mathbb{R}^{2\times 2}\times\mathbb{R}^{2\times 2}_{\operatorname{sym}}\to\mathbb{R}$ is a convex quadratic function bounded from below which is strictly convex on $\mathbb{R}^{2\times 2}_{\operatorname{sym}}\times\mathbb{R}^{2\times 2}_{\operatorname{sym}}$, there are constants $\bar{c},\bar{C}>0$ such that

[TABLE]

for all $E,F\in\mathbb{R}^{2\times 2}$. In particular,

[TABLE]

and consequently, by Poincaré’s inequality:

[TABLE]

We have then a subsequence (not relabeled) weakly converging in $W^{1,2}$ to some $z\in W^{1,2}$. In particular $\nabla z_{\varepsilon}\rightharpoonup\nabla z$ and $\operatorname{curl}z_{\varepsilon}\rightharpoonup\operatorname{curl}z$ in $L^{2}$. But also

[TABLE]

and therefore $\operatorname{curl}z=0$, i.e. $z\in Z$.

Now, for the sequence $u_{\varepsilon}$ we must work with $\overline{Q}_{2}^{\varepsilon}$ instead. First write

[TABLE]

and thus

[TABLE]

Since this applies pointwise, after (local) interpolation the estimate still holds:

[TABLE]

where in the firs step we have used that $\nabla_{s}u_{\varepsilon}$ is piecewise constant. So

[TABLE]

We claim now that $\|\hat{I}_{\varepsilon}(|z_{\varepsilon}|^{4})-|z_{\varepsilon}|^{4}\|_{0,1}=\mathcal{O}(\varepsilon)$. Indeed, by the local interpolation estimate (Lemma 1) and Hölder’s inequality for integrals and for sums:

[TABLE]

and this goes to zero as $\varepsilon\rightarrow 0$ by (33). But then $\int_{\omega}\hat{I}_{\varepsilon}(|z_{\varepsilon}|^{4})\leqslant C$ and by Korn-Poincaré’s inequality, the Sobolev embedding $W^{1,2}\hookrightarrow L^{4}$ and the previous bound, we have

[TABLE]

The sequence $(u_{\varepsilon})_{\varepsilon>0}$ is therefore also weakly precompact in $W^{1,2}(\omega;\mathbb{R}^{2})$ and the proof is complete. ∎

5.3 Discrete gradient flow

As a concrete example we specialize now to the prototypical example

[TABLE]

cf. (8). For each discrete problem, we compute local minimisers using gradient descent, for which the basic result is the following (see [5, §4.3.1]):

Theorem 9** **(Projected gradient descent)

Let $V_{\varepsilon}$ and $J^{\theta}_{\varepsilon}$ be given as in Problem 2 and let $(\cdot,\cdot)$ be the scalar product on $V_{\varepsilon}$. The map $F_{\varepsilon}:V_{\varepsilon}\times V_{\varepsilon}\rightarrow(V_{\varepsilon}\times V_{\varepsilon})^{\prime}$ given by

[TABLE]

is the Fréchet derivative of $J^{\theta}_{\varepsilon}$. Let $\pi_{u}:V_{\varepsilon}^{2}\rightarrow(V_{\varepsilon}\cap X_{u})^{2}$ be the linear orthogonal projection onto its image. The sequence defined as

[TABLE]

with $w_{\varepsilon}^{0}=(u_{\varepsilon}^{0},v_{\varepsilon}^{0})\in(V_{\varepsilon}\cap X_{u})^{2}$ and $d_{\varepsilon}^{j}\in V_{\varepsilon}\times V_{\varepsilon}$ such that

[TABLE]

and $\alpha_{j}$ determined with line search is energy decreasing. A line search means computing the maximal $\alpha_{j}\in\{2^{-k}:k\in\mathbb{N}\}$ such that

[TABLE]

where $\rho\in(0,1/2)$ is the proverbial fudge factor.

Proof.

The computation of $F^{\theta}_{\varepsilon}$ is straightforward. To see that the iteration is energy decreasing use (35) and the self-adjointness of $\pi_{u}=\pi^{2}_{u}$ to compute

[TABLE]

The existence of $\alpha_{j}>0$ is guaranteed as long as $J^{\theta}_{\varepsilon}\in C^{2}(V_{\varepsilon}^{2})$ because then we can perform a Taylor expansion and use again (35):

[TABLE]

∎

Remark 5** **(Caveat: local and global minimisers)

Even though we now know that the discrete energies correctly approximate the continuous one, as well as any global minimisers, gradient descent on each discrete problem is only guaranteed to converge to some local minimiser $w^{\star}_{\varepsilon}$. Lacking some means of tracking a particular $w^{\star}_{\varepsilon}$ as $\varepsilon\rightarrow 0$, there is not much one can do to prove that our method actually approximates the true global minimisers of $\mathcal{I}_{\operatorname{vK}}^{\theta}$. Unless $\theta\ll 1$, in which case we know local minimisers to be global (cf. Theorem 6).

5.4 Experimental results

For the implementation of the discretisation detailed above, we employ the FEniCS library [1] in its version 2017.1.0. The code is available at [10] and includes the model, parallel execution, experiment tracking using Sacred [16] with MongoDB as a backend and exploration of results with Jupyter [22] notebooks, Omniboard [35] and a custom application. Everything is packaged using docker-compose for simple reproduction of the results and one-line deployment.

We set $\omega=\hat{B}_{1}(0)$, a (coarse) polygonal approximation of the unit disc and test several initial conditions. The space $V_{\varepsilon}$ has $\sim$7000 dofs. We implement a general $Q_{2}$ for isotropic homogeneous material with the two (scaled) Lamé constants set to those of steel at standard conditions. We apply neither body forces nor boundary conditions, but hold one interior cell to fix the value of the free constants. We compute minimisers for increasing values of $\theta$ and $\mu_{\varepsilon}\sim 1/\sqrt{\varepsilon}$ via projected gradient descent (onto the space of admissible functions $V_{\varepsilon}\cap X_{u}$) and examine the symmetry of the final solution. The choice $\varepsilon^{-1/2}$ has shown to provide the fastest convergence results while keeping the violation of the constraint in the order of $10^{-4}$ (higher penalties have the expected effect of adversely affecting convergence). We track two magnitudes as measures of symmetry: on the one hand we compute the mean bending strain over the domain and on the other, as a second simple proxy we employ the quotient of the lengths of the principal axes.

The first initial configuration is the trivial deformation $y^{0}_{\varepsilon}=0$. Note that because the model is prestrained, the ground state is non-trivial and the plate “wants” to reach a lower energy state. In Figure 2 we depict the results of running the energy minimisation procedure for multiple values of $\theta$.

We further highlight the behaviour of the solution as a function of $\theta$ in Figure 3. In the first plot we compute the mean bending strains

[TABLE]

As mentioned, these act as an easy to compute proxy for the (mean) principal curvatures. We observe how as $\theta$ increases both strains decrease almost by an equal amount as the body gradually opens up and flattens out, while retaining its radial symmetry. However, around $\theta\approx 86$ a stark change takes place and one of the principal strains decreases while the other increases. This reflects the abrupt change of the minimiser to a cylindrical shape. We observe the same phenomenon with the quotient of the principal axes of the deformed disk in the right plot of the same Figure.

The second initial condition tested is an orthotropically skewed paraboloid. Basically, a spherical cap is pressed from the sides to obtain a “potato chip”. Testing this shape will highlight the effect of the initial configuration on the final curvature. We examine its strains and symmetry in Figure

5.

Again there is a critical value of $\theta\approx 50$ around which the shape of the minimiser drastically changes. Note however how the change is now gradual and we see intermediate shapes.

Acknowledgements

This work was financially supported by project 285722765 of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), “Effektive Theorien und Energie minimierende Konfigurationen für heterogene Schichten”.