## Abstract

An asymptotic study of the wrinkling of a pressurized circular thin film is performed. The corresponding boundary-value problem is described by two non-dimensional parameters; a background tension *μ* and the applied loading *P*. Previous numerical studies of the same configuration have shown that *P* tends to be large, and this fact is exploited here in the derivation of asymptotic descriptions of the elastic bifurcation phenomena. Two limiting cases are considered; in the first, the background tension is modest, while the second deals with the situation when it is large. In both instances, it is shown how the wrinkling is confined to a relatively narrow zone near the rim of the thin film, but the mechanisms driving the bifurcation are different. In the first scenario, the wrinkles are confined to a region which, though close to the rim, is asymptotically separate from it. By contrast, when *μ* is larger, the wrinkling is within a zone that is attached to the rim. Predictions are made for the value of the applied loading *P* necessary to generate wrinkling, as well as details of the corresponding wrinkling pattern, and these asymptotic results are compared to some direct numerical simulations of the original boundary-value problem.

## 1. Introduction

Many traditional buckling instabilities for thin elastic plates can be described by linear boundary-value problems that are easily amenable to analytical solutions (e.g. [1,2]). This is particularly true in the case of in-plane compressive loads when adopting the Kirchhoff–Love hypothesis and using inextensional plate theory. One of the main reasons for the lack of difficulties can be traced back to the very simple nature of the pre-buckling (or base state) which, more often than not, is simply constant. However, buckling becomes likely whenever compressive stresses develop in some region of the plate, and this does not require the presence of any compressive external loads. Indeed, this was the situation in some recent works [3,4] where the concerted effect of inhomogeneous pre-buckling deformations and tensile loads led to localized wrinkling patterns (when global tensile loads are involved, it is customary to use the words ‘buckling’ and ‘wrinkling’ somewhat interchangeably). Although the resulting linear equations were not solvable in closed form, the presence of initial tension introduced a large non-dimensional parameter that played a key role in obtaining analytic formulae for characterizing both the critical load and the buckling pattern. As explained shortly, the situation that we want to address in this study is somewhat different, in that the pre-buckling deformation will be described by a nonlinear boundary-value problem whose solution is not available explicitly.

The problem of interest here was inspired by a familiar situation encountered in the food industry. Figure 1 shows some pictures that illustrate the instability phenomenon that we have in mind. The upper picture contains a bird-eye's view of a cylindrical plastic food container, while closer snapshots can be seen in the other two photographs, in the lower part of the same figure. The container is covered by a stretched thin and transparent polyethylene lid that is glued to the circular rim of the cylinder. In this case, the product was packaged in a protective atmosphere, a process that involves flushing out the air from the container and replacing it with a mixture of different gases; as a result of this process, a negative pressure gradient is created inside the plastic chamber. In a first approximation, this scenario is equivalent to a uniformly stretched circular thin film subjected to a uniform pressure that acts perpendicular to its surface. The almost regular edge-wrinkling pattern that develops along the rim of the lid is clearly visible in figure 1. Such occurrence is undesirable, as it often indicates a high stress concentration and involves significant bending stresses that could eventually lead to loss of integrity in the packaged product.

Panov & Feodosiev [5] were the first to examine the possibility of buckling of pressurized thin plates using a special theory for the elastic buckling of shells and elementary Galerkin-type methods. Since then, there have been several related numerical studies that have been reviewed elsewhere [6], but for our purposes the most important is that by Adams [7], who considered a thin elastic plate uniformly stretched in the radial direction before a transverse load is applied. In his formulation, a nonlinear axisymmetric pre-bifurcation is assumed, and the Föppl–von Kármán equations are then solved for a concentrated load applied at the centre of the plate. The dimensionless equations contain just three parameters; Poisson's ratio of the elastic material, a re-scaled background tension and the applied load. Adams [7] found that buckling patterns tended to occur in the outer region of the disc, well away from the centre, and that as the background tension was increased so this effect became increasingly pronounced. The mechanisms driving this localization behaviour were not examined, which led Coman [6] to consider the related problem in which the concentrated load is replaced by a uniform pressure applied over the upper face of the disc. In that paper, a set of co-ordinate-free bifurcation equations were derived that govern the stability of an axisymmetric and nonlinear state; the calculation of the undisturbed fields and the non-axisymmetric disturbance equations required the numerical solution of a system of coupled ordinary differential equations of order 12. The localization of eigenmodes was clearly seen, and some of the main results are reproduced in figure 2.

The equations solved by Coman [6] contained three important parameters: the applied pressure *P*, a dimensionless measure of the background tension *μ*, and the wavenumber *m* of the buckling pattern. Although the results given in [6] are mostly computational, it was recognized that a careful analysis of the equations might provide further insight in to the mechanisms at work. As an example, it is immediately seen that many of the results summarized in figure 2 apply at large values of *P*, so presumably advantage of this can be taken to possibly unwrap some of the complications implicit within the system of bifurcation equations. Moreover, for each value of *μ* plotted in figure 2 there is a preferred wavenumber *m*_{c} which corresponds to the least loading *P*; Coman [6] conjectured that when *P*≫1 the bending effects are confined to a thin layer of size *O*(*P*^{−1/3}), that the buckling mode is confined to a region of extent *O*(*P*^{−1/6}) and that *m*_{c}=*O*(*P*^{1/4}).

The above observations provide the motivation for this paper; our main objective is to implement asymptotic methods and thereby gain insight into and explanation for the buckling characteristics noted by the first author in [6]. In order to flesh out our results, the remainder of the paper is organized as follows. To begin, in §2 we refresh quickly the important equations solved in [6]. The principal results of our analysis are then summarized in §§3 and 4. In the former, we show how it is possible to understand the behaviour of the lower curves in figure 2, for which *μ* is relatively modest while *P* is still large. Although it is not immediately clear as to why it could be expected, it transpires that some of the important features of the upper curves can be rationally described in terms of a large-*μ* analysis, even though *μ* is still roughly

## 2. Development of the key equations

Consider a circular thin film (or thin plate) of radius *a* and thickness *h*, with *h*/*a*≪1, under an initial in-plane tension *N*_{0}=*σh* acting along its circumference, and also subject to a uniform transverse pressure load *p*>0. Here *σ*>0 is the applied constant stress acting in the mid-plane of the plate, while the linearly elastic material is characterized by Young's modulus *E* and Poisson's ratio *ν*.

If *w* and ** u** denote the transverse and in-plane displacements of the plate, respectively, then, under certain kinematic assumptions (e.g. [8]), the Euler–Lagrange equations of the plate energy functional demand that these fields should satisfy the well-known Föppl–von Kármán equations

*D*=

*Eh*

^{3}/12(1−

*ν*

^{2}) represents the bending rigidity of the plate, and

**is the (symmetric) membrane stress tensor.**

*N*For a transverse displacement *w*=*w*(*r*,*θ*) and corresponding in-plane displacement ** u**(

*r*,

*θ*)=

*u*

_{r}(

*r*,

*θ*)

*e*_{r}+

*u*

_{θ}(

*r*,

*θ*)

*e*_{θ}, expressed in terms of the usual polar coordinates (

*r*,

*θ*), routine manipulations, such as those presented by Sheplak & Dugundji [9], lead to nonlinear partial differential equations for the components of the membrane stress tensor and the transverse displacement. The key balances may be cast in the forms

*a*

*b*

*c*where ∇

^{2}≡∂

^{2}/∂

*r*

^{2}+

*r*

^{−1}(∂/∂

*r*)+

*r*

^{−2}(∂

^{2}/∂

*θ*

^{2}) denotes the usual Laplacian in polar coordinates.

### (a) The base-state equations

In the axisymmetric base state, we have *u*_{r}=*u*_{r}(*r*), *u*_{θ}=0 and *w*=*w*(*r*), which imply that *N*_{rr} and *N*_{θθ} are functions of *r* and *N*_{rθ}≡0. Written in suitable dimensionless variables [6] the Föppl–von Kármán equations (2.2) may be cast as
*a*
*b*
*c*
where
*ρ*≡*r*/*a* so that the interval of interest is *ρ*∈[0,1] and the parameters *μ* and *P* are proportional to the in-plane tension *σ* and the applied transverse pressure load *p*, respectively; more precisely

It is worth noting at this early stage the definitions of the parameters *μ* and *P* as given in equation (2.5). Although our subsequent analysis deals with situations when at least one of these is large, the presence of factors of *a*/*h*≫1 in (2.5) implies that this regime does not require the uniform transverse pressure loading *p* nor the applied stress *σ* to be particularly excessive. Indeed for many practical situations concerned with a membrane, like the packaging example alluded to in the introduction, the parameters *μ* and/or *P* are very far from being small.

Asymptotic expansions of (2.3) with various boundary conditions, in the limit *P*≫1, have been known for some time; the case *μ*=0 is extensively reviewed in Stoker [10], which contains additional historical pointers. In more recent times, both [9] and [11] have examined the situation involving *μ*≠0. It is worth mentioning that these authors were concerned with an arbitrary scenario in which *P* and *μ* were unrelated, a fact that hampered their perturbation analysis. In fact, their results are of little use in our present stability problem, as in this case the basic state and the bifurcation equations are intimately connected, and require much more precise information about the pre-bifurcation fields.

### (b) The bifurcation equations

The bifurcation equations are obtained in the usual way; an infinitesimal increment is added to the basic state, the sum substituted into system (2.2), and then linearized. If the increments are written as
*a*
*b*
*c*
where the coefficients *A*_{ijk} were defined in [6], and are listed in appendix A for the sake of convenience. It is noted that these equations are of order eight, and need to be solved subject to appropriate regularity conditions at the centre of the thin film, clamped boundary conditions in the transverse direction at its edge (that is, *w*(*r*,*θ*)=∂*w*(*r*,*θ*)/∂*r*=0 at *r*=*a*) and prescribed radial stress distributions along the outer perimeter (so *N*_{rr}(*r*,*θ*)=*N*_{rθ}(*r*,*θ*)=0 at *r*=*a*). In dimensionless form, it is found that these boundary conditions are equivalent to
*a*
*b*

The numerical solution of the problem defined by (2.3) and (2.9) subject to (2.6), (2.7) together with (2.10) was performed using standard boundary-value problem solvers available in MATLAB. None of the computations required particularly innovative implementation, although in several cases it was found most helpful to track the forms of the solution curves using homotopy methods.

## 3. Asymptotic analysis for μ = O ( 1 )

We begin our consideration of the bifurcation problem when the in-plane tension parameter *μ* is not large; the precise meaning of this statement will become apparent as the structure of the solution is developed. The results illustrated in figure 2 suggest that even when *μ* is relatively small the corresponding *P* is not, and so we commence using the assumption that *P*≫1. It is immediate that these constraints fix the solution to the base-state problem (2.3) with (2.6) and (2.7) at the first few orders in powers of *P*.

### (a) Base-state solution

The structure of the base-state equations imply that across the majority of the plate the solution assumes the form

Equation (3.2*c*) guarantees the satisfaction of the regularity condition for *h*_{0}, but also indicates that the condition at *ρ*=1 cannot be fulfilled; this points the way to a boundary-layer structure for the base state near *ρ*=1 which we develop below. The elimination of *g*_{0} and *h*_{0} from (3.2) leads to the nonlinear equation
*Λ*≡(1−*ν*^{2})/2. This is the well-known nonlinear membrane equation derived by Föppl [12]. Heuristic power-series solutions of (3.4) were given in [13], although the existence and uniqueness of solutions for various boundary-value problems governed by this singular nonlinear differential equation were only later resolved (e.g. see [14] or [15] and references therein). The emergence of (3.4) in our asymptotic scheme is not at all surprising; indeed, as pointed out in [6], even when *μ*=0, the applied pressure *P* required to trigger the wrinkling instability is large enough so that the thin film is close to a membrane state (see fig. 5 of that reference).

To ensure the proper matching between the base state where 0≤*ρ*<1 to the boundary layer, it is necessary to understand the behaviour of *f*_{0} as *ρ*→1. In the view of the anticipated boundary condition that *f*_{0}(1)=0, it is elementary to determine that as *ρ*→1 so *f*_{0}∝(1−*ρ*)^{2/3}; more precisely, in this limit
*x*≡1−*ρ* satisfies 0<*x*≪1. This result contains one slight subtlety; the correction term to the leading-order behaviour satisfies a linear equation so that the constant *K* cannot be determined by the analysis local to *ρ*=1. Rather, it can only be settled by numerical solution of (3.4) solved subject to *f*′_{0}(0)=0 (an immediate consequence of equation (3.4)), the regularity condition *f*_{0}(0)=*g*_{0}(0) and behaviour (3.5). This calculation gives *K*≈−2.1755.

We now need to examine the structure of the base state in the boundary layer. Standard scaling analysis suggests that this layer is of extent *X* by
*Λ* into the definitions of variables. As *ν*. The alternative would be to omit these terms initially, but then all the independent and dependent variables would need to be subsequently re-scaled leading to more complexity in notation.)

#### (i) Base-state solution for X = O ( 1 )

We shall see later that to the precision that we require to isolate the structure of the bifurcation solution it is necessary to consider the first couple of terms in the forms of *N*_{rr}, *N*_{θθ} and *Θ*, where *a*
*b*
*c*
here, each quantity on the right-hand sides is a function of *X*. Unfortunately, none of the variables can be determined analytically, and require numerical solution. Thus, we just record the governing equations order by order: zeroth-order terms give
*X*, and at the next order it follows that:
*F*_{j}(0)=*H*_{j}(0)=0 for *j*=0,1 and suitable far-field conditions as *F*_{j}, *j*=0,1 can be inferred from appropriate re-scaling of the asymptotic behaviour (3.5), and they in turn fix the constraints on the *G*_{j} and *H*_{j}. In the interests of brevity, we do not list all these formulae as they can be routinely deduced.

### (b) Bifurcation solutions

The calculations of [6] confirmed that the bifurcation solutions reside near the rim of the thin film, and within the boundary layer for the basic state for which *X* scaling, then the necessary balance falls into place when
*X*=*X*_{0}, where *X*_{0} is thus far unknown and needs to be determined by the equations. Thus, it is convenient to define the scaled coordinate

Given this, we can now proceed with the analysis. First, we comment on the base state. This is easily written down for within (3.7) each quantity just needs to be expanded as Taylor series about *X*=*X*_{0}; hence we define constants *G*_{j} and *H*_{j}. Once *X*_{0} is known (see below), the constants *j*=0,1 can be obtained from the numerical solution of (3.4) and (3.9). The remaining coefficients may be deduced from (3.4) and (3.9) and their differential consequences.

Now we comment upon the structure of the disturbance solutions of (2.9). It may be demonstrated that
*m* is given by

If (3.13) and (3.14) are substituted in (2.9), a hierarchy of balances are obtained. At leading orders, the three equations give
*a*
*b*
Now the first two equations tell us the form of *u*_{0} and *v*_{0} once *w*_{0} is known; for the third equation to have a non-zero solution forces the constant inside the brackets to vanish. Next-order balances in (2.9c) leads to a similar type of equation, and the upshot is that for a non-trivial *w*_{0} then
*M*_{0} and the location *X*_{0}. As the constraints depend only on the properties of the zeroth-order basic-state equations in the *M*_{0}≈3.3179 and *X*_{0}≈0.4575.

The next non-trivial terms arising from (2.9c) yield a homogeneous differential equation for *w*_{0}, namely
*a*
*b*
*c*
It is at this stage that the reason behind the selection of the scaling (3.10) becomes apparent. With this choice the driving equation for *w*_{0} is a Gaussian-like differential equation which, for appropriate values of *Υ*_{j}, admits a solution such that *w*_{0}→0 exponentially quickly as both *α* if *X*=*X*_{0} and on the constant *M*_{6}, it is just a matter of calculation to conclude that *M*_{6}≈−1.0949.

Furthermore, we can also note the absence of the boundary conditions (2.10) from our considerations here. As the eigenfunctions are contained in this thin layer away from the rim of the disc at *ρ*=1, the bifurcation functions are exponentially tiny at *ρ*=1 so that the boundary conditions are automatically fulfilled.

Once we have established the pattern for determining the subsequent equations, the only obstacle that remains is that presented by the increasing complexity of the commensurate manipulations. It should be noted that the next non-trivial terms in (2.9c) throw up an algebraic equation for *M*_{8}, which may be solved to give

### (c) Numerical comparisons

Our analysis outlined above predicts that for values of the in-plane tension *m* and the loading *P* is

In figure 3, we show a direct comparison between this prediction and the direct numerical solution of the system (2.3) and (2.9) solved subject to the boundary conditions (2.6), (2.7) and (2.10). It is clear that the agreement between the numerical solution and the prediction is quite excellent. The three-term result seems to represent a significant improvement on the simpler two-term prediction, which is a partial justification for pursuing the additional asymptotic workings as far as described here. It might seem counterintuitive that the difference between even the three-term asymptotic result and the computations grows as *m* increases but this is no more than a reflection that the three-term expression neglects some growing, albeit slowly growing, terms as *m*=60 the three-term asymptotic prediction is within about 10% of the computed value; an error which decreases to about 6% when *m*=100 and less than 2% for *m*=200.

Of interest in a practical setting is the dependence of *m*_{c} and *P*_{c} (recall that is the critical wavenumber corresponding to the smallest *P*) for various values of *μ*; these results are noted in table 1. Unfortunately, the prediction (3.20) is only strictly applicable for *P* (or *m*) large; of course, it cannot be used to estimate the value of *m*_{c} for which *P* is not particularly large. By its very definition *m* should strictly be taken as an integer while the analysis outlined above does not impose this restriction. Thus, in table 1 minimum values of *P* are calculated for integer values of the wavenumber.

Last, we point out that, to the orders evaluated here, the results, in particular (3.20), hold unchanged whatever the value of *μ*.) In order to appreciate the potential usefulness of our predictions, in figure 4 we include the forms of the neutral stability curves for various values of *μ* up to *μ*=10. For relatively small values of *μ*, the result (3.20) continues to be a very good approximation to the numerical results, but becomes increasingly poor as *μ* grows. As *μ* becomes significant, the important balances in the problem must alter, and this is the subject of our next section.

## 4. Asymptotic analysis for large *μ*

Our results discussed above suggest that at relatively modest values of *μ* the derived results become forever less accurate and points to the eventual need for a modified analysis. In order to identify when results might be significantly changed, we refer back to the base-state equations (2.3). In §3a, we proposed that when *P*≫1 the stress component *N*_{rr} is *N*_{rr} that should be enforced at the rim of the disc *ρ*=1, constraint (2.6) suggests that once *μ*^{2}∼*P*^{2/3} the underlying base state will be significantly modified across the majority of the circular configuration. Thus, when *μ* is large, we are led to an expansion for the eigenvalue *P*, which takes the form
*P*, we adopt a modified strategy henceforth. In turns out to be mathematically convenient to seek *P* and the perturbation wavenumber *m* as functions of *μ*; of course, should results be preferred in terms of *P*, they can be deduced upon inversion of (4.1). However, before we consider the structure of the perturbation, let us start by developing the form of the base state to the accuracy required.

### (a) The base-state structure when *μ*≫1

Across the majority of the plate, the scalings proposed above suggest that a solution of the base-state equations (2.3) subject to (2.6) and (2.7) be sought in the form
*a*
*b*
*c*
At zeroth order, it follows that:
*a*
*b*
*c*
This nonlinear second-order system admits a solution which fulfils all the regularity conditions at the centre of the thin film if *N*_{rr}(1)=*μ*^{2}, so that *c*) imply
*ρ*=1 in which the base state must adjust to satisfy the rim conditions. Details of this layer will be noted later, but for the current purposes it is helpful to realize that *ρ*) and *K*_{0}.

At this juncture, it is helpful to anticipate some of the structure of the bifurcation equations to come. We have already noted that *N*_{θθ} is also *N*_{θθ} so large; the only resolution is to demand that *ρ*→1. By (4.4) this requires *K*_{0}≈5.1587 so that

At the first order, we find
*a*
*b*
*c*
from which it follows that
*ρ*=0 are automatically satisfied. Using (4.4) and (4.8) to eliminate *P*_{1}/*P*_{0}) and the quantity of importance for later matching turns out to be

It is helpful to note some key properties of the second-order problem for the base state. Now
*a*
*b*
*c*
These equations, in conjunction with the lower order solutions (4.4) and (4.8), may be used to show that

### (b) The structure of the bifurcated solution

In order to understand the nature of the bifurcated solution, we recall the calculations of [6] and the suggestion that, when *P*≫1, the wrinkling is confined to an *ρ*=1. In turns out that this zone is asymptotically smaller than our *a*
*b*
*c*
Anticipation of the important wavenumber regime prompts us to define
*a*
*b*
*c*
substitution into the bifurcation equation (2.9a) gives at zeroth order just

Equation (2.9b) is automatically satisfied at leading order because it simply yields the derivative of (4.18), but the third equation (2.9c) leads to the determining equation for *W*_{0},
*W*_{0}→0 as *Y* →0, and this is ensured if
*ξ*_{0}=−2.3381 is the first zero of the Airy function Ai. If the value of (4.10) is substituted into this result it follows that:

We now have the appearance of the first wavenumber-dependent term in the expansion of *P*. It is helpful to derive just one more term in this form, and this requires consideration of the next-order equations in the bifurcation system (2.9). Routine manipulations lead to
*a*
*b*
*c*
where
*a*
*b*
*c*

Recall that the leading-order term *U*_{0} was given by (4.18), once *W*_{0} is known; at this next order, we see that the equation that determines *V* _{0} is (4.25b). Given the solution *W*_{1}(0), which we record for future reference:
*Ai*′_{0} denotes the derivative of the Airy function

At this stage, we see the need for some sub-layer structure since at this order we cannot enforce all the boundary conditions (2.10). In particular, nothing as yet has been said about conditions (2.10b) that relate *U* and *V* at the rim *ρ*=1; since *U*_{0} was given completely by (4.18), and *V* _{0} by (4.25b), there is no freedom to ensure that (2.10b) holds. Consequently, we are led to examine both the base state and the bifurcated solution in a suitable inner layer.

### (c) The inner layer

Define the inner layer by *ρ*=1−*μ*^{−1}*η* with *η*=0.

For the bifurcation problem, a leading-order solution is sought in the form

At leading orders, the edge boundary conditions for the bifurcated solution reduce to

Matching with the solution to (4.21) and (4.27) forces *C*=*Ai*′_{0} and
*Γ*_{j} functions are substituted from (4.25c), and the results (4.10) and (4.13) used, it is then just a matter of some algebra to deduce that

### (d) Comparison with computations

In the analysis above, we have shown that when the tension *μ* is considered to be large, bifurcated solutions are restricted to the outer rim of the thin film, and occur when

In table 2, we exemplify the quality of this prediction by comparison of the critical values *P*_{c} and *m*_{c} with some direct computations of the full system over a range of values of *μ*. It is seen that the agreement is very respectable, even at relatively modest values of *μ*. Of course, away from the critical points it is unrealistic to expect excellent agreement for all wavenumbers, and eventually the result (4.38) begins to fail for values of *M* that are too small or too large; indeed if *m* is taken sufficiently large, the asymptotic result (4.38) becomes physically absurd with *P*<0! This is of course just an indication of the breaking down of the analysis, and it is not difficult to show that successive terms in (4.38) become disordered once

While the reconsideration of the problem might have some mathematical curiosity, it is likely to be of marginal interest for physically motivated problems. In practice, one of the important facets of wrinkling problems is the location and number of the corrugations that are likely to be observed. Fortunately, the workings above enable us to consider these aspects with minimal additional effort. We have already noted that the radial structure of the wrinkles is determined at leading order by the solution of a scaled Airy's equation; furthermore, the result (4.38) provides a simple determination for the preferred wrinkle wavelength. This is fixed by identifying that value of *M* leading to the smallest *P*; that is,
*μ* our simple formula for the identity of the most dangerous mode of wrinkling is extremely accurate, while the results of figure 6, although not quite as good, still show that the asymptotic methods provide a useful guide for the corresponding critical loads.

## 5. Concluding remarks and discussion

In this paper, we have employed asymptotic methods to explore the wrinkling mode of instability in a pressurized circular thin film, here modelled with the help of the nonlinear Föppl–von Kármán plate equations. All our results are expressed in terms of the dimensionless parameters *μ* and *P* but any of them can be recast in terms of the in-plane tension *σ*, transverse pressure load *p* and the important ratio of the plate thickness to radius *h*/*a* using the definitions (2.5). In particular, we see that the preferred wrinkle wavenumber *m*∝(*a*/*h*)^{3/4} so increases as the film becomes progressively thinner. Our work provides an analytical underpinning for the numerical results reported in [6], confirms many of the suggestions in that paper and shows how the structure of the governing equations can be used to infer the nature of the bifurcated modes. Of particular intrigue is how useful results seem to arise over a range of values of the background tension *μ*. For relatively small values of *μ*, the bifurcation consists of a mode confined to a thin zone which, although close to the rim of the plate, is asymptotically distinct from it. The zeroth-order mode satisfies a scaled parabolic cylinder equation that can admit solutions decaying in both directions away from the centre of the structure.

For larger values of *μ*, it is expected that the studies of §3 become ever more inaccurate; what is more surprising is how the problem at relatively modest values of *μ* is well described by the analysis of §4. Again, the modes are localized in a narrow zone at the rim, but the situation contrasts with the small-*μ* case, in so far as the zone now extends right to the edge of the plate. What is also of note is how the existence of the preferred wrinkle mode depends on a subtle coupling of the values of *μ* and *P*; the mode is supposed to be of small amplitude but its presence affects the structure of the base state across the whole of the plate (we saw how this provided a third boundary condition for equation (4.5)). It is noted that, while there are many examples in the literature describing problems for which the presence of some large nonlinear mode has ramifications for the underlying basic state, it is relatively rare that the existence of a mode of infinitesimal amplitude and confined to a small region plays a central role in fixing the basic state across the whole region of interest. Given this realization however, it is then a relatively standard exercise in perturbation analysis to determine the form the loading parameter *P* and the identity of the most dangerous wrinkle pattern.

Our analysis here provides a further example of how modern asymptotic methods can be used to unravel the structure of wrinkling present in a relatively complicated set of driving equations. While extremely accurate predictions for various physical quantities can only be determined via full numerical simulations, asymptotic techniques can detail the important mechanisms at play. In the problem studied in this paper, it is remarkable as to the range of parameters over which asymptotic predictions provide valuable insight; this is quite unusual and is a welcome bonus for the potential impact of our findings.

## Data accessibility

All the calculations in this paper were conducted using standard routines available in MATLAB.

## Author' contributions

All three authors contributed to the mathematical analysis, computations and the writing of this paper. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

There was no external funding supporting this work.

## Acknowledgements

The referees are thanked for their insightful comments which led to significant improvements and presentation of the paper.

## Appendix A. Coefficients for the bifurcation equations

The expression of the coefficients *A*_{ijk} that appear in (2.9) were originally derived in Coman [6] and are reproduced below for convenience. These quantities are labelled using the following convention: the first index *i* refers to the number of the equation in which it appears, the second index *j* corresponds to the position of the unknown function in the sequence (*U*,*V*,*W*) and the index *k* indicates the order of differentiation of the unknown function *U*, *V* or *W*. The coefficients for the first equation are given by

- Received July 10, 2015.
- Accepted September 3, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.