Spatial chaos as a governing factor for imperfection sensitivity in shell buckling

Shell buckling is known for its extreme sensitivity to initial imperfections. It is generally understood that this sensitivity is caused by subcritical (unstable) buckling, whereby initial geometric imperfections rapidly erode the idealized buckling load of the perfect shell. However, it is less appreciated that subcriticality also creates a strong proclivity for spatially localized buckling modes. The spatial multiplicity of localizations implies a large set of possible trajectories to instability—also known as spatial chaos—with each trajectory afﬁne to a particular imperfection. Using a toy model, namely a link system on a softening elastic foundation, we show that spatial chaos leads to a large spread in buckling loads even for seemingly indistinguishable random imperfections of equal amplitude. By imposing a dominant imperfection, the strong sensitivity to random imperfections is ameliorated. The ability to control the equilibrium trajectory to buckling via dominant imperfections or elastic tailoring creates interesting possibilities for designing imperfection-insensitive shells.


I. INTRODUCTION
It is well known that systems governed by subcritical bifurcations (unstable branching) are sensitive to initial imperfections.The unfolding of a subcritical pitchfork bifurcation in a two-parameter system reveals a cusp, oftentimes with a 2/3 power-law relationship between the driving and perturbing parameters [1].Hence, external perturbations or geometric imperfections round off the perfect bifurcation into limit points (folds) with rapidly decreasing magnitude of the critical parameter.
In engineering mechanics, shell buckling is quintessentially imperfection sensitive.Particularly axially compressed cylinders and externally pressurized spheres show a proclivity to buckle (and often collapse) at loads much smaller than the linear instability prediction (P cl ).Even for shells manufactured using the same process, controlled to the same geometric precision, and tested with the same apparatus, a large scatter in buckling loads is observed [2].Hence, the theoretically deterministic system turns out to be highly stochastic in practice.
Koiter [3] famously showed how geometric imperfections of the order of the shell thickness can drastically erode the idealized buckling prediction to levels observed in experiments (≈ 0.2P cl ).Geometric imperfections spanning the vector space of critical eigenvectors are particularly conducive to premature buckling [4].Within this framework, stochastic variations in observed buckling loads are explained by different imperfection modes and varying imperfection amplitudes.However, it is less clear what drives the large Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
scatter in buckling loads for similarly manufactured shells with seemingly indistinguishable random imperfections of comparable amplitude (see, e.g., 0.45P cl -0.8P cl in Ref. [5]).In other words, current imperfection sensitivity analyses use eigenmodes as imperfections and vary their amplitudes to modulate the onset of instability.The question of how random imperfections of the same amplitude lead to a large scatter in buckling loads is addressed here.
Recent applications of nonlinear systems theory to shell buckling [6][7][8] have revealed that the prebuckling equilibrium of an axially compressed, isotropic cylinder is surrounded by a large set of energetically favourable gateways to instability.On an energy landscape, the prebuckled state is described by a well, which is enclosed by a stability boundary (a mountain ridge) that can be traversed by larger (noninfinitesimal) perturbations.States on either side of the stability boundary evolve to the prebuckled state or into the postbuckling regime.Local minima along this stability boundary (saddle points in the multidimensional energy landscape) represent energetically favorable escape points from the prebuckling well and are known as mountain-pass points [6] or edge states [7].For the axially compressed cylinder, the smallest mountain pass corresponds to a fully localized single dimple in the cylinder wall [6].
Due to the multivalued nature of localizations, the single dimple can, with equal likelihood, develop anywhere across the spatial domain; that is, for an isotropic cylinder, the total strain energy is invariant to translations of a dimple of fixed size.Furthermore, a large set of many-dimpled edge states, with equally associated multiplicity, also exists on the stability boundary [8].Consequently, there exists a large set (for an infinitely long continuum, strictly, an infinite set) of localized escape routes from the prebuckling well over the stability boundary.Depending on the precise nature of initial imperfections and external perturbations, the cylinder can therefore transition out of the prebuckling well via one of these many routes.[13]).The stable prebuckling path is bounded by a set of localized unstable edge states.(b) Resistance of the cylinder to concentrated lateral probing (poking force F vs. indentation displacement w) for different levels of axial compression u.Plots adapted from Ref. [8].
This situation of an infinite set of escape routes-with each trajectory to buckling sensitively dependent on the initial conditions-has certain analogies with temporal chaos and is therefore often referred to as "spatial chaos" [9][10][11][12].Here we demonstrate how spatial chaos is a governing factor in the stochasticity observed in shell buckling, i.e., the large spread in buckling loads for seemingly indistinguishable imperfections of the same amplitude.

II. SPATIAL CHAOS IN SHELL BUCKLING
On an equilibrium diagram of applied end-shortening (u) vs. reaction force (P), an axially compressed, thin-walled, isotropic cylinder follows a quasilinear prebuckling path that destabilizes at a bifurcation [see Fig. 1(a)].Many unstable equilibrium paths of edge states embedded in the stability boundary exist immediately adjacent to the stable prebuckling path.As shown in Fig. 1(a), many of these edge states correspond to localizations.From a static equilibrium point of view, these single-dimple and multiple-dimple modes are circumferentially and, if away from the structure's boundary, axially invariant because the strain energy is invariant to movements of the dimples.Furthermore, many additional localized solutions exist that are not shown here for brevity.Due to the continuum nature of the cylinder, the prebuckling well of the perfect shell is thus bounded by a large set of localized edge states, each with an associated energy barrier that is eroded most easily by a specific affine imperfection or traversed by a stimulating perturbation.
These edge states do not exist for all levels of compression, as indicated by the limit points that bound the postbuckling paths in Fig. 1(a).The threshold for forming localizations arises due to the varying resistance of the cylinder to concentrated lateral perturbations (poking), as described in the stability landscape of Virot et al. [14].For low levels of global compression, the cylinder resists lateral poking with exclusively positive stiffness [path 1 in Fig. 1(b)].For greater levels of compression, the stiffness to lateral indentation reduces until limit points are traversed that lead to a region of negative stiffness [paths 2-3 in Fig. 1(b)].These two cases occur below the compressive threshold for localizations, i.e., in the region of Fig. 1(b) where poking force vs. indentation curves do not pass through the zero-force axis.For sufficiently large compression, the poking force passes through zero, and at this point, an edge state in the shape of a single dimple is encountered [path 4 in Fig. 1(b)].
Hence, uniaxial compression drives the cylinder's proclivity to form localizations by modulating the resistance to lateral indentation.Here, to model how the cylinder transitions out of the prebuckling well in the presence of imperfections, we analyze a simpler system with analogous lateral perturbation characteristics.The simplicity of the model allows us to rapidly experiment with a large set of random and dominant "engineered" imperfections and provides easier physical insight into the underlying physics.

III. REPRESENTATIVE MODEL
Due to its proclivity to localize, the buckling of a beam on a nonlinear elastic foundation has long been used as a model to investigate shell buckling (see, e.g., Ref. [10]).Here we choose an even simpler mechanical system comprising a series of pin-jointed links supported by nonlinear springs and compressed by displacement-controlled loading at one end [see Fig. 2(a)].This system has the unique advantage that the stretching and bending energies-representing the

FIG. 2. (a)
A schematic representation of a link system on an elastic foundation.Depending on the properties of the elastic foundation, the buckling response can either be (b) localized or (c) distributed.
equivalent contributions in the cylinder-are decoupled and separated in the links and springs, respectively.In addition, the system is readily translated into a small strain, large deformation finite-element model allowing for rapid analysis.Localizations typically occur for subcritical (unstable) bifurcations, whereas distributed (periodic) buckling dominates for supercritical (stable) bifurcations [15].In the link system, this transition from supercritical behavior [distributed buckling, see Fig. 2(c)] to subcritical behavior [localized buckling, see Fig. 2(b)] is readily controlled by tuning the elastic foundation from stiffening to softening.Although the link system has been used previously to explore spatial chaos [12], the interaction of spatial chaos and imperfection sensitivity, and the implications for shell buckling are addressed here for the first time.
As shown in Fig. 2(a), we consider a system of N horizontal links of length L 0 , supported at the joints by vertical, nonlinear springs of stiffness k nl and resting length 0 .The system is compressed by an actuating displacement, u a , at one end and constrained at the other end.The links are not rigid and can compress axially, but no bending moment can develop at the nodes or within the links.Links and springs are modelled as geometrically nonlinear, total Lagrangian finite elements (see Ref. [16], pp.64-70).The links are assumed to be linear elastic and isotropic with Young's modulus E = 5.56 GPa and cross-sectional area A 0 = 2.47 × 10 −4 m 2 .The spring stiffnesses are cubic and given by (1) In the following, we assume a link system comprising 20 links of total horizontal length NL 0 = 1 m.The linear and quadratic spring stiffness are constant at k 1 = 1 kN/m and k 2 = 0 N/m 2 , whereas the cubic stiffness k 3 is varied to modulate the behavior from stiffening (positive k 3 ) to softening (negative k 3 ).The finite-element model is solved using numerical continuation [17] that allows path following of equilibrium manifolds, pinpointing of bifurcations, and branch switching.
The perfect system first bifurcates (buckles) at an applied force of [18] For the chosen parameters, this gives P cl = 12.577 N. Accordingly, all results are presented in terms of the nondimensionalized quantities P/P cl and u a EA 0 /k 1 L 2 0 .

IV. RESULTS
We start with the mechanical response of the geometrically perfect (flat) link system for both a stiffening (k 3 = 100 GN/m 3 ) and softening (k 3 = −100 GN/m 3 ) elastic foundation.The imperfect behavior of the softening system for a number of random and engineered imperfections is discussed afterward.

A. Perfect response
The equilibrium manifolds of the stiffening and softening systems are shown in Figs.3(a) and 3(b), respectively.In both cases, the link system compresses axially in a flat state before losing stability at P/P cl = 1.
The system on a stiffening foundation bifurcates onto another stable equilibrium path corresponding to a distributed, zig-zagging buckling mode [mode 1 in Fig. 3(a)].The unstable bifurcation branches emanating from the second and third critical points are also shown for completeness and show higher-order buckling modes with the zig-zagging distribution interrupted by adjacent collinear segments [modes 2 and 3 in Fig. 3(a)].
In the case of the softening system, k 3 is tuned to produce a force-displacement response that mimics path 4 in Fig. 1(b), i.e., the cylinder lateral response to poking beyond the compressive threshold for edge states.The bifurcated paths now change from sloping upward-right (stable) to downward-left (unstable).The equilibrium paths branching from the first three critical points are again plotted [see Fig. 3(b)] and run parallel to the stable prebuckling path.Even though the eigenvectors associated with these three critical points are periodic and match mode shapes 1 -3 in Fig. 3(a), the specific nonlinearity of the system causes the deformation modes to localize immediately after bifurcating, as shown by modes iiii in Fig. 3(b).The localized nature of mode i means that a single buckle can be induced at any node to create an unstable equilibrium state.For example, the densely spaced unstable postbuckling curves in inset A of Fig. 3(b) (in gray) correspond to the single buckle of mode i translated to any of the other nodes [modes n in Fig. 3(b)].Notably, these additional paths are broken away from the prebuckling path and form limit points as indicated by the dots in inset A. Many additional unstable equilibrium paths also exist [not shown in Fig. 3(b)] that can be traced by forcing any node or combination of nodes laterally into an unstable edge state (until poking force reduces to zero).Each of these unstable edge states defines a possible escape route from the stable prebuckling well and different types of initial imperfections will trigger buckling via different trajectories and modes.The similarities between Fig. 1(a) (cylinder) and Fig. 3 (link system) are worth noting.For low levels of axial compression, the cylinder resists lateral poking with exclusively positive stiffness and localizations cannot form [see path 1 in Fig. 1(b)].Beyond a critical compressive threshold, the resistance to lateral indentation softens sufficiently to allow localizations to form as unstable equilibria [see path 4 in Fig. 1(b)].Hence, axial compression of the cylinder modulates the transition between the two qualitative behaviors described by Figs.3(a) and 3(b).Additionally, in the finitedimensional link system, the number of possible localized edge states scales exponentially with the number of links.In the limit, the link system is thus able to represent the infinitedimensional cylinder, featuring an infinite number of possible localizations.

B. Imperfection sensitivity
To assess how the softening system buckles in a nonperfect setting, we impose sets of random and engineered geometric imperfections.Due to the subcritical behavior, each imperfection rounds off the perfect bifurcation into a nonlinear prebuckling path that destabilizes at a limit point, subsequently leading to localization.Instabilities thus occur at lower loads than the perfect system with the knockdown of each imperfection proportional to the imperfection amplitude.
In the first instance, random geometric imperfections are imposed by defining a nonflat stress-free starting condition.A random vertical displacement in the range of δ y ∈ [−δ mag L 0 , δ mag L 0 ] (with δ mag = 5 × 10 −4 ) is applied to each free node.Random distributions are sampled from MATLAB's rand function and 10 4 different imperfection signatures are analysed.Stochastic variations in the foundation stiffness were also analyzed and produce qualitatively analogous results.
Figure 4(a) shows the normalized buckling loads for the set of random imperfections vs. the axial position (x/NL 0 ∈ [0.05, 0.95]) where individual localizations occur.Even though the mode shape on the unstable path immediately after the critical limit point often features several localized buckles, the unstable buckling mode quickly localizes completely at a single node, and this is the parameter plotted on the x axis in Fig. 4. The results clearly show a large spread in buckling loads.Similarly to other imperfection sensitivity studies on subcritical buckling [19], the P/P cl scatter is well described by a "generalized extreme value distribution"-here with parameters μ = 0.641 (location), σ = 0.0450 (scale), and ξ = −0.142(shape).The distribution has a median normalized buckling load of 0.657, a mode of 0.537, a standard deviation SD) of 0.0497 and a range of 0.308.
The large spread in buckling loads (range = 0.308) occurs even though the maximum random imperfection amplitude is constant and the imperfection signatures are seemingly indistinguishable.Furthermore, the link system buckles with equal likelihood at any of the nodes (apart from the nodes immediately adjacent to the boundaries).Interestingly, the buckling location is not determined by the maximum imperfection amplitude but generally by the location of the sharpest imperfection feature (greatest imperfection curvature).However, even for imperfections with the same local features (causing localization in the same place), the details of the imperfection signature over the rest of the structure strongly influence the nonlinear equilibrium path, thereby producing the observed spread in buckling loads.In a manner reminiscent of temporal chaos, the precise nature of the initial conditions thus defines the trajectory to buckling and thus buckling load and mode.
The second set of imperfections consists in local dimples of the same maximum amplitude δ y = δ mag L 0 as for the random imperfections.Due to the symmetry of the problem about the midspan, these "engineered" single-dimple imperfections are applied, in turn, only to the 10 nodes over half the system (x/NL 0 ∈ [0.05, 0.5]), while the remainder of the nodes are kept in the original (flat) state.This analysis leads to 10 buckling loads, one for each of the 10 single-dimple imperfections.
As expected, the link system always buckles at the location of the single-dimple imperfection and the dashed line in Fig. 4(b) shows the critical load for each of these 10 cases.As secondary random imperfections spanning the entire domain are introduced in addition to each single dimple, the buckling loads disperse around the dashed line; again, following a generalized extreme value distribution.To produce this behavior, we impose secondary random imperfections an order of magnitude smaller than the single dimples (i.e., δ y ∈ [−0.1δ mag L 0 , 0.1δ mag L 0 ]).For each of the 10 dimple imperfections, a total of 10 3 random imperfections are run.Away from the boundaries the median buckling load is 0.749, the mode 0.707, the SD 0.0164, and the range 0.0894.If the analysis is repeated for single-dimple imperfections and random imperfections scaled proportionally by a factor of 10 (dimple: δ y = 10δ mag L 0 ; random variation: δ y ∈ [−δ mag L 0 , δ mag L 0 ]), then the median buckling load decreases to 0.231 but the spread of buckling loads remains of similar magnitude (SD 0.0145 and range 0.0777).However, the closer the amplitude of random imperfections to single-dimple imperfections, the greater the stochasticity in buckling loads.
In conclusion, although imperfections of the same amplitude are imposed in Figs.4(a) and 4(b), dominant engineered imperfections lead to significantly less stochasticity than entirely random imperfections.The imposition of a dominant imperfection thus reduces the multiplicity of possible escape trajectories in spatially chaotic systems and improves the ability to make accurate predictions of the buckling load.

V. IMPLICATIONS FOR SHELL BUCKLING
The reduction in stochasticity due to a dominant imperfection accounts for the accurately predicted buckling loads of shells with precisely engineered imperfections [20] and the ability to control their postbuckling response [21].The dominant imperfection creates a localized prebuckling stress field that biases the shell to follow a specific equilibrium trajectory toward buckling.In effect, the dominant imperfection breaks the prebuckling symmetry, forcing the ensuing localization to form in one spatial region.This insight suggests possible avenues for designing less imperfection-sensitive shells.
Due to the imperfection insensitivity of supercritical systems, it is no surprise that stochasticity in the link system can be reduced by stiffening the elastic foundation.By changing the nonlinear stiffness component from k 3 = −100 GN/m 3 to k 3 = −1 GN/m 3 , the median buckling load, SD, and range over the set of 10 4 random imperfections change from 0.657, 0.0497, and 0.308 to 0.830, 0.0138, and 0.196, respectively.Hence, both the knockdown and spread of buckling loads due to random imperfections are reduced.Circumferential and axial stiffeners are often used in aerospace structures to prevent buckling of thin-walled shells.These stiffeners help to break the shell into smaller unsupported segments to increase buckling loads, but the stiffening effect also helps to reduce the observed imperfection sensitivity.
A second avenue toward reducing imperfection sensitivity is to induce nonuniform prebuckling stresses via elastic tailoring.To simulate a modulation in prebuckling stresses, a quadratic preforce distribution of F 0 = 0.08(x − 0.5) 2 N was applied to the springs in the link system and the imperfection sensitivity analysis was repeated over 10 4 random imperfections.In this case, buckling always occurred at the two most preloaded nodes (x/NL 0 = 0.05 and x/NL 0 = 0.95) with increased median buckling load of 0.740 and reduced SD and range of 0.0283 and 0.172, respectively.
With modern manufacturing technology, such blended means of elastic tailoring and controlled redistributions of prebuckling stresses are possible, e.g., via variable-angletow composites that steer fiber paths in curved trajectories across shells [22].The possibility of designing imperfectioninsensitive shells via elastic tailoring thus presents a promising avenue for future research.
Data are available at the University of Bristol data repository [23].

FIG. 1 .
FIG. 1. (a)Equilibrium paths of applied end-shortening u (normalized by cylinder radius R, length L, and wall thickness t) vs. the reaction load P (normalized by the linear stability prediction P cl[13]).The stable prebuckling path is bounded by a set of localized unstable edge states.(b) Resistance of the cylinder to concentrated lateral probing (poking force F vs. indentation displacement w) for different levels of axial compression u.Plots adapted from Ref.[8].

FIG. 4 .
FIG.4.A set of 10 4 buckling loads vs. the location along the link system where localized buckling occurs.(a) Results for a set of random imperfections and (b) results for single-dimple imperfections with additional superposed random imperfections.The dashed line shows the buckling loads for single-dimple imperfections without additional random imperfections.
Equilibrium diagrams of end-shortening u vs. reaction force P for (a) a stiffening foundation leading to supercritical behavior and distributed buckling and (b) a softening foundation leading to subcritical behavior and localized buckling.