## Abstract

The stability of the wrinkling mode experienced by a compressed half-space of neo-Hookean material is investigated using analytical and numerical methods to study the post-bifurcation behaviour of periodic solutions. It is shown that wrinkling is highly unstable owing to the nonlinear interaction among the multiple modes associated with the critical compressive state. Concomitantly, wrinkling is sensitive to exceedingly small initial imperfections that significantly reduce the compressive strain at which the instability occurs. The study provides insight into the connection between wrinkling and an alternative surface mode, the finite amplitude crease or sulcus. The shape of the critical combination of wrinkling modes has the form of an incipient crease, and a tiny initial imperfection can trigger a wrinkling instability that collapses into a crease.

## 1. Introduction

Surface instabilities are frequently observed when highly elastic soft materials are compressed (Tanaka *et al.* 1987; Gent & Cho 1999; Trujillo *et al.* 2008; Cai *et al.* submitted) and their importance has grown along with the steady increase in applications of soft materials (Crosby 2010). Biot (1963, 1965) appears to be the first to have demonstrated the existence of wrinkling instability modes at the surface of an incompressible neo-Hookean elastic half-space. These modes occur as a bifurcation from a state of uniform compression with the unusual feature that their wavelength is undetermined—the scale of wrinkle undulations is arbitrary as long as it is short compared with any other geometric dimension of the solid. Throughout this paper, the coordinate *x*_{1} is aligned with the direction of in-plane compression, *x*_{2} is aligned perpendicular to the free surface of the undeformed half-space and *x*_{3} is the out-of-plane coordinate (figure 1). The stretches in the uniform pre-bifurcation state are denoted by *λ*_{1},*λ*_{2} and *λ*_{3}. With *λ*_{3} imposed and fixed, Biot found that bifurcation into in-plane wrinkling modes occurs when the in-plane compression attains.
1.1The crease, or sulcus, surface mode, first analysed by Hohlfeld (2008) and Hohlfeld & Mahadevan (2011), is doubly unusual in that, in addition to having arbitrary wavelength, or depth, it does not emerge as a bifurcation but rather exists as a local state involving finite strain changes from the uniform compressive state. By carrying out a finite-element analysis of this state in a neo-Hookean half-space, Hong *et al.* (2009) have shown that for any fixed *λ*_{3} a crease is energetically favourable for compression in the fundamental state exceeding
1.2A crease will reduce the energy of the solid when the compressive strain exceeds equation (1.2), but the deformation pathway leading to crease formation was not determined by these authors. The mystery underlying these two modes is heightened by the fact that a crease can exist at smaller compressive (nominal) strain than that required for the onset of wrinkling, i.e.
Hohlfeld & Mahadevan (2011) explored the closing and opening pathways of a finite amplitude crease under a cycle of applied compression by attaching a very thin film with bending stiffness to the surface whose purpose is to regularize the numerical model by fixing the wavelength. As these authors emphasize, the free surface of any soft elastic solid is susceptible to wrinkling and creasing under compression because the mode wavelengths can be arbitrarily small and locally a surface will be effectively flat.

The present paper builds on the work cited above with a twofold objective: to determine the stability of wrinkling by carrying out a nonlinear post-bifurcation analysis, and to account for the role of initial imperfections in the form of slight surface undulations on the stability of the wrinkles. A clear pathway to crease formation emerges. The paper is organized into sections as follows. Section 2 develops the energy functional for the neo-Hookean half-space on which the analysis is based, and it briefly reviews Biot’s (1963, 1965); bifurcation results which form the starting point of the nonlinear analysis. In §3, Koiter’s post-bifurcation (Koiter 1945; van der Heijden 2009) approach is presented as relevant to wrinkling—the finite deformation of a nonlinear elastic solid with multiple modes associated with the critical bifurcation stress. The detailed analysis of stability and imperfection-sensitivity is executed in §4. The numerical analysis of the role of imperfections on the stability of wrinkling is presented in §5. Two types of imperfections are considered: a sinusoidal surface undulation in the shape of one of the classical wrinkling modes similar to that considered in the analytical study, and an isolated slight surface depression. The numerical analysis reveals connections between wrinkling and creasing which are summarized in the conclusions in §6.

## 2. Energy functional and the bifurcation solution

Let *x*_{i},*i*=1,3 be Cartesian coordinates defined above labelling material points in the undeformed body (figure 1). These coordinates will be used throughout the analysis. All tensor components will be referred to these coordinates. Let *u*_{i}(**x**) be the displacements of the material points in the deformed state with the Lagrangian strain tensor, *η*_{ij}, defined by 2*η*_{ij}=(*u*_{i,j}+*u*_{j,i})+*u*_{k,i}*u*_{k,j}. Denote the stretches in the fundamental uniform state by *λ*_{i} subject to a constraint of incompressibility, *λ*_{1}*λ*_{2}*λ*_{3}=1. The material points on the free surface are given by *x*_{2}=0 with a semi-infinite slab of neo-Hookean material below. An arbitrary uniform stretch *λ*_{3} is allowed but, once imposed, it is held fixed. The non-uniform wrinkling deformations associated with bifurcation are restricted to satisfy plane strain conditions in the (*x*_{1},*x*_{2}) plane. The applied load parameter is the stretch *λ*≡*λ*_{1} associated with the uniform solution, i.e. the average stretch in the *x*_{1} direction, *λ*, is imposed with *λ*_{2}=*λ*_{3}/*λ*_{1}. Solutions which are periodic with wavelength, *l*=2*π*/*k*, with respect to the reference *x*_{1} coordinate are sought, having wavelength *λl* in the deformed state.

Let *ϕ*(** η**) be the strain energy per unit volume of the strained material with

*μ*as the ground state shear modulus. Let 2.1be the energy change/wavelength (per unit depth of undeformed material) relative to the imposed uniform state having strain

*η*^{(0)}associated with imposed

*λ*(and

*λ*

_{3}). Let (no sum on

*i*) be the displacements associated with the uniform solution and denote the total displacements by 2.2where the additional displacements

*U*

_{i},

*i*=1,2 are restricted to have the periodicity noted above with zero average stretch in the

*x*

_{1}direction. For a neo-Hookean material: The incompressibility condition is

*C*(

*λ*,

*U*)=0, where The modified energy functional, including a Lagrangian multiplier function,

*q*(

*x*

_{1},

*x*

_{2}), to enforce incompressibility, is 2.3To eliminate the terms linear in

*U*

_{i}, let

*q*=

*r*+

*Q*(

*x*

_{1},

*x*

_{2}) with

*r*≡

*λ*

_{2}/

*λ*

_{1}, where

*Q*(

*x*

_{1},

*x*

_{2}) has the same periodicity as

*U*

_{i}; the linear terms in

*U*

_{2}cancel. By periodicity, the term linear in

*U*

_{1,1}integrates to zero. The modified functional becomes: 2.4and 2.5At prescribed

*λ*, conditions of equilibrium and incompressibility are given by the requirement that the first variations of with respect to

*U*

_{i}and

*Q*vanish subject to periodicity and such that the overall stretch

*λ*is not altered by

*U*

_{i}.

The lowest order terms in the functional are quadratic in the unknowns (*U*,*Q*). The eigenvalue problem for the critical stretch, *λ*_{W}, is the variational problem, , based on the quadratic terms in equation (2.4) with (*U*,*Q*) and their variations restricted to have the periodicity noted earlier and to decay as . This is the Biot wrinkling problem which is briefly outlined below. The Euler equations for the problem are
2.6with boundary conditions on *x*_{2}=0:*U*_{1,2}+*rU*_{2,1}=0 and 2*rU*_{1,1}+*λ*_{1}*Q*=0. This problem admits separated periodic solutions of the form
2.7The two characteristic solutions to equation (2.6), (*f*(*kx*_{2}),*g*(*kx*_{2}),*h*(*kx*_{2}))=(*F*,*G*,*H*)*e*^{ksx2}, which decay to zero as have *s*=*r* and *s*=1, with *F*=−*G*,*H*=*Gr*^{−1}(*r*^{2}−1) for *s*=*r* and *F*=−*Gr*^{−1}, *H*=0 for *s*=1. Satisfaction of the boundary conditions on *x*_{2}=0 gives the eigenvalue condition
2.8and
2.9where the normalization *g*(0)=1 has been enforced. The solution holds for any wave number *k* and prescribed *λ*_{3} with *λ*_{W} given by equation (1.1).

Multiple eigenmodes with the same periodicity exist associated with the critical stretch. To define them, identify the *n*th mode in the set using the notation and *Q*=*q*^{(n)}. With *k*=2*π*/*l*,
2.10for *n*=1,2,3,…. Here, the period, *l*, which is the only length scale in the problem, is employed as a dimensional normalizing factor. Attention has been restricted to modes that are symmetric about *x*_{1}=0. The normalization in equation (2.10) is such that on *x*_{2}=0 the modal displacement normal to the free surface is .

In the analysis which follows, the total displacement will be expanded in the form
2.11with *q*^{0}(*λ*)=*r* and *ξ*_{n} as the dimensionless amplitude of the *n*th mode. Higher order terms in the expansion are denoted by (Δ*u*_{i},Δ*q*_{i}).

## 3. Koiter’s initial post-bifurcation analysis for wrinkling

### (a) The perfect system

Here, a general result for the energy change in the vicinity of the bifurcation point will be derived using a compact notation. A general relation is sought between the prescribed overall stretch, *λ*, and the amplitudes of the bifurcation modes in the equilibrium post-bifurcation state which are denoted collectively by *ξ*. In the general notation, this relation has the form
3.1For wrinkling, it turns out that *a*≠0 and it will not be necessary to obtain *b* because it will be seen that *a*≠0 already implies that the bifurcation solution is highly unstable. With reference to expansion (2.11), the bifurcation modes are represented collectively as
In the same compact notation, let *ξ*^{2}*U*^{(2)} denote all the terms that are quadratic in the amplitudes of the bifurcation modes. In this compact notation, the initial post-bifurcation expansion takes the form
3.2The modified energy functional equation (2.4) is denoted by . Equilibrium in the bifurcated state and the constraint on volume change require satisfaction of the variational equation
for all admissible *δU*≡(*δU*_{i},*δQ*) satisfying periodicity with no average stretch. Expand this condition about (*λ*_{W},0), noting that and , obtaining
with notation such as
In increasing powers of *ξ*, this becomes
3.3The eigenvalue problem governing bifurcation solved earlier is obtained by setting the terms of order *ξ* to zero for all admissible *δU*, i.e.
Each of the higher terms in the expansion, *U*^{(m)}, is admissible and, thus,
3.4The variational problem for *U*^{(2)} is obtained from the terms of order *ξ*^{2}:
(Orthogonality conditions on *U*^{(2)} relative to *U*^{(1)} must also be imposed if one solves for *U*^{(2)}, but this will not be necessary.) With *δU*=*U*^{(1)} in the above equation and use of equation (3.1), one obtains the compact equation for the bifurcation coefficient, *a*:
3.5For most problems, this condition gives *a*=0, but for the wrinkle problem with multiple modes we will find *a*≠0 and it is not necessary to proceed further. Detailed information on mode coupling will also emerge.

The above results can be used to express the energy change, , from the fundamental state. Because and , it also follows that
Using expansions for *λ* and *U*, making use of the above, one finds
Accounting for the terms that vanish by virtue of the eigenvalue problem, gives
3.6Equation (3.6) allows one to identify the expression for as the cubic terms in *U* from equations (2.4) and (2.5), i.e.
3.7with
Based on the quadratic terms in equation (2.4) and noting ∂*λ*_{2}/∂*λ*=−*r* and ∂*r*/∂*λ*=−2*r*/*λ*_{1}, it is also straightforward to obtain
3.8In equations (3.7) and (3.8), *λ*_{1} and *r* are evaluated at *λ*_{W}.

### (b) The imperfect system: lowest order contribution of a geometric imperfection

A slight imperfection in the form of a periodic undulation of the surface of the undeformed half-space is assumed:
3.9with as the amplitude of the imperfection in the *n*th mode.

The objective is to obtain asymptotic results for the effect of very small imperfections on behaviour in the vicinity of the bifurcation point and, in particular, on the occurrence of the wrinkling instability. Only the lowest order influence of the imperfections is sought following an approach similar to that of Koiter (Koiter 1945; van der Heijden 2009).

For the initial undulation, *Φ* in equation (2.1) can be written as
3.10For very small *δ*(*x*_{1}),
such that the lowest order contribution to owing to the imperfection is
3.11Then, note that
3.12where, to lowest order in , and is the Piola–Kirchhoff stress in the fundamental state. With , the lowest order contribution of the imperfection to is
3.13

## 4. The instability of wrinkling

### (a) Evaluation of the post-buckling coefficients

In the notation of equations (2.10) and (3.2), we consider the first *N* modes:
4.1A direct evaluation of the integral in equation (3.8) gives
4.2while equation (3.7) gives
4.3with
Coefficients for *N*=6 are listed but some results below have been computed with *N* as large as 10. The nonlinear coupling of modes 1 and 2 is illustrated by the cubic term, , in equation (4.3). It arises owing to the fact that products of quadratic terms from mode 1, proportional to , and linear terms from mode 2, proportional to , are not orthogonal. On the other hand, the phasing of cubic terms proportional to and is such that they integrate to zero.

The contribution (3.13) from the initial imperfection is
4.4with *c*=0.95467.

### (b) The post-bifurcation equations

The modified energy functional in equation (3.6) plus the contribution owing to the imperfection is (illustrated for *N*=6)
4.5This result holds for any prescribed value of *λ*_{3}, with *λ*_{W} given by equation (1.1). Equilibrium requires for *i*=1,…,*N*. For *N*=6, the equations are:
4.6

These relations are asymptotically valid in the vicinity of the bifurcation point for sufficiently small imperfections.

### (c) Post-bifurcation solutions: the perfect system

#### (i) The two-mode approximation

Setting all the mode amplitudes to zero in equation (4.6) except for the first two modes and their imperfections, one has
4.7For the perfect system (), the solutions of interest are
4.8Subsequently, it will be evident why the solutions of interest are those associated with *overall compressive strains less than the bifurcation value* (i.e. *λ*−*λ*_{W}≥0 or, equivalently, *ε*−*ε*_{W}≤0 with *ε*=1−*λ* as the compressive strain). The solutions are shown in figure 2. The existence of the non-zero cubic term, , implies that *wrinkling bifurcation is unstable* at *λ*=*λ*_{W} because the energy change on the equilibrium path relative to the bifurcation state,
is negative. The shape of surface wrinkle for the combined two-mode approximation,
is plotted in figure 3. The wrinkle displays a deep-pronged penetration of the free surface into the material with relatively flat broad crests on either side.

#### (ii) N-mode approximations

By including the third mode in equation (4.6), one sees that the two-mode solution is indeed only an approximation—the term *b*_{123}*ξ*_{1}*ξ*_{2} in the third equation requires non-zero *ξ*_{3}. A finite number of modes can only generate an approximate solution to the order considered because the modes in the infinite set are all coupled through the quadratic terms in the equilibrium equations.

For a given *N*, the solution to the system (4.6) has the form *ξ*_{n}=*α*_{n}(*λ*−*λ*_{W}) for *n*=1,…,*N* and, thus, the normalized surface undulation, *u*_{2}(*x*_{1},0)/[(*λ*−*λ*_{W})*l*] in figure 3 is independent of (*λ*−*λ*_{W}) to the order considered in this paper. While any shape is possible according to the bifurcation solution, the post-bifurcation analysis identifies a definite shape, assuming periodicity. The normalized shape of the wrinkle at the surface in figure 3 has been determined by a sequence of calculations, each with *N* modes, for *N* ranging from 2 to 10. A standard numerical iterative algorithm for solving systems of coupled nonlinear algebraic equations was used to generate the *α*_{n}. The solution for the system of *N*−1 equations was employed as the initial assumption in the iteration for the solution for the system with *N* equations, thereby leading to the regular progression of shape approximations shown.^{1} The sequence is trending towards an incipient crease-like shape as more and more terms in the approximation are considered, although the shape for *N*=10 does not yet appear to have converged.

### (d) Post-bifurcation solutions: the imperfect system

#### (i) The two-mode approximation

Explicit results for the reduction of the compressive strain at which wrinkling becomes unstable are now given for an imperfection in the first mode , . The second part of equation (4.7) gives . Substituting this into the first part of equation (4.7), one finds , which provides the relation between *λ* and *ξ*_{1}. Denote the minimum of *λ* (i.e. the maximum compressive strain in the presence of the imperfection) by *λ**; it is associated with d*λ*/d*ξ*_{1}=0 (figure 2). Solving for *λ**−*λ*_{W}, one finds
4.9The stretch *λ** for the imperfect system corresponds to the maximum compressive overall strain on the equilibrium path. The solution at *λ** is unstable in the sense that it would snap dynamically and undergo a finite deformation into another configuration—almost certainly to a fully developed crease as will be seen later.

The effect of very small imperfections in lowering the compressive strain at which wrinkling becomes unstable is dramatic owing to the fact that it is proportional to , as seen in figures 2 and 4. The type of nonlinear coupling among simultaneous modes in wrinkling is rare but it is similar to that in two structural problems that also have multiple buckling modes and are notoriously imperfection-sensitive—the elastic buckling of cylindrical shells under axial compression (Koiter 1945; van der Heijden 2009), and spherical shells under external pressure (Hutchinson 1967).

#### (ii) N-mode approximations

Consider again the half-space with an initial imperfection in the first mode, , with for *n*>1. As in the perfect case, a sequence of calculations has been made with an increasing number of modes in the approximation. For any *N*, the solution to equation (4.6) at the point of the maximum overall compressive strain has the form and . The coefficient *c** is given in table 1, and the results for the reduction in the compressive stretch, *λ**−*λ*_{W}, at which wrinkling becomes unstable is plotted as a function of the imperfection amplitude in figure 4 for *N* ranging from 2 to 6. The results appear to be converging to a curve slightly above that for *N*=6. Equation (4.9) based on the two-mode approximation underestimates the reduction in the compressive strain at the wrinkling instability by about 35 per cent.

With *ξ*_{n}=*α*_{n}(*λ*−*λ*_{W}) for *n*=1,…,*N* as the *N*-mode solution for the perfect system in equation (4.6), consider an *imperfection in the shape of the solution for the perfect system*, i.e. for *n*=1,…,*N*, with as the single imperfection amplitude. With *ξ*_{n}=*α*_{n}*ξ* for *n*=1,…,*N*, it is straightforward to show that each of *N* equilibrium equations (4.6) reduces to the same equation: . The maximum compressive strain that can be imposed prior to instability is given by
4.10The imperfection-sensitivity implied by this result is similar to that predicted for an imperfection in the shape of the first mode.

An alternative instability condition for the imperfect half-space will be discussed in connection with the numerical solutions presented in §5.

## 5. Plane strain finite-element simulations of a half-space with initial imperfections

Two types of initial surface imperfections have been considered in the numerical simulations: a sinusoidal imperfection, , and a periodic array of non-interacting initial exponential depressions of the surface specified by
5.1A finite-element mesh conforming to the initial surface undulation was created on a rectangular region in the (*x*_{1},*x*_{2}) plane of width *L* and depth *D*=10*L* for the sinusoidal imperfection. The surface is traction-free, whereas the shear traction and *u*_{2} are taken to be zero on the bottom surface. The depth is sufficient to ensure that the boundary conditions on the bottom have no influence on the wrinkling behaviour. For the sinusoidal imperfection, the periodic boundary condition is imposed on the vertical sides of the region and *L* is taken to be *l*. For the exponential imperfection (5.1), periodic boundary conditions on the sides are also assumed for computational convenience. The imperfection is located within the centre of the region and *L* is chosen to be 20*l* with *D*=10*l* so as to ensure that there is essentially no interaction between neighbouring imperfections or the bottom—the results for this case can be regarded as that of an isolated imperfection of the form (5.1).

Plane strain (*λ*_{3}=1) finite-element simulations are performed via the commercial software, ABAQUS (2008). Considering that the instability and the wrinkling–creasing transition occur at the central region of the upper surface, a very fine mesh is used in this region with a ratio of *l* to the element size taken to be approximately 2000. In the finite-element simulations, the incompressible neo-Hookean material model is employed (ABAQUS 2008). The hybrid element (CPE6MH) suitable for simulations of incompressible materials is adopted. To introduce the initial surface imperfections, finite-element simulations are first run by specifying the boundary conditions on the upper surface for the sinusoidal imperfection and for the exponential imperfection. The function *IMPERFECTION in ABAQUS (2008) converts the displacements from this step to an initial stress-free geometric imperfection. This procedure is equivalent to directly meshing a block of stress-free material with the initial surface undulation. Simulations are performed to track the occurrence of the local instability and the formation of a crease, as reported below. Self-contact interaction is defined for the upper surface of the block. When a local wrinkling instability occurs, the global matrix of the system may be singular and the Riks method, the commonly used numerical method for dealing with limit points, will fail. In the present simulations, the pseudo-dynamic method has been adopted. A brief description of the key idea behind this nonlinear solution method is outlined as follows.

The nonlinear equations solved in a finite-element analysis can be written as
5.2where ** Ξ** is a nonlinear function of

**, symbolizing the displacements of the nodes and**

*u***(**

*Ξ***) and**

*u***denote internal forces and applied loads at the nodes, respectively. The pseudo-dynamic method regularizes the unstable problem by adding volume-proportional damping to the model such that equation (5.2) becomes 5.3where 5.4Here,**

*Ω***is an artificial mass matrix calculated with unity density,**

*M**c*is a damping factor,

**=Δ**

*v***/Δ**

*u**t*is the vector of nodal velocities and Δ

*t*is the increment of time. When the model is stable (quasi-static), viscous forces and viscous energy dissipation are very small such that the artificial damping does not perturb the solution significantly. When the structure goes dynamically unstable, however, nodal velocities increase and, consequently, part of the strain energy released is dissipated by the damping. In simulations of the wrinkling problem, pseudo-dynamic regularization, which is now a standard feature in ABAQUS (2008), allows solutions to be generated under prescribed overall compressive strain when the wrinkle becomes unstable and develops into a crease. The role played by the damping factor

*c*here is similar to that of the regularization factor in the Tikhonov regularization method (Tikhonov & Arsenin 1977), which is widely used to deal with ill-posed inverse problems.

### (a) Wrinkling instability for sinusoidal imperfections

Figure 5 presents the overall compressive strain at the point of wrinkling instability, *ε**=1−*λ**, as a function of the imperfection amplitude, . Included in this figure is the local compressive strain parallel to the surface at the deepest point of the wrinkle, *ε*_{A}=1−*λ*_{A}. As will be discussed in more detail below, the imperfect half-space becomes unstable with the local strain at the deepest point attains the Biot wrinkling strain, i.e. *ε*_{A}=*ε*_{W}. The simulations presented in figure 5 again reveal the extraordinarily strong imperfection-sensitivity of the overall compressive strain at instability, *ε**. Moreover, the numerical results confirm the accuracy of asymptotic result, , which has been included in figure 5, for reductions in overall strains larger than 0.1. In particular, imperfections larger than about reduce the overall strain at instability to below the level needed to sustain a crease in the perfect system, i.e. *ε**<*ε*_{C}. This is a tiny initial undulation with an amplitude less than 1 per cent of its wavelength. The crease that forms at this low value of *ε** is localized within the region of strain concentration at the bottom of the larger scale wrinkle.

Further details of the evolution of the instability are illustrated in figures 6 and 7 for the case . The difference in height between the highest and deepest points on the surface, Δ*u*_{2}, increases monotonically as the instability develops. Figure 6 displays the local compressive strain at the deepest point, *ε*_{A}, and the overall strain, *ε*, as functions of Δ*u*_{2}/*l* up to the onset of instability. At the onset, *ε*_{A}=*ε*_{W}, as already noted, and the overall compressive strain attains its maximum, *ε**. Upon attaining the onset condition at A, a small-scale wrinkle forms within a narrow region on either side of A (figure 7). This small-scale wrinkle evolves into a fully developed crease under conditions in which the overall strain remains essentially unchanged at *ε**. The wavelength of the small-scale wrinkle is comparable with the size of the finite elements, but once the crease develops, the crease depth is large when compared with element size, as seen in figure 7*b*. The crease relaxes the compressive strain in its vicinity as seen from the plot of the compressive strain at the surface as a function of horizontal distance from the crease in the deformed body in figure 8. For horizontal distances from the centre-line of the crease less than about 2.3*l*_{C}, the surface strain has been reduced below the creasing strain, *ε*_{C}=0.35, and at greater distances the strain does not exceed the wrinkling strain *ε*_{W}=0.456. The numerical simulations indicate that the characteristic point on the surface at which the strain attains *ε*_{C} is nearly coincident with the inflection point corresponding to the transition of the surface from being convex to concave. No attempt has been made to simulate behaviour at overall compressive strains beyond *ε** which would drive the crease even deeper than that shown in figures 7 and 8 and possibly nucleate new wrinkles and creases. Such calculations have recently been performed and compared with experimental observations by Cai *et al.* (submitted).

In connection with figure 6, it was noted that the onset of wrinkling instability is associated with (almost) simultaneous satisfaction of two conditions: (i) attainment of a maximum in the overall compressive strain, and (ii) *ε*_{A}=*ε*_{W}. The analytical modelling of wrinkling instability in §4*d* is based on condition (i). Motivated by the numerical findings related to the role of condition (ii), the analytical approach in §4*d* has been used to compute the overall stretch *λ** at which the local wrinkling condition, *ε*_{A}=*ε*_{W}, is met at the deepest point of the wrinkle on the surface at *x*_{1}=*l*/2. The details of this calculation will not be given because they involve only a minor extension of the analysis in §3. The result has precisely the same functional form, , as in the case of condition (i), where, as before, the coefficient *c** depends on the number of modes, *N*, in the approximation. The coefficient is presented in table 1 along with that computed earlier based on condition (i). According to the analytical approximation, the local wrinkling condition (ii) is attained slightly before the maximum overall strain is reached. However, the difference between the overall critical stretch *λ** from two conditions is nearly negligible when six modes are used in the calculation (table 1). Thus, both the analytical and numerical predictions indicate that the two conditions, (i) and (ii), are attained nearly simultaneously at the onset of wrinkling instability and crease formation.

### (b) Wrinkle instability for exponential imperfections

Simulations with the exponential imperfection (5.1) have also been carried out with results presented in figure 9. Wrinkling instability and the formation of a crease again occurs when the local strain at the deepest point of the surface wrinkle attains *ε*_{W}. A slight local depression on the surface of the half-space reduces the overall strain at the wrinkling instability to levels similar to that seen for the sinusoidal imperfection, based on comparable values of the normalized imperfection amplitudes that have been defined.

## 6. Conclusions

The post-bifurcation analysis of Biot’s wrinkling problem reveals that wrinkling is extremely unstable and highly imperfection-sensitive. Wrinkling is also seen to be one pathway to the finite amplitude creasing mode. Wrinkling is so unstable and imperfection-sensitive that well-developed wrinkles are not likely to be observed—a slight wrinkle will become dynamically unstable and trigger the formation of a crease. In this sense, the crease can be regarded as the collapse state of a wrinkle. The wrinkle/crease connection has an analogue in the elastic buckling of cylindrical shells under axial compression and spherical shells under external pressure which, like wrinkling, are characterized by multiple bifurcation modes associated with the critical stress. Buckling of these shell structures is also so unstable and imperfection-sensitive (Hutchinson 1967; van der Heijden 2009) that their short wavelength bifurcation modes are almost never observed because they become unstable at very small amplitudes and snap dynamically into a collapse state. Buckling modes observed in the collapsed state of the shell typically have much larger wavelengths than those of the bifurcation modes. A few experiments have employed high-speed cameras to capture the bifurcation modes right after they are triggered (Brush & Almroth 1975) or have used a mandrel to arrest the buckles immediately after they have formed (Carlson *et al.* 1967).

In addition, these shell structure/loading combinations are so imperfection-sensitive that, of the large number of shells tested over many years, none has reached a buckling load greater than about one-half of the buckling load of the perfect shell when the radius to thickness ratio exceeds 1000. In this respect, as well, there may be a close analogue to wrinkling/creasing, i.e. the imperfection-sensitivity of wrinkling may be so strong that the maximum compressive strain of any actual realization of an elastomer layer will always lie well below Biot’s wrinkling strain (1.1) owing to inevitable surface imperfections.

There is one important respect in which wrinkling of a uniform half-space differs from cylindrical and spherical shell buckling—the wrinkling wavelength is undetermined and can be arbitrarily small. In principle, a perfectly flat surface should reach the Biot wrinkling strain but in practice, as noted above, it seems reasonable to assume that imperfections will always be present at some scale to trigger creases at strains just above the creasing strain. As Hohlfeld & Mahadevan (2011) have noted, wrinkles and creases are always lurking to destabilize a smooth surface of a compressed elastomer because their scale can be arbitrarily small. Examples have been presented in this paper for both sinusoidal and isolated imperfections, wherein small-scale wrinkles form in the vicinity of the point on the surface of maximum compression—wrinkles within wrinkles—and then spontaneously collapse into a local crease. An open question concerns the lower limit on the size of these instabilities. A continuum representation such as the neo-Hookean material, with the elastomer represented by a constitutive model having no material length dependence, provides no lower limit on the scale of the instabilities. Surface effects such as a stiff thin layer of oxidized material would provide a limit. Strain gradient strengthening associated with deformation gradients that become comparable with scale of the polymeric microstructure would also place a lower limit on the size of the instabilities, but such effects have not yet been quantified for elastomers.

Finally, we note that the unstable wrinkling behaviour of the uniform neo-Hookean half-space is in sharp contrast to the highly stable wrinkling behaviour of a system comprising a thin stiff film bonded to compliant half-space substrate. The film–substrate system buckles into wrinkling modes at very small compressive strains (Allen 1969). These systems can be compressed well beyond the critical bifurcation strain with the buckled state remaining stable (Cai *et al.* 2011). It is not at all unusual for experimental systems to display stable wrinkling behaviour at an overall compressive strain 10 times the bifurcation strain. Imperfections play a secondary role in the behaviour of these systems. It remains for future work to explore the full parameter space of film–substrate systems to uncover the range of parameters in which a transition occurs from the highly stable buckling behaviour associated with very stiff films to the highly unstable behaviour associated with wrinkling of the uniform half-space.

## Acknowledgements

The authors acknowledge the input from a reviewer who suggested the form of the incompressibility condition used in this paper. An earlier version of the paper employed an equivalent but less compact expression. Y.P.C. acknowledges the financial support from Tsinghua University (2009THZ02122).

## Footnotes

- Received June 16, 2011.
- Accepted July 27, 2011.

- This journal is © 2011 The Royal Society