## Abstract

We compute nonlinear force equilibrium solutions for a clamped thin cylindrical shell under axial compression. The equilibrium solutions are dynamically unstable and located on the stability boundary of the unbuckled state. A fully localized single dimple deformation is identified as the *edge state*—the attractor for the dynamics restricted to the stability boundary. Under variation of the axial load, the single dimple undergoes homoclinic snaking in the azimuthal direction, creating states with multiple dimples arranged around the central circumference. Once the circumference is completely filled with a ring of dimples, snaking in the axial direction leads to further growth of the dimple pattern. These fully nonlinear solutions embedded in the stability boundary of the unbuckled state constitute critical shape deformations. The solutions may thus be a step towards explaining when the buckling and subsequent collapse of an axially loaded cylinder shell is triggered.

## 1. Introduction

Natural and engineered structures ranging from egg shells to air- and spacecraft are built from curved thin elastic shells that offer exceptional structural rigidity at minimal weight. Accurately predicting failure loads of those shells is thus an important yet not fully resolved classical problem of structural engineering. By far, most experimental and theoretical studies have focused on the prototypical problem of a cylindrical elastic shell under axial compression. A straightforward linear stability analysis finds the unbuckled state to lose stability at a critical buckling stress *λ*_{C}. Experimentally, however, buckling is observed at stresses as low as 20% of the critical value. This significant reduction in buckling load is quantified by the so-called knockdown factors [1] and has been linked to an extreme sensitivity of the buckling load on unavoidable imperfections of a shell. As first realized by Koiter [2], the extreme influence of imperfections can be related to nonlinear couplings of linear eigenmodes suggesting that weakly nonlinear analysis can capture the reduction in buckling load due to imperfections [3]. However, these approaches have not been fully successful in predicting failure loads of compressed shells.

More recently, fully nonlinear post-buckling solutions of the cylindrical shell equations were constructed. Those states are non-trivial force-equilibria that are speculated to be significant for understanding the shells’ route to collapse. The fully nonlinear equilibria may hence provide a path towards predicting realistic buckling loads and the shock sensitivity of shells [4,5]. For a cylinder shell, two types of fully nonlinear solution families have been found that differ from commonly studied eigenmodes in that they are localized to a part of the shell:

First solutions were presented by Hunt, Lord, Peletier, Champneys and co-workers [6–11]. Based on a Galerkin approximation in azimuthal (circumferential) direction, the equilibrium shell equations were converted into a set of coupled ODEs in the axial space coordinate. Assuming strict periodicity in the circumferential direction together with asymptotic decaying boundary conditions ensuring that the solutions smoothly approach the unperturbed cylinder on the top and bottom, the equilibria were captured as periodic orbits that can be computed and followed using the continuation package AUTO. The seminal work yields solutions that are localized in the axial direction and consist of a ring of dimples arranged periodically around the circumference. By carefully varying the load, the solutions are shown to undergo homoclinic snaking. This means they undergo a sequence of saddle-node bifurcations that create states with additional rings of dimples. Thereby the dimple pattern grows in the axial direction, but it remains localized in that direction and strictly periodic along the circumference.

A second type of equilibrium solutions that appear fully localized not only in axial but both in axial and circumferential directions were computed by Horák, Lord & Peletier [12,13] using methods to find *mountain-pass points* [14,15]. Specifically, a single dimple solution was constructed. However, instead of asymptotically decaying boundary conditions as in the previous work, symmetry boundary conditions in the axial direction were chosen, equivalent to periodic boundary conditions with a restriction to even functions. These boundary conditions significantly simplify the analysis, but as a result, the single dimple appears localized in two dimensions, but is actually supported by a column-like deformation extending over the whole cylinder length (see also appendix A). The extended column deformation is permitted by the boundary conditions, which do not force the deformation to vanish at the loaded ends. The single dimple state was found by exploiting the variational structure of elasticity, namely the fact that the dynamics can be derived from an energy functional. The single dimple is a mountain pass, a saddle point of the energy landscape located between the unbuckled state and the post-buckling regime. Additional mountain pass points, consisting of multiple dimples, have been computed with different initial conditions of the algorithm. But contrary to the rows of dimples in the azimuthally periodic subspaces discussed above, no connection between the different states is observed; in particular, homoclinic snaking of the single dimple was not observed in [12].

To further investigate the relation of fully nonlinear equilibrium solutions and experimentally observed buckling loads, such states need to be constructed for realistic loading scenarios and boundary conditions. Instead of experimentally hard to realize periodic boundary conditions, we therefore consider an axially loaded cylinder shell of finite length subject to clamped boundary conditions. Assuming overdamped dynamics, we construct attracting equilibrium solutions on the basin boundary of the unbuckled state’s basin of attraction. These so-called edge states are dominated by a fully localized single dimple solution related to the mountain pass state constructed by Horák *et al.* Continuing the solution branch for varying loads, we observe homoclinic snaking in the azimuthal direction: the single dimple undergoes a sequence of saddle-node bifurcations leading to an increasing number of dimples localized in the circumferential direction. Once the whole circumference is filled by a periodic ring of dimples, the onset of axial snaking as reported by Hunt, Lord, Peletier, Champneys and co-workers [6–11] is observed.

The paper is organized as follows: §2 introduces the equations describing the compressed shell as well as the numerical representation of the equations and the algorithm used to find the edge state. The results are presented in §3, starting in §3a with a description of the edge state—a single localized dimple—and its properties. The effect of load variations on the single dimple state is discussed in §3b. Finally in §3c, we restrict the domain to a periodic section of the cylinder and document homoclinic snaking of the state in the axial direction. The significance of the results is discussed in §4. In appendix A, we provide a verification of our numerical scheme by recomputing the mountain pass solution from [12,13].

## 2. Equations and numerical implementation

### (a) Equations and boundary conditions

We consider a thin cylindrical shell of thickness *t*, radius *R* and length *L* subject to uniform axial load *λ*. The Young modulus *E*, and Poisson ratio *ν* characterize the material of the elastic shell. Let *x* and *y* denote the axial and the azimuthal coordinate on the domain *Ω*=(−*L*/2,*L*/2)×(−*πR*,*πR*), respectively. The elastic response of the shell is faithfully captured by the Donnell–Mushtari–Vlasov (DMV) equations, also known as the von Kármán–Donnell equations. These small strain and moderate rotation asymptotic reductions of three-dimensional elasticity are approximations for thin shells in the nonlinear regime up to rotations of 15–20^{°} [16]. The state of the shell is characterized in terms of the outward displacement normal to the surface *W*(*x*,*y*) and the Airy stress function *Φ*(*x*,*y*). The latter is a potential of the stresses in the shell and related to them by *Φ*_{,xx}=*N*_{yy}, *Φ*_{,yy}=*N*_{xx} and *Φ*_{,xy}=−*N*_{xy} (subscript comma denotes differentiation), where *N* is the in-plane stress integrated over the shell thickness. Under load, the shell responds with a uniform axial compression and radial expansion described by the axisymmetric base state with constant normal displacement *W*^{0}=*λRν*/*tE* and stress potential *w*=*W*−*W*^{0} and *ϕ*=*Φ*−*Φ*^{0}, the DMV equations take the form [16]
^{2}*A*=*A*_{,xxxx}+2*A*_{,xxyy}+*A*_{,yyyy} and the bracket operator
*D*=*Et*^{3}/12(1−*ν*^{2}), Young’s modulus *E* and Poisson’s ratio *ν* characterize the elastic shell material. *λ* is the axial compressive load per unit length. For small values of *λ*, the unbuckled configuration *w*(*x*,*y*)=0 is linearly stable. A linear stability analysis of the perfect, unbuckled, infinitely long cylindrical shell predicts a critical load of (see e.g. [1])

The equations (2.1) and (2.2) are nonlinear differential equations of fourth order and thus need to be supplied with two boundary conditions for each field on each side of the domain. To most faithfully reflect experimental conditions where a cylinder shell is placed between two compressing plates, we use clamped boundary for the normal displacement
*R*+*W*^{0}. Following [17] for the stress function at the loaded ends, we enforce

Equations (2.1) and (2.2) together with the boundary conditions are invariant with respect to reflections *Ω*_{1/4}=(0,*L*/2)×(0,*πR*) and enforce the symmetry conditions at the planes *x*=0 and *y*=0. The symmetry conditions are equivalent to demanding odd derivatives of the solutions to vanish at the symmetry planes. Exploiting these discrete symmetries reduces the numerical cost of the simulations and improves convergence properties of root-finding algorithms. The full set of boundary conditions in the symmetric subspace is as follows:
*ϕ* nor Δ*ϕ* is fixed in their absolute value by these boundary conditions, additional integral conditions need to be imposed

### (b) Numerical implementation

To compute equilibrium solutions of (2.1), perform linear stability analyses of the solutions and follow solutions as a function of load, we consider a first-order dynamical system resulting from the overdamped evolution of the shell.

The r.h.s. of (2.1) indicates unbalanced forces in the shell. If the forces vanish identically, an equilibrium solution is found. To characterize the stability of an equilibrium, we consider the temporal evolution of an infinitesimal perturbation. A non-zero r.h.s. will cause the shell to move. In general, elastic forces are balanced by both acceleration *w*_{,ττ} in the normal displacement and an additional damping term −*δw*_{,τ} capturing dissipation. As any acceleration from rest causes a velocity in the same direction, the stability of an equilibrium with respect to infinitesimal perturbations is not altered by inertia. Thus, for characterizing the stability of an equilibrium, we may neglect acceleration and consider an overdamped dynamics, where −*δw*_{,τ}≫*w*_{,ττ}. This overdamped situation may physically be realized by submerging the shell in a highly viscous liquid. It is the simplest dynamics that allows computing equilibria and characterizing their stability but for describing inertial effects such as shell vibrations, the acceleration term would need to be included. The normal displacement evolves according to the nonlinear dynamical system
*w*, obtained by solving the boundary value problem (2.2) for a known *w*.

Computing solutions to the equilibrium equations (2.1) is equivalent to finding fixed points *f* linearized around a fixed point *w** have negative real part, the fixed point is stable and small deviations will exponentially return to *w** under time evolution. The undeformed state *w*=0 is a solution *f*(*w*=0)=0 and all eigenvalues of *f* linearized around 0 have negative real part for *λ*<*λ*_{C}.

The equations are discretized on a two-dimensional, rectangular grid and derivatives are approximated by high-order finite differences. This discretization approach combines the flexibility of finite differences when modifying boundary conditions with high accuracy of the derivatives comparable to spectral methods. The latter is of particular importance for capturing the fourth-order spatial derivatives in the shell equations. To avoid Runge oscillations at the non-periodic and non-symmetric domain boundaries, a non-uniform grid with points more densely clustered towards the boundaries is chosen. Specifically, with *M* the stencil size, the first *M* points at the boundary are distributed proportional to the roots of the 2*M*th Chebyshev polynomial. Interior points, with a distance larger than *M* from the boundaries, are uniformly distributed.

High-order finite difference stencils for one-dimensional derivatives are computed on the non-uniform grid using Fornberg’s algorithm [18]. The matrix representations of the two-dimensional derivative operators Δ, Δ^{2} and ∂_{xy} are computed by sums and products of the matrix representations of the one-dimensional operators.

The first step in computing *f*(*w*) is to compute the Airy stress function *ϕ*(*w*) for given normal displacement *w*. To solve the corresponding boundary value problem in (2.2) with the boundary conditions (2.6), an auxiliary variable *ψ*=Δ*ϕ* is introduced. Thereby, the fourth-order equation is converted into two second-order equations

Let *w*, *ϕ* and *ψ* on the discrete grid and let *D*_{ab} and *D*_{Δ} be matrix representations of the derivative operators of ∂_{ab} and Δ, respectively. To solve the Poisson equation (2.12) for *D*_{Δ}. Technically, lines corresponding to boundary points are replaced with the finite difference stencil for first derivatives and the entries of *D*_{Δ} and following the same outlined procedure. With

The bi-Laplacian Δ^{2} contains fourth-order spatial derivatives. As a result, the equation *w*_{,τ}=*f*(*w*) is stiff and explicit time integration requires impractically small time steps *δτ*. Owing to nonlinearities, fully implicit time integration is, however, very expensive. Therefore, we choose a semi-implicit integration scheme, where the linear parts of the equation are treated implicitly and the nonlinear parts explicitly. Specifically, semi-implicit backwards differentiation (SBDF) of order three [19] with first order SBDF for initialization proved to be accurate. Note that the matrix *D*_{Δ} that needs to be inverted for computing

Our code is fully object-oriented and written in C++. Basic vector- and matrix operations as well the the direct solvers are based on the C++ library ‘Eigen’ [20].

### (c) The edge tracking algorithm

Starting from a given shell deformation, the time evolution follows a unique path described by the first-order dynamical system (2.11) representing the overdamped dynamics of the shell. The current state of the system is, therefore, completely described by the shell’s deformation; the systems’ state space—the space of all possible states—is the space of all possible deformations *w*(*x*,*y*) that satisfy the boundary conditions.

If the compressive load is below *λ*_{C}, small deformations will decay towards the unbuckled, straight state. The set of all these deformations that return to the unbuckled state is its *basin of attraction*, a region in state space that surrounds the unbuckled state. States that lie outside this basin will *not* return to the unbuckled state and almost all of them will evolve towards collapse. Consequently, there is a boundary between these states and the basin of attraction of the unbuckled state—states in this boundary are the objective of this study as they may guide the transition to collapse.

The chosen overdamped dynamics can be derived from a variational principle such that the evolution follows the gradient of an energy functional. Consequently, the basin boundary can also be interpreted as a stucture in the system’s energy landscape: In the energy landscape, the linear stability of the unbuckled state is reflected by the fact that it is located at the bottom of a valley. The basin boundary corresponds to a mountain ridge surrounding it. States on the boundary will slide down neither on one side of the ridge towards the unbuckled state nor on the other side towards collapse. They may, however, still evolve *on* the basin boundary until they reach a saddle point. Such a saddle point is an equilibrium and was termed *mountain pass* by Horák *et al.* [12,13]; in terms of dynamical systems theory, it is an *edge state* [21,22] of the first-order system (2.11). A defining property of those edge states is that they only have one unstable direction transversal to the basin boundary, while all directions within the boundary are attractive. That implies the stable manifold of an edge state to be of co-dimension one. Thus, it locally divides state space.

An algorithm to compute the edge state, termed *edge tracking*, has been developed and successfully employed in fluid mechanics [23,24]; it allows to follow the time evolution of a state in the basin boundary, which will evolve towards the attracting state in the boundary, i.e. the edge state. The same algorithm is used here to compute the edge state of the axially compressed cylinder: Let *w*_{0} be a random deformation of the shell strong enough to trigger the evolution towards collapse. Now, consider the set of scaled deformations *α*⋅*w*_{0}. For a prefactor *α*=1, the collapsing deformation outside the basin boundary is recovered, while *α*=0 corresponds to the stable unbuckled state. Consequently, there must be a value 0<*α**<1 such that *in* the basin boundary and neither leads to collapse nor evolves towards the unbuckled state—the trajectory starting from *α** can only be approximated. Using bisection, the edge tracking algorithm, therefore, computes two values *α*_{L} and *α*_{H} that bracket the ‘true’ value *α**, so that *w*_{0,L}=*α*_{L}*w*_{0} returns to the unbuckled state, while *w*_{0,H}=*α*_{H}*w*_{0} collapses. Following the time evolution of both *w*_{0,L} and *w*_{0,H}, the true trajectory in the basin boundary is well approximated as long as the distance between *w*_{0,L}(*τ*) and *w*_{0,H}(*τ*) remains small. Once the distance between both approximating trajectories becomes too large at *τ*_{0}, time integration is stopped and a new pair of states that bracket the edge is computed by scaling *w*_{1}=*w*_{0,H}(*τ*_{0}). The two states that result from this bisection, *w*_{1,L} and *w*_{1,H}, are then again integrated forward in time until their separation becomes too large. Iterating this process, the trajectory *in* the basin boundary can be approximated for arbitrarily long times, eventually converging towards the edge state.

### (d) Root finding and numerical continuation

Equilibrium solutions to the DMV equations (2.1) are computed using a Newton root finding method. We specifically use the optimized Newton–Krylov–Hookstep algorithm from John Gibson’s ‘channelflow’ package [25]. The method operates in lower dimensional Krylov-subspaces of the full space and is matrix-free. Consequently, large systems of nonlinear algebraic equations in millions of coupled degrees of freedom can be solved. To increase the radius of convergence of the Newton method, we further employ advanced trust region optimization techniques [26]. Based on the experience in fluid mechanics and following the methodology of channelflow, we do not compute roots of (2.1) directly. Instead, roots of *F*^{T}(*w*)−*w*=0 are constructed, where *F*^{T} denotes the forward evolution for time *T*. Obviously, finding roots of the finite time difference equation and of the direct equation, *f*(*w*)=0 is equivalent. We have, however, observed that the Krylov subspace methods converge faster for the finite time difference equation, presumably because stable eigendirections are damped by the finite time integration. The Newton–Krylov–Hookstep algorithm applied to the finite time difference equation provides a robust method for finding both stable and unstable equilibrium solutions from initial guesses such as those constructed by edge tracking.

To continue equilibrium solutions under parameter variation, we use numerical continuation techniques. The predictor corrector method from channelflow combines a prediction step for a new parameter value obtained by quadratic extrapolation with correction by Newton–Hookstep. Thereby, the full bifurcation behaviour of equilibrium solutions can be studied.

## 3. Localized post-buckling states

We choose parameters corresponding to the seminal experiments of Yamaki [27],
*et al.* [7]. In the following, we non-dimensionlize the problem by choosing the thickness *t* as a length scale and *Et*^{3} as an energy scale. The unit of time is set by the damping constant *δ*, which can be chosen arbitrarily.

Owing to symmetries, the simulations are carried out on a subdomain consisting of one quarter of the cylinder. Finite differences of order 14 and a spatial resolution of *N*_{x}×*N*_{y}=76×300 for the quarter domain are chosen. Doubling the spatial resolution and varying finite difference orders, we confirm that all observables are fully converged with relative errors below 10^{−13}.

### (a) The edge state—the attractor on the basin boundary

In the following, we consider the cylinder at a fixed load of *λ*/*λ*_{C}=0.7 and compute the equilibrium deformation corresponding to the edge state. The edge tracking algorithm is initalized with a deformation consisting of a random superposition of sine functions up to wavenumbers *k*_{x}=7 and *k*_{y}=15, filtered to satisfy the clamped boundary conditions. The approximating trajectories starting from *w*_{L} (*w*_{H}), indicated by blue (red) lines in figure 1, bracket a trajectory in the basin boundary. This trajectory settles to a constant energy value indicating that it approaches an attracting fixed point: the edge state. With the deformation reached by edge tracking as an initial guess, the Newton algorithm quickly converges, indicating that a locally attracting equilibrium solution of (2.1) embedded in the basin boundary is found.

The edge state is visualized in figure 2*a* where the normal deflection is shown on the unrolled cylinder. Figure 2*b* shows a three-dimensional rendering of the edge state. The attractor in the basin boundary corresponds to a single dimple inward indentation sandwiched between two small outward buckles on its right and left. The deformation field is quickly damped to zero towards the clamped boundaries in axial direction and also tends to zero in the azimuthal direction; the dimple is, thus, a fully localized equilibrium deformation of the shell. Figure 2*c* shows a cut through the dimple in the axial direction and figure 2*d* a cut in the azimuthal direction, confirming the localization of the state. In azimuthal direction, the tails show damped spatial oscillations, while the decay is monotonic in the axial direction.

Linearizing (2.11) around the single dimple solution, we perform a linear stability analysis of perturbations to that state. In the symmetric subspace, the edge state has a single unstable eigenvalue with a positive real part of 0.0243, all other eigenvalues have strictly negative real parts. Relaxing the symmetry restrictions does not affect the unstable eigenvalue, but gives rise to a neutrally stable mode that corresponds to rotating the deformation around the cylinder axis. Owing to the clamped boundary conditions, there is no translational invariance in the axial direction even when no symmetry restrictions are enforced. As a consequence, the eigenmode corresponding to shifts in the axial direction has a negative eigenvalue of −3.25×10^{−5}, causing a dimple displaced in the axial direction to return to its central position under time evolution. Though computed with symmetry restrictions, the single dimple is hence also an edge state for the full cylinder, as it has a single unstable direction irrespective of symmetry constraints. Perturbations in the direction defined by this unstable eigenmode will return to the unbuckled state or—depending on the sign—lead to collapse. The single unstable direction, thus, points out of the basin boundary and is controlled by the edge tracking algorithm. All remaining directions are stable, rendering the edge state attractive for the dynamics within the basin boundary.

Note that in addition to the single dimple solution, other attracting equilibria are embedded in the basin boundary. Which attractor the edge tracking algorithm converges to depends on the choice of the (arbitrary) initial deformation. These alternative edge states consist of several isolated dimples distributed over the cylinder, similar to the multi-dimple mountain pass points discussed in [12]. The single dimple configuration, is however, the attractor in the basin boundary that most initial perturbations converge to. Among all states we found, it is also the state with smallest elastic energy and smallest norm ∥*w*∥_{2}.

### (b) Varying the load: homoclinic snaking

The single dimple state has been found by edge tracking at 70% of the critical load. How the equilibrium solution changes under variation of load is investigated by numerical continuation of the solution branch. The bifurcation diagram obtained by this procedure is presented in figure 3, where both the elastic energy and maximum deflection are used to characterize the states. The single dimple state appears to bifurcate subcritically from the unbuckled state at *λ*=*λ*_{C}. Thus, both energy and maximum deformation decrease when approaching the critical load. When lowering the load, the solution branch turns around in a saddle node bifurcation at *λ*/*λ*_{C}=0.48. Subsequently, a sequence of further saddle-node bifurcations give rise to altogether five folds in the solution branch. These folds, peaked on one end and more rounded on the other, tilt and align with the *y*-axis, indicating that large changes in the solution are observed for decreasing intervals in load.

The physical mechanism underlying the folds in the bifurcation diagram is revealed by visualizing the deformation field. In figure 4, the equilibrium configuration is shown at the subsequent minima of the load (indicated by the red circles in the bifurcation diagram). At the first minimum, the equilibrium shows a single localized dimple. After one fold at the second minimum, two additional dimples have formed right and left of the original dimple. Subsequently, after each fold, an additional pair of dimples appears, until finally the whole cylinder circumference is filled with 10 dimples. The marginal eigenmode associated with the saddle-node bifurcation between the three and the five dimple state is shown in figure 4*f*. It is weighted most heavily at the fronts of the three dimple state so that adding a small component of the eigenmode creates a new pair of dimples and thereby widens the structure. This bifurcation structure in which a sequence of saddle-node bifurcations leads to the growth of localized states by subsequent addition of structures at the fronts of a localized state is known as homoclinic snaking [28]. The mechanism underlying homoclinic snaking has been most carefully studied in pattern-forming PDE models such as the Swift–Hohenberg equation [29], but the phenomenon can be observed in many nonlinear systems (see [28] for a recent review). Examples include flows driven by convection [30] or shear [31,32] as well as forced nematic crystal layers [33]. Homoclinic snaking was also observed for the axially localized, circumferentially periodic solutions of the compressed cylindrical shell from [6–11] discussed in the introduction. Based on studies of the Swift–Hohenberg equation and various fluid problems, we expect the existence of additional solution branches corresponding to states composed of an even number of dimples. Such states with even number of dimples will not be connected to the single dimple edge state, but together with asymmetric configurations connecting states of even and odd number of dimples may form a full snakes-and-ladders structure observed in other systems.

Figure 5 presents profiles of the single dimple and the 10 dimple state from figure 4. The amplitude of the single dimple is smaller than the 10 dimple state reflecting the higher load at which it has been computed. While the inner region of the 10 dimple state appears periodic, a small defect in the periodic pattern leads to weakening of the outer two dimples. As a consequence, the solution branch does not connect to the circumferentially periodic state of wavenumber 10 discussed below. In other systems showing homoclinic snaking on a periodic domain, reconnections of the localized branch and the periodic state are observed only for some parameter values, while other system parameters lead to the formation of defects [34] as observed here for the cylinder with parameters from Yamaki’s experiments. For which parameters the branch reconnects is related to the compatibility of the internal wavelength of the localized structure with the imposed 2*πR* period in azimuthal direction. As the naturally chosen azimuthal wavelength of the localized patterns is close to 2*πR*/10, we expect a reconnection for the slightly modified *R*/*t*.

To characterize the stability of the states, spectra of the linearized evolution operator are calculated along the branch. Colours in figure 3 represent the number of unstable eigenvalues. The single dimple has a single unstable eigenvalue from its initial bifurcation at the critical load down to the saddle-node bifurcation point at *λ*/*λ*_{C}=0.48. On the upper branch of the saddle-node, the single dimple state is stabilized, even if the symmetry restrictions are relaxed. Thus, in a small parameter region, the state can be found without edge tracking or the mountain pass algorithm by standard time integration. Further following the branch, the state picks up additional unstable directions. Most multi-dimple states have more than one unstable eigenvalue and thus cannot be constructed by edge tracking.

### (c) Axially localized post-buckling in a periodic subspace

Under continuation in load, the single dimple solution grows in the circumferential direction. The dimple structure approaches a state of 10-fold azimuthal symmetry but for the parameters from Yamaki’s experiment, a developing defect (figure 5) indicates that the snaking branch does not fully connect to a periodic state. To analyse such a 10-fold azimuthally periodic state, we follow Hunt *et al.* [8] and consider the invariant subspace S10 of 10-fold circumferential symmetry. Technically the computational domain is restricted to an azimuthal slice of size 2*πR*/10 subject to periodic boundary conditions in *x*.

Edge tracking in the symmetry subspace at *λ*/*λ*_{C}=0.7 followed by a Newton search yields the state shown in figure 6*a* with the same colour coding as in figure 2. The individual dimples of the periodic state appear geometrically almost identical to the localized single dimple and the amplitude measured in terms of maximum deflection is virtually indistinguishable.

The bifurcation behaviour in subspaces of discrete rotational symmetry has been discussed in great detail in [6–9], where asymptotically vanishing boundary conditions in axial direction are considered. The cylinder was shown to support symmetric and so-called cross-symmetric dimple patterns. The latter ones are obtained if a shift-and-reflect symmetry is imposed at the equator instead of the reflection symmetry. Under continuation in load, the azimuthally periodic states undergo homoclinic snaking in the axial direction. Starting from a single, axially localized ring of dimples, additional rings subsequently are added in axial direction, leading to the growth of a checker-board pattern.

For the clamped boundary conditions considered here, a similar behaviour is observed as the load is varied. The bifurcation diagram in figure 7 shows the solution branch in the S10 subspace together with the results obtained for the localized dimple. In the S10 subspace, one axial snaking step can be observed. The visualization of the states at the minima of *λ* in figure 6*b*,*c* reveals a pair of additional rings of dimples being added above and below the state found by edge tracking. The additional rings are displaced by half a wavelength in the azimuthal direction so that eventually a checker-board configuration is reached. Owing to the clamped boundary conditions and the finite length of the cylinder, no further snaking steps can be observed.

Comparing the bifurcation curves for the localized and the periodic state is particularly enlightening for the maximum deflection. Both states have identical deflection *min*(*w*) for

## 4. Discussion

We have identified edge states—fully nonlinear equilibrium states on the basin boundary of the unbuckled state’s basin of attraction—for an axially loaded cylindrical shell subject to realistic boundary conditions. For the parameters used in Yamaki’s experiments [27] and a load of 70% of the critical load, this edge state is a fully localized single dimple solution. The depth of the dimple is approximately 1.5 times the thickness of shells.

Under continuation in load, the single dimple undergoes homoclinic snaking in the circumferential direction. Next to the single dimple, additional dimples are subsequently created until the circumference is filled with 10 dimples. In the subspace of 10-fold azimuthal symmetry, the state—now only localized in the axial direction—then shows the onset of axial snaking as reported previously by Hunt, Lord, Peletier, Champneys and co-workers [6–8]. Thus, similar to the two-dimensional Swift–Hohenberg equation [35], the equations for a loaded cylindrical shell support homoclinic snaking of localized dimple patterns in two dimensions.

For most parameters, all equilibrium solutions constructed by edge tracking and parameter continuation are unstable and located on the basin boundary of the unbuckled state. This implies that a slightly smaller deformation will return to the unbuckled state, while a slightly larger one leads to collapse. In this sense, the unstable equilibria define the most ‘dangerous’ deformations of the cylinder. This is especially relevant for edge states that are attracting for the dynamics confined to the basin boundary. The edge states coincide with mountain pass points, which are (locally) the energetically lowest points that allow leaving the unbuckled states attraction basin. As all equilibria constructed physically correspond to an arrangement of multiple dimples, each requiring elastic energy, we may conjecture that the single dimple edge state is not only a local but also the global energy minimum on the basin boundary.

With the single dimple being an edge state and the energy minimum on the basin boundary of the unbuckled state, the state’s energy is an indication of how much a perfect shell needs to be deformed before it buckles and collapses. The edge state’s energy approaches zero as the applied load approaches the critical value (figure 3). This reflects the observation that smaller deflections are sufficient to destabilize a shell at higher load. Below loads of *λ*/*λ*_{C}≃0.48, the single dimple does not exist as an equilibrium configuration and the edge state is instead (at least for a certain range of loads) a configuration with three dimples. This suggests that depending on the load, a deformation in the form of a single dimple might be less efficient in triggering a shell to buckle than a deformation of more complex geometric shape.

The fully nonlinear equilibrium states and specifically edge states may define the shape and amplitude of shell deformations that guide the transition to buckling. This interpretation is backed by the obvious similarity between the dimple patterns presented here and the onset of post-buckling dynamics in experiments [27]. How to relate the fully nonlinear states to experimentally observed buckling thresholds is not yet entirely obvious and will require additional parameter studies, the construction of many more unstable equilibria and potentially a statistical approach. Nevertheless, we expect unstable fully nonlinear localized equilibrium states to be an important element of a fully nonlinear description of when and how elastic cylinder shells buckle.

## Data accessibility

This article has no additional data.

## Authors' contributions

T.M.S. conceived the work. T.K. developed the numerical tools and carried out the simulations in consultation with T.M.S. Both authors interpreted the results, wrote the paper and gave their final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the Swiss National Science Foundation under grant nos. 200021-165530 and 200021-160088.

## Acknowledgements

We thank Hecke Schrobsdorff for computational support and John Hutchinson for very valuable discussions and comments on the manuscript.

## Appendix A. Verification of the numerical code against previous results

The only non-trivial solutions of the DMV equations on a full cylinder without azimuthal symmetry constraints in the literature are, to our knowledge, the solutions of Horak *et al.* [12,13]. The authors employ a mountain pass algorithm to find saddle points on a cylinder with symmetric boundary conditions in axial direction.
*et al.* employ different representations of the derivatives ∂_{xy} including a computation in Fourier space that yields the best results.

To verify our code, we compute the single dimple solution from §3a on the same domain and with the same set of boundary conditions as in [13] at an axial load of *λ*/*λ*_{C}=0.7. The best quantitative data in [13] are found in figure 7, where the normal displacement *w*(*x*=0,*y*) (defined as the inward displacement by the authors) is shown. Note that the figure reveals a strong dependence on the method by which the derivative ∂_{xy} is computed, indicating that the grid is very coarse; for the comparison we chose the solution where this derivative is computed in Fourier space, assuming to be the most accurate. In figure 8, we plot the solution of [13] together with the single dimple state computed with our code (with a *N*_{x}=*N*_{y}=100 and finite differences of order 10 on the quarter-domain). The figure shows a near-perfect agreement between the two states, confirming that we are able to accurately reproduce the former results.

- Received March 15, 2017.
- Accepted August 11, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.

## References

Sign in for Fellows of the Royal Society

Fellows: please access the online journals via the Fellows’ Room

Not a subscriber? Request a free trial

### Log in using your username and password

### Log in through your institution

Pay Per Article - You may access this article or this issue (from the computer you are currently using) for 30 days.

Regain Access - You can regain access to a recent Pay per Article or Pay per Issue purchase if your access period has not yet expired.