Simulating Thin Sheets: Buckling, Wrinkling, Folding and Growth

Numerical simulations of thin sheets undergoing large deformations are computationally challenging. Depending on the scenario, they may spontaneously buckle, wrinkle, fold, or crumple. Nature's thin tissues often experience significant anisotropic growth, which can act as the driving force for such instabilities. We use a recently developed finite element model to simulate the rich variety of nonlinear responses of Kirchhoff-Love sheets. The model uses subdivision surface shape functions in order to guarantee convergence of the method, and to allow a finite element description of anisotropically growing sheets in the classical Rayleigh-Ritz formalism. We illustrate the great potential in this approach by simulating the inflation of airbags, the buckling of a stretched cylinder, as well as the formation and scaling of wrinkles at free boundaries of growing sheets. Finally, we compare the folding of spatially confined sheets subject to growth and shrinking confinement to find that the two processes are equivalent.


Introduction
Thin sheets are omnipresent in nature, technology and everyday life, appearing at virtually all length scales. Being much thinner in one than in the other two dimensions, they can develop an unparalleled, rich variety of deformation modes when subjected to external forces, spatial constraints, or intrinsic growth. They buckle, wrinkle, fold, and crumple. Numerical simulations are an often-indispensable approach for studying the complex interplay of these modes. The folding and crumpling of a piece of paper [1][2][3][4][5][6][7] and metal sheets wrinkling and crumpling in vehicle collisions [8][9][10][11] are two out of many examples. For many such problems, the finite element method (FEM) has shown to be amongst the most efficient and flexible tools, especially in cases with strong material nonlinearity, complex geometry, or anisotropy. Even though the Kirchhoff-Love theory [12] provides a simple kinematic description, numerically sound finite element implementations of thin sheets have turned out to be difficult and cumbersome in the past. These problems can be successfully overcome since the subdivision surface paradigm was introduced to the FEM [13,14].
Large deformations of soft thin tissue such as insect wings, plant leaves, cell membranes, or flowers are often induced by growth (or shrinkage) [15,16], inevitably leading to the development of residual stress [17,18]. In this paper, we present an extension of the Kirchhoff-Love theory to allow for anisotropic inplane growth, which we implement with Loop subdivision shape functions. The combination of these two concepts grants access to a very straightforward and highly efficient, yet powerful and flexible numerical tool for the simulation of nonlinear thin sheet mechanics. Our approach accounts for the change of reference curvature when the surface grows, generalizing recently developed tethered mass-spring models [19,20]. The next section summarizes the mentioned extension. It is followed in the subsequent sections by a series of thin sheet problems that we solve using the developed FEM implementation. Special attention is paid to the formation of self-similar wrinkles along a plastically stretched free edge as well as on the scaling of single-wavelength wrinkles of growing cylinders similar to flowers.

The Kirchhoff-Love Sheet with Anisotropic Growth
Let Ω ⊂ E 3 be the stress-free undeformed ("reference") middle surface of a sheet with small thickness h. Under the action of external forces or growth, the sheet deforms into a new configuration with middle surface Ω ⊂ E 3 . In the following, let Greek indices α, β, γ, δ ∈ {1, 2}, and Latin indices i, j ∈ {1, 2, 3}. Lower (upper) indices will denote covariant (contravariant) components. Moreover, let {θ 1 , θ 2 , θ 3 } be a curvilinear coordinate system, and let x(θ 1 , θ 2 ) and x(θ 1 , θ 2 ) be parametrizations of Ω and Ω, respectively (see Fig. 1). The material points p and p = χ(p) in the reference and deformed sheet are parametrized as where θ 3 ∈ [−h/2, h/2]. χ is a diffeomorphism that maps from the reference to the deformed material positions. The tangent spaces of Ω and Ω are spanned by the respective vector fields By virtue of the Kirchhoff assumption, straight material lines normal to the middle surface retain these properties as well as their length. They are determined by the unit normal vectors Figure 1: Reference, grown and deformed configurations of the sheet's middle surface.
The covariant components of the first fundamental forms follow as a αβ = a α · a β and a αβ = a α · a β , while those of the second fundamental forms are given by b αβ = a 3 · a α,β and b αβ = a 3 · a α,β .
Assuming that the thin sheet obeys the St. Venant-Kirchhoff law of linear elasticity, the connection between its kinematics and energetics is provided by the Koiter energy density functional [21,22]. Let the sheet be characterized by Young's modulus E and Poisson's ratio ν. The elastic energy U e of the Koiter sheet is obtained by integrating the energy per unit area over the middle surface: where dΩ = |a 1 × a 2 | dθ 1 dθ 2 . The Einstein summation applies to repeated indices. H is often referred to as the "elastic tensor", and is given component-wise by α = (a − a)/2 and β = b − b are the in-plane (2×2) membrane and bending strain tensors, respectively. The Koiter shell (6) can be extended to incorporate anisotropic growth through the multiplicative decomposition of the geometric deformation gradient ∇χ = F e F g [23,24] into a growth tensor F g and a purely elastic response F e , that ensures continuity and compatibility of the body. Owing to the Kirchhoff constraints, we may write and the growth-modified elastic strains then simply read [25] We further augment the elastic energy (6) with an inertial term to capture the dynamics of the thin sheet. The kinetic energy reads where ρ is the mass density of the sheet, andẋ = ∂x/∂t is the velocity. Our aim is to find the minimizer x of the total energy U = U e + U k for given growth tensors F g or external driving forces.

Finite Element Implementation with Subdivision Surfaces
To account for out-of-plane bending rigidity, the bending term in Eq. (6) integrates the Gaussian and mean curvatures, which comprise second derivatives of the displacement field u = x − x, over the middle surface. For boundedness of the integral in the weak formulation, continuously differentiable finite element shape functions (C 1 -continuity) are needed. This requirement has proven very challenging in the history of shell finite elements, until Cirak et al. [13,14] have introduced Loop subdivision surfaces to the FEM. A fundamental difference to traditional finite elements is that subdivision surfaces gain C 1 -continuity at the expense of a larger support of the individual shape functions. Details on their implementation are given in Refs. [13,25]. Aside from guaranteeing convergence, subdivision surfaces allow a classical Rayleigh-Ritz formulation of the sought finite element deformation: no rotational variables are needed and the only unknowns are the three nodal displacements. Moreover, a single quadrature point per element is sufficient [13,14,25], rendering this finite element approach computationally highly efficient and flexible. Of course, increasing the number of quadrature points may assist in resolving strongly anisotropic growth fields or constitutive relations. We employ Loop subdivision surface shape functions here to minimize the total energy U numerically by solving Newton's equation of motion with a standard predictor-corrector scheme. Subcritical viscous damping is added for numerical stability and equilibration.

Inflated Pillows and Airbags
As a first instance of folding and wrinkling, we reproduce the inflation of pillows and airbags from Ref. [14]. A square sheet with diagonal length d = 120 cm and thickness h = 1 mm and a circular sheet with radius R = 35 cm and thickness h = 0.4 mm are instantaneously pressurized with 5 kPa to buckle out of their flat initial configuration. The elastic moduli are given by E = 588 MPa, ν = 0.4 and E = 60 MPa, ν = 0.3, respectively. We exploit the reflection symmetry by simulating only the upper half of the geometry. The quasi-static equilibrium configurations are shown in Fig. 2, and a movie of the dynamic deformation is provided in the supplementary material. The square pillow features a distinct folding pattern in the middle of the four edges, resulting from the non-uniform distribution of Gaussian curvature: The edges are closer to the point of maximal uplift than the corners, thus getting curved more and pulled towards the center. They are laterally stretched and longitudinally compressed, yielding the observed folds. The circular airbag, on the other hand, is initially axisymmetric and therefore behaves differently. Wrinkles develop similarly to elastic plates stamped into curved cavities [26] or ultrathin films placed on fluid drops [27]. On top of these low-amplitude wrinkles, axisymmetry of the stress field is broken and large crumples occur that localize the geometrically imposed Gaussian curvature.

Buckling of a Stretched Cylinder
A frequently studied buckling problem is the laterally stretched open cylinder [28,29]. Two opposite point forces of equal increasing magnitude F cause the cylindrical sheet with radius R = 4.953 cm, length L = 10.35 cm, thickness h = 0.94 mm, and free edges to first bend before snapping through at F c = 11.836 Eh 3 /R to the post-buckling regime, where further deformations are dominated by stretching. The elastic moduli are fixed to E = 10.5 MPa, ν = 5/16. In Fig. 3, we plot the radial displacements of the points A, B and C, which are degenerate at the buckling threshold F c . The corresponding movie can be found in the supplementary material.

Boundary Instabilities and Wrinkling
An interesting feature observed in plant growth is the occurrence of self-similar wrinkles along free tissue boundaries, such as the edges of flowers and leaves [19,20,[30][31][32][33]. The morphology is very similar to the shape of torn plastic sheets, and apparently, both phenomena are characterized by a plastic longitudinal metric profile g l (z) = 1/(1 + z/l), where l > 0 is a characteristic length scale and z ≥ 0 is the coordinate perpendicular to the growing edge. Audoly and Boudaud [34] were able to show that the solution of the Föppl-von Kármán equations on the edge of a free rectangular sheet with such growth profiles consists of self-similar wrinkles governed by odd integral scaling factors. The Föppl-von Kármán equations are, however, geometrically limited as they don't allow reentrancy. The present growing Koiter shell model allows us to numerically solve the problem with its full geometric nonlinearity taken into account. Consider a flat rectangular sheet of thickness h = 10 −4 , length L = 4, and width W = 1, initially lying in the xz plane. We clamp the long edge at z = W and constrain the short edges to stay at x = −L/2, L/2, leaving them free to move in other directions. Plastic growth is imposed by setting the growth tensor to F g = diag 1 + g l (z), 1, 1 (12) in Cartesian coordinates (x, y, z), and we choose a characteristic length l = 40h for the growth field in this example. Fig. 4 shows the resulting equilibrated configuration after growth (or tearing), and a movie showing the equilibration is provided in the supplementary material. The self-similarity of the free edge at z = 0 is apparent, clearly resembling the wrinkling cascades observed in experiments [20,30]. We have measured the fractal dimension of the grown edge depicted in Fig. 4(a) using the box counting method [35] and the self-similarity method [36]. In the former, the curve length L in contained in a cubic box is determined as a function of the edge length L box of the box. A fractal curve is expected to scale as L in ∼ L D f box . Such scaling is indeed observed with fractal dimension D f = 1.15(1) (Fig. 5(a)). The scaling breaks down due to the influence of the clamped opposite edge, which introduces a global straightening effect when the box size is large (L box ≈ W ). The second method is more robust to global orientation and is thus better suited here. The length L s of a piecewise linear path along the curve with segment size s is measured and expected to scale as L s ∼ s 1−D . We find a self-similarity dimension D = 1.196(5) (Fig. 5(b)), which is very close to the Hausdorff dimension of a 60-degree Fibonacci word fractal ( Fig. 4(b)), D H = 1.2083 [37], and a bit lower than that of a triadic Koch curve (Fig. 4(c)), D H = 1.2619 [38].  The metric profile used above is not the only one yielding wrinkled edges. When it comes to wavy flowers like certain orchids for instance, single-wavelength undulations instead of self-similar edges are not uncommon. The feature causing wrinkle cascades is the presence of a non-constant geometric length scale defined by l geo (z) = −g(z)/g (z) [33]. The following families of growth fields will thus produce similar boundary instabilities: On the other hand, an exponential growth field yields only a single wavelength [19] because l geo (z) ≡ l in this case. On some flowers, these undulations may be forced to integral wavenumbers n by angular periodicity of a single petal. Let's hence significantly increase the characteristic length l and thickness h such that only a single wavelength λ prevails even for growth in the form of Eqs. (13,14), and let's consider a cylinder with height H and radius R instead of a flat plate. This change in geometry delivers dramatic consequences: A thin cylindrical sheet growing in circumferential direction according to g(z), where z is the cylinder axis, only breaks its axisymmetry if growth leads to a circumference that changes faster than the sheet's metric can account for [31]. (Note that the excluded linear case p = 1 of Eq. (14) produces the excess cone [16,39] in the limit R → 0, which doesn't wrinkle at the boundary because g l,1 (z) ≡ 0.) A direct consequence of the Gauss-Bonnet theorem is that the axisymmetry is preserved as long as and broken otherwise. The origin for this instability is the (non)existence of embeddings of the surface: According to Gauss's Theorema Egregium, the creation or elimination of Gaussian curvature must be accompanied by in-plane stretching, which is traded for out-of-plane buckling if the sheet is sufficiently thin (see Fig. 6). (c) A relatively short cylinder (small H/l) with free boundaries also buckles away from the wrinkled edge, breaking C n symmetry further to C 2 .
How does the number n of boundary waves scale with l, h, R when axisymmetry is broken, and is it universal for all positive, monotonically decreasing and strictly convex growth profiles satisfying lim z→0 l geo (z) = l ?
Since the preferred wavelength λ is a local feature independent of global geometry and topology (independent of R), we may use the ansatz [33] λ ∼ h α l 1−α geo , i.e., On the other hand, geometry implies that λ = 2πR n (1 + g(0)) . (19) since z = 0 is where Ineq. (16) is violated first, given that g < 0 and g > 0. After combining Eqs. (16)(17)(18)(19), one thus finds a scaling for the number of wrinkles Up to the first term, which accounts for the mean curvature of the cylinder, this coincides with the scaling law reported in [34] for the wrinkling hierarchies in initially flat sheets, where α = 2/5 is found for the family g l (z) ∝ 1/(1 + z/l). Our numerical data, which best fits Eq. (20) with α = 0.39 (2), shows that the wrinkles studied here fall in the same category, see Fig. 7. Moreover, the data collapse of all employed growth profiles on a single line indicates that the scaling is universal in this respect. For this numerical study, we used F g = diag 1, 1 + g(z), 1 in cylindrical coordinates (r, ϕ, z) with various growth profiles proportional to Eqs. (13)(14)(15), and with slowly increasing proportionality prefactors. A movie showing such spontaneous wrinkling for different n is provided in the supplementary material.  Figure 7: Scaling of the wrinkle number n along the edge of a circumferentially growing cylinder. The mode with the lowest energy is the integer n closest to the power law (20), but excited modes also randomly occur and are metastable (data points lying significantly above or below the straight line).

Confined Growth and Crumpling
Many numerical simulations of thin sheets getting folded and crumpled inside of shrinking hollow spheres have been carried out recently [4][5][6]40]. An important result is that the high bulk stiffness of crumpled sheets is due to a network of vertices and lines carrying large mean curvature and bending energy [1]. But what if instead of externally forced compression, the thin sheet intrinsically grows inside of a fixed spatial container? Are the processes that crumple a plant leaf or petal growing inside a bud the same as for a piece of foil crumpled by hand? Indeed they are in the elastic limit, as the following simulation demonstrates.
A thin circular sheet (radius R) is placed inside of a spherical cavity (radius R). In the first setup, the container is shrunk, folding and crumpling the sheet into a ball of the size of the container. In the second setup, the container sustains its size while the sheet undergoes uniform isotropic growth, both in plane and in thickness. The only important parameter for this problem is the Föppl-von Kármán number γ ∼ (R/h) 2 = 10 4 . Equivalent time scales are obtained by shrinking the sphere according to R(t) = R/(1+g(t)), where g(t) = λt is the growth factor of the growing sheet. The growth rate λ is chosen small enough to keep inertial effects negligible. We add repulsive contact forces penalizing volumetric overlap between any two pieces of the sheet. Initially, both sheets buckle to form a developable cone with a single vertex (see Fig. 8, first column) that starts to nucleate at R/R ≈ 0.53. The emerging ridge network of focused mean curvature (third column) and bending energy (fourth column) is the same in both scenarios. Compared to similar measurements [5], the cross correlation r = 0.89 of the mean curvature ridge patterns is very high. The differences are only local and of the order of the mesh resolution. A movie showing the crumpling process is contained in the supplementary material.

Conclusions
Simulating the plastic growth of thin sheets can be numerically demanding. The finite element method, being perhaps the best-suited technique for anisotropies, lacked an expedient C 1 -continuous discretization, which is indispensable from a theoretical viewpoint, until subdivision surfaces were ported to it. We have extended the Kirchhoff-Love theory by arbitrary volumetric in-plane growth [25] and used the Loop subdivision surface paradigm to build a highly flexible and efficient numerical tool for the simulation of nonlinear thin sheet mechanics. Requiring no rotational variables, a thin sheet representation of this kind is remarkably simple to implement and superior to traditional approaches in terms of computational costs. A series of example simulations were carried out to demonstrate these strengths. In particular, we have quantified the self-similarity and scaling of wrinkles along plastically stretched free edges of thin sheets, and found that the problem of a growing confined sheet is equivalent to the narrowing confinement of a static sheet.