The secondary buckling transition: wrinkling of buckled spherical shells

# The secondary buckling transition: wrinkling of buckled spherical shells

Sebastian Knoche Department of Physics, Technische Universität Dortmund, 44221 Dortmund, Germany    Jan Kierfeld Department of Physics, Technische Universität Dortmund, 44221 Dortmund, Germany
###### Abstract

We theoretically explain the complete sequence of shapes of deflated spherical shells. Decreasing the volume, the shell remains spherical initially, then undergoes the classical buckling instability, where an axisymmetric dimple appears, and, finally, loses its axisymmetry by wrinkles developing in the vicinity of the dimple edge in a secondary buckling transition. We describe the first axisymmetric buckling transition by numerical integration of the complete set of shape equations and an approximate analytic model due to Pogorelov. In the buckled shape, both approaches exhibit a locally compressive hoop stress in a region where experiments and simulations show the development of polygonal wrinkles, along the dimple edge. In a simplified model based on the stability equations of shallow shells, a critical value for the compressive hoop stress is derived, for which the compressed circumferential fibres will buckle out of their circular shape in order to release the compression. By applying this wrinkling criterion to the solutions of the axisymmetric models, we can calculate the critical volume for the secondary buckling transition. Using the Pogorelov approach, we also obtain an analytical expression for the critical volume at the secondary buckling transition: The critical volume difference scales linearly with the bending stiffness, whereas the critical volume reduction at the classical axisymmetric buckling transition scales with the square root of the bending stiffness. These results are confirmed by another stability analysis in the framework of Donnel, Mushtari and Vlasov (DMV) shell theory, and by numerical simulations available in the literature.

## I Introduction

When spherical shells, such as sports and toy balls or microcapsules, are deflated, they always go through the same sequence of shapes, see fig. 1. At small deflation, the capsule remains spherical. Upon reducing the volume, an axisymmetric dimple forms in an abrupt transition. For sufficiently thin shells, this dimple finally loses its axisymmetry upon further volume reduction, resulting in a polygonally buckled shape. This is shown by daily life experience, microcapsule experiments Quilliet2008 (); Datta2010 (); Datta2012 (), and computer simulations Quilliet2008 (); Vliegenthart2011 (); Quilliet2012 () based on triangulated surfaces (e.g. with the surface evolver) or finite element methods Vaziri2008 (); Vaziri2009 (); Vella2011 (). This generic behaviour also applies to other types of deformation, for example when a dimple is formed by indenting the capsule with a point force Vaziri2008 (); Vaziri2009 (); Pauchard1998 (), when the capsule is pressed between rigid plates Vaziri2009 () or when the capsule adheres to a substrate Komura2005 ().

The first buckling transition, where an axisymmetric dimple forms, is a classical problem in continuum mechanics, which is very well understood. Linear shell theory can be used to calculate the onset of instability of the spherical shape, see references Landau1986 (); Ventsel2001 (). Furthermore, nonlinear shell theory has been used to investigate the postbuckling behaviour Koiter1969 (), which revealed that the buckled shape is unstable with respect to further volume reduction if the pressure is controlled. Numerical analyses of nonlinear shape equations have been used to obtain bifurcation diagrams for the axisymmetric deformation behaviour Bauer1970 (); Knoche2011 (). In this paper, we will approach this classical axisymmetric buckling using two different models for axisymmetric deformations. The first model is based on nonlinear shell theory, which leads to a complete set of shape equations which are to be solved numerically Knoche2011 (). The second model is an approximate analytical model which has been proposed by Pogorelov Pogorelov1988 ().

In contrast, the secondary buckling transition, where the dimple loses its axisymmetry, has been merely observed experimentally or in computer experiments so far. A theoretical approach which explains the mechanism underlying this instability and which can predict the corresponding critical volume for the secondary buckling transition is still lacking. Here we will offer an explanation within continuum elasticity theory. This also demonstrates that polygonal capsule shapes can also occur in the absence of any discretisation effects. For crystalline elastic capsules, defects in the triangulation gives rise to additional faceting effects upon deflation Yong2013 ().

Analysing the results of the axisymmetric models, we get a hint on the physical mechanism of the secondary buckling transition: we observe a region of compressive hoop stress, which is located in the inner neighbourhood of the dimple edge, just in the place where the secondary buckling occurs in experiments and simulations. In order to release the compressive stress, the circumferential fibres buckle out of their circular shape if the hoop stress reaches a critical value; quite comparable to the Euler buckling of straight rods Landau1986 (). A quantitative investigation of the secondary buckling transition therefore consists of two steps: Firstly, quantifying the stress distribution in the axisymmetric buckled configuration, and secondly, finding the critical compressive stress at which the axisymmetric configuration loses its stability.

The first task is readily accomplished by the two axisymmetric models, the approach by nonlinear shape equations or Pogorelov’s approximate analytical model. For the second task, we will develop a simplified model which captures the essential features of the geometry and stress distribution of the axisymmetric buckled shape, and derive a critical compressive stress by using the stability equations of shallow shells Ventsel2001 ().

We verify our results by a second approach in which we conduct a stability analysis of the full shape, where the exact geometry and shape distribution is taken into account. The framework to do this is the nonlinear DMV (Donnel, Mushtari and Vlasov) theory of shells Niordson1985 (); Ventsel2001 ().

The resulting critical volumes for the primary and secondary buckling transition are presented in a phase or stability diagram in fig. 12, from which the state of a deflated capsule with given bending stiffness and volume difference can be read off and which is the main result of the paper. A short account of these results focusing on the secondary buckling mechanism and the parameter dependencies of the critical buckling volumes has been given elsewhere Knoche2014 ().

## Ii Axisymmetric primary buckling of shells

### ii.1 Nonlinear shape equations

Figure 2 shows the parametrisation of the capsule midsurface. The undeformed shape is parametrised in cylindrical polar coordinates by the functions and , where is the arc length. The slope angle , defined by the two relations

 dr0ds0=cosψ0anddz0ds0=sinψ0, (1)

see fig. 2 b), permits simple calculation of the principal curvatures of the midsurface. These are found in meridional and circumferential direction, respectively, and read

 κs0=dψ0ds0andκφ0=sinψ0r0. (2)

In our case, the spherical reference configuration is explicitly given by and , and the principal curvatures reduce to . The geometrical relations (1) and (2) also hold for the deformed midsurface when all indices “0” are omitted.

Upon deformation into a different axisymmetric shape, the midsurface undergoes stretching and bending. We measure the stretches in meridional and circumferential direction by

 λs=ds/ds0andλφ=r/r0, (3)

respectively. The function defined in this context determines the position at which a shell element originally located at can be found after deformation. The strains that correspond to these stretches are centred around and are defined as and . To measure the change in curvature, the bending strains

 Ks=λsκs−κs0andKφ=λφκφ−κφ0 (4)

are defined Libai1998 (); Pozrikidis2003 ().

The deformation results in an elastic energy which is stored in the membrane. This elastic energy can be calculated as the surface integral over an elastic energy density (measured per undeformed surface area), which we assume to be of the simple Hookean form

 wS=12EH01−ν2(e2s+2νeseφ+e2φ)+12EB(K2s+2νKsKφ+K2φ). (5)

In this expression, is the (three-dimensional) Young modulus, the membrane thickness, is the (three-dimensional) Poisson ratio which is confined to , and is the bending stiffness. It can be shown Libai1998 () by the principle of virtual work that the meridional tension and bending moment can be derived as

 τs =1λφ∂wS∂es=EH01−ν21λφ(es+νeφ), (6) ms =1λφ∂wS∂Ks=EB1λφ(Ks+νKφ).

from the energy density. The corresponding relations for the circumferential tension and bending moment are obtained by interchanging all indices and in these equations. Note that the tensions and bending moments are measured per unit length of the deformed midsurface, whereas the energy density is measured per unit area of the undeformed midsurface – this is the reason why the factors occur in (6).

To close the problem of determining the deformed shape, we need equations of equilibrium. These read, for the tangential force, normal force and bending moment,

 0 =−cosψrτφ+1rd(rτs)ds−κsq, (7) 0 =−p+κφτφ+κsτs+1rd(rq)ds, (8) 0 =cosψrmφ−1rd(rms)ds−q. (9)

In these equations, is the transversal shear force, and is the applied normal pressure, which can also be interpreted as a Lagrange multiplier to control the capsule volume.

For a numerical treatment, the equilibrium equations should be written as a system of first-order differential equations, which are called shape equations. They follow from (1), (2) and (7) - (9), see Knoche2011 (). Boundary conditions must be imposed which assure that the capsule is closed and has no kinks at the poles. In ref. Knoche2011 (), a detailed discussion of this issue is given, as well as numerical procedures that are suitable for solving the shape equations.

For the further analysis, it is convenient to introduce dimensionless quantities by measuring tensions in units of and lengths in units of . Specifically, this results in a dimensionless bending stiffness which is the inverse of the Föppl-von-Kármán-number ,

 ~EB≡EBEH0R20andγFvK=1~EB. (10)

### ii.2 Analytic Pogorelov model

The second approach to axisymmetric buckled shapes is based on a model of Pogorelov Pogorelov1988 (). The basic idea is that for small bending stiffness, the shape of the axisymmetric dimple will be close to an isometric deformation of the sphere: a shape where a spherical cap is mirror inverted (see fig. 3, grey lines). For vanishing bending stiffness, , this shape has vanishing elastic energy and, thus, represents the stable equilibrium shape, since it is free of stretching. It only involves bending to invert the curvature of the cap. However, the bending strain at the edge of the inverted cap is infinitely large, which gives rise to infinitely large bending energy for . Thus, switching from to , the sharp edge of the dimple has to be smoothed out.

In order to describe the deformation from the isometric shape to the final smooth shape, we follow the ideas of Pogorelov and introduce displacements and in - and -direction, respectively, see fig. 3. Assuming that and are small (and some further simplifications), we use linear shell theory to calculate the bending and stretching energies in the final shape. This technique is quite remarkable because it enables us to describe large deformations with linear shell theory by choosing not the undeformed shape as reference state, but the isometric buckled shape. In the following, we will summarise this procedure and results from Pogorelov1988 () regarding the axisymmetric primary buckling because the results will provide the basis for a quantitative theory of the secondary buckling transition.

First, there are geometric relations for isometric deformations of spheres, which are obtained by mirror reflection of a spherical cap, see fig. 3 a). The resulting dimple is characterised by its opening angle , from which we can calculate the dimple radius , depth and volume difference . For later use, we already introduce a first order approximation in , since we will assume that the dimple is small compared to . The exact and approximated relations then read

 rD =R0sinα≈αR0 (11) h =R0−R0cosα≈α2R0/2 (12) ΔV =2h2π(3R0−h)/3≈πα4R30/2. (13)

In order to evaluate the elastic energies in the final shape, we split it into different regions which are investigated separately. We define the inner neighbourhood and outer neighbourhood of the dimple edge (see fig. 3) as the regions in which the displacements are significant. This is only the case in a narrow strip to both sides of the dimple edge, see (20) below. Outside these regions, the displacements are negligible; these regions are referred to as (for the rest of the dimple) and (for the rest of the undeformed part).

The elastic energies that we have to evaluate consist of the bending needed to invert the curvature in , and the bending and stretching due to the displacements in and . The uniform compressive strains in and , which result from the negative inner pressure and which are already present in the spherical (pre-buckled) shape, are neglected.

The bending energy in is readily written down. From the spherical shape to the isometric deformed shape, the curvatures change from to , thus giving bending strains , where the stretches in the definition (4) are neglected. The energy density (5) is integrated over the area of the inverted spherical cap. We neglect that the area of the dimple should be reduced by the area of since we assume that is much smaller than , and hence the area of is . Thus, the bending energy in reads

 UB(G2) =12EB∫G2dA{K2s+2νKsKφ+K2φ} =4πα2EB(1+ν). (14)

In the region , the deformation energy is governed by the displacements and (where is the undeformed arc length). From “graphic considerations” Pogorelov1988 () Pogorelov concludes that a meridian will not stretch or compress very much: . Using our nonlinear shape equations, we checked that this assumption holds at the secondary buckling transition if the Poisson ratio is not too large; in fact, the magnitudes of and (in a root mean square sense) are found to behave roughly as at the secondary buckling transition. At the first buckling transition, we find that depends on , and that our assumption holds for small reduced bending rigidities. The stretching of circumferential fibres (cf. eq. (3)) results in the strain , where the approximation holds when is sufficiently narrow. The integration of the stretching energy is performed over the area element with . Here the arc length coordinate was centred around the dimple edge and runs up to a point where the displacements have decayed sufficiently to be neglected. The stretching energy of the outer neighbourhood is thus given by

 US(Go) =12EH01−ν2∫Goe2φdA =πEH0(1−ν2)rD∫ε0u(s0)2ds0 (15)

Pogorelov approximates the bending strains as and . In the integral of the bending energy, a term occurs, which can be readily integrated to , and we need to specify boundary conditions for the displacement to proceed. We require the dimple to have a horizontal tangent at the dimple edge , hence . At the other end , the displacement shall have decayed and we enforce . The resulting expression for the bending energy is

 UB(Go)=πEBrD∫ε0v′′(s0)2ds0−παEBνrD/R0. (16)

In total, the elastic energy of the outer neighbourhood is therefore given by

 U(Go)= ∫ε0ds0{πEBrDv′′(s0)2+πEH0(1−ν2)rDu(s0)2} −παEBνrD/R0. (17)

Analogously, the elastic energy of the inner neighbourhood can be calculated. The stretching energy is the same as for the outer neighbourhood. But note that the sign of is negative in this case because we have hoop compression in the inner neighbourhood and hoop stretching in the outer neighbourhood (we already anticipate that the result for will be positive in the outer neighbourhood, see fig. 3). The bending strains have to be modified because we have to take into account that the curvature of the inverted cap is already inverted, and we get and . The resulting elastic energy of the inner neighbourhood is given by

 U(Gi)= ∫0−εds0{πEBrDv′′(s0)2+πEH0(1−ν2)rDu(s0)2} −4παEB(1+ν)rD/R0+παEBνrD/R0. (18)

The integrand coincides with the result (17) for the outer neighbourhood; only the constant terms differ.

To find the functions and which represent the final shape, we have to minimise the functional of total elastic energy with respect to and . During the variation and we keep the parameters , , and of the isometric shape constant. Since the volume change due to and can be neglected (compared to the large due to the isometric deformation), this corresponds to a variation at constant capsule volume, and we do not need to introduce a Lagrange multiplier controlling the volume. As the integrals in the elastic energies of the inner and outer neighbourhood are identical, we expect a symmetric shape of the dimple, with an odd function and an even function . It is thus sufficient to determine the solution on the interval by minimising (17).

During the minimisation we have to impose a constraint on and , because the energy functional was set up under the assumption of vanishing (or negligible) meridional strain . The final solution must satisfy this condition, which can be written as

 u′(s0)+αv′(s0)+12v′(s0)2=0. (19)

Furthermore, the variation has to respect the boundary conditions, which are evident from geometrical considerations: so that the capsule is not ripped apart at the dimple edge, for a horizontal tangent at the dimple edge, and because the displacements must have decayed at .

The number of parameters in the problem can be greatly reduced with a suitable nondimensionalisation by introducing characteristic length and energy scales. Inspection of the integrand in (17) and the constraint (19) shows that the substitutions

 s0 =ξ¯s0,ξ≡[~EB(1−ν2)]1/4(R0rDα)1/2 (20) u =ξα2¯u, (21) dv/ds0 =α¯w, (22)

with a typical arc length scale prove useful. For small , the length scale is also small, which proves that the regions and are indeed narrow strips. The substitutions lead to a dimensionless form of the energy (17),

 U(Go) =Uξ∫¯ε0d¯s0{¯w′(¯s0)2+¯u(¯s0)2}+const, (23) Uξ ≡πα5/2r1/2DR3/20EH0(1−ν2)1/4~E3/4B, (24)

with an energy scale . Using the geometric relations (11) - (13), the scaling parameters , , and can be expressed as functions of the elastic moduli, the reduced volume difference , and the capsule radius ,

 α ≈(83ΔVV0)1/4 (25) ξ ≈[~EB(1−ν2)]1/4R0 (26) Uξ ≈π(83)3/4EH0(1−ν2)1/4(~EBΔVV0)3/4R20. (27)

where we used the first-order approximations in .

Both the arc length scale and the energy scale emerge from the competition of stretching and bending energies in (17) (under the constraint (19)): gives the typical arc length size of the neighbourhoods and and gives the typical energy of the buckled configuration. The final result for the Pogorelov buckling energy, which is obtained after minimisation of the total energy with respect to and will differ from only by a numerical prefactor, see eq. (34) below.

For the minimisation of the total energy, the integral term in eq. (23) has to be minimised. Note that the limit of the integral has been rescaled, too, according to (20). Following Pogorelov, we consider the case of small , where because according to (26). Thus our task for the calculus of variations is to minimise the functional

 J[¯u,¯w]=∫∞0d¯s0{¯w′2+¯u2} (28)

subjected to the constraint (19) which reads

 ¯u′+¯w+12¯w2=0 (29)

in rescaled variables and with boundary conditions

 ¯u(0)=0,¯w(0)=−1,¯u(∞)=0,¯w(∞) =0. (30)

Pogorelov solved the constrained variational problem analytically. His results for the minimising functions and are presented in appendix A. These functions are defined piecewise, due to some simplifications, on two intervals and , where the optimal choice for is . The minimal value of the functional is found to be .

Now we can switch back from the nondimensionalised quantities to physical quantities in order to analyse the features of Pogorelov’s model, plot solutions and compare them to our results from the nonlinear shape equations. The rescaling of the functions and describing the capsule shape is obviously given by (20) - (21). We also have to take into account that the origin of was shifted to the dimple edge. In the coordinate system of the nonlinear shape equations, the origin of starts at the south pole of the capsule, and the dimple edge (of the isometric deformation) is located at . When comparing these solutions, the functions from the Pogorelov model must be shifted.

The displacements and can be used to plot the deformed capsule shape and to calculate further properties, like curvatures and tensions. With the strain definitions and , we obtain tensions from the linearised Hookean law, see (6),

 τφ=EH01−ν2u(s0)rDandτs=ντφ. (31)

The definitions of the bending strains imply curvatures

 κs ={−1/R0+v′′(s0),s0

The total elastic energy is obtained by adding . We see that the constant terms cancel each other, and only the integral terms of and survive. Each integral term gives . Thus, our final result for the elastic energy of an axisymmetric buckled capsule with a given volume difference is

 UPog =2JminUξ (34)

### ii.3 Comparison of the two models

Figure 4 shows a plot of the capsule shape and tension distribution in a solution of the shape equations and a solution of the analytic Pogorelov model. The shape is very well captured by the analytic model. Only in the close-up, deviations from the solution of the shape equations can be recognised. The behaviour of the hoop tension also agrees well in both models. We see a characteristic region of strong compression, which is located in the inner neighbourhood of the dimple edge, in which the models even agree quantitatively. In this region, the compressive hoop tension has a parabolic profile to a good approximation. In the outer neighbourhood of the dimple edge, a corresponding region of strong hoop stretch is present. The meridional tension is quite badly captured by the Pogorelov model, which is a consequence of the strong simplification of vanishing meridional strain .

A simple explanation why hoop compression arises in the inner neighbourhood and hoop tension in the outer neighbourhood of the dimple edge is evident from the close-up in fig. 4 a): after smoothing of the dimple edge, the inner neighbourhood is displaced horizontally to the left, and the outer neighbourhood to the right. Thus, the circumferential fibres of the inner neighbourhood are compressed and those of the outer neighbourhood are stretched.

Now, let us turn to the comparison of the elastic energies of the deformed shapes, from which the bifurcation behaviour can be deduced (see Knoche2011 () for bifurcation diagrams of spherical shells). First of all, there is the trivial (spherical) solution branch, which can be handled analytically. If the deformed capsule is a sphere with radius , the strains are uniform and homogeneous, . The bending strains vanish. The elastic energy density (5) must be integrated over the undeformed surface, which has an area . Hence, the elastic energy for a spherical shape is given by

 Usph=4πEH01−νR20[(VV0)1/3−1]2. (35)

This function is plotted in fig. 5 and increases rapidly with decreasing capsule volume .

The energy of a buckled solution of the shape equations can be calculated by numerical integration of the surface energy density of eq. (5). The data points for numerical solutions are shown in fig. 5. At a critical volume , the classical buckling volume, this solution branch starts to separate from the spherical branch (see inset in fig. 5). At first, it runs to the right and lies at slightly higher energies than the spherical solution branch (see the inset). These shapes correspond to unstable energy maxima. But after a return, the buckled branch crosses the spherical branch at a volume . From there on, the axisymmetric buckled configuration is energetically favourable to the spherical shape, representing the global energy minimum. The spherical shape is still metastable between and and represents a local energy minimum. Koiter’s stability analysis Koiter1969 () suggests that the buckling transition of real (imperfect) shells occurs somewhere in this region, depending on the severity of the imperfections. For , however, the spherical branch is unstable. This bifurcation scenario with metastable spherical and buckled branches below and above , respectively, is typical for a discontinuous shape transition. This summarises briefly the discussion of the buckling behaviour presented in ref. Knoche2011 ().

In Knoche2011 (), we also showed that our findings for the classical buckling volume coincide with the literature results for the classical buckling pressure (see, for example, Timoshenko1961 (); Landau1986 (); Ventsel2001 ()). To show that, we need to convert the pressure into a volume, by using the pressure-volume relation of the spherical solution branch. This can be derived from the elastic energy (35) via

 p=∂Usph∂V=2EH01−ν1R0[(VV0)1/3−1](V0V)2/3. (36)

Inverting this relation between and and inserting results in the classical buckling volume

 VcbV0=(12+√14+2(1−ν)√~EB)−3 (37)

or equivalently, for small Quilliet2012 (),

 ΔVcbV0≈6(1−ν)~E1/2B. (38)

This volume coincides very well with the point where the branch of axisymmetric buckled shapes separates from the spherical solution branch (see inset in fig. 5).

The elastic energy derived in the Pogorelov model, eq. (34), is also incorporated in the bifurcation diagram fig. 5. For volumes smaller than , it agrees well with the data points from the shape equations. Deviations start to develop for large deformations (), which is due to the simplification that the dimple was assumed to be small in the Pogorelov model. For shapes with too small dimples, , the model is also inaccurate, for two reasons. Firstly, it was assumed that the neighbourhood of the dimple which is deformed, and , is narrow. For too small dimples, this condition is not satisfied since the size of the inner neighbourhood becomes as large as the dimple itself. Secondly, for shapes with small dimples, the pre-buckling deformation (that is, the uniform contraction due to the negative inner pressure) cannot be neglected as it was done in the Pogorelov model. This is also the reason why the tensions and do not reach the correct limits far away from the dimple edge (see fig. 4 b).

However, the inset in fig. 5 suggests that the Pogorelov model can be successfully used to calculate the first buckling volume by setting . Doing so using (34) and (35), we obtain for small the critical volume difference for the first buckling,

 ΔV1stV0∣∣∣Pog=6J4/5min(1−ν)4/5(1−ν2)1/5~E3/5B. (39)

Equations (38) and (39) define two lines in a phase or stability diagram, see fig. 12, which is set up in the --plane. In the phase diagram, the line of the first buckling transition is represented by a continuous line according to the Pogorelov model (39) and data points which were derived from the shape equations by requiring equal energies of the spherical and buckled shape. Both approaches are in good agreement. The data points can be fitted with a power law (see fig. 12),

 ΔV1stV0∣∣∣shape eqs.=(4.78±0.03)~E0.6127±0.0006B (40)

This must be compared to (39) (evaluated at as used in fig. 12), which is and thus very close. The classical buckling transition is also represented by a continuous line, according to (37), and data points from the shape equations in the phase diagram fig. 12. In the space between the lines of first and classical buckling, spherical shapes and axisymmetric buckled shapes can both exist, since the spherical shapes are metastable and the buckled shapes stable.

With that we close our investigation of the axisymmetric shapes. In conclusion, we have shown that the hoop tension has a negative peak in the inner neighbourhood of the dimple edge. This is the region where wrinkles will develop when the compressive tension exceeds a critical value. Furthermore, we have established two critical volumes for the transition from the spherical shape to the axisymmetric buckled shape. These are the first buckling volume , at which the buckled shape becomes energetically favourable to the spherical shape, and the classical buckling volume where the spherical shape loses its (meta)stability.

## Iii Secondary buckling as wrinkling under locally compressive stress

### iii.1 Simplification of geometry and stress state

One main result of the previous section is that close to the dimple edge, a region of compressive hoop tension with a parabolic profile occurs, see fig. 4 b). This motivates our analysis of wrinkling or buckling of an elastic plate under a locally compressive parabolic stress in this section: We expect wrinkles to occur in the region of compressive hoop stress because the formation of wrinkles can release the compressive stress, which is energetically favourable Cerda2003 (); Wong2006 (). Within this relevant region, the capsule is shallow, and we can approximate it as a shallow shell or curved plate.

The specific plate geometry and the state of stress we impose are shown in fig. 6. The key features of the stress distribution and midsurface geometry in the wrinkling region can be reduced to the following simple functions: The stress distribution can be approximated by a parabola (dashed red line), and the section through the midsurface by a cubic parabola (dashed blue line). Note that these approximations only have to hold in the wrinkled region, i.e. the region of compressive , or, at the end of our analysis, the region with a large wrinkle amplitude. We need three parameters to describe the state of stress and the plate geometry in the following: the parameters and characterise the depth and width of the parabolic stress profile, and the parameter the curvature of the section through the midsurface.

The cubic parabola describing the midsurface is fitted to the point where the exact midsurface has vanishing curvature (see fig. 6). In the vicinity of this point, the real midsurface shows an approximately linear increase in curvature, . To obtain a cubic parabola with the same slope of curvature, we choose a height profile

 z(x,y)=16acy3 (41)

to describe the approximated midsurface. Here we have oriented the -coordinate along the original -coordinate and centred it at the point of vanishing curvature. By writing down (41), we have neglected the radius of curvature of the axisymmetric solution, because we assume that the wrinkles have shorter wavelength. Therefore, we consider a curved rectangular plate with the -coordinate corresponding to the uncurved -coordinate. Since the relevant portion of the shell is shallow, differences in the metric for the description with or as a coordinate can also be neglected Ventsel2001 (). So we can calculate by simple differentiation of with respect to instead of , which can be done numerically for a given axisymmetric solution,

 ac=dκsds0∣∣∣sc. (42)

The parabola to approximate the hoop tension is chosen to have the same minimum value and the same integral over the compressive part as the exact numerical function . Let denote the exact numerical integral, which has the physical interpretation of the net force in the compressive region . We can evaluate it by numerical integration for a given solution. A parabola of the form

 τx(x,y)=−τ0(1−apy2) (43)

has the roots . The integral over the parabola between its roots is and must equal . Thus we have

 ap=(4τ0/3F)2withF=∫s2s1τφ(s0)ds0 (44)

to determine the parameter of the parabola from a given axisymmetric shape. Note that the parabola is centred at , corresponding to the point of vanishing meridional curvature. This point does not exactly agree with the minimum of the exact hoop tension , but is very close (see fig. 6).

In the following, we will neglect the meridional tension , since it is small compared to , see fig. 4. Furthermore, there are no shear tensions in the axisymmetric configuration. Hence the stress state to which the plate is subjected reads , and according to (43).

Figure 7 summarises our simplified model. We want to investigate the wrinkling (or buckling) of a curved plate in the plane with height profile (41), subjected to a locally compressive stress (43). The appropriate tool for this task are the stability equations of shallow shells. The result of this analysis will be a wrinkling criterion in the form of a critical tension such that for compressive stress levels the plate wrinkles. The critical tension will depend on the parameters characterizing the width of the parabolic stress profile and characterizing the curvature of the midsurface section, and the elastic moduli of the shell.

### iii.2 Wrinkling criterion for the curved plate

Before presenting a more detailed stability analysis, we start with a scaling argument. Here we neglect curvature effects completely and approximate the compressed region of the plate by a rectangular strip of width (identical to the compressive part of the parabola) under a homogeneous compressive stress . For clamped long edges, the wrinkling wave length is given by the width, Timoshenko1961 (), and the resulting critical Euler buckling stress is . This result turns out to give the correct parameter dependence in leading order, see eq. (52) below.

The stability equations are partial differential equations for the normal displacement and Airy stress function and are given by Ventsel2001 ()

 Δ2ϕ =−EH0∇2κw (45) EBΔ2w =∇2κϕ+τx∂xxw+2τxy∂xyw+τy∂yyw, (46)

see appendix B for a derivation. Here, is the Laplacian and is the Vlasov operator. Note that the tensions and curvatures occurring in these equations are the tensions and curvatures of the initial state, prior to wrinkling. The Airy stress function, or stress potential, permits the calculation of the additional stresses in the plate, which arise because of the wrinkling, by the relations , and Ventsel2001 (). The existence of a non-trivial solution of these stability equations indicates the existence of an unstable deformation mode for the axisymmetric buckled solution (i.e. a negative eigenvalue of the second variation of the elastic energy).

In the present geometry and stress state, some terms vanish in the stability equations. We solve the equations using an ansatz

 w(x,y)=W(y)sinkx,ϕ(x,y)=Φ(y)sinkx. (47)

The wrinkle shape function represents a wrinkle pattern consisting of wrinkles extending in -direction with an amplitude function for each wrinkle, which are arranged in a periodic pattern in -direction with a wave number . Wrinkle shape and wave number are to be determined. Inserting this ansatz as well as the expressions for tensions and curvatures into the stability equations results in two coupled linear ordinary differential equations for the amplitude functions,

 0 =(∂4y−2k2∂2y+k4−k2EBτ0(1−apy2))W+ack2yEBΦ 0 =(∂4y−2k2∂2y+k4)Φ−(acyk2EH0)W. (48)

For a numerical solution, it is necessary to nondimensionalise the equations, which gives useful information on the relevant parameters. At first, we choose a length unit , which is the root of the parabolic stress profile, so that we can expect to decay on this scale. Substituting and in (48) induces further substitutions for the parameters of the differential equations such that they can finally be written in the form

 0 =(∂4^y−2^k2∂2^y+^k4−^k2^τ0(1−^y2))^W+(^ac^k2^y)^Φ 0 =(∂4^y−2^k2∂2^y+^k4)^Φ−(^ac^k2^y)^W. (49)

The substitutions included here are

 ^y =√apy, ^k =k√ap, ^ac =√EH0EBaca3/2p, (50) ^τ0 =τ0apEB, ^Φ =ΦEB, ^W =√EH0EBW.

A shooting method Stoer2010 (); numrec () can be applied to solve the differential equations numerically, when boundary conditions are provided. Because of the symmetry of the problem, we expect to be an even function. Then, from (49), it follows that is an odd function. Thus, it is sufficient to solve the differential equations on an interval . The boundary conditions at the left end of this interval, the start-point of integration, follow from the symmetry conditions, , and . The latter choice is arbitrary, since the differential equations are homogeneous. We imagine the plate to be infinitely large, so that the wrinkles are confined by the local nature of the compression rather than plate edges. For , the wrinkle amplitude has to approach , as well as the slope of the stress potential because the additional tension derived from shall approach . In practice, we integrate up to a sufficiently large and impose the boundary conditions and . When the shooting method must satisfy four boundary conditions at the endpoint of integration, it needs to vary four independent shooting parameters, which are typically the starting values of integration. But in the present case, due to the homogeneity of the differential equations, there are only three free starting conditions, , and ; the choice of is arbitrary and cannot serve as a shooting parameter. Thus, the shooting method must be allowed to vary one of the additional parameters of the differential equations, , or .

In fact, the differential equations can be interpreted as an eigenvalue problem, where one of the three parameters , or plays the role of an “eigenvalue” and must be chosen so that the equations do have a non-trivial solution. In our case, we choose as the eigenvalue, because this has the simplest physical interpretation: we increase the stress on the plate until a non-trivial solution in form of wrinkles exists. Then, this value is the critical stress which the plate can bear; for larger loads it will wrinkle. Obviously, the critical stress will depend on the other two parameters, . The value of is known when an axisymmetric solution is given, see (44) and (42). The wave number , however, is unknown. We are interested in the wrinkling mode which becomes unstable first, i.e. for the smallest possible load. To find this critical mode and the corresponding critical stress, we can minimise with respect to ,

 ^τc(^ac)=min^k^τ0(^ac,^k), (51)

so the non-dimensional critical stress only depends on . The wavenumber of the minimum is related to the wavelength of the wrinkles, . Due to boundary conditions in -direction, the wavelength is quantised, and thus the minimisation also has to be performed on a quantised space for . For simplicity, however, we will neglect the quantisation, which corresponds to an infinitely long plate (in -direction) on which all wavelengths are permitted.

A numerical analysis along the lines just presented reveals the function , which is shown in fig. 8 a) together with the critical wavelength. The critical tension increases with increasing curvature parameter . This reflects the well known fact that bent surfaces (like corrugated cardboard) are harder to bend in the perpendicular direction than flat surfaces. In addition, our numerical procedures return the shape of the wrinkles, that is, the amplitude function . As shown in fig. 8 b), the amplitude indeed decays rapidly outside the compressive region, i.e. for , as it is required for the approximation to be accurate.

The main result of this section is the formula for the critical compressive stress on which a curved plate will start to wrinkle,

 τc=apEB^τc(^ac)with^ac=√EH0EBaca3/2p. (52)

The function is known numerically (see fig. 8 a). For a given axisymmetric solution of the shape equations (or the Pogorelov model), the parameters and can be calculated according to (42) and (44), respectively, and inserted into (52). This gives the critical tension which the capsule can bear in the axisymmetric state. If the minimum value in the compressive region exceeds the critical stress, i.e. if , the axisymmetric capsule shape is unstable with respect to a wrinkling mode with wavelength

 λc=1√ap^λc(^ac), (53)

where the function is also known numerically (see fig. 8 a). This is our secondary buckling criterion according to the curved plate model.

Our further analysis also shows that the secondary buckling transition is a continuous transition, see appendix B. We checked that the fourth order terms in the elastic energy are positive, and that the energy change upon wrinkling reads for a given wrinkle amplitude (all numerical prefactors have been omitted). Thus, if the compressive tension exceeds the critical tension, , the formation of wrinkles lowers the elastic energy in second order, which yields a wrinkle amplitude growing like when the critical tension is exceeded. At the secondary buckling transition the system thus undergoes a supercritical pitchfork bifurcation. The continuity of the secondary buckling is also observed in the numerical simulations Quilliet2012 (). This is in contrast to the primary buckling transition, which is a discontinuous transition with metastability above and below the transition as discussed above and with an axisymmetric dimple of the buckled state which always has a finite size.

### iii.3 Secondary buckling transition threshold

We will now apply the secondary buckling criterion to buckled shapes from the shape equations and from the Pogorelov model in order to quantify the threshold for secondary buckling.

For the shape equations, we exemplarily discuss a capsule with elastic parameters and (see fig. 9). Our numerical analysis shows that the critical volume difference for the secondary buckling is , for the dimensionless curvature parameter we find , and the wrinkle wavelength is . The wrinkle wavelength can be converted into a wrinkle number: Since the wrinkles are centred at where the radius of the axisymmetric midsurface is , the number of wrinkles is given by , so either or . These results are illustrated in fig. 9 in a meridional section and in a three-dimensional view. In fig. 9 a), the same section as in fig. 6 is shown, but additionally the normal displacement, or wrinkle amplitude, is plotted. The results show that the wrinkle amplitude is very small outside the compressive region (compare with fig. 6), and that the approximated midsurface is quite accurate in the region of large wrinkle amplitude. This indicates that our approximations are justified, because their inaccuracies lie in regions where the wrinkles do not develop.

We have applied this analysis to a whole range of different bending stiffnesses , and thus generated a new line for the phase diagram fig. 12 which represents the critical volume of the secondary buckling transition (red data points). The data points can be fitted by the power law

 ΔV2ndV0∣∣∣shape eqs.=(2550±50)~E0.946±0.002B. (54)

Now, the secondary buckling criterion will be applied to the Pogorelov model. Because we have an approximate analytic expression for the compressive hoop stress profile within the Pogorelov model, we can obtain an analytic result for the secondary buckling threshold. At first, the parameters and must be determined, so that we can evaluate the critical hoop tension (52). In the Pogorelov model, the hoop tension is given by (31). Its minimum value is

 τ0=σmin3EH01−ν2ξαR0, (55)

and the integral between its roots can be calculated as

 F=EH01−ν2ξ2αR0[524σ