## Abstract

Multistable shells are thin-walled structures that have more than one stable state of self-stress. We consider isotropic axisymmetrical shallow shells of arbitrary polynomial shapes using a Föppl–von Kármán analytical model. By employing a Rayleigh–Ritz approach, we identify stable shapes from local minima in the strain energy formulation, and we formally characterize the level of influence of the boundary conditions on the critical geometry for achieving bistable inversion—an effect not directly answered in the literature. Systematic insight is afforded by connecting the boundary to ground through sets of extensional and rotational linear springs. For typical cap-like shells, it is shown that bistability is generally enhanced when the extensional spring stiffness increases and when the rotational spring stiffness decreases, i.e. when boundary movements in-plane are resisted but when their rotations are not; however, for certain other shapes and large in-plane stiffness values, bistability can be enhanced by resisting but not entirely preventing edge rotations. Our predictions are furnished as detailed regime maps of the critical geometry, which are accurately correlated against finite-element analysis. Furthermore, the suitabilities of single degree-of-freedom models, for which solutions are achieved in closed form, are evaluated and compared to our more accurate predictions.

## 1. Introduction

Multistability in shells arises from the interaction of out-of-plane bending and in-plane stretching. A bistable example from nature is the Venus Flytrap, which can rapidly ‘snap-through’ and close from an open equilibrium position to ensnare prey [1]; artificial applications range from familiar objects such as flick bracelets and thermostatic switches in kettles [2] to semiconductors in the field of microelectronics [3], and to recent innovations in energy-harvesting devices [4] or deployable spacecraft elements [5].

Fundamental research on bistable structures was first undertaken by Timoshenko [6] for curved beams, and by Wittrick *et al.* [2] as well as Hyer [7] for shells. As multistability, in general, is concomitant with large deflections, the nonlinear and coupled Föppl–von Kármán (FvK) theory of shallow shell elastic deformation [8] provides a suitable theoretical framework for ascertaining their behaviour. For more intricate shells with bespoke support conditions, the finite-element (FE) method enables their analysis straightforwardly but not without computational expense or the complete assurance that all possible stable states have been detected. Correspondingly, the FvK framework remains competitive and affords direct insight into the factors which govern multistability. Often, its governing equations can be intractable unless simplifying assumptions about the shape of the shell are made. One successful approach describes the dominant deflections by *uniform* curvatures (UC) and twisting curvatures for the sake of simplicity. Approximating the shape globally this way, however, fails to capture local edge effects [9], leading typically to a non-vanishing bending moment at the boundary. Despite this, the UC approach provides a good approximation (≈ ±10%) of the ratio of bending-to-stretching energies for shells with constant thickness [2,3,7,10–13], which is suitable for most design purposes. For a lenticular shell, where a vanishing thickness precludes edge effects, Mansfield [14] obtains an exact solution for ‘extra’, i.e. bistable states.

More detailed numerical approaches use higher-order polynomial basis functions for describing the large deflection shape of shells with greater accuracy [15,16]. Other elaborate analytical models also consider higher-order curvatures, for example, Vidoli [17] demonstrates the increased accuracy of models with three degrees of freedom and a quadratically varying curvature (QVC). Other recent approaches show that the snap-through processes of magnetically actuated caps [18] or the multistable behaviour of a one-sided clamped plate [19] can be captured analytically with higher-order models. To some extent, these models show superior accuracy because a violation of the out-of-plane boundary condition (BC) is minimized at least or completely prevented.

In this paper, we consider the bistable behaviour of linear elastic, axisymmetrical shallow shells using polynomial displacement fields. One aim is to establish the trade-off between performance accuracy and model complexity as the polynomial order is increased. By enforcing axisymmetry and thus displacement fields that only vary radially, our analysis is simplified without compromising the generality of results. Our main aim, however, is to explore the effect that support conditions have upon the bistable character because the capability of any multistable structure is realized only by connecting it externally. Not only do we couple to the first aim insofar as capturing these BCs exactly (where the UC model is obliged to neglect the edge moment), we may observe if bistability is directly enhanced or eroded by the nature of the external connection—in this case—to ground. This has not been formally addressed although some researchers have dealt with the ‘limits’ of connectivity, where an edge is pinned or clamped as part of their particular study, and others have noted the importance of particular support conditions in stability problems experimentally [20]. We offer a broader sense of performance by implementing linear elastic connections of arbitrary stiffness, which behave as spring-like foundational elements.

Our approach derives a strain energy potential for the deformed shell whose stable equilibria are identified by familiar calculus of variations. We seek load-free, statical configurations without considering the unstable transition between them, and we use three examples to highlight the effectiveness of our methodology. The first is a familiar spherical cap, which is mounted on ‘pinned’ rollers at its edge so that bistable inversion takes place without constraint. This example can be directly compared to similar caps in the literature, and we determine the critical geometry at which bistable inversion *just* becomes possible in terms of the minimum height of the cap relative to its thickness.

In the second example, linear springs connect the roller supports of the first cap to ground. These springs apply radial forces to the cap that resist in-plane displacements at the edge, where in-plane radial stresses are no longer zero. Despite this simple modification in how the cap is supported, the critical height is found to decrease in a nonlinear fashion as the stiffness increases. In the final example, rotations of the shell boundary are resisted by torsional springs that locally apply a radial bending moment. We consider a non-uniformly curved shell initially with no boundary gradient but having some height because this allows us to solve efficaciously for the deformed shape by superposing two deflection fields with specific edge characteristics: one with zero bending moment, the other with zero slope. The contributions from each field are obtained by matching the actual bending moment and rotation of the cap edge to the spring element. This approach is novel and affords a reasonably compact working compared to one that uses a general polynomial displacement field, e.g. as in Wittrick *et al.* [2]. This shell is also supported on in-plane extensional springs for a completely sprung boundary and, hence, broader insight. A nonlinear performance of the critical geometry is revealed where we can observe, for example, the effect of the limits of connectivity for stiffness extrema. We constrained our investigations to closed shells consisting of isotropic materials with external connections at the *edge*, but point out possible extensions for more complex constitutive laws [21,22], support conditions [23,24] or geometries [25] (e.g. polar-orthotropic annuli on elastic foundations).

The structure of this paper is as follows. In the next section, we derive the general and particular solutions of the FvK governing equations. All results are then compared to FE simulations in §3, where we emphasize the stark relevances of different BCs; the FE model is briefly described in this section. A summary and outlook follows in §4.

## 2. Analytical formulation

We undertake a higher-order Rayleigh–Ritz method that extends previous UC models for load-free shells, e.g. as found in Seffen & McMahon [26]. A systematic overview of the approach is given in figure 1, which is elaborated next.

By assuming a polynomial expression of the transverse deflection field, *w*, on a domain *Ω* over the shell with boundary, *Γ*, we describe the curvatures, *κ*, the deflection gradients, *φ*_{r}, the bending moment, *m*, the Gaussian curvature, *g*, as well as the bending energy, *Π*_{B}, in terms of *n* degrees of freedom (d.f.). Out-of-plane curving is coupled to the in-plane stresses by relating the change in Gaussian curvature to the Airy stress function, *Φ*, using Gauss’s Theorema Egregium [9]. From *Φ*, the in-plane stress resultants, *σ*, the mid-plane strains, *ε*, as well as the displacement, *u* and stretching energy, *Π*_{S}, are derived, while satisfying the in-plane BCs. Finally, stable configurations are identified by minimizing the total strain energy, *Π* equal to *Π*_{B}+*Π*_{S}. The general solution is given next; the particular solution, which elucidates satisfaction of the BCs on *Γ*, is given afterwards.

### (a) General solution

Our shell is axisymmetrical and can be connected to ground everywhere along its circumference by linear elastic springs. Because deformations are also assumed to be axisymmetrical, the springs resist in-plane extensions and out-of-plane rotations of the shell boundary purely in the radial direction without circumferential variation. A shallow spherical segment is depicted in figure 2*a*, with a through-thickness *t*, and an outer planform radius *a*. We choose this segment only to illustrate more easily the cylindrical coordinate field given by *r* and *θ*, and by *z* in the direction of transverse *w*. For a radial variation alone and no initial deflections, the FvK plate equations can be written as [28,29]
*a*and
*b*where ∇^{2} is the familiar Laplacian operator, equal to d^{2}(..)/d*r*^{2}+(1/*r*)d(..)/d*r*.

These expressions are nonlinear and coupled, and are valid only when the radial curvatures *κ*_{r}≪1/*t*, and the gradient of the shell d*w*/d*r*≪1, everywhere. The external transverse pressure load is denoted by *q*, and the flexural rigidity is *D*=*Et*^{3}/[12(1−*ν*^{2})], with Young’s modulus, *E*, and Poisson’s ratio, *ν*. Even though large out-of-plane deflections are considered, the shallow shell condition arises because the strain formulation
*w*, and later define a surjective and not injective stress function, *Φ*(*w*), which implies that there is not a unique inverse function, *w*(*Φ*).

### Remark 2.1

The equilibrium equation, (2.1a), is in general not satisfied for polynomial approaches of *w* and *Φ* with a finite number of terms. For an alternative state of self-stress, *q* is absent, and hence, for a deflection field, *w*, of order *p*>1, the first term in equation (2.1a) is of order *p*−4, while the second terms order is *p*−3+deg(d*Φ*/d*r*). For matching orders, the Airy stress function, *Φ*, has a logarithmic component to enforce the compatibility for the highest order term in equation (2.1a); however, such an Airy stress function invokes the incompatible behaviour of infinite strain energy due to infinite stress peaks at the centre of the plate (cf. equation (2.6)). This problem cannot be overcome by increasing the number of terms, since every additional term also requires an additional correction term, and a mismatch always remains. By contrast, polynomial approaches are suitable for the second FvK equation, (2.1b), with 2deg(*w*)=deg(*Φ*). As equilibrium has the same axiomatic nature as work, we ignore the equilibrium equation, (2.1a), by assuming some (as yet unknown) deflection instead, and find an approximate *yet* accurate solution via energy minimization.

#### (i) Out-of-plane curving

The initially stress-free shape, *w*^{0}, is distinguished from the current shape, *w*, by a superscript 0; it is specified by a function *change* between the two configurations is denoted by a hat, so the deflection reads *f*^{0}(*r*), is an arbitrary polynomial expression of order *p*, and both are specified in §2b. For axisymmetrical structures, the change in gradient, *m*_{r} and *m*_{θ}, as well as the shear force, *q*_{r}, are independent of the initial configuration, so the change in these values corresponds to the current value, i.e.

#### (ii) In-plane stretching

The Airy stress function, *Φ*, is a bi-potential function for describing the mid-plane stress resultants as

#### (iii) Bending–stretching interaction

The in-plane and out-of-plane responses are connected by equating the extrinsic definition of Gaussian curvature to in-plane strains according to Gauss’s Theorema Egregium [9]:
*change* in Gaussian curvature
*p*−4 in *r*. By introducing the coefficients *α*_{i}, this term can be sorted by order:
*r*:
*C*_{1},*C*_{2} and *C*_{3} are now used to satisfy the in-plane BCs.

### (b) Particular solution

While the initial shape can have the form of any arbitrary polynomial expression, we choose approximations of two very well-known examples for evaluation in §3 and enforce zero displacement at *r*=*a*:
*a*or
*b*where *ρ* the dimensionless radius *r*/*a*. Shape (a) resembles a uniformly curved shallow cap and is inspired by a study of initially curved bistable beams [6] using a half cosine-wave; (b) represents a full cosine-wave akin to the deflection field of a plate with a clamped edge [30–32], where the gradient is zero. These shells are subjected to two possible types of BCs: Dirichlet BCs, which impose displacements or rotations on a boundary *Γ*_{D}, and Neumann BCs, which prescribe generalized stresses on *Γ*_{N}. When springs are connected to the boundary, these conditions become coupled, and the specific stresses (from in-plane forces or bending moments) are related to an in-plane displacement or rotation via the spring stiffness, *K*_{u} or *K*_{φ}, respectively, by
*a*and
*b*The minus sign in both expressions accords a resistive force or bending moment for positive directions of *u*_{r} and *a*
*b*
*c*
*d*
*e*and
*f*

The first three supports are directly employed in the examples of §3. The fixed pinned support and the clamped support are equivalent to having infinite spring stiffness, so by letting *K*_{u} and *K*_{φ} tend to large values in equations (2.16b) and (2.16c), we may observe how the ideal support conditions are approached from a ‘spring’ perspective. The final condition ensures that we deal with plates without central holes: if we wished to consider planform annuli, we may substitute this final condition with one from any of the five expressions above it but written in terms of the inner radius of the hole; see Sobota & Seffen [25] for details.

#### (i) Out-of-plane bending

In order to describe the transverse deflections of a shell bounded by rotational springs, *n* degrees of freedom, *η*_{i}, which by itself is the solution for cases where the boundary of the shell is pinned for all time. The rotational spring equation, equation (2.15b), prescribes the compatibility between these two fields and is later used to express the clamped midpoint deflection, *n* degrees of freedom, *η*_{i}; cf. equation (2.26). Correspondingly, the total midpoint deflection is

*Hinged subset*. The hinged subset has to satisfy the BCs of a vanishing bending moment and deflection at *ρ*=1. To ensure this, we combine equation (2.4) with equation (2.5) and choose a polynomial series of even powers for the radial bending moment, i.e.
*p*=2*n*+2, as:
*η*_{1}, these curvature and bending moment expressions have a quadratic variation in *ρ*.

*Clamped subset*. The clamped subset is chosen to be [30,32]
*ρ*=1. Note also that the two conditions concerning the deflection in equation (2.16e) are satisfied for arbitrary values of the midpoint deflection *w*^{0}, and the hinged deflection,

*Resulting deflection field* *ρ*=1, which enables us to find *ρ*=1, *m*_{r}=0 for the hinged subset, so the bending moment is solely evoked by the clamped part in equation (2.23),
*n* hinged degrees of freedom, *η*_{i}. For a vanishing spring stiffness, *K*_{φ}=0, the clamped part vanishes; for *K*_{φ} tending to infinity, the edge rotation becomes very small but is not allowed to vanish because the solution is achieved with respect to the hinged degrees of freedom. The consequences of this assumption are addressed in §3.

### Remark 2.2

For roller-supported caps, a UC approach is also investigated in this paper by setting

### Remark 2.3

The midpoint deflection of the clamped mode is expressed in terms of *η*_{i} from equation (2.26) for *n* degrees of freedom but where the corresponding deflection variation in *ρ* is always of order four; cf. equation (2.21). For the hinged mode in equation (2.18), the variation in *ρ* increases with *n*, up to order 2*n*+2. We do not impose a similar polynomial expansion in *ρ* for the clamped mode because we find that it is simply not needed in view of producing accurate results.

#### (ii) Bending–stretching interaction

The change in Gaussian curvature is calculated according to equation (2.11) with *κ*_{r}, *κ*_{θ}, *α*_{i} terms, equation (2.12), are given for a three degree-of-freedom model with rotational edge springs; lower-order models are straightforward to derive from these expressions by setting ‘unused’ degrees of freedom to zero. Since the product of initial curvatures also features in *α*_{i} terms are given for both initial shapes expressed by equation (2.14a) and (2.14b). For instance, the curvatures for a hinged (*K*_{φ}=0), uniformly curved cap initially *n*=1) are equal to
*ρ*, leading to a change in Gaussian curvature:
*α*_{i} terms, equation (2.12), with 2*p*−4=4*n* (+4) now read
*K*_{φ}=0, which ought to be expected.

#### (iii) In-plane stretching

From the Airy stress function in equation (2.13), the in-plane stresses are calculated via equation (2.6) to be
*ρ*=0, the constants *C*_{2} and *C*_{3} have to be zero and, thus, equation (2.30) is simplified to
*σ*_{αr} and *σ*_{αθ} denote the summation terms with *α*_{i}. These expressions also satisfy the condition *u*_{r}|_{ρ=0}=0, which is equivalent to having a biaxial stress state, *σ*_{r}=*σ*_{θ}, at the centre of the plate. The remaining constant *C*_{1} is determined for a linear elastic in-plane extensional spring by using *u*_{r}=*rε*_{θ} with
*σ*_{r} from equation (2.31) into the extensional spring relationship, equation (2.15a), with *ρ*=1 gives
*C*_{1} can be substituted into the equations for the stresses and strains, in order to calculate and then minimize the strain energy.

### (c) Identifying stable configurations

All stress- and strain resultants are known and, hence, the bending and stretching energy components, *Π*_{B} and *Π*_{S}, respectively, which include the spring contributions, can be calculated via
*Π*=*Π*_{B}+*Π*_{S}, where
_{η} is the differential operator with respect to the *n* degrees of freedom, *η*_{i}. These configurations are stable if, and only if, all eigenvalues of the strain energy function are positive, which is assured by a positive definite Hessian matrix of stiffness, **H**, where

## 3. Results

Recall that the first example is a spherical cap supported on rollers according to the BCs in equation (2.16a). The second is another spherical cap connected to ground by in-plane springs, equation (2.16b), and the third is a non-uniformly curved cap with initial deflection given by equation (2.14b) and supported by both extensional and rotational spring supports, equation (2.16c). For all cases, the variation with key physical parameters is established as a function of the number of degrees of freedom while comparing their accuracy to FE results.

FE simulations were conducted with the commercial package ABAQUS [27] using a quasi-static implicit dynamic analysis. Only one-quarter of the shell was modelled, with biaxial symmetry applied to over 600 elements and the following parameters: *E*=10^{7}, *t*=0.01, *a*=1 and density equal to 10^{−5}; cf. Sobota & Seffen [33] for details. Mesh refinement of randomly picked samples did not lead to any changes in the critical properties.

### (a) Roller-supported spherical cap

The critical dimensionless height, *ν* is governed by naturally stable limits of *ν*=−1 and *ν*=0.5 for isotropic materials: for homogeneous shells, the Young’s modulus does not affect the critical geometry. Figure 4*a* first compares results from previous studies to FE for the sake of revision before divulging our predictions in figure 4*b*. Although we do not plot them, the results for the UC assumption in Vidoli [17] and in Seffen & McMahon [26] are therefore identical to the current case using *a*,*b*, where an initial height-to-thickness ratio of at least six guarantees bistable inversion irrespective of *ν*. The QVC model of Vidoli [17] shows a better approximation of the FE solution for positive Poisson’s ratios in figure 4*a* while our QVC results in figure 4*b* are superior for negative ratios, which may be exhibited, say, by auxetic materials. Wittrick *et al.* [2] deduce that the critical height must exceed *ν*=1/3 and *ν*=1/4; our QVC model returns lower values of 4.13 and 3.92, respectively. Even though Mansfield [14] provides a closed-form solution of *a* because, astonishingly, his solution predicts the critical initial deflection precisely (±0.39%) for *ν*≥0, but for negative Poisson’s ratios the deviations rise up to 12.8%. When we increase the number of degrees of freedom to three, we find a better approximation of *b* for all values of *ν*.

### (b) Spherical cap with extensional spring supports

We first set the Poisson’s ratio to be 0.5 because this correlates to the *worst* case in figure 4 in terms of the accuracy of our QVC model for a roller-supported cap. The performance of this model in the present case is indicated in figure 5 along with UC predictions (cf. remark 2.2), which is also a single degree-of-freedom approach. Closed-form solutions of the critical bistable height, *K*_{U}=0, and to 16/(7−*v*) when

The reason for including the simpler UC model is now apparent in figure 5 because it yields a *better* approximation of the critical geometry than the QVC model and is surprisingly close to the FE trend. However, Sobota & Seffen [25] show that this holds only for the *ratio* of bending- to stretching energies—that the close accuracy of the UC model is fortuitous, and that our single degree of freedom QVC model gives more accurate variations of the stress resultants. The predictions from using two and three degrees of freedom are also plotted, and are almost identical to the FE solutions with maximum absolute deviations of 0.78 and 0.73%, respectively.

When *K*_{U} is small, in-plane displacements of the boundary are largely unrestricted, and the cap tends towards being roller-supported, equation (2.16a). The most accurate prediction of the associated critical height is 5.5 times the thickness, which matches that in figure 4. As *K*_{U} increases, the critical height decreases; moreover, at *K*_{U} ≈ 1 we observe that *lowest* Poisson’s ratio value: by stiffening the cap edge, the critical height can be smaller than that of the first case in view of promoting bistability. A discernible lower bound of *Et*/*a* acting on a length of 2*πa*. When *K*_{U} becomes larger, in-plane boundary displacements are increasingly resisted, and the cap is supported in a manner akin to a fixed pinned support, equation (2.16d).

The detailed dependency on Poisson’s ratio is given in figure 6*a*, where a similar accuracy is confirmed from using three degrees of freedom. The variation in critical height is now plotted as a ‘landscape’ with respect to *ν* and to *K*_{U}. Its discrete contours in figure 6*b* allow values of

We impose the same range for *K*_{U} from figure 5, noting little variation in *K*_{U}=10. We plot against the logarithm of spring stiffness for compactness but this also allows us to infer robustly the asymptotic performances. For the case of a roller support (*K*_{U}→0), the critical heights are largest and the variation with Poisson’s ratio is the most distinct. As *K*_{U} increases, *beam* [6]. While the beam is bistable if the ends are pinned to the ground and *quantified* for shells: depending on the Poisson’s ratio, the critical initial height for a roller-supported cap must be between two and four times larger than that of a fixed pinned cap.

Out of interest, we have performed separate FE analysis of other shells with rhomboidal and elliptical planforms of aspect ratios up to 7 : 1. Their initial shapes conform to bi-directional sinusoidal profiles with fixed pinned supports along their boundary plane. For both types of shell, we find that the critical bistable height lies in the range *ν*>0 for the rhombuses and independently of *ν* for the ellipses. These limits are similar to those found in figure 5 as *K*_{U} becomes large, which suggests a universal height for bistability of around two thicknesses for this type of boundary connection irrespective of the planform shape.

### (c) Non-uniformly curved shell with extensional and rotational spring supports

The UC approach is no longer suitable since neither the initial shape nor the deformation field have UC, recall equation (2.14b). Sample results for different degrees of freedom are displayed in figure 7 and compared to FE predictions for the case of *K*_{U}=1 and *ν*=0.3. The range of rotational stiffness, *K*_{ϕ}, is specified from a very small value up to *K*_{ϕ} ≈ 100, where each variation in *K*_{ϕ}<1, the error in *ν*=0.3.

Note that the axes variation for *K*_{U} has been reversed in figure 8*a* for a more open perspective of the solution landscape, and that larger values of *K*_{U} have been used compared to previous figures, in order to reveal the pertinent variation. The landscape topography is clearly nonlinear, for example, *K*_{U} increasing for a given *K*_{ϕ} but the trend in general reverses for *K*_{ϕ} varying. The behaviour for *K*_{U} varying mirrors that of the second example despite it having a different initial shape and no rotational support: that by resisting in-plane edge displacements, bistability can prevail. A rotational spring tends to revert the shell by definition since it applies an edge-bending moment in the sense of the initial shape. Hence, by reducing the applied moment by decreasing *K*_{ϕ}, we also promote bistability. Closer inspection reveals an anomalous ‘dip’ in the downward trend of *K*_{ϕ} for the largest values of *K*_{U}. Importantly, the critical height is lowest for a non-zero value of *K*_{ϕ} (≈0.39), which suggests that a moderate value of rotational spring stiffness benefits bistability—here, this stiffness is equivalent to 4.26*D*/*a*. We surmise that this attests to the highly nonlinear interaction between stretching and bending in shells, and to their sensitivity to edge effects, which we hope to quantify in further study. A discrete contour plot is also given in figure 8*b*, where the asymptotic support conditions are also indicated in the four corners of the plot. In the sense of in-plane *versus* rotational freedoms, these are fixed–fixed (*K*_{U} & *K*_{ϕ}→0).

### (d) Final remarks

Because our deflection field comprises a clamped deformation mode and a hinged (pinned) mode, we cannot capture an ideal clamped support outright because the hinged mode is inherently coupled to the overall solution. We can only *approach* an ideal clamp by setting *K*_{ϕ} to be a very large finite value. In the earlier derivation, the total strain energy, *Π*, ultimately expresses *K*_{ϕ} to the fourth power (but not so for *K*_{U}), which compounds further the effect that large values of *K*_{ϕ} can have upon the numerical solution procedure, especially when solving for the eigenvalues of the Hessian matrix, equation (2.37). A sensible limit for *K*_{ϕ} should be correlated to the range of numerical floating-point precision levels of the particular software package. For basic packages, this range is 16 decimal places (or ‘16-point’), hence setting *K*_{ϕ}≤10^{4}) to be indistinguishable from those of Sobota & Seffen [25] who deal with a clamped support shell from the outset.

We have chosen specific initial shapes because of their ubiquity in other studies. Other shape profiles are indeed amenable to our method, but we have to consider the nature of our approach, which differs from FE analysis: both methods can theoretically lead to ‘false negatives’ or to ‘false positives’, but the first is more likely in FE analysis, while our choice *not* to linearize the work equation increases the possibility of the latter. Using FE, we cannot exclude the possibility of alternative equilibrium states since the linear stiffness matrix would require us to evaluate an infinite number of possible combinations. However, because of linearization, the problem becomes amenable to a numerical solution process even for discretizations with a large number of degrees of freedom. The current method considers a nonlinear stiffness matrix, **H**, which implies that computational efficiency limits us to a few degrees of freedom; hence, our set of shape functions is not necessarily accurate enough. This problem can be illustrated by the example of a uniformly curved cap now clamped on its edge, where the gradient remains non-zero for all time. Although we never observed bistable behaviour in FE simulations for a wide range of initial geometries, our theoretical analysis suggests its feasibility. In the FE simulations, we observe that the radial curvature becomes highly concentrated just before the clamped edge, which gives way to a large restoring moment despite the ‘holding’ effects of significant in-plane circumferential tension. Our lower-order models fail to capture the strongly curved domain next to the edge, but simply increasing the order is not a viable solution approach because the midpoint deflections and, more importantly, deflection gradients (recall d*w*/d*r*≪1) increase beyond the limits acceptable for shallow shell behaviour and thus beyond the scope of this study.

## 4. Summary

We have addressed the bistable properties of an initially stress-free shell supported on its circular edge by extensional and rotational springs. Our analytical solutions have been derived from a Rayleigh–Ritz method with a polynomial series up to four degrees of freedom. The FvK compatibility equation has been satisfied by making use of Gauss’s Theorema Egregium, and the connection between the Gaussian curvature and the Airy stress function has been established independently of the in-plane BCs. Solutions for the uniform- and QVC approaches have been captured in closed form and are accurate for most design processes when the edge is hinged. In all results, having a radial spring stiffness promotes bistability; in the example of a roller-supported cap, the critical initial midpoint deflection for bistability has been shown to be between two and four times higher than that of the same pinned cap but where radial edge displacements are now prevented. While the response of the former was highly dependent on the Poisson’s ratio, the latter was almost independent of *ν*. Incorporating a rotational spring stiffness has shown that the response of the structure produces a strongly non-uniform Gaussian curvature field, and therefore, higher-order models with at least three degrees of freedom are needed to capture the bistable response. It has been demonstrated that a horizontally immovable boundary in combination with a specific value of rotational spring stiffness, which is not an extreme case, can minimize the critical initial height needed to achieve bistability. Compared to structures with hinged edges, fully clamped boundaries tend to prevent bistability.

## Data accessibility

This work does not have any physical experimental data. The essential codes for finite-element analysis [33] via ABAQUS [27] and for solving the analytical models [35] via Mathematica [34] are available in the cited repositories and for open consultation, which also contain the data for generating the results presented in figures 4–8.

## Authors' contributions

P.M.S. and K.A.S. conceived the methodology, interpreted the results and wrote the paper in equal measure. P.M.S. performed the finite-element simulations and wrote the Mathematica code, all in consultation with K.A.S. Both authors gave their final approval for publication.

## Competing interests

The authors declare that they have no competing interests.

## Funding

P.M.S. is grateful to the Friedrich-Ebert-Foundation for their financial support; this work was otherwise not supported by funding from any other agency.

## Acknowledgements

Both authors are grateful to Prof. Stefano Vidoli for explaining his QVC model at the start of the study.

## Appendix A

For a *n* degree of freedom model, the current radial and circumferential curvatures, *κ*_{r} and *κ*_{θ}, respectively, read
*α*_{i}-terms for equation (2.12), where 2*p*−4=4*n*, for three degrees of freedom and a uniformly curved initial shape (cf. equation (2.14a)) read
*α*-terms of a model with three degrees of freedom are calculated to be

- Received April 10, 2017.
- Accepted June 21, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.