Morphing shell structures: A generalised modelling approach

Morphing shells are nonlinear structures that have the ability to change shape and adopt multiple stable states. By exploiting the concept of morphing, designers may devise adaptable structures, capable of accommodating a wide range of service conditions, minimising design complexity and cost. At present, models predicting shell multistability are often characterised by a compromise between computational efficiency and result accuracy. This paper addresses the main challenges of describing the multistable behaviour of thin composite shells, such as bifurcation points and snap-through loads, through the development of an accurate and computationally efficient energybased method. The membrane and the bending components of the total strain energy are decoupled by using the semi-inverse formulation of the constitutive equations. Transverse displacements are approximated by using Legendre polynomials and the membrane problem is solved in isolation by combining compatibility conditions and equilibrium equations. This approach provides the strain energy as a function of curvature only, which is of particular interest, as this decoupled representation facilitates efficient solution. The minima of the energy with respect to the curvature components give the multiple stable configurations of the shell. The accurate evaluation of the membrane energy is a key step in order to correctly capture the multiple configurations of the structure. Here, the membrane problem is solved by adopting the Differential Quadrature Method (DQM), which provides accurate results at a relatively small computational cost. The model is benchmarked against three exemplar case studies taken from the literature. ∗corresponding author: ettore.lamacchia@bristol.ac.uk Preprint submitted to Composite Structures June 20, 2015


Introduction
Multistable structures may play an important role in future engineering designs due to their potential for reconfiguration between different states. Shape changing concepts may allow the designer to devise morphing structures able to adapt to a wide range of service conditions, minimising the need for complex multiple interacting components.
A large variety of morphing concepts have been developed. These often utilise stiffness-tailored nonlinear structures to provide shape changing features. Notable examples include non-symmetric composite plates or prestressed laminates [1,2,3,4,5,6,7]. A variety of numerical (e.g. based on Finite Element Analyses) and semi-analytical nonlinear models to describe and design these structures can be found in the literature. This paper focuses on the latter kind, which are either bespoke to specific problems, and of limited applicability, or more general which often necessitate a compromise between computational efficiency and accuracy. For example, Hyer [1] and Hamamoto et al. [8] proposed energy-based models for bistable plates that assume laminates to be initially flat. Subsequently, Galletly et al. [9] and Guest et al. [10] relaxed this assumption, introducing different restrictions on permissible constitutive equations or initial Gaussian curvatures. More recently, Eckstein et al. [11] extended Hyer's work to describe the multimode morphing of thermally-actuated initially cylindrical shells, in which temperature-dependent material properties were also modelled. In all of these works, simple shape functions were employed to keep computational cost to a minimum, with the limiting assumption that transverse displacements could only capture constant curvatures, whereas in reality, satisfaction of zero forces at free boundaries creates significant nonlinear variation of curvatures into the interior of the structure [12,7]. Conversely, Pirrera et al. [13] developed a shell model employing higher order polynomials to approximate displacements more accurately than had been done previously. The large number of degrees of freedom allowed the restriction to constant Gaussian curvature to be relaxed. The authors could then demonstrate that refined shape functions are required for the accurate prediction of bifurcations, buckling, snap-through loads and, in general, of the deep post-buckling behaviour of morphing panels. However, the high order of the set of complete polynomials implied a large computational cost, thus a loss of efficiency of the method.
Guest et al. [10] were the first to decouple the membrane and bending components of the total strain energy to describe the equilibria and stability of a composite cylindrical shell with constant curvature. The advantage of this approach is that the resulting membrane and the bending problems can be solved in independently, with a consequent reduction of total computational cost. Based on this observation, Seffen [14] applied the approach proposed in [10] to plates with constant curvatures and elliptic planforms, whilst Vidoli [15] added quadratically varying curvatures.
By assuming constant curvature over the surfaces with elliptic planform, Seffen [14] could find the closed form solution for the membrane problem by combining the compatibility condition with the in-plane equilibrium equations. He solved the compatibility condition in terms of the Airy stress function by expressing the membrane stress as a second order polynomial, which upholds the boundary conditions of zero force and zero shear, respectively normal and tangent to the free boundary. Vidoli [15] relaxed the assumption of constant curvature over the surface and extended Seffen's work to model rectangular planforms and composite materials. To solve the governing equations, Vidoli adopted a Rayleigh-Ritz approach and deliberately imposed symmetry whilst forcing the bending natural boundary conditions to be satisfied on average.
Our approach builds upon [10,14,15] by considering thin anisotropic laminated shells possessing general planforms and initial curvatures. In this paper, following the work of Guest et al. [10], Seffen [14] and Vidoli [15], the total strain energy is expressed as the sum of the stretching and bending contributions, which are decoupled via the semi-inverse formulation of the constitutive equations [16]. Upon decoupling, one can solve a set of independent membrane problems by combining the membrane equilibrium equations and the compatibility condition that links strains and curvatures according to Gauss' Theorema Egregium. The resulting system of equations gives a strain field which is a function of the curvatures only, which are subsequently used as Lagrangian parameters, i.e. degrees of freedom. By substituting the strain field into the stretching energy, the total strain energy also becomes a function of the curvatures only. As with the Rayleigh-Ritz method, structural equilibria are found by minimisation of the total potential energy with respect to the Lagrangian parameters. The stability of the equilibria is then evaluated assessing the positive or negative definiteness of the system's Hessian.
Although we consider structures that display shell-like behaviour, the mathematical assumptions behind the present model are those of initially curved plates. In other words, we consider the kinematics of flat plates, with initial curvature, whereby bending and membrane behaviour is uncoupled for linear response but becomes coupled for the nonlinear case, as facilitated by large displacements. Simultaneously, we consider that membrane and bending coupling can occur through material anisotropy via linear response or nonlinear or indeed a combination of the two. These assumptions are valid if the shell is sufficiently geometrically and kinematically shallow, as used, for example, by Mansfield in describing the bending behaviour of initially curved strips [17].
Vidoli and Maurini [18], suggested that the solution of the membrane problem over general planforms calls for a numerical solution. They also noted that the accuracy of such a solution is crucial for the correct estimation of the relative importance of the bending and stretching contributions to the strain energy. Following these observations, the Differential Quadrature Method (DQM) is adopted here for its proven accuracy and computational efficiency [19]. In contrast to Guest et al. [10], Seffen [14] and Vidoli [15], the speed and robustness of DQM, combined with the generality of Legendre polynomials, allows a general description of the displacement field over the computational domain to be modelled.
The Differential Quadrature Method was introduced by Bellman and Casti [20]. It is based on the premise that any continuous function can be approximated with a high-order polynomial and that the function's derivatives can be expressed as a linear combination of functional values at given grid points distributed over the domain. Owing to the high-order polynomial approximation, DQM usually requires fewer degrees of freedom in comparison to other discretisation methods, such as Finite Element Method (FEM) or Finite Difference Method, to achieve accurate results [21]. In this sense, and unlike other semi-analytical approaches [22,23,7], DQM guarantees accuracy at a relatively computational cost. This paper is structured as follows. In sections 2 and 3, we derive a generalised model for initially shallow, anisotropic shells of aribitrary planform that undergo large displacements and can, but not necessarily, change shape from one stable configuration to another. In section 4 the systems of nondimensionalised equations that have been developed are validated against available benchmark results, thereby showing our approach's characteristic features and advantages. These are summarised as follows.
Accuracy of results. The transverse displacements are approximated using high-order Legendre polynomials. This discretisation allows the membrane and the bending components of the total strain energy to be captured with fewer Lagrangian parameters, yet accurately, in comparison to using non-orthogonal polynomials [7,11].
Symmetry and boundary conditions. As in [15] and in contrast to conventional Rayleigh-Ritz approaches, the boundary conditions on the inplane stress resultants are satisfied exactly and point-wise. Ordinarily, the real value of the bending moments at the boundary is approached asymptotically. However, as a novel contribution, by making use of DQM we are able to adopt a systematic approach for the solution of the structural problem that does not require any assumption on the symmetry of the displacement field and captures natural bending boundary conditions without additional approximate constraining relations.
Computational efficiency. Our model is computationally inexpensive, because, given a geometry, the relative domain and computational grid, and a load, the membrane problem is solved once and for all inverting one small, sparse matrix of DQM weighting coefficients. This advantage is illustrated in section 2.3 and is in contrast to previous works, in which a set of matrix inversions had to be carried out. The efficiency of the proposed model is benchmarked against results in the literature and becomes particularly evident, for instance, trying to recover the snap-through diagram of bistable plates. We achieve greater accuracy with the additional benefit of fewer degrees of freedom when compared with previous work [7].

Kinematic Equations and Total Strain Energy
In a Cartesian reference system x-y-z, the strains are defined in the usual form as a nonlinear combination of linear deformations [24]. Under the assumption of small strains and moderate rotations, the nonlinear deformations for thin shallow shells are fully described by the Föppl-von Kármán formulation of the strain field ε, defined as [25]: (2) are the mid-plane strains and curvatures, representing flat plate kinematics, respectively. Here, u, v and w represent the displacements along x, y and z.
According to Claperyon's theorem [26] the system's total potential energy is defined as where, following Classic Laminate Analysis (CLA) and conventional notation, the stress tensor is defined as whereQ is the reduced stiffness matrix, Ω is the surface area of the plate and h its thickness. The subscript 'r' indicates contribution of non-mechanical deformations, such as thermal, piezoelectric or hygroscopic. By using (1) and (4) and integrating through the thickness, equation (3) can be recast as where A, B and D are the in-plane, coupling and bending stiffness matrices, ∆k = k − k 0 is the change in curvature from the reference/undeformed configuration, having curvature k 0 , to the current/deformed configuration, 6 and N r and M r are the in-plane stress and moment resultants due to nonmechanical deformations. For a complete derivation the reader is referred to references [25,16]. By substituting the semi-inverse formulation of the constitutive relations [16] into (5), the total strain energy becomes where is the in-plane stress resultsnt and M = M xx M yy M xy is bending moments resultant.
Notably, equation (6) allows the elastic component of total strain energy to be decoupled into a membrane and a bending part. This can be seen observing that the stiffness matrix in (7) is now diagonal, with bendingstretching couplings embodied in D * .
In other words, we now have an energy formulation in which the terms N A * N and ∆k D * ∆k fully and independently capture the membrane and bending behaviour of the panel. Essentially, this formulation reflects the linear distribution of strain through the thickness that is assumed to exist in a thin-walled structure but where the neutral axis does not necessarily lie at the mid plane. As usual, structural equilibria may be found by minimisation of the total potential energy. In the following sections, we describe the membrane problem and its use in the determination of overall equilibrium and stability.

The Membrane Component of the Total Potential Energy 2.2.1. Theorema Egregium and Equilibrium Equations
In order to fully define the so-called membrane problem, two key ingredients are required. They are the compatibility condition and the in-plane equilibrium equations. These are now derived and cast in a notation compatible with DQM.
Compatibility condition.. A generic Riemannian manifold is completely characterised by its first and second fundamental forms [27]. The components of the first and second fundamental forms are not independent, since they must satisfy the compatibility (or integrability) conditions. In two dimensions, this conditions reduce to three non-vanishing differential equations [28], i.e. the Gauss Theorema Egregium and the Codazzi-Mainardi equations. These equations can be obtained from the Riemann tensor recast for a twodimensional space. Physically, the compatibility relations between strains and curvatures represent the conditions of uniqueness of the normal vector at each point of the surface.
In the adopted Cartesian reference system and assuming small strains and moderate rotations, the compatibility condition for a thin shell is Equation (8) links the in-plane variation of the strain field with the change in Gaussian curvature, given by the determinant of the curvature tensor. 1 Introducing the differential operator and using (6), the compatibility condition (8) becomes In-plane equilibrium.. The membrane equilibrium equations for a free-free plate [25] are 1 note that when writing det k, k is intended in its 2 × 2 form, i.e. k = k xx k xy k xy k yy .

8
where ∂Ω is the boundary of the plate and n its normal. Introducing the differential operators the equilibrium equations become

Formulation of the Membrane Problem
By combining equation (10) and (13), we define the membrane problem as where the differential operator L is defined as It is worth noting that equation (14) is linear in N and therefore a bijective (one-to-one) relation exists between curvatures and in-plane stress resultants. This means that, given a suitable spectral discretisation for the displacement field, the membrane problem can be solved independently and using methods of linear superposition.
To solve equation (14), herein, we discretise the matrix of differential operators defined in (15) by using the matrix formulation of DQM. A description of DQM in matrix form can be found in Lamacchia et al. [29]. By applying the DQM to equation (15), the differential operators L 1 , L 2 and L 3 become matrices of weighting coefficients. The weighting coefficients have been evaluated by following the approach presented by Shu [19]. In particular, the matrices of weighting coefficients for the first order derivatives in x and y have been evaluated by using Chebyshev polynomials, which have been shown to be numerically efficient and well behaved. The matrices of the higher order derivatives have been evaluated by using the recursive formula presented by Shu [19], which allows the weighting coefficients for the derivative of order p to be evaluated by simply applying p times the matrix multiplication rule to the weighting coefficients of the first order. Moreover, Runge's phenomenon in which oscillations occur at the edges of the domain has been avoided by evaluating the weighting coefficients on a Chebyshev-Gauss-Lobatto mesh grid [30,31,19].

Overall Equilibrium and Stability
The curvature field to be used in (14) is defined here by approximating the transverse displacements w with Legendre polynomials. Hence, where w 0 represents the reference (undeformed) shell geometry, P i (x) and P j (y) are Legendre polynomials of order i and j defined as [32] and here used as shape functions; and where, finally, q ij are unknown parameters to be determined. Legendre polynomials are a set of complete, orthogonal functions which have been shown to provide great robustness [33], both in terms of accuracy and flexibility in describing the deformed shape of laminates. Using (16), w is approximated with a polynomial of order 2n, composed of (n + 1) 2 shape functions and degrees of freedom, namely the P i (x) P j (y) products and the q ij parameters.
Assuming a displacement field via (16), the in-plane stress resultants are obtained as function of curvatures only by solving equation (14): with N · n = 0 on ∂Ω. Computationally, the boundary conditions are imposed by converting the expression N ·n = 0 on ∂Ω in DQ form and appending the resulting rows of ones and zeros to L. N is then found by taking the Moore-Penrose pseudo-inverse [34,35] of the extended matrix, which is rectangular.
It should be noted that the size of L depends on the size of the DQ grid only, whereas increasing the order of the Legendre polynomials only increases the number of terms and the complexity of f , which has a marginal effect on the computational cost. At this point in the procedure, f is a function of q ij and, therefore, so are the in-plane stress resultants. The solution, N , of the membrane problem can now be substituted into (7) to obtain an expression of the total potential energy that is function of q ij terms only.
Equilibrium configurations correspond to extrema of Π(q ij ), therefore satisfying the expression which results in a set of nonlinear equilibrium equations of the kind ∂Π ∂q ij = 0. This system is solved with Newton-Raphson schemes. The stability of the solutions of (19) is assessed with the Hessian of the total potential energy. Equilibria are stable if and only if

Non-dimensional Formulation
In order to reduce possible numerical ill-conditioning of the nonlinear model and to analyse the relative importance of each term in the governing equations, the following scaling scheme is adopted, in which the symbol tilde (˜) denotes non-dimensionalisation: where the dimensional scaling parameters • 2L x and 2L y are the side lengths of the plate along the Cartesian axes, • U d , V d and W d are coefficients defined as [7,36,37] • E and K are used to scale the in plane strains and curvatures, respectively, and are defined as • Σ scales the membrane stress resultant N and it is defined as Σ = AE.
For illustrative purposes we assume that the non-mechanical contributions in the total strain energy are due to a thermal load. The thermoelastic in-plane stress resultants and bending moments are then Following [7], where ∆T 0 = T − T ref is the thermal gradient from cure temperature,τ is a non-dimensional parameter defined from ∆T 0τ = T − T ref and it is noted that, in general,Q k (T ) andᾱ k (T ), the coefficients of thermal expansion, may both be functions of the temperature.
Substituting the parameters defined in equation (21) and (25) into equation (7), one obtains the non-dimensional form of the total strain energỹ where the non-dimensional matrices are defined as and is the parameter, defined as in [7], used to scale the total strain energy, so that Π = Π dΠ . For completeness, the membrane problem must be scaled and recast in terms of non-dimensional parameters too. Hence, the compatibility condition (8) becomes and by definingL and substituting the non-dimensional form of equation (6), equation (29) is recast as 13 Similarly, the equilibrium conditions (13) become Finally, combining (32) and (33), the membrane problem in non-dimensional form is given by The solution of the membrane problem is simply given bỹ which is a function ofx,ỹ and of the Lagrange parameters q ij . Analogously to section 2.3, by substituting equation (37) into (26) we obtain the total potential energy as a function of q ij terms only. The equilibrium positions are obtained by minimising equation (26) with respect to the Lagrange parameters.
It is interesting to observe that the operator expressed in equation (30) is a second order differential operator. Hence, in equation (32) the effects of any constant or linearly varying external field in x and y vanish. In other words, as an interesting observation, any external non-mechanical loading field contributes to the development of curvatures if and only if it is at least of the second order in x or y.
In order to visualise this property, we consider a constant and a quadraticallydistributed thermal field applied to the cylindrical shell depicted in Figure 1a. For the sake of simplicity, and to avoid any membrane-bending coupling due to anisotropic material properties, we consider an isotropic aluminium (Al  Figure 1a shows the isotropic aluminium (Al 6061-T6) cylindrical shell, having initial curvature of 3 m −1 about the x-axis. The shell is subjected to a uniform and quadratically varying thermal field over the surface. Figure 1b and 1c show respectively the displacements of the initially curved cross section (blue line in Figure 1a) and of the initially flat cross section (red line in Figure 1a) due to the thermal fields. From Figure 1c it is evident that the initially flat cross section deforms just in-plane under a uniform thermal load, whereas it deforms both with v and w due to the quadratically-distributed thermal field over the shell surface. From Figure 1b is evident that the initially curved cross section deforms with u and w independently of the distribution of the thermal field. Here u, v and w are the displacements along x, y and z respectively. 6061-T6) shell, having in plane dimensions 1×1 m and thickness 0.005 m. The initial curvature is 3 m −1 about the x-axis. The magnitude of the applied thermal gradient is 100 • C and the shell is clamped at the centre. Figure 1b and 1c show the deformations of the two orthogonal cross sections of the cylindrical shell, which are represented by the blue and red lines in Figure 1a. In particular, Figure 1b shows that the originally curved cross section exhibits both u and w components of deformation, independently of the applied thermal field. Conversely, by observing Figure 1c, it is evident that the originally flat cross section deforms just in-plane along the y-axis under a uniform thermal load, whereas it deforms with v and w components when a quadratically distributed thermal field is applied. Hence, the Gaussian curvature of the deformed configuration does not change with respect to that of the reference shape if the distribution of the external thermal field is uniform (or linear) over the surface. Conversely, the initially cylindrical shell deforms by changing its Gaussian curvature to accommodate the nonmechanical stress if the external field is at least of the second order in x or y.

Results
In this section, we present three case studies whose collective aim is to showcase the main features and advantages of our model in comparison to the state-of-the-art. For ease of comparison and to obtain meaningful benchmarks, we reproduce and complement results available in the literature.

Case-Study 1: Snap-Through Load
The computational efficiency of the proposed model is evaluated in terms of the number of degrees of freedom, or Lagrangian parameters, required to obtain converged results. Here, the work done by Pirrera et al. [7] on the thermally-induced bistability of thin unsymmetric plates is taken as a benchmark. The reason for this choice is twofold. Firstly, Pirrera et al. [7] have shown that a large number of degrees of freedom is necessary to capture the regions of the design space in which cross-ply laminates are bistable. Secondly, they have also shown, together with Diaconu et al. [38], that an even larger number is required to determine the snap-through load of the panel accurately. Naturally, the necessity of many degrees of freedom comes at the expenses of the computational cost.
Here, we focus on the ability of our model to capture the snap-through load of the plate, i.e. the load required to snap from one stable shape to the other. Pirrera et al. [7], adopted a Rayleigh-Ritz model, based on CLA, employing complete linearly independent polynomials to discretise the displacement field. For a correct estimation of the snap-through load, within 5% from converged FE results, polynomials of order 11 had to be used. This entailed 62 degrees of freedom, after symmetric displacements were imposed. As stated by the authors, this level of accuracy required an extremely high computational cost, and any further increment of the order of the polynomial approximation was impractical. Figure 2 shows the snap-through load for the [0 2 /90 2 ] plate studied in reference [7]. Material properties are listed in table 1. The plate is free at the boundaries and clamped in the middle. It has a 0.25×0.25 m square planform and the ply thickness is 0.131 mm. Four corner forces have been applied along the z-axis. The expression for the total potential energy, equation (26), is amended accordingly by adding the contribution of the external force. Initially, the plate is flat and unloaded. The cool-down from cure to room temperature is simulated by applying a ∆T = −180 • C. At room temperature, the plate is in its bistable region. The corner forces are applied progressively up to the snapping point, where the structure loses stability. Figure 2 shows the force-displacement diagram at ∆T = −180 • C. The curves therein are obtained with a Chebyshev-Gauss-Lobatto DQ grid consisting of 29 × 29 nodes, which guarantees accuracy. The Legendre polynomials' order increases up to convergence, as shown in Table 2. Legendre polynomials of seventh order reach convergence, with 64 degrees of freedom and a 1.8% difference with respect to the FE analysis. The solution for n = 6 (49 degrees of freedom) is closer to the FE result than that with 62 degrees of freedom in Pirrera et al. [7]. This result shows that the proposed model achieves greater accuracy, using fewer degrees of freedom, in comparison to [7]. This, still retaining all of the advantages of [7] with respect to conventional FE analyses.

Case-Study 2: Multi-mode Morphing
Eckstein et al. [11] studied the deformation of a thermally actuated cylindrical [0 4 /90 4 ] non-symmetric laminate. The combination of initial curvature and thermally-induced moments was found to give rise to multiple deformation modes. The laminate indeed shows two orthogonal bending modes at curing and room temperature, and an intermediate twisting mode.   The authors formulated a Rayleigh-Ritz analytical model assuming a second order polynomial approximation for the out-of-plane displacements, with a resulting constant curvature tensor as in Dano et al. [39]. The model was validated against a FE simulation consisting of 256 Abaqus S8R elements, using a uniformly distributed mesh. The FEM and the Rayleigh-Ritz model showed excellent agreement in predicting the two orthogonal cylindrical shapes. However, the analytical model overpredicted the values of curvature in the range of temperature where twisted configurations are present. The authors attributed this limitation to the simplicity of the shape functions used to describe the in-plane strains, which resulted in an overestimation of the membrane stiffness. This constraint makes it preferential for the laminate to store strain energy in the form of bending, hence the overprediction of the curvatures. The authors also observed that the approach adopted in Fernandes et al. [40] would likely yield more accurate results, because the membrane problem is, first, approximated with more degrees of freedom and, then, solved separately using FEM.
The model in [40], as that in [11], relies on the simplifying hypothesis of uniform curvature throughout the morphing process. As discussed in section 1, our work follows an approach similar to that of Fernandes et al. [40], in that the membrane problem is solved numerically and independently. However, differently from Eckstein et al. [11] and Fernandes et al. [40] the model proposed here does not rely on any restricting hypothesis on the curvature field. The multi-mode morphing presented in Eckstein et al. [11] is then used as a benchmark for the ability of our model to predict any kind of deformed configuration accurately by dividing the bending and stretching components of the total potential energy more precisely. In particular, in this section we show the benefits provided by the Legendre polynomials in capturing the membrane and the bending components compared to full polynomials. This choice results in further insight as we are able to predict bifurcation points and transitions from one deformation mode to the other accurately.
For the sake of simplicity, we have adapted the work done in [11] by considering material properties which are constant with temperature. These are reported in table 3. Shell dimensions and curvatures are unchanged from reference [11]. An FE analysis has been performed by using Abaqus. The shell is modelled using 8100 S4R elements with a total of 8281 nodes, which provide converged results. The cooling-down process has been simulated by using both Static-General and Riks steps with Nlgeom on. An imperfection in the form of a corner force has been introduced in order to drive the deformation towards the stable twisted configuration. The effect of the imperfection is negligible on linear behaviour because the shell assumes a cylindrical (non-twisted) configuration outside of the bifurcation region. For the nonlinear analyses to run successfully, it has been necessary to introduce numerical dissipation. The dissipation factor has been chosen to be sufficiently small so that the ratio between dissipated energy and total strain energy was smaller than 10 −3 [13]. Figure 3 shows the results given by the proposed model in terms of Legendre polynomials of increasing order. Convergence of the DQM grid is reached using a 31×31 Chebyshev-Gauss-Lobatto mesh. In particular, Figure 3 shows the displacement w of two adjacent corners (solid and dashed lines) of the shell as a function of the temperature. The shell is initially cylindrical in its stress free configuration at 180 • C (note that the solid and dashed lines are hence coincident). Upon cooling down, the shell approaches the first bifurcation point at 124 • C when the twisted configuration appears. The shell remains in its twisted configuration until the second bifurcation point at 54 • C, when a new cylindrical configuration appears, orthogonal to the initial one. By decreasing the temperature down to 0 • C, the shell keeps its cylindrical configuration and no further bifurcation points are present. Convergence has been reached by using sixth-order Legendre polynomials. These results are in excellent agreement with those given by FE Static-General runs. It is interesting to observe that the FE simulation performed using the Riks method was unable to capture the stable twisted configurations, despite the presence of the imperfection. Conversely, the analysis carried out using Static-General steps captures the twisted shape, although without apparent bifurcation points, as clear from the insets of Figure 3. Riks and Static-General analyses are then both necessary and complementary, because Static-General alone does not provide accurate solutions for values of temperature nearby the two bifurcation points. Figure 4 shows the difference between results from our model and FE analyses. The colour maps plotted onto deformed shells represent the norm of the error between the displacements obtained by using the Abaqus' Static-General solution and Legendre polynomials of order six, for the three characteristic configurations at 0 • C, 90 • C and 180 • C. The norm of the error on the displacement field is defined in percentage as The largest difference in terms of displacements is ∼ 0.22%, localised at the edges of the shell and any further increments of n do not lead to any significant improvement of results.
In conclusion, the shape functions used in Eckstein et al. [11] describe accurately the cylindrical shapes of multi-mode morphing panels. However, they are shown to be not sufficiently flexible to describe the twisted configurations. Hence, the overestimation of the values of curvature for the range of temperature where the twisted configurations are expected. Conversely, the use of higher-order Legendre polynomials in the present model allows us to describe the deformed shapes accurately for all values of temperature.

Case-Study 3: Tristable Shell
As a final benchmark for our model we now consider the work of Vidoli [15] on initially doubly-curved shells having orthotropic  As mentioned in section 1, Vidoli [15] investigates the tristability of thin shells, assuming constant, linearly and quadratically-varying curvatures to approximate the displacement field. As stated by the author, the aim of that study was to provide an efficient tool for the qualitative analysis and quick parametric design of multistable shells, with an intrinsic incapability to obtain highly accurate results.
The model is extremely efficient as it captures the three stable configurations using just two unknowns. However, to limit the number of unknowns and hence the computational cost, the author imposed symmetric displacements and purposely exploited some relationships between degrees of freedom, derived assuming that the boundary conditions on the bending moments are satisfied on average.
Vidoli also states: "To further enhance the model predictive capabilities, higher order expansions can be considered for the transverse displacement field, following the same reasoning lines. However, it is useful to remark that the proposed procedure, despite not being limited to polynomial ansatz only, does not scale efficiently when a large number of degrees of freedom has to be considered" [15].
Having shown in our first two case studies that by using DQM the procedure does scale efficiently and without the need to impose symmetries, we now intend to demonstrate that, in addition, overcoming this limitation also does not require the boundary conditions on the bending moments to be satisfied on average.
In order to do this, the proposed model is compared to the results in [15], reproduced numerically by FEA. The relative position in the energy landscape of the three stable equilibria can be shown with a stability plot similar to those shown in Vidoli and Maurini [18]. Such information is reported in figure 5 that is obtained using a second order Legendre approximation. This choice allows the energy to be plotted against the only two non vanishing Legendre parameters q 02 and q 20 , which play the role of constant curvatures in this case. The contour plot of the energy shows three points of local minima that correspond to the three stable equilibria. Knowing where the stable equilibria are is useful information because it drives the FE analysis as follows. Our FE model consists of 8100 S4R Abaqus elements with a total of 8281 nodes. The snapping process from the undeformed configuration to the two deformed, but stable, states is modelled separately with two different simulations. In the first simulation, the shell deforms from the initial undeformed configuration, denoted as stable state 1 (SS1), into the second stable configuration denoted as stable state 2 (SS2). In the second simulation, the shell transitions from SS1 to stable state 3 (SS3). In each of these analyses, we set three Static-General steps, with Nlgeom on. In the first step, the deformed shape (SS2 or SS3), as predicted by our model, is imposed over SS1 as an applied displacement field. In the second step, the forcing displacements are released at all nodes but the mid points on the four edges, thus letting the shell partially relax. In the third step, all applied displacements are removed, so that the shell is completely free to deform, reaching the closest stable equilibrium shape. This process is important as it avoids the need of introducing numerical damping. These FE results are hence highly accurate and reliable. Table 5 shows the displacement w of a corner of the shell for Legendre polynomials of different orders. Results are in excellent agreement with FE and convergence is reached for n = 5 with a 29×29 Chebyshev-Gaus-Lobatto mesh grid.
The three stable configurations obtained with fifth order Legendre approximations are shown in full in Figure 6. The contour plots represent the position error as per equation (38). Results are in excellent agreement, with a maximum difference of ∼ 0.8% at the corners of SS3.
These results show again that the efficiency of DQM allows the order of the polynomial approximation to be increased up to convergence, yet at relatively low computational cost. Furthermore, we have been able to reproduce Vidoli's results without introducing assumptions on the symmetry of the displacements and on the satisfaction on average of the bending boundary conditions. Table 5: Out-of-plane displacement w of a corner of the shell for the three stable configurations SS1, SS2 and SS3 as function of the order n of the Legendre polynomials.
[m] (c) Position error for SS3. Figure 6: Difference in terms of displacements between FE Static-General and the order 5 of the present model, for the three stable states SS1, SS2 and SS3.

Conclusions
By exploiting anisotropy of composite materials along with the inherent geometric nonlinearity of curved plates, these characteristically shell-like structures display unique and interesting behaviours, including multistability and temperature-triggered shape-morphing capability. At present, models predicting shell multistability present solution methods that are often a compromise between computational efficiency and result accuracy. In this paper, we presented an efficient, yet accurate, energy-based semi-analytical model to describe the multistability of thin anisotropic shells for morphing applications. Key steps of the proposed method are (1) the decoupling of the total potential energy into the stretching and bending contributions via the semi-inverse formulation of the constitutive equations and (2) the solution of the membrane problem by using DQM, which provides a high degree of accuracy at a relatively small computational cost.
We benchmarked the proposed model reproducing results available in the literature, complementing them and providing further insight into the numerical intricacies required to capture characteristic features of kinematically nonlinear morphing structures. The main features and advantages of our model in comparison to the state-of-the-art are summarised as follows.
Accuracy of results. By using a set of complete, orthogonal shape functions, in this case Legendre polynomials, to approximate transverse displacements, the values of the membrane and bending components of the total strain energy are captured with higher accuracy and fewer degrees of freedom, in comparison to previous works using non-orthogonal polynomial basis functions [7,11].
Symmetry and boundary conditions. By virtue of the DQM we are able to adopt a systematic approach for the solution of the governing equations that does not require any limiting assumption on the symmetry of the displacement field. Furthermore, no additional constraints on the boundary conditions are imposed. Despite this added complexity, the robustness of DQM facilitates converged solutions.
Computational efficiency. The membrane problem is, first, recast in terms of differential operators and then solved, once and for all, inverting one small (as compared with FE) sparse matrix of DQ weighting coefficients.
By virtue of the efficiency of DQM, higher-order polynomial approximations can be used up to convergence, yet keeping the computational cost at relatively low levels. The utility of the proposed model proves particularly useful when computing snap-through loads of bistable plates. We achieve greater accuracy, compared to the state-of-the-art [7], using fewer degrees of freedom.