Orthotropy as a driver for complex stability phenomena in cylindrical shell structures. Composite Structures

By applying a quasi-static transverse load, shallow cylindrical shells can be snapped between two inverted stable states. This dynamic transition, initiated at either a limit point (fold) or a subcritical bifurcation, traverses a region of instability. The ensuing large displacements are attractive for shape adaptation of multi-functional engineering structures. Previous research on isotropic cylindrical shells has indicated the presence of isolated regions of stability in the otherwise unstable transition region. While these regions are theoretically attractive for multi-stage snap-morphing, they are di ﬃ cult to attain in practice by means of a single control parameter. In this paper, we study the e ﬀ ect of orthotropic material properties on the structural stability and snap-through be- havior of shallow cylindrical shells. In addition, we explore complex stability phenomena created by material orthotropy, which are attractive for multi-stage morphing. The problem is analyzed in a robust manner using a technique known as generalized path-following, which combines the mathematical domains of ﬁ nite element analysis and numerical continuation. The present study shows, in particular, that laminates comprising layers of di ﬀ erent directional material properties provide the ability to tailor the elastic bifurcation behavior of cylindrical shells. For example, the lamination scheme can be varied to add isolated regions of stability or entirely remove them; to induce snaking that allows transitioning between the two inverted cylindrical shapes through a series of snaps; and ﬁ nally, to create groupings of multiple unique stable con ﬁ gurations that can be attained via addi- tional shape control.


Introduction
In engineering design, structural instabilities are generally regarded as unwanted sources of failure rather than beneficial design features. For example, structural optimization often calls for a concurrency of failure modes, but this approach can make the structure imperfection sensitive or lead to post-critical collapse [1,2]. Therefore, the traditional philosophy of designing engineering structures is to prevent buckling or at least make its effects benign, e.g. by guaranteeing stable post-critical behavior. On the other hand, when additional functionality is to be embedded within structures-e.g. to allow for large elastic shape reconfigurations [3], to create meta-materials with innovative properties [4], or to design self-encapsulating devices [5]-instabilities are often encouraged.
Here, we focus on shape-morphing, which is a biomimetic design strategy that aims to adapt structures to conform to different loading conditions. Particularly in the aerospace industry, shape-morphing structures are viewed as a promising technology to enable more structurally efficient designs (see e.g. [6][7][8][9][10][11][12][13]). Strictly speaking, high-lift devices on aircraft wings, such as leading-edge slats and trailingedge flaps, are examples of this technology, but current solutions typically rely on heavy hydraulic or electric actuators to move rigid components through large displacements. A multi-stable structure can, in theory, be used to the same effect, as snapping from one stable shape to another creates the large displacements necessary for meaningful shape adaptation. The simplest example of a multi-stable structure is a post-buckled Euler strut, which can be snapped from one sinusoidal post-buckled shape to another by means of a transverse actuation load (see Fig. 1). An additional advantage of these multi-stable snap-through devices is that sensing, actuation and control may be embedded within the nonlinear mechanics of the structure, and this attractive passive control strategy has recently been demonstrated successfully for fluid flow control [14,15].
Multi-layered, fiber-reinforced composite materials have played a predominant role in the research on shape-morphing because their orthotropic material properties and laminated nature facilitates a union of materials with dissimilar directional material properties. The underlying physical principle is the same as that of commonplace thermally https://doi.org/10.1016/j.compstruct.2018.05.013 Received 11 February 2018; Accepted 2 May 2018 actuated bimetallic strips or piezoelectric bimorphs used for actuation, sensing and energy harvesting applications. For example, when two layers of fiber-reinforced plastic with reinforcing fibers running in perpendicular directions are assembled into a laminate, and then cured at high temperature and cooled, the differences in thermal expansion between layers cause warping into curved configurations [16]. This manufacturing-induced warping can then be exploited to snap the laminate between two inverted self-equilibrated shapes. As such, the orthotropic nature of laminated composites allows the creation of bistable [17][18][19] or even tristable shell structures [20,21].
Of course, the same bistable behavior can also be obtained from an isotropic shallow cylindrical shell, which can be snapped from one cylindrical shape to its geometric inverse. Furthermore, the onset of snapping in an isotropic cylindrical shell is not as straightforward as simple loss of stability at a limit point. An isotropic cylindrical shell subjected to a quasi-static transverse dead load contains a multitude of symmetry-breaking bifurcation branches within the unstable transition region between the two inverted shapes. Although most of these additional branches are statically unstable, a recent study has shown that an isolated region of stability exists between the two geometrically inverted cylindrical states [22]. While this region is theoretically attractive for creating morphing devices capable of multi-stage snapping, it is impossible to attain this additional stable configuration by means of a single control parameter from the unloaded state.
The aim of this work is to investigate additional bifurcations and isolated regions of stability that arise when orthotropy is introduced in cylindrical shell structures. In particular, we show that certain lamination schemes lead to snaking [23,24] that enables multi-stage morphing.
The remainder of the paper is structured as follows. Section 2 briefly summarizes the chosen modeling framework, which is based on a combination of finite element discretization and numerical continuation algorithms. Section 3 outlines the problem definition and the details of the isotropic shell behavior. Section 4 then studies the bifurcation behavior of the orthotropic shell and highlights the effect of lamination sequence on the isolated regions of stability (Section 4.1); the occurrence of snaking (Section 4.2); and the grouping of multiple unique stable configurations (Section 4.3). Finally, a summary of the work and conclusions are presented in Section 5.

Modeling framework
Due to its geometric versatility the finite element method is the preferred technique for solving complex problems arising in engineering mechanics. In the applied mathematics community, the mathematical methods of multi-parameter analysis, branch-switching and bifurcation tracking implemented in numerical continuation codes are well established [25][26][27]. However, these are not classically used for structural mechanics (statics) applications, where specialized arclength techniques restricted to a single parameter-the applied load-are predominantly used [28,29]. The application of broader numerical continuation methods within structural finite element solvers often go by the name of "generalized path-following" [30][31][32], and it is such a framework that we use herein to analyze the bifurcation behavior of cylindrical shell structures.
The present formulation considers a discretized model of a slowly evolving, conservative and elastic structure, where the internal forces, f u ( ), and tangential stiffness, K u ( ) T , are uniquely defined from the current displacements, u, by means of the first and second variations of the total potential energy. Thus, non-conservative loading and historydependent problems such as plasticity are not included in the current formulation.
An accurate stability analysis using these numerical methods depends on suitable choices of nonlinear finite elements that have sufficient fidelity in their kinematic assumptions and strain definition to accurately represent the higher-order energy terms of a Taylor series expansion of the total potential energy [33]. Throughout this work we use so-called "degenerated shell elements" [34] based on the assumptions of first-order shear deformation theory [35] formulated using the full Green-Lagrange strain tensor and a total Lagrangian reference system. Hence, the analysis allows for finite displacements and rotations but is restricted to the small strains implicit of a classic Hookean constitutive law.
In statics, an equilibrium state is expressed as a balance between internal and external forces, where, in a displacement-based finite element setting, this balance is written in terms of n discrete displacement degrees-of-freedom, u, and a scalar loading parameter, λ, The vectors p λ ( ) and f u ( ) are the external (non-follower) load and internal force vector, respectively. In the case of linear and proportional loading we have≡ , wherep is a constant reference loading vector (dead loading) and the subscript comma notation is used to denote partial differentiation.
For generalized path-following, Eq. (1) is adapted to incorporate any number of additional parameters, such that, is a vector containing p control variables, with Λ 1 corresponding to parameters that influence the internal forces (e.g. material properties, geometric dimensions, temperature and moisture fields) and Λ 2 relating to externally applied mechanical loads (e.g. forces, moments, tractions).
The n equilibrium equations in (2), correspond directly to the n displacement degrees of freedom in the system. Because the structural response is parametrized by p additional parameters, a p-dimensional solution manifold in    + ( ) exists. Specific solution subsets on the pdimensional manifold are defined by incorporating additional auxiliary equations, g. Hence, we wish to evaluate solutions to the augmented system For r auxiliary equations, the solution to Eq. (3) becomes − p r ( )-dimensional, and hence − p 1 auxiliary equations are required to define a one-dimensional curve-a so-called subset curve of the multi-dimensional solution manifold. These subset equations can define fundamental equilibrium paths (the fundamental load parameter is varied); parametric equilibrium paths (a non-load parameter is varied); branchconnecting paths that determine the total number of additional branches emanating from a bifurcation point; paths that track bifurcations in parameter space; etc. To track bifurcations, we need to simultaneously enforce the fulfillment of a criticality condition, e.g. . In the most general form, not limited to but including the previous criticality condition, a vector of q auxiliary variables, v, may be added to the auxiliary equations g, .
Hence, Eq. (4) describes n equilibrium equations and r auxiliary equations in + + n p q ( ) unknowns leading to a + − p q r ( ) -dimensional solution. To determine a one-dimensional subset curve of singular points, we thus require = + − r p q 1 auxiliary equations to constrain the system. Following the criticality example for bifurcation tracking referenced above, when the n-dimensional null vector at the critical state is introduced as the auxiliary variable, v, a singular subset curve in two parameters, = p 2, is appropriately constrained by the associated = + r n 1 auxiliary equations = K v 0 When evaluating one-dimensional curves ( = + − r p q 1), one additional equation is needed to uniquely constrain the system to a solution point = y u v Λ ( , , ). Hence, , N where N is a scalar equation which plays the role of a multi-dimensional arc-length constraint. A specific solution to Eq. (5) is determined by a consistent linearization coupled with a Newton-Raphson algorithm, where the superscript denotes the jth equilibrium iteration and the subscript denotes the kth load increment. The iterative correction cycle is typically started by a predictive forward Euler step. For most problems, the inversion of the iteration matrix, is significantly simplified by partitioning the system into blocks such that only the symmetric tangential stiffness matrix, K T , needs to be factorized in the solution process (see e.g. [36]). Within this framework, the meaning of the term generalized pathfollowing becomes clear. It refers to the fact that any arbitrary curve can be traced on the equilibrium manifold, as long as a pertinent auxiliary equation is defined that constrains the equilibrium equation to the locus of points required. For implementation-specific details on the modeling framework and finite elements in general, the interested reader is directed to the papers by Eriksson and co-workers [31,37,38].

Islands of stability in isotropic shells
The hinged, thin, isotropic, shallow cylindrical shell section drawn schematically in Fig. 2 is a classical benchmark problem in the nonlinear, large deformation finite element literature [39]. The radius of the shell is = R 2540 mm with hinged longitudinal edges of length = L 508 mm and free circumferential edges of equal arc-length = B 508 mm. The thickness, Young's modulus and Poisson's ratio of the shell are = = t E 6.35 mm, 3102.75 MPa and = υ 0.3, respectively. A transverse point load P c (dead load) is applied at the center C of the roof plan-form and acts as the fundamental loading parameter. Equally, under a rigid loading regime the dead load P c may be replaced by a controlling displacement parameter w c . The nature of the fundamental loading parameter (dead vs rigid loading) leaves the equilibrium paths unchanged-i.e. the roots of the first variation of the energy potential-but alters the stability of these equilibria-i.e. the sign of the determinant of the second variation of the energy potential.
For many years the full bifurcation behavior of this transversely loaded shell structure remained elusive. While the equilibrium path of symmetric deformation modes (fundamental equilibrium path) was successfully determined in its entirety, including the existence of two load limit points on this fundamental equilibrium path, a symmetrybreaking pitchfork bifurcation before the first limit point was often ignored (see the two modes in Fig. 2). After this asymmetric equilibrium path was established in the literature by means of initial imperfections applied to the structure or asymmetric meshing techniques [40], two further pairs of bifurcation points and connecting secondary paths were recently determined [22]. This latter study showed that a small interval of one of these secondary equilibrium paths is stable, whereas all other portions of the secondary paths are unstable. Indeed, as shown later in this section, there exist three further tertiary equilibrium paths, branching from the recently determined secondary paths, that are unstable under dead loading (load control), one of which can be stabilized under rigid loading (displacement control). Throughout the analyses in Section 3 and 4 the shell is discretized into a × 31 31 node mesh of 100 fully integrated 16-node, total Lagrangian shell elements using the shell director parametrization of Ramm [41] (two rotations per node, no drilling around the shell director). The effects of shear and membrane locking are minimized by using bi-cubic isoparametric interpolation functions (16-node elements) and by refining the mesh sufficiently until convergence with respect to a high-fidelity × 68 68 element mesh of reduced integration S4R elements of the commercial finite element software ABAQUS is obtained. Fig. 3a shows the full equilibrium behavior of the isotropic shell under dead loading, plotted in terms of the transverse displacement w c versus the applied transverse load P c at the center C of the roof. Blue segments denote stable equilibrium paths (all eigenvalues of the tangential stiffness matrix are positive), whereas red segments denote unstable equilibrium solutions (at least one negative eigenvalue). Black points correspond to critical points, i.e. equilibrium solutions where at least one eigenvalue of the tangential stiffness matrix is exactly zero (to a numerical tolerance on the equilibrium Eq. (2) of − 10 ) 6 . These critical points are either symmetry-preserving limit points or symmetrybreaking bifurcation points, where at least one secondary branch intersects the primary branch.
The equilibrium path beginning at the origin (path 1) in Fig. 3a is the typical fundamental path reported in numerous studies in the literature and corresponds to the symmetric mode shape shown in Fig. 2. This symmetric mode becomes unstable at the first bifurcation point, where another equilibrium path 2 (corresponding to the asymmetric mode in Fig. 2) branches from the fundamental path 1. Because this first pitchfork bifurcation is subcritical, in an experimental dead loading scenario the shell would deform symmetrically along fundamental path 1 and then snap across to the inverted shape once the first bifurcation point is reached. The fundamental path 1 and bifurcation branch 2 traced by ABAQUS' arc-length solver using S4R shell elements match closely with the present model. In total, there are three secondary bifurcation branches emanating from the fundamental path (paths 2-4), which are mostly unstable under the imposed dead loading regime, apart from the isolated region of stability on bifurcation branch 3 as previously established by Zhou et al. [22]. Identifying this region experimentally is not a trivial task, as this region cannot be reached through a continuous stable path from the unloaded state. Rather, a load of around 200 N could be applied from the unloaded state, and the roof then perturbed into the mode shape corresponding to the isolated region of stability. Additional tertiary branches (paths 5-7), to the authors' knowledge shown here for the first time, emanate from the three secondary branches (paths 2-4), which either connect back to the secondary branches from which they emanated (paths 5 and 6), or connect to other secondary branches, such as tertiary branch 7 connecting branch 4 to branch 3.
Although these additional tertiary equilibrium paths 5-7 are all unstable under dead loading, an interval of path 5 can be stabilized by switching to rigid loading. As shown in Fig. 3b, under rigid loading the equilibrium solutions remain unchanged but their stability changes-the load limit points in Fig. 3a are replaced by displacement limit points in Fig. 3b (symmetry-breaking bifurcation points do not change), and previously unstable equilibrium paths are now partly stable. Fig. 3b shows that the asymmetric bifurcation branch 2 is stable under rigid loading, such that the shell does not snap into its inverted state upon reaching the first critical load but smoothly transitions between the two inverted cylindrical shapes via the asymmetric mode. Furthermore, the isolated region of stability on bifurcation branch 3 covers a larger range of the control parameter, such that a rigid loading regime makes it easier to verify this isolated region of stability experimentally. For example, the shell could be loaded to = w 15 mm c , or close to the region where the asymmetric path 2 and the isolated region on path 3 intersect, and the shape then perturbed from the path 2 mode (shown in Fig. 3b) to the path 3 mode (shown in Fig. 3a) via additional control points. In addition, inset A in Fig. 3b shows that a small portion of the asymmetric bifurcation branch 2 temporarily loses stability causing the structure to snap to the adjacent equilibrium path 5. The mode shape of path 5 maintains the asymmetry in the curved direction of path 2, but adds an orthogonal asymmetry in the uncurved direction of the cylindrical shell. Hence, the equilibrium diagram in Fig. 3b suggests that, in a rigid loading experiment, the shell would briefly snap into and out of the path 5 mode shape while being loaded along path 2.
Ideally, the isolated region of stability on path 3 could be exploited for multi-snap morphing between the two cylindrical shapes via an intermediate stage. Because this region cannot be reached through a continuous stable path from the unloaded state, this is not possible unless additional points of control, such as actuating piezoelectric patches, are attached to the shell. While this would still be a useful solution for a potential engineering application, the isolated region may be more readily induced as certain design parameters are changed. For example, if for some value of shell thickness the isolated stable region intersected the zero load axis, then the associated mode shape would be self-equilibrated and could be observed in an unloaded state, thereby converting the shell from bi-to tristable. This investigation can be performed parametrically by running the full solution shown in Fig. 3 for a number of models with different thickness. Given that the stable region is bounded by two limit points, a computationally more efficient method is to trace a foldline of these two critical points with respect to variations in the shell thickness via bifurcation tracking.
The circular black curve in Fig. 4 shows the locus of two limit points, which delimits the boundary of the stable region for changes in shell thickness. The direction of decreasing thickness is indicated by the direction of the arrows with the −t label. As an equilibrium manifold plotted in w vs P vs t c c space, the stable region interior to the black foldlines visually represents an island of stability surrounded by unstable equilibrium solutions. This region is bounded at = t 7.12 mm, where the stable region ceases to exist at a cusp catastrophe, and the limiting case of zero thickness. On the left-hand side, the foldline approaches a locus of pitchfork bifurcation points, showing that the bifurcation point and minimum limit point asymptotically converge as the thickness of the shell decreases. The foldline also suggests that the length of the stable portion of path 3 can be maximized for a thickness of ≈ t 4.75 mm. Similar analyses were also conducted for varying length, width and rise of the cylindrical panel but in all cases it was not possible to force the stable region to intersect the zero load axis. This suggests that it is not possible to convert this structure into a selfequilibrated tristable shell.

Othotropy as a driver for more complex stability phenomena
We now turn to the effect of changing the material system of the shell from isotropic to orthotropic. The material properties are highlighted in Table 1 and correspond to a fiber-reinforced plastic with different directional material properties along and perpendicular to the reinforcing fibers.
Multiple orthotropic layers of equal thickness with different fiber orientations are laminated into a layered structure. Most of the laminates considered herein are quasi-isotropic, meaning that the effective axial modulus of the laminate is the same in all directions. One way to achieve this is by guaranteeing equal amounts of°°0 ,90 and ±°45 layers in the laminate, where the°0 direction is defined along the uncurved ydirection in Fig. 2, and the°90 direction corresponds to the curved xdirection. Here we consider quasi-isotropic, symmetric and balanced laminates with eight layers, e.g. − [0/90/45/ 45] s , where subscript s denotes a laminate that is symmetric about its mid-plane. This means that rearranging the order of the individual layers does not affect the effective axial rigidity of the shell. It does, however, affect the directional bending rigidity because the bending rigidity of individual layers increases with distance from the laminate mid-plane. The material properties in Table 1 have been chosen so that the effective axial modulus of the quasi-isotropic shell-constant in all directions-is equal to the isotropic Young's modulus (3102.75 MPa) considered in the previous section.

Stabilizing and destabilizing effects
By using quasi-isotropic laminates with the same effective modulus as the isotropic shell discussed in Section 3, we can investigate the effect of bending rigidity on the snap-through behavior of the shell while controlling for axial rigidity (total thickness remains 6.35 mm). A comparison of six quasi-isotropic laminates is shown in Fig. 5 For each of the equilibrium path diagrams in Fig. 5 (dead loading) only the branches emanating from the first two bifurcation points occurring on the fundamental path are computed, as all additional branches are unstable and only unnecessarily complicate the figures. Even though comparing the different lamination schemes is not straightforward due to the complexity of the stable and unstable equilibria, some general observations can be made: • Maximizing the bending stiffness in the uncurved direction and minimizing it in the curved direction (Fig. 5a − [0/45/ 45/90] s ) removes the isolated stable equilibrium in the unstable transition region, and decreases the snap-through load compared to the equivalent isotropic panel. By increasing the bending rigidity in the curved direction slightly, i.e. moving the°90 layer outwards (Fig. 5b − [0/90/45/ 45] s ), the snap-through load is increased and two isolated regions of stability are created in the unstable transition zone. For both cases, the secondary branches are considerably more entangled than for the isotropic shell.
• Maximizing the twisting rigidity removes the isolated stable regions in the unstable transition zone and significantly simplifies the nature of the secondary equilibrium branches. In this case, i.e. with the ± 45 layers outermost, increasing the bending rigidity in the curved direction by placing the°90 layer further from the mid-plane ( • Maximizing the bending rigidity in the curved direction ( Fig. 5f − [90/45/ 45/0] s ) significantly increases the snap-through load compared to the isotropic shell by stabilizing the fundamental path. Indeed, in this case the first bifurcation point on the fundamental path occurs after the first limit point, such that the shell snapsthrough upon reaching a limit point but snaps back in the reverse direction when reaching a bifurcation point. By decreasing the bending rigidity in the curved direction slightly (Fig. 5e − [90/0/45/ 45] s ), the first bifurcation point again occurs before the limit point.
To summarize, we can conclude that increasing the bending rigidity in the uncurved direction destabilizes the shell and leads to more entangled bifurcation branches, whereas increasing the bending rigidity in the curved direction stabilizes the shell on the fundamental path as well as on the secondary bifurcation branches in the transition region. From fundamental mechanics of one-dimensional arch structures loaded transversely, we know that the snap-through buckling instability is driven by destabilizing compressive stresses that are induced as a reaction to the outwards thrust against the supports. The classical buckling load of a two-hinged, isotropic, circular arch with incompressible center line subjected to a uniform radial pressure is directly proportional to the bending rigidity EI [42]. Because Hurlbrink's equation can also be used for cylindrical shells if EI is replaced with it is no surprise that increasing the bending rigidity in the curved direction stabilizes the shell.

Snaking and multi-snap morphing
The observations in the previous section suggest that placing the°90 layer outermost in the stacking sequence is most effective at creating isolated regions of stability in the unstable transition region. We next consider the − [90/0/45/ 45] s laminate under rigid loading (Fig. 6), but with half the original thickness, i.e. = t 3.175 mm, to increase the extent of the snaking behavior described in this section. In all equilibrium diagrams of Fig. 6 the fundamental path has been curtailed at point A and then restarted at point B to reduce the complexity of the figures.  space the stable region interior to the black foldlines forms an island of stability surrounded by unstable equilibrium solutions.

Table 1
Orthotropic material properties used for the layered shell, where the 1-direction points along the fibers of a reinforced plastic, the 2-direction is perpendicular to the fibers, and the 3-direction is normal to the orthotropic layer.  Fig. 6a shows the first secondary branch connecting the two bifurcation points denoted by P. Note that the first bifurcation point on the fundamental path no longer connects to the final one, as was previously the case for the thicker shell in Figs. 3 and 5, but rather to the second bifurcation point on the fundamental path. Similarly, Fig. 6b shows that another bifurcation branch connects the final bifurcation point on the fundamental path with the third one. An interesting feature of this shell's behavior is that especially the first secondary equilibrium path (Fig. 6a) shows a visual example of snaking behavior, whereby the structure loses stability at a series of limit points to snap through a series of mode shapes where no further lines of symmetry are broken. Hence, upon increasing rigid loading the shell initially transitions into the asymmetric mode I at the first bifurcation point and then experiences two discrete snaps to higher-order mode shapes II and III.
However, due to the existence of the second bifurcation branch (Fig. 6b), the exact snapping behavior throughout the loading history can only be accurately determined via a dynamic analysis. Hence, an implicit dynamic analysis was conducted in the commercial FE solver ABAQUS with a mesh of × 50 50 4-noded, reduced integration S4R elements. The density of the material was assumed to be = ρ 0.0016 g/mm 3 (typical of a fiber-reinforced plastic) and Rayleigh damping coefficients proportional to mass ( = α 0.05) and stiffness ( = β 0.05) were applied to dampen out vibrations after snapping occurs. To simulate a quasi-static loading regime, the snap-through displacements from 0 to 30 mm were applied in 1500 s, and equally the unloading snap-back behavior was also modeled through a time period of 1500 s. Additionally, a small geometric imperfection had to be applied to allow the structure to easily follow the first bifurcation branch off the primary path. The solution from ABAQUS Dynamic for the snap-through and snap-back behavior is shown in Fig. 7a and Fig. 7b, respectively, and the correlation between the dynamic solution from ABAQUS and the stable regions of the static solution from our in-house solver is excellent.
Upon increasing rigid loading from zero displacement (Fig. 7a), the shell initially breaks the symmetry of the fundamental symmetric mode and transitions into the asymmetric mode shape I (for diagrams of all mode shapes referenced hereon, please refer to Fig. 6). At the following limit point the shell then snaps into mode shape II. Surprisingly, the shell does not snap into mode shape III at the next limit point but snaps directly into the S-shaped mode V. When the shell snaps, internally stored strain energy is converted into kinetic energy or dissipated as heat by the imposed damping factors. Thus, the velocity of the shell during snapping is governed by the density of the material (fixed) and the damping applied in the model. It is possible that for increased damping (reduced snapping velocities) the shell may stabilize in mode shape III instead, but because the applied Rayleigh damping coefficients ( = = α β 0.05) are already relatively high, this is unlikely to occur in reality. Upon further loading in mode V the shell then undergoes a small snap into mode IV and a final smooth transition into the full    inverted cylindrical shape. Upon decreasing rigid loading (Fig. 7b), the shell reverses the transition into mode IV and V. When mode V cannot sustain further decreases in displacement, the shell snaps into mode VI and shortly thereafter into mode I. From this point the structure can be smoothly unloaded into its original cylindrical shape.
This multi-snap behavior is interesting from a shape-morphing perspective because it creates the opportunity for multi-stage morphing. Rather than immediately snapping from one cylindrical shape to another, the transition can be designed in a more gradual fashion such that different shapes can be exploited at different stages of the operating environment. One potential application would be a tristable configuration, whereby mode V is used as a baseline configuration from which the inverted modes I and IV are readily induced. Hence, the shell could be designed to snap between these three modes under a specific loading regime to, for example, alter the topological profile.
To show that this multi-snap behavior is driven by orthotropy, rather than by the thickness reduction from = t 6.35 mm to = t 3.175 mm, we have plotted the fundamental primary path and first two secondary bifurcation branches of the isotropic shell with reduced thickness ( = t 3.175 mm) in Fig. 8. The equilibrium paths are very similar to those observed for the thicker isotropic shell in Fig. 3b, whereby the   Fig. 9. Plots of the primary equilibrium path, secondary and tertiary bifurcation branches that form a collection of closely spaced stable regions towards the righthand side of the diagrams. The full behavior is shown in (f) and individual curves are shown in (a)-(e).
secondary branches connect the same bifurcation points and no snaking behavior occurs. Hence, we conclude that the complex stability phenomena shown in Fig. 6 are indeed driven by orthotropic material properties.

Grouping of isolated stable states
As a final case study, we consider a fully orthotropic (not quasiisotropic) shell with the stiffer material axis aligned with the curved direction, and the unreinforced direction along the cylindrical shell axis. Hence, referring to the stacking sequence notation used previously for the quasi-isotropic multi-layered shells, we are analyzing a [90] s laminate comprising one or multiple°90 layers. The overall thickness of the shell is = t 6.75 mm, and a dead loading regime is considered. Fig. 9a shows the fundamental primary path (stable regions in black and unstable regions in gray), and the first bifurcation branch connecting the first and last bifurcation points ( − P P 1 1 ) on the fundamental path. Similarly, Fig. 9b shows the primary path and the second bifurcation branch connecting the second and penultimate bifurcation points − P P 2 2 on the fundamental path. Finally, Fig. 9c shows the primary path and the third bifurcation branch connecting the third and third-last bifurcation points − P P 3 3 on the fundamental path. What is immediately apparent from these three figures is that the four equilibrium paths depicted all group stable equilibria towards the right-hand side of the plots. Fig. 9d shows a tertiary path which connects the secondary path − P P 1 1 (Fig. 9a) with secondary path − P P 3 3 (Fig. 9b) from points P 4 to P 4 . Similarly, Fig. 9e shows another tertiary path which connects the same two secondary paths via a second set of points − P P 5 5 . Again, both these tertiary equilibrium paths feature isolated stable regions on the righthand side of the plot close to the primary stable path. Fig. 9f then aggregates these six equilibrium paths on a single diagram with stable regions shown in blue and unstable regions shown in gray. This figure clearly shows how four different stable equilibrium paths aggregate together towards the right-hand side of the plot. This bunching together is, of course, partly only a feature of the chosen displacement metric w c . If a different displacement on the shell planform was plotted on the horizontal axis, the shape and position of the equilibrium paths would change. Indeed, Fig. 10 shows that the mode shapes of these four stable equilibrium solutions are indeed quite distinct.
Generally speaking, there are a number of striking features about the behavior of this particular shell structure. First, the equilibrium diagram (Fig. 9f) features multiple entangled post-buckling paths (additional fully unstable paths are not depicted for clarity), and this behavior is indicative of spatial chaos [44]. If this shell were to be cyclically loaded to snap from one cylindrical shape into its inverted form, the interactions between stable and unstable equilibria would likely lead to chaotic behavior.
Second, it is not entirely clear into which mode the structure would naturally snap into once the first bifurcation branch loses stability at the limit point. The solution from ABAQUS Dynamic (Fig. 10a) shows that the shell always snaps-through to the fully inverted shape IV and does not stabilize in any of the three intermediate mode shapes I-III (Fig. 10b). It is also worth noting that each of these intermediate mode shapes I-III only exists over a finite range of the fundamental loading parameter delimited by limit and bifurcation points (insets A, B, D and E of Fig. 9). This means that beyond a load > P 1800 N c the only deformation mode shape that is statically stable is the full inverted shape IV.
However, for a load ⩽ ⩽ P 545 N 879 N c (inset F of Fig. 9f) the shell is indeed stable in five different mode shapes-the fundamental mode starting from the unloaded state at the origin, and mode shapes I-IV (Fig. 10b). Thus, this behavior creates another interesting opportunity for shape adaptation. In the range ⩽ ⩽ P 545 N 879 N c the primary structural shape to be exploited may be the first fundamental mode, but by means of additional control points (e.g. piezoelectric patches) distributed over the plan-form of the shell, the structure could be forced to morph into one of four possible shapes at constant load. The advantage here is that, for all four shapes, the associated displacement at the load application point remains the same but the structure may adapt its shape around this point. How, and if, such an elaborate control scheme is indeed possible is a potentially interesting topic for future work.

Conclusions
Shallow cylindrical shells are attractive structures for shape-adaptive morphing structures because their multi-stability is driven by large displacements which lead to appreciable changes in shape. When deformed by a transversely applied point load, isotropic cylindrical shells are often only bistable, and any further isolated regions of stability are difficult to induce using a single load-controlling parameter. We have shown that introducing orthotropy significantly opens the design space for many other interesting stability phenomena. The case studies considered here reveal the surprising complexity of stable and unstable equilibrium solutions that can be induced simply by varying the elastic moduli along the two principle axes of curvature of a cylindrical shell.
For example, by varying the bending rigidity (at constant axial rigidity) using different lamination sequences of orthotropic layers, the behavior of the shell can either be stabilized or further destabilized. Particularly by stiffening the bending response in the curved direction, additional stable regions can be added to the otherwise unstable  Fig. 10. The dynamic snapping behavior of the shell analyzed quasi-statically in Fig. 9. The dynamic analysis in (a) shows that the shell does not stabilize in any of the closely spaced modes I-III, shown in (b), but snaps straight into the fully inverted mode IV. transition region between two inverted cylindrical shapes. This phenomena arises as a result of snaking, whereby the structure loses stability at a series of limit points, thereby sequentially snapping through a series of higher-order mode shapes. Finally we have shown that in extreme cases, the stability behavior of the shell can become incredibly complex, featuring a plethora of entangled stable and unstable postbuckling paths that are indicative of spatial chaos. For example, for a cylindrical shell comprising orthotropic fiber-reinforced layers aligned purely with the curved direction, four structural stable shapes can be achieved for almost the same combination of centrally applied load and resulting displacement at the control point. This means that if additional control points are introduced, the cylindrical shell can theoretically be snapped into one of four possible shapes depending on which control points are actuated.

Data statement
Data are available at the University of Bristol data repository, data.bris, at https://doi.org/10.5523/bris.j48ferikjn4k205v9o38zt42z.