Faceted wrinkling by contracting a curved boundary

Single-mode deformations of two-dimensional materials, such as the Miura-ori zig-zag fold, are important to the design of deployable structures because of their robustness; these usually require careful pre-patterning of the material. Here we show that inward contraction of a curved boundary produces a fine wrinkle pattern with a novel structure that suggests similar single-mode characteristics, but with minimal pre-patterning. Using finite-element representation of the contraction of a thin circular annular sheet, we show that these sheets wrinkle into a structure well approximated by an isometric structure composed of conical sectors and flat, triangular facets. Isometry favours the restriction of such deformations to a robust low-bending energy channel that avoids stretching. This class of buckling offers a novel way to manipulate sheet morphology via boundary forces.

Single-mode deformations of two-dimensional materials, such as the Miura-ori zig-zag fold, are important to the design of deployable structures because of their robustness; these usually require careful pre-patterning of the material.Here we show that inward contraction of a curved boundary produces a fine wrinkle pattern with a novel structure that suggests similar single-mode characteristics, but with minimal pre-patterning.Using finite-element representation of the contraction of a thin circular annular sheet, we show that these sheets wrinkle into a structure well approximated by an isometric structure composed of conical sectors and flat, triangular facets.Isometry favours the restriction of such deformations to a robust low-bending energy channel that avoids stretching.This class of buckling offers a novel way to manipulate sheet morphology via boundary forces.
The flourishing field of extreme mechanics [4][5][6][7] seeks ways in which large, spontaneous material deformations yield robust, reproducible structures and motions.Many of these structures exploit the selective deformations shown by thin sheets, which bend or fold easily but resist internal stretching strongly [8][9][10][11][12][13][14][15][16][17][18][19] .The large resistance to stretching creates strong, non-local constraints that restrict the sheet to selected modes of deformation.A classic example is the Miura-ori folding pattern 1,20,21 .This periodic pattern enables one to contract a spread-out flat sheet of paper into a compactly folded state with one single mode of deformation dictated by the pattern (see Fig. 1a).Here we report similar robust deformation of a flat sheet without the position-dependent processing needed for an origami fold.Instead, we specify an initial annular shape and a simple mode of edge deformation: contracting the inner boundary radially.The result is a strongly reproducible three-dimensional pattern of high-amplitude wrinkles.We attribute the reproducibility to a qualitative feature of the simulated, finite-element simulation-each wrinkle approaches an isometric shape, with no stretching.This unstretched morphology is made possible by a distinct pattern of flat, triangular facets alternating with smoothly-curving conical segments.
We consider a minimal form of radial contraction that we call "inner Lamé" contraction.We displace the inner boundary of the annulus with no further constraints.Our interest in this geometry arose from experiments on the curling of edge-connected polymeric sheets with unequal solvent swelling.Ref. 22 shows the † Electronic Supplementary Information (ESI) available: [details of any supplementary information available should be included here].
See DOI: 10.1039/cXsm00000x/ result for two half sheets with different swelling separated by a straight boundary.We simulated the radial analogue of this situation, where the straight boundary was replaced by a circular boundary separating an inner region with little swelling from an outer one with more swelling.The swelling resulted in fine radial wrinkling like that shown in Fig. 1.The present simplified system was developed to understand this wrinkling.
As shown in Fig. 1c, the inner Lamé contraction produces fine wrinkling like the conventional Lamé contraction of Fig. 1b.However the two wrinkled states are radically different.The conventional Lamé wrinkle wavelength is explained in terms of the applied tension 2 .However, in the inner Lamé, the exterior tension of the Lamé state is absent, and so the established mechanism for the wrinkle wavelength cannot apply.Indeed, our wrinkling contrasts with most wrinkled systems, where the bending energy favoring longer wavelengths is countered by the energy associated with some external forcing 9,23 .The wrinkle wavelength then depends explicitly on the strength of this forcing.The inner Lamé system has no such strength parameter.Instead, the deformation is defined in terms of the geometric displacement of the inner edge only.
Even among geometrically-defined forms of wrinkling, the inner-Lamé geometry is anomalous.A priori, it doesn't conform to typical one-dimensional buckling treatable as an Euler elastica system 24 .Other forms in which the constraint is applied at one boundary 14,25 do not give a well-defined wave number as we find.
Furthermore, we note that the forms of wrinkling mentioned above are explained as equilibrium ground-states of the system.In our system the wrinkles emerge from a quasistatic deformation.That is, the inner displacement is applied gradually, and the configuration is allowed to relax at every stage.This procedure Fig. 1 Geometry-controlled robust radial wrinkling.a) Miura-ori patterns guide a flat sheet to contract into a folded state with a single, continuous motion (adapted from 1 ).b) Smooth radial (i.e., Lamé) wrinkling in an annulus under both inner and outer tension, indicated by arrows (adapted from 2 ).In this paper, we study Miura-ori like robust deformation of an annulus in a radial Lamé setting.This deformation can be generated in multiple ways.c) Finite-element calculated shape under uniform radial contraction of the inner boundary, described in Methods (inset shows detail of emergent zigzag inner boundary).d) Circular-cone-triangle geometric construction defined in Sec. 2, with triangles shown in blue and cones in red.e) Ruff collar garment (reproduced from 3 ) made by manually pleating the inner boundary of an annular piece of cloth.
corresponds to how one would perform the deformation experimentally.In this limit, the system follows a given energy minimum as the deformation is imposed.It need not be an energy ground state.
In this paper, we address the structure of the individual wrinkles resulting from this unusual mode of deformation * .We show that a strain-free isometric shape quantitatively accounts for observed features of our simulated shape, notably the elastic energy shown in Fig. 5.This energy contrasts strongly with the much higher energy expected with a simple sinusoidal wrinkling pattern (Sec.3.3).Moreover, the simulated deformation matches an explicitly isometric deformation of the original annulus into a cone-triangle pattern.We describe the robustness of the structure using a variety of initial and boundary conditions for the contracted inner boundary.We argue that the observed robustness of the cone-triangle structure is a natural consequence of its isometric property.We note how the cone-triangle structure can be induced via different paths, as illustrated in Figs.1c, d, and e, and discuss generalisations to a broader range of structures.Using the findings of this paper, we address the question of wavelength determination in a separate work 26 .
Sections 1 and 2 define the inner Lamé deformation that we simulate, and the isometric model shape used for comparison.In Section 3, we show that the shape and energy of the simulated sheet are consistent with the model.Section 4 addresses

Inner Lamé deformation
Given the circular annulus of Fig. 2a, we define the inner Lamé deformation as the equilibrium shape induced by drawing the inner boundary inward by a distance ∆, so that it is forced to live on a cylinder whose radius is reduced by ∆. † We define the initial inner radius as our unit of length.The pull distance ∆ thus acts as the sole loading parameter in the system.The system can thus be conveniently defined using only three geometric parameters: thickness t, the width w, and displacement 'loading' ∆.The internal forces determining the shape arise from the bending modulus B and the in-plane stretching modulus Y .Here, we seek the asymptotic behaviour of our annulus in the thin limit: t ≪ w;t ≪ 1.
The conventional "FT" reasoning used to explain classical Lamé wrinkling 2 suggests that fine wrinkling like that of Fig. 1. would require negligible radial stretching.In the FT regime, the wrinkling slopes grow to the order of unity, and the associated azimuthal stress becomes negligible compared to the pre-buckling stress.Only the bending azimuthal stress, owing to wrinkling, proportional to the bending modulus and the curvature, remains.The same is expected for our inner Lamé case: consider the radial stress in a narrow sector of the annulus (see Fig. 2c).Equilib- When the bases of the triangles are moved inward, they must tilt to retain their shape and connectivity.The bridging sectors are squeezed laterally to form cones.A typical deformed example is shown in Fig. 1c.c) Top view of a narrow sector of the sheet, with detail of an area element at the outer edge of Lamé sheet showing radial (horizontal) and azimuthal (almost vertical) boundary forces.In the fully buckled (FT) regime 2 , the azimuthal forces nearly vanish (as indicated by the upper and lower X's).In the the present "inner Lamé" system, the outer radial force also vanishes, indicated by the X at right.Thus, the inward radial force must also nearly vanish for force equilibrium.All elements to the left of the pictured one are subject to the same argument.Thus except for bending stresses, all radial stresses vanish.
rium requires that the net radial force on the sector be zero.In the unbuckled state, the inner tension driving the deformation is balanced by the net outward component of the azimuthal stresses on the sides of the sector.But in the FT regime, this stress is negligible, as argued above.In addition, there is no radial tension at the outer boundary.We conclude that the only forces to balance an inner tension are the weak bending stresses.
The absence of any deformation stress except for bending stress would imply that the sheet deforms only by bending in the thin limit.That is, it is isometric to the undeformed flat annulus.Such an isometric limit can be attained by wrinkles of diverging fineness 16,17 .But if a limiting smooth shape is to exist, then it must obey the geometrical property of "developability".That is, in view of Gauss's Theorema Egregium 27 , there must be a straight generator line through each point of the surface that extends to the boundary without curving in space 28 .In the next sections, we show that the contracted annulus is indeed consistent with a specific developable surface, suggested by the numerically determined shape.

A developable cone-triangle model
The faceted inner Lamé morphology of Fig. 1 suggests a further hypothesis: that the wrinkled morphology can be seen as a piecewise union of alternating flat triangles and curved conical sectors.We take this as the starting point for our geometric model.For a configuration with wavenumber m, we consider a flat configuration (see Fig. 2b) where the inner boundary is a 2m-sided polygon.Each segment of this boundary is the base of an isosceles triangle of height w.To complete the annulus, the apexes of adjacent triangles are joined with an arc.The constructed surface then consists of these triangles and the sectors joining them.
One may readily specify an isometric deformation that allows the inner boundary of this structure to contract.We deform the initial surface by translating the base of each triangle inward toward the centre, specifying that the triangles remain rigid.In order to maintain their shape and connectivity, the triangles must tilt (i.e.rotate about their mid-lines) increasingly as they are drawn inward; adjacent triangles tilt in opposite directions.The sectors connecting two triangles must then also bend in order to remain connected to their triangles.A simple isometric choice for this deformation is to bend the flat sectors into circular conical arcs as the boundaries of each sector squeeze together.These arcs do not in general join smoothly with their adjacent triangles, as further discussed in Sec.3.3.

Comparing simulated deformation with model
In this section, we test whether the isometric/developable conetriangle construction is able to predict the quantitative features of our simulated inner Lamé deformation from Sec. 1.Our finiteelement analysis (FEA) software Abaqus uses standard methods to represent the deformed equilibrium shape as the (cylindrical) inner boundary is gradually contracted inward.The inner boundary is allowed to buckle by moving along the cylinder (see Methods).For simplicity we restrict this motion to the axial direction, and allow no azimuthal motion of the boundary points.We describe the effect of this restriction in Sec.4.3.4.In comparing this calculation to the cone-triangle construction, we must supply two free parameters: the boundary displacement ∆ and the wavenumber m.Thus, we record the ∆ and m of the finite-element deformation, and then construct the corresponding geometric solution with the same ∆ and m.The first comparison we make is purely qualitative -we superpose the two solutions and observe how well they match up.Fig. 3 shows such a comparison, where we have rotated the geometric solution about the cylinder axis to obtain the best superposition.As the figure shows, one may align the two figures so that they overlay well over a significant part of the annulus, with matching zigzags at the inner boundary and matching waves at the outer boundary.Between the inner and outer boundaries, the surfaces match closely enough that the two surfaces alternately hide each other.
The qualitative resemblance seen in the figure encourages a more quantitative comparison.In particular, we propose that for a given wave number m, the shape of the sheet approaches an isometric cone-triangle structure in the limit of zero thickness.Further, approximating each cone as a circular cone sector that con-J o u r n a l N a me , [ y e a r ] , [ v o l .] , 1-16 | 3 Fig. 3 Qualitative comparison between simulations and our cone-triangle model for a representative sample (w = 1.67, t = 1.33 × 10 −3 ).(Left) Given w, ∆ and the wavenumber m from the numerical solution, we construct a model cone-triangle solution (cones in red, triangles in blue) using the procedure described in Section 2. (inset) Zoom-in of a single conical sector (in red) with the entire cone shown (in cream).S marks the length of the outermost arc (as in Fig. 2b) after deformation.For details of construction, see SI. (Right) For a qualitative comparison, we superpose the two solutions, with azimuthal rotation as the only degree of freedom.We find close agreement between them; with this top-down view, areas with colour indicate where the model solution is above the numerical solution, while those in grey indicate the opposite.In particular, the zoomed-in boxed portion shows close matching at both inner and outer boundaries.Specifically, the two inner boundaries show a regular zigzag shape of the same amplitude, which extends to qualitatively flat facets and terminates in a smoothly undulating outer boundary that also has similar amplitude in the two sheets.Zooming out to the whole annulus, we see that there is a 'phase mismatch' between the two solutions, with part of the solution in-phase (the boxed region) and part of it out-of-phase.This is due to the non-uniform size of the wrinkles in the numerical solution.
nects its neighboring triangles gives a useful match to the shape and energy of the sheet.
Accurate comparison is hindered, however, by the variability of the simulated wrinkles.Though regular in general, the simulated wrinkle pattern is not strictly periodic on the annulus.To optimise the comparison, we perturb the initial flat state of the numerical sheet to promote periodicity.Specifically, we deform the initial state with a low-amplitude, smooth, periodic wrinkle with a wavelength similar to the original simulations.‡

Inextensibility and bending strains
We may roughly verify the anticipated isometry of the sheet as announced in the Introduction.To this end, we measure changes in lengths of the material circles at the inner (r = 1) and outer (r = 1 + w) boundaries.
First, we consider radial inextensibility.Here, we use approximations appropriate for fine buckling and narrow wrinkles as in the far-from-threshold regime of Ref. 2 .A material circle at radius r in the annulus has fixed initial arc length L 0 (r) = 2πr.After the contraction, this circle deforms into a wavy line, whose horizontal projection is approximately circular under the observed fine ‡ The initial wrinkle shape chosen is an eigenmode obtained from linear stability analysis performed using the Abaqus software.
buckling.The radius of this contracted circle is the original r plus some radial displacement e r due to the inward contraction.The fractional change in the projected length L xy is thus given by L xy /L 0 = 1 + e r /r. ( At the inner boundary (r = 1), e r = −∆ because of the imposed boundary conditions.However, at the outer boundary (r = 1 + w), e r is free to assume any value.If it's different from −∆, this implies the presence of radial strain.In Fig 4a, we use Eq. 1 to provide a simple estimate of the radial strain.For the outerboundary displacement we infer e r from Eq. 1 and subtract the inner boundary displacement ∆.This gives the change in the horizontal component of the radial length.We interpret this small negative apparent strain in Sec.4.3.2.
We may gauge azimuthal extensibility by measuring the full (not projected) change of arc length for the circle at radius r.
Here we consider only the outer boundary region.As reported in Fig 4b, there is a nonzero compression of these circles.This compression reflects the Elastica bending stress of pure bending mentioned in Sec. 1. Specifically, bending at a scale λ produces a residual bending stress, σ resid θ θ ∼ −B/λ 2 , where B is the bending modulus.Thus, from wrinkling of wavenumber m, we expect the outer boundary to have an azimuthal bending strain given by: where we have used λ = 2π(1 + w)/m and B/Y ∼ t 2 (t is the thickness).The strain ε resid θ θ is calculated for our numerical samples by directly measuring the change in azimuthal bond lengths along a material circle slightly inside r = 1 + w § , and then averaging.
ε resid θ θ is negative as expected; moreover, it is roughly independent of ∆ beyond a transient initial phase.Fig. 4b shows a log-log plot of |ε resid θ θ | for multiple samples with differing values of m, t and w.We see data collapse over the range of a decade, with the predicted slope of unity.We attribute the large variability of the points to variations in the specific wrinkle pattern, which change from one simulation to the next.Based on the linear fit, the estimated pre-factor for Eq. 2 is 0.24.We note that this is roughly 50 times greater than the pre-factor calculated theoretically for purely rectangular Elastica bending, which gives σ resid = −2B/λ 2 17 .Such a larger value, however, is not unexpected given the inherently two-dimensional nature of our conical Elastica arcs.
Figs. 4a and 4b together provide suggestive evidence for radial and azimuthal inextensibility of the numerical solutions.

Evolution of wrinkle amplitude
We now compare specific features of the wrinkle shape to the predictions of the cone-triangle construction.In Fig. 4c, we consider the amplitude of the observed wrinkles.A dimensionless way of measuring the amplitude is by determining the average slope angle α of the wrinkles at the inner and outer boundaries.Though the cone-triangle construction can be implemented for a wide range of wavenumbers m for any given ∆ and w, we considered only the narrow range of m that appeared spontaneously in the unperturbed numerics, where 7 ⪅ m ⪅ 35 26 .
We define the average slope angle α to be the slope angle of the line joining the zero-crossing of the wrinkle with the adjacent maximum (illustrated in inset of Fig. 4a).For angle α in at the inner boundary, the figure explicitly relates α in to L xy /L.For fine wrinkles, we may use Eq. 1 to get: where we have introduced the reduced variable ∆(r) ≡ ∆/r.At the outer edge, L xy /L is likewise the ratio of the horizontal line to the arc length of the circular sector (see inset of Fig. 4a).For our circular cone model, it is related to the slope angle α out by In Fig. 4c, we compare these cone-triangle slope angles to the simulated ones as a function of ∆(r) for both inner (r = 1) and outer (r = 1+w) boundaries.The upper red curve and data points compare α in .The agreement of the data points with the curve confirms that each triangle segment tilts to preserve its length § This is done to avoid possible strain originating in Gaussian curvature boundary layers at the free outer boundary of the annulus; see 29 .
as it is drawn towards the axis.The solid black curve gives the cone-triangle prediction of Eq. 4 for α out .The data points on and beneath the curve represent annuli of different widths w and wavenumber m, coloured by their splay: η ≡ π(1 + w)/(mw).In general, our circular cone model is consistent with the prediction for narrowly-splayed wrinkles, but over-estimates the slopes for widely-splayed ones.The pictures below the graph show the annulus shape for different amounts of splay η.They suggest the source of the over-estimated slopes: for the widely splayed example with the greatest discrepancy, the boundary curve is significantly flatter than a circular arc.Thus the simplifying assumptions of our circular-cone construction are not satisfied in this widely-splayed regime.

Comparing energies
Perhaps the main significance of this isometric deformation is that its energy behaves so differently from that of a generically deformed sheet.It is different in two respects.Firstly, the energy depends only on the bending modulus B, and is independent of the stretching modulus Y for given bending modulus.Secondly, the energy in this limit is qualitatively smaller than for a generic wrinkling deformation with curvature and thickness comparable to our simulations.Below we verify that the energy of the contracted annulus has these distinctive features.The elastic energy is a sum of strain energy and bending energy: Thus, for a strain-free isometric configuration, we have U elastic ≈ U bend .For an arbitrary conical deformation extending between lengths L min and L max along the cone axis, the bending energy is given by: where B is the bending modulus and κ N (s) is the curvature of the outermost arc of the cone normal to the surface, parameterised by the arc length s.For a circular cone as in our model, κ N (s) = const..For fine wrinkles ¶ , we can write κ N ≈ 1/ρ c , where ρ c is the radius of curvature of the outer circle of the cone (see SI Fig. 7b).Then, adding up the contributions of 2m cones for a wrinkled solution of wavenumber m (see SI for details), we get: Also, in order to avoid high-curvature regions near the inner boundary, where the thin-sheet approximation κt ≪ 1 becomes questionable, we choose to compare energies for only the outer 3/4 th part of the annulus, so that log( L max L min ) = log(4).In Eq. ( 6) only the quantity ρ c (∆) is not explicitly known.Our circular cone construction (see SI) permits us to get an expression for ρ c (∆) in terms of the outer slope angle α out : We recall that α out depends only on w and ∆, and not on m.Thus ρ c ∼ 1/m, and so U model elastic ∼ m 2 , which is the quadratic scaling with wavenumber (at fixed amplitude) expected for a bending energy.
Our primary subject of interest, however, is the scaling of U model elastic with ∆ for fixed m.Since α out (∆) ∼ √ ∆ (see Fig. 4c), we expect the squared curvature and the bending energy to scale linearly with We note that this linear form is the same as for the Elastica energy of a rod under weak Euler buckling 30 , which is expected since each transverse arc of the cone undergoes Euler buckling with this same ∆ dependence.
We may compare this isometric energy with that of conventional wrinkles 17 whose height h(r, θ ) has the separable form h = f (r) sin(mθ ).As discussed in the SI, radial wrinkles of this form have substantial strain ε c of order of ε c ≃ ∆/m 2 , much larger than the bending strain of Eq. 2, viz.≃ m 2 t 2 .Accordingly, for fixed ∆ and m, the conventional strain ε c and its associated elastic energy are much larger than their bending counterparts as t → 0. The SI shows that ε c is also large for the observed m values: In Fig. 5, we plot three energy curves for a typical sample.These curves compare the simulated energy for a typical example, U sim elastic (∆) ( lower curve), with two theoretical estimates: U model elastic (∆) ( upper curve; see Eq. 6) calculated from our conetriangle model, and U sim ⟨κ 2 ⟩ (∆) ( middle curve) calculated indirectly from the simulated solution through measurement of the mean squared curvature at the outer boundary: ⟨κ 2 ⟩ out ≡ 1 2π(1+w) ds κ 2 || , and by using Eq. 5.The three curves are similar, differing by a roughly constant factor of order unity.Notably, the full simulated energy including any strain energy, is smaller than the model bending energies.This discrepancy is natural for the circular cone model, since this circular ansatz does not minimise the bending energy.The assumed uniform, circular curvature does not minimize the energy of the cone shape needed in order to bridge between the adjacent triangles.We expect that the system chooses a curvature profile κ(s) that minimizes this energy.Thus, || Here, κ is the arc curvature, without reference to the surface.As before, for fine wrinkles, we have κ ≈ κ N .as expected, the use of the measured κ(s) improves the agreement with the measured energy greatly.But a clear discrepancy remains.
We discuss the possible source of this discrepancy in SI Sec.6.7.It compares the simulated energy profile with that of a conical shape.The mismatch between the two profiles implies that the real surface does not have the straight director lines needed for a developably isometric surface.Such a mismatch requires strain energy.Accordingly, the mismatch should lessen with decreasing thickness.On the other hand this mismatch enables a reduction in bending energy below that of a cone.Thus the SI suggests how a small strain distortion can reduce the total energy below our isometric estimate, as the figure shows.
We conclude that the elastic energy of deformation for the wrinkles encountered in this study is strongly consistent with our claim of an isometric conical deformation.This minimal geometric picture gives a quantitative account of the structure and energy for a given wave number m. Residual discrepancies seem consistent with expected residual stretching effects.Furthermore, the explicit circular-cone-triangle construction gives a workable quantitative estimate of the full elastic energy.This approximation is least accurate when the splay parameter η is large.

A developable wrinkling morphology
As shown above, the Lamé radial contraction system is qualitatively altered by removing the outer tension from that system.Specifically, the purely geometric definition of our deformation results in a faceted shape that is in clear contrast to prior forcebased deformations 2,31,32 .Our work demonstrates a strong resemblance between the observed shape and the exactly isomet-ric cone-triangle shape defined above.For the numerically feasible thicknesses we report, the cone-triangle hypothesis allows us to reproduce structural and energetic features to good accuracy.Still, detailed analysis shows clear-cut departures from this shape that remain unexplained.
Developable approximations such as the cone-triangle model appear to be a novel approach to describing periodic buckled structures.For the annular contraction deformation studied here, this developable picture seems entirely adequate as an initial quantitative account of the shape.This contrasts with other faceted shapes 33,34 such as twisted ribbons 35 , where a comparable understanding of the energy requires explicit consideration of non-developable geometry and non-bending energy.Our conetriangle decomposition also provides an instructive example of a buckled shape in a radially symmetric system that does not have the conventional separable form -i.e. the buckling displacement is not a product of a radial factor and azimuthal factor.
The developability of the faceted wrinkling structure makes it fundamentally geometric, so that it can be created easily on a computer.This provides a way to programmatically construct smooth radially wrinkled patterns in materials using intrinsically flat components, much in the spirit of origami and other designed materials 7 .In fact, a well-known real-world example of such programmed radial wrinkling is the Elizabethan ruff collar.One way of making a ruff is exactly using the inner Lamé boundary conditions: we take a cloth annulus, and manually make pleats in it such that the inner boundary of the annulus lives on a contracted radius.Figs.1e and 1d show, respectively, a ruff collar and its reconstruction using our circular-cone-triangle model; the visual similarity between the two is striking.The conical outer solution seems to capture the ruff wrinkles, in particular their radial straightness and their shape at the outer boundary.
As noted above, accounting for the wrinkles of Fig. 1c is not equivalent to finding a ground-state configuration consistent with the imposed boundary condition.Our deformation was performed quasi-statically, and the wrinkles appeared continuously without abrupt jumps.Thus the system followed a single energy minimum as it deformed.This observed deformation may well have significantly higher energy than the ground state has.

Nature of the thin-sheet limit
The wrinkles investigated here closely approximate the isometric cone-triangle structure, in both shape and energy, as shown above.That is, many numerical examples with a broad range of thicknesses are well described by a model with no dependence on thickness.We take this as strong evidence that, for a given fine wrinkle number m ≫ max{1, 1/w}, the sheet approaches this isometric shape as thickness goes to zero.This leaves open the question of how the system selects this m, which clearly does depend on thickness.Given the present understanding of the energy and shape, the groundwork is laid to explain how this wrinkle number is selected.A separate paper by one of us 26

Approximations and limitations 4.3.1 Circular cone Ansatz
The key modelling assumption used in this paper is the circular Ansatz for the bent cones.As we have shown, this Ansatz serves as a serviceable approximation to the true energy-minimising Elastica solution, while being simple enough to allow for explicit calculation of several key observables, namely, the slope α out and the energy U elastic .Sec. 3 shows that this approximation is accurate only for relatively narrow cones (small η) and large w.This can be partially explained as follows.We expect the circular Ansatz to hold true only for weak Euler buckling."Weak" here means that the angular squeezing of the cones and the resultant slope, are both small.In the SI, we show that this squeezing experienced by a conical sector during deformation increases sharply with decreasing w, such that for our narrowest annuli, it can be as high as 20% for our maximum ∆.At such large values, one no longer expects the weak Euler regime to be valid, and the resultant shape manifests its non-linearity through the flattening apparent in Fig. 4c.We note, however, that even at such large contractions, both the circular Ansatz and the numerical simulations give a linear energy profile, although their slopes are no longer close.

Radial Strain
Our main quantitative evidence for isometric wrinkling comes from the azimuthal strain of Fig. 4b and the total energy of Fig. 5.We were not able to obtain comparable quantitative evidence by measuring the radial strain.The nominal radial strain reported in Fig. 4a is of the same order as the azimuthal strain, but it is negative, contrary to expectations.These small strains are significantly affected by further effects that we did not include.One such effect is the vertical slope of the material radial lines.For example a material line extending radially from the tip of an inner zigzag to the peak of an outer wrinkle undergoes a change of height z from z in to z out .Now, the estimate of Fig. 4a amounts to approximating this line by its horizontal projection, ℓ xy .The full displacement between inner and outer ends of the line is thus For the data point in Fig. 4a at ∆ = 0.2, this effect adds a correction to the strain of +.0002 to +.0003, greatly reducing the reported value of strain.In addition, the material line is necessarily longer than displacement between its inner and outer ends.Further corrections come from localised strain effects near the inner boundary, as discussed in the next section.
Because of such corrections Fig. 4a shows only a rough consistency with the asymptotic isometry we claim here.For stronger evidence we rely on the energy measurements of Fig. 5.

Inner boundary approximation
Despite this strong resemblance between the simulated wrinkles and the cone-triangle structure, the match cannot be exact.Our construction of the cone-triangle shape in Sec 2 required that the inner boundary be a polygon.It is only a good approximation to the simulated circular annulus when the wavelength (i.e.2π/m) is sufficiently small.In particular it must be small in relation to the width w if the triangles are to resemble those of Fig. 1.Details are provided in the SI, Sec.6.5.
One clear disagreement concerns the length of the inner boundary.In the simulation, the undeformed inner boundary is a circle.Under contraction, each tilted circular arc segment must stretch if it is to remain on the surface of a cylinder.Indeed, we observed such stretching along the inner boundary nodes.However, this stretching zone was confined to this first ring of nodes, with a much weaker stretching in the adjacent ring.In addition, we verified that the associated elastic energy was much smaller than that reported in Fig. 5.
The inexactness of the triangle construction at the inner surface makes it clear that the construction is not a realistic representation of a circular annulus contracted onto a cylinder for arbitrary wavenumber m.We demonstrated the isometric shape for only a limited range of geometries; both the displacement ∆ and the width w were confined to moderate ranges.For these moderate conditions the system appears to choose wavelengths fine enough that the inaccuracies of the cone-triangle construction are not important.For more extreme ranges, we expect that the cone-triangle picture shown here will need modification and that the t → 0 limit will be more subtle.In particular, the locally stable isometric structure observed here may give way to macroscopic folding 23 .Finding the conditions in which the conetriangle wrinkles occur is a fruitful subject for future work.

Lateral constraint boundary condition
As noted in Methods, we prohibit any lateral movement of the nodes on the inner boundary, in order to stabilise the faceted wrinkled regime against folding instabilities.However, even if we do allow lateral movement, we see the same faceted wrinkling at early simulation times, before azimuthal movement of the boundary eventually leads to folding 23,36 .That is, the, excess material concentrates at a single sector along the circumference.Our observations suggest that faceting can be seen for initial buckling, even in the absence of lateral constraints.But our evidence is still only suggestive and needs to be confirmed through physical experiments.

Application
The deformation geometry treated here amounts to a strategy for avoiding the typical fate of deformed sheets.Such sheets adopt structures that require significant elastic energy beyond the needed bending energy.Often these structures produce uncontrolled disorder, e.g. as seen in a crumpled sheet.By contrast, the radial contraction used here guides the sheet into a locally lowenergy channel resulting in only bending energy.This low-energy state is moreover attainable by continuous deformation from a flat reference state.Since these deformations lack the additional deformation costs of competing structures, they are intrinsically (locally) stable.They provide a general means to fold and compact a sheet continuously and reversibly with minimal external guidance.We note that the curved inner boundary appears important in stabilising this structure; straight boundaries with no radial contraction 14 do not produce analogous regular wrinkling.
These guiding channels should appear more generally than in the simple circular geometry studied here.The same principles control other inwardly curved boundaries, and other directions of contraction, e.g.onto a constraining cone instead of a cylinder.

Conclusion
While the emergence of developability in wrinkled systems has been studied intensively over the past decade, the simplicity and symmetry of the inner Lamé geometry makes it a particularly clean realisation of such a phenomenon.This faceted morphology is qualitatively different from previously known forms of twodimensional wrinkling.Significantly, it creates sharply dictated, strongly distorted features with negligible stretching and without a separable radial/azimuthal shape.Moreover, it has many variants to be explored, with real world applications.Finally, this work sets the stage to explain the selection of the wavenumber m 26 .

Materials and methods
To investigate the wrinkling morphology of the inner Lamé system, we used the commercial finite-element (FE) solver Abaqus 2018 (SIMULIA, Dassault Systèmes).Using both implicit and explicit dynamic methods ** in a quasi-static regime, we performed simulations of the inner Lamé system using its defining boundary conditions: radial displacement e r (r = 1) = −∆ at the inner boundary mesh points, with the outer boundary nodes (at r = 1 + w) kept free.The maximum value of ∆ applied was ∆ max = 0.267.
The simulation seeks the configuration that would be attained by adiabatically contracting the inner boundary (i.e., its nodes).Accordingly, it gradually decreases e r from zero to −∆ as time t increases from 0 to t f .The value of t f is made sufficiently large that the kinetic energy is always a small fraction of the potential energy.Then the final configuration was not noticeably affected by further increase of t f , as detailed in the SI.
To account for the possibility of high strains, the sheet was initially modelled as a Neo-Hookean hyperelastic material with coefficients equivalent to the linear moduli: Young's modulus, E = 0.907125 MPa, and Poisson ratio, ν = 0.475 (c.f.Eq. 6), corresponding to a rubber-like material.The associated bending stiffness B with thickness t is then given by B = Et 3 12(1−ν 2 ) .To verify that results are independent of the material model, we also reperformed several simulations with a linear material model with these same moduli.Further details on our FEA methods are given in the SI.
To test the validity of our results over a range of parameters, we kept the inner radius fixed and varied the other two parameterswidth w and thickness t -over the range of a decade.The width w of the annulus was varied between w = 0.33 (very narrow) to w = 1.67 (moderately wide) -a factor of almost 10.The thickness t was varied between t = 2.67 × 10 −3 (thick) to t = 1.33 × 10 −4 (very thin), a factor of 20.We performed consistency checks to ensure that the final morphology was independent of the choice of any simulation parameters.For data extraction, we used the ** 'Dynamic' means through time integration of the equations of motion.'Implicit' (using a modified Newton-Raphson method) and 'explicit' refer to the method of integration.
We used two principal deformation protocols to ensure that the wrinkles we observed were characteristic of the quasi-static limit.Our baseline protocol was done using explicit dynamics, in which the finite elements are accelerated with a constant damping factor.The other deformation protocol used an implicit scheme for integrating the forces, with the damping factor automatically selected by the software to favour quasi-static behaviour.In both cases, to avoid gross kinetic effects, we increased ∆ from 0 to ∆ max at a rate slow enough such that the kinetic energy of the system always remained ⪅ 5% of the elastic energy.This is standard procedure for quasi-static analyses in finite-element simulations (see SI for details).
The specifics of the wrinkling morphology depend on how the contracted boundary is constrained.For example, the clamped boundary studied by Mora and Boudaoud 11 leads to a highly compressed boundary layer whose relation to the sheet is difficult to discern.The other extreme is to allow each point on the boundary circle to lie anywhere on the confining cylinder; such boundaries lead to wrinkling only in a transient regime, leading ultimately to collapse into a fold 23,38 .Here we study an intermediate constraint in which points on the bounding circle may move only axially (i.e.vertically) on the bounding cylinder; azimuthal (θ ) motion is not allowed.Doing so automatically prohibits folding (since this requires lateral motion), but without interfering with the wrinkling.
6 Supplementary Information

Construction of circular conical sections
In the main text, we describe how we divide the flat annulus into two regions -an inner region composed of isosceles triangles, and an outer region composed of arc sectors.The main text describes how the inner triangles can isometrically deform under an inner displacement (or contraction) of length ∆ by tilting about the normal to the inner boundary -this gives an inner 'buckled' solution.
Here, we describe how to construct the outer solution whereby the arc sectors isometrically bend into sections of right-circular cones.In other words, each arc bends into a circular profile.The requisite variables for our geometric construction are described in Fig. 6.In the flat state, each arc sector, lying between two adjacent flat triangles, subtends an angle φ f at the common vertex of these two triangles.If the edges of the triangles have length R f , then the length of the largest arc is S = R f φ f .On tilting, the edges of these same triangles approach each other, so that they subtend a reduced angle φ t (see SI Fig. 7a).This reduction in angle (∆φ = φ t − φ f ≤ 0) constitutes the effective compressive strain ε cone = ∆φ /φ f for the arc sector.SI Fig. 7a shows the evolution of ε cone with boundary contraction ∆.We see clearly that narrower sheets feel a greater ε cone for the same ∆.
The excess material can be accommodated by bending the sector into a cone.The exact shape of this cone can be determined by solving the Elastica equation with the appropriate boundary conditions.However, for simplicity, in this paper we choose to approximate the exact shape by a circular cone.The aim then is to find the circular cone that fits in between the two adjacent triangles while exactly accommodating the excess angle ∆φ .Fig. 6b shows all the variables needed for this operation.The isometrically bent conical section has to fit inside the triangle △AOB (drawn in teal) defined by the apex 0 and the tilted edges OA and OB (see also Fig. 3 in main text).The conical section itself is drawn in bold black.Let the circular cone have angular extent φ c , and let its maximum radius be ρ c , centred at point C. Then the conical shape can be determined from the following set of geometric constraints.
The first constraint is that of inextensibility, i.e. length conservation.This gives us: The second constraint involves specifying the end-to-end distance S 1 between points A and B. Thus, using the blue triangle △AOB and the grey triangle △ACB in Fig. 6b, we have: Eqns. ( 8) and ( 9) constitute a set of simultaneous equations for the cone variables ρ c and φ c , since R f and φ f are fixed by the flat geometry, while φ t is set by the tilt of the neighbouring angles.Finally, using basic trigonometry, we have for the cone's apex angle: In sum, we have three independent equations for three unknowns: φ c , ρ c and β (all positive-definite), which fully specify the cone in space (including the centre C).SI Fig. 7b shows a representative solution of such a constructed cone.SI Fig. 7c shows the evolution of this construction with increasing ∆.

Correcting for scalloped arc sectors
The flat arc sectors as defined above -circular arcs with their origin at the inner boundary -lead to the creation of a flat shape that is in fact a scalloped annulus, whose outer circumference 2mR f φ f is greater than the expected 2π(1 + w).
We can correct for this length discrepancy when constructing the bent cones, by changing R f → χR f in Eq. ( 8) where 9) remains unchanged since it specifies the end-to-end-distance, which is set by the neighbouring tilted triangles.We note that this factor χ is smallest for samples with small w and m (wide splay η), and becomes ≈ 1 for large w and m (narrow splay).

Calculating the conical strain ε cone (∆)
The inset of Fig. 7a shows the evolution of the conical strain ε cone with contraction ∆ for our cone-triangle construction.The observations (plotted with dots) show that ε cone (∆) is linear, with the slope depending significantly on the width w of the sheet, but independent of the wavenumber m.To understand this behaviour, we analytically derive an expression for ε cone (∆) = (φ t − φ f )/φ f .Since φ f is fixed by the initial configuration, to find ε cone , we only need to calculate φ t (∆).The angle φ t (∆) can be calculated given any one of the inner triangles, the flat arc angle φ f , and the tilt angle α.We use the triangle △OPQ from Fig. 7a, located in the x − y plane.The annular width w gives the height of the triangle, while h gives its base.Now, instead of rotating about O, we choose to rotate the triangle about its mid-line.For this, we let X be the mid-point of OP, and we set it to be the origin (0, 0), as shown in Fig. 8.In similar fashion, we define Z to be the  midpoint of the edge QR, so that the bisector of the cone is the line OZ, shown as a dotted blue line in Fig. 8.The angle between OQ and the blue dotted line is then φ f /2.This triangle is then tilted clockwise out-of-plane about − → XQ (parallel to the y-axis).Under this rotation, P → P ′ and O → O ′ , while X and Q remain unchanged.In this tilted configuration, φ t /2 is given by the angle between the rotated vector In what follows, we use vector notation to denote both 3d and 2d vectors, precising their nature where necessary.The steps for determining φ t are as follows: 1. Find the 3-vector where R α is the standard 3x3 rotation matrix of angle α about the y-axis (positive α is considered anti-clockwise), and − → XO = (h/2, 0, 0) .

Then
3. We only need its orthogonal projection on to the x-y plane: 4. We need the vector angle between this and the vertical blue plane defined by − → OZ.Here, Calculating this in a symbolic software like Mathematica, we get a complicated trigonometric expression that depends on three variables: α, φ f , and the dimensionless triangle aspect ratio h/w.
To simplify this expression, we consider the regime of large m, which corresponds to both small h/w and small φ f .Expanding to first order in h/w, and using the exact relation cos α = 1 − ∆ to replace α by ∆, we get: Thus, we get for the cone contraction factor: Substituting h/w ≈ π mw and φ f ≈ π(1+w) mw (valid for large m), we get: Eq. 13 predicts that, for large m, the slope of ε cone vs. ∆ should be given by the ratio of inner to outer radius of the annulus.We plot this prediction as solid lines alongside the measurements in the inset of SI Fig. 7a.The two are fully consistent.While Eq. 13 for the slope has been shown only in the regime of narrow cones, in practice, we find it to be approximately valid even for wider cones.The qualitative conclusion here is that, for given ∆, wider annuli effectively feel less squeezed than narrower annuli in a cone-triangle deformation, independent of the number of triangles/cones.

Calculating bending energy for a cone
In this section, we continue with the variables introduced in Sec.6.1.A right-circular cone is made up of a series of circles whose radii ρ increase linearly with distance (say, ζ ) along the cone axis, measured from the cone tip.For a cone of axial length L c (see SI Fig. 7a) and maximum radius ρ c , we thus have: We can also write ρ c L c = tan β , where β is the cone's vertex angle.Now, consider an infinitesimally wide circular band on the cone, of radius ρ, angular extent φ c , and width dζ cos β .The band has bending energy dU = B 2 × κ 2 N × (φ c ρ) × dζ cos β , where B is the bending modulus and κ N is the local normal (i.e.out-of-plane) curvature.Given the geometry of the cone, we have κ N = κ cos β , where κ = 1/ρ is the arc curvature of the circular band.Thus dU = B 2 × (cos β /ρ) 2 × (φ c ρ) × dζ cos β .Simplifying and integrating over an entire conical sector gives: The singularity as ζ → 0 means that the conical shape must be modified there.Thus we consider only the region beyond some L core .
where β = tan −1 ρ c L c .For the majority of numerical samples discussed in this paper, where m ≫ max{1, 1/w}, it is sufficient to use small-angle approximations for β .Thus, setting cos β ≈ 1 and sin β ≈ β , Eq. ( 16) gets reduced to : Here, φ c , ρ c and L c are all functions of the contraction ∆.But φ c and ρ c are related through the constraint: (see Eq. 8), and L c ≈ w (the width of the annulus), so U cone bend can be reduced to a function of the single dynamic variable ρ c (∆).For the entire annulus, we need to multiply this by the number of cones 2m.Thus, we get: More generally, for a cone extending between axial limits L min and L max , we have: This is the expression given in the main text.While Eqs. ( 18) and ( 19) seem independent of m, it is not so.Eq. ( 7) in the main text shows that ρ c ∼ (1 + w)/m =⇒ U theory bend ∼ m 2 , as expected of a bending energy.We note that this ρ c ∼ 1/m scaling could have been predicted in another way.In the limit of maximum possible contraction (∆ → 1), the point C approaches the xy-plane, and so the diameters of the 2m circles must approximately equal the reduced projected outer perimeter: 2π(1 + w − ∆) → 2πw.This gives us: 4mρ c → 2πw =⇒ ρ c → πw 2m (the dotted line in the inset of SI Fig. 7c).Indeed, even if C is off the x-y plane (e.g. for smaller ∆), the diameters of the 2m circles must still equal the reduced perimeter 2π(1 + w − ∆) up to some factor.Thus, we have the scaling relation ρ c ∼ 1/m as expected.
Finally, we note that the data presented in the paper represents the full expression 16, without any approximation.

Limits of the cone-triangle construction
To the best of our knowledge, the above conical construction works as long as the initial (flat) angle φ f < π, i.e. as long as the edges of two adjacent triangles define a triangle.The value of φ f depends on the width w and the wavenumber m, and increases as m decreases.This defines a minimum wavenumber m cone min (w) below which a conical solution is invalid.For a sufficiently wide annulus, we find that m cone min = 2, which is the minimum possible value for any wrinkled solution.However, for very narrow annuli, this value goes up.Thus, for w = 0.2, we find m cone min = 4. Fig. 9 shows two such contrasting geometries.m cone min (w) defines a geometric limit beyond which we expect our cone-triangle model to fail.However, the main text shows that we already see significant deviations from our model for Abaqus solutions with wavenumber m significantly higher than m cone min (w).

Comparison with separable, sinusoidal wrinkling
The cone-triangle shape described in the main text sharply reduces the elastic strain and energy relative to conventional wrinkling with height profile h c (r, θ ) e.g. of the form 17 as we now illustrate.
In the far-from-threshold regime, this h c profile must relax the azimuthal strain ε θ θ arising from the inward displacement u r = −∆.In general this azimuthal strain has the form 2 : For purposes of comparison we may evaluate this ε at a zero of h c (r, θ ) e.g.mθ = π/2.There by symmetry the azimuthal displacement u θ = 0 and ∂ θ u θ = 0. Thus the second term in Eq.21 vanishes.In the last term ∂ θ h c = ±m f (r).We use this expression to estimate ε θ θ .Choosing h c (r, θ ) to make ε θ θ vanish implies So that f (r) ≃ 1 m √ 2r∆ This non-constant f (r) entails a radial strain ε rr given by 2 Here u r is constant as noted above and ∂ r h ≃ ∂ r f (r) ≃ 1 m 1 2 2∆/r, giving an estimated conventional wrinkle strain ε c rr of This radial strain of conventional wrinkling is to be compared to the bending strain ε B in the cone-triangle model.This strain arises from the bending stress σ B θ θ where bending modulus B is related to the thickness t, the Poisson ratio ν and the bulk Young's modulus E by 39 Sec.12 Using ε = (Et)σ , then Eq. 25 yields For the annuli simulated above the conventional wrinkles have much greater elastic strain than the bending strain we report: their ratio is given by Evidently for fixed m the ratio diverges as t → 0. For the specific annulus used for Fig 5 with r = 1, ∆ = 0.2, ε c /ε B ≃ 50.By these estimates the sinusoidal wrinkling carries substantially larger strain and thence elastic energy relative to the isometric wrinkling we report.

Departures from conical shape
Here we gauge the impact of the stretched zigzag vertices on the total energy by examining how the elastic energy varies with distance r − 1 from the inner rim.Fig. 10a shows this energy for the right-most data points of Fig. 5.The red simulation point (lower curve) in Fig. 5 is found by summing the energies in each azimuthal ring of finite elements over the range of r indicated by the hashed region.The simulated ring energies for selected r are shown as black rectangles.
The corresponding cone-triangle energy, i.e. the right-most point on the (yellow) middle curve in Fig. 5, is found using the cone ring energy plotted as a solid curve, using Eq. 5.This equation implies a ring energy varying as 1/(r − 1).Since the yellow point is found by matching the simulated curvature at the outer boundary, the cone model gives a ring energy that matches the simulated ring energy there, as Fig. 10a shows.
By comparing the simulated ring energies with the cone model, we can gain insight into the difference in energies seen in Fig. 5.The inner simulation point shows the expected large discrepancy with the model, which unrealistically diverges at the inner boundary.The middle simulation point is at the inner boundary of the region treated in Fig. 5.This ring energy is 40 percent smaller than the cone model prediction.This difference is consistent with the 15 percent differences between the energies of Fig. 5. Any cone that would cure this discrepancy would have to extrapolate to a vertex beyond the inner boundary.
The azimuthal strain profile shown in Fig. 10b also shows departures from the model.The plot shows the azimuthal dependence of the azimuthal strain ε θ θ for a radial position r equidistant between the inner and outer boundaries..The plotted strain is maximal at peaks and troughs of the wrinkles.Though the averaged strains,are consistent with bending strains, as discussed in Sec.3.1, there are strong oscillations adjacent to these peaks and troughs with a period of two finite elements.This suggests that our simulations have limited reliability for predicting these weak strains near the peaks or troughs.Away from the peaks and troughs the strains vary smoothly.Of special interest is the point where the model cone and its adjacent triangle would meet.These points, marked by dashed lines show no sign of discontinuity.
These detailed features of Fig. 10 show shortcomings of the cone-triangle model.However they underscore the relevance of this model for understanding this form of buckling.

Finite-element method (FEM) simulation details
For our simulations, we used the commercial finite-element software Abaqus 2018 (Simulia, Dassault-Systèmes, Providence, RI).This section describes the different steps for generating a typical simulation of our inner Lamé system, in the order typical of a finite-element software.
The assembly consisted of only a single annulus, with inner radius fixed and taken to be unity, and with varying width w and thickness t in order to test our system over a wide range of system parameters.For width, we used values w = 0.20, 0.33, 0.67, 1.0, 1.67 (a factor of almost 10, ranging from very narrow to moderately wide), and for thickness, we used values t = 2.67 × 10 −3 , 1.33 × 10 −3 , 6.67 × 10 −4 , 2.67 × 10 −4 , 1.33 × 10 −4 (a factor of 20, ranging from moderately thick to very thin).While these thickness values vary over a decade, the values still fall well within the thin sheet limit.The annular part was made of 2d shell quad (S4R) elements 40 .This choice was made mainly to optimise speed, since we used a fine enough mesh to ensure that doubling the linear mesh size change the energy by a negligible amount (⪅ 1 − 2%).For comparison, the coarsest mesh we used was for the w = 1.67 annulus, with 60 elements across the radius and 1400 elements across a circle, giving a maximum linear size for an element ≈ 0.01.For consistency checks however, we also ran some simulations with annuli made of 3d volume cubic (C3D8R) elements 40 , which gave the same morphology (with the same wavenumber), but which much longer running times.
When discussing the material properties, for concreteness, we will use SI units (and thus take the inner radius to be 1 m).For the material properties, we mostly used a standard neo-Hookean hyperelastic model 40  The elastic modulus can be used, along with the material density ρ and the average linear mesh size l e , to determine the average integration time scale (i.e. the 'stable time increment') in the simulation, as follows.The elastic bulk modulus K and the density ρ determine the speed of sound in the material, c s = K/ρ † † .The stable time increment in the simulation is on the scale of the time required for elastic information to traverse an average mesh element: ∆t ∼ l e /c s .For our values of K ∼ 10 5 Pa, ρ = 10 3 Kg/m 3 and l e ∼ 10 −2 m, we get ∆t ∼ 10 −3 s.This in turn determines the dimensionless number of iterations n iter performed by the solver in a simulation running over time period T : n iter = T /∆t.Below, we discuss the typical values of T used in our simulations, and how increasing T allows us to reach a quasi-static limit in dynamic integration methods (i.e.where kinetic energy is present but negligible compare to elastic energy).
For the radial displacement loading at the inner boundary, we applied velocity and displacement boundary conditions (BCs) interchangeably.Typically, we applied velocity BCs with a linear amplitude profile, and displacement BCs with a smooth-step profile 40 , in order to assure a smooth (i.e.zero velocity) pull at the † † Alternately, one can use the Young's modulus E instead of the bulk modulus K in the definition of c s .But this does not qualitatively alter our argument above.
beginning.These were applied so that the maximum displacement amplitude ∆ max = 0.267 is attained within a time period T (defined in the same units as the ∆t given above).The value of T was chosen to be large enough to ensure small kinetic energy and give a T -independent configuration, as defined below.Typically, we used T = 20 for the thicker sheets, and T = 80 for the thinnest sheets.
For the simulation protocol, we employed both 'dynamic explicit' and 'dynamic implicit' 40 integration schemes in the quasistatic limit (as compared to a fully 'static' energy minimisation scheme).The word 'dynamic' refers to the presence of inertia, while 'explicit' and 'implicit' refer to the solution scheme.'Explicit' means explicit time-integration of Newton's second law, while 'implicit' refers to implicit integration (viz.through iterative root-finding) of Newton's law, using a modified Newton-Raphson method.The mixture of these two methods was done partly as a consistency check, partly for convenience, and partly by necessity.While the implicit method in the quasi-static limit was faster for most jobs, the explicit solver was indispensable for the thinnest samples, where the static solver ran into convergence problems.Ensuring the quasi-static limit is also easier for 'dynamic implicit' than for 'dynamic explicit'.In dynamic implicit, the quasi-static option is in-built, but for dynamic explicit, it has to be ensured manually by applying the loading slowly enough so that further slowing has no effect on the final shape and energy.
For this, we need to look at the available energy modes.The energy balance equation in Abaqus 40 is given by (ignoring possible terms coming from viscosity, friction, heat, contact and constraint penalties, etc.), is: where E I is the internal energy, E KE is the kinetic energy, and E W is the work done by externally applied loads.For us, the internal energy is just the elastic energy (by design, there are no other energy modes).Thus, for quasi-static loading, one generally requires the kinetic energy (E KE ) to be < 10% of the elastic energy (E I ).In practice, we kept the ratio to ⪅ 5%.Since the elastic energy is thickness-dependent (always increasing with increasing thickness), thinner sheets required slower applications of the loading.E.g., for the thinnest samples (t = 1.33 × 10 −4 ), this involved applying the contraction ∆ over a time period T = 80−120, viz.using n iter ∼ 10 5 solver iterations.As a result, the slowest simulations, for the thinnest and widest sheets, lasted ≈ 150 corehours; the average simulation however, lasted between ≈ 30 corehours.For reference, the validation case (described blow) of folding a flat sheet into a cylinder, albeit with a much coarser mesh, was accomplished using ≈ 0.2 core-hours.
For data extraction, we used field output for the displacement variables, and history output for the energy 40 .For the elastic energy, we used the "ALLIE" (internal energy) variable 40 , equivalent to the E I variable in Eq. 29.Since the simulations are done using a dynamic time-integration scheme with inertia, there is an inherent noise in the energy values arising from imprecision in finding the exact energy minimum (viz.due to inertial oscillations).We cannot estimate this noise precisely, but a rule-ofthumb estimate is ⪅ 5%, i.e. of the same order of magnitude as the ratio E KE /E I .However, in reality, it might well be less.Significantly, this noise cannot account for the discrepancy between energy measurements and model predictions in Fig. of the main text.Finally, post-processing was done using Abaqus2Matlab 37 .

Testing the numerics for known cases
We used the above procedures (albeit with a static energy minimisation scheme) to calculate an analytically solvable case, to verify that the shape and energy agreed with the known results.The example was a rectangular sheet of width w = 1, length L = 2π, and thickness t = 1 × 10 −3 , in which we prescribed boundary conditions on position and orientation of the short edges to make them curve up and in, to form a circular cylinder of unit radius.
We verified the circularity of the cross-section by projecting onto the plane and measuring the distance from the centre.We found that no point differed in its axial distance by more than .001%.The measured elastic energy differed only slightly from the analytic result, U cylinder bend = Bπ ≈ 3.07 × 10 −4 Joule, where B = Et 3 /12(1 − ν 2 ) is the bending modulus, obtained using the thickness t, Young modulus E and Poisson ratio ν quoted above.The simulation gave an energy ≈ .01%larger than this.A discrepancy of this sign is expected because the analytic form neglects the small strain energy owing to the nonzero thickness of the sheet simulated. -----------------

Fig. 2 a
Fig. 2 a) Sketch of the inner Lamé deformation of a horizontal annulus.(Top) Top view showing the undeformed annulus, of inner radius 1 and width w. (Bottom) Vertical cross-section of the deformed annulus through a line where the vertical displacement of the wrinkle is maximal.The undeformed state is shown in light shading; the dark lines are a section of the deformed annulus, contracted horizontally by a distance ∆.These lines are also displaced upward and tilted downward.b) Sketch of the circular-cone-triangle construction in its initial, flat state.Here we choose wavenumber m = 6, hence it has 12 triangles (shown in colour).The mid-line and height w of one triangle are indicated.The sectors bridging the triangles are shown in white; S is its outermost arc length.When the bases of the triangles are moved inward, they must tilt to retain their shape and connectivity.The bridging sectors are squeezed laterally to form cones.A typical deformed example is shown in Fig.1c.c) Top view of a narrow sector of the sheet, with detail of an area element at the outer edge of Lamé sheet showing radial (horizontal) and azimuthal (almost vertical) boundary forces.In the fully buckled (FT) regime 2 , the azimuthal forces nearly vanish (as indicated by the upper and lower X's).In the the present "inner Lamé" system, the outer radial force also vanishes, indicated by the X at right.Thus, the inward radial force must also nearly vanish for force equilibrium.All elements to the left of the pictured one are subject to the same argument.Thus except for bending stresses, all radial stresses vanish.

Fig. 4
Fig.4Comparisons of geometric quantities between theory and simulations (using the explicitly periodic initial conditions described at the beginning of Sec. 3).(a) Testing radial inextensibility.The graph shows an average radial strain inferred from Eq. 1, i.e. the small relative difference between the inward displacement at the outer boundary and that at the inner boundary, further discussed in Sec.4.3.2This example had w = 0.33, m = 25, and t = 2.67 × 10 −4 .(inset) Sketches of the wrinkle profile at the two boundaries, showing the definition of lengths L 0 and L xy , and slope angles α out and α in .The blue arrows indicate the compression direction.(b) Testing azimuthal inextensibility at the outer boundary.The magnitude of the residual azimuthal strain |ε resid θ θ | recorded slightly inside r = 1 + w, plotted logarithmically against the collective variable m 2 t 2 /(1 + w) 2 at large ∆, for multiple values of m, t and w.The black line shows the linear fit: |ε resid θ θ | ≈ 0.24 m 2 t 2 (1+w) 2 .(c) Wrinkle slope angle vs. applied contraction ∆ ≡ ∆/r at the inner and outer boundaries (see Sec. 3.2).(Top row) We separately test Eq.(3) for α in (in red) and Eq.(4) for α out (in black).The simulation data points are plotted alongside the model predictions.Simulated samples of differing width w and wavenumber m are coloured using the dimensionless cone splay parameter, η ≡ π(1+w) mw .The error bars are standard deviations from averaging over all the wrinkles.Markers: w = 0.33 (circles), w = 0.67 (squares), w = 1.0 (diamonds), w = 1.67 (triangles).(Bottom row) Example wrinkle profiles for small and large η: (left) for η = 0.25 (w = 1.0, m = 25), and (right) for η = 1.8 (w = 0.33, m = 7).Colouring denotes height.
addresses this question.J o u r n a l N a me , [ y e a r ] , [ v o l .] , 1-16 | 7

Fig. 6
Fig. 6 Isometrically constructing a sector of a cone from a flat arc sector.(a) The flat configuration.(b) The curved configuration; the actual conical sector and its variables are depicted in black, the triangle △AOB formed by the edges of the two adjacent tilted triangles, along with its variables, is in blue, and the putative cone with its variables is in grey.

Fig. 7
Fig. 7 Constructing a single cone (a) Two adjacent flat triangles (in yellow) subtending an angle φ f tilt towards each other by angle α.The tilted triangles (in deep blue) subtend a smaller angle φ t .Thus, the region between the triangles (shown approximately as a violet triangle) gets squeezed by an angle ∆φ = φ t − φ f .(inset) Evolution of the conical strain ε cone = (∆φ )/φ f with ∆ for three different widths w.Model measurements are plotted alongside the predicted lines of slope = −1/(1 + w) (see Eq. 13).(b) To fit within the contracted boundary defined by OA and OB, the region must bend into a circular cone of radius ρ c , angular extent φ c , axial length L c and tilt angle β .The only ingredients needed for the construction are △AOB and the tilt angle α(∆).(inset) Cross-section of the cone (in magenta) through △OAC, showing that the arc curvature κ and the normal curvature κ N are related by κ N = κ cos β .(c) The evolution of the conical solution with ∆ (here, 0.05 ≤ ∆ ≤ 0.6).Greater opacity means larger ∆; the cones are translated with respect to each other for viewing clarity.(inset) As ∆ increases, the radius of curvature ρ c decreases, and for ∆ → 1, attains an asymptotic value of ρ c ≈ πw 2m (dotted line).

Fig. 8
Fig.8Diagram supporting the calculation in Sec.6.3.Triangle OPQ is the same as in Fig.7a.Here, we give the points explicit coordinates ( shown in black text).For more details, see main text.

Fig. 9
Fig. 9 Two different flat state geometries for w = 0.2.(Left) For m = 4, the flat angle φ f < π, which allows the conical construction described in Sec.6.5.(Right) For m = 2, the flat angle φ f > π, which means that our conical construction is not valid here.

Fig. 10
Fig. 10 a) Energy distributions for the annulus of Fig. 5. Horizontal axis is radial distance from inner rim measured in finite element widths.Vertical axis is ring energy described in the text.Hashed region shows the region treated in Fig. 5. Its area represents the simulated energy with largest ∆ on the lower curve of Fig. 5. Black marks show simulated ring energies at three selected radii.Solid curve shows the ring energy profile calculated from cone-triangle model used for the middle curve in Fig. 5. b) Simulated azimuthal strain ε θ θ (in red) as a function of azimuthal finite-element coordinate θ , for a material circle located at the middle of an annulus with (w = 0.33, t = 1.33 × 10 −3 , m = 10, ∆ = 0.27).Light curve in blue shows height profile.Region spanning one trough and peak are shown.Variability of other troughs and peaks is similar.Vertical black lines indicate the location of two consecutive cone-triangle boundaries.