## Abstract

This work addresses a generalization of Dean's classical problem, which sought to explain how an annular thin elastic plate buckles under uniform shearing forces applied around its edges. We adapt the original setting by assuming that the outer edge is radially stretched while the inner rim undergoes in-plane rotation through some small angle. Boundary-layer methods are used to investigate analytically the deformation pattern which is set up and localized around the inner hole when this angle reaches a well-defined critical wrinkling value. Linear stability theory enables us to identify both the critical load and the preferred number of wrinkles appearing in the deformed configuration. Our asymptotic results are compared with a number of direct numerical simulations.

## 1. Introduction

Shear-induced buckling is a form of structural instability whose relevance to the failure of rectangular plates was first explained by Timoshenko (1921) and Southwell & Skan (1924). Shortly afterwards, Dean (1924) showed how shear-induced buckling can occur in an annular thin plate whose boundaries are subjected to uniform but opposite shearing forces. The key to Dean's classical work was a simplification of the pre-buckling stress distribution that led to an ordinary differential buckling equation of Eulerian type, and hence solvable in closed form. The conclusion of his paper makes it clear that Dean realized that his analysis could not be extended for when tension is applied on one (or both) of the circular boundaries. A review of the vast literature on elastic instabilities appearing since suggests that the first attempt to revisit the problem left open by Dean was some preliminary work described in Coman & Bassom (2007*b*). There we supposed that the annulus was stretched by applying a uniform displacement field to the outer boundary while the inner rim was subjected to a small torque. This situation falls within a broader context related to *partial wrinkling*/*buckling* instabilities due to non-compressive edge loads, a situation which has attracted considerable interest over the past few decades.

There are a number of practical instances where non-compressive edge loading heralds the onset of localized and regular buckling. Examples include the azimuthal shearing of stretched annular thin plates (Li & Steigmann 1993; Miyamura 2000), the in-plane bending of a uniformly stretched rectangular thin plate (Stein & Hedgepeth 1961) and annular-shaped plates, which are uniformly stretched by imposed displacements on the two boundaries (Haughton & McKay 1995). Traditionally, such problems have been analysed using various tension field theories, which share the property that they all disregard the obvious bending stiffness involved. The seminal work of Steigmann (1990) has proved to be particularly versatile for handling a large number of interesting situations (see, for example, Haseganu & Steigmann (1994); Li & Steigmann (1995) or Roxburgh *et al.* (1995)). Tension field theory was formulated with the objective of determining the extent of the wrinkled region and the direction of the wrinkles themselves; the price of tackling these relatively large-scale questions is that the details of the fine structure of the wrinkling pattern are lost. In many cases, the theory has been shown to be a very good approximation to reality in the advanced post-buckling regime, although recently Iwasa *et al.* (2004) have questioned whether this is always so. Most of the problems for which the tension field approach is relevant have the particular feature that there is no ‘branching’ in the wrinkling process as the load increases; although wrinkles do grow in amplitude and lengthen, their number does not change.

The genesis of the tension field theory can be traced to the attempts of the German engineer Wagner in the 1920s to account for the strength of thin metal webs and spars, which can carry a shear load well in excess of the initial buckling threshold. Since it is an approximate theory valid for severe loading conditions, it is likely to be poor for describing the weakly nonlinear regime that sets in immediately after wrinkling is first triggered. This is not a major handicap, though, as the incipient behaviour is very well suited for classical buckling/post-buckling analyses in which the bending rigidity is fully accounted for. It is such a linearized, buckling regime which will be the main focus of our study in what follows.

Recent developments in cell biomechanics (Burton & Taylor 1997; Boal 2002; Bernal *et al.* 2007) suggest that the quantitative analysis of wrinkled patterns produced by living cells crawling on polymer nanomembranes can be applied to give an estimate of the force applied by the cell cytoskeleton. This requires a delicate modelling that must clearly take into account the bending rigidity of material surfaces. Using classical plate theory, Géminard *et al.* (2004) addressed several numerical and experimental aspects of this problem for a simplified model consisting of an annular thin film with prescribed displacements on both boundaries. Coman & Haughton (2006) and Coman & Bassom (2007*a*) extended the Géminard *et al.* (2004) study by investigating the possible relevance of singular perturbation methods and showed how systematic analytical progress can be made. The developments pursued in these two works have proved to have a broader scope and subsequently motivated further studies (Coman 2007; Coman & Bassom in press).

It is helpful to review briefly the rationale behind some of these papers. As we have mentioned already, when an elastic film is acted upon by suitable edge loads, it may become susceptible to a wrinkle instability. If, for the moment, we allow ourselves to lapse into relative vagueness, and we shall tighten up our meanings in due course, the articles mentioned consider situations in which a loading parameter, say *λ*, is related to the number of wrinkles generated, *n*∈. The form of *λ*=*λ*(*n*) can be determined via the numerical solution of a suitable ordinary differential eigenproblem, and, for any particular film, it would be expected that there is a preferred wrinkled configuration. In other words, of all the possible *n*, there is one for which the corresponding *λ* is the smallest; this mode number then requires least loading to excite it and is likely to be the one provoked first in practice. It can be shown that, in the limit where the thickness of the elastic film tends to zero, the underlying ordinary differential problem becomes increasingly singular and numerical experiments predict that the preferred wrinkling number grows. Moreover, asymptotic analysis then becomes viable and the first theories in this direction appealed to standard WKB arguments. Coman & Haughton (2006) demonstrated that WKB methods can provide accurate approximations to the form of *λ*(*n*) when *n*≫1, but it is not easy to use their results to infer the identity of the most-favoured wrinkling mode. In contrast, the underlying differential system can also be studied using boundary-layer theory and matched asymptotic expansions. While such techniques are commonplace in many fields of fluid mechanics, their power and versatility seem to have been rather less recognized for problems involving the stability of solids. In Coman & Bassom (2007*a*), we demonstrated how these methods can be used to determine the preferred wrinkle mode for a prestressed annular thin film in tension. Although the analysis necessary to isolate this instability pattern can be carried out relatively quickly, it was found that the resulting predictions are in excellent agreement with numerical simulations. The picture that has emerged from these studies looks seductively general: in each of the examples studied, wrinkling has been shown to be controlled by a hierarchy of boundary-layer problems each solvable in terms of Airy functions. In passing, we remark that it should not be concluded that the boundary-layer approach is superior to the WKB results for in many ways the two techniques yield complementary, but distinct, information. The relative merits of the two techniques have been discussed at length in Coman & Bassom (2007*a*, in press), hence will not be repeated here.

Here, our specific concern is with a description of the wrinkling mode that can be initiated by the circular shearing of an annular thin film. The relevant eigenproblem, which is presented in §2, has already been the subject of a numerical study by the present authors. The simulations show that a torque applied to the inner rim of the annulus triggers a localized wrinkle mode concentrated in the vicinity of this inner edge, and this problem was analysed using WKB methods in Coman & Bassom (2007*b*). For small wrinkle wavelengths, it was found that the WKB predictions correlate well with the numerical results, but it was far from clear how the preferred wrinkle pattern could be isolated; this is tackled here via the boundary-layer approach. Moreover, the analysis described in §2 is rather more than a trivial modification of studies like Coman & Bassom (2007*a*) for reasons we describe later.

We proceed as follows. Section 2 provides the background to the partial wrinkling problem that forms the main part of our study. Based on the linearized Donnell–von Kármán bifurcation equation for thin plates, an eigenproblem governed by a partial differential equation is derived. With the help of a normal-mode solution, this is then recast as a singularly perturbed Sturm–Liouville problem for an ordinary differential equation with complex-valued coefficients. The main part of our analysis is detailed in §3, where a boundary-layer technique is used to obtain the relationship between the number of wrinkles and a dimensionless form of the bending stiffness of the plate. The attractive feature of the boundary-layer method adopted is that approximations for the critical wrinkling load and the asymptotic structure of the corresponding localized eigenmodes are obtained naturally along the way. Comparisons with direct numerical simulations are included in §4, and we conclude in §5 with a summary and brief discussion.

## 2. The annular model

The model adopted in this paper was introduced and discussed in Coman & Bassom (2007*c*), so here we restrict ourselves to highlighting only the main features. The general setting is depicted in figure 1*a*: a clamped annular film of inner radius *R*_{1}, outer radius *R*_{2} and thickness *h* (*h*/*R*_{2}≪1) is stretched by imposing the uniform displacement field *U*_{0}>0 around the outer edge while the inner boundary is rotated through some (small) angle by the application of a torque *M*. Classical plate theory is used to describe the statics of the thin film and the notation used is standard.

Assuming an axisymmetric deformation prior to the onset of instability, the pre-bifurcation state of stress is easily deduced by solving the system of equations for plane stress elasticity. This information is then coupled to the linearized Donnell–von Kármán buckling equation, and so leads to the governing eigenproblem in the usual (*r*, *θ*) polar coordinates for the (infinitesimal) out-of-plane film displacement *w*≡*w*(*r*, *θ*). This system can be written aswhere the dimensionless radial co-ordinate *ρ*≡*r*/*R*_{2} and denotes the aspect ratio of the annulus. Here we have also introduced the differential operatorstogether with the two auxiliary constants(2.1)where *ν* is the Poisson ratio of the material. The analysis that leads to (_{μ}), and which is detailed in Coman & Bassom (2007*b*), also identifies two significant combinations of physical parameters, which play key roles in the description of the wrinkling instability. In particular,(2.2)where *E* is the appropriate Young's modulus. We see that *λ* is a dimensionless quantity proportional to the ratio between the shear stress and the initial tension in the annulus, while *μ* measures the dominance of membrane action over bending effects. The assumptions used tacitly in the asymptotic analysis included in §3 are: (i) *η*=_{S}(1) and (ii) *μ*≫1. The second assumption is quite natural as *μ*∝*R*_{2}/*h* and reflects our interest in the thin-film limit of the configuration described in figure 1*a*. We mention in passing that Dean's original work and a related numerical contribution by Bucciarelli (1969) dealt with the complementary case *μ*→0.

If (*r*, *λ*, *η*, *ν*) denotes the second-order plane stress tensor that characterizes the pre-bifurcation state due to the applied loads, it is known that this tensor has exactly one negative eigenvalue within the annular region(2.3)as long as(2.4)The expression that appears on the right-hand side of the inequality (2.4) represents the loading parameter threshold that marks the onset of compressive stresses in the film; in the case of a true membrane (i.e. *μ*=∞), this corresponds to the critical wrinkling load. For the sake of completeness, a sketch of the pre-buckling principal stress distribution in the annular domain is included in figure 1*b*.

Solutions of the eigenproblem (_{μ}) are sought in the form(2.5)with the understanding that the real part of (2.5) represents the physical quantity, and where *n*∈ is the mode number (the number of identical half-waves of the wrinkling pattern). Using the normal-mode ansatz (2.5), it is found that the complex amplitude *W*≡*W*(*ρ*) satisfies the differential equation(2.6)whereand Equation (2.6) is supplemented with the boundary conditions(2.7)which follow directly from the original form of (_{μ}).

It is the system (2.6) and (2.7) that was tackled numerically in Coman & Bassom (2007*b*). A sample of those results is illustrated in figure 2, which shows the form of the eigenvalue *λ* as a function of mode number *n* and aspect ratio *η*. The graphs displayed are quite representative of the general behaviour of *λ* as *n*∈ varies. If we think of some chosen *η*, then for smallish *n* the corresponding *λ* is relatively large. As *n* grows, *λ* falls before reaching a minimum value, and then it rises again with further increase of *n*. In a physical situation, it is probable that the mode number with corresponding least *λ* is the most susceptible to excitation, on the obvious grounds that it requires the smallest applied torque for wrinkling to be induced.

The stability characteristics of elastic plates are usually deduced from the plots of the neutral stability curves akin to those shown in figure 2. (The terminology is borrowed from the literature on hydrodynamics, but has a slightly different connotation in the present context.) Although the numerical solution of the normal-mode system (2.6) and (2.7) is not difficult (but the problem does become increasingly tricky once *μ*⪆500), there remains the task of identifying that value of *n*∈ for which *λ* is least; in other words, we need to find the envelope of the neutral stability curves. While WKB methods can be used to approximate this envelope when *μ*≫1, and can lead to really quite accurate results for minimal effort (see Coman & Bassom 2007*b*), it is not at all obvious how one can adapt such analysis to infer the identity of the most dangerous azimuthal mode. The results discussed in Coman & Bassom (2007*b*) suggest that as *μ*→∞, the mode number of the most important wrinkling mode grows, the corresponding eigenvalue *λ* appears to tend to a constant, and the eigenfunction of the system (2.6) and (2.7) seems to be concentrated near the inner rim *ρ*=*η*. The objective of the rest of this paper is to put these statements on a firmer footing, and we do so using boundary-layer analysis combined with asymptotic matching.

Before we launch into a detailed discussion of the properties of eigensolutions of (2.6) and (2.7), it is worth reconsidering some of the findings described in Coman & Bassom (2007*a*). The eigenproblem under discussion there, that describing the wrinkle instability of a purely tensioned annulus, has, at first glance, a very close parallel to (2.6). In fact, the only clear distinction is that the imaginary parts of the coefficients in (*ρ*) and (*ρ*) of (2.6) are absent. (We remark that the suppression of these terms would seem to remove the eigenvalue *λ* completely from the eigensystem! However, in the pure tension problem, the corresponding forms of the auxiliary constants *A* and *B* defined in (2.1) are in fact functions of the eigenvalue.) In Coman & Bassom (2007*a*) it was shown that, when the geometrical parameter *μ*≫1, the most important modes reside in the scaling regime *n*=(*μ*^{3/4}). For such values *n*∈, the corresponding eigenfunctions are compressed into a thin region of depth (*μ*^{−1/2}) attached to *ρ*=*η*. Here, the leading-order solutions are governed by a scaled form of an Airy-like equation, which can be solved so that the eigenfunction vanishes away from the inner rim of the annulus and satisfies the first equation of the boundary conditions (2.7). However, the derivative parts of (2.7) cannot be immediately enforced on *ρ*=*η*, and this points to the presence of a second layer embedded within the main (*μ*^{−1/2})-depth layer. With this solution structure, it is quite straightforward to deduce the form of the eigenvalue *λ*∼_{S}(1) as a function of *n* and see that the *n*∼*μ*^{3/4} regime contains the most dangerous mode sought.

Given the apparent close similarity between the problem in Coman & Bassom (2007*a*) and that discussed here, it would seem reasonable to assume that the underlying solution structure detailed earlier should transfer across with minimal difficulty; at worst, only a slight reworking of that theory, as suggested by the various minor modifications, is needed. Unfortunately, it soon becomes clear that the inclusion of the imaginary terms in the coefficients (*ρ*) and (*ρ*) of (2.6) is rather more than small perturbations and, indeed, if we naively adopt the scalings suggested by the paper mentioned, then these terms are actually the dominant ones in the underlying equation. The presence of these relatively large extra terms suggests that the eigenfunctions of the problem oscillate quickly on a short scale: this short scale is much less than the principal (*μ*^{−1/2}) Airy-equation layer identified in Coman & Bassom (2007*a*, in press)—although it is long compared with the extent of the embedded layer. Rapidly oscillating eigenfunctions are consistent with the numerical simulations discussed in Coman & Bassom (2007*b*; see, in particular, fig. 2) and means that the asymptotic description of wrinkles when *μ*≫1 is somewhat more than a trivial extension of calculations performed before. With these various comments in mind, we now embark on a detailed analysis of (2.6) and (2.7).

## 3. The boundary-layer analysis

The rapid spatial oscillations in eigensolutions of our problem leave us with the potential difficulty of matching together solutions over three distinct length scales: the principal region attached to the inner rim; its small embedded sub-zone; and the intermediate scale of the oscillations. Rather than try to account for all of these simultaneously, it appears easiest to strip out the spatial oscillations right from the outset and search for solutions of the particular type(3.1)with *κ*∈, a constant to be tied down in due course. The cost of writing the solution in this way is that the new equation for the (slowly varying) amplitude *V*(*ρ*) is rather more involved than that for *W*(*ρ*). Nevertheless, it can be arranged conveniently in a form similar to (2.6), so that(3.2)where the more complicated coefficients , , and are listed in appendix A. We note that the required boundary conditions for (3.2) follow immediately from (2.7), so that *V*(*η*)=*V*(1)=0 and *V*′(*η*)=*V*′(1)=0.

Guided by the analysis discussed at the end of §2, we seek the form of the solution in an (*μ*^{−1/2}) neighbourhood of *ρ*=*η*, where the coordinate *X* is defined byThe size of this region, together with the (*μ*^{−3/4}) length scale on which the eigenfunction oscillates, suggests that we seek a solution in which *V*(*ρ*), the eigenvalue *λ* and the mode number *n* all expand in series of the types(3.3a)(3.3b)and(3.3c)here, the expectation that *λ*=_{S}(1) and *n*=(*μ*^{3/4}) at leading orders is suggested in the work of Coman & Bassom (2007*a*). On substituting these expressions into (3.2), we find at the leading (zeroth) order that non-trivial solutions for *V*_{0}(*X*) can exist only if a certain algebraic constraint between *N*_{0}, *λ*_{0} and *κ* holds. First-order terms throw up a simple first-order linear differential equation for *V*_{0}(*X*); it is possible to obtain bounded non-trivial solutions only if the coefficients of the *V*_{0} and d*V*_{0}/d*X* terms individually vanish. Taken together, these three requirements lead to(3.4)(3.5)and an expression for *λ*_{1}. The second of these relations can be used to express *κ* in terms of *N*_{0} and *λ*_{0}, whence it follows that(3.6)It is remarked that this value is precisely the one appearing in the inequality (2.4), denoting the loading parameter threshold for an idealized membrane. It then follows that(3.7)with *κ*<0 (from (3.5)) and, in turn, *λ*_{1}=0.

It is at the following order that the equation which fixes the structure of the eigensolution appears. In particular, it is found thatwhere *Δ*_{j} (*j*=0,1,2) areThis is just a rescaled Airy equation, hence, in order to ease its solution, the dependent variable is changed according to(3.8)Elementary manipulations show that the amplitude *Q*_{0}(*X*) in (3.8) must satisfy the equationwhere(3.9)In order to obtain solutions confined to the *X*=_{S}(1) region, it is necessary that *Q*_{0}(*X*)→0 as *X*→∞. Thus, by choosing the proportionality constants conveniently, *Q*_{0} can be given in terms of the Airy function of the first kind,(3.10)Furthermore, since the inner rim of the annulus is restrained against out-of-plane displacements, we also need , where −*ζ*_{0} is the first zero of *A*i (*ζ*_{0}≈2.3381). This relation can then be rewritten with the aid of (3.9), thus(3.11)Now, we see the way in which the minimum possible loading parameter can be determined. Whereas the leading order *λ*_{0} satisfies (3.6) and is independent of *N*_{0}, and *λ*_{1}≡0, the correction term *λ*_{2} does vary with *N*_{0} and grows without bound as both *N*_{0}→0 and *N*_{0}→∞. It is a simple matter to check that *λ*_{2} is minimized when(3.12)and then(3.13)

### (a) Higher-order terms

Given the initial scalings for our various quantities, the analysis summarized thus far is neither particularly long winded nor complicated. Our predictions for the critical mode number and corresponding eigenvalue should be formally valid as *μ*→∞, but we would be fortunate indeed if these effectively leading-order expressions proved to be accurate approximations for all but very large *μ*. In order to extend the probable usefulness of our results, it would be convenient to have a little more information at hand. Although it can be anticipated that the associated algebraic manipulations will become increasingly tedious as we move to further orders, modern symbolic algebra tools make the task quite tractable.

#### (i) Third order

If we seek a solution with(3.14)then expressed in terms of *Z*=*α*^{1/3}*X*, the amplitude *Q*_{1} is found to satisfy(3.15)where the constants are defined according toProperties of the Airy operator can be used to find the particular integral of (3.15) analytically. Routine calculations show thatand the boundary conditions for *V* require that *Q*_{1}(*Z*)→0 as *Z*→0. Hence, we deduce that(3.16)This is a complex-valued constraint, which, when simplified, shows that at critical conditions the correction termNote also that when , the coefficient of the terms proportional to *N*_{1} in (3.16) automatically vanishes so that we can derive no information at this stage about the optimum value . At first sight, since the previous order problem fixed , we might have expected that knowledge of *N*_{1} would appear here but this is not the case, solely because our interest is in the neighbourhood of the minimum of the hypersurface *λ*=*λ*(*N*_{j}) (*j*=0, 1, 2, …).

#### (ii) Fourth order

The method of solution at this stage follows the now-established pattern. If we relate *V*_{2}(*Z*) to *Q*_{2}(*X*) in the obvious way given in (3.8) and (3.14), this function satisfieswhere the numerous *c*_{ij}∈ are complicated constants, which can be expressed in terms of quantities introduced at previous orders. Despite the large number of terms, it is purely a mechanistic task to construct a particular solution for *Q*_{2}(*Z*) and it follows that this function tends to (some constant) as *Z*→0. Rather than give the general form of the constant, which is complicated and not particularly illuminating, it is worthwhile commenting on some of its properties. Given the pattern established so far, it can be anticipated that the expression for *Q*_{20} should contain the eigenvalue correction term *λ*_{4}, a quadratic function of *N*_{1} as well as various powers of *N*_{0}; this is already indicated in the above notation for the constant. In order to determine the smallest possible value for the loading, it follows that the mode number correction *N*_{1} needs to be chosen to minimize the quadratic; in fact, the details of the calculation for *Q*_{20} show that the term linear in *N*_{1} is absent so that the desired minimization occurs when we choose . Then, the substitution of the remaining critical values yields the result(3.17)where .

### (b) The embedded layer

The solutions obtained so far in the *X*=_{S}(1) layer have the property that although *Q*_{0}(*X*) and *Q*_{1}(*X*) vanish as *X*→0, *Q*_{2}(*X*) does not and, furthermore, the constraints on the derivative of *V*(*X*) at the inner rim of the annulus are not fulfilled. All these matters can be rectified by considering the form of the solution in a thin embedded layer, where we haveand where the eigenfunction is expanded as . To the order of our working we need to compute only the leading order terms in this inner layer. Routine manipulations show that satisfies the equationBy demanding on *Y*=0, the resulting solution leads to the asymptotic behaviourand then matching with the main solution gives the critical value of *λ*_{4},(3.18)

## 4. Numerical results

With the calculations outlined previously, it follows that for large *μ* the critical *λ=λ*^{*} is(4.1)where (*j*=0, 2, 4) are given by (3.6), (3.13) and (3.18). The corresponding critical value of the mode number is available directly from (3.12), so that(4.2)In order to assess the usefulness of our predictions, it is necessary to compare these asymptotics with some direct numerical simulations. The neutral stability curves for two typical cases are illustrated in figure 3, in which we show the results for *μ*=200 and 500. The behaviour depicted here is identical to that indicated in figure 2, but a few more remarks are helpful. The dashed lines correspond to those pairs (*n*, *μ*) that satisfy *n*≪(*μ*^{1/2}) or, following the terminology adopted in Coman & Bassom (2007*b*), they fall within the so-called *membrane-like regime* in which the response curves depend monotonically on both *η* and *n*. When *μ*=200, the mode numbers up to 10 lie in the membrane-like regime, but, once *μ*=500, it is only the first 14 modes that comprise this regime. The presence of a small bending rigidity in the film does not allow the number of wrinkles to grow indefinitely, and it appears that within the range (*μ*^{1/2})<*n*≪*μ* the film identifies a neutrally stable energy configuration with *n* wrinkles. Then, the film is in the *plate-like regime* which is distinguished by the continuous lines in figure 3.

The stability envelopes are indicated by the white circles, and it is clear that in both cases the asymptotic formula (4.1) faithfully reproduces the numerical simulations; moreover, as would be expected, the accuracy improves as *μ* grows.

Another view of the usefulness of the asymptotic predictions is provided in figure 4. Figure 4*a* shows the form of the critical mode number for three values of inner rim radius *η*; the solid line denotes the asymptotic result (4.2) while the squares show the outcome of direct numerical solutions. Of course, there is a slight inconsistency in our description of *n*_{cr}; in practice, this quantity must be an integer, but the asymptotic expansions do not impose this restriction. The numerical work sought to minimize the predicted loading parameter *λ* over *n*∈ and so, to facilitate a fair comparison, what is shown by the continuous line is the nearest integer function applied to the prediction (4.2). This gives rise to the staircase-like form of *n*_{cr}, which is observed to be in excellent agreement with the numerical results. We remark that although the theory can be relied upon only for formally large values of *μ*, it appears to provide surprisingly accurate predictions for even quite modest-sized *μ*.

Figure 4*b* shows how the critical load prediction compares with its actual value based on numerical solutions. While the inclusion of the term gives a reasonable approximation to *λ*^{*}, the next (*μ*^{−1}) contribution improves the agreement quite substantially. Now, especially at values of *μ* greater than a couple of hundred, the asymptotic result is really quite excellent and justifies the effort involved in the higher-order work described in §3. Of course, the accuracy deteriorates for smaller *μ*, but, if it is remembered that the asymptotic analysis was conducted under the implicit assumption that *μ*^{1/4}≫1, it is perhaps surprising that the agreement is as good as it is. The relative accuracy (RA) for the sets of data shown is best for the larger values of *η*, although this is not immediately obvious on inspection of figure 4*b*. When *μ*=500 and if *η*=0.3, then RA≈0.6 and 3.2% for the three- and two-term approximations (4.1), respectively; these results change to RA≈9.4 and 22%, respectively, when *μ*=50. To give an indication of how these results vary with *η*, it suffices to mention that the counterparts for *η*=0.1 of figure 4*a*,*b* are RA≈2.9 and 8.7%, respectively, when *μ*=500.

## 5. Remarks

In this work, we have revisited a problem left open by Dean (1924), which is relevant to the buckling instability experienced by annular thin elastic plates under circular shearing. It has long been known that such plates are susceptible to a wrinkling, which takes the form of a large number of equiangular spiral wrinkles that are concentrated around the inner rim of the annulus. Figure 5 shows the type of pattern that is seen in practice; it is observed how a large number of wrinkles are focused in a small neighbourhood of the inner boundary and the individual wrinkles spiral around the centre of the hole. This spiralling effect is the physical manifestation of the intermediate-scale oscillation in the solution eigenfunction. The presence of tensile loads at the outer boundary gives rise to a situation in which the linearized bifurcation equation is not amenable to closed-form solution so that Dean's approach becomes inoperable. The existence of the natural large parameter *μ*, which is related to the initial state of prestress, means that boundary-layer methods can be used to determine both the critical load and the preferred number of wrinkles. It is surprising how the computation of the first few terms in the asymptotic series for the loading *λ* and the number of wrinkles *n* give rise to expressions which are in very good accord with direct numerical simulations. We have already mentioned that, in theory at least, symbolic algebra methods combined with matching could be used to calculate further terms in (4.1) and (4.2) but, given the impressive accuracy of the existing terms as evidenced by the results in figure 4, it is by no means certain that such extra work would be commensurate with the gain.

The situation explored has the tantalizing prospect that it could be related to other physical scenarios. Lindsay (1992) and Haughton & Lindsay (1993) investigated the deformation of a hyperelastic annular slab caused by the rotation of a rigid shaft passing through the hole and bonded to the annulus. Assuming that the lateral curved surface is kept fixed, a second-order elasticity theory was used to describe the deformation of the top surface of the slab. For a very broad class of constitutive equations, a ‘pinching’ deformation in the vicinity of the shaft was noted, and fig. 3 of Haughton & Lindsay (1993) reveals a boundary-layer behaviour akin to that in the present work. Those authors conducted an axisymmetric analysis and it would be of significant interest to incorporate azimuthal terms. This would both increase substantially the range of possible deformations and enable us to determine which is preferred in practice. Moreover, the studies of Lindsay (1992) and Haughton & Lindsay (1993) demonstrate, to a certain extent, a finite elasticity version of the classical *Weissenberg effect* for highly viscous fluids (e.g. Reiner 1958). In its simplest form, this phenomenon appears when a rod, which is rotated in a non-Newtonian fluid, leads to the latter climbing up the rod; a comprehensive mathematical analysis was given by Joseph & Fosdick (1973) and Joseph *et al.* (1973). We intend to return later to probe this apparently unexplored link between our work described here and these other mathematical problems.

## Acknowledgements

C.D.C. acknowledges with gratitude the financial support received from the Royal Society for a visit to the University of Western Australia. We would like to thank the referees for their comments on this work.

## Footnotes

- Received June 25, 2007.
- Accepted August 13, 2007.

- © 2007 The Royal Society