## Abstract

We consider the point indentation of a pressurized elastic shell. It has previously been shown that such a shell is subject to a wrinkling instability as the indentation depth is quasi-statically increased. Here we present detailed analysis of this wrinkling instability using a combination of analytical techniques and finite-element simulations. In particular, we study how the number of wrinkles observed at the onset of instability grows with increasing pressurization. We also study how, for fixed pressurization, the number of wrinkles changes both spatially and with increasing indentation depth beyond onset. This ‘Far from threshold’ analysis exploits the largeness of the wrinkle wavenumber that is observed at high pressurization and leads to quantitative differences with the standard ‘Near threshold’ stability analysis.

This article is part of the themed issue ‘Patterning through instabilities in complex media: theory and applications.’

## 1. Introduction

The buckling instability of elastic objects is a staple of much research in engineering, applied mathematics and physics. Traditionally, the focus has been on determining the properties of these instabilities close to the onset of instability with a view to avoiding them: where does instability occur and how sensitive to imperfections is this threshold [1,2]? However, more recently this ‘buckliphobia’ has, to a certain extent, been replaced by ‘buckliphilia’ [3]: the patterns that are commonly observed in instability (notably wrinkling, folding and creasing) are not only common in the natural world (for example, in skin [4], the cortical convolutions in the human brain [5] and in growing plants [6]) but are also useful in generating regular patterns in technological applications including photonics [7] and stretchable electronics [8]. In these applications, interest is focused not so much on the onset of instability, but rather the behaviour far from onset, where novel behaviour, such as period doubling [9,10], may emerge.

A classic example of elastic instability is the wrinkling of a thin elastic object (a beam) under compression on a liquid bath or elastic foundation [9,11,12,13]. As well as being visually striking, these wrinkle patterns allow, in principle, for a simple assay of a sheet’s mechanical properties [14]. The essential mechanics of the pattern selection are clear: there is a competition between destabilizing effects (here a compression) and stabilizing effects, such as a bending stiffness (as in Euler buckling). However, the selection of an intermediate length scale, different from the system size, requires a competition between restoring forces. For example, in the compression of an elastic beam (of bending stiffness *B*) on a foundation (of stiffness *K*), the resistance to bending is minimized through the selection of the largest wavelength possible while the foundation prefers the smallest wavelength possible. The wavelength that is observed emerges from a compromise between these two [15]; in particular,
1.1

In uniaxial compression, the result (1.1) applies both at the onset of instability and far beyond it. However, in other scenarios, the behaviour of the system changes qualitatively beyond onset. For example, attempts to determine the location of wrinkles based on determining where compression would occur ordinarily under-predict the wrinkle extent dramatically [14,16,17]. It is now well known that to correctly predict the location of wrinkles it is necessary to account for the fact that wrinkling relaxes the compressive stress that caused instability in the first place [18,19]—what is usually referred to as tension-field theory [20] or the relaxed energy approach [21]. However, an understanding of how this stress relaxation affects the wrinkle wavelength has only come more recently [15,22,23]. In particular, Cerda & Mahadevan [15] showed that tension along the direction of the wrinkles gives rise to an effective ‘tensional’ stiffness, *K*_{tens}, that may be substituted into (1.1) in place of *K*. Generalizing this work, Paulsen *et al.* [23] recently showed that curvature along the wrinkle direction can also give rise to a ‘curvature-induced’ stiffness, *K*_{curv}.

In this paper, we give a new perspective on wrinkling ‘Far from threshold’ (FT) for the indentation-induced wrinkling of a pressurized elastic shell [17,18]. We use both asymptotic analysis of the shallow shell equations, together with ABAQUS simulations. In particular, we perform the FT expansion [22] by exploiting the large number of wrinkles, *m*≫1, that are frequently observed in such scenarios, and determine results for the wrinkle number *m* as the indentation depth, pressurization and also radial position within the shell vary. This also allows us to generalize the curvature-induced stiffness [23] to the case in which the object has a natural curvature.

A key theme throughout this work will be comparing and contrasting the results at the onset of instability, or ‘Near threshold’ (NT), with those FT. Indeed, our study is the first that observes the transition from NT to FT in direct numerical simulations. We, therefore, begin with a discussion of the governing shallow shell equations in §2 before performing the standard NT analysis in §3. In §4, we present our ABAQUS simulations, which motivate the FT analysis presented in §5. We compare these results in §6 before summarizing our results in §7.

## 2. Shallow shell formulation

### (a) Governing equations

Figure 1 shows a schematic of a pressurized, spherical elastic shell, subject to an indentation at its pole. Provided that the horizontal length scale over which the shell is deformed by indentation is much smaller than the radius of curvature of the shell, *R*, a convenient framework within which to analyse the deformations of such thin elastic shells is ‘shallow-shell theory’ [24,17]: the intrinsic geometry of the shell is assumed to be identical to the geometry of the plane of its projection. We shall discuss the circumstances under which this assumption is valid later, but for now we introduce a cylindrical polar coordinate system (*r*,*θ*,*z*) and orient the shell with its planform lying in the (*r*,*θ*) plane. The key quantities of interest are the deflection of the shell in the normal direction, *w*(*r*,*θ*), and the stress components within the shell, *σ*_{rr},*σ*_{θθ} and *σ*_{rθ}. For the shell to be in equilibrium, the in-plane stress balances, i.e.
2.1must be satisfied, together with the normal force balance. For a spherical elastic shell of radius *R*, subject to a constant internal pressure *p* and a localized indentation force *F*, this normal force balance reads [24,17,25]
2.2Here *B*=*Et*^{3}/[12(1−*ν*^{2})] is the bending stiffness, with *E* the Young modulus, *t* the thickness and *ν* the Poisson ratio of the shell. The localized force *F* is implemented as a Dirac *δ*-function but may also be understood in the first integral of (2.2) as a shear force *Q*_{r}=−*F*/(2*πr*) present only outside a cylindrical region of radius *r*=*ϵ*→0 [24]. Note that in (2.2) the Laplacian operator is that used in plane polar coordinates, i.e.

The problem is closed by combining the linear stress–strain relationships (Hooke’s Law) with the geometrically nonlinear relationships between strain and the in-plane and normal displacements, and *w*(*r*,*θ*). For a spherical shell, these read
2.3
2.4
2.5(Here *ϵ*_{ij}, with *i*,*j*=*r*,*θ* denotes the entries of the symmetric strain tensor.) We note that in many problems, it is often convenient to introduce an Airy stress function *Φ* to ensure the in-plane stresses automatically balance; this then requires the solution of an equation for *Φ* that ensures the compatibility of strains. This approach is less convenient for the perturbation analysis that we shall undertake here.

The problem is to be solved subject to appropriate boundary conditions. At the point of indentation, we require that the displacement takes an imposed value, *δ*, that the shell does not have a cusp, and that the radial displacement *u*_{r} vanishes, i.e.
2.6Note also that, for convenience, we imagine imposing a given indentation depth, *δ*; the force *F* required to impose that displacement can, in principle, be determined as part of the solution, though we do not discuss the behaviour of this force here.

Far from the indenter, we require the shell to return to its undeformed state. We measure normal displacements from this state so that
2.7We also note that the pressure jump across the shell requires a uniform, isotropic tension within the shell *σ*_{rr}=*σ*_{θθ}=*pR*/2 prior to indentation (by analogy with Laplace’s Law). We, therefore, also require that
2.8

### (b) Non-dimensionalization

In previous studies [17,18], it has proved useful to introduce the shell ‘capillary length’
2.9as a natural horizontal length scale. (The shell ‘capillary length’ is the horizontal distance over which the pressure-induced tension, ∼*pR*, together with a restoring force from the curvature of the shell, *Et*/*R*^{2}, return the shell to flat following a small vertical deformation; it is somewhat analogous to the capillary length of a liquid–vapour interface.) By then introducing dimensionless variables
2.10With this non-dimensionalization, the normal force balance equation (2.2) becomes
2.11where is the dimensionless indentation force and
2.12is the dimensionless pressure. Alternatively, *τ*^{−2} is the dimensionless bending stiffness, so that *τ*^{2} may be thought of as the ‘bendability’ of the shell [22,26]. (Note now that the operator ∇ is written in terms of the dimensionless radius *ρ* in the obvious way.) The dimensionless in-plane stress balances become
2.13and
2.14while the stress–strain–displacement relations become
2.15
2.16
2.17

The non-dimensionalization also introduces a second dimensionless parameter,
2.18which is the dimensionless indentation depth. Generally, we shall be interested in the limit of highly pressurized (and hence highly bendable), and highly indented shells, i.e. . These limits may seem to be in contradiction with the additional requirements that the typical strains be small, *pR*/*Et*≪1, (for Hooke’s Law to remain valid) and that the lateral scale of deformation ℓ_{p}≪*R* (for the shallow shell theory to be valid). For a given indentation depth, shell radius and elastic modulus, these restrictions can in fact be satisfied in the limit of vanishing pressure and thickness, *p*,*t*→0, provided that *p*/*t*→0 but , as discussed elsewhere [18].

### (c) Axisymmetric behaviour

Axisymmetric solutions of (2.11)–(2.17) have been studied in some detail [17,25]. The key observation is that as increases, the hoop stress is non-monotonic: close to the indenter the stress becomes large, while far away it takes its undisturbed far-field value, . Between these two tensile regions, however, there is an annulus in which the hoop stress becomes relatively compressive; this behaviour can be understood as follows: indentation acts to pull material circles in to a position where their circumference is too long for their final position and so they become relatively compressed. As grows, this hoop stress becomes more and more compressive (relative to the far-field value, *σ*_{θθ}→1/2) until eventually it becomes absolutely compressive, . In figure 2, we see that this transition has already occurred by : as *τ* increases, the hoop stress in some region, even as . A more detailed calculation shows that with this transition occurs at [17].

In the limit of large *τ* the shell becomes a membrane-shell, i.e. it has zero resistance to compression and we expect that it should buckle for any . At finite values of *τ*, the classic problem is to try and understand *when* this buckling first occurs: what is the smallest indentation depth for which a buckled solution exists? What buckling mode is adopted at this buckling point? These questions may be addressed by a standard linear stability analysis, which we refer to as the NT analysis (as it is only valid very close to the threshold of instability).

## 3. ‘Near threshold’ analysis

To understand the onset of buckling in this problem, we begin by making a small perturbation, with wavenumber *m*, of the axisymmetric state. We, therefore, take the ansatz
3.1where superscript (0) indicates an axisymmetric quantity and a (1) indicates a small perturbation.

We saw above that in the limit there is a critical indentation depth, , at which first becomes compressive. We expect that for this will be a lower bound on the indentation required to cause buckling. We determine by linearizing the governing equations for the perturbation quantities, i.e. those with a (1) superscript. The result is a quadratic eigenvalue problem [27] for *m*^{2} with a given indentation depth ; this eigenvalue problem does not have a solution for , allowing the critical indentation depth, together with the associated mode number, to be determined numerically.

Our numerical results are shown in figure 3 and suggest that in the highly pressurized limit, *τ*≫1,
3.2These results can be rationalized by a scaling analysis of the governing equations, along the lines of the analysis presented in [22]; the details of this scaling calculation are to be presented elsewhere for a thin sheet floating at a liquid interface [28]. Here, however, we focus instead on first comparing these results with the results of more detailed numerical simulations, before probing what happens beyond the onset of instability.

## 4. Finite-element analysis

To provide an alternative perspective on the wrinkling of indented, pressurized elastic shells, we used finite-element simulations, following previous work [29,30,18]. These simulations were performed using the commercial finite-element code ABAQUS 6.14 (Dassault Systémes Simulia Corp., Providence, RI, USA). A linearly elastic (Young’s modulus *E*=7 GPa and Poisson’s ratio *ν*=0.3) spherical cap of radius *R*=1 m, base chord diameter 2*L*=1.6 m and thickness *t*=0.05 mm is discretized using quadrilateral linear shell elements with reduced integration and finite membrane strain (S4R). The number of elements is chosen to ensure that the resulting mesh is able to properly resolve azimuthal variations with a large wavenumber. Here, we investigate pressurizations up to *τ*=1888 (for which the number of wrinkles at onset is *m*_{NT}(1888)=70) in ABAQUS. We, therefore, ensure that there are enough elements to resolve azimuthal variations with a substantially larger wavenumber; specifically we ensure that our simulations are able to resolve wrinkles with wavenumber 3*m*_{NT}(1888)=210, and find that 94 891 elements is sufficient for this purpose. The same number of mesh elements is used in all simulations reported here, and no initial imperfection is imposed.

Our simulations are conducted in two steps: first, the shell is inflated to the desired pressure (*p*=1 kPa, *p*=5 kPa and *p*=10 kPa are investigated here, corresponding to *τ*=188, *τ*=944 and *τ*=1888, respectively). This inflation is performed using the built-in procedure *Static, General. The shell is then indented by imposing the desired displacement within a circular area of radius *L*_{A}=5 mm at the apex of the shell; this indentation mimics that caused by a flat punch whose size still remains much smaller than the shell width *L* and the typical horizontal length scale though we note that the limit of a point indenter is well defined and regular [25]. The indentation step is of the type *Dynamic, Implicit and stabilizes the solution by introducing a Rayleigh damping *α* proportional to the mass. Here the damping coefficient *α* and the density of the material of the shell are chosen arbitrarily; however, it is verified *a posteriori* that the amount of kinetic energy in the simulation remains less than 10% of the strain energy. The simulations may, therefore, still be considered to be in the *quasi-static* regime. Our numerical approach falls within the class of *dynamic relaxation techniques* in the classification proposed by Taylor *et al.* [31]: both a mass and a damping (proportional to the mass) are introduced and a dynamic problem is solved. In contrast with the *kinetic damping approach* used by Taylor *et al.* [32], however, we did not optimize the problem (reducing the dynamics to be dependent only on the density parameter).

The shapes of a deformed shell with *τ*=188 held fixed and increased, as found in our ABAQUS simulations, are shown in figure 4*a*–*d*. These images show clearly that the simulations predict a wrinkling instability, as expected. However, two features of these images are immediately inconsistent with the NT account of wrinkling. Firstly, the pattern evolves with the indentation depth, . Secondly, the wrinkle number is not spatially uniform: wrinkles appear at different angular positions, depending on the indentation depth (see particularly figure 4*c*,*d*).

Understanding the two features of the wrinkling pattern shown in figure 4, namely the variation of the wrinkle number in space and with indentation depth, will be our focus for much of the remainder of this paper. However, such a study requires a reliable way of quantifying the number of wrinkles at a given radial position, *ρ*. While the presence of wrinkles in a shell is visually evident in simulations, manually counting wrinkles is neither practical nor reliable. Other methods, including the use of Fourier transforms give broad ranges of the wrinkle number obscuring the trends evident in images of the shells, see e.g. figure 4*a*–*d*. Here, we examine the behaviour of the minimum principal stress (MPS): since wrinkling occurs when the hoop stress becomes sufficiently compressive (*σ*_{θθ}<0) and the radial stress remains tensile (*σ*_{rr}>0), the MPS should highlight both the direction of the compression and its frequency. Although the normal displacement might be a more natural candidate through which to detect wrinkles, the MPS is expected to vary sinusoidally in the same way. Empirically, we find that the MPS calculated by ABAQUS gives a less noisy signal than the normal displacement, making the identification of the wrinkle wavenumber this way more robust.

Colour maps of the MPS, with increasing and *τ*=188 fixed are shown in figure 4*e*. These plots show that the MPS highlights the wrinkled region: outside the wrinkled region *σ*^{(p)}_{min}=0, which is shown in grey. For a given dimensionless indentation depth , the MPS field at each node of the original mesh is mapped onto a regular polar grid so that the shell may then be cut into slices at increasing radial distance from the indentation point. (This interpolation is performed using the Matlab function *griddata*, equipped with the *nearest* method.) Once this has been done, the wavenumber at a given radial distance from the indentation point may be extracted by counting the troughs in the azimuthal variation of the MPS (figure 4*f*).

The raw results of this analysis for the simulations reported in figure 4 are shown in figure 5. As expected based on the images from the simulation, we see that the wrinkle number varies both with radial position and with indentation depth, . This is different from the result of the NT (or linear stability) analysis, which instead predicts a particular wavenumber at the onset of instability. We would like to be able to explain these two features of the patterns quantitatively and note that the wrinkle number *m*, although varying, is typically large, *m*≫1. We, therefore, seek a theoretical approach that exploits this largeness to describe the behaviour of the system FT.

## 5. ‘Far from threshold’ analysis

To understand the evolution of the wrinkle pattern as the indentation depth increases beyond the threshold value , we must go beyond the linear stability analysis, and account for the relevant nonlinearities. However, to make analytical progress we adapt the approach used by Davidovitch *et al.* [33] and Hohlfeld & Davidovitch [26], to expand the governing equations in terms of a suitable small parameter. In the case of the NT analysis, this small parameter is the wrinkle amplitude and hence the distance beyond the onset of the wrinkling threshold. To understand the behaviour of FT analysis, we exploit the observation from the NT analysis that the wrinkling wavenumber *m* diverges as : the parameter 1/*m* seems to be a candidate small parameter.

To make progress, we neglect, for the moment, the observation that the wrinkle number *m* varies in space and consider a single-mode ansatz for the normal displacement *W*(*ρ*,*θ*) within the wrinkled zone
5.1It is then entirely natural to seek expansions for the displacements and components of the stress tensor as
5.2where superscript (0) indicates the base solution and the superscripts (1) and (2) the higher order corrections. However, we emphasize that the base solutions will not, in general, be solutions of the fully axisymmetric problem [18].

The general form of the energy functional in the problem under analysis and in dimensionless variables reads
5.3where the terms
5.4are the bending energy density, the stretching energy density due to the radial and shear stresses, the stretching energy due to the hoop stress contribution and the work per unit area in compressing the gas, respectively. The FT approach divides this energy additively as , where is the portion of the system’s energy that does not depend on the (small) wrinkle amplitude, and hence dominates the portion that does depend on the wrinkle amplitude, , or the sub-dominant energy [22]. In this approach, the macroscopic properties of the wrinkle pattern (such as *where* the shell is wrinkled) is determined by minimizing , while the wrinkle number is determined by minimizing . The position of the wrinkles has been considered in detail before [17,18]. Here, we shall instead focus on the number of wrinkles *m* and so it is interesting to understand the general structure of , before formally deriving the order by order balances. Let us insert the expansion given in (5.1) and (5.2) and the relations in (2.15), (2.16) and (2.17) into the various energy terms presented in equation (5.4). Moreover, let us also anticipate that radial derivatives are small in comparison to azimuthal derivatives, and that, in particular, ∂/∂*ρ*∼1 while ∂/∂*θ*∼*m*. Since is, by definition, the contribution depending on the wrinkle amplitude and considering that the perturbation is sinusoidal, i.e. only the even powers remain after integration in the azimuthal direction, we find that the leading order contributions are
5.5The most general form of the energy depending on the wrinkle amplitude is, therefore,
5.6Although how *m* scales with *τ* is not known in advance, we see from (5.6) that the quantity in the first set of square brackets in the integral is minimized when the terms *τ*^{−2}*m*^{2} and *m*^{−2} are of the same order, so that *τ*=*O*(*m*^{2}). This is consistent with other FT problems in which *m*∼*ϵ*^{−1/4} with *ϵ*^{−1} the bendability, since here *ϵ*^{−1}=*τ*^{2} [15,22]. We also note that we expect (since compression is required for wrinkling); the third term in the integrand of (5.6), therefore, decreases the energy of the system and hence is the destabilizing term.

### (a) Expansion of the governing equations

The expression for the sub-dominant energy leads us to expect that *τ*=*O*(*m*^{2}). We, therefore, now proceed to expand the stress equilibrium equations in (2.11), (2.13) and (2.14) together with the stress–strain relations in (2.15), (2.16) and (2.17) order by order in *m*^{−1} using the expansions given in equations (5.1) and (5.2). At this stage, we do not know precisely how the wrinkle number *m* is related to the pressurization *τ*; indeed, a key aim of our analysis is to determine more precisely how *m* and *τ* are related. We shall see that the bending stiffness does not enter the first two orders in this FT expansion.

#### (i) Order *m*

The highest order involves only the stresses in the axisymmetric base state (superscript (0)), though we shall see that this base state is *not* the same as that determined from the axisymmetric membrane theory. Vertical force balance reveals that
5.7while in-plane stress equilibrium gives
5.8In principle, one could imagine that a term due to bending, ∝*m*^{2}*τ*^{−2}, could enter (5.7). However, as we have concluded that *τ*=*O*(*m*^{2}), we neglect this term at *O*(*m*) and see immediately that (noting that *W*^{(1)}(*ρ*)≠0 by the assumption of wrinkling): at leading order the hoop stress is relaxed by wrinkling. This result is variously known as the tension-field [20], or relaxed energy [21] limit, and is well-established as the appropriate limit for wrinkled membranes. We also see from this order that .

#### (ii) Order 1

Assuming from the results of the previous order, the normal force balance at *O*(1) reads
5.9while the equilibrium of in-plane stresses requires
5.10As we saw at *O*(*m*), a term from bending ∝*τ*^{−2}*m*^{3} could, in principle, enter the normal force balance (5.9) but is neglected because *τ*=*O*(*m*^{2}). Since , equation (5.9) reveals that the r.h.s. can be balanced only by the first term of the l.h.s., i.e.
5.11Moreover, we find that since no other terms are present in the normal force balance equation. When this latter result is used in the second equation of (5.10), we see immediately that , where the constant *A*_{1} vanishes on requiring the shear stress to be continuous across the unwrinkled/wrinkled inner interface (noting that the unwrinkled domain is axisymmetric). The result, equation (5.11), is consistent with the tension-field approach for the shape of a wrinkled pressurized shell adopted previously [18]. Moreover, since is only radius dependent, the first equation of (5.10) then shows that (recall that must remain 2*π*-periodic). The first equation of (5.10) then gives the classic tension-field result .

Recalling that , the *O*(1) axisymmetric and asymmetric expansions of the stress–strain relation (2.16) read, respectively
5.12and
5.13From the first of these, we find the so-called slaving condition [22]
5.14Equation (5.13) then allows us to determine , though this is of less interest here. Note that, in contrast with the NT analysis, the wrinkle shape, *W*^{(1)}, is fixed in the FT expansion: a certain amount of length of material has to be ‘wasted’ by wrinkling and the product of the amplitude and the wavenumber is fixed by this condition. Recalling that , the *θ*-dependent expansion of (2.17) at this order reveals
5.15so that
5.16whereas the *θ*-independent expansion shows
5.17and hence that , for *C*_{1} a constant. For the hoop displacement , to match the axisymmetric state outside the wrinkled region, we must have that *C*_{1}=0, and hence that .

#### (iii) Order *m*^{−1}

Assuming from the results of the previous order, the force balance at *O*(1/*m*) reveals that
5.18Note that the term due to the bending stiffness (terms proportional to *τ*^{−2}) is finally included at this order since energy minimization led us to expect *τ*^{−2}*m*^{4}=*O*(1). The equilibrium of stresses in the other two directions gives
5.19

Recalling that , we see from the second equation of (5.19) that (to ensure 2*π*-periodicity). We therefore see that the first non-vanishing component of the hoop stress is axisymmetric, as noted previously [22]. We also find that , where *A*_{2} vanishes upon requiring the stress to vanish at the edges of the wrinkled region.

To say more about wrinkling requires a more detailed analysis of (5.18), which in turn requires that we determine more precisely . We can note that (5.18) suggests that and so we write
5.20where the radial function *g*(*ρ*) has yet to be identified. We may then write (5.18) as
5.21This has been written in such a way that the r.h.s. shows the destabilizing hoop stress (the compression that leads to wrinkling), while the l.h.s. contains the various restoring forces, including the effects of the bending stiffness (first term), a geometrical stiffness due to the curvature of the shell (second term) [23] and a tensional stiffness due to the tension along the wrinkles (third term) [15]. To make further progress, therefore, we must determine *g*(*ρ*), and to do so we note that, recalling and equation (5.16), the only *θ*-dependent terms in the *O*(1/*m*) expansion of (2.15) are
5.22immediately giving
5.23We, therefore, may rewrite (5.21) as
5.24Here we note that the second term on the l.h.s. of (5.24) may be written in dimensional terms as
5.25this term is a generalization of the curvature-induced stiffness introduced recently [26,23] for naturally flat sheets, i.e. 1/*R*=0, to objects with an intrinsic curvature.

### (b) Spatial variation of the wrinkle number

Equation (5.24) is the key result with which we can begin to understand the wrinkle patterns displayed by our ABAQUS simulations. Having generalized the effect of wrinkle curvature studied by Paulsen *et al.* [23] to include the effect of natural curvature *together* with wrinkle curvature, we shall follow their study to develop a prediction for the wrinkle number.

For any value of the wrinkle number *m*, (5.24) gives a solution for the value of that will give that wrinkle number. However, in the FT approach, we instead ask the question: what value of *m* minimizes the energy associated with this deformation? In what follows, we assume that the wrinkle number *m*(*ρ*) varies only slowly with *ρ*, and so we treat *m*(*ρ*) as the value of *m* that minimizes the energy locally.

Considering the general expression for the elastic energy of the wrinkled state, (5.6), we note that the l.h.s. of (5.24) corresponds to the functions
5.26Since the slaving condition (5.14) fixes *W*^{(1)}(*ρ*), the energy in (5.6) is to be minimized over the wrinkle number *m*, with *f*_{B}(*ρ*) and *f*_{K}(*ρ*) as in (5.26). Assuming that the wrinkle number varies only slowly with radial position *ρ* (so that the wrinkle number may be determined locally), we find that
5.27where
5.28is the curvature-induced stiffness due to the difference between the natural curvature of the shell and the curvature along the wrinkles caused by indentation, while
5.29is the stiffness due to the tension along the wrinkles [15]. Equation (5.27) is the appropriate expression of the local *λ*-law introduced by Paulsen *et al.* [23] in a setting with circular symmetry, where it is more natural to talk about wavenumber than wavelength. There are two important differences here with the local *λ*-law as previously expressed: first, our expression for *K*_{curv} (5.28) includes the effect of natural curvature (rather than imposed curvature, as studied previously [26,23]). Second, our problem does not contain a stiffness due to an elastic foundation, though this could be added to (5.27) in the obvious way.

To make the final step, we require the base state of the wrinkled shell *W*^{(0)}; the calculation of *W*^{(0)} using tension field theory, is outlined by Vella *et al.* [18]. In general, a multipoint boundary value problem must be solved numerically to find *W*^{(0)}(*ρ*). However, in the limit of Vella *et al.* [18] showed that wrinkles occupy an increasing portion of the shell, and that, within this region, the base state of the shell is given by
5.30Using this basic shape, together with (5.28), it is then easy to see that
5.31

The slaving condition (5.14) can be used with the tension-field profile for *W*^{(0)}(*ρ*), (5.30), to determine the wrinkle amplitude *W*^{(1)}(*ρ*); this has been calculated previously [18]. Using this previous result, we find that in the limit
5.32It is clear that this expression diverges both at the indenter position, *ρ*→0, and at the edge of the wrinkled zone, . Our assumption in the analysis leading to (5.27) was that variations in *m* occur only slowly with *ρ*; while this remains plausible as *ρ*→0 (since then *m*∼*ρ*^{1/4}), it is clearly not the case as . In fact, near the edge of the wrinkles there is a spatial boundary layer in which this divergence is resolved; we discuss this further in §7 but for now proceed by approximating *K*_{tens} as
5.33Combining the expression for the two stiffnesses with (5.27) the radial variation of the wavenumber is then approximately described, for , by
5.34We now compare this result with the results of our finite-element analysis.

## 6. Comparison between ‘Far from threshold’ analysis and ABAQUS

The preceding analysis resulted in a (relatively) simple prediction for the spatial variation of the number of wrinkles, *m*(*ρ*), in the limit of very large indentation depths . The scaling behaviour *m*∼*τ*^{1/2} has been predicted previously [17,18]. However, we emphasize that the dependence on indentation depth and radial position, *ρ*, are novel. Furthermore, the presence of two different stiffnesses, *K*_{curv} and *K*_{tens}, together with the very different behaviours of these two stiffnesses with radial position and indentation depth, , mean that different scaling laws may hold in different regions of the shell. In particular, we see that the two tensions balance one another when
6.1Then, we expect that for , *K*_{tens} dominates and the wrinkles are ‘tensional’; for , it is *K*_{curv} that dominates and the wrinkles are ‘curvature-dominated’. Finally, we note that wrinkles only occupy the region and that the wrinkled shell adopts a universal shape, given by (5.30), in this case. It is, therefore, natural to rescale *ρ* to exploit this universal shape, which we do by letting . In this rescaled coordinate, wrinkling in the majority of the shell is dominated by *K*_{curv} (since as ) and so we have
6.2where the dependence on Poisson’s ratio *ν* emerges from the ratio of the bending and stretching moduli, *B*/*Et*=*t*^{2}/12(1−*ν*^{2}). We note that written in this way, the wrinkle number is independent of the modulus and the pressure, depending only on the geometrical properties of the shell and the deformation. This result is somewhat surprising but has also been observed in a thin membrane floating on a liquid [23].

A particular instance of (6.2) is that the wrinkle number at the edge, *m*=*m*_{edge}, may be found by substituting *ξ*=1.19, to give
6.3We note that this corresponds to a constant wavelength at the wrinkles’ edge as indentation progresses. Results for other radial positions are summarized in table 1. We note, in particular, that when viewed at a fixed radial position, *ρ*_{fix}, the number of wrinkles observed initially decreases (according to (6.2)), while for very large indentation depths it will increase once .

Figure 6 presents a comparison for the wrinkle numbers predicted theoretically and those observed in our finite-element simulations. In particular, figure 6*a* investigates the effect of varying while holding *τ*=188 fixed, while figure 6*b* investigates the effect of varying *τ* while holding fixed. In each case, the raw dimensionless data showing the variation of *m* against the radial coordinate *ρ* is presented in the inset, while the main figure shows the data normalized according to (6.2).

The agreement between theory and numerics shown in figure 6 is generally very good. However, this comparison also highlights that the agreement begins to break down at the largest indentation depths , possibly in a manner similar to the crumpling observed in a film on a liquid drop [34]. Indeed, we note that at the external boundary of the wrinkled zone, the simulation results shown in figure 4*d* are highly reminiscent of the crumples observed in the sheet-on-drop problem (see e.g. fig. 3*c* of King *et al.* [34]).

Table 1 summarizes four different scaling behaviours, depending on which region of the shell we examine. Figure 7 shows tests of predictions for the region in which *K*_{curv} dominates (as increases, this region fills an increasing proportion of the shell). We emphasize that apparently complex behaviour can observed as increases: for a given radial coordinate *ρ*, the wrinkle wavenumber decreases with but if we instead focus on a fixed universal position, fixed, we see that *m* increases.

Despite the different spatial variations of wrinkle number, we always find that *m*∼*τ*^{1/2}, which is the typical scaling found in FT analyses [22,23]. It is particularly important to emphasize that this scaling law is different from the prediction of the NT analysis, which yielded *m*∼*τ*^{2/3} and is reproduced as the black curve in figure 7*b*. Our finite-element simulations show both the NT regime and the FT regime; we believe that this is the first time both of these have been recovered in finite-element simulations, though both regions have been studied in the theoretical analysis of Davidovitch *et al.* [33]. Indeed, from a numerical point of view, the comparison provided in Taylor *et al*. [31] suggests that different numerical methods are better suited to either the NT or the FT regimes, not both. We also note that, in our results, it seems that the transition between the two regimes is not clear-cut.

It is interesting to note that the scaling of the wavenumber with *τ*, i.e. *m*∼*τ*^{1/2}, has already been suggested in two previous works [17,18]. Vella *et al.* [17] used a scaling argument along the lines of [15] to explain the wrinkle number observed at onset (the NT behaviour). Such an argument should only be used to determine the FT behaviour, however; the apparently good agreement between the two is, we believe, an accident of the fact that the wrinkle number at onset only reaches the expected *τ*^{2/3} scaling when while, for the values investigated previously [17], the NT behaviour was still transitioning between the *τ*≪1 and *τ*≫1 regimes. A later consideration of the global energy consistent with the FT approach was given by Vella *et al.* [18]. This result predicted some variation in wrinkle wavenumber with indentation depth (though did not observe such variation), but did not take into account the spatial variation of the wrinkle pattern. As we have seen in this work, the radial variation means that there are no simple scalings for the evolution of the wrinkle pattern with applied indentation: one can see that the wrinkle number increases or decreases with indentation depth depending on whether one remains at fixed radial position, or moves with the edge of the wrinkled zone.

## 7. Conclusion

In this paper, we have presented a combined numerical and asymptotic study of the well-developed wrinkling of a pressurized shell subject to point indentation. We have shown that the scaling for the wrinkle number depends on the dimensionless pressurization *τ* according to *m*∼*τ*^{1/2}; a key result of our analysis is that this FT scaling is different to the scaling *m*∼*τ*^{2/3} that is observed at the onset of instability. Furthermore, we have shown that at a fixed indentation depth, the wrinkle number varies spatially within the shell, according to which of two effective stiffnesses (one dominated by tension, the other by curvature) is more significant. Finally, we generalized the concept of a curvature-induced stiffness proposed by Paulsen *et al.* [23] to situations in which this curvature occurs together with an intrinsic curvature.

Generally speaking, we found very good agreement between finite-element simulations and our asymptotic analysis. Nevertheless, a number of features of our results warrant further investigation. One such feature is the discrepancy found for very large indentations at the outer boundary of the wrinkled area (the crumples shown in the snapshots of deformed shapes shown in figure 4*d*, together with the decrease in wrinkle number observed there, compared with the FT prediction). We note that while strikingly similar crumples have been observed experimentally in related systems [34], one possible confounding factor here is our use of shallow shell theory. Shallow shell theory requires that the horizontal scale of the deformation is much smaller than the radius of curvature of the shell itself. In the largest indentation depths probed here (, ℓ_{p}≈5 cm) the extent of the wrinkled region, *r*_{O}≈0.5 m, which is close to the scale of the shell itself (recall that here *L*=0.8 m). As a result, future work should focus on understanding whether this drop in wrinkle number at large indentation depths is a finite size effect, or a robust feature due to crumpling. The numerical simulations also show a second interesting behaviour. At large indentations, for example or in figure 5, some wrinkles are visually more apparent than others. This is reminiscent of the wrinkle-to-fold transition that has already been noted in related problems, such as the confinement of a thin film over a liquid substrate [35]. This transition is still to be fully understood in the relatively simple two-dimensional setting of a sheet on a liquid bath, and so is beyond the scope of the present work; nevertheless, this remains an interesting possibility for future work.

Another pressing theoretical concern exposed by our analysis is the apparent divergence of the tensional stiffness at the edge of the wrinkled zone, . Simple experiments with a beachball (figure 8) confirm our finding from finite-element simulations that the wrinkle number does not diverge at the edge of the wrinkled zone and so we rule out the presence of a cascade of wrinkles as seen at the boundaries of wrinkled zones in related problems [36,37]. From (5.29) it is clear that the divergence in wrinkle number is caused by the divergence in *K*_{tens} that occurs when *W*^{(1)}→0 with d*W*^{(1)}/d*ρ* finite. In reality, this divergence is regularized by accounting for the term (d*W*^{(1)}/d*ρ*)^{2} in the computation of the slaving condition, as shown in a similar problem [33]. (This term is neglected in our analysis because of our assumption that radial derivatives are *O*(1) while azimuthal derivatives are *O*(*m*). This assumption breaks down in a boundary layer near the outer edge of the wrinkles.) For our purposes, we note merely that the existence of this smoothing out justifies our assumption that *K*_{tens} can be neglected near the edge of the wrinkled zone (as was also observed in the Lamé problem [31]). In turn, our neglect of *K*_{tens} at the edge leads to the prediction of a constant wrinkle wavelength at the edge of the wrinkled zone, in qualitative agreement with figure 8.

Finally, we note that much of the experimental interest in understanding well-developed wrinkle patterns, particularly the wrinkle number, has come from the desire to use wrinkles as a metrological assay for determining the properties of thin sheets [14,38]. Our results here, particularly the prediction of a curvature-induced stiffness that depends on both natural and imposed curvature, may lead to a reassessment of previous experimental work that used the wrinkling of deflated bubbles to infer the two-dimensional modulus of complex interfaces [38]. Our observation of a spatially variable wrinkle wavelength also suggest a possible mechanism for ‘chirp’ in photonics applications.

## Data accessibility

Data associated with this paper may be found at https://doi.org/10.5287/bodleian:DOeqDMwqG.

## Author' contributions

M.T. carried out the numerical analysis. D.V. conceived of and designed the study. Both authors drafted and approved the manuscript.

## Competing interests

We declare that we have no competing interests.

## Funding

The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 Programme/ERC grant agreement no. 637334 (D.V.).

## Acknowledgements

We are grateful to Benny Davidovitch for many discussions about the wrinkling of pressurized shells, and for thoughtful comments on an earlier draft of this paper.

## Footnotes

One contribution of 13 to a theme issue ‘Patterning through instabilities in complex media: theory and applications.’

- Accepted December 20, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.