Scholarly article on topic 'Morphing shell structures: A generalised modelling approach'

Morphing shell structures: A generalised modelling approach Academic research paper on "Materials engineering"

CC BY
0
0
Share paper
Academic journal
Composite Structures
OECD Field of science
Keywords
{"Morphing structures" / "Differential Quadrature Method" / "Semi-inverse formulation" / "Differential operator"}

Abstract of research paper on Materials engineering, author of scientific article — E. Lamacchia, A. Pirrera, I.V. Chenchiah, P.M. Weaver

Abstract 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 energy-based 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.

Academic research paper on topic "Morphing shell structures: A generalised modelling approach"

ELSEVIER

Contents lists available at ScienceDirect

Composite Structures

journal homepage: www.elsevier.com/locate/compstruct

Morphing shell structures: A generalised modelling approach

E. Lamacchia a,*> A. Pirreraa, I.V. Chenchiahb, P.M. Weaver3

a Advanced Composites Centre for Innovation and Science, University of Bristol, Queen's Building, University Walk, Bristol BS8 1TR, UK b School of Mathematics, University of Bristol, University Walk, Bristol BS8 1TW, UK

CrossMark

ARTICLE INFO

ABSTRACT

Article history:

Available online 4 July 2015

Keywords:

Morphing structures Differential Quadrature Method Semi-inverse formulation Differential operator

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 energy-based 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 bench-marked against three exemplar case studies taken from the literature.

© 2015 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://

creativecommons.org/licenses/by/4.0/).

1. 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 pre-stressed laminates [1-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

* Corresponding author. E-mail address: ettore.lamacchia@bristol.ac.uk (E. Lamacchia).

example, Hyer [1] and Hamamoto and Hyer [8] proposed energy-based models for bistable plates that assume laminates to be initially flat. Subsequently, Galletly and Guest [9] and Guest and Pellegrino [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 multi-mode 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,

http://dx.doi.org/10.1016/j.compstruct.2015.06.051 0263-8223/© 2015 The Authors. Published by Elsevier Ltd.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

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 and Pellegrino [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 aniso-tropic laminated shells possessing general planforms and initial curvatures. In this paper, following the work of Guest and Pellegrino [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 ani-sotropy 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 and Pellegrino [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 in-plane 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].

2. Methodology

2.1. 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 Foppl-von Karman formulation of the strain field e, defined as [25]:

e = e + zk, (!)

gu i 1 ¡gw\2 gx ^ 2 I gx)

gv i 1 (gw^' gy ^ 2 \ gy I

_ gv 2 \ gx

gw gw ' gy gx

and k =

~kxx~ g2w gx2

kyy — g2w ~w

.kxy. 2 g2w _ gxgy

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

- rTe dzdX,

-h/2 2

where, following Classic Laminate Analysis (CLA) and conventional notation, the stress tensor is defined as

r = Q(e - er),

where Q is the reduced stiffness matrix, X 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, Eq. (3) can be recast as

e T rA B e [Nr! T e

Ak B D Ak Mr Ak

e A* B* N + Nr"

M + Mr. -BT D* Ak

where A, B and D are the in-plane, coupling and bending stiffness matrices, Ak = k - k0 is the change in curvature from the refer-ence/undeformed configuration, having curvature k0, to the current/deformed configuration, and Nr and Mr are the in-plane stress and moment resultants due to non-mechanical 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

r -lTrfl * D* -irW"l\

where A* = A"1, B* = -A-1B, D* = D - BA-1B, and N = [Nxx Nyy Nxy ]T is the in-plane stress resultsnt and M = [Mxx Myy Mxy ]T is bending moments resultant.

Notably, Eq. (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 bending-stretching couplings embodied in D*.

In other words, we now have an energy formulation in which the terms NTA*N and AkTD*Ak 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

u N T A* 0 N Nr T A* B* N r ■

2 Ak 0 D* Ak 2Mr 0 I Ak

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.

2.2. 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.

2.2.1.1. 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 two-dimensional 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

= det k - det k0.

5y2 " dxdy

Eq. (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

L1 _ [ay7 aF 2âxy7] (9)

and using (6), the compatibility condition (8) becomes

L1 [A*(N + Nr) + B*Ak] = det k - det ko. (10)

2.2.1.2. In-plane equilibrium. The membrane equilibrium equations for a free-free plate [25] are

gx + gy — on X;

gy — 0; on X;

gNxy SN gx

on gx,

where dX is the boundary of the plate and n its normal. Introducing the differential operators

gx 0 0

0 0 gy

0 it, 0 gy

the equilibrium equations become (¿2 + ¿3)N — 0.

2.2.2. Formulation of the membrane problem

By combining Eqs. (10) and (13), we define the membrane problem as

det k — det k0 — ¿1 0

Nn = 0,

on X, on gx,

Note that when writing det k, k is intended in its 2 x 2 form, i.e. k =

where the differential operator L is defined as

L2 + L3

It is worth noting that Eq. (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 Eq. (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 Eq. (15), the differential operators L1, L2 and L3 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].

2.3. 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,

It should be noted that the size of L depends on the size of the DQgrid 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 qj 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 qj terms only.

Equilibrium configurations correspond to extrema of n(q,j), therefore satisfying the expression

Vn = 0;

which results in a set of nonlinear equilibrium equations of the kind

—P = 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

—p > 0, for all qj. (20)

3. 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:

x = Lxx, y = Lyy, u = Ud U, v = VdU, w = WdW, e = EU, Dk = KAfc, N = RN, where the dimensional scaling parameters

w(x, y) = Wo(x, y) + Y/qijPi(x)Pj(y),

i=0 j=0

where w0 represents the reference (undeformed) shell geometry, Pi(x) and Pj (y) are Legendre polynomials of order i and j defined as [32]

Pi(/) = £

l\i-l- 1V1 -

for l = i,j and / = x, y (17)

and here used as shape functions; and where, finally, qj 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 Pi(x)Pj(y) products and the q¡j parameters.

Assuming a displacement field via (16), the in-plane stress resultants are obtained as function of curvatures only by solving Eq. (14):

N = L—

det k — det k0 — L1 A*Nr — L1B*Dk 0

= L—1f,

with N ■ n = 0 on IX. Computationally, the boundary conditions are imposed by converting the expression N ■ n = 0 on IX 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.

2Lx and 2Ly are the side lengths of the plate along the Cartesian axes,

Ud, Vd and Wd are coefficients defined as [7,36,37]

1ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Ud = Fx 2A1lA22D1lD22, V d = i^WÄ

wd = 4 A11A22D11D22,

E and K are used to scale the in plane strains and curvatures, respectively, and are defined as

1 Wd 0

2 L2 0

2 l2 0

— W1 0

_ Wd — l2

LxLy ■

• R scales the membrane stress resultant N and it is defined as R = 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

Nr = Rn S and Mr = RM U. Following [7],

Npy r Zk+1 -Rn — £/ (lk(T)Xk(T)AT 0 dz,

k—1 zk

Nil* fZk+1 _ Rm — Y Qk (T )«k(T )AT0zdz,

k—1 zk

where AT0 — T - Tref is the thermal gradient from cure temperature, ~ is a non-dimensional parameter defined from AT0s — T - Tref and it is noted that, in general, Q.k(T) and ak(T), the coefficients of thermal expansion, may both be functions of the temperature.

Substituting the parameters defined in Eqs. (21) and (25) into Eq. (7), one obtains the non-dimensional form of the total strain energy

n — / i (1ntA*N+1 AkD*Ak -^TÀA*! - xb*Ak- xd*Ak) dxdy,

where the non-dimensional matrices are defined as

A*— ^ RT A*R — ^ ETAE, nd nd

A — lxTV RN a rN ,

b* = D *

nd — tr

lxly rTb*k,

LxLy KD*K, d*— P RM K, nd

E 0 0 K

A B B D

E 0 0 K

is the parameter, defined as in [7], used to scale the total strain energy, so that n = nd 1!.

For completeness, the membrane problem must be scaled and recast in terms of non-dimensional parameters too. Hence, the compatibility condition (8) becomes

1 g2èxx 1 g2è„

2 g~2 2 g~2

and by defining

^ — kjiyy - 4% - (~0xx~0yy - 4k2x^,

¿! =1 h JL sL -4S— ]

1 2 L sv2 a*ak J '

B* = R BK,

and substituting the non-dimensional form ofEqs. (6), (29) is recast as

LW JV - B*Ak + R-1Rn~ ) — kxxkyy - 4fc?v - ( k0xxk0yy - 4E

Similarly, the equilibrium conditions (13) become

^ri (cL2R + 1L3R)N —

Finally, combining (32) and (33), the membrane problem in non-dimensional form is given by

(LiV — /, on X e [-1,1]x[-1,1], \n • n — 0, on gx, where

^t;^ (Lx ^ 2 + LyL 3)R

.V —

L1B *k

kxxkyy - 4k2y -0

k0xxk0yy - 4k2„, I

The solution of the membrane problem is simply given by JV = L,

which is a function of x, and of the Lagrange parameters qij.

Analogously to Section 2.3, by substituting Eq. (37) into (26) we obtain the total potential energy as a function of qij terms only. The equilibrium positions are obtained by minimising Eq. (26) with respect to the Lagrange parameters.

It is interesting to observe that the operator expressed in Eq. (30) is a second order differential operator. Hence, in Eq. (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 quadratically-distributed thermal field applied to the cylindrical shell depicted in Fig. 1a. For the sake of simplicity, and to avoid any membrane-bending coupling due to anisotropic material properties, we consider an isotropic aluminium (Al 6061-T6) shell, having in plane dimensions 1 x 1m 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.

Fig. 1b and c show the deformations of the two orthogonal cross sections of the cylindrical shell, which are represented by the blue and red lines in Fig. 1a. In particular, Fig. 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 Fig. 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 non-mechanical stress if the external field is at least of the second order in x or y.

4. 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.

4.1. 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 unsym-metric 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.

Fig. 1. Fig. 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. Fig. 1b and c show respectively the displacements of the initially curved cross section (blue line in Fig. 1a) and of the initially flat cross section (red line in Fig. 1a) due to the thermal fields. From Fig. 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 Fig. 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. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

-6 -4 -2 0 2 4 6

Corner displacement w [m] -10-2

Fig. 2. Force-displacement diagram at various orders of the Legendre polynomial, compared to the results in reference [7]. The comparison in the proximity of the snapping load is shown in the insets. Solid and dotted lines represent the stable and the unstable solutions, respectively. The applied thermal gradient is AT = -180 °C.

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.

Table 1

Material properties for the unsymmetric [02/902] laminate.

E11 [Gpa] E22 [Gpa] G,2 [Gpa] vn [-] a„ [10-6 /°q «22 [10-6 /°C] 161 11.38 5.17 0.32 -1.8 3

Table 2

Snapping load at various orders of Legendre polynomial compared with the results of reference [7] (Error = L°'dL'„'adLE°

: 100).

Polynomial order Number of unknowns Load [N] Error [%]

Legendre 2 9 4.98 48.6

Legendre 3 16 4.32 28.9

Legendre 4 25 3.97 18.5

Legendre 5 36 3.65 8.96

Legendre 6 49 3.42 2.1

Legendre 7 64 3.41 1.8

Complete 11 (symmetric), Ref. [7] 62 3.51 5

FE, Ref. [7] 15606 3.35 -

Table 3

Material properties for the unsymmetric [04/904] laminate.

En [GPa] E22 [GPa] Gn [GPa] vn [-] aii [10-6/°C] «22 [10-6/°C]

129.55 8.85 5.28 0.33 -2.26 23.4

Fig. 2 shows the snap-through load for the [02/902] 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 x 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, Eq. (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 DT = -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. Fig. 2 shows the force-displacement diagram at DT = -180°C. The curves therein are obtained with a Chebyshev-Gauss-Lobatto DQ grid consisting of 29 x 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.

4.2. Case-study 2: multi-mode morphing

Eckstein et al. [11] studied the deformation of a thermally actuated cylindrical [04/904] 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 and Hyer [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].

Fig. 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 x 31 Chebyshev-Gauss-Lobatto mesh. In particular, Fig. 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

•10"2

50 100

Temperature [°C]

Fig. 3. Displacements w of two adjacent corners of the cylindrical shell as function of temperature. The solid and dashed line in the figure correspond to any two adjacent corners of the shell. From cure temperature to the first bifurcation point at 124 °C, the shell has a cylindrical configuration and any two adjacent corners displays the same out-of-plane displacement. In the region between the first bifurcation point and the second bifurcation point located at 54 °C, the shell assumes a twisted configuration, hence two adjacent points display different out-of-plane displacements. By decreasing the temperature below the second bifurcation point, a new cylindrical configuration appears, orthogonal to the initial one. The figure shows the results given by various order of the Legendre polynomials compared to the FE analysis, and their convergence. Results provided by using the sixth order of the Legendre polynomials (n — 6) are in excellent agreement with the FE, and any further increment of n does not lead to any significant improvement.

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 Fig. 3. Riks and Static-General analyses are then both necessary and complementary, because

Static-Generai alone does not provide accurate solutions for values of temperature nearby the two bifurcation points.

Fig. 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

Error [%] 0.20

(a) Position error for the deformed shape (b) Position error for the deformed shape at 0°C at 90 °C

-o.i -o.i

(c) Reference shape at 180 °C

Fig. 4. Difference in w between FE Static-General and the order 6 of the present model, for the three characteristic configurations at 0 °C,90 °C and 180 °C.

(ULeg - UFE)2 + (VLeg - Vfe)2 + (WLeg - WFE)2

100. (38)

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.

4.3. 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 layup [45, -45, -45,45, -45,45,45, -45]. These laminates are of interest because they may feature three stable shapes. The panel studied here has free-free boundary conditions, is clamped in the middle and has a 0.25 x 0.25 m planform. Material properties are reported in Table 4.

As mentioned in Section 1, Vidoli [15] investigates the tristabil-ity 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

Table 4

Material properties for the [45, -45, -45,45, -45,45,45, -45] laminate.

E11 [GPa]

E22 [GPa]

G12 [GPa]

tpiy [mm]

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 Fig. 5 that is obtained using a second order Legendre approximation. This choice allows

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.

Polynomial order

SS1 w [m]

SS2 w [m]

SS3 w [m]

-0.1148 -0.1148 -0.1148 -0.1148 -0.1148

-0.08721 -0.08725 -0.08730 -0.08731 -0.08732

0.05533 0.05591 0.05623 0.05624 0.05671

-0.04 -0.02

0.02 0.06

-0.1 -0.05 0 0.05 0.1 %2

Energy [J]

Fig. 5. Contour plot of the total strain energy as function of the curvatures q02 and q20. The insets show the three points of local minima, corresponding to the three stable states SS1, SS2 and SS3.

the energy to be plotted against the only two non vanishing Legendre parameters q02 and q20, 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 x 29 Chebyshev-Gaus-Lobatto mesh grid.

The three stable configurations obtained with fifth order Legendre approximations are shown in full in Fig. 6. The contour plots represent the position error as per Eq. (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.

5. 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 aniso-tropic 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

Error [%]

x [ra]

(b) Position error for SS2.

(c) Position error for SS3.

Fig. 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.

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.

Acknowledgements

This work was supported by the Engineering and Physical Sciences Research Council through the EPSRC Centre for Doctoral Training in Advanced Composites for Innovation and Science [Grant No. EP/G036772/1].

References

[1] Hyer MW. Some observations on the cured shape of thin unsymmetric laminates. J Compos Mater 1981;15:175-94.

[2] Daynes S, Potter KD, Weaver PM. Bistable prestressed buckled laminates. Compos Sci Technol 2008;68(15-16):3431-7.

[3] Thill C, Etches J, Bond IP, Potter KD, Weaver PM. Morphing skins. Aeronaut J 2008:117-39.

[4] Diaconu CG, Weaver PM, Mattioni F. Concepts for morphing airfoil sections using bi-stable laminated composite structures. Thin-Walled Struct 2008;46:689-701.

[5] Daynes S, Weaver PM. A morphing composite air inlet with multiple stable shapes. J Intell Mater Syst Struct 2011;22(9):961-73.

[6] Lachenal X, Weaver PM, Daynes S. Multistable composite twisting structure for morphing applications. Proc R Soc A 2012;468:1230-51.

[7] Pirrera A, Avitabile D, Weaver PM. Bistable plates for morphing structures: a refined analytical approach with high-order polynomials. Int J Solids Struct 2010;47(25-26):3412-25.

[8] Hamamoto A, Hyer M. Non-linear temperature-curvature relationships for unsymmetric graphite-epoxy laminates. Int J Solids Struct 1987;23:919-35.

[9] Galletly DA, Guest SD. Bistable composite slit tubes. ii. a shell model. Int J Solids Struct 2004;462:4503-16.

[10] Guest SD, Pellegrino S. Analytical models for bistable cylindrical shells. Proc R Soc A: Math Phys Eng Sci 2006;462:839-54.

[11] Eckstein E, Pirrera A, Weaver PM. Multi-mode morphing using initially curved composite plates. Compos Struct 2014;109:240-5.

[12] Coburn BH, Pirrera A, Weaver PM, Vidoli S. Tristability of an orthotropic doubly curved shell. Compos Struct 2013;96:446-54.

[13] Pirrera A, Avitabile D, Weaver PM. On the thermally induced bistability of composite cylindrical shells for morphing structures. Int J Solids Struct 2012;49:685-700.

[14] Seffen KA. Morphing bistable orthotropic elliptical shallow shells. Proc R Soc A: Math Phys Eng Sci 2007;463:67-83.

[15] Vidoli S. Discrete approximation of the Foppl-Von Karman shell model: from coarse to more refined models. Int J Solids Struct 2013;50:1241-52.

[16] Mansfield EH, The bending and stretching of plates. 2nd ed. Cambridge University Press; 1989.

[17] Mansfield EH. Bending, buckling and curling of a heated thin plate. Proc R Soc A: Math Phys Eng Sci 1962;268:316-27.

[18] Vidoli S, Maurini C. Tristability of thin orthotropic shells with uniform initial curvature. Proc R Soc A: Math Phys Eng Sci 2008;464:2949-66.

[19] Shu C, Differential quadrature and its application in engineering. Springer; 2000.

[20] Bellman R, Casti J. Differential quadrature and long-term integration. J Math Anal Appl 1971;34(2):235-8.

[21] Striz AG, Chen W, Bert CW, Static analysis of structures by the quadrature element method (qem), International Journal of Solids and Structures 1994;2807(31).

[22] Mattioni F, Weaver PM, Friswell MI. Multistable composite plates with piecewise variation of lay-up in the planform. Int J Solids Struct 2009;46(1):151-64.

[23] Aimmanee S, Hyer MW. Analysis of the manufactured shape of rectangular thunder-type actuators. Smart Mater Struct 2004;13(6):1389-406.

[24] Koiter WT. On the nonlinear theory of thin elastic shells. Proc. Koninklijke Nederlandse Akademie van Wetenschappen 1966;38(B69):1-54.

[25] ReddyJN, Mechanics of laminated composite plates and shells. 2nd ed. CRC Press; 2004.

[26] Love AEH, A treatise on the mathematical theory of elasticity. 4th ed. Cambridge. 1927. p. 173.

[27] Novozhilov VV, The theory of thin shells - A translation by Lowe P.G. and Radok J.R.M., P. Noordhoff LTD 1 1959. p. 27-29.

[28] McConnel AJ. Application of the absolute differential calculus, 1. Blackie and Son Limited; 1931 [p. 193-207].

[29] Lamacchia E, Pirrera A, Chenchiah IV, Weaver PM. Non-axisymmetric bending of thin annular plates due to circumferentially distributed moments. Int J Solids Struct 2014;51(3-4):622-32.

[30] Runge C. Über empirische funktionen und die interpolation zwischen äquidistanten ordinaten. Zeitschrift für Mathematik und Physik 1902;46:224-43.

[31] Hille E. Analytic function theory, vol. I. AMS Chelsea Publishing; 1962.

[32] Koepf W, Hypergeometric summation: an algorithmic approach to summation and special function identities, Braunschweig, Germany: Vieweg; 2014.

[33] Wu Z, Raju G, Weaver PM. Comparison of variational, differential quadrature, and approximate closed-form solution methods for buckling of highly flexurally anisotropic laminates. J Eng Mech 2012;139:1073-83.

[34] Moore EH. On the reciprocal of the general algebraic matrix. Bull Am Math Soc 1920;26(9):394-5.

[35] Penrose R. A generalized inverse for matrices. Proc Cambridge Philos Soc 1955;51:406-13.

[36] Diaconu CG, Weaver PM. Postbuckling of long unsymmetrically laminated composite plates under axial compression. Int J Solids Struct 2006;43(22-23):6978-97.

[37] Nemeth MP, Nondimensional parameters and equations for nonliner and bifurcation analyses of thin anisotropic quasi-shallow shells, NASA/TP-2010-216726. 2010. p.589-606.

[38] Diaconu CG, Weaver PM, Arietta AF. Dynamic analysis of bi-stable composite plates. J Sound Vib 2009;322:987-1004.

[39] Dano M-L, Hyer MW. Thermally-induced deformation behavior of unsymmetric laminates. Int J Solids Struct 1998;35:2101-20.

[40] Fernandes A, Maurini C, Vidoli S. Multiparameter actuation for shape control of bistable composite plates. Int J Solids Struct 2010;47:1449-58.