Semi-implicit methods for the dynamics of elastic sheets
Silas Alben, Alex A. Gorodetsky, Donghak Kim, Robert D. Deegan
SSemi-implicit methods for the dynamics of elastic sheets
Silas Alben , ∗ , Alex A. Gorodetsky , Donghak Kim , and Robert D. Deegan Department of Mathematics, Department of Aerospace Engineering, Department of Physics & Center for the Study of Complex SystemsUniversity of Michigan, Ann Arbor, MI 48109, USA ∗ Recent applications (e.g. active gels and self-assembly of elastic sheets) motivate the need toefficiently simulate the dynamics of thin elastic sheets. We present semi-implicit time steppingalgorithms to improve the time step constraints that arise in explicit methods while avoiding much ofthe complexity of fully-implicit approaches. For a triangular lattice discretization with stretching andbending springs, our semi-implicit approach involves discrete Laplacian and biharmonic operators,and is stable for all time steps in the case of overdamped dynamics. For a more general finite-difference formulation that can allow for general elastic constants, we use the analogous approachon a square grid, and find that the largest stable time step is two to three orders of magnitudegreater than for an explicit scheme. For a model problem with a radial traveling wave form of thereference metric, we find transitions from quasi-periodic to chaotic dynamics as the sheet thicknessis reduced, wave amplitude is increased, and damping constant is reduced.
I. INTRODUCTION
In recent years there have been various studies of how spatial variations in the composition of a thin sheet canproduce global conformational changes. Examples include the appearance of spontaneous curvature due to strainvariations across the thickness of the sheet [1–3] or non-Euclidean reference metrics induced by in-plane strain varia-tions [4–8]. A theory of incompatible elastic plates [9, 10] has been developed to determine equilibrium configurationsof such sheets. Related approaches have been used to develop self-folding origami gels [11]. A catalog of respon-sive materials are now available for investigating the mechanics of thin sheets, including non-uniform responsive gelsheets [12–14], sheets of nematic elastomers with a non-uniform director field [15, 16], responsive gels combined withoriented micro-rods [17], and sheets in confined geometries [18–21].The dynamics of responsive gel sheets were the focus of work by Yoshida and collaborators, who synthesized a gelthat locally swells in response to chemical waves propagating entirely within the gel [22]. They used self-oscillatinggels in various narrow strip geometries to make a variety of soft machines [23–27]. The Balazs group used poroelasticsimulations of self-oscillating gels to demonstrate additional examples of soft machines with uniaxial bending orisotropic swelling [28–31]. Here we will focus on simulating a simple but extensively-studied model, a thin elasticsheet driven by changes in its equilibrium metric. We will present an efficient semi-implicit time-stepping algorithm ∗ Electronic address: [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] A p r for the case of overdamped dynamics, a representative case with the same form of numerical stiffness as more detailedfluid-elastic and fluid-poroelastic models. The approach is also useful in other applications where the dynamics ofthin elastic sheets are important, such as the self-assembly of thin sheets under magnetic forces [32, 33], and the therolling of actuated bilayers [34, 35].A related problem is simulating the dynamics of fluid membrane vesicles with surface tension [36]. Here a similartime-step constraint arises for bending forces, while the thin elastic sheets considered in the present work also havestiffness due to stretching forces. Like [36], we develop a semi-implicit time-stepping approach for computationalefficiency, though our formulation differs due to the different mechanical forces.In general, a semi-implicit (or implicit-explicit) time discretization writes some of the terms (typically those withthe highest spatial derivatives) implicitly, to improve time-step constraints for stability [37–45]. The implicit termsare typically linear in the unknowns at the current time step, so they can be solved directly at each time step, avoidingsome of the complication and computational expense of nonlinear iterative solvers (e.g. Newton-type methods) infully implicit discretizations [46]. If the linearized implicit term is sufficiently large in comparison to the explicit terms,the semi-implicit method may be stable for a wide range of time steps. For nonlinear PDEs, a somewhat empiricalapproach to formulating schemes, based on analogies with time-stepping for simpler linear PDEs, is often necessary.In this work we begin with Seung and Nelson’s discretization of elastic sheets by bending and stretching springson a triangular lattice [47]. We use the approach of [40, 43] to split the stretching force into an implicit linear termcorresponding to zero-rest-length springs, and a nonlinear remainder. The implicit stretching term is proportionalto a discrete Laplacian matrix multiplying the current sheet position. To stabilize the bending force, we add animplicit bending term that is proportional to a discrete biharmonic matrix multiplying the current sheet position.The resulting method appears to be stable for all time steps. To allow for general elastic constants, we formulatea finite-difference discretization of the elastic energy with the analogous semi-implicit approach. We validate andcompare the methods on test problems with internal in-plane stretching forces (nontrivial equilibrium metrics) andstudy the effects of basic physical parameters on the sheets’ dynamics. II. ELASTIC SHEET
We consider a thin sheet or bilayer that undergoes large time-dependent deformations due to internal forces (froma prescribed, time- and space-varying reference metric). We assume the sheet obeys linear (Hookean) elasticity, butthat the midsurface (the set of points located midway through the sheet in the thickness direction) can be an arbitrarysmooth time-dependent surface, so elastic forces depend nonlinearly on its position. The extension of the Kirchhoff-Love (and F¨oppl-von-K´arm´an) models of elastic plates to nonflat reference metrics and/or large deformations withsmooth midsurfaces has been called the Koiter shell theory [48–50] or the theory of non-Euclidean plates [9, 51]. Herewe mainly follow the latter’s notation.The elastic energy involves stretching and bending energy terms determined by the position of the body, r ( x ).Here r lies in R , as does the material coordinate x = ( x , x , x ). In the classical situation, the sheet has a zero-energy flat state r ( x ) = x , where x lies in the interval [ − h/ , h/
2] ( h is the sheet thickness, much smaller thanthe other dimensions), and x and x lie in a planar region, the same for each x . For a small line of material ∆x = ˜ x − x connecting two material points ˜ x and x , its squared length on the undeformed surface is dl = ∆ x i ∆ x i (with summation over repeated indices). Denote the squared length on the deformed surface r ( x ) by dl (cid:48) . Using theTaylor series of r ( x ) up to first derivatives we have dl (cid:48) − dl = 2 (cid:15) ij ∆ x i ∆ x j . (1)Here (cid:15) ij is the strain tensor, (cid:15) ij = 12 ( g ij − δ ij ) . (2)where g ij = ∂r k ∂x i ∂r k ∂x j (3)is the metric tensor. For curved shells or active materials, the rest state may be curved and/or time-varying, in whichcase the reference metric, δ ij in (2), becomes ¯ g ij ( x , t ), a time- and space-varying function determined for example bychemical activity [12, 22, 28]: (cid:15) ij = 12 ( g ij − ¯ g ij ( x , t )) . (4)The reference metric is assumed to take the form¯ g = ¯ g ¯ g g ¯ g
00 0 1 (5)with upper 2-by-2 reference metric ¯ g αβ . The dependence on out-of-plane components (¯ g α , ¯ g α ) is trivial, so shearingthough the plate thickness is not imposed and expansion/contraction in the thickness direction occurs only passively,due to the Poisson ratio effect. The components of the metric g i for i = 1 , , g α = g α = 0 for α = 1 ,
2, and g = 1 [9].We may write the energy in terms of the midsurface deformation by expanding g αβ about the sheet midsurface x = 0. We obtain, at leading order, its thickness-average a αβ and its thickness-gradient [9]: g αβ = a αβ − x b αβ + O ( h ) . (6)Here a αβ is the upper 2-by-2 part of the metric tensor (3) evaluated at the sheet midsurface, x = 0. The thickness-gradient is written in terms of b αβ , the second fundamental form b αβ = ∂ r k ∂x α ∂x β n k (7)also evaluated at the sheet midsurface, x = 0, with n k the components of the midsurface unit normal vector n . Forthe reference metric ¯ g , we write ¯ a and ¯ b for the corresponding terms in the expansion about the midsurface.For an isotropic sheet with Young’s modulus E , Poisson ratio ν , and thickness h , the elastic energy per unit volumeis w = 12 ¯ A αβγδ (cid:15) αβ (cid:15) γδ , (8)a quadratic function of the in-plane components of the strain tensor in (4) with elasticity tensor¯ A αβγδ = E ν (cid:18) ν − ν ¯ g αβ ¯ g γδ + ¯ g αγ ¯ g βδ (cid:19) (9)where the entries of ¯ g αβ are those of ¯ g − αβ . Integrating w over the sheet thickness, the energy per unit midsurface areais w D = (cid:90) h/ − h/ wdx = w s + w b + h.o.t., (10)a sum of stretching energy per unit area w s = h A αβγδ ( a αβ − ¯ a αβ ) ( a γδ − ¯ a γδ ) (11)and bending energy per unit area w b = h A αβγδ (cid:0) b αβ − ¯ b αβ (cid:1) (cid:0) b γδ − ¯ b γδ (cid:1) . (12)in terms of the midsurface elasticity tensor A αβγδ = E ν (cid:18) ν − ν ¯ a αβ ¯ a γδ + ¯ a αγ ¯ a βδ (cid:19) , (13)and higher order terms in h . The total elastic energy is, to leading order in h , W = W s + W b , W s = (cid:90) (cid:90) w s (cid:112) | ¯ a | dx dx , W b = (cid:90) (cid:90) w b (cid:112) | ¯ a | dx dx . (14) W is a function of the midsurface configuration, r ( x , x ,
0) and midsurface reference metric ¯ a . III. SHEET DYNAMICS
The sheet midplane evolves according a force balance equation, where the elastic force per unit area f acting at apoint on the midplane is a sum of stretching and bending forces per unit area: f = f s + f b = δw s /δ r + δw b /δ r (15)given by taking the variation of w s + w b with respect to r . For a sheet moving in Stokes flow (i.e. at zero Reynoldsnumber) the elastic forces would be balanced by external fluid forces which depend linearly on the sheet velocity [36]: ∂ r ∂t = S [ f s + f b ]( r ) . (16)where S is the Stokes operator S [ f ]( r ) = (cid:90) G ( r , r (cid:48) ) f ( r (cid:48) ) dA ( r (cid:48) ) , G ( r , r (cid:48) ) = 18 πµ (cid:18) (cid:107) ρ (cid:107) I + ρ ⊗ ρ (cid:107) ρ (cid:107) (cid:19) , ρ ≡ r − r (cid:48) . (17)In this work, we consider a simplification of (16) which has similar numerical stiffness issues: overdamped dynamics,in which the Stokes operator is replaced with a multiple of the identity: µ ∂ r ∂t = f . (18)Here µ is a parameter that can be used to model the effects of internal and external damping. The equation can alsobe used to identify equilibria as a gradient descent method [52] with time playing a nonphysical role. Extensions of(18) including the inertia of the sheet or a surrounding fluid will add a dependence on ∂ r /∂t . When discretized,the ∆ t − dependence of such terms will improve the time-step constraint compared to the overdamped problem, forexplicit schemes [53]. In this work, we develop a semi-implicit approach for (18), which should also apply with moredetailed forms of internal and external damping/forcing, e.g. from a fluid. In the subsequent computational resultswe nondimensionalize the sheet lengths by the radius or half-width R (e.g. for hexagonal and square sheets), energyby the bending energy scale Eh /
12, and time by the period of ¯ g αβ ( x , t ) (assumed to be time-periodic). These choicesdefine dimensionless versions of all the parameters (e.g. µ ). IV. TRIANGULAR LATTICE WITH STRETCHING AND BENDING SPRINGS r j r i n l n k FIG. 1:
Triangular lattice mesh with examples of elements in the elastic energy: adjacent vertices r i and r j (red) and adjacent facenormals n k and n l (blue). We first consider a simple model of an elastic sheet with material points connected by an equilateral triangularlattice mesh. Nearest neighbor points are connected by Hookean springs and the total stretching energy is a sum ofthe squares of nearest neighbor distances minus the local spring rest length d ij : U s = K s (cid:88) i,j ( (cid:107) r i − r j (cid:107) − d ij ) , (19)with K s a stretching stiffness constant. A bending energy is applied to adjacent triangular faces based on the anglesbetween the normals to the faces. The total bending energy is a sum over nearest neighbor pairs with bending stiffnessconstant K b : U b = K b (cid:88) k,l (cid:107) n k − n l (cid:107) = K b (cid:88) k,l − n k · n l . (20)Seung and Nelson used this model to study buckling due to defects in elastic membranes [47], and it was used by manyother groups to study other deformations of thin sheets and shells due to defects and/or external forces [35, 54–61],as well as polymerized and fluid membranes [62].Seung and Nelson showed that for a lattice with d ij ≡ d , a constant, as d tends to 0 the stretching energy U s tends tothat of an isotropic thin sheet with stretching rigidity Eh = 2 K s / √ ν = 1 /
3. The continuum limitof the bending energy contains two terms, one proportional to the mean curvature squared and the other proportionalto the Gaussian curvature. The term involving mean curvature tends to that of an isotropic thin sheet with bendingrigidity Eh / − ν ) = √ K b /
2. However, with this bending rigidity, the prefactor of the Gaussian curvatureterm is too large for ν > − / ν = 1 /
3) [63]. Nonetheless, for many problems, theGaussian curvature term plays a negligible role because it can be integrated to yield only boundary terms. For closedshells, the Gaussian curvature term integrates to a constant and thus does not affect the elastic forces [56]. For opensheets (with boundaries), the equilibrium shape could be insensitive to the boundary shape or boundary conditions,particularly if the external or internal actuation is not localized at the boundary. In this work we will compare themodel to a finite-difference discretization with ν = 1 / { r i } . Thegradient of (19) with respect to r i is ∇ r i U s = K s (cid:88) j ∈ nhbrs( i )( (cid:107) r i − r j (cid:107) − d ij ) ( r i − r j ) (cid:107) r i − r j (cid:107) . (21)where nhbrs( i ) is the set of vertex neighbors to i . Following [40, 43] we write the summand as a linear term plus aterm with constant magnitude, ∇ r i U s = K s (cid:88) j ∈ nhbrs( i ) r i − r j − d ij ( r i − r j ) (cid:107) r i − r j (cid:107) . (22)To write the algorithms we define r = (cid:2) r (cid:124) x , r (cid:124) y , r (cid:124) z (cid:3) (cid:124) (23)as the vector of 3 N vertex coordinates, with r x , r y , and r z the N -vectors of x -, y -, and z -coordinates. Thus r i = (cid:2) ( r x ) i , ( r y ) i , ( r z ) i (cid:3) (cid:124) . Also note that the italicized r ∈ R N is different from the position function r ( x ) and avertex on the discretized surface r i , both of which take values in R . It turns out that the linear term in (22) can bewritten as the product of K s and a block diagonal matrix L with r , where L has three blocks (along the diagonal),each of which is a discretized Laplacian on the triangular mesh with free-edge boundary conditions (see examples ofstencils in figure 2, top row), a negative semidefinite matrix. Each block multiplies r x , r y , and r z , respectively. Wetreat the linear term implicitly and the constant-magnitude term explicitly. Collecting the terms (22) for all vertices i , we obtain the total stretching force. A semi-implicit first-order temporal discretization of (18) with stretching forcesonly is: µA p r n +1 − r n ∆ t = K s L r n +1 + f SE ( r n ) , [ f SE ( r ) i , f SE ( r ) N + i , f SE ( r ) N + i ] (cid:124) ≡ − K s (cid:88) j ∈ nhbrs( i ) d ij ( r i − r j ) (cid:107) r i − r j (cid:107) . (24)Here A p = n p √ d /
12 is the area per point on the undeformed lattice, with n p the number of triangles of which thepoint is a vertex, 6 for interior points and fewer for boundary points. f SE ( r ) is the nonlinear term in (22), with 3 N entries, given on the right side of (24) for i = 1 , . . . , N .Now assume the spring rest lengths are bounded for all time: d ij ≤ ¯ d , a constant, and each vertex has at most p neighbors (6 for the triangular lattice). Rearranging (24) and using the boundedness of f SE ( r ), we have an upperbound at time step n + 1: (cid:13)(cid:13) r n +1 (cid:13)(cid:13) ≤ (cid:13)(cid:13) ( I − ∆ tK s L / ( µA p )) − r n (cid:13)(cid:13) + (cid:13)(cid:13) ( I − ∆ tK s L / ( µA p )) − f SE ( r n ) (cid:13)(cid:13) (25) ≤ (cid:107) r n (cid:107) + K s p ¯ d (26) ≤ (cid:13)(cid:13) r (cid:13)(cid:13) + ( n + 1) K s p ¯ d. (27)So (cid:107) r n +1 (cid:107) grows at most linearly in time with this discretization. Empirically, the iteration appears to be boundedfor all time steps for spring rest lengths d ij that are bounded in time. For comparison, a forward Euler discretizationof (18) with stretching forces only results in r n +1 = ( I + ∆ tK s L / ( µA p )) r n + ∆ tµA p f SE ( r n ) . (28)Neglecting the rightmost term in (28), we have a 2D diffusion equation, and stability is possible only when the largesteigenvalue of I + ∆ tK s L / ( µA p ) is bounded in magnitude by 1. Since the eigenvalues of L ∼ / ∆ x for lattice spacing∆ x , this requires ∆ t < C s µA p ∆ x /K s for a constant C s .The gradient of the bending energy (20) with respect to a lattice vertex r i is ∇ r i U b = K b (cid:88) k,l sin θ kl ∇ r i θ kl , (29)using n k · n l = cos θ kl . The dihedral angle θ kl depends on the four points in the union of the neighboring triangles k and l (see figure 1). Two of these points are the endpoints of the edge shared by the triangles. At each of the othertwo points, ∇ r i θ kl is directed along the normal to the triangle in which it lies, with magnitude equal to the reciprocalof its distance from the shared edge. At the endpoints of the shared edge, ∇ r i θ kl can be found by requiring that thenet force and torque due to θ kl is zero, which gives six equations (for the three components of net force and torque)in six unknowns (the forces on the two endpoints of the shared edge). The bending force f B is a 3 N -vector withcomponents [ f B ( r ) i , f B ( r ) N + i , f B ( r ) N + i ] (cid:124) ≡ −∇ r i U b , i = 1 , . . . , N, (30)similar to f SE in (24). Our linearized approximation to the bending force is similar to that of [36]. They write theterms with the highest spatial derivatives in the form B ( r n ) r n +1 . Here B involves fourth derivatives with prefactorsthat include lower-order derivatives extrapolated to time step n + 1 from previous time steps. In fact, we use a simpler FIG. 2:
Examples of the stencils corresponding to a discrete Laplacian operator with free edge boundary conditions (one of the diagonalblocks of L defined below (22)) (top row), and ˜∆ x ,x , a discrete biharmonic operator with free edge boundary conditions (bottom row),at different mesh points (circled) away from and near a boundary on the triangular lattice. expression: B r n +1 , where B is a block diagonal matrix with each of the three blocks equal to D ˜∆ x ,x . Here D isthe bending modulus ( Eh / (12(1 − ν )) in dimensional form, 1 / (1 − ν ) in dimensionless form) and ˜∆ x ,x is thediscretized biharmonic operator on the triangular lattice in the orthogonal material coordinates x and x . If in-planestrain (shearing and dilation) are not very large, x and x are close to orthogonal arclength coordinates s , s alongthe midplane surface. For any surface X ( s , s ) parametrized by orthogonal arclength coordinates s , s we can write∆ s ,s X ( s , s ) = ∆ s ,s ( κ + κ )ˆ n + N. (31)where κ + κ is twice the mean curvature and N = ( κ + κ )∆ s ,s ˆ n involves derivatives of X that are of lower orderthan those in the first term on the right hand side. The highest-derivative term in the continuum bending force isalso the first term on the right hand side of (31) when the equilibrium metric is the identity (see [36, 65]). Thus theleft hand side of (31) is a reasonable linear (and constant-coefficient) approximation to the bending force. A moreaccurate linear approximation to the bending operator could include corrections that take into account the nontrivialinverse equilibrium metric (¯ g αβ in (13)) and in-plane strain, e.g. by including nonuniform prefactors extrapolated fromprevious time steps. We find however that the constant-coefficient biharmonic is a sufficiently good approximation inthe sense that it damps out spurious mesh-scale bending oscillations (as occurs with a fully explicit bending term) upto large-amplitude variations in the reference metric, ≈ . A defined in (34)–(36),below. With larger variations in the reference metric, the position vector remains bounded in time, but there is verylarge deformation and self-intersection even for the case of purely planar deformations, without bending forces. Weexplain how the discrete biharmonic operator ˜∆ x ,x is calculated in appendix A.Semi-implicit (or implicit-explicit) schemes for cloth animation have sometimes left bending forces explicit, whenthey are much smaller than stretching forces [41]. Considering a forward Euler discretization of (18) with bendingforces only, and approximating the bending force by B r n , stability requires ∆ t < C b µA p ∆ x /K b for a constant C b ,since the eigenvalues of B ∼ / ∆ x for lattice spacing ∆ x . The ratio of bending to stretching time step constraints ∼ K s ∆ x /K b ∼ ∆ x /h , the square of the ratio of lattice spacing to sheet thickness. For the parameters used in thiswork, the stretching time step constraint is typically smaller, but not much smaller than the bending time constraint,so a semi-implicit approach needs to address both terms.Our first-order semi-implicit discretization for sheet dynamics with both stretching and bending forces is: µA p r n +1 − r n ∆ t = K s L r n +1 + f SE ( r n ) + B r n +1 − B r n + f B ( r n ) . (32)The last two terms on the right hand side approximately cancel in the highest-derivative term, leaving the implicitbending term B r n +1 as the dominant one. The second-order version with uniform time stepping (an approximatebackward differentiation formula) is: µA p r n +1 − r n + r n − t = K s L r n +1 + 2 f SE ( r n ) − f SE ( r n − ) + B r n +1 − B r n + B r n − + 2 f B ( r n ) − f B ( r n − ) . (33)For test problems that model the deformation and dynamics of active gel sheets driven by internal swelling [12, 22,28, 30], we construct a uniform triangular lattice with spacing d . We then set the dilatation factor d ij /d ≡ η tocorrespond to one of three examples: a static radial distribution, a unidirectional traveling wave, or a radial travelingwave of isotropic dilation/contraction: η ( x , x , t ) = 1 + A cos (cid:18) π (cid:18) k (cid:113) x + x (cid:19)(cid:19) (34) η ( x , x , t ) = 1 + A sin (2 π ( kx − t )) (35) η ( x , x , t ) = 1 + A sin (cid:18) π (cid:18) k (cid:113) x + x − t (cid:19)(cid:19) . (36)The corresponding equilibrium metrics are¯ a αβ ( x , x , t ) = η ( x , x , t ) δ αβ , (37)and we take zero reference curvature (¯ b αβ = 0) for simplicity. Before presenting results, we describe in the nextsection a more direct finite difference simulation that allows for more general elastic parameters than the triangularlattice model. However, the triangular lattice simulations have the advantage of remaining bounded in time for alltime steps across wide ranges (several orders of magnitude) of values for the parameters ( h , A , k , µ ) with large meshsizes (hexagonal domains with N = 3367 and 9241, for example). Simulations appear smooth up to strain amplitudes A ∼ .
3. Above this value, sheet shapes become jagged and self-intersect, but remain bounded in time. Hence thetime step and lattice spacing are constrained only by the need to resolve the dynamics at a given parameter set. Themethod can also be used to converge to static equilibria, in which case the time step becomes a step length for agradient descent algorithm. Here a large step length may be used initially to rapidly approach the neighborhood of0an equilibrium, and then a smaller step length allows for convergence. The resulting convergence is geometric (notsuperlinear) but generally quite fast, and due to the simplicity of the formulation it is a good alternative to Newtonand quasi-Newton methods in the static case (as well as the dynamic case).
V. FINITE DIFFERENCE ALGORITHM
We now propose a second algorithm, inspired by that for the triangular lattice, but based on a finite differencediscretization of the continuum elastic energy (14), and which therefore allows the full range of values of E , h , and ν (that are physically reasonable). We use a square grid with grid spacing ∆ x for simplicity and define second-orderaccurate finite-difference operators: D α ≈ ∂ α , D αβ ≈ ∂ αβ , α, β = { , } (38)In the energy (14) we use a αβ ≈ D α r x (cid:12) D β r x + D α r y (cid:12) D β r y + D α r z (cid:12) D β r z , (39)where (cid:12) denotes a componentwise (Hadamard) product of two vectors. We use the analogous expression for b αβ anda trapezoidal-rule quadrature for the integrals in (14). For equilibrium metrics in the form (37) the discrete form ofthe stretching energy ( W s in (14)) is˜ W s = h A αβγδ (cid:0) q (cid:12) η (cid:1) (cid:124) [( a αβ − ¯ a αβ ) (cid:12) ( a γδ − ¯ a γδ )] (40)where q is the vector of quadrature weights for the trapezoidal rule on the rectangular mesh, and we use the usualsummation rule for repeated indices. We compute ∇ r x ˜ W s by using the chain rule with (39) and (40): ∇ r x ˜ W s = h D (cid:124) α (cid:2) A αβγδ (cid:0) q (cid:12) η (cid:12) ( a γδ − ¯ a γδ ) (cid:1) D β r x ) (cid:3) + h D (cid:124) β (cid:2) A αβγδ (cid:0) q (cid:12) η (cid:12) ( a γδ − ¯ a γδ ) (cid:1) D α r x ) (cid:3) (41)and the same expressions for ∇ r y ˜ W s and ∇ r z ˜ W s , with r x in (41) replaced by r y and r z respectively. Our linearizedapproximation to the stretching force, to compute r n +1 semi-implicitly, is M ns r n +1 , where M ns is a block diagonalmatrix with three blocks, each given by h D (cid:124) α (cid:2) A αβγδ (cid:0) q (cid:12) η (cid:12) a nγδ (cid:1) D β ) (cid:3) + h D (cid:124) β (cid:2) A αβγδ (cid:0) q (cid:12) η (cid:12) a nγδ (cid:1) D α ) (cid:3) . (42) M ns r n +1 is the discrete stretching force with zero reference metric, analogous to the discrete Laplacian on the triangularlattice, which gives the stretching force with zero-rest-length springs. The use of a nγδ instead of a n +1 γδ makes M ns independent of r n +1 . An extrapolation that is higher-order in time can also be used.We also compute ∇ r W b using the chain rule, resulting in a similar (though somewhat lengthier) expression for thebending force. Our linearized approximation to the bending force is M b r n +1 , where M b is the product of the bendingmodulus and the discrete biharmonic operator, the same as for the triangular lattice but now on a rectangular mesh.Unlike the triangular lattice algorithm, the finite difference algorithm is only stable for first-order time-stepping,i.e. µA p r n +1 − r n ∆ t = M ns r n +1 − M ns r n − (cid:110) ∇ r ˜ W s (cid:111) n + M b r n +1 − M b r n − (cid:110) ∇ r ˜ W b (cid:111) n , (43)1 n ∆ t SImax ∆ t FEmax
33 0.04 1.8 × −
44 0.03 1.0 × −
55 0.025 5 × −
66 0.02 3 × − TABLE I: Approximate upper bounds on stable time step for the semi-implicit finite difference scheme (SI) compared toforward Euler (FE), with different grid sizes n . The number of grid points is n and the number of unknowns is 3 n . the analogue of (32), with A p = ∆ x the area per point (multiplied by 1/2 at points along the sides and 1/4 at thecorners). It is unstable for the second-order version, µA p r n +1 − r n + r n − t = M [ n +1] s r n +1 − M [ n ] s r n + M [ n − s r n − − (cid:110) ∇ r ˜ W s (cid:111) n + (cid:110) ∇ r ˜ W s (cid:111) n − (44)+ M b r n +1 − M b r n + M b r n − − (cid:110) ∇ r ˜ W b (cid:111) n + (cid:110) ∇ r ˜ W b (cid:111) n − . (45)due to the extrapolated gradient terms (i.e. (44)–(45) becomes stable when the extrapolation reverts to first order forthe gradient terms in (45)). The superscripts in brackets denote a second-order extrapolation to the indicated timestep using the two preceding time steps. We note that second-order accuracy can be obtained from the first ordermethod via Richardson extrapolation.While the triangular lattice algorithm is essentially unconditionally stable, the first-order finite-difference algorithmis only stable up to moderately large time steps. However, the semi-implicit operators do yield orders-of-magnitudeimprovements in the largest stable time steps compared to an explicit scheme (forward Euler). As an example, weset η = η in (36), with A = 0 . h = 0 . µ = 1000, and ν = 1 /
3. In table I we give approximate values for themaximum stable time step, defined to be a time step such that the sheet deflection remains bounded (below 10 inmaximum norm) up to t = 5. The sheet is a square of side length 2, with n grid points in each direction. As n isincreased from 33 to 66, the maximum stable time step decreases by a factor of 2 for the semi-implicit method, versusa factor of 6 for the explicit method. Because the cost of solving the linear systems in (32), (33), and (43) is onlyslightly larger than the cost of the rest of the algorithm for the largest n in table I, the orders-of-magnitude differencein time step translates to an orders-of-magnitude difference in overall computational cost for these mesh sizes. VI. RESULTS
We now present a sequence of simulation results to display basic aspects of the algorithms and parameters. First,we compare the triangular lattice and finite difference methods in two situations. The first is the equilibrium sheetdeformation with a nontrivial, but static (time-independent) reference metric given by η = η in (34), with A = 0 . k = 1, h = 0 . µ = 1000, and ν = 1 /
3. The sheets are initially nearly flat squares ( z = 0 . x + x ), − ≤ x = x, x = y ≤ -0.301 z y x y -0.300.3 z -1 0 1 y -0.300.3 z -1 0 1 y -0.300.3 z -0.301 z y x z y x FIG. 3:
Static equilibria for the triangular lattice algorithm (left), and the finite difference algorithms with different Poisson ratios in oneof the bending energy terms (center and right). The reference metrics are given by η = η in (34), with physical parameters A = 0 . k = 1, h = 0 . µ = 1000, and ν = 1 / − / row) and side view (bottom row). The leftmost sheet is computed with the triangular lattice algorithm, with jaggededges along one pair of sides. The center sheet is the result of the finite difference algorithm, with a different Poissonratio value in one of the bending energy terms, to match those of the triangular lattice algorithm. For equilibriummetrics (37) we may write (12) as w b = Eh − ν ) η (cid:0) ( b + b ) − − ν )( b b − b ) (cid:1) . (46)In the stretching energy term (11) and the prefactor of the right hand side of (46), we set ν = 1 /
3, but we set ν to-1/3 in the second appearance of ν (within 1 − ν ) on the right hand side of (46). This gives an elastic energy consistentwith that of the triangular lattice model [63]. The rightmost sheet is also given by the finite difference algorithm butwith ν = 1 / z values, ∆ z , are within 1% for the center and leftmost sheets, and thedifference is about 2% for the center and rightmost sheets. In this case at least, the error in Gaussian bending rigidityin the triangular lattice algorithm has a modest effect, as was found by [64] in a different problem.We now consider an unsteady reference metric, η = η in (35), a unidirectional traveling wave, in both algorithms.We again take square sheets, with the same spatial grids and initial conditions as before, with (smaller) A = 0 . k = 1, h = 0 . µ = 1000, and ν = 1 / t = 0 . z versus time for the triangular lattice (black) andfinite difference (red) algorithms. The dynamics after buckling of the initial state are quite complex, and include theappearance of a higher buckling mode during 6 ≤ t ≤
8. For both algorithms the solutions have a strong oscillatorycomponent with the same frequency as the reference metric, as one would expect. The solutions have other componentsthat evolve on much longer time scales (tens to hundreds of periods), and even with modest deflections (here, about7% of the square side length), the solutions can evolve in complicated ways over long time scales. The two algorithmsshow rough qualitative agreement, though the discrepancy in deflection reaches 10–15% at certain times in panel A.3 -1 0 1 x -101 y B) -0.0500.050.1 -1 0 1 x -101 y C) -0.0500.050.1 Time z A) -1 0 1 x -101 y D) -0.0500.050.1 -1 0 1 x -101 y E) -0.0500.050.1 FIG. 4:
A comparison of sheet dynamics with the triangular lattice and finite difference algorithms with reference metric factor η = η in(35), a unidirectional traveling wave, starting from a nearly flat sheet. The physical parameters are A = 0 . k = 1, h = 0 . µ = 1000,and ν = 1 /
3. A) A comparison of the z -deflection over time for the triangular lattice (black) and finite difference (red) algorithms.Distributions of z deflection at time t = 29.5 and 30 are shown in panels B and C for the triangular lattice and D and E for the finitedifference method, respectively. The z deflection at t = 29.5 and 30 is shown in panels B and C for the triangular lattice and D and E for the finitedifference method, respectively. Both algorithms show a similar distribution of deflection on the right side of thesheets. On the left side, the triangular lattice has an additional peak (near x = − . , y = 0) that is not present inthe finite difference method.We have illustrated the behavior of the algorithms in two simple cases. We’ve seen close agreement in the finalequilibrium after buckling with a steady reference metric, and somewhat less agreement in a dynamical problem witha unidirectional traveling wave metric. The disagreements may be attributed to the sensitivity of the dynamicalproblem to different prefactors in the Gaussian curvature term of the bending energy and to the different spatialdiscretizations. We now proceed to illustrate some other basic features of the buckling behavior and dynamics, withthe triangular lattice approach only for brevity.If a static reference metric cannot be realized by a surface in R , the sheet has nonzero stretching energy, and ingeneral becomes unstable to out-of-plane buckling for sufficiently small h [66]. Intuitively, buckling allows a reductionin stretching energy at the expense of bending energy. The relative cost of bending decreases with decreasing h ,making buckling more favorable. For the static metric with η = η in (34), we plot a computational estimate of the“buckling threshold”—for various A ∈ [10 − , − ], the values of h at which buckling occurs. Here buckling is defined4 -6 -4 -2 -3 -2 -1 -10 -9 -10 -8 -6 -4 -2 FIG. 5:
Change in out-of-plane deflection from t = 0 to t = 2 for a unit hexagon with a small initial deflection z = 10 − x x , and imposedreference metric (34) with k = 1, varying metric factor amplitude A (horizontal axis) and varying sheet thickness h (vertical axis). Thedamping constant µ = 12 . d = 1/33. The solid red line shows the locus of zero change in deflection, essentially the“buckling threshold.” The red crosses show, at several values of h , the values of A on the buckling threshold with a finer mesh spacing,1/66. The white line shows the scaling A ∼ h . by whether a small initial deflection ( z = 10 − x x ) grows after two time units, with a certain damping constant( µ = 12 . d = 1 /
33. The mesh spacing becomes more important at smaller h , where bucklingdeformations may occur with a smaller wavelength due to the decreased bending energy. To check the effect of meshrefinement, we repeat the computations with d = 1 /
66 at selected points and obtain the red crosses. The white lineshows the scaling A ∼ h , which approximately matches the buckling threshold data presented in [66] for a differentreference metric. Our data appear to follow this trend for h ≥ .
01. Deviations at smaller h are likely due to thefinite mesh spacing.Next, we study the postbuckling behavior, again with the static metric factor η = η in (34) but with A = 0 . k = 1, and different initial perturbations of the form z = cr m sin( mφ ), where r is the initial distance of sheet pointsfrom the hexagon center, φ is the azimuthal angle (with respect to the hexagon center), and m is the azimuthalwavenumber, an integer from 2 to 6. The amplitude c ranges from 0 to 0.1, and µ ranges from 1 to 1000. We obtainbuckling into various equilibria akin to those studied by [66] with other radially symmetric reference metrics. Weshow six equilibria in figure 6, in two views (top and bottom rows). The first three (starting at the left) have a twofoldazimuthal rotational symmetry, but are distinct configurations. The second has an additional bilateral symmetry andthe third has an additional pair of local maxima along the sheet edges. The fourth and fifth have a threefold rotationalsymmetry, and the fourth has bilateral symmetries in addition. The sixth has a sixfold rotational symmetry. These5 FIG. 6:
A selection of six buckled equilibria of the unit hexagon. Each is shown in two views: the top row gives views from 60 degreeswith respect to the z -axis, and the bottom row gives views along the z -axis. The sheet thickness is h = 0 .
03 and the imposed referencemetric factor is η = η given in (34) with A = 0 . k = 1. The equilibria are obtained by starting from different initial perturbations,with lattice spacing d = 1/33. are only a selection of local equilibria for a given static reference metric, and illustrate the complexity of the energylandscape for such sheets, even without dynamics and/or a time-dependent reference metric. A. Parameter sweeps
We next study the effects of parameters on the sheet dynamics, using the reference metric factor η = η in (36),a radial traveling wave. We perform a sequence of parameter sweeps, varying one of A , h , or µ while keeping theothers fixed. We fix the wavenumber k at 1 in all cases. If k is much smaller, the sheet does not buckle becausethe reference metric is almost uniform, and planar dilation/contraction is preferred energetically. If k is much larger,buckling is also inhibited because the reference metric averaged over a local region approaches the identity tensor.As a base case, we take A = 0.1, µ = 1000, and h = 0.03, and vary each physical parameter in turn with numericalparameters d = 1 /
33 and ∆ t = 0 . z = 0 . x + x ), − ≤ x = x, x = y ≤ t = 19 to 20 in time increments of 0.2 (from left to right), forfour different values of µ , with A = 0.1 and h = 0.03. When µ is large, the sheet responds more slowly to the referencemetric. The tendency is to oscillate about a mean state of no deformation (because the long-time-average referencemetric is the identity tensor), and at sufficiently large µ the sheet does not buckle out of plane. At the largest µ (3000), the dynamics are essentially periodic (with period 2, so only a half-period is shown; the position at time t + 1is that at time t but reflected in the z = 0 plane). The out-of-plane deflection is smaller and more symmetric thatat smaller µ , where the dynamics become more chaotic and asymmetric. At smaller µ , the sheet moves more rapidlythrough the complicated, time-dependent elastic energy landscape. At µ = 1000, the sheet deflection ∆ z ( t ) has alarge component with the same period as the reference metric, but also large components that are not periodic. At µ = 300, bilateral symmetry is lost, while at µ = 100, the sheet almost assumes a threefold symmetry.Next, we vary h , with A = 0.1 and µ = 1000. Snapshots are shown in figure 8. With larger h , bending is relativelymore costly, so the deformation is smoother, and at sufficiently large h the sheet relaxes back to a planar state6 FIG. 7:
Snapshots of a unit hexagon, from t = 19 to 20 in time increments of 0.2 (from left to right) with different damping constants µ (top to bottom), metric factor amplitude A = 0.1 in (36) and h = 0.03. The colors show the z coordinate value, and the color scale isscaled to the minimum and maximum z coordinate value of each sheet. FIG. 8:
Snapshots of sheet dynamics with various sheet thicknesses h (top to bottom), from t = 19 to 20 in time increments of 0.2 (fromleft to right). Here A = 0.1 and µ = 1000. ( h = 0 . h , the sheet motion is more chaotic and asymmetric, as for decreasing µ . Notably, at smaller h fine wrinkling features appear, and the deformation is far from any kind of symmetry. Eventually the sheet mayintersect itself (not shown) because forces due to self-contact are not included.We vary A next, with h = 0 .
03 and µ = 1000. Here (figure 9) there is a sequence of dynamics from slow relaxationto planar motions ( A = 0.03), to periodic and symmetric ( A = 0.06, with period 2), to aperiodic but still bilaterallysymmetric ( A = 0.1 and 0.2). At larger A the deflection amplitude increases but the deformation remains relativelysmooth, and in this sense the dynamics are more similar to smaller µ than to smaller h .Finally, we consider the effect of the sheet shape on the dynamics. We consider the two cases with periodic dynamicsabove, ( A = 0 . µ = 3000) in figure 7 and ( A = 0 . µ = 1000) in figure 9, both with h = 0 .
03. In the top tworows of figure 10, we compare the hexagonal sheet snapshots with those of a square sheet with the same parameters.7
FIG. 9:
Snapshots of sheet dynamics with various metric factor amplitudes A (in (36)), from t = 19 to 20 in time increments of 0.2 (fromleft to right). Here h = 0.03 and µ = 1000. FIG. 10:
Comparisons of the dynamics of hexagonal and square sheets with A = 0 . µ = 3000 in the top two rows and A = 0 .
06 and µ = 1000 in the bottom two rows. All other parameters and initial conditions are the same as in figures 7 and 9 (e.g. h = 0 . t = 19 to 20 in time increments of 0.2 (from left to right). Both sheets have width 2 when the reference metric is the identity. The bottom two rows make the same comparisonwith A = 0 .
06 and µ = 1000. At each time, the hexagonal and square sheets have qualitative similarities. At thefirst time (leftmost snapshots), the sheets have a central hump. At the second time, an upward ring (yellow) appears.Third, the ring breaks up into an array of humps with sixfold symmetry for the hexagons, and either six- or fourfoldsymmetry for the squares. Fourth, the humps reach the sheet boundary. Fifth, a central depression forms, and sixth,it widens and is surrounded by an upward ring. The main differences between the squares and hexagons are that thesquares more often have a fourfold rather than sixfold symmetry, the maximum deflections of the squares are larger,and the squares have a significant nonperiodic component in the dynamics. The comparison illustrates that periodicdynamics can be somewhat sensitive to sheet shape, which is perhaps not surprising given that chaotic dynamics arecommon and small changes in various parameters can shift periodic dynamics to chaotic dynamics.8 VII. CONCLUSION
We have presented semi-implicit algorithms that can simulate the dynamics of thin elastic sheets with large timesteps. We focus on the case of elastic sheets with nontrivial steady or time-varying reference metrics in overdampeddynamics, but the methods should apply to other problems with internal or external forcing. The first algorithmsimulates a triangular lattice mesh, and uses a splitting of the stretching force that has been used previously forcomputer graphics simulations of hair and cloth [40, 43]. The semi-implicit bending force uses a biharmonic operatorwith free-edge boundary conditions. The algorithm appears to be unconditionally stable for a periodic referencemetric. The second algorithm simulates a rectangular grid with finite-difference derivatives of the energy, and allowsfor general values of the Poisson ratio (and unlike the triangular lattice, a consistent value in all terms of the bendingand stretching energies). The semi-implicit finite difference algorithm is analogous to that for the triangular lattice,involving a stretching force operator for zero-rest-length springs, and the biharmonic operator for an approximatebending force. The finite difference algorithm is not unconditionally stable, but typically has a maximum stable timestep two to three orders of magnitude greater than that of an explicit scheme.The two algorithms agree very closely for the deformation of a square sheet under a static reference metric, evenallowing a different Poisson ratio ( − / Acknowledgments
We acknowledge support from the Michigan Institute for Computational Discovery and Engineering (MICDE).
Appendix A: Bending force
We denote the discrete biharmonic operator on the triangular mesh with free-edge boundary conditions as ˜∆ x ,x ,and because it is independent of the sheet configuration we may derive it for a sheet that is a portion of a flattriangular lattice with edge length d . ˜∆ x ,x is a mapping from (small) out-of-plane displacements at each point tothe bending force at that point. The mapping is a sum of the force-displacement mappings for each pair of adjacent9equilateral triangles, because the bending energy is also a sum over such units. For the four vertices in a given pairof neighboring triangles, a unit upward displacement of one of the two outer vertices (those not on the shared edge)yields a downward bending force at each of the outer vertices, equal by symmetry, and an upward force at each of theinner two vertices (those on the shared edge), equal and opposite to those at the outer vertices, to main net force andtorque balance. The forces are proportional to the displacement by a constant that gives the desired bending modulusof the sheet. The bottom row of figure 2 shows examples of the stencils for this biharmonic operator at a few meshpoints: an interior point (analogous to the 13-point biharmonic stencil on a rectangular mesh), a next-to-boundarypoint, and a boundary point. The discrete biharmonic operator ˜∆ x ,x has these stencil values multiplied by thetriangle altitudes raised to the − (cid:0) √ d/ (cid:1) − , while the linearized bending force operator has instead thesame stencil values multiplied by (cid:0) √ d/ (cid:1) − because the bending force is the gradient of the bending energy, whichis approximately the biharmonic of the deflection multiplied by the sheet area per vertex. [1] C. R. Calladine. Theory of shell structures . Cambridge University Press, 1983.[2] L. B. Freund and S. Suresh.
Thin film materials: stress, defect formation and surface evolution . Cambridge UniversityPress, Cambridge, 2004.[3] Shahaf Armon, Efi Efrati, Raz Kupferman, and Eran Sharon. Geometry and Mechanics in the Opening of Chiral SeedPods.
Science , 333(6050):1726–1730, 2011.[4] B Audoly and A Boudaoud. Self-similar structures near boundaries in strained systems.
Physical Review Letters , 91(8),2003.[5] C. D. Santangelo. Buckling thin disks and ribbons with non-Euclidean metrics.
Europhysics Letters , 86(3), 2009.[6] John A. Gemmer and Shankar C. Venkataramani. Shape selection in non-Euclidean plates.
Physica D , 240(19):1536–1552,2011.[7] Eran Sharon. Swell Approaches for Changing Polymer Shapes.
Science , 335(6073):1179–1180, 2012.[8] Eran Sharon and Hillel Aharoni. Frustrated shapes.
Nature Materials , 15(7):707–709, 2016.[9] E. Efrati, E. Sharon, and R. Kupferman. Elastic theory of unconstrained non-Euclidean plates.
Journal of the Mechanicsand Physics of Solids , 57(4):762–775, 2009.[10] Efi Efrati, Eran Sharon, and Raz Kupferman. The metric description of elasticity in residually stressed soft materials.
SoftMatter , 9(34):8187–8197, 2013.[11] Jun-Hee Na, Arthur A. Evans, Jinhye Bae, Maria C. Chiappelli, Christian D. Santangelo, Robert J. Lang, Thomas C. Hull,and Ryan C. Hayward. Programming Reversibly Self-Folding Origami with Micropatterned Photo-Crosslinkable PolymerTrilayers.
Advanced Materials , 27(1):79–85, 2015.[12] Yael Klein, Efi Efrati, and Eran Sharon. Shaping of elastic sheets by prescription of non-Euclidean metrics.
Science ,315(5815):1116–1120, 2007.[13] Jungwook Kim, James A. Hanna, Myunghwan Byun, Christian D. Santangelo, and Ryan C. Hayward. Designing ResponsiveBuckled Surfaces by Halftone Gel Lithography.
Science , 335(6073):1201–1205, 2012.[14] Zi Liang Wu, Michael Moshe, Jesse Greener, Heloise Therien-Aubin, Zhihong Nie, Eran Sharon, and Eugenia Kumacheva.Three-dimensional shape transformations of hydrogel sheets induced by small-scale modulation of internal stresses.
Nature Communications , 4:1586, 2013.[15] T. H. Ware, M. E. McConney, J. J. Wie, V. P. Tondiglia, and T. J. White. Voxelated liquid crystal elastomers.
Science ,347(6225):982–984, 2015.[16] Timothy J. White and Dirk J. Broer. Programmable and adaptive mechanics with liquid crystal polymer networks andelastomers.
Nature Materials , 14(11):1087–1098, 2015.[17] A. Sydney Gladman, Elisabetta A. Matsumoto, Ralph G. Nuzzo, L. Mahadevan, and Jennifer A. Lewis. Biomimetic 4Dprinting.
Nature Materials , 15(4):413, 2016.[18] Yen-Chih Lin, Ji-Ming Sun, Jen-Hao Hsiao, Yeukuang Hwu, C. L. Wang, and Tzay-Ming Hong. Spontaneous emergenceof ordered phases in crumpled sheets.
Phys. Rev. Lett. , 103:263902, 2009.[19] B Roman and J Bico. Elasto-capillarity: deforming an elastic structure with a liquid droplet.
Journal of Physics: CondensedMatter , 22(49):493101, 2010.[20] Hillel Aharoni and Eran Sharon. Direct observation of the temporal and spatial dynamics during crumpling.
NatureMaterials , 9(12):993–997, 2010.[21] Joseph D. Paulsen, Vincent Demery, K. Bugra Toga, Zhanlong Qiu, Thomas P. Russell, Benny Davidovitch, and NarayananMenon. Geometry-Driven Folding of a Floating Annular Sheet.
Physical Review Letters , 118(4), 27 2017.[22] R Yoshida, T Takahashi, T Yamaguchi, and H Ichijo. Self-oscillating gel.
Journal of the American Chemical Society ,118(21):5134–5135, 1996.[23] Shingo Maeda, Yusuke Hara, Takamasa Sakai, Ryo Yoshida, and Shuji Hashimoto. Self-walking gel.
Advanced Materials ,19(21):3480, 2007.[24] O Tabata, H Hirasawa, S Aoki, R Yoshida, and E Kokufuta. Ciliary motion actuator using self-oscillating gel.
Sensorsand Actuators A-Physical , 95(2-3):234–238, 2002.[25] O Tabata, H Kojima, T Kasatani, Y Isono, and R Yoshida. Chemo-mechanical actuator using self-oscillating gel forartificial cilia. In
MEMS-03: IEEE The Sixteenth Annual International Conference on Micro Electro Mechanical Systems ,Proceedings: IEEE Micro Electro Mechanical Systems, pages 12–15. IEEE, 2003.[26] Shingo Maeda, Yusuke Hara, Ryo Yoshida, and Shuji Hashimoto. Peristaltic motion of polymer gels.
Angewandte Chemie-International Edition , 47(35):6690–6693, 2008.[27] Yusuke Shiraki and Ryo Yoshida. Autonomous Intestine-Like Motion of Tubular Self-Oscillating Gel.
Angewandte Chemie-International Edition , 51(25):6112–6116, 2012.[28] VV Yashin and AC Balazs. Modeling polymer gels exhibiting self-oscillations due to the Belousov-Zhabotinsky reaction.
Macromolecules , 39(6):2024–2026, 2006.[29] Olga Kuksenok, Victor V. Yashin, and Anna C. Balazs. Mechanically induced chemical oscillations and motion in responsivegels.
Soft Matter , 3(9):1138–1144, 2007.[30] Victor V Yashin and Anna C Balazs. Theoretical and computational modeling of self-oscillating polymer gels.
The Journalof chemical physics , 126(12):124707, 2007.[31] Olga Kuksenok, Debabrata Deb, Pratyush Dayal, and Anna C. Balazs. Modeling chemoresponsive polymer gels.
AnnualReview of Chemical and Biomolecular Engineering , 5(1):35–54, 2014. PMID: 24498954.[32] Mila Boncheva, Stefan A Andreev, L Mahadevan, Adam Winkleman, David R Reichman, Mara G Prentiss, Sue Whitesides,and George M Whitesides. Magnetic self-assembly of three-dimensional surfaces from planar sheets.
Proceedings of theNational Academy of Sciences , 102(11):3924–3929, 2005.[33] Silas Alben and Michael P Brenner. Self-assembly of flat sheets into closed surfaces.
Physical Review E , 75(5):056113, NanoLetters , 11(6):2280–2285, 2011.[35] Silas Alben. Bending of bilayers with general initial shapes.
Advances in Computational Mathematics , 41(1):1–22, 2015.[36] Shravan K Veerapaneni, Abtin Rahimian, George Biros, and Denis Zorin. A fast algorithm for simulating vesicle flows inthree dimensions.
Journal of Computational Physics , 230(14):5610–5634, 2011.[37] HH Rosenbrock. Some general implicit processes for the numerical solution of differential equations.
The Computer Journal ,5(4):329–330, 1963.[38] David Gottlieb and Steven A Orszag.
Numerical analysis of spectral methods: theory and applications , volume 26. Siam,1977.[39] Thomas Y Hou, John S Lowengrub, and Michael J Shelley. Removing the stiffness from interfacial flows with surfacetension.
Journal of Computational Physics , 114(2):312–338, 1994.[40] Mathieu Desbrun, Peter Schr¨oder, and Alan Barr. Interactive animation of structured deformable objects. In
GraphicsInterface , volume 99, page 10, 1999.[41] Bernhard Eberhardt, Olaf Etzmuß, and Michael Hauth. Implicit-explicit schemes for fast animation with particle systems.In
Computer Animation and Simulation 2000 , pages 137–151. Springer, 2000.[42] Thomas Y Hou, John S Lowengrub, and Michael J Shelley. Boundary integral methods for multicomponent fluids andmultiphase materials.
Journal of Computational Physics , 169(2):302–362, 2001.[43] Andrew Selle, Michael Lentine, and Ronald Fedkiw. A mass spring model for hair simulation. In
ACM Transactions onGraphics (TOG) , volume 27, page 64. ACM, 2008.[44] Silas Alben. An implicit method for coupled flow–body dynamics.
Journal of Computational Physics , 227(10):4912–4933,2008.[45] Silas Alben. Simulating the dynamics of flexible bodies and vortex sheets.
Journal of Computational Physics , 228(7):2587–2603, 2009.[46] Hsiao-Yu Chen, Arnav Sastry, Wim M van Rees, and Etienne Vouga. Physical simulation of environmentally induced thinshell deformation.
ACM Transactions on Graphics (TOG) , 37(4):146, 2018.[47] HS Seung and DR Nelson. Defects in flexible membranes with crystalline order.
Phys. Rev. A , 38(2):1005 – 1018, 1988.[48] Warner Tjardus Koiter. On the nonlinear theory of thin elastic shells.
Proc. Koninkl. Ned. Akad. van Wetenschappen,Series B , 69:1–54, 1966.[49] Philippe G Ciarlet. Un mod`ele bi-dimensionnel non lin´eaire de coque analogue `a celui de wt koiter.
Comptes Rendus del’Acad´emie des Sciences-Series I-Mathematics , 331(5):405–410, 2000.[50] Roman Vetter, Norbert Stoop, Thomas Jenni, Falk K Wittel, and Hans J Herrmann. Subdivision shell elements withanisotropic growth.
International Journal for Numerical Methods in Engineering , 95(9):791–810, 2013.[51] Efi Efrati, Yael Klein, Hillel Aharoni, and Eran Sharon. Spontaneous buckling of elastic sheets with a prescribed non-Euclidean metric.
Physica D , 235(1-2):29–32, 2007.[52] Jorge Nocedal and Stephen Wright.
Numerical optimization . Springer Science & Business Media, 2006.[53] Ernst Hairer, Syvert P Norsett, and Gerhard Wanner.
Solving Ordinary Differential Equations II: Stiff Problems . Springer,1987.[54] Gerhard Gompper and Daniel M Kroll. Fluctuations of polymerized, fluid and hexatic membranes: continuum models andsimulations.
Current opinion in colloid & interface science , 2(4):373–381, 1997. [55] Alexander E Lobkovsky and TA Witten. Properties of ridges in elastic membranes. Physical Review E , 55(2):1577, 1997.[56] Jack Lidmar, Leonid Mirny, and David R Nelson. Virus shapes and buckling transitions in spherical shells.
Physical ReviewE , 68(5):051910, 2003.[57] GA Vliegenthart and G Gompper. Forced crumpling of self-avoiding elastic sheets.
Nature materials , 5(3):216, 2006.[58] Eleni Katifori, Silas Alben, Enrique Cerda, David R Nelson, and Jacques Dumais. Foldable structures and the naturaldesign of pollen grains.
Proceedings of the National Academy of Sciences , 107(17):7635–7639, 2010.[59] Etienne Couturier, J Dumais, E Cerda, and Elenei Katifori. Folding of an opened spherical shell.
Soft Matter , 9(34):8359–8367, 2013.[60] Chloe M Funkhouser, Rastko Sknepnek, and Monica Olvera De La Cruz. Topological defects in the buckling of elasticmembranes.
Soft Matter , 9(1):60–68, 2013.[61] Duanduan Wan, David R Nelson, and Mark J Bowick. Thermal stiffening of clamped elastic ribbons.
Physical Review B ,96(1):014106, 2017.[62] David R Nelson, Tsvi Piran, and Steven Weinberg.
Statistical mechanics of membranes and surfaces . World Scientific,2004.[63] Bernd Schmidt and Fernando Fraternali. Universal formulae for the limiting elastic energy of membrane networks.
Journalof the Mechanics and Physics of Solids , 60(1):172–180, 2012.[64] Brian Anthony DiDonna. Scaling of the buckling transition of ridges in thin sheets.
Physical Review E , 66(1):016601, 2002.[65] Joel L Weiner. On a problem of Chen, Willmore, et al.
Indiana University Mathematics Journal , 27(1):19–35, 1978.[66] Efi Efrati, Eran Sharon, and Raz Kupferman. Buckling transition and boundary layer in non-Euclidean plates.