Wrinkling composite sheets

Abstract

We examine the buckling shape and critical compression of confined inhomogeneous composite sheets lying on a liquid foundation. The buckling modes are controlled by the bending stiffness of the sheet, the density of the substrate, and the size and the spatially dependent elastic coefficients of the sheet. We solve the beam equation describing the mechanical equilibrium of a sheet when its bending stiffness varies parallel to the direction of confinement. The case of a homogeneous bending stiffness exhibits a degeneracy of wrinkled states for certain lengths of the confined sheet; we explain this degeneracy using an asymptotic analysis valid for long sheets, and show that it corresponds to the switching of the sheet between symmetric and antisymmetric buckling modes. This degeneracy disappears for spatially dependent elastic coefficients. Medium length sheets buckle similarly to their homogeneous counterparts, whereas the wrinkled states in large length sheets concentrate the bending energy towards the soft regions of the sheet.

---

1 Introduction

From surfactant monolayers to tectonic plates, the deformation of elastic sheets underlies a vast landscape of problems in science and engineering. Not only do their patterns evoke a great aesthetic appeal, but they underlie both function and form in the world around us.^1–4

The buckling of a homogeneous sheet on top of a supporting foundation has been the object of soft matter studies for decades see e.g. ref. 5–7 and references therein. In a free-standing confined sheet (Euler's elastica), the compressive load (which we refer to as “compression” throughout the paper) at which buckling occurs, and the scale of this buckling, are set by the length L* of the confined sheet. However, in the presence of a supporting foundation, an intrinsic length scale [small script l]* is set by the mismatch between the elasticity of the sheet and that of the foundation. When a resting sheet with clamped ends is confined longitudinally, a uniaxial compression builds up and wrinkles of the surface emerge.⁸ Such wrinkle patterns have been observed in, among many other settings, thin polymer sheets resting on the surface of water,⁹ bacterial biofilms,¹⁰ granular rafts,¹¹ the mineral veins of rocks,³ strained epitaxial films,¹² and in glaciology.^13,14 From a theoretical standpoint, wrinkles of confined elastic sheets are predicted by a simple set of geometric rules.¹⁵ The theory for supported sheets also maps onto a model to describe the elasticity of an unsupported epithelial monolayer.¹⁶

In many settings the deforming material is treated as a homogeneous solid, despite a potentially important intrinsic composite structure. Indeed, whether or not the composite structure can be treated as effectively homogeneous requires appropriate quantitative analysis. When the host material is stiff, such as ice, and the inclusions are soft, the theory of Eshelby¹⁷ shows that the effective Young's modulus, and hence the bending stiffness is reduced relative to the case of the stiff host material alone. When the host solid is soft and the inclusions are uniformly distributed,¹⁸ two possibilities exist—composite softening or stiffening—depending on the size of inclusions relative to the elastocapillary length, , where γ* is the surface tension of the inclusion-host interface and is the Young's modulus of the host material. However, by controlling the distribution of inclusions we can, for example, prepare a gradient of bending stiffness, and thereby manipulate the failure properties of such elastic composites.

Here we examine the deformation of composite elastic sheets floating on a liquid foundation. We treat the case with a gradient of bending stiffness parallel to the direction of confinement. This gradient is envisaged to be created by controlling the distribution of liquid inclusions.

First, we present key results for the critical compressive load and wrinkle patterns in homogeneous floating sheets. The interplay between the intrinsic length scale, [small script l]*, and the undeformed size of the confined sheet, L*, determines the wrinkled state that has the smallest compression. The observed wrinkles are then either symmetric or antisymmetric about the centre of the sheet. Importantly, at certain sheet lengths the wrinkled states are degenerate: both symmetric and antisymmetric modes exist at the same compression.

While the homogeneous case provides important insights into the wrinkle pattern, our main interest is the study of inhomogeneous sheets. Taking advantage of the effective medium behaviour of composites,^17–19 we examine the wrinkles in composite sheets.

We impose a spatial distribution of liquid inclusions that translates into spatially varying elastic moduli (Young's modulus and Poisson ratio) and hence a spatially varying bending stiffness. The inhomogeneous stiffness is responsible for breaking the symmetry of the wrinkle patterns and eliminating the degeneracies. For sheets large compared to the intrinsic length scale, this gradient shifts the position of the maximum amplitude of wrinkles. This displacement happens at the onset of wrinkling, and hence it does not require any non-linear behaviour. Inhomogeneous sheets are technologically appealing for they offer a new means to control fracturing via the spatial variation of the wrinkle amplitude.

2 Mechanics of a floating composite sheet

We start with a presentation of the governing equations for a composite thin sheet floating on a liquid, including the effects of in-plane forces. For more detail the reader is referred to the book by Mansfield.⁴

2.1 Problem formulation

We denote all dimensional quantities with a superscript (·)* and let x*,y*,z* be the Cartesian coordinates in the horizontal direction, into the page, and in the vertical direction, respectively. We consider a sheet with thickness h* and longitudinal undeformed length L*, such that h* ≪ L*. The sheet (of density ) rests on a liquid of density ρ*, and in the absence of all forces apart from gravity, it floats with its mid-line at a height above the free surface. We measure all vertical displacements relative to this equilibrium level,† as in Fig. 1 where the origin of the vertical axis is set at the liquid level. The sheet is confined lengthwise by a distance d*, and has clamped edges located at x* = ±a* in the longitudinal direction, with a* ≡ (L* − d*)/2. However, we shall impose boundary conditions at x* = ±L*/2, which are the Lagrangian coordinates of the sheet's edges; this avoids confusion with their position at the onset of buckling, x* = ±a*, and is correct within small deformation elasticity. A quantity of interest from our analysis is the value of a* at the onset of buckling for a given sheet length L*.

Fig. 1 Schematic diagram showing a thin elastic sheet floating on a liquid with a lateral compressive force τ* due to confinement.

The geometric incompatibility between the sheet's undeformed length and the confinement length yields a mechanical instability producing a vertical displacement w*(x*,y*). This displacement causes, and is influenced by, the restoring pressure p* = −ρ*g*w* that the liquid foundation exerts on the sheet. An in-plane compression τ* (with dimensions of force per unit length) is calculated as an emergent property, rather than being imposed.

The static equilibrium of a confined floating sheet is determined by the balance between elastic forces and the hydrostatic pressure. In this paper, we are concerned with the effects of a variable bending stiffness, or flexural rigidity, B*(x*), which may be achieved using composite materials whose elastic properties are inhomogeneous. Such inhomogeneous sheets can be made by embedding one material in the other to control the resulting composite structure, as illustrated schematically in Fig. 2. Examples include ionic liquid droplets in silicone,¹⁸ hydrogel particles inside elastomeric matrices,²⁰ 3D printing materials using silicone double networks²¹ and fibre–silicone mixtures,²² or hydrogel substrates with a photo-sensitive cross-linker and a gradient photo-mask.²³ Despite a varying composition, we neglect any variation of the density of the sheet. This is justified for stiff composites in the dilute regime, and for soft composites, because the matrix and the inclusions have approximately the same density. (A varying bending stiffness would also correspond to a single component material sheet with varying thickness, but we do not treat this case here because it would complicate the free-floating equilibrium.)

Fig. 2 Strip of a composite sheet of length L, with liquid inclusions and the linear volume fraction profile of eqn (23) in the x direction, parallel to the confinement.

In the absence of inclusions, namely for a homogeneous elastic sheet, the elastic coefficients (Young's modulus, , and Poisson ratio, ν₀; we denote the elastic constants of the homogeneous sheet with subscript “0”) are constant, and hence so too is the bending stiffness; in particular, the baseline bending stiffness is

The intrinsic length scale of the displacements, or wrinkles,plays a central role in the problem. Therefore, we use [small script l]* to rescale all lengths, while we rescale the compressive force τ* with . Dimensionless quantities are unstarred so that x = x*/[small script l]*, and . This non-dimensionalization also introduces a typical scale for the size of the stresses induced by compression, which we define aswhich is also a measure of the relative ease of bending to stretching of the sheet, or the stretching-stiffness. Therefore, [scr S, script letter S] is important in determining how much confinement is required to induce buckling, as we describe later. Because so that , for thin sheets we expect that [scr S, script letter S] ≫ 1, and hence deformation can be accommodated more easily by bending than by stretching. For example, a 10 cm thick sheet of ice, for which and ν₀ ≈ 0.3, has [scr S, script letter S] ∼ 1000.² A soft material, such as PDMS, for which and ν₀ ≈ 0.5, has [scr S, script letter S] ∼ 100 for h* on the order of cm. However, thin sheets do not always bend more easily than they stretch. Provided that [scr S, script letter S] is finite, sheets may accommodate some of the imposed deformation by compressing (negative stretching), particularly when the confinement is small and the supported sheet resists buckling. We therefore refer to [scr S, script letter S] as the ‘inextensibility’ of the sheet. Before returning to this idea, we introduce buckling in two dimensions where the equations of mechanical equilibrium simplify.

2.2 Two-dimensional buckling of a sheet

Our focus here is on two-dimensional deformations of the sheet under uniaxial compression. A simple argument (see Appendix A) shows that the in-plane stress τ is constant and that small out-of-plane displacements of the sheet satisfywhich is Euler's linearized elastica equation (see e.g. Section 20 of ref. 24) with a lateral load due to the hydrostatic pressure in the liquid foundation and varying elastic properties along the axis of confinement. This is analogous to the small-deflection equilibrium of a compressed inhomogeneous column.²⁵ The boundaries of the thin sheet at x = ±L/2 are clamped so that w_,x(±L/2) = w(±L/2) = 0.

The elastic response of a composite is characterized by the elastocapillary length l* discussed in the Introduction (Section 1). When inclusions of radius R* are smaller (larger) than l*, the composite is difficult (easy) to deform because surface tension can maintain (cannot maintain) their sphericity. For example, for a typical liquid inclusion with γ* = O(0.01 N m⁻¹),¹⁸ and soft materials with a range of E* = O(1 − 100 kPa), l* ranges from 1 to 100 μm. Thus, surface tension effects are important for micron-sized inclusions. A summary of different theoretical models that predict the bulk mechanical properties of a composite is given in Appendix B.

3 Buckling of a homogeneous sheet

The force balance of an elastic sheet with homogeneous properties, namely B(x) = 1 in eqn (4), is given by

(w₀)_,4x + τ₀(w₀)_,2x + w₀ = 0, (5)

where we use subscript “0” to denote the displacement in the homogeneous case. This is to be solved with boundary conditions corresponding to a clamped sheet:

w₀(±L/2) = 0, (6a)

and

(w₀)_,x|_x=±L/2 = 0. (6b)

3.1 Buckling profiles

Equations (5) and (6) define an eigenvalue problem whose detailed solution is given in Appendix C. As expected from the reflection symmetry about x = 0, there are two distinct types of solutions with opposite parity: even functions of the form w^(s)₀(x) = A_s cos kx, which correspond to symmetric buckling profiles, and the odd counterpart w^(a)₀(x) = A_a sin kx, giving antisymmetric profiles. For both cases the wavenumber k satisfies

k⁴ − τ₀k² + 1 = 0. (7)

Therefore, there are two possible wavenumbers, k_±, given byImportantly, the domain of k_±(τ₀) (such that k_±(τ₀) ∈ [Doublestruck R]) is bounded from below, where τ₀ = 2 gives the value where the two branches of the solutions meet, k₊(τ₀) = k₋(τ₀). We show in Appendix D that there are no solutions of the eigenvalue problem (eqn (5) and (6)) for complex k_±.

In general, w₀(x) will contain both of the wavenumbers given in eqn (8). The condition of zero vertical displacement at x = ±L/2 is satisfied by the symmetric and antisymmetric combinations viz.,

where A_s, A_a are, yet to be determined, constants. The remaining boundary condition that (w₀)_,x(±L/2) = 0 leads to the relations between k₊ and k₋;which are associated with the symmetric (9a) and antisymmetric (9b) modes respectively. Since k_± = k_±(τ₀), the solutions of eqn (10a) determine the compression τ^(s)₀ required to produce the corresponding symmetric profile w^(s)₀(x) from eqn (9a). Analogously, the roots τ^(a)₀ of eqn (10b) are associated with the odd functions w^(a)₀(x) from eqn (9b). For a given value of L, each relation (10a) and (10b) has an infinite number of solutions, the smallest of which is τ₀ = 2. However, when τ₀ = 2, k₊ = k₋ and eqn (9) gives the trivial solution w₀(x) ≡ 0. Therefore, each value of τ₀ > 2 that solves either of eqn (10) corresponds to a different, and non-trivial mode of buckling in the sheet. We are concerned with determining only the lowest mode of buckling, which corresponds to the smallest value of τ₀ > 2 that solves either eqn (10a) or eqn (10b) for a given sheet length L.

3.2 Asymptotic results for large sheets

Fig. 3 shows numerical results for the two smallest values of τ₀ at the onset of instability as a function of the sheet length L. The quantity of principal interest is the smallest buckling load, τ₀ > 2. These numerical results suggest that the critical compressive load at the onset of wrinkling corresponds to a symmetric or antisymmetric mode depending on the precise value of L. To understand this behaviour, we also note that the critical buckling load appears to approach 2 from above as the sheet length L → ∞. We use this as motivation to seek an asymptotic solution of eqn (10a) and (10b) for ε = τ₀ − 2 ≪ 1. This analysis (see Appendix E) reveals that:for symmetric modes, whilefor antisymmetric modes. These results go to higher order in L⁻¹ than the equivalent result by Rivetti and Neukirch.²⁶

Fig. 3 Minimum compressive force τ₀ as a function of the size of the confined sheet L. (a) τ₀(L) for medium length sheets on a log–log plot. (b) log–log plot of τ₀(L) − 2 for large sheets showing the good agreement between numerical solutions of (10a) (black solid) and (10b) (gray solid) and the asymptotic relationships (11) (green dashed) and (12) (red dashed), respectively.

Note that, in either case, the result that τ₀ → 2 suggests that k_± → 1 and hence that the natural dimensionless wavelength is λ = 2π. It is therefore often useful to measure the sheet length in terms of the number of half wavelengths and so we introduce

L̃ = L/π. (13)

We use L and L̃ interchangeably, for convenience.

The asymptotic predictions (11) and (12) show excellent agreement with our numerical solutions of (10a) and (10b) for L̃ ≳ 8 (see Fig. 3(b)). They also illustrate a L̃⁻² power-law (red triangles in both Fig. 3(a) and (b)) and how the mode with smallest critical compression oscillates between the symmetric and antisymmetric modes as sin L oscillates. Asymptotically we have that the mode at onset should be antisymmetric when sin L > 0 (i.e. when 2nπ < L < (2n + 1)π), while it should be symmetric when sin L < 0, (i.e. when (2n − 1)π < L < 2nπ), for integer n. A more detailed discussion of this feature, for general values of L, is given in Appendix F.

The asymptotic results for the critical loads also allow us to give simpler expressions for the symmetric and antisymmetric mode shapes in the limit L ≫ 1 (see Appendix E):

and

These results show how the mode shapes consist of a sinusoid of (short) wavelength equal to 2π, whose amplitude is modulated by another sinusoid with a large wavelength, equal to 2L. This result also makes clear that the amplitude of the wrinkles varies spatially as a result of a beating phenomenon—this spatial variation of wrinkle amplitude is similar to, but has a simpler underlying cause than, the spatial variation observed by Tovkach et al.²⁷

Fig. 4 shows a comparison between the asymptotic expressions of (14) with the exact mode shapes predicted by eqn (9), for a moderately large sheet (L̃ = 8). The agreement between the asymptotic predictions and the shapes given by eqn (9a) and (9b) is good, and small deviations only become visible near the ends, which is more pronounced for the symmetric mode, shown in the green dashed line of Fig. 4. Note that the boundary condition eqn (6b) is only satisfied at O(L⁻¹).

Fig. 4 Mode shapes corresponding to the two smallest compressions for a sheet of size L̃ = 8. Results are computed from the exact solutions, eqn (9) plotted with solid lines, and the asymptotic expressions, eqn (14) in dashed lines. The amplitudes are computed assuming the condition (a_c − a) = 10⁻³ ≪ 1, which ensures that the sheet bends close to the onset of buckling, imposed in eqn (18).

Our asymptotic results are related to those reported by Pocivavsek et al.⁹ Indeed, the linear combination of A_s = sin(L/2) in eqn (14a) and A_a = cos(L/2) in eqn (14b), yields the asymptotic solution of Pocivavsek et al.⁹

3.3 Critical confinement for buckling

Having determined the critical compressive force required to obtain buckling, τ₀(L), we now determine the properties of the compressed sheet around this buckling threshold. The non-zero lower bound for the compression reveals that the sheet is subject to a uniform compression prior to buckling. We consider a sheet that is confined to a domain of size 2a slightly shorter than its undeformed length (L) such that its associated compression is τ_a < τ₀(L). (Here, the subscript “a’’ is used to distinguish the applied compression needed to confine the sheet to a length 2a from the values required for buckling, i.e. the eigenvalues of the mechanical equilibrium from eqn (5).) Since τ_a < τ₀(L) the sheet is stable in its flat configuration, i.e. w(x) = 0, and is uniformly strained in the x-direction. The horizontal strain ε_xx may be linked to the horizontal displacement u(x) and also to the compressive stress τ_avia Hooke's law, givingwhich, upon integration from −L/2 to L/2 givesso that the critical value of a at which the sheet buckles is

The sheet can accommodate any imposed compressive stress τ_a < τ₀(L), or equivalently any confinement a > a_c, through in-plane compression. Below that length, the sheet buckles. As further confinement is imposed (a < a_c) the sheet's stress remains at the eigenvalue τ₀(L) and it tends to accommodate this further confinement by bending out of plane, rather than compressing in plane. We now consider how this amplitude grows just beyond the onset of buckling, i.e. for a < a_c.

3.4 Amplitudes

Having determined τ₀ by solving eqn (10), the solution for the shape of the buckled sheet is given by eqn (9) up to the multiplicative constants A_s and A_a. Once the sheet has buckled out of plane, the expression for the strain in eqn (15) is modified to (see Appendix A). Integrating this from x = −L/2 to x = L/2 and using symmetry of about x = 0, results in

Using w₀(x) from eqn (9a) and (9b) we find

orwhich can be substituted into eqn (18) to give A_s and A_a. Since k_± and τ₀ are determined by the numerical solution of eqn (10a) or (10b), the amplitudes depend only on the value of L. Therefore, for given values of a and L, the shape of the sheet can be completely determined numerically, including the amplitude, up to a sign.

We classify the sheets according to their size: Short length sheets are shorter than the intrinsic length scale of the wrinkles, [small script l]*, (L̃ < 1), medium length sheets are a few length scales in size [L̃ = O(1,10)], and large length sheets are many length scales in size (L̃ ≫ 1). This classification is relevant to the choice of the displacement parallel to the direction of confinement.

3.5 Switching between symmetric and antisymmetric modes

For large L, the asymptotic analysis of Section 3.2 shows how the sheet deformation switches between antisymmetric and symmetric modes whenever L increases by π. A similar oscillation is observed in the numerical results for the smallest roots of eqn (10), plotted in Fig. 3(a) for medium length sheets, and Fig. 3(b) for large length sheets. However, this period of π is only attained asymptotically as L → ∞: for finite L the length between mode switches must be calculated separately, as shown in Appendix F.

When the modes switch, or their compressions cross as L varies, the two relations (10) are satisfied by the same values of τ₀ and L̃. The crossing points that we call type I are located at

with l ∈ [Doublestruck N]*; index l labels the crossings – starting at l = 1 for the crossing at the smallest L̃, which is of type I, and increasing with the length of the sheet. Similarly, for what we call type II crossings we haveThe index l also explicitly determines the corresponding compressions ‡. The second index in eqn (20) and (21) evaluates to 1 because here we are interested in the first pair of compressions. The crossings for higher modes are included in Appendix F.

For short length sheets (L̃ < 1) the compressions are inversely proportional to L̃² (see Fig. 3a). This power-law behavior can be understood by dimensional analysis; as the size of the confined sheet decreases it becomes the smallest length scale in the problem and hence . However, a fundamental assumption for the floating beam eqn (4) is that the size of the sheet is larger than its thickness, and hence we will not explore the limit in which this condition is violated.

4 Buckling of inhomogeneous sheets

The results of the homogeneous case show that the amplitude of the buckling mode at the onset of buckling varies spatially, with a maximum close to the centre of the sheet. One natural question is: can we control this amplitude variation, for example by introducing a stiffness gradient? To make this idea more concrete, we imagine adding liquid inclusions to our soft host such that there is a stiffness gradient within the beam, and solve the resulting problem numerically.

4.1 Theoretical setting

Imposing a gradient of the liquid volume fraction parallel to the direction of confinement, ϕ(x), means that the elastic coefficients also vary spatially, i.e. E = E(x) and ν = ν(x). Therefore, we have a variable bending stiffness, viz.,Note that while B* may vary due to changes in the thickness h*, we do not consider this possibility here, which would modify the isostatic floating requirement as noted above.

We anticipate that a gradient in ϕ will break the symmetries about x = 0 exhibited by homogeneous sheets. Naively, one might expect the portion of the buckled sheet with the largest amplitude to reside where the modulus is smallest, since bending is easiest here, but this poses further questions, such as where should such a sheet will experience maximum bending stress, and hence be most likely to fracture?

We consider for the numerical analysis a sheet that is embedded with stiffening and incompressible liquid inclusions (γ′ = 10, with γ′ ≡ l*/R* as per Appendix B; R* denotes the size of the inclusions), which are linearly distributed parallel to the direction of confinement:

where ϕ_end:= ϕ(x = −L/2) is the concentration of inclusions at the left-hand end of the sheet. The scaled bending stiffness B(x) equals the dimensionless effective Young's modulus for stiffened soft composites given by Style et al.:^18,28where ϕ(x) is given by eqn (23). The control parameters here are ϕ_end, and the size of the sheet, L.

4.2 Spatial variation of amplitude

For medium length sheets [L̃ = O(1,10)], we numerically solve the equation of mechanical equilibrium (4) using the Chebfun package³⁶ to obtain the buckling profiles corresponding to the smallest compression. The contour plot of w(x) for different lengths L̃ in Fig. 5(a) shows that the maximal wrinkle amplitude (the trajectory in green) is shifted to the softer side of the sheet (x > 0). Indeed, the buckling profile at L̃ = 8, in Fig. 5(b), is no longer symmetric, unlike the modes of the homogeneous sheet (L̃ = 8) in Fig. 4. This intuitive feature is also observed in the buckling profiles of an inhomogeneous column.²⁵ We also see that the wrinkle wavelength is not noticeably altered by the change in elastic modulus. To understand this, we note that our asymptotic solution for large homogeneous sheets, leading to (14a) and (14b), shows that the (short) wrinkle wavelength depends on the sum k₊ + k₋, while the large wavelength of amplitude modulations is controlled by the difference k₊ − k₋; as such the fine scale wrinkle wavelength is affected only at higher order in spatial variation of properties than the amplitude modulation itself, as is observed. (To plot the function w(x), we use the integral constraint eqn (18), with (a_c − a) = 10⁻³ ≪ 1; this is close to the onset of buckling where our theory is valid.)

Fig. 5 (a) Contour plots for w(x) with varying L̃ and the largest vertical displacement (highlighted in green), for a medium length stiffened sheet. Here, ϕ_end = 0.1, γ′ = 10. The amplitudes are computed assuming the condition (a_c − a) = 10⁻³ ≪ 1, that ensures that the sheet bends close to the onset of buckling, imposed in eqn (18). (b) Buckling profile in a sheet with L̃ = 8.

Despite the lack of symmetry, the buckling profiles of a medium length inhomogeneous sheet are similar to the fundamental symmetric/antisymmetric modes characteristic of the homogeneous sheets. We project the numerical solutions for the inhomogeneous sheet, w(x), onto the modes of the homogeneous sheet (eqn (9)), i.e. the basis of eigenfunctions of the eigenvalue problem given by eqn (5).§ More concretely, we expand the solution of eqn (4) as the infinite linear combination

where we use the index j to denote each pair of solutions—symmetric (s) and antisymmetric (a). We redefine the eigenfunctions (9a) and (9b) asandand choose to normalize to unity the coefficients in the linear combination eqn (25).

The coefficients a^(s)₁ and a^(a)₁ are plotted in Fig. 6a. The first pair of modes—either the symmetric or the antisymmetric—are dominant in the linear combination that describes w(x). This dominant term changes in the vicinity of the crossing points (depicted by vertical dashed lines in Fig. 6). However, for larger and more inhomogeneous sheets the residue that is not captured by the lowest symmetric and antisymmetric modes, 1 − [(a^(s)₁)² + (a^(a)₁)²], plotted in Fig. 6b, grows. (Those cases with larger values of ϕ_end are plotted in lighter gray in Fig. 6.) This indicates a reduction in the projection of w(x) onto the first pair of modes.

Fig. 6 Smallest pair of squared coefficients in the expansion (25) for a stiffened inhomogeneous sheet (γ′ = 10). (a) Solid curves denote [a^(s)₁]², and dotted curves [a^(a)₁]². (b) Residual that is not captured by the lowest symmetric and antisymmetric modes. Different shades of gray for trajectories corresponding to each value of ϕ_end, from ϕ_end = 0.02 (black), to ϕ_end = 0.2 (lightest gray), and equally spaced jumps for the curves in between (increasing in the direction of the arrow in (b)). The dashed vertical lines denote the crossing points as in eqn (20) and (21) (j = 1, l = {1,2,3,4}).

Higher order modes are necessary to account for the buckling profile of longer and more inhomogeneous sheets. Thus the difference between the wrinkles in homogeneous and inhomogeneous sheets increases. This superposition of modes also explains the concentration of the bending energy towards the softer end of the sheet, as seen in Fig. 7 for the long inhomogeneous sheet.^9,29–31

Fig. 7 Displacement of the spatial variation in amplitude in the buckling profile w(x) for a large length stiffened inhomogeneous sheet (ϕ_end = 0.1, γ′ = 10) as L varies. The sheet is stiffer at x = −L/2 and softer at x = L/2. (a) Contour plots for w(x) for varying L̃ and the largest vertical displacement in green. (b) Buckling profile in a sheet with L̃ = 30. (a_c − a) = 10⁻³ to compute the amplitudes from eqn (18).

4.3 The failure of hard composites

Our study of wrinkles thus far, leads to further questions. For example, when we consider brittle materials, such as ice, in the homogeneous case we expect failure (e.g. fracture) to happen at the highest bending stress, i.e. at the centre. If we introduce a stiffness gradient, we have shown above how the position of the maximum amplitude is displaced from the centre. How does this position compare with that of the maximum bending stress? In other words, where is fracture more likely to occur in an inhomogeneous sheet?

Because ice is a hard material, we model the composite structure of water inclusions with the effective Young's modulus and the Poisson's ratio of Eshelby¹⁷ as

and(see Appendix B). A model for fracture in ice floes is given in Vella and Wettlaufer;² when a crack is initiated, the stresses exceed a critical value .¶

For elastic sheets, stresses vary linearly through the thickness of the sheet, and so the maximum stress is achieved at the surface of the sheet. This stress is related to the maximum bending moment per unit length acting on an element of the sheet (see p. 5 of ref. 4):

where the maximum bending moment per unit length iswhere is the curvature of the sheet in the small slope approximation (which is implicit in the derivation of the floating beam eqn (4)⁴). This maximum represents a trade-off between the bending stiffness, which peaks at the stiffer end of the sheet, and the curvature of the buckling profile, which is larger towards the softer side of the sheet.

We rescale the bending moment per unit length by , and the stress by . Now, considering a sheet with constant thickness h* and varying Young's modulus, the dimensionless failure criterion is

σ_m ≤ 6|M_max|. (32)

We plot in Fig. 8(a) the buckling profiles of a thin floating sheet with a decreasing volume fraction of softening inclusions (eqn (23)) and elastic moduli given by eqn (28) and (29) along the direction of confinement (L̃ = 19.6). The vertical axis in the contour plot denotes an increasing volume fraction at x = −L/2, ϕ_end. The dimensionless stretching-stiffness that models the behaviour of a fresh ice sheet, with , , h* = 1 mm and ν₀ = 0.33,^32,33 floating on water, is a large value of [scr S, script letter S] ≈ 10⁴. To compute the amplitudes of the buckling profiles using eqn (18), we assume (a_c − a) = 10⁻⁵.

Fig. 8 Bending profiles w(x) and maximum bending stress of a thin, wide softened sheet (L̃ = 19.6) close to the onset of buckling ((a_c − a) = 10⁻⁵). (a) Contour plots for w(x) with varying volume fraction ϕ_end are computed numerically.³⁶ The largest vertical displacement for every profile of a given ϕ_end is marked by green dots. Red dots denote the position of the maximum bending stress. (Greyscale is used to show the value of w(x), as indicated by the bar to the right.) (b) The maximum bending stress (dots) and the rescaled version of the yield strength ^32,33 (horizontal dashed line).

In addition to the spatial variation of the wrinkle amplitude towards the softer end of the sheet, which is more prominent at larger values of the volume fraction ϕ_end, there is a deviation between the position of the largest deflection (green dots) and the coordinate of |M_max| (red dots). To wit, the sheet's maximum bending stress is not at the wrinkle of largest amplitude, but at wrinkles on the stiffer side.

The maximum bending stress and the largest deflection in Fig. 8(a) move stepwise towards the soft end as ϕ_end increases, and the separation between their positions changes with ϕ_end. Now, because the largest curvature occurs at the largest deflection, and |M_max| is proportional to the curvature (see eqn (31)) then |M_max| is a discontinuous function of ϕ_end. In Fig. 8(b), the discontinuities correspond to when the maximum bending stress jumps to the left and hence approaches the location of the largest deflection.

Finally, we note that the maximum bending stress, plotted in Fig. 8(b), increases in inhomogeneous sheets.

4.4 Critical compression

The buckling profiles of an inhomogeneous sheet have no symmetry for any finite L. Thus, the switching between symmetric and antisymmetric modes observed in the homogeneous case as L varies has no analogue in the inhomogeneous case. We recall that those sheet lengths at which the homogeneous modes switch correspond to degenerate values of the associated compressions, and hence these degeneracies should disappear in the inhomogeneous case. This is a well established result in quantum theory, where the jargon for degenerate eigenvalues is ‘level crossing’, the simplest case of which was analyzed by von Neumann and Wigner.³⁷ From a mathematical perspective, the problem here is very similar: the linear differential operator defining the mechanical equilibrium, eqn (4), can be parametrized as [scr A, script letter A](ϕ_end) = [scr A, script letter A]_H + ϕ_end [scr A, script letter A]_I. When ϕ_end = 0, this parametrization yields the homogeneous case, eqn (5): [scr A, script letter A]_Hw₀ + τ₀(w₀)_,2x = 0. The problem of finding the eigenvalues of [scr A, script letter A](ϕ_end) by perturbation theory, which expresses each eigenvalue as a power series in ϕ_end, then shows the unification of the eigenvalues into a single multi-valued analytic function of ϕ_end, and the avoidance of crossing for ϕ_end ≠ 0.^38,39

We examine in Fig. 9 the two smallest compressions of inhomogeneous sheets as the sheet size varies. We use the rescaling (Ē denotes the spatial average of eqn (24)) in the inhomogeneous case and include for comparison the smallest pair of critical buckling compressions from the homogeneous sheet. We note that this scaling collapses the inhomogeneous case onto the homogeneous case. However, by blowing up the region around L̃ = L̃^(II)_1,1 (see the inset in Fig. 9, in which the inhomogeneous curves are computed every ΔL̃ = 0.032), we see that the compressions of the inhomogeneous sheet do not cross. In other words, they are ordered from smaller to larger for any finite L̃. We thus denote the two smallest compressions τ⁽¹⁾(L̃) and τ⁽²⁾(L̃), with τ⁽¹⁾(L̃) < τ⁽²⁾(L̃). Importantly, τ⁽¹⁾(L̃) is larger than the smallest compression in the first pair (τ^(s)₀,τ^(a)₀), and hence an inhomogeneous stiffness induces larger compressive loads; at most, compressions that match the homogeneous counterpart are observed for certain values of L̃.

Fig. 9 Critical compression of a stiffened inhomogeneous sheet (ϕ_end = 0.2, γ′ = 10) as a function of the sheet size L̃, with numerical results using Chebfun³⁶ (τ⁽¹⁾, τ⁽²⁾). The quantities τ^(s)₀, τ^(a)₀ denote the smallest pair of compressions of the corresponding homogeneous sheet (ϕ_end = 0). Here, the critical compression τ has been rescaled by the square root of the mean Young's modulus, , to facilitate comparison with the results from the homogenous case (solid curves).

In a thought experiment in which we infinitesimally increase ϕ_end, starting at ϕ_end = 0, i.e. a homogeneous sheet, the intersecting compressions at L̃ given by eqn (20) and (21) split when ϕ_end > 0. We measure this separation between consecutive compressions, and we examine its growth with ϕ_end close to the corresponding crossing points of the homogeneous sheet, in Fig. 10.

Fig. 10 Minimal separation between the two smallest compressions Δτ⁽²⁻¹⁾_min of a stiffened inhomogeneous sheet (γ′ = 10). Different symbols denote those minima corresponding to the crossings of type I (circles), and of type II (triangles) when ϕ_end/L = 0. Shades of gray correspond to the crossing point index l: lighter gray for l = 1, darkest gray for l = 4. The red right triangle is one decade in length on both legs, showing the linear relation Δτ⁽²⁻¹⁾_min ∝ ϕ_end/L.

We fix ϕ_end and then we find numerically the local minima of Δτ⁽²⁻¹⁾ with respect to L̃, which we denote Δτ⁽²⁻¹⁾_min. The values of L̃ that minimize Δτ⁽²⁻¹⁾ start at the crossing points when ϕ_end = 0, and increase with ϕ_end. The first eight results Δτ⁽²⁻¹⁾_min (corresponding to the first eight crossings in the homogeneous case) are plotted in Fig. 10 for different values of ϕ_end/L, where L is the corresponding sheet size that minimizes Δτ⁽²⁻¹⁾. All curves Δτ⁽²⁻¹⁾_min collapse into one. Therefore we conclude that Δτ⁽²⁻¹⁾_min grows linearly with ϕ_end/L, which is the gradient from eqn (23).

5 Conclusions

We have examined the two-dimensional buckling and wrinkle patterns in floating homogeneous and inhomogeneous thin elastic sheets. With two control parameters, the size of the confined sheet and the gradient of bending stiffness, we quantified the wrinkled states using the Föppl–von Kármán theory of thin sheets. A central test of the results is to vary the bending stiffness by varying the volume fraction of inclusions in the host solid.

In homogeneous sheets, the only control parameter is the sheet size L. The buckling profile is determined by the smallest compressive load, and so we expect that the mode that will be observed in a buckling experiment corresponds to the smallest compression. However, for some particular sizes of the confined sheet, the same compression is associated with two different wrinkling modes. We gave asymptotic results for the shape of the sheet at the onset of buckling, together with the critical loads and the critical sheet sizes for degeneracy in the limit of large sheets, L ≫ 1.

In contrast, this degeneracy is not observed in inhomogeneous sheets. Indeed, the otherwise crossing compressions of the homogeneous case split when a gradient of stiffness is applied parallel to the direction of confinement. The size of this splitting grows linearly with the magnitude of the gradient of the volume fraction at all the crossing points. Importantly, the wrinkled states of confined inhomogeneous sheets depend sensitively on their size. While medium length sheets buckle very much like their homogeneous counterparts, the wrinkled states in large length sheets are a superposition of many modes. This feature of large length sheets allows for the bending energy to be spatially concentrated, which is crucial in establishing a failure criterion, with particular relevance in glaciology.

Finally, the results presented here for floating sheets are also relevant for sheets on a linear soft elastic foundation, see e.g. ref. 40. A more complex behaviour is expected in the more general case of a linear elastic foundation, whose response is expected to be geometrically nonlinear, in which localization of buckling occurs, see e.g., ref. 41 and 42. In addition to the range of applications of interest, from soft composites of biological relevance⁴³ to hard composites of engineering or geophysical importance, a thorough mathematical analysis, rather than the numerical study given here, may provide additional insights.

Author contributions

MS, DV and JSW conceived and designed the research. MS, CA and AFB developed and implemented the theory. MS and CA conducted and analysed the simulations. MS, CA, AFB, DV and JSW wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.

Appendices

Appendix A Details of theoretical formulation

In the main text, we gave the equation governing the out-of-plane displacement of the beam without a formal derivation. Here, we expand upon the derivation of this. A key detail concerns the state of stress within the sheet—to ensure that this stress satisfies in-plane equilibrium, ∇·σ = 0, we introduce a force function χ*: the internal in-plane forces per unit length are obtained by double differentiation of χ* so that , and σ_xy = −∂²χ*/∂x*∂y*.⁴ Following the standard derivation of the plate equation, see ref. 4, but accounting for the possibility that ν = ν(x*), we find that normal displacements of the sheet satisfy:

∇²(B*(x*)∇²w*) − [B*(x*){1 − ν(x*)},w*] + ρ*g*w* = [w*,χ*], (33)

where B*(x*) is the bending stiffness (or flexural rigidity) of the sheet and the von Kármán operator is

For a midplane displacement with components (u*,v*,w*), the sheet's in-plane strains, ε_ij, are given by

where, for example, (w*)_,x* denotes the partial derivative of w* with respect to x*. Note that the displacements u* and v* may be eliminated from these relationships by cross-differentiation, see p. 13 of ref. 4. Relating these derivatives of strains to the in-plane forces, and hence to the derivatives of the force function χ*, one finds the compatibility equationwhich gives the stress in the plane of the sheet induced by the stretching of the sheet's mid-plane.

Appendix B Effective stiffness of composite materials

The foundational Eshelby theory of solid composites¹⁷ describes the elastic behaviour of rigid composites with a dilute dispersion of noninteracting incompressible inclusions. Stiff-matrix materials such as ice, glass, ceramics and steel have E* = O(GPa) and ν ∼ 0.3, and thus have subnanometric elastocapillary length. Therefore, for typical inclusion sizes the effect of surface tension is negligible and we can use Eshelby theory¹⁷ to compute the effective elastic moduli of compression, κ*, and rigidity, μ*, which arewhere α ≡ (1 + ν₀)/[3(1 − ν₀)] and β ≡ 2(4 − 5ν₀)/[15(1 − ν₀)]. We denote the host matrix's elastic constants with subscript “0”, those corresponding to the inclusions with subscript “1”, and symbols with no subscript denote the solid composite. The volume fraction of inclusions is ϕ.

The incompressible liquid inclusions have zero shear modulus, , and infinite bulk modulus, (due to incompressibility), so that the Young's modulus and Poisson's ratio of the composite are

andwhere we have used eqn (37) and expanded to first order in ϕ. For a stiff-matrix composite with ν₀ = 0.3 with liquid inclusions, the composite Young's modulus (Poisson's ratio) is less than (greater than) the host matrix; and ν > ν₀.

For soft composites, micron sized inclusions create non-negligible interfacial stresses with an effective Young's modulus given by

where γ′ ≡ l*/R* captures the size regime where surface tension operates.^18,28Eqn (24) assumes that the inclusion concentration is dilute and hence we refer to it as the dilute theory (DT). However, we note that this approach is quantitatively accurate up to ϕ ≈ 0.2 when γ′ > 2/3, which is in the stiffening regime where .⁴⁴ The constituents of the composite are incompressible and hence ν = 1/2 throughout.

In the softening regime, where γ′ < 2/3, we use the expression for the effective Young's modulus of Mancarella et al.¹⁹ (MSW):

Appendix C The eigenvalue problem for a homogeneous sheet

Here we solve the linear fourth-order differential equation for the vertical displacement w₀,where the eigenvalue τ₀ is determined by the requirement to have a non-trivial that satisfies the homogeneous boundary conditionsSince eqn (42) has constant coefficients, we seek solutions of the form w₀(x) ∝ e^ikx, yielding

k⁴ − τ₀k²+ 1 = 0, (44)

and the solutions for k² areor

Thus for the wavenumber k to be real, we need τ₀ ≥ 2. The eigenfunctions w₀(x) are linear combinations of e^±ik₊x and e^±ik₋x, but we are interested in real solutions. Thus, we write the general solution of eqn (42) as

w₀(x) = C₁ cos(k₊x) + C₂ cos(k₋x) + C₃ sin(k₊x) + C₄ sin(k₋x), (47)

where the k_± are the positive roots of eqn (45) and the C_i (i = 1,2,3,4) are real constants. The problem is linear and has homogeneous boundary conditions, so those constants can only be determined up to a multiplicative factor. Enforcing eqn (43), we obtainThe vanishing determinant of the system (48) can be factorized asHence, for a given L we have two possible relations involving k_± and L:

k₊tan(k₊L/2) = k₋tan(k₋L/2) (50)

or k₊cot(k₊L/2) = k₋cot(k₋L/2), (51)

where ‘or’ means that either relation (50) or relation (51) is fulfilled, or both. Note that when only eqn (50) is fulfilled, we must set C₃ = C₄ = 0 for the last two equations in (48) to be satisfied, and hence the resulting eigenfunction is even. Similarly, when only eqn (51) is fulfilled, we must set C₁ = C₂ = 0 and the eigenfunction is odd.

Clearly when τ₀ = 2, k_± = 1 and the relations (50) and (51) are both satisfied simultaneously. However, the resulting solution is trivial, w₀(x) = 0. In Appendix F we ask for which values of L both relations (50) and (51) are satisfied simultaneously with τ₀ > 2 and show that there are certain values of L, the ‘crossing points’, for which both odd and even solutions emerge with the same value of τ₀. More generally, however, one of the relations (50) and (51) has a smallest value of τ₀ > 2. We therefore expect that as the compressive stress τ₀ is increased from 0, the mode with the smallest value of τ₀ will be obtained; this emergent buckling mode will therefore be symmetric or antisymmetric depending on which of the relations (50) and (51) is solved by the smaller value of τ₀. In Appendix E we determine asymptotic expressions, valid for L ≫ 1, for the smallest τ₀ > 2 that satisfies each of the relations (50) and (51); this allows us also to determine which is the smaller compression and hence which mode, symmetric or antisymmetric, should be expected at the onset of wrinkling.

We are generally interested in the dependence of the eigenvalue τ₀ on the natural length of the sheet L. We note that we can rewrite the characteristic eqn (44) as

which holds for k being either k₊ or k₋. We also note, from the original quartic, that k₊²k₋² = 1 and hence, taking k_± to be positive, we have

k₊k₋ = 1. (53)

Appendix D The buckling wavenumber for homogeneous sheets is real

In Section 3 we assumed that the wavenumber k observed in buckling is purely real. Here, we demonstrate this is the case by supposing instead that eqn (44) has complex roots. Since the tension τ₀ is real and eqn (44) contains only even powers, there must be two complex conjugate pairs of solutions. We may write one pair as k_± = k_r ± ik_i, with k_r ≥ 0, and hence the other pair will be −k_±. We extend sine and cosine to the complex plane using analytic continuation, and thereby extend the boundary conditions detailed in Appendix C to obtain the counterparts of eqn (50) and (51) as

(k_r + ik_i)tan(k_r + ik_i)L/2 = (k_r − ik_i)tan(k_r − ik_i)L/2 (54)

and

(k_r + ik_i)cot(k_r + ik_i)L/2 = (k_r − ik_i)cot(k_r − ik_i)L/2, (55)

where for both relations (54) and (55) the right-hand side is the complex conjugate of the left-hand side, so the imaginary part of either side must be zero. This condition takes the form

f(k_rL/2) = −g(k_iL/2) and f(k_rL/2) = g(k_iL/2), (56)

for eqn (54) and (55), respectively, where

By plotting the functions f, g and −g one finds that their ranges do not overlap (though f and g have the same limit as x tends to zero) and hence there is no solution of eqn (56). Therefore, k is real.

Appendix E Asymptotic solution for large sheet sizes, L ≫ 1

Having shown that the roots of eqn (50) and (51) are necessarily real, we consider in this appendix the behaviour of these roots for large sheets, L ≫ 1. Our starting point is the observation, from numerical simulations, that as L → ∞ it appears that τ₀ → 2 for both symmetric and antisymmetric modes. We therefore let

τ₀ = 2 + ε, (58)

with ε ≪ 1. From (8) we then immediately have that

We consider first the case of symmetric modes, rewriting eqn (50) as

which can then be written in terms of ε asThe quantity ε^1/2L appears frequently, and so we letfor some α which is an O(1) quantity to be determined. We find that eqn (61) becomesA non-trivial solution requires sinα ≪ 1, and hence that α ≈ nπ for some integer n. The smallest τ₀ > 2 corresponds to the smallest positive α so that the relevant root is α ≈ π. A simple calculation of the correction α–π from eqn (63) then yields the asymptotic expression for the value of τ₀ for the even mode, τ^(s)₀, that is given in eqn (11) of the main text. Precisely the same calculation, with minor modifications of signs on the right hand side of eqn (60)–(63), follows through for the asymmetric mode and yields eqn (12) for τ^(a)₀.

We note further that the asymptotic expressions for τ^(a)₀ and τ^(s)₀ agree to leading order in L⁻¹ when sin L = 0; therefore we expect the crossing points to occur at L = nπ with n ≫ 1 integer.

The asymptotic expressions for symmetric and antisymmetric mode shapes, given in eqn (14) of the main text, follow from expanding eqn (9) to leading order in L⁻¹.

Appendix F Crossings of symmetric and antisymmetric modes

Here, we determine expressions for the values of L for which both eqn (50) and (51) are satisfied simultaneously; this corresponds to the sheet lengths for which a symmetric mode and an antisymmetric mode have the same eigenvalue. We refer to such a point as a ‘crossing’.

We multiply eqn (50) and (51) to obtain k₊² = k₋² and hence k₊ = k₋ = 1. However, this implies τ₀ = 2, which, as we have already seen, corresponds to the trivial solution. Instead, we must have that either tan k₊L/2 = tan k₋L/2 = 0 or cot k₊L/2 = cot k₋L/2 = 0, thereby allowing simultaneous solutions of eqn (50) and (51) with k₊ ≠ k₋. Thus, there are two families of non-trivial common solutions to eqn (50) and (51):

(1) k₊L = π + 2mπ and k₋L = π + 2nπ, with m,n ∈ [Doublestruck N]. Eqn (53) then implies that

We refer to the eigenvalues τ₀ at these crossing points as type I; they are given by (see eqn (52))

(2) k₊L = 2m̃π and k₋L = 2ñπ, with m̃,ñ ∈ [Doublestruck N]* (the set of non-zero natural numbers). Eqn (53) implies that

We refer to the eigenvalues τ₀ at these crossing points as type II; they are given by

Each symmetric and antisymmetric mode comes in a pair, which we label with the index j = m − n ∈ [Doublestruck N]* such that j = 1 corresponds to the pair with the smallest eigenvalues. (Note that since k₊ > k₋, m > n.) We find that the crossings within the pair j occur for

(m_j,n_j) = (l − 1,l − 1 + j) and (m̃_j,ñ_j) = (l,l + j), l ∈ [Doublestruck N]*. (68)

The index l labels the crossings within a given pair, such that l = 1 corresponds to the smallest size for which a crossing of either type occurs. Hence, we obtain two sets of crossing points given byThe corresponding eigenvalues areFor very small values of L, within a pair the symmetric mode always has the smallest eigenvalue. As L increases, the first pair crosses at (type I), corresponding to τ₀ = 10/3. The next crossing (type II) occurs at , corresponding to τ₀ = 5/2. Between those crossing points, the antisymmetric mode has the smallest eigenvalue. As L increases this pattern repeats infinitely many times (see Fig. 3). Indeed, as l → ∞, we note that:

L^(I)_l,1 ∼ 2πl, L^(II)_l,1 ∼ π(2l + 1), (71)

reproducing the result of the asymptotic analysis for L ≫ 1 that followed on from eqn (11) and (12), namely that the system should switch between symmetric and asymmetric modes (and vice versa) each time L increases by a multiple of π.

Acknowledgements

We thank A. Souslov and D. Mitra for helpful comments. MS, AFB and JSW acknowledge the support of Swedish Research Council Grant No. 638-2013-9243. CA acknowledges the support of the Swedish Research Council Grant No. 2018-04290. DV acknowledges the Leverhulme Trust for support through a Philip Leverhulme Prize. Nordita is supported in part by NordForsk.

Notes and references

1. D. Vella and J. S. Wettlaufer, Phys. Rev. Lett., 2007, 98, 088303.
2. D. Vella and J. S. Wettlaufer, J. Geophys. Res., 2008, 113(C11), C11011.
3. P. J. Hudleston and S. H. Treagus, J. Struct. Geol., 2010, 32, 2042–2071.
4. E. H. Mansfield, The Bending and Stretching of Plates, Cambridge University Press, 2nd edn, 1989.
5. S. T. Milner, J. F. Joanny and P. Pincus, Europhys. Lett., 1989, 9, 495–500.
6. E. Cerda and L. Mahadevan, Phys. Rev. Lett., 2003, 90, 074302.
7. D. Vella, P. Aussillous and L. Mahadevan, Europhys. Lett., 2004, 68, 212–218.
8. L. Bourdieu, J. Daillant, D. Chatenay, A. Braslau and D. Colson, Phys. Rev. Lett., 1994, 72, 1502–1505.
9. L. Pocivavsek, R. Dellsy, A. Kern, S. Johnson, B. Lin, K. Y. C. Lee and E. Cerda, Science, 2008, 320, 912–916.
10. M. Trejo, C. Douarche, V. Bailleux, C. Poulard, S. Mariot, C. Regeard and E. Raspaud, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 2011–2016.
11. E. Jambon-Puillet, C. Josserand and S. Protière, Phys. Rev. Mater., 2017, 1, 042601.
12. W. T. Tekalign and B. J. Spencer, J. Appl. Phys., 2007, 102, 073503.
13. A. D. Kerr, J. Glaciol., 1976, 17, 229–268.
14. N. B. Coffey, D. R. MacAyeal, L. Copland, D. R. Mueller, O. V. Sergienko, A. F. Banwell and C.-Y. Lai, J. Glaciol., 2022, 68, 867–878.
15. I. Tobasco, Y. Timounay, D. Todorova, G. C. Leggat, J. D. Paulsen and E. Katifori, Nat. Phys., 2022, 18, 1099–1104.
16. U. C. V. Andrenšek, P. C. V. Ziherl and M. Krajnc, Phys. Rev. Lett., 2023, 130, 198401.
17. J. D. Eshelby, Proc. R. Soc. London, Ser. A, 1957, 241, 376–396.
18. R. W. Style, R. Boltyanskiy, B. Allen, K. E. Jensen, H. P. Foote, J. S. Wettlaufer and E. R. Dufresne, Nat. Phys., 2015, 11, 82–87.
19. F. Mancarella, R. W. Style and J. S. Wettlaufer, Proc. R. Soc. London, Ser. A, 2016, 472, 20150853.
20. S. Moser, Y. Feng, O. Yasa, S. Heyden, M. Kessler, E. Amstad, E. R. Dufresne, R. K. Katzschmann and R. W. Style, Soft Matter, 2022, 18, 7229–7235.
21. T. J. Wallin, L.-E. Simonsen, W. Pan, K. Wang, E. Giannelis, R. F. Shepherd and Y. Mengüç, Nat. Commun., 2020, 11, 4000.
22. C. Mo, H. Long and J. R. Raney, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2123497119.
23. J. R. Tse and A. J. Engler, Curr. Protoc. Cell Biol., 2010, 47, 10.
24. E. Lifshitz, A. Kosevich and L. Pitaevskii, Chapter II - The Equilibrium Of Rods And Plates, Butterworth-Heinemann, Oxford, 3rd edn, 1986, pp. 38–86.
25. M. Suñé and J. S. Wettlaufer, Phys. Rev. Mater., 2021, 5, 055603.
26. M. Rivetti and S. Neukirch, J. Mech. Phys. Solids, 2014, 69, 143–155.
27. O. Tovkach, J. Chen, M. M. Ripp, T. Zhang, J. D. Paulsen and B. Davidovitch, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 3938–3943.
28. R. W. Style, J. S. Wettlaufer and E. R. Dufresne, Soft Matter, 2015, 11, 672–679.
29. H. Diamant and T. A. Witten, Phys. Rev. Lett., 2011, 107, 164302.
30. B. Audoly, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 011605.
31. O. Oshri, F. Brau and H. Diamant, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 052408.
32. P. Hobbs, Ice Physics, OUP, Oxford, 2010.
33. E. M. Schulson, JOM, 1999, 51, 21–27.
34. W. F. Weeks and D. L. Anderson, Trans., Am. Geophys. Union, 1958, 39, 641–647.
35. R. J. Evans and N. Untersteiner, J. Geophys. Res., 1971, 76, 694–703.
36. Chebfun Guide, ed. T. A. Driscoll, N. Hale and L. N. Trefethen, Pafnuty Publications, 2014.
37. J. von Neumann and E. Wigner, Phys. Z., 1929, 30, 467–470.
38. C. Bender and S. Orszag, Advanced Mathematical Methods for Scientists and Engineers I, Springer-Verlag, New York, 1999.
39. P. Lax, Linear Algebra and its Applications, John Wiley and Sons Inc, 1996.
40. G. W. Hunt, M. K. Wadee and N. Shiacolas, J. Appl. Mech., 1993, 60, 1033–1038.
41. M. Potier-Ferry, in Amplitude Modulation and Localization of Buckling Patterns, ed. J. Thompson and G. Hunt, Cambridge University Press, 1983.
42. F. Brau, H. Vandeparre, A. Sabbah, C. Poulard, A. Boudaoud and P. Damman, Nat. Phys., 2011, 7, 56–60.
43. A. Goriely, The Mathematics and Mechanics of Biological Growth, Springer, New York, NY, 2017, vol. 45.
44. F. Mancarella, R. W. Style and J. S. Wettlaufer, Soft Matter, 2016, 12, 2744–2750.

---

Footnotes

† We treat sheets floating on water with density ρ* ∼ 10³ kg m⁻³, and the standard value of the gravitational acceleration g* ∼ 9.8 m s⁻².

‡ At the crossing points, the compressions are and for type I and type II crossings respectively. See Appendix F.

§ Integration by parts shows that the operators and are Hermitian, the latter is also positive definite, and hence is a Hermitian definite pencil. Therefore, the problem is one of standard Hermitian eigen-theory and we can take advantage of the completeness of the known eigenfunctions, eqn (9a) and (9b), of the homogeneous sheet problem, eqn (5).

¶ For ice floes the yield strength, , is 1–3 MPa for fresh ice,^32,33 and 0.1–0.4 MPa,³⁴ or 0.4 MPa,³⁵ for sea ice.

---