Finite element simulation of nonlinear bending models for thin elastic rods and plates

TL;DR

This paper develops and analyzes finite element methods for simulating nonlinear bending in thin elastic rods and plates, addressing complex constraints and applications in nanotechnology and biological structures.

Contribution

It introduces discretization techniques for nonlinear bending models based on recent mathematical frameworks, including isometry constraints and gradient flows, with rigorous analysis and practical experiments.

Findings

Effective numerical schemes for nonlinear bending problems

Successful simulation of elastic knots and singularities

Incorporation of torsion and self-avoidance in models

Abstract

Nonlinear bending phenomena of thin elastic structures arise in various modern and classical applications. Characterizing low energy states of elastic rods has been investigated by Bernoulli in 1738 and related models are used to determine configurations of DNA strands. The bending of a piece of paper has been described mathematically by Kirchhoff in 1850 and extensions of his model arise in nanotechnological applications such as the development of externally operated microtools. A rigorous mathematical framework that identifies these models as dimensionally reduced limits from three-dimensional hyperelasticity has only recently been established. It provides a solid basis for developing and analyzing numerical approximation schemes. The fourth order character of bending problems and a pointwise isometry constraint for large deformations require appropriate discretization techniques which are discussed in this article. Methods developed for the approximation of harmonic maps are adapted to discretize the isometry constraint and gradient flows are used to decrease the bending energy. For the case of elastic rods, torsion effects and a self-avoidance potential that guarantees injectivity of deformations are incorporated. The devised and rigorously analyzed numerical methods are illustrated by means of experiments related to the relaxation of elastic knots, the formation of singularities in a M\"obius strip, and the simulation of actuated bilayer plates.

Figures

Tables

Table 1. Table 1. Frequently used notation.

$(0,L)$, $\omega$, $\Omega$	one-, two-, and three-dimensional domains
$L^p(A;{ℝ}^{\ell })$	Lebesgue functions with values in ${ℝ}^{\ell }$
$W^{k,p}(A,{ℝ}^{\ell })$	Sobolev functions with values in ${ℝ}^{\ell }$
$H^k(A,{ℝ}^{\ell })$	Sobolev space with $p=2$
$x=(x_1,x_2,x_3)$	spatial variable
$x^{\prime }=(x_1,x_2)$	planar component of spatial variable
$|\cdot |$	Euclidean or Frobenius norm of a vector or matrix
$(\cdot ,\cdot )$, $\parallel \cdot \parallel$	scalar product and norm in $L^2$
$x\cdot y$, $A:B$	scalar products of vectors and matrices
$I_{\ell }$	identity matrix in ${ℝ}^{\ell \times \ell }$
$\mathrm{sym}\left(A\right)$, $\mathrm{tr}\left(A\right)$	symmetric part and trace of a matrix
$SO\left(3\right)$	orthogonal matrices with positive determinant
$y^{\prime }$	one-dimensional derivative
$\nabla y=[{\partial }_{1}y,{\partial }_{2}y,{\partial }_{3}y]$	gradient of a vector field
${\nabla }^{\prime }y$	planar component of gradient
$G^{\prime }$, $I^{\prime }$	total or Fréchet derivative
$P_k\left(A\right)$	polynomials of degree at most $k$ on a set $A$
$h$, $h_{\min }$	maximal and minimal mesh-sizes
${𝒯}_h$	triangulation with intervals or triangles
${𝒩}_h$, ${𝒮}_h$	nodes and sides in a triangulation
$z$, $z_T$, $z_S$	vertices and midpoints of elements, midpoints of sides
${𝒮}^{k,\ell }\left({𝒯}_h\right)$	elementwise degree $k$ polynomials in $C^{\ell }$
${ℐ}_h^{1,0}$, ${\widehat{ℐ}}_h^{1,0}$	global and elementwise nodal $P1$ interpolants
$Q_h$	elementwise averaging operator
$\parallel \cdot {\parallel }_{L_h^p}$, $\parallel \cdot {\parallel }_h$, ${(\cdot ,\cdot )}_h$	discrete $L^p$ norms, case $p=2$, discrete $L^2$ product
$\tau$	step size
$d_t{a}^k=\left(a^k-a^{k-1}\right)/\tau$	backward difference quotient for step size $\tau >0$
$\delta$	small thickness parameter
$I\left[y\right]$	energy functional
$L\left[y\right]$	linear operator
$W$	energy density
$Q_3$, $Q_{\mathrm{rod}}$, $Q_{\mathrm{plate}}$	quadratic forms
$\lambda$, $\mu$	Lamé parameters
$c_b$, $c_t$	bending and torsion rigidity
$𝒜$, ${𝒜}_h$	sets of admissible deformations
$ℱ\left[y\right]$, ${ℱ}_h\left[y_h\right]$	(shifted) tangent spaces
${(\cdot ,\cdot )}_{\star }$, ${(\cdot ,\cdot )}_{†}$	metrics used to define gradient flows
$c$, $c^{\prime }$, $c^{\prime \prime }$, …	generic constants

Table 2. Table 2. Average runtimes in seconds and average iteration numbers of the preconditioned conjugate gradient scheme (pcg) per iteration in parentheses for solution strategies (1)-(4) for the constrained linear systems of equations arising in the iterative approximation of harmonic maps on triangulations 𝒯 ℓ subscript 𝒯 ℓ \mathcal{T}_{\ell} of the unit cube with mesh-sizes h ℓ ≈ 2 − ℓ subscript ℎ ℓ superscript 2 ℓ h_{\ell}\approx 2^{-\ell} , ℓ = 2 , 3 , … , 6 ℓ 2 3 … 6 \ell=2,3,\dots,6 .

	saddle, direct	reduced, direct	reduced, pcg	reduced, pcg
$\#{𝒩}_{\ell }$	(backslash)	(backslash)	(diagonal $C^{𝖳}AC$)	(incompl. Chol.)
$125$	$\mathrm{0.000\hspace{0.17em}738\hspace{0.17em}8}$	$\mathrm{0.000\hspace{0.17em}324\hspace{0.17em}6}$	$\mathrm{0.000\hspace{0.17em}083\hspace{0.17em}5}$  (11.9)	$\mathrm{0.000\hspace{0.17em}030\hspace{0.17em}0}$  (6.0)
$729$	$\mathrm{0.015\hspace{0.17em}046\hspace{0.17em}6}$	$\mathrm{0.006\hspace{0.17em}483\hspace{0.17em}6}$	$\mathrm{0.000\hspace{0.17em}190\hspace{0.17em}6}$  (25.4)	$\mathrm{0.000\hspace{0.17em}074\hspace{0.17em}3}$ (10.2)
$4913$	$\mathrm{0.596\hspace{0.17em}117\hspace{0.17em}8}$	$\mathrm{0.243\hspace{0.17em}627\hspace{0.17em}3}$	$\mathrm{0.000\hspace{0.17em}920\hspace{0.17em}4}$   (52.5)	$\mathrm{0.000\hspace{0.17em}388\hspace{0.17em}4}$ (19.9)
$35937$	$\mathrm{34.472\hspace{0.17em}256}$	$\mathrm{8.263\hspace{0.17em}754\hspace{0.17em}2}$	$\mathrm{0.008\hspace{0.17em}106\hspace{0.17em}3}$ (105.9)	$\mathrm{0.003\hspace{0.17em}208\hspace{0.17em}2}$ (38.8)
$274625$	–	–	$\mathrm{0.109\hspace{0.17em}453\hspace{0.17em}9}$ (212.1)	$\mathrm{0.058\hspace{0.17em}556\hspace{0.17em}4}$ (76.1)

Equations

$$I_{3d}^{\delta }[y]={\int }_{{\Omega }_{\delta }}W(\nabla y)\mathrm{d}x$$

$$\underset{y\in \mathcal{A}}{\min }{I}_{3d}^{\delta }[y]=\mathcal{O}({\delta }^3)$$

$$y:(0,L)\to {\mathbb{R}}^3$$

$$|y^{\prime }(x_1)|=1$$

$$I_{\mathrm{rod}}[y]=\frac{1}{2}{\int }_0^L|y^{\prime \prime }(x_1){|}^2\mathrm{d}{x}_1.$$

$$(y^{\prime \prime },w^{\prime \prime })=(\lambda y^{\prime },w^{\prime })⟺y^{(4)}=(\lambda y^{\prime }{)}^{\prime },$$

$$0=z_0<z_1<\cdots <z_N=L$$

$$|y_h^{\prime }(z_i)|=1$$

$$y_h\mapsto I_{\mathrm{rod}}^h[y_h]=\frac{1}{2}{\int }_0^L|y_h^{\prime \prime }(x_1){|}^2\mathrm{d}{x}_1,$$

$$y_{h}\equiv\big{(}y_{h}(z_{i}),y_{h}^{\prime}(z_{i})\big{)}_{i=0,\dots,N}.$$

$$({\partial }_{t}y,w{)}_{\star }=-(y^{\prime \prime },w^{\prime \prime }{)}_{L^2},y(0)=y_0,$$

$${\partial }_t{y}^{\prime }(t,x_1)\cdot y^{\prime }(t,x_1)=0,w^{\prime }(x_1)\cdot y^{\prime }(t,x_1)=0.$$

$$|y^{\prime }(s,x_1){|}^2-1={\int }_0^s\frac{d}{dt}|y^{\prime }(t,x_1){|}^2\mathrm{d}t=2{\int }_0^s{\partial }_t{y}^{\prime }(t,x_1)\cdot y^{\prime }(t,x_1)\mathrm{d}t=0,$$

$$(d_t{y}_h^k,w_h{)}_{\star }=-([y_h^k{]}^{\prime \prime },[w_h{]}^{\prime \prime }{)}_{L^2},y_h^0=y_{0,h},$$

$$[d_t{y}_h^k{]}^{\prime }(x_i)\cdot [y_h^{k-1}{]}^{\prime }(x_i)=0,[w_h{]}^{\prime }(x_i)\cdot [y_h^{k-1}{]}^{\prime }(x_i)=0$$

$$I_{\mathrm{rod}}^h[y_h^k]+\tau \Vert d_t{y}_h^k{\Vert }_{\star }^2\le I_{\mathrm{rod}}^h[y_h^{k-1}].$$

$$\begin{array}{rl}|[y_h^k{]}^{\prime }(z_i){|}^2-1& =|[y_h^{k-1}](z_i{)}^{\prime }{|}^2+{\tau }^2|[d_t{y}_h^k{]}^{\prime }(z_i){|}^2-1\\ & =\cdots ={\tau }^2\sum _{\ell =1}^k|[d_t{y}_h^{\ell }{]}^{\prime }(z_i){|}^2.\end{array}$$

$$y:\omega \to {\mathbb{R}}^3$$

$$(\nabla y\big{)}^{\sf T}\nabla y=I_{2}$$

$$I_{\mathrm{plate}}[y]=\frac{1}{2}{\int }_{\omega }|D^{2}y{|}^2\mathrm{d}{x}^{\prime }.$$

$$(D^{2}y,D^{2}w)=0$$

$$(\nabla y{)}^{\mathsf{T}}\nabla w+(\nabla w{)}^{\mathsf{T}}\nabla y=0.$$

$$\big{(}\nabla y_{h}(z)\big{)}^{\sf T}\nabla y_{h}(z)=I_{2}.$$

$$I_{\mathrm{plate}}^h[y_h]=\frac{1}{2}{\int }_{\omega }|D_h^2{y}_h{|}^2\mathrm{d}{x}^{\prime }.$$

$$\big{(}y_{h}(z),\nabla y_{h}(z)\big{)}_{z\in\mathcal{N}_{h}}.$$

$$({\partial }_{t}y,w{)}_{\star }=-(D^{2}y,D^{2}w),y(0)=y_0,$$

$$L_{\nabla y}^{\mathrm{iso}}[{\partial }_t\nabla y]=0,L_{\nabla y}^{\mathrm{iso}}[\nabla w]=0,$$

$$L_A^{\mathrm{iso}}[B]=A^{\mathsf{T}}B+B^{\mathsf{T}}A.$$

$$(d_t{y}_h^k,w_h{)}_{\star }=-(D_h^2{y}_h^k,D_h^2{w}_h)$$

$$L_{\nabla y_h^{k-1}}^{\mathrm{iso}}[d_t\nabla y_h^k]=0,L_{\nabla y_h^{k-1}}^{\mathrm{iso}}[\nabla w_h]=0.$$

Full text

1. 1 Introduction
2. 2 Formal dimension reductions
3. 3 Convergent finite element discretizations
4. 4 Iterative solution via constrained gradient flows
5. 5 Linear finite element systems with nodal constraints
6. 6 Applications, modifications, and extensions
7. 7 Conclusions

1. Introduction

Thin elastic structures occur in various practical applications and in fact truly three-dimensional objects are hardly ever used. Important reasons for this are the reduction of weight and cost but also the special mechanical features of rods and plates. Correspondingly, their numerical treatment is expected to be more efficient when such structures can be described as lower-dimensional objects. Because of the different mechanical behavior they cannot be treated like three-dimensional objects and new discretization techniques are needed. Typical large bending deformations of rods and plates are distinct from those of three-dimensional objects and are illustrated in Figure 1.

In this article we address the numerical approximation of dimensionally reduced models for describing large deformations of thin elastic rods and plates. These models result from rigorous limiting processes of general three-dimensional hyperelastic material descriptions when the diameter of a circular rod or the thickness of a plate is small compared to the length or diameter and when the acting forces lead to deformations with energies comparable to the third power of the diameter or thickness. Examples of such situations are the bending of a springy wire or sheet of paper.

Characteristic for large nonlinear bending phenomena is that nearly no shearing or stretching effects of the object occur and that curvature quantities define the amount of energy required for particular deformations. These aspects become explicitly apparant in the dimensionally reduced models: the energy functionals depend on curvature quantities and an isometry condition arises in the vanishing thickness or diameter limit. In particular, this condition implies that length and angle relations remain unchanged by a deformation.

The employed models for elastic rods and plates result from dimension reductions of general descriptions for hyperelastic material behavior. We thus consider an energy density $W:\mathbb{R}^{3\times 3}\to\mathbb{R}$ and a corresponding energy minimization of

[TABLE]

in a set $\mathcal{A}\subset W^{1,p}(\Omega_{\delta};\mathbb{R}^{3})$ of admissible deformations $y:\Omega_{\delta}\to\mathbb{R}^{3}$ that includes boundary conditions. The parameter $\delta>0$ indicates a small diameter or thickness of the reference configuration $\Omega_{\delta}\subset\mathbb{R}^{3}$, e.g., $\Omega_{\delta}=(0,L)\times\delta S$ for a thin rod with cross-section $\delta S\subset\mathbb{R}^{2}$ containing zero or $\Omega_{\delta}=\omega\times(-\delta/2,\delta/2)$ for a thin plate with midplane $\omega\subset\mathbb{R}^{2}$. Assuming that the minimal energies are comparable to $\delta^{3}$, i.e.,

[TABLE]

as $\delta\to 0$, and following the contributions [38, 36, 48, 43], it is possible to identify limiting, dimensionally reduced theories that determine the corresponding limits of solutions as $\delta\to 0$. The particular cubic scaling characterizes bending phenomena of the elastic body and excludes membrane effects, we refer the reader to [37, 29, 34] for discussions of models corresponding to other scaling regimes. We outline the numerical methods that have been developed for simulating nonlinear bending behavior of rods and plates after a discussion of the related literature.

Throughout this article we use energy minimization principles to determine deformations subject to boundary conditions and external forces. For other approaches to the modeling of rods and plates via equilibria of forces or conservation of momentum we refer the reader to [2, 27, 3]. Only a few numerical methods have been discussed mathematically for the numerical solution of nonlinear bending models with inextensibility or isometry constraint. The articles [64, 21] devise various methods to compute discrete curvature quantities. The focus of this article is on the reliability of methods, i.e., the accuracy of finite element discretizations and the convergence of iterative solution methods for the discrete problems. The methods discussed here use techniques developed for the approximation of harmonic maps into surfaces in the articles [1, 6, 11]. We review the numerical treatment of rods following [9, 18] and plates as proposed in [8, 10]. For the efficient iterative solution we adopt ideas from [45, 41]. We discuss the treatment of bilayer plates following [14, 13], illustrate a method that enforces injectivity of deformations in the case of rods following [19, 17], and propose methods for the numerical solution of bending deformations with shearing effects following ideas from [12]. The problems considered in this article have similarities with problems related to the length-preserving elastic flow of curves and the surface area and volume preserving Willmore–Helfrich flow of closed surfaces but require different numerical methods. For contributions related to those problems we refer the reader to [33, 32, 4, 5, 52, 49, 24]; for examples of modern applications of nonlinear bending phenomena including the construction of micromachining fingers, the fabrication of nanotubes, the occurence of wrinkling in plastic sheets, and the description of certain properties of DNA molecules, we refer the reader to [58, 55, 57, 61].

1.1. Bending of elastic rods

We consider an elastic rod, e.g., a springy wire, which in its reference configuration occupies the region $(0,L)\times\{0\}\subset\mathbb{R}^{3}$. A low energy deformation

[TABLE]

leaves distances of pairs of points on the rod unchanged. This is described by the inextensibility (and incompressibility) condition

[TABLE]

for almost every $x_{1}\in(0,L)$. For appropriate boundary conditions a deformation then minimizes the bending energy

[TABLE]

The inextensibility condition implies that $y$ defines an arclength parametrization of the deformed rod and hence its curvature is given by the second derivative of $y$. The simple energy functional $I_{\rm rod}$, which has been proposed by Bernoulli in 1738, ignores torsion effects and arises as a special case of the dimension reduction from three-dimensional hyperelasticity.

It is interesting to see that the dimension reduction leads to significant changes in the nature of the energy functionals. The three-dimensional model depends on strains, is not constrained, and often provides existence of unique solutions. The dimensionally reduced functional $I_{\rm rod}$ depends on curvature, is constrained, and is singular in the sense that the set of admissible deformations may be empty, e.g., for extensive boundary conditions, and that solutions may be non-unique, e.g., for simple compressive boundary conditions. These aspects are related to the presence of a critical nonlinearity via a Lagrange multiplier for the inextensibility constraint in the Euler–Lagrange equations for critical points of $I_{\rm rod}$, i.e.,

[TABLE]

where the scalar function $\lambda$ depends nonlinearly on $y$. The explicit presence of a Lagrange multiplier can be avoided if only test functions are considered that satisfy the linearized inextensibility condition $y^{\prime}\cdot w^{\prime}=0$. This corresponds to normal, i.e., non-tangential perturbations of a curve in the energy minimization.

The inextensibility condition requires a suitable numerical treatment to avoid locking phenomena or other artifacts. For a partitioning of the interval $(0,L)$ with nodes

[TABLE]

and a subordinated conforming finite element space $\mathcal{A}_{h}\subset H^{2}(0,L;\mathbb{R}^{3})$ we impose the inextensibility condition only at these nodes, i.e.,

[TABLE]

for $i=0,1,\dots,N$. Since $y_{h}\in H^{2}(0,L;\mathbb{R}^{3})$ we obtain linear convergence with respect to the meshsize $h$ of the constraint violation error away from the nodes. The discrete minimization problem then seeks a minimizer $y_{h}\in\mathcal{A}_{h}$ for the functional

[TABLE]

subject to the nodal constraints $|y_{h}^{\prime}(z_{i})|=1$ for $i=0,1,\dots,N$. A possible choice of a finite element space uses piecewise cubic, continuously differentiable functions. This space has the advantage that its degrees of freedom are the positions and tangent vectors at the nodes, i.e.,

[TABLE]

The discretized inextensibility condition can thus be explicitly imposed on certain degrees of freedom, cf. Figure 2.

We iteratively solve the discrete minimization problem by using a gradient flow, i.e., on the continuous level we consider a family $y:[0,T]\times(0,L)\to\mathbb{R}^{3}$ of deformations that solve the evolution equation

[TABLE]

subject to the linearized inextensibility conditions

[TABLE]

Provided that we have $|y_{0}^{\prime}(x_{1})|^{2}=1$ it follows that

[TABLE]

i.e., the inextensibility condition is satisfied. We use an implicit discretization of the evolution equation and a semi-implicit treatment of the linearized constraint, i.e., with the backward difference quotient operator $d_{t}$ we consider the time-stepping scheme

[TABLE]

subject to the linearized constraints evaluated at the nodes, i.e.,

[TABLE]

for $i=0,1,\dots,N$. The scheme is unconditionally energy-decreasing and convergent to a stationary configuration, i.e., choosing the admissible test function $w_{h}=d_{t}y_{h}^{k}$ directly shows

[TABLE]

The inextensibility constraint will not be satisfied exactly at the nodes but its violation is controlled by the step size $\tau$ and the initial energy. A proof in the discrete setting imitates the continuous argument given above. Using the orthogonality $[d_{t}y_{h}^{k}]^{\prime}(z_{i})\cdot[y_{h}^{k-1}]^{\prime}(z_{i})=0$ and the property $|[y_{h}^{0}]^{\prime}(z_{i})|^{2}=1$ we have

[TABLE]

Because of the unconditional energy stability the term on the right-hand side is of order $\mathcal{O}(\tau)$. We will show below that these properties are also valid if torsion effects are taken into account.

1.2. Elastic plates

The mathematical description and numerical treatment of elastic plates generalizes that of elastic rods. In the dimensionally reduced model we consider deformations of a two-dimensional midplane

[TABLE]

that leave angle and area relations unchanged, i.e., they satisfy the isometry condition

[TABLE]

almost everywhere in $\omega\subset\mathbb{R}^{2}$ with the identity matrix $I_{2}\in\mathbb{R}^{2\times 2}$. This is equivalent to saying that the tangent vectors $\partial_{1}y$ and $\partial_{2}y$ of the deformed plate and the normal vector $b=\partial_{1}y\times\partial_{2}y$ define an orthonormal basis for $\mathbb{R}^{3}$ in almost every point $x^{\prime}\in\omega$. The actual deformation for appropriate boundary conditions minimizes the bending energy proposed by Kirchhoff in 1850,

[TABLE]

Because of the isometry condition, the integrand coincides with the mean curvature of the deformed plate while its Gaussian curvature vanishes. For a minimizing or critical isometry $y$ we have that

[TABLE]

for all test fields $w$ satisfying appropriate homogeneous boundary conditions and the linearized isometry condition

[TABLE]

A finite element discretization uses a possibly nonconforming finite element space such as so-called discrete Kirchhoff triangles and imposes the isometry condition in the set of nodes $\mathcal{N}_{h}$, i.e., for all $z\in\mathcal{N}_{h}$ we have

[TABLE]

With a discrete Hessian $D_{h}^{2}$ the numerical minimization is then realized for the functional

[TABLE]

In the case of the discrete Kirchhoff triangle, which may be seen as a natural generalization of the space of one-dimensional cubic $C^{1}$ functions, the degrees of freedom are the deformations and the deformation gradients in the nodes, i.e., the quantities

[TABLE]

An image of a discrete Kirchhoff deformation is depicted in Figure 3.

The isometry constraint is thus imposed directly on certain degrees of freedom. The iterative numerical minimization follows closely the approach used in the one-dimensional situation. For an initial isometry $y_{0}$, we consider the continuous evolution problem

[TABLE]

for appropriate test functions $w\in H^{2}(\omega;\mathbb{R}^{3})$ subject to the linearized isometry condition

[TABLE]

with the linearized isometry operator

[TABLE]

A semi-implicit discretization of this constrained evolution problem leads to a sequence of linearly constrained problems: given an admissible $y_{h}^{0}\in\mathcal{A}_{h}$ compute the sequence $(y_{h}^{k})_{k=0,1,\dots}$ via $y_{h}^{k}=y_{h}^{k-1}+\tau d_{t}y_{h}^{k}$, where $d_{t}y_{h}^{k}$ solves

[TABLE]

subject to the conditions

[TABLE]

Again, straightforward calculations show that the iteration is energy decreasing and convergent, and that the constraint violation is of order $\mathcal{O}(\tau)$.

1.3. Outline of the article

The article is organized as follows. In Section 2 we discuss the arguments that lead to dimensionally reduced models for elastic rods and plates in the case of small energies. We review partial $\Gamma$-convergence results and explain the occurrence of the inextensibility and isometry constraints. Section 3 is devoted to the convergent and practical finite element discretization of the one- and two-dimensional minimization problems describing the elastic deformation of rods and plates. The rigorous justification of the finite element methods will be established via showing $\Gamma$-convergence of the discretized functionals to the continuous one as discretization parameters tend to zero. Difficulties arise in the appropriate treatment of nonlinear constraints and higher order derivatives. The practical minimization of the discretized energy functionals is addressed in Section 4. We use appropriately discretized gradient flows that lead to sequences of linear systems of equations together with guaranteed energy decrease. Moreover, we verify that they converge to stationary configurations. In view of the nonuniqueness and limited additional regularity properties of solutions this appears to be the best attainable result if no further assumptions are made. The special saddle-point structure of the linear systems of equations that arise in the time steps is investigated in Section 5. It turns out that the nodal constraints can be incorporated in the solution space leading to reduced linear systems with symmetric and positive definite system matrices. Section 6 is concerned with extensions and modifications of the models and solution methods. In particular, we discuss the numerical treatment of bilayer bending problems, the inclusion of a self-avoidance potential, and a bending problem allowing for the formation of wrinkles. The final section provides a summary and conclusions of our considerations. We end the introduction with an overview of employed notation.

1.4. Notation

Throughout this article we use standard notation for derivatives and integrals, matrices and inner products, Lebesgue, Sobolev, and finite element spaces. The list in Table 1 provides an overview of the most important symbols.

2. Formal dimension reductions

Following the articles [38, 36, 43] we illustrate in this section how the dimensionally reduced minimization problems can be obtained from a general three-dimensional hyperelastic energy minimization problem:

[TABLE]

The set of admissible deformations $\mathcal{A}$ is assumed to be a weakly closed subset of a Sobolev space $W^{1,p}(\Omega;\mathbb{R}^{3})$ and required to include appropriate boundary conditions which imply a coercivity property. In an abstract way these are defined by a bounded linear operator

[TABLE]

and given data $\ell_{\rm bc}\in Y$ for a suitable linear space $Y$, e.g., traces of functions in $W^{1,p}(\Omega;\mathbb{R}^{3})$ restricted to a subset ${\Gamma_{\rm D}}$ of $\partial\Omega$. We assume that the energy density $W\in C^{2}(\mathbb{R}^{3\times 3})$ satisfies the following standard requirements:

• •

$W$ is frame-indifferent, i.e., for all $F\in\mathbb{R}^{3\times 3}$ and $Q\in SO(3)$ we have

[TABLE]

• •

$W$ vanishes at the identity $I_{3}\in\mathbb{R}^{3\times 3}$ and grows at least quadratically away from $SO(3)$, i.e., for all $F\in\mathbb{R}^{3\times 3}$ we have

[TABLE]

• •

$W$ is isotropic, i.e., for all $F\in\mathbb{R}^{3\times 3}$ and $R\in SO(3)$ we have

[TABLE]

From the first two conditions we have that $W^{\prime}(I_{3})=0$ and a Taylor expansion yields

[TABLE]

where $Q_{3}(G)=W^{\prime\prime}(I_{3})[G,G]$ is the quadratic form defined by the second variation of $W$ at the identity matrix. Incorporating the implicitly assumed homogeneity of the underlying material we have

[TABLE]

where the linear operator $\mathbb{C}:\mathbb{R}^{3\times 3}\to\mathbb{R}^{3\times 3}$ is given by

[TABLE]

with the Lamé parameters $\lambda,\mu>0$, the symmetric part $\operatorname{sym}(A)=(A+A^{\sf T})/2$, and the trace $\operatorname{tr}(A)=A:I_{3}$.

2.1. Elastic rods

We assume that the deformation $y:(0,L)\to\mathbb{R}^{3}$ of an elastic rod of vanishing thickness and length $L$ preserves distances, i.e., satisfies $|y^{\prime}(x_{1})|=1$, and complement the vector field $y^{\prime}$ by normal vector fields $b,d:(0,L)\to\mathbb{R}^{3}$ to an orthonormal frame

[TABLE]

We then consider a three-dimensional deformation $y_{\delta}$ obtained from extending the deformation of the centerline $(0,L)$ to the three-dimensional body $\Omega_{\delta}=(0,L)\times\delta S$ with scaled cross section $S\subset\mathbb{R}^{2}$ containing zero, i.e.,

[TABLE]

with a correction function $\beta:\Omega_{\delta}\to\mathbb{R}^{3}$. Inserting this deformation into the three-dimensional energy functional and considering the limit as $\delta\to 0$, we expect to identify a dimensionally reduced functional minimized by $y$ and the normal fields $b$ and $d$. We note that $x_{2},x_{3}=\mathcal{O}(\delta)$ and that we expect $\partial_{2}\beta,\partial_{3}\beta=\mathcal{O}(\delta^{-1})$. We therefore set

[TABLE]

with the matrix $R=[y^{\prime},b,d]\in SO(3)$. The matrix $R^{\sf T}\nabla y_{\delta}$ is thus a perturbation of the identity matrix $I_{3}$ and a Taylor expansion of the energy density yields with

[TABLE]

that we have

[TABLE]

Letting $\alpha=\delta^{2}R^{\sf T}\beta$ and noting that $R$ does not depend on $x_{2}$ and $x_{3}$ we have

[TABLE]

We note that since $R^{\sf T}R=I_{3}$ we have $(R^{\sf T})^{\prime}R=-R^{\sf T}R^{\prime}$ so that $R^{\sf T}R^{\prime}$ is skew-symmetric. A minimization of the integral of the energy density over the cross section $\delta S$ motivates defining for skew-symmetric matrices $A=[a_{1},a_{2},a_{3}]\in\mathbb{R}^{3\times 3}$ the reduced quadratic form $Q_{\rm rod}$ via

[TABLE]

With the particular representation of $Q_{3}$ by the Lamé parameters one finds with the entries $a_{ij}=-a_{ji}$ of $A$ for a circular cross section $S=B_{1/\pi}(0)$ that

[TABLE]

The constant factors on the right-hand side define the bending and torsion rigidities and are abbreviated by

[TABLE]

We always assume $\lambda,\mu>0$ so that $c_{\rm b}\geq 2c_{\rm t}$. For the particular matrix $A=R^{\sf T}R^{\prime}$ and $R=[y^{\prime},b,d]\in SO(3)$ we have that

[TABLE]

Noting that $y^{\prime\prime}\cdot y^{\prime}=0$ we have that

[TABLE]

is the squared curvature of the deformed rod and that

[TABLE]

is its squared torsion. We eliminate the variable $d$ via the identity $d=y^{\prime}\times b$ in what follows. We thus expect that the deformation $y:(0,L)\to\mathbb{R}^{3}$ of the centerline of a thin rod and the unit normal vector field $b:(0,L)\to\mathbb{R}^{3}$ solve the following dimensionally reduced problem:

[TABLE]

Here, we abbreviate

[TABLE]

The second part for justifying the dimensionally reduced model consists in showing that for any sequence $(y_{\delta})_{\delta>0}$ of three-dimensional deformations with $I_{\rm 3d}[y_{\delta}]\leq c\delta^{3}$ there exists an appropriate limit $(y,b)\in\mathcal{A}$ such that

[TABLE]

This so-called compactness property is proved in [43] which provides the complete rigorous dimension reduction in a more general setting. We refer the reader to [42] for further aspects of the description of elastic rods.

Remark 2.1**.**

To illustrate that the inextensibility or isometry condition $|y^{\prime}|=1$ arises naturally in the dimension reduction we consider the planar deformation of a two-dimensional thin beam $\Omega=(0,L)\times(-\delta/2,\delta/2)$ with the simple energy density

[TABLE]

We assume that the deformation is given by

[TABLE]

for a deformation $y:(0,L)\to\mathbb{R}^{2}$ of the centerline and a corresponding normal field $b:(0,L)\to\mathbb{R}^{2}$, i.e., we have $y^{\prime}(x_{1})\cdot b(x_{1})=0$. Noting that

[TABLE]

we find that

[TABLE]

We insert this expression into the energy functional and carry out the integration in vertical direction, i.e.

[TABLE]

For a cubic scaling of the elastic energy we need $A=0$, i.e., $|y^{\prime}|^{2}=1$ and $|b|^{2}=1$. This implies $b^{\prime}\cdot b=0$, hence $|B|^{2}=4|y^{\prime}\cdot b^{\prime}|^{2}$, and shows that up to terms of order $\delta^{5}$ we have

[TABLE]

where we used the identities $y^{\prime}\cdot b=-y^{\prime\prime}\cdot b$ and $y^{\prime\prime}\cdot y^{\prime}=0$ in combination with the fact that $(y^{\prime},b)$ is an orthonormal basis in $\mathbb{R}^{2}$ so that $(y^{\prime}\cdot b)^{\prime}=0$.

2.2. Elastic plates

To explain the derivation of the bending model for elastic plates we consider an isometry $y:\omega\to\mathbb{R}^{3}$ and a corresponding unit normal field $b:\omega\to\mathbb{R}^{3}$, i.e., we have

[TABLE]

for $1\leq i,j\leq 2$ and

[TABLE]

for almost every $x^{\prime}\in\omega$ and $j=1,2$. We define a deformation $y_{\delta}$ of the three-dimensional body $\Omega_{\delta}=\omega\times(-\delta/2,\delta/2)$ by extending $y$ in normal direction, i.e.,

[TABLE]

with a quadratic correction term $\beta:\omega\to\mathbb{R}^{3}$. Using the planar gradient $\nabla^{\prime}=[\partial_{1},\partial_{2}]$ we have that

[TABLE]

and

[TABLE]

This implies that

[TABLE]

We insert the deformation $y_{\delta}$ into the hyperelastic energy functional, use the Taylor expansion $W(I+x_{3}G)=Q_{3}(x_{3}G)+o(x_{3}^{2}|G|^{2})$, with $G=R^{\sf T}\nabla y_{\delta}$, and carry out the integration in vertical direction. This leads to

[TABLE]

The correction field $\beta:\omega\to\mathbb{R}^{3}$ is eliminated via a pointwise minimization, i.e., for $M\in\mathbb{R}^{2\times 2}$ extended by a vanishing third row to a matrix $\widehat{M}\in\mathbb{R}^{3\times 2}$, we define

[TABLE]

Since we assume a homogeneous and isotropic material one obtains for a symmetric matrix $M\in\mathbb{R}^{2\times 2}$ that

[TABLE]

For the matrix $\widehat{M}=R^{\sf T}\nabla^{\prime}b\in\mathbb{R}^{3\times 2}$ the third row vanishes and its uppper $2\times 2$ submatrix coincides with the second fundamental form $II$ of the surface parametrized by $y$, i.e.,

[TABLE]

for $i,j=1,2$, where in fact $b=\pm\partial_{1}y\times\partial_{2}y$. Using that $y$ is an isometry we have that the squared mean curvature is up to a fixed factor given by the identical expressions

[TABLE]

Hence, the dimensionally reduced problem seeks an isometric deformation that minimizes the integral of the squared Hessian:

[TABLE]

The bending rigidity is defined by $c_{\rm b}=2\mu+\lambda 2\mu/(2\mu+\lambda)$. As in the case of rods, a rigorous derivation additionally requires showing that the functional defines a general lower bound, i.e., establishing a lim-inf inequality, and we refer the reader to [38, 36] for details. Analogously to the functional, also the boundary conditions change their nature in the dimension reduction. A fixed part of the lateral boundary leads to a clamped boundary condition in the reduced model which imposes a condition on the deformation and its gradient.

3. Convergent finite element discretizations

We discuss in this section the discretization of the dimensionally reduced nonlinear bending models using appropriate finite element methods. Challenges are the treatment of higher order derivatives and a nonlinear pointwise constraint. We establish the correctness of the discretizations by showing that the discrete functionals $I^{h}$ converge in the sense of $\Gamma$-convergence with respect to weak convergence on a space $X$, cf., e.g., [31], to the continuous, dimensionally reduced functional $I$. We use the terminology almost-minimizing for a sequence of objects that are minimizers of a sequence of functions up to tolerances that converge to zero with $h$. This follows from verifying the following three conditions:

• (a)

Well posedness or equicoercivity: The discrete functionals are uniformly coercive, i.e., if $I^{h}[y_{h}]\leq c$ then it follows that $\|y_{h}\|_{X}\leq c^{\prime}$ with $h$-independent constants $c,c^{\prime}\geq 0$, and admit discrete minimizers.

• (b)

Stability or lim-inf inequality: If $(y_{h})_{h>0}\subset X$ is a bounded sequence of discrete almost-minimizers then every weak accumulation point $y$ belongs to the set of admissible deformations $\mathcal{A}$ and we have

[TABLE]

• (c)

Consistency or lim-sup inequality: For every $y\in\mathcal{A}$ there exists a sequence $(y_{h})_{h>0}$ of admissible discrete deformations such that $y_{h}\rightharpoonup y$ in $X$ and

[TABLE]

It is an immediate consequence of (a)–(c) that sequences of discrete almost-minimizers accumulate at minimizers of the continuous problem. Well posedness typically follows from coercivity properties of the functional $I$ while the stability is established with the help of lower semicontinuity properties of $I$. If the union of discrete sets of admissible deformations is dense in the set of admissible deformations then consistency is obtained via continuity properties of $I$. We specify these concepts for the finite element approximation of elastic deformations of rods and plates in what follows. We always use a regular triangulation $\mathcal{T}_{h}$ of the domain $A=(0,L)$ or $A=\omega$ into intervals or triangles, respectively, with a set of nodes (vertices of elements) denoted $\mathcal{N}_{h}$, i.e.,

[TABLE]

We often use numerical integration or quadrature, defined with the elementwise applied nodal interpolation via

[TABLE]

for elementwise continuous functions $v,w:A\to\mathbb{R}^{\ell}$ with $A\subset\mathbb{R}^{d}$ and

[TABLE]

If $p=2$ we write $\|v\|_{h}$ instead of $\|v\|_{L^{2}_{h}}$. We note that these expressions define equivalent scalar products and norms on function spaces containing elementwise polynomials of bounded degree, cf., e.g., [10].

3.1. Elastic rods

For a discretization of the bending-torsion model (2.1) for elastic rods we first derive a suitable reformulation of the minimization problem. Recalling that for an admissible pair $(y,b)\in\mathcal{A}$ and the vector $d=y^{\prime}\times b$ we have that $[y^{\prime},b,d]\in SO(3)$ almost everywhere in the interval $(0,L)$, we deduce that

[TABLE]

Since $|b|^{2}=1$ the last term on the right-hand side vanishes while the orthogonality $b\cdot y^{\prime}=0$ implies that $b^{\prime}\cdot y^{\prime}=-b\cdot y^{\prime\prime}$. We thus have that

[TABLE]

This identity leads to the equivalent representation

[TABLE]

The dimension reduction of Section 2.1 shows that we have $c_{\rm b}\geq 2c_{\rm t}$ so that the last term is controlled by the first one and the coercivity of $I_{\rm rod}^{h}$ becomes explicit. Another advantage of this representation is that the last term is separately concave which allows for an effective iterative treatment. To define the discrete funtional $I_{\rm rod}^{h}$ we consider a partitioning of the reference interval $(0,L)$ defined by sets of nodes $\mathcal{N}_{h}$ and elements $\mathcal{T}_{h}$. For this partitioning we define the linear and cubic finite element spaces with different differentiability requirements via

[TABLE]

with sets of polynomials of degree at most $k$ on $T$ given by $P_{k}(T)$. The degrees of freedom of the finite element spaces are are depicted in Figure 4.

It is straightforward to verify that there exist nodal bases $(\varphi_{z})_{z\in\mathcal{N}_{h}}$ and $(\psi_{z,j})_{z\in\mathcal{N}_{h},j=1,2}$ such that for all $\phi_{h}\in\mathcal{S}^{1,0}(\mathcal{T}_{h})$ and $v_{h}\in\mathcal{S}^{3,1}(\mathcal{T}_{h})$ we have

[TABLE]

The right-hand sides define nodal interpolation operators $\mathcal{I}_{h}^{1,0}$ and $\mathcal{I}_{h}^{3,1}$ on $C^{0}([0,L])$ and $C^{1}([0,L])$, respectively. For ease of notation we use the abbreviation

[TABLE]

For efficient numerical quadrature we introduce the elementwise averaging operator $Q_{h}$ defined for a vector field $v\in L^{1}(0,L;\mathbb{R}^{3})$ and every element $T\in\mathcal{T}_{h}$ via

[TABLE]

With the product finite element space $V_{\rm rod}^{h}$ and the operator $Q_{h}$ we consider the following discretization of the minimization problem (2.1) in which the pointwise orthogonality relation $y^{\prime}\cdot b=0$ is approximated via a penalty term:

[TABLE]

Note that the constraints are imposed on particular degrees of freedom which makes the method practical. We have the following existence and convergence result.

Proposition 3.1** (Convergent approximation).**

For every pair $(h,\varepsilon)>0$ there exists a minimizer $(y_{h},b_{h})\in\mathcal{A}_{h}$ for $I_{\rm rod}^{h,\varepsilon}$ satisfying

[TABLE]

As $(h,\varepsilon)\to 0$ we have that every accumulation point of a sequence of discrete almost-minimizers is a minimizer for $I_{\rm rod}$ in $\mathcal{A}$.

Proof (sketched).

We outline the main arguments of the proof and refer the reader to [9, 18] for details.

(a) Let $(y_{h},b_{h})\in\mathcal{A}_{h}$ with $I_{\rm rod}^{h,\varepsilon}[y_{h},b_{h}]\leq c$. The coercivity of the discretized functional follows from the fact that $c_{\rm b}\geq 2c_{\rm t}$ and the identity

[TABLE]

with the positive semi-definite matrix

[TABLE]

The discrete coercivity and the continuity properties of the functional $I_{\rm rod}^{h,\varepsilon}$ imply the existence of discrete minimizers.

(b) Given a sequence of discrete almost-minimizers $(y_{h},b_{h})_{h,\varepsilon>0}$ one first checks that accumulation points $(y,b)$ as $(h,\varepsilon)\to 0$ belong to the continuous admissible set $\mathcal{A}$. Noting that $b_{h}\to b$ strongly in $L^{\infty}$ we find that

[TABLE]

(c) It remains to show that $I_{\rm rod}[y,b]$ is minimal. For this, we choose a smooth almost-minimizing pair $(\widetilde{y},\widetilde{b})\in\mathcal{A}$ obtained from an appropriate regularization of a minimizing pair and verify that the sequence of interpolants $(\widetilde{y}_{h},\widetilde{b}_{h})$ satisfies $\lim_{(h,\varepsilon)\to 0}I_{\rm rod}^{h,\varepsilon}[y_{h},b_{h}]=I_{\rm rod}[\widetilde{y},\widetilde{b}]$. ∎

Remark 3.2**.**

We note that the result of the proposition can be also be established if the orthogonality relation $y_{h}^{\prime}\cdot b_{h}=0$ is imposed exactly in the nodes of the triangulation. For an efficient numerical solution of the minimization problem the approximation via a separately convex term is advantageous as this allows for a decoupled treatment of the variables.

3.2. Elastic plates

Constructing finite element spaces that provide convergent second order derivatives and which are efficiently implementable is significantly more challenging in two space dimensions. Among the various possibilities is the discrete Kirchhoff triangle, cf., e.g., [25], which defines a nonconforming finite element method in the sense that its elements do not belong to $H^{2}$. The space can be seen as a natural generalization of the space of one-dimensional cubic $C^{1}$ functions since the degrees of freedom are the deformations and deformation gradients at the nodes of a triangulation which are appropriately interpolated on the individual elements. To define this finite element space we choose a triangulation $\mathcal{T}_{h}$ of $\omega$ into triangles and set

[TABLE]

Here, $P_{3-}$ denotes the subset of cubic polynomials on $T$ obtained by eliminating the degree of freedom associated with the midpoint $z_{T}$ of $T$, i.e., we have

[TABLE]

The degrees of freedom in $\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})$ are the function values and the derivatives at the vertices of the elements. It is interesting to note that a particular basis for $\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})$ will not be needed. A canonical interpolation operator $\mathcal{I}_{h}^{\rm dkt}:C^{1}(\overline{\omega})\to\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})$ is defined by requiring that the identities

[TABLE]

hold at all nodes $z\in\mathcal{N}_{h}$. The employed finite element spaces are depicted in Figure 5.

Crucial for the finite element discretization of the bending problem is the definition of a discrete gradient operator

[TABLE]

which allows us to define discrete second order derivatives of functions $w_{h}\in\mathcal{S}^{\rm dkt}$ via

[TABLE]

Here we make use of the fact that $\nabla_{h}w_{h}\in H^{1}(\Omega;\mathbb{R}^{2})$. The discrete gradient operator $\nabla_{h}$ is for given $w_{h}\in\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})$ defined as the unique piecewise quadratic, continuous vector field $q_{h}\in\mathcal{S}^{2,0}(\mathcal{T}_{h})^{2}$ that satisfies the condition

[TABLE]

for all $z\in\mathcal{N}_{h}$ while the degrees of freedom associated with the sides of elements are defined by the two conditions

[TABLE]

for all sides $S=[z_{S}^{1},z_{S}^{2}]\in\mathcal{S}_{h}$ with normals $n_{S}$, tangent vectors $t_{S}$, and midpoints $z_{S}=(z_{S}^{1}+z_{S}^{2})/2$. For $w\in C^{1}(\overline{\omega})$, we set

[TABLE]

With the discrete second derivatives we are in a position to state the finite element discretization of the plate bending model:

[TABLE]

To prove the correctness of this discretization we show that existing finite element minimizers accumulate at admissible isometries, incorporate that the bending energy is weakly lower semicontinuous, and use that isometries can be approximated by smooth isometries which is a result proved in [47, 40].

Theorem 3.3** (Convergent approximation).**

For every $h>0$ there exists a minimizer $y_{h}\in\mathcal{A}_{h}$. If $(y_{h})_{h>0}$ is a sequence of almost-minimizers, then $\|\nabla y_{h}\|\leq c$, for all $h>0$, and every accumulation point $y\in H^{1}(\omega;\mathbb{R}^{3})$ of the sequence is a strong accumulation point, belongs to the continuous admissible set $\mathcal{A}$ and is minimal for $I_{\rm plate}$.

Proof (sketched).

We follow the typical steps for establishing a $\Gamma$-convergence result provided in [7, 10].

(a) Using the boundary conditions included in the set $\mathcal{A}_{h}$ it follows that the mapping $z_{h}\mapsto\|D_{h}^{2}z_{h}\|$ is a norm on the subset of $\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})^{3}$ with functions satisfying corresponding homogeneous boundary conditions. This leads to a coercivity property and the existence of discrete solutions.

(b) The uniform discrete coercivity property implies that the sequences $(D_{h}^{2}y_{h})_{h>0}$ and $(\nabla y_{h})_{h>0}$ have weak accumulation points $\xi$ and $\nabla y$ in $L^{2}$ which are compatible in the sense that $\xi=D^{2}y$. Moreover, we have that $y\in\mathcal{A}$ and weak lower semicontinuity of the $L^{2}$ norm shows that

[TABLE]

(c) Let $y\in\mathcal{A}$ be a minimizer for $I_{\rm plate}$. The continuity of the functional $I_{\rm plate}$ with respect to the strong topology in $H^{2}$ in combination with the density results for smooth isometries established in [40] allow us to assume that $y$ is smooth. We may thus define an approximating sequence of finite element functions by setting $\widetilde{y}_{h}=\mathcal{I}_{h}^{\rm dkt}[y]$. Approximation properties of the interpolation operator lead to the inequality

[TABLE]

which proves the statement. ∎

4. Iterative solution via constrained gradient flows

The practical solution of the finite element discretizations of the nonlinear bending problems is nontrivial due to the presence of nonlinear pointwise constraints and the corresponding lack of higher regularity properties. To provide a reliable strategy that decreases the energy we adopt gradient flow strategies. Our estimates show that these converge to stationary, low energy configurations. We will always use a linearized treatment of the constraints which is then discretized semi-implicitly. This makes the iterative scheme practical. To illustrate the main idea, consider the following abstract minimization problem in a Hilbert space $X\subset L^{2}(\Omega;\mathbb{R}^{\ell})$:

[TABLE]

Here, we assume that the constraint is understood pointwise with a function $G:\mathbb{R}^{\ell}\to\mathbb{R}$. The Euler–Lagrange equations for the problem are then formally given by the identity

[TABLE]

for all $w\in X$ with a Lagrange multiplier $\lambda\in L^{1}(\Omega)$. Note that the term involving $\lambda$ disappears if $w$ satisfies $G^{\prime}[y;w]=0$ and that this is sufficient to characterize a stationary point subject to the constraint. The corresponding gradient flow is formally defined via

[TABLE]

Our corresponding time-stepping scheme uses the backward difference quotient operator $d_{t}a^{k}=(a^{k}-a^{k-1})/\tau$ for a step-size $\tau>0$ and determines iterates via the linearly constrained problems

[TABLE]

We specify the meaning of the iterative scheme in the following algorithm.

Algorithm 4.1** (Abstract constrained gradient descent).**

*Let $y^{0}\in X$ be such that $G[y^{0}]=0$ and $I[y^{0}]<\infty$ and choose $\tau>0$, set $k=1$.

(1) Compute $y^{k}\in X$ such that*

[TABLE]

for all $w\in X$ under the constraints $d_{t}y^{k},\,w\in\ker G^{\prime}[y^{k-1}]$, i.e.,

[TABLE]

(2) Stop the iteration if $\|d_{t}y^{k}\|_{X}\leq\varepsilon_{\rm stop}$; otherwise increase $k\to k+1$ and continue with (1).

We remark that it is useful to regard $d_{t}y^{k}$ rather than $y^{k}$ as the unknown in the iteration steps. In particular, we may eliminate $y^{k}$ via the identity $y^{k}=y^{k-1}+\tau d_{t}y^{k}$ with the known function $y^{k-1}$. A geometric interpretation of the iteration is that given an iterate $y^{k-1}$ the correction $d_{t}y^{k}$ is computed in the tangent space of the level set $\mathcal{M}^{k-1}$ of $G$ defined by the value $G[y^{k-1}]$. Note that we do not use a projection step onto the zero level set of $G$ since this will in general not be energy stable. Figure 6 illustrates the conceptual idea of Algorithm 4.1.

The following theorem states the main features of Algorithm 4.1, i.e., its unconditional energy stability with the resulting convergence to a stationary configuration of lower energy and the control of the constraint violation by the step size.

Theorem 4.2** (Convergent iteration).**

Assume that $I$ is convex, coercive, continuous, and Fréchet differentiable on $X$ and $G:\mathbb{R}^{\ell}\to\mathbb{R}$ is twice differentiable with uniformly bounded second derivative, i.e., we have $G^{\prime\prime}[r;s,s]\leq 2c_{G}|s|^{2}$ for all $r,s\in\mathbb{R}^{\ell}$. Then, the iterates of Algorithm 4.1 are uniquely defined and satisfy for $L=0,1,2,\dots$ the energy estimate

[TABLE]

Moreover, if $\|\cdot\|_{H}$ is a norm on $X$ with the property $\||z|^{2}\|_{H}\leq c_{H}\|z\|_{X}^{2}$ for all $z\in X$ then we have the constraint violation bound

[TABLE]

where $e_{0}=I[y^{0}]$. In particular, we have that $\|d_{t}y^{k}\|_{X}\to 0$ and Algorithm 4.1 terminates within a finite number of iterations. The output $y^{L}\in X$ satisfies the residual estimate

[TABLE]

Proof.

(a) The existence of the iterates follows by applying the direct method of the calculus of variations to the minimization problems:

[TABLE]

The solution is unique and the corresponding Euler–Lagrange equation coincides with the equation that defines the iterates in Algorithm 4.1.

(b) Choosing the admissible test function $w=d_{t}y^{k}$ and using the convexity of $I$ leads to

[TABLE]

A multiplication by $\tau$ and summation over $k=1,2,\dots,L$ prove the energy stability.

(c) With a Taylor expansion of $G$ about the iterate $y^{k-1}$ and the imposed identity $\delta G[y^{k-1};d_{t}y^{k}]=0$ we find that

[TABLE]

Repeating this argument and noting that $G[y^{0}]=0$ leads to the estimate

[TABLE]

Applying the norm $\|\cdot\|_{H}$ to the estimate, using the triangle inequality, and incorporating the assumed bound $\||z|^{2}\|_{H}\leq c_{H}\|z\|_{X}^{2}$ as well as the energy bound proves the estimate for the constraint violation error.

(d) The estimate for $I^{\prime}[y^{L};w]$ is an immediate consequence of the bound $\|d_{t}y^{L}\|_{X}\leq\varepsilon_{\rm stop}$. ∎

Examples of pairs of norms $\|\cdot\|_{H}$ and $\|\cdot\|_{X}$ that satisfy the assumed estimate are the $L^{1}$ norm in combination with the $L^{2}$ norm or the $L^{\infty}$ norm together with a Sobolev norm in $H^{s}$ with $s$ sufficiently large. It is remarkable that the violation of the constraint is independent of the number of iterations. An explanation for this is that the updates $d_{t}y^{k}$ converge quickly to zero in the gradient flow iteration.

4.1. Elastic rods

We apply the abstract framework for constrained minimization problems to the energy functional describing the bending-torsion behavior of elastic rods. For a vector field $y_{h}\in\mathcal{S}^{3,1}(\mathcal{T}_{h})^{3}$ we set

[TABLE]

while for a vector field $b_{h}\in\mathcal{S}^{1,0}(\mathcal{T}_{h})^{3}$ we define

[TABLE]

The functionals $L_{{\rm bc},y}^{\rm rod}$ and $L_{{\rm bc},b}^{\rm rod}$ are the components of $L_{{\rm bc}}^{\rm rod}$ corresponding to the variables $y$ and $b$, respectively. We generate a sequence $(y_{h}^{k},b_{h}^{k})_{k=0,1,\dots}$ that approximates a stationary configuration for $I_{\rm rod}^{h,\varepsilon}$ with the following algorithm.

Algorithm 4.3** (Gradient descent for elastic rods).**

*Choose an initial pair $(y_{h}^{0},b_{h}^{0})\in\mathcal{A}_{h}$ and a step size $\tau>0$, set $k=1$.

(1) Compute $d_{t}y_{h}^{k}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ such that for all $w_{h}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ we have*

[TABLE]

(2) Compute $d_{t}b_{h}^{k}\in\mathcal{E}_{h}[b_{h}^{k-1}]$ such that for all $r_{h}\in\mathcal{E}_{h}[b_{h}^{k-1}]$ we have

[TABLE]

(3) Stop the iteration if

[TABLE]

otherwise, increase $k\to k+1$ and continue with (1).

Again, it is useful to view $d_{t}y_{h}^{k}$ and $d_{t}b_{h}^{k}$ as the unknowns in Steps (1) and (2) instead of $y_{h}^{k}=y_{h}^{k-1}+\tau d_{t}y_{h}^{k}$ and $b_{h}^{k}=b_{h}^{k-1}+\tau d_{t}b_{h}^{k}$. The algorithm exploits the fact that the penalty term is separately convex while the nonquadratic contribution to the torsion term is separately concave. Therefore, the decoupled semi-implicit treatment of these terms is natural and unconditionally energy stable.

Proposition 4.4** (Convergent iteration).**

Assume that we have

[TABLE]

for all $(w_{h},r_{h})\in V^{h}_{\rm rod}$ with $L_{\rm bc}^{\rm rod}[w_{h},r_{h}]=0$. Algorithm 4.3 is well defined and produces a sequence $(y_{h}^{k},b_{h}^{k})_{k=0,1,\dots}$ such that for all $L\geq 0$ we have

[TABLE]

The iteration controls the unit-length violation via

[TABLE]

where $e_{0,h}=I_{\rm rod}^{h,\varepsilon}[y_{h}^{0},b_{h}^{0}]$. In particular, the algorithm terminates within a finite number of iterations.

Proof.

(a) To prove the stability estimate we note that the functional

[TABLE]

is separately convex, i.e., convex in $y_{h}$ and in $b_{h}$. Therefore, we have that

[TABLE]

which by summation leads to the inequality

[TABLE]

Similarly, the functional

[TABLE]

is separately convex and we have

[TABLE]

By choosing $w_{h}=d_{t}y_{h}^{k}$ and $r_{h}=d_{t}b_{h}^{k}$ in the equations of Steps (2) and (3) of Algorithm 4.3 we thus find that

[TABLE]

Since

[TABLE]

we deduce the asserted estimate.

(b) The nodal orthogonality conditions encoded in the spaces $\mathcal{F}_{h}[y_{h}^{k-1}]$ and $\mathcal{E}_{h}[b_{h}^{k-1}]$ lead to the relations

[TABLE]

By induction and incorporation of the stability estimate we deduce the asserted estimates for the constraint violation. ∎

We illustrate the performance of Algorithm 4.3 via a numerical experiment showing the relaxation of a twisted initially flat curve which is clamped at both ends. We plotted in the bottom of Figure 7 the total energy and the torsion contribution defined by

[TABLE]

We observe that the curve quickly releases its large energy and becomes a spatial curve attaining a stationary configuration with equilibrated curvature after approximately 2000 iterations. The model parameters used in the simulation were set to $c_{\rm b}=2$ and $c_{\rm t}=1$. We used a partition into 1006 subintervals corresponding to a mesh-size $h=1/80$. The step-size $\tau$ and the penalty parameter $\varepsilon$ were chosen proportional to $h$.

4.2. Elastic plates

We next apply the conceptual approach for solving nonlinearly constrained minimization problems outlined above to the case of approximating bending isometries. For this we recall that the discrete bending energy is given by

[TABLE]

with the set of admissible discrete deformations defined as

[TABLE]

Using the linearized isometry operator

[TABLE]

we define the (shifted) tangent space of the set of discrete isometric deformations $\mathcal{A}_{h}$ at a deformation $y_{h}$ via

[TABLE]

We then decrease the bending energy for given boundary conditions by iterating the steps of the following algorithm.

Algorithm 4.5** (Gradient descent for elastic plates).**

*Choose an initial $y_{h}^{0}\in\mathcal{A}_{h}$ and a step size $\tau>0$, set $k=1$.

(1) Compute $d_{t}y_{h}^{k}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ such that for all $w_{h}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ we have*

[TABLE]

(2) Stop if $\|d_{t}y_{h}^{k}\|_{\star}\leq\varepsilon_{\rm stop}$; otherwise, increase $k\to k+1$ and continue with (1).

The iterates $(y_{h}^{k})_{k=0,1,\dots}$ will in general not satisfy the nodal isometry constraint exactly, but the violation is again independent of the number of iterations and controlled by the step size $\tau$.

Proposition 4.6** (Convergent iteration).**

The iterates $(y_{h}^{k})_{k=0,1,\dots}$ of Algorithm 4.5 are well defined and satisfy for every $L\geq 0$ the energy estimate

[TABLE]

Moreover, if $\|\nabla w_{h}\|_{h}\leq c_{\star}\|w_{h}\|_{\star}$ for all $w_{h}\in\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})^{3}$ with $L_{\rm bc}^{\rm plate}[w_{h}]=0$, then we have the constraint violation bound

[TABLE]

where $e_{0,h}=I_{\rm plate}^{h}[y_{h}^{0}]$.

Proof (sketched).

(a) Since for any $y_{h}^{k-1}\in\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})^{3}$ we have that $\mathcal{F}_{h}[y_{h}^{k-1}]$ is a nonempty linear space the Lax-Milgram lemma implies that there exists a unique solution $y_{h}^{k}\in\mathcal{F}_{h}[y_{h}^{k-1}]$.

(b) The energy decay property is an immediate consequence of choosing $w_{h}=d_{t}y_{h}^{k}$ and using the binomial formula

[TABLE]

(c) The error bound for the isometry violation follows from the orthogonality defined by the linearized isometry condition and the energy decay property as in the proof of Proposition 4.2. ∎

Figure 8 illustrates the discrete evolution defined by Algorithm 4.5 via snapshots of different iterates. The clamped boundary conditions imposed at the ends $\gamma_{\rm D}=\{0,L\}\times[0,w]$ of the strip $\omega=(0,L)\times(0,w)$ with length $L=10$ and width $w=1$ are defined via the operator

[TABLE]

and functions $y_{\rm D}\in L^{2}({\gamma_{\rm D}};\mathbb{R}^{3})$ and $\widehat{y}_{\rm D}\in L^{2}({\gamma_{\rm D}};\mathbb{R}^{3\times 2})$. These functions are constructed in such a way that the segment $\{L\}\times[0,w]$ is rotated and mapped onto the fixed opposite segment $\{0\}\times[0,w]$. In this way the formation of a Möbius strip is enforced, we included a small linear forcing term to avoid certain nonuniqueness effects. We observe that the nonsmooth choice of the starting value with large bending energy does not influence the robustness of the iteration and that within less than $10.000$ iterations a stationary configuration for the stopping parameter $\varepsilon_{\rm stop}=5\cdot 10^{-3}$ and evolution metric with $\|\cdot\|_{\star}=\|D_{h}^{2}\cdot\|$ is attained. We also observe that we obtain a satisfactory stationary shape for coarse triangulations and that the energy decreases monotonically as predicted, with a small violation of the isometry constraint indicated by the quantity

[TABLE]

which appears to be nearly independent of the iteration. We finally remark that we observe concentrations of curvature at boundary points corresponding to certain singularities discussed in [15].

5. Linear finite element systems with nodal constraints

The discrete gradient flows devised in the previous sections lead to linear systems of equations in the time steps that have a special saddle point structure. In particular, they are given in the form

[TABLE]

with a fixed positive definite symmetric matrix $A\in\mathbb{R}^{n\times n}$ and a matrix $B\in\mathbb{R}^{p\times n}$ that changes in the iteration steps. The matrix $B$ is of full rank and block diagonal, i.e.,

[TABLE]

with vectors $b_{i}\in\mathbb{R}^{\ell}\setminus\{0\}$ for $i=1,2,\dots,p$ and $n=p\ell$. We note that the solution $x\in\mathbb{R}^{n}$ of the linear system of equations above is equivalently characterized by the system

[TABLE]

for all $z\in\mathbb{R}^{n}$ belonging to the kernel of $B$, i.e., subject to the conditions

[TABLE]

To obtain a simpler, unconstrained system of equations we choose for each $i=1,2,\dots,p$ a set of orthonormal vectors $(c_{2}^{i},c_{3}^{i},\dots,c_{\ell}^{i})\subset\mathbb{R}^{\ell}$ that are orthogonal to $b_{i}$, i.e., we have

[TABLE]

We may then represent the kernel of $B$ by the image of the matrix $C\in\mathbb{R}^{p\ell\times p(\ell-1)}$ defined by

[TABLE]

The matrix $C$ defines an isomorphism $C:\mathbb{R}^{p(\ell-1)}\to\ker B$ and we may thus reformulate the linear system of equations as

[TABLE]

With the solution $\widehat{x}$ we obtain the solution $x$ of the orginal system (${\rm S}$) via $x=C\widehat{x}$. Since the columns of $C$ are linearly independent the matrix $C$ has full rank. Hence, the symmetric matrix $\widehat{A}=C^{\sf T}AC$ is positive definite and the linear system of equations $\widehat{A}\widehat{x}=\widehat{f}$ can be solved efficiently with a preconditioned conjugate gradient scheme. The construction of suitable preconditioners has been discussed in [45] and the main issue in their justification is to understand how accurate the approximation

[TABLE]

is. Straightforward manipulations show that the solution $x$ of the saddle-point system (${\rm S}$) satisfies

[TABLE]

With the equvialent formulation $(C^{\sf T}AC)\widehat{x}=C^{\sf T}f$ and $x=C\widehat{x}$ we find that

[TABLE]

Since $C$ is orthogonal in the sense that $C^{\sf T}C=I_{p(\ell-1)}$ it follows that

[TABLE]

This justifies the approximation $(C^{\sf T}AC)^{-1}\approx C^{\sf T}A^{-1}C$ if the second term on the right-hand side is small. Another relation is obtained by formally applying the Neumann series

[TABLE]

to the product $T=\varrho(C^{\sf T}AC)(C^{\sf T}A^{-1}C)$ with a suitable parameter $0<\varrho\leq 1$ to deduce that

[TABLE]

Accepting the approximation $(C^{\sf T}AC)^{-1}\approx C^{\sf T}A^{-1}C$ the next step is to choose a preconditioner $P$ for $A$, i.e., $P$ approximates $A^{-1}$ and the multiplication by $P$ is inexpensive, and use the matrix $\widehat{P}=C^{\sf T}PC$ as a preconditioner for $\widehat{A}$. Noting that $C$ is orthogonal this is justified by the spectral norm estimate from [45],

[TABLE]

A further discussion of related preconditioners in the context of micromagnetics can be found in [41]. We illustrate the construction of the bases of the null space and the performance of different solution strategies in the context of harmonic maps into spheres.

5.1. Application to harmonic maps

Harmonic maps are vector fields with values in a given manifold which are stationary for the Dirichlet energy. In the case of the unit sphere, the vector field satisfies a pointwise unit-length constraint. Already this simple case illustrates analytical difficulties when dealing with constrained partial differential equations. In particular, harmonic maps are nonunique even for fixed boundary values and may be discontinuous everywhere, cf. [50]. Partial regularity results can be proved if a harmonic map is energy minimizing, cf. [56]. These observations underline the importance of computing harmonic maps with low energy. We aim at applying the concepts of the previous sections to the approximation of harmonic maps and consider the following model problem defined with a function $u_{\rm D}\in C(\partial\Omega;\mathbb{R}^{d})$ with $|u_{\rm D}(x)|=1$ for all $x\in\partial\Omega$ which is the trace of a function $\widetilde{u}_{\rm D}\in H^{1}(\Omega;\mathbb{R}^{d})$:

[TABLE]

Following [1, 6] we discretize the admissible set $\mathcal{A}$ by piecewise affine vector fields and impose the unit-length constraint in the nodes of the underlying triangulation. This leads to the following discrete minimization problem:

[TABLE]

For a function $u_{h}\in\mathcal{S}^{1,0}(\mathcal{T}_{h})^{d}$ with nonvanishing nodal values we define

[TABLE]

With a gradient flow approach and a linearized treatment of the nodal constraint we are led to the following algorithm.

Algorithm 5.1** (Gradient descent for harmonic maps).**

*Choose $u_{h}^{0}\in\mathcal{A}_{h}$ and $\tau>0$, set $k=1$.

(1) Compute $d_{t}u_{h}^{k}\in\mathcal{F}_{h}[u_{h}^{k-1}]$ such that for all $w_{h}\in\mathcal{F}_{h}[u_{h}^{k-1}]$ we have*

[TABLE]

(2) Stop if $\|d_{t}u_{h}^{k}\|_{\star}\leq\varepsilon_{\rm stop}$; otherwise, increase $k\to k+1$ and continue with (1).

The Lax–Milgram lemma shows that the iteration is well defined. Because of the orthogonality condition we have that

[TABLE]

which by induction gives $|u_{h}^{k}(z)|\geq\dots\geq|u_{h}^{0}(z)|=1$ for $k=0,1,\dots,K$. If the restriction of the inner product $(\cdot,\cdot)_{\star}$ to $\mathcal{S}^{1,0}(\mathcal{T}_{h})^{d}$ is represented by the matrix $M$ then the linear systems in the iteration can be written as

[TABLE]

with the vectorial finite element element stiffness matrix $S$ and the constraint matrix

[TABLE]

To define the matrix $C$ that provides an isomorphism from $\mathbb{R}^{p(d-1)}$ onto the kernel of $B$ we need to compute for a given vector $b\in\mathbb{R}^{d}\setminus\{0\}$ an orthonormal basis $(c_{2},\dots,c_{d})$ of its orthogonal complement. We proceed as follows: if $b$ is parallel to $e_{1}$ then we choose the remaining canonical basis vectors $(e_{2},\dots,e_{d})$; otherwise, we set

[TABLE]

where $b^{\perp}$ denotes the rotation of $b$ by $\pi/2$, followed by the normalization $c_{j}=\widetilde{c}_{j}/|\widetilde{c}_{j}|$ for $j=2,\dots,d$.

Example 5.2**.**

Let $\Omega=(-1/2,1/2)^{3}$ and define $u_{\rm D}$ on ${\Gamma_{\rm D}}=\partial\Omega$ for $x\in\partial\Omega$ via

[TABLE]

The employed triangulations $\mathcal{T}_{\ell}$ result from $\ell$ uniform refinements of a reference triangulation of $\Omega$ into six tetrahedra. The starting value $u_{h}^{0}\in\mathcal{S}^{1,0}(\mathcal{T}_{\ell})^{3}$ is defined by generating random nodal values of unit length.

Table 2 compares the following general strategies to solve the linear problems in the iterative solution of harmonic maps:

• (1)

Solve the original indefinite saddle-point system (${\rm S}$) with a direct solution method for sparse systems.

• (2)

Solve the reduced positive definite system (${\rm R}$) with a direct solution method for sparse linear systems.

• (3)

Solve the reduced positive definite system (${\rm R}$) with the preconditioned conjugate gradient method using a diagonal preconditioning of the full system matrix $C^{\sf T}AC$.

• (4)

Solve the reduced system (${\rm R}$) with the preconditioned conjugate gradient method using the preconditioner $C^{\sf T}(L_{\rm ic}L_{\rm ic}^{\sf T})^{-1}C$ with the incomplete Cholesky factorization $A\approx L_{\rm ic}L_{\rm ic}^{\sf T}$ computed once.

We measured the average time needed for different discretizations to solve the linear systems of equations using strategies (1)-(4). We always chose the step size $\tau=h$, the stopping parameter $\varepsilon_{\rm pcg}=10^{-8}$ for the relative residuals in the preconditioned conjugate gradient scheme, and the stopping criterion $\varepsilon_{\rm stop}=h_{\ell}/10$ in Algorithm 5.1. We employed the $H^{1}$ inner product to define $(\cdot,\cdot)_{\star}$ and used Matlab’s backslash operator as a model direct solver for sparse linear systems. From the obtained numbers we see that the preconditioner particularly designed for the structure of the constrained problems with changing constraints clearly outperforms all other approaches. However, it does not lead to mesh-size independent iteration numbers for this example with a nonsmooth solution. We also remark that we did not observe a further improvement when additional terms in the Neumann series were used with a damping factor $\varrho=1$. In a two-dimensional setting with a smooth solution the complete Cholesky factorization led to mesh-independent iteration numbers.

6. Applications, modifications, and extensions

We address in this section the application of the developed methods to the simulation of extended problems. In particular, we devise a numerical method for simulating bilayer plates, discuss how injectivity can be enforced in the case of elastic rods, and investigate a model that involves the thickness parameter and thereby allows for describing situations in which low energy membrane and bending effects occur simultaneously.

6.1. Bilayer plates

An important class of applications of nonlinear plate bending is related to composite materials, i.e., structures such as bilayer plates which are made of two sheets with slightly different mechanical features that are glued on top of each other. The difference in the material properties, e.g., concerning the response to environmental changes such as temperature, allows for externally controlled large bending effects. In bimetal strips one of the metals contracts while the other one expands leading in combination to the formation of rolls with prescribed curvature.

The preferred curvature is defined by the difference in elastic properties of the involed materials and is explicitly visible in the dimensionally reduced model identified in [53]. Following that work we consider a thin plate $\Omega_{\delta}=\omega\times(-\delta/2,\delta/2)$ and the inhomogeneous energy density

[TABLE]

with a material parameter $\alpha\in\mathbb{R}$. The energy density $W$ is assumed to satisfy the standard requirements stated in Section 2 so that $W_{\rm bil}$ is minimal for deformation gradients $F$ which are multiples of rotations with $\det F=(1+\delta\alpha)^{3}>1$ if $x_{3}>0$ and $\det F=(1-\delta\alpha)^{3}<1$ for $x_{3}<0$ corresponding to expansive and contractive behavior in the upper and lower sheets, respectively. It has been shown in [53] that for deformations with energy proportional to $\delta^{3}$ a limiting, dimensionally reduced problem is given by the following minimization problem:

[TABLE]

Here, $Q_{\rm plate}$ is the reduced quadratic form introduced in Section 2.2 which is applied to the difference of the second fundamental form $II$ and a multiple of the identity matrix with a factor given by half of the jump in the material difference divided by $\delta$. This material mismatch introduces a term that is similar to so-called spontaneous curvature terms in the context of biomembranes described via the Helfrich-Willmore model. Note that in the definition of $W_{\rm bil}$ the difference is $2\delta\alpha$ and its proportionality to the thickness $\delta$ of the plates is important to avoid delamination effects. The reduced model prefers deformations $y$ that define surfaces with mean curvature $\alpha$ and vanishing Gaussian curvature due to the isometry constraint. Hence, rolling effects in one direction occur. In [54] global minimizers for $I_{\rm bil}$ have been shown to be given by cylinders with radius $r=1/\alpha$. Recalling that the quadratic form $Q_{\rm plate}$ is given by the scalar product

[TABLE]

we do not have the property that the integrand is proportional to $|D^{2}y|^{2}$ as in the case of a single layer or two identical materials corresponding to the case $\alpha=0$. Instead we infer

[TABLE]

Because of the isometry condition the first term is proportional to $|D^{2}y|^{2}$ while the third term is a constant $c_{\alpha}$ which is irrelevant for the minimization process. For the second term on the right-hand side we have

[TABLE]

Therefore, using $II_{ij}=-\partial_{i}\partial_{j}y\cdot(\partial_{1}y\times\partial_{2}y)$ we have

[TABLE]

Convergent finite element discretizations and iterative strategies for the numerical solution of the bilayer plate bending problem have been devised in the articles [14, 13] extending the ideas described in the Sections 3.2 and  4.2. The iterative scheme of [14] is unconditionally stable but requires a subiteration, i.e., the solution of a nonlinear system of equations in every time step, which limits its practical applicability, while the scheme used in [13] is explicit and efficient but can only be expected to be conditionally stable. We devise here a new semi-implicit scheme with improved stability properties. The set of admissible deformations $\mathcal{A}_{h}$ and the tangent spaces $\mathcal{F}_{h}[y_{h}]$ are defined as in Section 3.2 and we set

[TABLE]

for a function $v_{h}\in\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})$ with a componentwise application in the case of a vector field. The discretized energy functional describing actuated bilayer plates is therefore given by

[TABLE]

A uniform discrete coercivity property of this discretization has been established in [14]. The new iterative minimization from [16] uses an explicit treatment of the spontaneous curvature term.

Algorithm 6.1** (Gradient descent for bilayer plates).**

*Choose $y_{h}^{0}\in\mathcal{A}_{h}$ and $\tau>0$, set $k=1$.

(1) Compute $d_{t}y_{h}^{k}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ such that for all $w_{h}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ we have*

[TABLE]

(2) Stop if $\|d_{t}y_{h}^{k}\|_{\star}\leq\varepsilon_{\rm stop}$; otherwise increase $k\to k+1$ and continue with (1).

We show that the iteration of Algorithm 6.1 is energy stable on finite time intervals under a moderate condition on the step-size resulting from the use of an inverse estimate. To formulate it, we first note that we assume that the boundary conditions imply a Poincaré type estimate

[TABLE]

for all $w_{h}\in\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})^{3}$ with $L_{\rm bc}^{\rm bil}[w_{h}]=0$. With this estimate one verifies that the mesh-dependent estimate

[TABLE]

holds for all $w_{h}\in\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})^{3}$ with $L_{\rm bc}^{\rm bil}[w_{h}]=0$ and with the minimal mesh-size $h_{\rm min}$.

Proposition 6.2** (Convergent iteration).**

The iterates of Algorithm 6.1 are well defined. Assume that we have

[TABLE]

for all $w_{h}\in\mathcal{S}^{\rm dkt}(\mathcal{T}_{h})^{3}$ with $L_{\rm bc}^{\rm bil}[w_{h}]=0$. Then there exist constants $c_{\rm bil},c_{\rm bil}^{\prime}>0$ such that if $\tau|\log h_{\rm min}|^{2}\leq c_{\rm bil}^{\prime}$ for $L=0,1,\dots,K\leq T/\tau$ we have

[TABLE]

and

[TABLE]

where $e_{0,h}=I_{\rm bil}^{h}[y_{h}^{0}]$.

Proof.

We argue by induction and assume that the estimate has been established for $L-1$. Choosing $w_{h}=d_{t}y_{h}^{k}$ and using Hölder’s inequality shows that for $k\leq L$ we have

[TABLE]

Absorbing terms involving $d_{t}y_{h}^{k}$ on the left-hand side and noting that terms involving $y_{h}^{k-1}$ are bounded because of the bounds for $k\leq L-1$, we find that

[TABLE]

Here, we also used that the discrete energy is coercive. With this intermediate estimate and the orthogonality condition encoded in the definition of $\mathcal{F}_{h}[y_{h}^{k-1}]$ we infer with $y_{h}^{k}=y_{h}^{k-1}+\tau d_{t}y_{h}^{k}$ that

[TABLE]

and hence with the inverse estimate we obtain the suboptimal estimate that

[TABLE]

To obtain the energy bound we note that a discrete product rule shows that the discrete time derivative of the spontaneous curvature term is given by

[TABLE]

Comparing this expression with the right-hand side of the equation in Step (1) of Algorithm 6.1 with $w_{h}=d_{t}y_{h}^{k}$ leads to

[TABLE]

We use the estimates $\|\partial_{j}d_{t}y_{h}^{k}\|_{L^{\infty}_{h}}\leq c|\log h_{\rm min}|\|d_{t}y_{h}^{k}\|_{\star}$ and $\|\partial_{j}y_{h}^{k-\ell}\|_{L^{\infty}_{h}}\leq c$ for $j=1,2$ and $\ell=0,1$, to bound the right-hand side by $\tau|\log h_{\rm min}|c^{\prime\prime\prime}\|d_{t}y_{h}^{k}\|_{\star}^{2}$. This leads to

[TABLE]

and proves the asserted energy estimate. With this bound we also obtain the improved bound for the constraint violation error. ∎

The good stability properties of the newly proposed numerical scheme are confirmed by a numerical experiment whose outcome is visualized in Figure 10. The setup uses a rectangular strip of length $L=10$ and width $w=4$ that is clamped at one end. The spontaneous curvature parameter is $\alpha=-1$ and we have $c_{\rm b}=1$ and $c_{\rm sc}=1$. The figure shows snapshots of the evolution on a grid with medium mesh-size and stationary states for different triangulations with mesh-sizes $h$ proportional to $2^{-\ell}$, $\ell=1,2,3$. The step-size was always set to $\tau=h/20$. The stopping criterion was chosen as $\varepsilon_{\rm stop}=10^{-3}$. Because of the extreme geometry and the large deformation the asymmetry of the underlying triangulations is reflected in the numerical solutions but this effect disappears for smaller mesh-sizes. The theoeretical results about the energy monotonicity and controlled constraint violation are confirmed by the bottom plot of Figure 10.

6.2. Selfavoiding curves and elastic knots

A strategy for finding useful representatives of knot classes, i.e., closed curves within a given isotopy class with particular features, is to minimize or decrease the bending energy within the given class via continuous evolutions defined by a gradient flow. To ensure that the flow does not change the topological properties of the curve, appropriate terms have to be included in the mathematical formulation. To this end we add a self-avoidance potential to the bending energy that prevents the curve from self-intersecting or pulling tight, i.e., we consider flows determined by functionals

[TABLE]

in sets of curves satisfying periodic or, e.g., clamped boundary conditions. We refer the reader to [46] for a general discussion of appropriate functionals. A functional that turned out to have several advantageous features is the tangent-point functional proposed in [39]. It is for a $C^{1}$ curve $y:[0,L]\to\mathbb{R}^{3}$ defined via the tangent-point radius $r_{y}(x,z)$ which is the radius of the circle that is tangent to $y$ in $y(x)$ and intersects the curve in the point $y(z)$. For curves that are parametrized by arclength we have

[TABLE]

which follows from investigating the geometrical configuration sketched in Figure 11. If $z\to x$ then we have that $r_{y}(x,z)$ converges to the inverse of the curvature of the curve $y$ at $x$. If otherwise $y(x)\to y(z)$ for fixed $x\neq z$ then the radius converges to zero as depicted in Figure 12. An alternative choice is the use of the Menger curvature which is discussed in [59].

The tangent-point functional ${\rm TP}$ is for an exponent $q\geq 1$ defined via

[TABLE]

The exponent has to be suitably chosen so that singularities corresponding to a vanishing tangent-point radius are sufficiently strong to lead to an infinite value of the functional. Thereby, an energy barrier is defined that separates different isotopy classes. The functional ${\rm TP}$ has important and remarkable features.

Remarks 6.3**.**

*(i) If the arclength-parametrized curve $y\in C^{1}([0,L];\mathbb{R}^{3})$ is injective then ${\rm TP}[y]$ is finite if and only if $y\in W^{2-1/q,q}(0,L;\mathbb{R}^{3})$, cf. [22].

(ii) For all $M>0$ and $q>2$ there exists a constant $c_{M,q}>0$ such that for all arclength-parametrized curves $y\in C^{1}([0,L];\mathbb{R}^{3})$ with ${\rm TP}[y]\leq M$ we have the bi-Lipschitz estimate $|x-z|\leq c_{M,q}|y(x)-y(z)|$ for all $x,z\in[0,L]$, cf. [23]. In particular, ${\rm TP}$ is a knot energy in the sense that for a family of curves $(y_{k})_{k\in\mathbb{N}}$ converging pointwise to a curve with self-intersection, the values ${\rm TP}[y_{k}]$ blow up, cf. [63, 46, 60, 23].

(iii) The functional ${\rm TP}$ is twice continuously differentiable with bounded variations and $L^{1}$ integrands, cf. [23, 17].*

To decrease the total energy $I^{\rm tot}$ of a given initial curve $y_{0}$ within its isotopy class we use a discretization of the evolution

[TABLE]

subject to the linearized arclength conditions $\partial_{t}y^{\prime}\cdot y^{\prime}=0$ and $w^{\prime}\cdot y^{\prime}=0$. The discretization uses the ideas outlined in Sections 3.1 and 4.1 with an explicit treatment of the potential.

Algorithm 6.4** (Gradient descent for selfavoiding curves).**

*Choose an initial $y_{h}^{0}\in\mathcal{A}_{h}$ and a step-size $\tau>0$, set $k=1$.

(1) Compute $d_{t}y_{h}^{k}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ such that for all $w_{h}\in\mathcal{F}_{h}[y_{h}^{k-1}]$ we have*

[TABLE]

(2) Stop the iteration if $\|d_{t}y_{h}^{k}\|_{\star}\leq\varepsilon_{\rm stop}$; otherwise, increase $k\to k+1$ and continue with (1).

The explicit treatment of the potential has several advantages. First, we obtain linear problems in the time steps. Second, we avoid the inversion of fully populated matrices. Third, the assembly of the vector on the right-hand side can be fully parallelized without communication costs. Finally, we do not change the stability of the iteration compared to a fully implicit time stepping scheme since the potential does not have obvious convexity properties. Using the differentiability properties of the functional ${\rm TP}$ and assuming that the flow metric is the $H^{2}$ scalar product, it has been shown in [17] that the conditional energy decay property

[TABLE]

holds for all $0\leq L\leq K\leq T/\tau$ with a finite time horizon $T>0$. This inequality implies the convergence of the iteration to a stationary configuration. Since the total energy is bounded uniformly, the values of the potential remain controlled so that in case of a sufficient resolution no self-contact or topology changes can occur.

The good stability properties and corresponding selfavoidance behavior of the numerical scheme are illustrated by the snapshots of configurations and the corresponding energy decay shown for a generic setting in Figure 13. We observe that the curve relaxes to a curve with equilibrated curvature and which is nearly flat. While the total energy decreases the (unscaled) tangent-point functional increases. In the experiments we used the model parameters $c_{\rm b}=10$, $\varrho=10^{-3}$, $q=3.9$. The initial curve belongs to the knot class $8_{10}$, cf. [51], and has a total length $L\approx 45.45$ which is accurately preserved during the evolution. To compute the evolution we used a partitioning of the reference interval into $551$ subintervals. The step-size $\tau$ was chosen proportional to the mesh-size $h$.

6.3. Föppl–von Kármán model

The nonlinear bending model discussed in the previous sections is not suitable to describe certain phenomena such as the formation of wrinkles which occur for very small energies. In order to describe these effects in a dimensionally reduced model it is essential to involve the thickness parameter $\delta>0$ as this number determines the oscillating or nonsmooth behavior of deformations. A model that captures such effects is the Föppl–von Kármán model which can be identified via different scalings of the in-plane and out-of-plane components of a three-dimensional deformation when $\delta$ is small. It determines a planar displacement $u:\omega\to\mathbb{R}^{2}$ and a deflection $w:\omega\to\mathbb{R}$ as a minimizing pair for the energy functional

[TABLE]

in a set $\mathcal{A}$ of admissible pairs contained in the product space $H^{1}(\omega;\mathbb{R}^{2})\times H^{2}(\omega)$. Here, we use the symmetric gradient

[TABLE]

and the dyadic product

[TABLE]

where $\nabla w$ is understood as a column vector. We refer the reader to [26, 27, 35, 37] for justifications of this minimization problem as a simplification of three-dimensional hyperelastic material models. It has been shown rigorously in [44, 62] that the presence of the parameter $\delta$ leads to minimizers with wrinkling patterns whose geometry is determined by the parameter $\delta$. This is done by identifying optimal scaling laws of the energy in terms of $\delta$, cf. [20, 30] for further examples. Only a few numerical methods for minimizing the Föppl–von Kármán energy are available, cf. [28] for an abstract investigation.

To minimize the energy for given boundary conditions we follow [12] and again adopt a gradient flow appraoch defined by the system

[TABLE]

Here, $(\cdot,\cdot)_{\star}$ and $(\cdot,\cdot)_{\dagger}$ are inner products on $H^{2}(\omega)$ and $H^{1}(\omega;\mathbb{R}^{2})$, respectively, and we used the identities $|a\otimes a|^{2}=|a|^{4}$ and $M:(a\otimes b)=(Ma)\cdot b=(Mb)\cdot a$ for a symmetric matrix $M\in\mathbb{R}^{2\times 2}$ and $a,b\in\mathbb{R}^{2}$. The temporal discretization of the system decouples the equations via a semi-implicit evaluation of different terms. In particular, given $(u^{k-1},w^{k-1})$ we compute $(u^{k},w^{k})$ such that

[TABLE]

for all $(v,z)$ satisfying appropriate homogeneous boundary conditions. The average $w^{k-1/2}$ is defined for subsequent approximations $w^{k}$ and $w^{k-1}$ via

[TABLE]

The unconditional stability and energy decay of the iteration follows from choosing $v=d_{t}w^{k}$ and $z=d_{t}u^{k}$ and exploiting a discrete product rule, cf. [12] for details. The problems in the time steps determine minimizers of certain functionals and if $0<\tau\leq\tau_{0}$ then these functionals are strongly convex and minimizers are uniquely defined. Correspondingly, we expect the Newton scheme to converge provided that the step-size $\tau$ is sufficiently small.

For a full discretization we use the set of admissible pairs

[TABLE]

and the corresponding homogeneous space

[TABLE]

To avoid an unnecessary small uniform step-size $\tau$ we instead use variable step-sizes $(\tau_{k})_{k=0,1,\dots}$ that are adjusted according to the performance of the Newton scheme using the following rules:

[TABLE]

The parameter $r\geq 0$ defines an upper bound for the step-sizes and avoids a numerical overflow. The ideas lead to the following algorithm.

Algorithm 6.5** (Gradient descent for Föppl–von Kármán functional).**

*Choose $(u_{h}^{0},w_{h}^{0})\in\mathcal{A}_{h}$, an integer $N_{\rm max}>0$, stopping tolerances $\varepsilon_{\rm stop},\varepsilon_{\rm N}>0$, and an initial step-size $\tau_{1}>0$, set $k=1$.

(1a) Repeatedly decrease $\tau_{k}$ until the Newton scheme terminates within $N_{\rm max}$ steps and tolerance $\varepsilon_{\rm N}$ to determine $d_{t}w_{h}^{k}\in\mathcal{F}_{h,w}^{0}$ such that*

[TABLE]

*for all $v_{h}\in\mathcal{F}_{h,w}^{0}$.

(1b) Compute $d_{t}u_{h}^{k}\in\mathcal{F}_{h,u}^{0}$ such that*

[TABLE]

*for all $z_{h}\in\mathcal{F}_{h,u}^{0}$.

(2) Stop if $\|d_{t}w_{h}^{k}\|_{\star}+\|d_{t}\widetilde{\varepsilon}(u^{k}_{h})\|_{\dagger}\leq\varepsilon_{\rm stop}\min\{1,\tau_{k}\}$; otherwise, define*

[TABLE]

increase $k\to k+1$, and continue with (1).

Note that in the algorithm a stopping criterion is used that is proportional to the step-size which is important in the case of small step-sizes. To illustrate the performance of the algorithm and features of the mathematical model we follow [12] and consider the compression of a plate along one of its sides. Particularly, we let $\omega=(-1/2,1/2)\times(0,1)$ and compress the side $\{0\}\times[-1/2,1/2]$ by 10 percent which defines the boundary data for the in-plane displacement $u$. We use homogeneous clamped boundary conditions along the same side for the deflection $w$. The results shown in Figure 14 confirm the formation of wrinkling structures. These depend on both the numerical resolution and the thickness parameter. Moreover, the energy decay shown in the bottom plot of Figure 14 indicates that these are only obtained within a reasonable number of iterations if an adaptive step-size strategy is used.

7. Conclusions

We have addressed in this article the accurate finite element discretization and reliable iterative solution of nonlinear models for describing large deformations of elastic rods and plates. To avoid unjustified regularity assumptions we have adopted the concept of $\Gamma$-convergence and thereby proved convergence of approximations. To compute stationary configurations of low and possibly minimal energy we have used gradient flow discretizations that are guaranteed to decrease the elastic energy and preserve inextensibility and isometry constraints appropriately. These two concepts leave a theoretical gap in the sense that the framework of $\Gamma$-convergence is concerned with global energy minimizers while gradient flows can only be expected to determine stationary configurations. However, $\Gamma$-convergence goes beyond global energy minimization and gradient flows typically avoid unstable critical configurations. These properties are convincingly confirmed by our numerical experiments which did not employ any a priori knowledge about expected solutions. The state of the art of reliable methods for computing nonlinear bending phenomena leaves open several aspects such as the development of optimal preconditioners, adaptive refinement of stress concentrations, determination of convergence rates, effective combination with membrane effects, or simulation of realistic dynamic models open. We believe that the methods presented in this article can be useful in investigating them in future research.

Acknowledgments. The author wishes to thank his coworkers in various collaborations that led to the results presented in this article.