Geometrically incompatible confinement of solids

# Geometrically incompatible confinement of solids

Benny Davidovitch Department of Physics, University of Massachusetts, Amherst, MA 01003    Yiwei Sun Department of Physics, University of Massachusetts, Amherst, MA 01003    Gregory M. Grason Department of Polymer Science and Engineering, University of Massachusetts, Amherst, MA 01003
July 19, 2019
###### Abstract

The complex morphologies exhibited by spatially confined thin objects have long challenged human efforts to understand and manipulate them, from the representation of patterns in draped fabric in Renaissance art to current day efforts to engineer flexible sensors that conform to the human body. We introduce a mathematical principle, broadly generalizing Euler’s elastica – a core concept of continuum mechanics that invokes the energetic preference of bending over straining a thin solid object and has been widely applied to classical and modern studies of beams and rods. We define a class of geometrically incompatible confinement problems, whereby the topography imposed on a thin solid body is incompatible with its intrinsic (“target”) metric and, as a consequence of Gauss’ Theorema Egregium, induces strain. Focusing on a prototypical example of a sheet attached to a spherical substrate, numerical simulations and analytical study demonstrate that the mechanics is governed by a variational principle, which we call the “Gauss-Euler elastica”. This emergent rule states that – despite the unavoidable strain in such an incompatible confinement – the ratio between the energies stored in straining and bending the solid may be arbitrarily small. The Gauss-Euler elastica underlies a theoretical framework that greatly simplifies the daunting task of solving the highly nonlinear equations that describe thin solids at mechanical equilibrium. This development thus opens new possibilities for attacking a broad class of phenomena governed by the coupling of geometry and mechanics.

## I Introduction

The spatial confinement of thin, solid objects has been the target of theoretical inquiry since the very early formulations of continuum mechanics. The best known example is “Galileo’s beam”, a problem that is central to solid mechanics, and moreover, to the foundations of variational calculus. Considering the buckling of a wooden beam, Euler recognized that the problem is described by the variational principle Levien (2008),

 Euler′s \emphelastica:  δUbend=0; subject to Ustrain=0 , (1)

where and , respectively, are the contributions to the elastic energy of a solid rod due to bending and straining, and represents the first variation of with respect to shape. Equation (1), known as “Euler’s elastica”, is borne out of two principles that intertweave mechanics and geometry. First, thin solids are energetically far less costly to bend than to stretch, owing to the fact that the ratio of bending to stretching moduli for a rod (or sheet) of thickness varies as   Landau and Lifschitz (1986). Second, confining a rod in space – say, by forcing its ends to be closer than its contour length – can generically be done without straining its centerline. Thus, by restricting configurations to isometric deformations (i.e. ), one guarantees that they acquire minimal energy as , thereby greatly reducing the space of possible equilibria states, making the problem analytically tractable.

Euler’s rule can be extended to two-dimensional (2D) solids under particular circumstances in which the deformed mid-surface can be embedded in 3D isometrically (i.e. without strain) Mansfield (1989), for example, when analyzing the “developable cones” realized by gently pushing a xerox paper into a ring Cerda and Mahadevan (1998). However, when addressing the generic problem of confinement of a thin solid in 3D space, namely, confinements that require a change in Gaussian curvature of the mid-surface, Gauss’ Theorema Egregium implies that a perfectly isometric deformation of the mid-surface is impossible Witten (2007), revoking the applicability of Euler’s variational rule (1) and requiring, instead, minimization of the fully-fledged Föppl-von Kármán (FvK) energy Landau and Lifschitz (1986),

 FvK:   δ[Ubend+Ustrain]=0 . (2)

The unavoidable strain inherent to such geometrically incompatible confinement (GIC), underlies both the impossibility of drawing an accurate map of the globe as well as the frustration of gift wrapping a ball. A quintessential example of GIC is given by the “spherical stamping” experiments of Hure et al. (Fig. 1A), in which a thin sheet (of size and thickness ) is placed in a narrow gap () between two rigid, concentric spherical shells (of radius ) Hure et al. (2012). Here, the positive Gaussian curvature of the confining topography gives rise to complex patterns that cannot be described by Euler’s variational rule, necessitating consideration of the highly non-linear interplay between bending and strain in the FvK equations.

In the article, we aim to bridge the gap between the bending-dominated mechanics of Euler and strain-dominated geometry of Gauss through the “Gauss-Euler elastica” – a generalization of the variation rule (1) that characterizes the mechanical equilibrium of GIC problems. We define GIC as the imposition of a smooth “substrate”, whose shape has a Gaussian curvature, , onto a thin solid sheet or shell, whose “target” metric is characterized by a different curvature, . In addition to , a GIC problem is equipped with two dimensionless parameters, denoted and , which characterize, respectively, the confinement strength, and the bendability of the confined solid. In what follows, we will show how to define these parameters for other GIC problems (e.g. Fig. 1B-F), but to begin, their meaning is best introduced through the spherical stamping problem, where . Considering , the confinement strength and bendability are conveniently defined as:

 χ−1=(W2/2Rδ)2   ;   ϵ−1=(W2/2Rt)2 , (3)

When (i.e. ), the wide gap requires no deflection and the sheet remains flat (i.e. no confinement). If the gap is sufficiently small (i.e. ), the stamp imposes perfect spherical shape, requiring a “bare” geometric strain , due to the elongation of its longitudes and shortening of latitudes. Our study focuses on the intermediate regime, , in which the sheet shape is strongly sensitive to confinement, but can relieve the energetically costly strain down to a residual value, , through deflections within the gap that are only penalized by their bending cost. We propose that these elastic deformations are described by a variant of the variational rule (1):

 Gauss−Euler \emphelastica: δUbend=0; subject to  Ustrain/Ubend→0 , (4)

where “” refers to the doubly asymptotic limit, of high bendability and strong confinement, .

Eq. (4) describes a variational problem that is dominated by bending energy, while being constrained by the high cost of strain, much like the original elastica problem, but with two key and nontrivial distinctions. First, the minimization principle, , means that remarkably, and despite the presence of Gaussian curvature, the elastic energy is dominated by bending rather than strain; thus, somewhat surprisingly, Eq. (4) is closer to Euler’s elastica  (1) than to the generic FvK problem (2). Second, the suppression of strain that is expressed as an exact isometry () in Euler’s elastica, is now imposed in a weaker, asymptotic manner (). Again, this latter distinction derives from the unavoidable distortion of the confined solid from its intrinsic metric imposed by a confining topography with . We will show that the ratio , and the deformed shape of the confined solid, are controlled by the confinement strength and bendability parameters, and , yielding distinct types of energy minimizers at various sectors of the parameter regime, .

We commence with a simpler model of a 2D sheet attached to a ball of stiff, radial springs Grason and Davidovitch (2013); Hohlfeld and Davidovitch (2015); Bella and Kohn (2017). The relative simplicity of this “minimal” GIC problem enables a pedagogic exposition and an explicit derivation of the variational principle (4), which we then employ to study spherical stamping.

## Ii Spherical Winkler foundation

Consider a circular solid sheet of thickness and radius , made of Hookean solid with stretching and bending moduli, and , respectively, attached to a stiff spherical surface of radius ; deflections away from the substrate are penalized by radial, uniformly spaced harmonic springs, each spring with constant (Fig.1G-H). The energy density (per area), , of the system is written schematically as:

 u=ustrain+ubend+usub , (5)

where: and The strain and curvature tensors, and , respectively, are given in terms of the in-plane displacements, , and normal deflection of the sheet from its rest, planar state (see Appendix A). Here, the two dimensionless groups that quantify the strength of confinement and bendability of the sheet are defined as Grason and Davidovitch (2013); Hohlfeld and Davidovitch (2015):

 χ−1=KsubR2Y  ;  ϵ−1=YW4BR2∼(W2tR)2 . (6)

The physical meaning of these parameters can be grasped by considering an ideal, axisymmetric deformation of the system. If is very small, the sheet remains nearly planar and the energetic cost of flattening the substrate beneath it is , whereas if is large the substrate is barely deformed and the sheet conforms to its spherical shape with characteristic energies, and (where ). In this light, the dimensionless parameters in Eq. (6) are seen as the ratios , . Here we focus on the regime:

 ϵ−1≫χ−1≫1    ⟹ χ≪1 , (ϵ/χ)≪1  , (7)

in which is large, and the substrate closely retains its shape to avoid a high energy . However, rather than conforming axisymmetrically to the substrate, which would generate a bare strain and an energy density, , the high bendability of the sheet enables a substantial suppression of strain, down to a residual level, , via formation of small-wavelength wrinkles that “absorb” the excess length of compressed latitudes at the expense of only slight deviation of the substrate from its spherical shape.

The deformed shape and the residual strain can be found via a numerical simulation of the Seung-Nelson model of a 2D elastic sheet bound to spherical Winkler foundation Seung and Nelson (1988); Vernizzi et al. (2011), whose discrete elements capture the full geometrical nonlinearity underlying the FvK energy (Appendix B). In our simulations (Fig. 2) we explore the relevant parameter regime (7) by fixing the confinement strength at a large value, or , and gradually increasing the bendability, (from to ) by reducing the sheet’s thickness. As exceeds a critical value (dashed vertical in Fig. 2G,2H), radial wrinkles emerge in an annular zone, . Increasing further, the wrinkles cover a larger and larger portion of the sheet, and their characteristic wavelength, , becomes smaller (Fig. 2G). We find that:

 λ/W≈2π(χϵ)1/4  ;  L/W≃24/3(ϵ/χ)1/6 (8)

Analyzing the radial and hoop components of the in-plane stress tensor, (which are linear combinations of the respective components of the strain tensor, see Appendix C), we find that unlike the axisymmetric case, for which , here:

 σrr∼σθθ∼Y⋅εres  ,where:  εres≈√BKsub/Y , (9)

(Figs. 2D-2F), such that in the asymptotic limit (7). Let us recall that uniaxial compression of supported sheets gives rise to wrinkles with wavelength of the form (8), reflecting the suppression of a compressive stress component (here ) down to residual level (obtained by balancing sheet bending and substrate stiffness)Bowden et al. (1998); Cerda and Mahadevan (2003). Strikingly, Figs. 2D-2F show that not only the compressive hoop stress, , but also the purely tensile radial stress, , is characterized by this same residual scale.

Further analysis of simulations reveals that the wrinkle-assisted suppression of stress (9) reshuffles the energetic hierarchy in the asymptotic limit (7). Fig. 2H shows that the strain energy is dominant at the vicinity of threshold, but becomes negligible in comparison to the bending and substrate energies in the limit (7). This reversal of the energetic hierarchy can be understood through estimates of respective energy densities which, at first pass, neglect spatial variation. Since Eq. (9) shows that all components of strain are proportional to the residual scale, we have , while , where is the characteristic amplitude of wrinkles, the wavelength, , is given by (8), and . To estimate , we employ a “slaving” condition that links the wrinkle amplitude and wavelength, , expressing the fact that the excess arclength “wasted” by azimuthal wrinkles matches the inward radial displacement of latitudes on the substrate Hohlfeld and Davidovitch (2015); Bella and Kohn (2017). Together with Eq. (8) we thus obtain:

 ubend∼usub∼u(bare)√ϵ/χ  ;  ustrain/ubend∼√ϵ/χ . (10)

This result exhibits two critical features. First, the elastic energy is asymptotically vanishing in comparison to the bare energy, , required for perfectly conforming the flat sheet to a spherical topography. Second, the energetic hierarchy is reversed in comparison to the bare, asymptotic deformation – being governed by bending (balanced by substrate deformation) rather than strain. Taken together, these two features manifest the Gauss-Euler variational principle (4). Notwithstanding the fact that imposing Gaussian curvature generates some strain in the sheet, its energetic cost is negligible and the energy-minimizing state is found by minimizing bending energy (for this case, in balance with the substrate) over a family of “asymptotically isometric” configurations, in a surprisingly analogous manner to Euler’s elastica (1).

## Iii Inverted tension field theory (ITFT)

Motivated by the numerical results, and by classical tension field theory (TFT), which describes an asymptotic, compression-free stress field in confinement problems that are governed by tensile boundary loads Mansfield (1989), we introduce here an “inverted” tension field theory (ITFT) that describes an asymptotic, strain-free stress field, in confinement problems governed by geometrical incompatibility. For the spherical Winkler problem, this theory is formulated as a doubly-asymptotic expansion of FvK equations in and , around the singular limit (7), to which we refer below by the symbol “”. The succinct exposition below is supported by a detailed analysis of the displacement field, force balance equations and energy balance in Appendix C.

1. Hoop confinement: Elimination of radial strain, namely, preserving length of longitudes, requires a finite radial displacement: . This follows directly from the condition , where . Consequently, the projection of sheet’s latitudes onto the sphere is contracted by: .

2. Wrinkle wavelength and hoop stress: Considering each latitude as an elastic hoop attached to a substrate of stiffness under confinement , wrinkles relax the bare compressive strain, , namely: , such that the amplitude, , and wavelength, , satisfy the aforementioned “slaving” condition Davidovitch et al. (2011). Minimizing bending and substrate energies one obtains:

 λ≈2π(B/Ksub)1/4=W(ϵχ)1/4 ; σθθ≈−2B/λ2 . (11)

3. Radial stress: Turning now to force balance in the radial direction, we address the FvK equation:

 ∂r(rσrr)−σθθ=0 (12) ⟹σrr→2√BKsub(−1+W/r) , (13)

where we employed from Eq. (11) as a non-homogeneous source for along with the free boundary condition, . Notably, Eq. (13) manifests the counter-intuitive concept of “bending-induced” tension along wrinkles, which was envisioned already by Hure et al. Hure et al. (2012).

4. Strained core: The spatial divergence of radial stress, in (13), must be alleviated by a strained, unwrinkled core, , whose radius vanishes asymptotically as . The characteristic strain in this core is , yielding a radial (and similarly, hoop) tension, . The continuity of at is required by radial force balance and yields the scaling relation: , in accord with our simulations (8). Notably, the core-assisted regularization of radial stress (13) modifies our estimate of the strain energy (10) by only a logarithmic factor, .

5. Vanishing shear: The arguments above describe the vanishing of hoop and radial stress components as . However, an asymptotically-isometric configuration requires also the vanishing of shear stress, , which is accommodated by way of small oscillations of radial displacement in registry with wrinkles, Paulsen et al. (2016). The energetic cost of these “shear cancelling” radial oscillations leads to curvature-induced stiffness, , that suppresses further the wrinkle amplitude. Taken together with the substrate stiffness, this leads to a generalized version of Eq. (11):

 λ=2π(B/Keff)1/4 ; Keff=Ksub+Kcurv . (14)

(Since , another contribution to , due to tension along wrinkles Paulsen et al. (2016), is negligible here). Since , curvature-induced stiffness is also negligible for the spherical Winkler problem. However, we will show below that the elastic cost of shear suppression is critical for other GIC problems.

The solution of our ITFT equations yields quantitative expressions for the wavelength (Fig. 2G) and residual stresses (solid curves in Figs. 2D-2F), as well as the energies, , and (Fig. 2H), in excellent agreement with simulations. This substantiates the Gauss-Euler elastica (4) and the ITFT equations, whose polar representation is Eqs. (11,14,12), as valuable tools for solving GIC problems.

## Iv ITFT versus TFT

We elucidate how ITFT embodies the Gauss-Euler elastica (4) by comparing the radial force balance equation and radial stress, Eqs. (12,13), with their counterparts in two standard studies of a sheet attached to a spherical substrate (Fig. 3A).

The “bare”, axisymmetric (unwrinkled) state results from solving Eq. (12) by expressing the stress tensor through the axisymmetric displacement field (Fig. 3A, top row) Grason and Davidovitch (2013), yielding a order equation for . The consequent mean strain energy is .

A TFT analysis of this problem (Fig. 3A, middle row) Grason and Davidovitch (2013); Hohlfeld and Davidovitch (2015) is suitable when a tensile load, , is sufficiently smaller than , yet much larger than the residual hoop compression in a wrinkled state, (11). Under such conditions, Eq. (12) becomes a order homogenous equation for the radial stress, and the asymptotic stress tensor has a single non-vanishing component: , which underlies the dominant part of the energy.

In contrast, the ITFT stress yields a subdominant contribution to the energy (). Equation (12) becomes an inhomogeneous equation for , whose source () is set a priori by minimization of the bending energy subject to an (asymptotic isometry) constraint, disregarding the explicit contribution of to the total energy. The asymptotic stress field has the form , acting as a tensorial Lagrange multiplier that results from a constraint (asymptotic isometry) and generalizes the scalar Lagrange multiplier of Euler elastica.

## V The spherical stamping problem

Spherical stamping differs from the Winkler foundation problem in that there is no energetic penalty for deforming the spherical “substrate”. Instead, the energy consists only of bending and straining the sheet, subjected to the constraint that its deflection from the sphere is bounded: . A qualitative understanding of this problem can be obtained by assuming that wrinkles fill the gap, such that the slaving condition, , implies a wrinkle wavelength (up to -dependent pre-factors). The “local law” (Eq. 11 then suggests that the gap effect is akin to a spherical Winkler substrate: (see Appendix DPaulsen et al. (2016). The ITFT equations (11,12) thus yield . Using the results of the preceding analysis and recalling the dimensionless groups defined in Eq. (3) we readily find that for a given ratio , radial wrinkles cover the sheet barring a core of radius . Evaluating the hoop and radial stresses and the various energies, once finds that: and . Hence, the deformation satisfies the Gauss-Euler elastica (4) in the doubly asymptotic limit, .

Note that if (), the gap-induced stiffness, , falls below . In this parameter regime Eq. (14) restricts the wavelength to: ; consequently, wrinkles do not fill the gap (i.e. ). This reflects the energetic cost of supressing shear strain when wrinkling about a curved shape; wrinkles with amplitude larger than would not satisfy , and thus, do not constitute an asymptotic solution of the Gauss-Euler elastica (4).

In order to obtain the actual ITFT solution for the spherical stamping problem, one has to consider a “3-zone” pattern: a “gap limited” zone (), where the amplitude , exploiting the possible deflection at the vicinity of the edge where confinement is strong; a “curvature-limited” zone () – where the wavelength and , so that wrinkles do not fill the gap; and finally a strained, unwrinkled core at (Appendix D). A central result of this analysis is the emergence of two “sub-phases” of ITFT (see Fig. 4B): For (), the pattern comprises only a gap-limited zone, which terminates sharply at an unwrinkled core ( in the above notations); the intervening, curvature-limited zone appears when , and expands (with and ) upon increasing the gap until it reaches the outer edge at (, at which stage a wrinkled state loses contact with the confining shells. Figure 4A presents the energy of our ITFT solution (normalized by ), and the related stamping force , as a function of , for some values of the inverse bendability, , exhibiting two noteworthy features. First, each plot terminates at a finite value, , at which the wrinkle pattern is fully curvature-limited and detaches from the confining shells (signified by a vanishing stamping force). Second, as bendability increases, the energy relaxation becomes more efficient – the residual energy, at scales as .

## Vi Discussion

We introduced the Gauss-Euler elastica (4) for GIC problems – a variational rule where bending is minimized subject to the constraint that strain energy is negligible. We demonstrated the applicability of this principle in two variants of a spherically-confined sheet through ITFT – an asymptotic expansion of FvK equations around a singular, strain-free limit. We showed that the ITFT may be regarded a generalization of Euler elastica, whereby the stress field acts as a tensorial, spatially-varying Lagrange multiplier that enforces an asymptotic isometry constraint. The schematic in Fig. 3B shows how the general, highly nonlinear FvK problem of simultaneous minimization of bend and strain for a confined sheet (2), reduces to three distinct types of simplified variational problems as the sheet thickness becomes small: Euler elastica is applicable for developable deformations (no tension, no Gaussian curvature); TFT is applicable in the presence of tensile loads whose work dwarfs bending energy; ITFT applies for GIC problems, where non-developable deformations emerge in the absence of tensile loads.

The GIC problem we addressed – confinement of a planar sheet onto a spherical topography – corresponds to , where confinement is preferentially azimuthal and the emerging wrinkle pattern is splayed. The GIC problem shown in Fig. 1B is characterized by ; there, confinement is preferentially radial and wrinkles are bent. We expect that the Gauss-Euler elastica (4) is valid also for such problems, and even for more complex GIC problems that result in a mix of splay and bend of wrinkle textures (Fig. 1D).

We emphasize that our definition of GIC problems, and thereby the applicability of Gauss-Euler elastica, precludes a large class of problems in which the confinement is not sufficiently strong to impose a finite, thickness-independent difference, . Well known examples for such “weak confinement” problems are the biaxial confinement of a stiff skin attached to a (planar) compliant substrate Bowden et al. (1998); Audoly and Boudaoud (2008); Kohn and Nguyen (2013), and the contact of curved liquid surfaces with a planar sheet King et al. (2012). The complexity of such problems can be anticipated by considering the two problems we addressed here in the parameter regimes complementary to (7), namely, . In the spherical Winkler problem, this corresponds to the regime , where the substrate is very soft and can deform substantially – flattening uniformly to reduce sheet’s strain and thereby developing stress-focusing, crumpled zones King et al. (2012), or alternatively, enabling the formation of localized folds Pocivavsek et al. (2008); Diamant and Witten (2011); Paulsen et al. (2017).

For the spherical stamping problem, a parameter regime that corresponds to “weak confinement” is, , where radial wrinkles predicted by ITFT no longer fill the gap, and another deformation is required to relax the residual energy, , stored in a wrinkled state (Fig. 4A). To wit, let us consider a very weak confinement, namely, , where (light purple in Fig. 4B); here, sheet is barely touching the confining shells, and the problem seems qualitatively similar to confining a sheet in a ball whose radius is only slightly smaller than its own Aharoni and Sharon (2010). There, the deformation is governed by a few “stretching ridges” that separate the sheet into stress-free facets Witten (2007). Such a deformation is distinguished from a wrinkled state in two key aspects. First, the ridges, which dominate the energy, are characterized by a balance of bending and strain Witten (2007), thereby obeying the general FvK rule (2), but not the restrictive version invoked in Gauss-Euler elastica (4). Second, in contrast to wrinkle patterns, which are described by a mapping of the plane (namely, continuous curvature), for which the stress is well-defined everywhere, ridges become infinitely narrow as , such that the asymptotic deformation is only , “origami-like” mapping (i.e. the sheet’s tangent is discontinuous).

Furthermore, one may envision an intermediate regime (heavy purple in Fig. 4B), namely, , with , at which the asymptotic deformation is a mapping of the planar disk – such that the tangent remains continuous everywhere (even asymptotically), and only the curvature becomes discontinuous. Such a nontrivial isometric mapping is suggested by Nash embedding theorem Nash (1954), whose relevance for the deformations of solid sheets has very recently begun to be explored Gemmer et al. (2016). Whether such Nash mappings do characterize an asymptotically-isometric response of confined solids, and whether their regularized versions (for small, but finite and ) are subjected to the Gauss-Euler elastica (4), are fascinating questions that we hope will inspire future studies.

###### Acknowledgements.
We thank P. Bella, R.V. Kohn, and N. Menon for many discussions, B. Roman, J. Bico and J. Hure for discussing with us their experiment and model Hure et al. (2012), R. Sknepnek for his assistance with numerical simulations, D. Vella, D. O’Kiely, J. Bico and J. Paulsen for useful comments, and participants of the program “Geometry and elasticity of 2D soft matter”, KITP Santa-Barbara 2016. This work was funded by the ACS Petroleum Research Fund Grant No. 54513-ND and NSF Award CBET-14-38425. Simulations were performed using the UMass Cluster at the Massachusetts Green High Performance Computing Center.

## Appendix A A sheet on a slightly deformable sphere

Here we describe in detail the “spherical Winkler” model presented in the main text.

We consider a circular elastic film of radius attached to a “spherical Winkler” substrate, or equivalently a ball made of radial, non-interacting harmonic springs of stiffness and rest length (Fig. 1 of main text). We assume that the sheet to be perfectly attached to the substrate (i.e. an infinite adhesion energy); however, the force exerted by the deformed substrate on the sheet is purely normal to the sheets midplane, with no tangential component (such that the sheet slides freely on the substrate surface). This model captures the essential mechanics of adhesion onto a curved isotropic solid, here simplified by assuming that the substrate response is local and without tangential forces. This “spherical Winkler” model was described in Sec. II of Hohlfeld and Davidovitch (2015); here we will briefly review its main features, focusing on the case of a sheet with stress-free boundary.

Since , the displacement can be expressed through the Monge parametrization:

 u(r,θ)=ur(r,θ)^r+uθ(r,θ)^θ+ζ(r,θ)^z , (15)

with the components of the strain tensor, , given by

 εrr=∂rur+12(∂rζ)2 , (16a) εθθ=1r∂θuθ+1rur+12r2(∂θζ)2 , (16b) εrθ=εθr=12(1r∂θur+∂ruθ+1r∂rζ∂θζ) , (16c)

The stress in the sheet is given by the Hookean relationship Landau and Lifschitz (1986); Timoshenko and Goodier (1970); Mansfield (1989):

 σrr=Y1−Λ2(εrr+Λεθθ) , (17a) σθθ=Y1−Λ2(εθθ+Λεrr) , (17b) σrθ=Y1+Λεrθ , (17c)

where the Poisson ratio of the sheet. Furthermore, anticipating the shape to be characterized by small slopes (even in the wrinkled zones), the curvature tensor can be approximated as:

 κrr=∂2rrζ  ;  κθθ=1r∂rζ+1r2∂2θθζ  ;  κrθ=2r∂2rθζ (18)

The energy () is conveniently expressed through its areal density ():

 u=uFvK+usub , (19a) where uFvK=ustrain+ubend , and: usub=Ksub2(ζ+r2/2R)2 , ustrain=12σijεij  ;  ubend=B2Tr(κ)2 (19b)

(Eq. 5 of main text). In Eq. (19b) we once again employed the inequality, , approximating the rest (spherical) state of the substrate as . Furthermore, we will see later that the bending energy is governed by the rapid, azimuthal undulations of the shape, and hence can be well approximated as (i.e. variable shape along radial direction leads to negligible bending energy).

Choosing appropriate units of energy and length, one finds that the model energy (19) depends on two dimensionless parameters (Eq. 6 of main text):

 ϵ=BR2YW4∼(√tRW)4  ;  χ−1=KsubR2Y , (20)

whose physical meaning was discussed in the main text. As we explained there, the focus of our study is the asymptotic regime:

 ϵ−1≫χ−1≫1 , (21)

where the sheet is highly bendable and the substrate is barely deformable. A direct counting of dimensionless groups yields a third parameter, , but in the Föppl ván Kárman (FvK) framework it can be absorbed into a global energy scale Hohlfeld and Davidovitch (2015).

In order to find the displacement field it is useful introduce the Euler-Lagrange equations of the energy functional (19), which are essentially the FvK equations for the mechanical equilibrium of the sheet, subjected to the restoring force :

 BΔ2ζ−σrr∂2rζ−2rσrθ(∂r−1r)∂θζ −1r2σθθ(∂2θζ+r∂rζ)=−K[ζ(r,θ)+r2/2R] , (22a) ∂rσrr+1r(∂θσrθ+σrr−σθθ)=0 , (22b) ∂rσrθ+1r(∂θσθθ+2σrθ)=0 , (22c)

where the Laplacian . Eq. (22a) is the FvK equation and expresses force balance in the normal direction; whereas Eqs. (22b,22c) form the FvK equation, and express force balance on each infinitesimal piece of the film in the two directions locally tangent to the sheet Landau and Lifschitz (1986); Timoshenko and Goodier (1970); Mansfield (1989).

It is useful to consider the axially symmetric state, which is the solution for the case of an infinitely rigid substrate (). In this case the effect of the substrate can be expressed as a constraint – the sheet must conform perfectly to a sphere.Clearly, the whole displacement field is axially symmetric, such that , , and . Since any explicit dependence on vanishes, the only non-vanishing components of the stress tensor remain and . Furthermore, for an infinitely stiff substrate, the normal force balance (Eq. 22a) is replaced by the geometrical constraint, , and the only remaining equation is the radial force balance, supplemented by the strain-displacement Eq. (16), akin to the compatibility condition Azadi and Grason (2014a). The radial displacement , as well as the radial and hoop stresses for the axisymmetric state are given in Azadi and Grason (2014b); Azadi16. Crucially, the only stress scale that governs the magnitude of both stress components, and , of the axisymmetric state, is , and the energy of such a state – which is obviously dominated by strain – is therefore .

## Appendix B Simulation on sheets on a spherical Winkler foundation

### b.1 Numerical model

We employ the “bead-spring” model of Seung and Nelson Seung and Nelson (1988) implented in the MEMBRANE code developed by Skenpnek et al. Vernizzi et al. (2011). The sheet is modeled by an initially planar, equitriangular array of vertices , connected by nearest neighbor springs of rest length (all length scales are measured in units of the stress-free lattice spacing). Upon deformation, in-plane strain energies are captured by spring stretching energy,

 Estrain=k2∑⟨αβ⟩(|xα−xβ|−a)2, (23)

where the sum is taken over nearest neighbor bonds and is a spring constant. To model the bending energy of the sheet, normal vectors are associated to each triangular plaquette, and the bending energy derives from the dihedral rotation of normals between adjacent plaquettes

 Ebend=κ2∑⟨μν⟩∣∣nα−nβ∣∣2, (24)

where the sum is taken over neighbor (edge-sharing) plaquettes and is an angular stiffness. In the limit , this discrete model tends towards the continuum limit of an elastic membrane with stretching and bending moduli, eq. 19b, with modulii:

 Y=2√3k; B=√32κ, (25)

which we use to compare numerical simulations to analytical predictions. Finally, we consider an elastic interaction between the sheet and the spherical substrate of preferred radius of the form of a Winkler foundation that penalizes local deflection of the substrate in the radial (normal) direction,

 Esub=K2∑α(|xα|−R)2 (26)

where the origin is taken as the sphere center and is a stiffness parameter for the substrate. In the limit of small triangles (and small strains) this energy tends towards the continuum form of eq. 19b, with

 Ksub=2√3a2K. (27)

Given an initial configuration of the sheet, the simulation proceeds via energy minimization performed using the BFGS. The convergence is determined by requiring that the sum over absolute value of derivative of energy on all vertices should be smaller than a threshold value. A threshold level of was established empirically by confirming that two benchmarks of the shape profile, the extent of the wrinkled zone the wavelength of wrinkles, do not change by more than a few % upon further decrease of the threshold.

Initial configurations are constructed by cutting planar triangular meshes from within a circular “stencil” of radius . To ensure that the discretization has minimal impact on simulated wrinkle patterns, we require sufficiently large ratios of so that the ultimate patterns satisfy the criteria over the full range of simulated parameters. Practically, this was achieved using sheets with corresponding to 326,467 vertices. Energy relaxation begins from a state with the planar sheet projected orthographically onto the sphere of radius . Bending modulus and substrate stiffness were varied over a range of parameter values (corresponding to ) and (corresponding to ).

### b.2 Simulation analysis

Simulation results are analyzed from two morphological feature: the radius of the unwrinkled core and the wrinkle wavelength , as well as the profiles of radial and hoop stress. To analyze the shape, vertex positions are converted into the deflection profile: , the normal deflection of a vertex from the sphere, its arc radius and azimuthal angle around the pole, . The vertex positions are then binned into radial annuli of width , and then the mean number of wrinkles at a given arc radius is determined from counting the number of amplitude oscillations (crests or trothes) in the deflection profile plotted vs. . The length is determined via two algorithms. In the first, which is implemented numerically for , the radially averaged bending energy density of simulated sheets is plotted as a function of radial position. The start of the wrinkled zone is taken as a radius at which the bending energy density has the largest second derivative (computed via interpolation), and the error bars as estimated as the full width, half maximum of this second derivative. When the edge of the wrinkled zone approaches the edge of the sheet (i.e. when ), this numerical algorithm fails, and is determined from plots of vs. , as the extrapolation of the apparent wrinkle envelope at to .

The stress profiles in simulated sheets are computed using the virial formula for a vertex position

 (28)

where the sum is taken over nearest (bonded) neighbors of , is the area per vertex, and force from the bond is

 f(α,β)=k(|xα−xβ|−a)^xαβ, (29)

where . The radial and hoop components (plotted in Fig. 2 in the main text) are computed by projecting the stress tensor onto the local radial direction that points along the longitude from to , where .

## Appendix C Basis of ITFT

In main text we introduced the inverted tension field theory (ITFT) equation for the spherical Winkler problem. In this exposition we explicitly assumed that, despite the strong oscillations (wrinkles) of the shape in the azimuthal direction, the residual strains and consequently stress field are governed by force balance equations that depend only on the radial direction, . Here we justify this assumption by considering also the leading oscillatory terms in the displacement field and consequent strain components. We will show that the ITFT equations (11-13 in the main text) are essentially the leading order in an expansion of the FvK equations around an asymptotically strainless state, obtained in the singular limit (21), where the sheet becomes asymptotically strainless in the sense of the Gauss-Euler variational principle, Eq. 4 of main text.

We further employ this asymptotic expansion to evaluate the total energy of the wrinkled state, , in the limit (21). Normalizing energies by the factor , which is the energy (up to a a numerical pre-factor) of the axisymmetrically deformed (unwrinkled, strained) state, we find that the leading contribution to the energy originates from bending the sheet and substrate deformation:

 Ubend+Usub→YW6R4⋅π6√(ϵ/χ) [1+O(√(ϵ/χ))] , (30)

and the strain energy contributes only a sub-leading term:

 Ustrain→YW6R4⋅2π3(ϵ/χ) [−log(ϵ/χ)+O(1)] . (31)

The asymptotic expansion of the energy, given by Eqs.(30,31), provides a self-consistent confirmation to the Gauss-Euler elastica principle (Eq. 4 of the main text): The minimum energy of a GIC problem (in the limit (21)) is characterized by a dominant energy associated with bending and substrate deformation, and a sub-dominant energy due to residual strain.

We adopt a variational method to approximate the energy via a parametric family of states that assumes radial wrinkles develop in an annular zone, , anticipating that the sheet becomes fully covered with a diverging number of wrinkles () in the singular limit (21). Hence, we have

 L

where . Our objective is to find the radial profile , the number of wrinkles (which may vary with ), the core radius , and the other components of the displacement field, and , that minimize the total energy . Central to our reasoning is the high energetic price that comes with any level of strain due to the large contrast between the stretching and bending moduli. The overarching principle of our variational construction in this section is thus to keep strain at the smallest possible level, while requiring the sheet to remain as close as possible to the spherical substrate. Implementation of this principle leads to a much stronger result, shown in the next subsection, that there exists a wrinkle pattern for which the energy is negligible in the limit (21) in comparison with and . Finally, we will show that, despite this overall negligibility of strain, its minimization has a dominating effect on the separation between wrinkled and unwrinkled zones. The spatial distribution of becomes concentrated in the limit (21) in the (internal) margin of the wrinkled zone, yielding a new type of “stress focusing” phenomenon, which underlies the core size, .

### c.1 Asymptotically vanishing strain

We start by requiring the radial strain to vanish in the limit (21). Using Eqs. (16a, with ), the boundary conditions , and recalling our assumption that , we obtain:

 εrr→0  ⇒  ur→−16r3R2  ⇒  ~Δ(r)≡−urr→16r2R2 , (33)

as was mentioned already in the main text.

The above limit characterizes the leading term in the radial displacement, which remains finite in the limit (21). Establishing the hierarchical structure of the energy, as described in the next subsection, requires consideration of the displacement beyond leading order. Anticipating this necessity, we generalize the above expression for as

 ~Δ(r)→16r2R2+gΔ(r) , (34)

where the function is required to vanish in the limit (21).

#### Azimuthal strain

Inserting the ansatz (32), into Eq. (16b), and requiring hoop strain to vanish we obtain from the non-oscillatory (i.e. -independent) part of :

 εθθ→0  ⇒  14r2f2m2→~Δ(r) , (35)

where the azimuthal confinement is given by Eq. (33). Anticipating again the analysis in the next subsection, we complement the above expression with a higher-order term:

 14r2f2m2→~Δ(r)+gf(r) , (36)

where the function is required to vanish in the limit (21).

At a leading order approximation, where , Eq. (36) is a “slaving condition” between the amplitude and wavelength () of wrinkles Davidovitch et al. (2011); Hohlfeld and Davidovitch (2015): both and become indefinitely small in the limit (21), but their ratio must remain finite in order to allow collapse of hoop and radial components of the strain tensor.

The oscillatory (i.e. -dependent) part of Eq. (16b) shows that an asymptotically vanishing hoop strain implies an additional constraint on the azimuthal displacement :

 εθθ→0  ⇒  uθ→+r2m~Δ(r)sin(2mθ) , (37)

where we used Eq. (36). Once again, we allow for higher-order contributions to as we elaborate in the following subsection:

 uθ → rgθ1(r)sin(mθ) + r2m[~Δ(r)+gθ2(r)]sin(2mθ) , (38)

where the functions , are required to vanish (faster than ) in the limit (21).

#### Shear strain

The last component of the strain tensor is the shear. Inserting the ansatz (32), into Eq. (16c), and using the above result, Eq. (35), one may be tempted to conclude that the shear component approaches a finite value: . However, an asymptotically vanishing strain is retained (at an ultimately nominal energetic cost) through an oscillatory addition to the radial displacement:

 εrθ→0⇒  ur→−r~Δ(r) + 1m2r2R [√~Δ(r)+gr1(r)] cos(mθ) , (39)

where is yet another function which is required to vanish in the limit (21). Note that the oscillatory correction to the radial strain, , induced by the necessity to eliminate shear, scales inversely with , and is thus sufficiently small such that the radial strain still vanishes in the limit (21). As we will show in the next subsection, this oscillatory contribution to underlies the “geometric stiffness”, , Eq. (14) of the main text Paulsen et al. (2016); Hohlfeld and Davidovitch (2015).

Reconsidering Eq. (16c), one may jump to the conclusion that the residual shear strain is . Since in the limit (21)), such a contribution seems harmless; however, we will see later that it is inconsistent with an asymptotic vanishing of other strain components. Inspection of the various terms in (16c) reveals that the actual shear strain can be reduced to terms of order or higher, by including a higher-order contribution to the :

 ur→−r~Δ(r)+1m2r2R [√~Δ(r)+gr1(r)]cos(mθ) +1m2r[−34~Δ(r)−14rddr~Δ(r)−gf(r)+gr2(r)]cos(2mθ) , (40)

where we introduced a final function, , that must also vanish in the limit (21).

### c.2 Dominant and sub-dominant energies

We showed above that requiring all strain components to asymptotically vanish provides strict constraints on the form of the displacement field , as expressed in Eqs. (34, 36,38,40). The ability to eliminate strain suggests that in the limit (21), the total energy is not “equipartitioned” between the various deformation modes, but instead – similarly to buckling and wrinkling phenomena in one-dimensional (1D) systems – becomes dominated by the bending and substrate deformation costs, and , rather than by . We will prove this nontrivial claim in a self-consistent manner, evaluating first the limit value of