Initial post-buckling of variable-stiffness curved panels

Abstract Variable-stiffness shells are curved composite structures in which the fibre-reinforcement follow curvilinear paths in space. Having a wider design space than traditional composite shells, they have the potential to improve a wide variety of weight-critical structures. In this paper, a new method for computing the initial post-buckling response of variable-stiffness cylindrical panels is presented, based on the differential quadrature method. Integro-differential governing and boundary equations governing the problem, derived with Koiter׳s theory ( Koiter, 1945 ), are solved using a mixed generalised differential quadrature (GDQ) and integral quadrature (GIQ) approach. The post-buckling behaviour is determined on the basis of a quadratic expansion of the displacement fields. Orthogonality of the mode-shapes in the expansion series is ensured by a novel use of the Moore–Penrose generalised matrix inverse for solving the GDQ–GIQ equations. The new formulation is validated against benchmark analytical post-buckling results for constant stiffness plates and shells, and compared with non-linear finite-element (FE) analysis for variable-stiffness shells. Stability estimates are found to be in good agreement with incremental FE results in the vicinity of the buckling load, requiring only a fraction of the number of variables used by the current method.


Introduction
Curved panels are used extensively in modern aerospace, civil, mechanical and naval engineering structures. Their excellent structural efficiency, based on high buckling loads and low mass, make them attractive solutions for lightweight structures which develop large compressive stresses. In order to take their full advantage, non-linear aspects of their behaviour are an important consideration in the design phase. In compression they often have unstable buckling modes, causing unwanted jumps in load. Furthermore, one can expect to see buckling occurring well below the calculated critical load under axial compression (Cox and Clenshaw, 1941;Jackson and Hall, 1947;Peterson, 1969;Thornburgh and Hilburger, transition region exists in the design space, where stable post-buckling stops and unstable post-buckling begins. The numerical results of this study were confirmed by Pope (1965Pope ( , 1968) who used the Principle of Virtual Work to determine approximate equilibrium paths for panels with finite-displacements and fixed mode-shapes. Fig. 1 shows Pope's results for an isotropic panel. Based on these curves it may be concluded that the effect of the curvature is to (a) significantly increase the buckling load and (b) reduce the tangent stiffness at the buckling load and hence, the stability. In practice, a panel under displacement controlled edge loading, exhibits dynamic jumping from the unstable bifurcation to the nearest stable state. This effect results in a significant drop in load.
Studies into the post-buckling behaviour of composite curved panels are much less numerous. Bauld and Khot (1982) used a Rayliegh-Ritz and finite-difference approach to compute equilibrium paths for ½0; 90 2s laminates. This was done in an attempt to replicate experimental results, but with relatively little success. Zhang and Matthews (1983) investigated the effect of a small family of stacking sequences using a non-linear Rayleigh-Ritz Fourier series formulation. Within the group of laminates analysed, bend-twist coupling was found to have a significant effect on the drop in load following bifurcation. Madenci and Barut (1994) investigated composite panels under compression using non-linear finite-element analysis, but were primarily interested in the effect of cut-outs rather than stacking sequences. They did, however, find extremely good agreement with previous experiments. More recent publications, such as Hilburger et al. (1999Hilburger et al. ( , 2004, have also focused on the effect of boundary conditions and cut-outs. There have, therefore, been no firm conclusions made with regard to the general effects of the stiffness constants on the post-buckling behaviour of curved panels (presumably due to the large size of the design space). If these effects are to be explored more fully, the development of efficient design tools is essential.
The design space of these structures may be increased further by assuming that the shell is manufactured with variableangle tows (VAT) and has variable stiffness-constants. One of the first papers to introduce this concept was concerned with the optimisation of perforated composite plates for increased buckling capacity (Hyer and Lee, 1991). It was demonstrated that by forcing well supported regions of the structure to be under the highest stresses, the maximum buckling load could be increased considerably over an optimal constant-stiffness structure. It was later shown by Gürdal et al. (2008) that the use of VAT may not only increase the buckling load, but the conflicting requirements for both high in-plane stiffness and high buckling loads may be resolved. In addition to buckling optimisation, it has also been shown that VAT laminates can be used to delay first ply failure of plates under compression (Lopes et al., 2007); tune the fundamental frequency of shells (Blom et al., 2008) and increase global structural stiffness of cylindrical shells (Wu, 2008). Furthermore, advances in manufacturing techniques, such as the development of the continuous tow shearing technique by Kim et al. (2012), are making VAT designs a practical possibility. Current investigations into the post-buckling behavior of VAT structures have focused on flat plates (Rahman et al., 2011;Raju et al., 2013;Wu et al., 2013aWu et al., , 2013b. The development of closed form solutions for these structures is hindered by the generality of the variable laminate and the complexities of the shell governing equations. Numerical methods are, therefore, the only practically viable option for computing the dependent field variables and the stability characteristics. In the context of constant-stiffness curved panels, previous authors have mostly applied incremental procedures to the post-buckling problem, such as the Newton-Raphson or Riks algorithms (Panda and Ramachandra, 2010;Thornburgh and Hilburger, 2006). Such methods are, of course, useful when a 'deep' non-linear solution is required, they do however, produce large quantities of data, and have relatively large computational requirements: two characteristics which are not well suited to optimisation problems. Asymptotic methods, on the other hand, have the advantage that they may project an approximation of the non-linear solution from a point of interest (i.e. a bifurcation point) and hence only require a small number of variables. They are useful, therefore, for estimating the behaviour of a system without actually following an equilibrium path. Previous implementations of asymptotic methods for shell problems have mostly been produced within the finite-element framework and for general structures (Casciaro et al., 1998;Garcea et al., 2009;Kheyrkhahan and Peek, 1999 Pope (1965)].
In this work, we explore the use of Koiter's original theory to compute approximate equilibrium paths specifically for variable-stiffness cylindrical panels. The Differential Quadrature (DQ) method has been used as a solution scheme for determining pre-buckling, buckling and second-order displacement fields. It has been found by previous investigations that the DQ method is an efficient means of solving thin-plate equilibrium equations where variable stiffness constants occur (Raju et al., 2012). The method has, therefore, been used here in conjunction with asymptotic expansion of the non-linear problem in an effort to minimise the total number of variables in the model.
The governing equations of various displacement vector fields are converted from their algebraic form into DQ form (i.e. matrix form), along with their boundary conditions and further constraints. Equations for the equilibrium paths are then derived in terms of the discrete displacement vectors for fast calculation. We compare the numerical solutions to benchmark algebraic results for orthotropic plates and isotropic panels. Finally, numerical solutions for a small group of representative variable-stiffness shells are compared with results from finite-element based Riks analyses.

Theory
In the present model, we consider a cylindrical panel with a VAT laminate as shown in Fig. 2. It is assumed that the tow orientations of the kth ply are given by a general function θ k ðx; yÞ and that the laminate has no significant coupling between membrane and bending deformation (i.e. B ¼ 0Þ. The panel has length l x , width l y , thickness h and a radius of curvature R. The variation of the reinforcement direction results in a laminate with the following constitutive properties: where A and D are the membrane and bending stiffness matrices, and and M ¼ ½M x M y M xy T are the membrane strains, changes in curvature, and membrane and bending stress resultants, respectively. As an approximation, shear deformation of the shell's section in the xz and yz planes are assumed to be negligible (i.e. plane sections remain plane (Kirchhoff, 1850)). The laminate is, therefore, represented by an equivalent single layer.
In our model we have adopted Donnell's (1933) shallow-shell kinematic equations. We write the strains and curvatures as The vector u ¼ ½u v w T contains the mid-surface displacement fields corresponding to the x, y and z directions. Eq.
(2) has been derived assuming that all terms involving 1=R 2 and quadratic terms involving the first derivatives of u and v are negligible (Donnell, 1933). It is explicitly assumed, therefore, that the shell is shallow, and the displacements in the zdirection are larger than those in the x and y directions. Fig. 2. The panel geometry with the current co-ordinate system.

Governing equations
Non-linear solutions for the displacement fields u, v and w are determined approximately using Koiter's (1945) method. In this scheme the loading parameter and the displacements are expanded in series of the following form: where a is a perturbation parameter. The fields u 0 , u 1 and u 2 are the pre-buckling and linear and quadratic buckling modes.
The coefficients λ c , λ 1 and λ 2 are the critical buckling load, the slope and curvature of the equilibrium path, respectively. Prior to buckling, the shell's displacements are assumed to increase linearly in proportion to a loading parameter λ; thus u ¼ u 0 ¼ λû 0 . At some critical load value, λ c , the shell bifurcates with a corresponding linear mode-shape u 1 . In the analysis, the stability of the linear mode-shape is determined by computing the second-order field u 2 such that the shell is in equilibrium for a small value of a. Evaluating the changes in energy due to a change in a (based on the solutions for u 1 and u 2 ), it is possible to predict the tangent of the load-deflection path and hence the post-buckling stiffness of the structure.
The governing equations forû 0 , u 1 and u 2 have been derived in A.1 using the principle of stationary potential energy. Introducing the following notation for the stress-resultants: with m A f0; 1; 2g, we write Pre-buckling: Buckling: Post-buckling: where A is the shell's area. In the literature Eqs. (5)-(7) are usually called Donnell's shallow-shell equations. Eq. (8) is not an equilibrium equation, rather a functional expression, which is used to enforce orthogonality of the linear and quadratic buckling modes u 1 and u 2 . It is required to ensure uniqueness of the solution to u 2 (Koiter, 1945).

Boundary conditions
All external forces in the model are applied at the boundary (see Fig. 3). The edges corresponding to x ¼ f0; l x g are termed ends and given the symbols Φ 1 x and Φ 2 x . Similarly, the edges corresponding to y ¼ f0; l y g are termed sides and have the symbols Φ 1 y and Φ 2 y . Bending boundary conditions are identical for each of the displacement fields u 0 , u 1 and u 2 (due to the linearity of the changes in curvature k). They are written as where m is the state (i.e. 0 pre-buckling, 1 buckling and 2 post-buckling). The fields Q x and Q y are the transverse shear stress resultants (for a shallow-shell) and are given by the expressions: In our model the following (optional) membrane constraints are also common to u 0 , u 1 and u 2 , they are The displacement functions U k x and V k x , and U l y and V l y , are the fixed axial and circumferential displacements at the ends and sides, respectively. The scalars u k xm , v k xm , u l ym , and v l ym are the axial and circumferential displacements of control points at the ends and the sides in each state m (see Fig. 3). In writing the edge displacements as they are in Eq. (10) we have reduced each sub-boundary to a discrete point having two degrees of freedom.

Pre-buckling membrane boundary conditions
External forces are introduced into the model in the pre-buckling state and are all multiplied by the loading factor λ.

Fixed stress boundary conditions are
where F k x and S k y , and S l x and F l y are prescribed stress-resultants at the ends and the sides. Eqs. (10) and (11) are mutually exclusive (at each edge), that is, one may either prescribe known displacements or known loads.
In addition to the standard boundary conditions given in Eq. (11), we have also incorporated force boundary conditions into our formulation. These have been implemented in a similar manner to multi-point constraints (MPCs) found in finiteelement formulations. Forces are applied to the edge control points, each of which has two degrees of freedom. We write the

Control point
Sides Φ l y : where the variables f k x , p k x , f l y and p l y are the concentrated external forces at the end and the side control points (see Fig. 3). Using Eqs. (10) and (12) one can enforce known forces on edges with unknown stress distributions. These constraints are useful in variable-stiffness shell problems since the non-uniform stiffness results in non-uniform stress-resultant distributions.

Buckling and post-buckling membrane boundary conditions
Due to the non-linearity in the strains ε, the membrane force boundary conditions are slightly different for each of the displacement fields u 0 , u 1 and u 2 . At edges where standard stress boundary conditions (11) have been applied, we have the following Buckling: Post-buckling: Alternatively, at edges where force boundary conditions (12) have been applied, we have Buckling: Sides Φ l y : Post-buckling: Sides Φ l y : Eqs. (15) and (16) must be used in conjunction with Eq. (10).

Path coefficients
Assuming that solutions for u 0 , u 1 and u 2 satisfying the governing equations (5)-(7) and the boundary conditions (9)-(16) have been found, it is possible to compute the change in load factor λ for a given change in the perturbation parameter a. Expanding λ as a Taylor series in terms of a gives Eq. (3).
The path parameters for the variable stiffness shell have been derived in A.2 and are given by the expressions: The stability of a shell in its initial post-buckling range may be determined from the path parameters. The slope indicates asymmetry of λ about the point a¼0. A shell with a non-zero slope must have at least one unstable branch originating from the bifurcation point. The curvature indicates the asymptotic stability of the adjacent state for both positive and negative values of a. In earlier works, it has been represented by the so-called b-imperfection sensitivity parameter (b ¼ λ 2 =λ c ), due to its relationship with buckling load imperfection sensitivity. An excellent explanation of the dependency of stability on the path coefficients can be found in El Naschie's (1991) text.

Differential quadrature implementation
In this paper a mixed differential-integral quadrature (GDQ-GIQ) technique has been used to solve the integrodifferential governing and boundary equations for the problems (5)-(16). The benefit of this is that the orthogonality (8) and force (12) constraints, as well as the equilibrium equations, may be satisfied in a single matrix inversion operation.
In differential quadrature (DQ) the gradient of a discretised function is assumed to be equal to a weighted sum of the function at all discrete values. The standard DQ approximation (Bellman et al., 1972) for a function f(x) is written in matrix form as where the matrix T x is a weighting coefficient matrix, f is a vector of grid-point variables, and N is the number of grid points in the x direction. Higher-order derivative weighting arrays are derived on the basis of Eq. (18) (i.e. by matrix multiplication) or by recursive formulae (Shu, 2000). Numerical values of the weighting coefficients T x ij are dependent on the DQ basis functions and the spacing of the grid points. In structural stability problems it has been established that using Lagrange polynomials as basis functions and a Chebyshev grid spacing yields fast convergence of both the statics and buckling problems (Sherbourne and Pandey, 1991;Wang et al., 2003;Raju et al., 2012).
In the present problem, the dependent fields are functions of two variables. We, therefore, discretise both axes, resulting in an N Â M grid of points ðx i ; y j Þ and a two-dimensional function-array f ðx i ; y j Þ ¼ f ij . As a matter of computational convenience it is useful to hold the discretised function f ij as a vector rather than as an array; thus permitting the use of standard matrix multiplication and inversion. We therefore extend the definition of f into a second dimension by writing it as where it now has dimensions ðNM Â 1Þ (Zeng and Bert, 2001). The numbering system of the elements in f is shown graphically in Fig. 4. Based on Eq. (18) and the relationship between f ij and f it is possible to write the partial derivatives of the twodimensional field f ðx; yÞ in DQ form as where the matrices F q and F qr both have dimensions ðNM Â NMÞ, and the dummy variables q and r are assigned to a particular co-ordinate direction x or y.
In a similar manner to GDQ, it is possible to replace the integration operations in the discretised equations by a matrix multiplication operation. In contrast to Gaussian quadrature, in which the grid spacing and the weighting coefficients are both unknown variables, GIQ uses a pre-defined grid. In this work, we use the Penrose's (1955) pseudoinverse operator to compute these integral weightings. This approach has not been reported previously and is an alternative to the integral weighting coefficient formulae developed in Shu (2000). Defining an indefinite-integral we write the discrete indefinite-integral of f in vector form as This is computed approximately by the expression: where ð Þ þ is the pseudoinverse operator and t is a vector of integration constants. We define the new integral weighting matrix in the q direction as where I q has the same dimensions as T q . The use of the pseudoinverse operator is required because T q is singular (reflecting the fact that constant terms are consumed by differentiation). On the basis of Eq. (23) we define a definite-integral weighting vector which gives the integral of the functional over a single grid direction, hence which have dimensions ð1 Â NÞ and ð1 Â MÞ, respectively. The vectors j x and j y are used to compute the definite integral of f along lines of constant y or x respectively. Noting that the function f is two-dimensional and that f has a structure given by Eq. (19), the elements of vectors j x and j y are cast into the matrices J x and J y , such that where J x and J y have dimensions ðM Â NMÞ and ðN Â NMÞ, respectively. Multiplying the vector of function values f by the matrix J x results in an ðM Â 1Þ vector of definite x-direction integrals of f at different co-ordinates y j . Likewise, multiplication of f by J y yields a vector of y-direction integrals at different co-ordinates x i .

Variable discretisation
Making use of the field discretisation in Eq. (19), we represent the displacement fields u, v and w at the DQ grid points by the ðNM Â 1Þ column vectors u, v and w, respectively. We do the same for the linear strains e and curvatures k, and the stress-resultants N and M, giving their DQ counterparts lower-case boldface symbols. Substituting the DQ approximations (20) into the equations of kinematics (2) where I is an ðNM Â NMÞ identity matrix, E and K are the strain and curvature co-efficient matrices, and Substituting Eq. (26) into Eq. (1), and noting that the membrane and bending stiffness coefficients are functions of x and y, the stresses and moments at the grid points are and m x m y m xy 2 6 4 where ðN x ; N y ; N xy Þ and ðM x ; M y ; M xy Þ are the membrane and bending stress-resultant coefficient matrices. They each have dimensions ðNM Â 3NMÞ.

Common boundary conditions
Bending boundary conditions (9) and kinematic membrane boundary conditions (10) are common to each of the vector fields u 0 , u 1 and u 2 . The bending boundary conditions are expressed as and and where the Φ superscript represents rows only corresponding to points on the relevant boundary. Collecting the coefficients of U from Eqs. (29)a to h, we may assemble a system of equations: where B B is the bending boundary coefficient matrix which has dimensions ðð2N þ 2MÞ Â NMÞ. Next, by sampling the prescribed membrane boundary displacements U k x , V k x , U l y and V l y at the grid-points and equating them to the difference between the slave and control point displacements, we write the following system of discrete equations: where C is the membrane kinematic constraint matrix and V is the constraint displacement vector. The dimensions of C depend on the number of displacement boundary constraints which have been imposed. They may, therefore, be anywhere between ð0 Â NMÞ and ðð2N þ 2MÞ Â NMÞ.
3.2. Linear algebra problem forÛ 0 Substituting the discrete strains and curvatures from Eqs. (26a,b) into the discrete stresses and moments from Eq. (28), and applying further numerical differentiation, we write the equilibrium equations (5) as and where Eq. (32) is the membrane equilibrium equation, Eq. (33) is the bending equilibrium in the equation, andÛ 0 is the prebuckling displacement vector. Combining Eqs. (32) and (33) and collecting the various coefficients ofÛ 0 we may write the following system of linear equations: where the matrices L ij represent the contribution of the displacement field j to the equilibrium of the shell in the i direction. This method of developing the field equations for variable-stiffness shells differs from other DQ implementations in that matrix manipulations have been used to express them in terms of the dependent variables. In a variable stiffness context the advantage of this technique is that the stiffness coefficient matrices are substituted directly into L, thus avoiding the need to evaluate their differentials. The membrane boundary conditions (11) are converted into DQ form by sampling the functions F k x , S k x , F l y and S l y at the boundary grid-points, and equating the values to the relevant rows of Eq. (28a). These equations are used to build the system of equations: where B M is the membrane boundary coefficient matrix and F is a vector of the prescribed boundary stresses. Where force boundary conditions have been used, Eq. (12) applies, and the DQ equations for the control point are Sides Φ l y : j xT ○ðN y Þ Φ ÁÛ 0 ¼ f l y and ð36cÞ The integration vectors j x and j y have been used to integrate the stresses along each edge based on the global displacement vectors. Coefficients of U 0 0 in Eq. (36) are collected and assembled into B M . The forces f k x , p k x , f l y and p l y are contained in the vector F 0 .
The solution of the pre-buckling problem is now a solution of the over-determined linear algebraic system of Eqs. (34), (29), (31) and (35). In this work the matrix pseudoinverse operator ðÞ þ (Penrose, 1955) has been used to solve forÛ 0 , hencê The advantage of this approach is that all of the DQ governing and constraint equations are satisfied simultaneously, without the need for any modification to the stiffness matrix. It has been found to be a robust method of solving the shell static problem, finding solutions even when equilibrated rigid-body modes exist, for example when p 1 y ¼ p 2 3.3. Linear algebra problem for U 1 The left-hand sides of Eqs. (6a-c) have the same form as those in the pre-buckling equations. We may, therefore, use L as the stiffness matrix for the buckling problem. On the right-hand side of Eq. (6a), a pressure term exists which is linear in w 1 . This is the basis of the geometric stiffness matrix. Converting this term into DQ form and substituting the stresses from the pre-buckling problem we have where G zw is the so-called geometric stiffness matrix. The unconstrained eigenvalue problem for U 1 is, therefore, given by where λ c is the critical buckling load. The array G has only one non-zero sub-matrix G zw . Eq. (39) must be constrained by boundary conditions before it is solvable. Unlike the pre-buckling problem, the inclusion of boundary conditions into the linear algebra problem is done via substitution. This is done because an eigenvalue problem must have square coefficient matrices.
Using identical methods to those used in Section 3.2 we may assemble the buckling boundary conditions as In this work, Shu and Du's (1997) method was used to introduce the membrane and bending boundary conditions into the system of equations. Rows of Eqs. (39) were replaced directly by discretised boundary conditions from Eq. (40). The constrained eigenvalue problem is therefore where the asterisk denotes matrices which have been modified by the addition of boundary condition equations. Interestingly, because we have ignored small pre-buckling rotations of the section in deriving the problem for u 1 (see A.1), the eigenvalue problem only concerns the surface normal displacement components w 1 . This approximation has reduced the total number of variables in the numerical eigenvalue problem by two-thirds. The modified and condensed eigenvalue problem is written as with where P is the compatibility matrix and K is the modified stiffness matrix. It should be noted that the removal of the membrane variables from the eigenvalue problem of shells is not a new concept. Previously, in an analysis of complete cylindrical shells in torsion, Batdorf (1947) recognised that the compatibility equation may be eliminated from the problem by using an inverse differential operator E=R, À 4 f g, the so-called Batdorf operator. The results were a modified equilibrium equation of the same order as the original Donnell equation. The matrix P performs a similar task to the Batdorf operator, albeit numerically. When multiplied by the normal displacements, or when operating on the normal displacements, it returns a compatible set of membrane displacements. The Batdorf operator operates on the normal displacement field to return a compatible set of membrane stresses.

Linear algebra problem for U 2
Converting both sides Eqs. (7) into DQ form, we find that the equilibrium equation for U 2 is given by where Q is a ð3NM Â 1Þ vector of loads. One will notice that the stiffness matrix in Eq. (43) has already been determined to be singular in Section 3.3. We therefore require Eq. (8) to render U 2 unique. Converting into DQ form, using the integration vectors and matrices, we may write it as which is a single row equation. Collecting coefficients of U 2 , Eq. (44) simplifies to where S is the orthogonality matrix with dimensions ð1 Â NMÞ. In the same manner as in the pre-buckling and buckling problems, we convert the boundary conditions into DQ form and express them as The problem for the post-buckling displacements U 2 is, therefore, given by the overdetermined system of Eqs. (43), (44) and (46). We solve this system of equations using the pseudoinverse operator:

Path coefficients
Assuming that the displacement vectors have been determined we compute the path coefficients directly using the GDQ and GIQ weighting matrices. First, we define the non-linear strains E½u i and E 11 ½u i ; u j at the grid-points as and respectively. Replacing the terms in Eqs. (17) with the GDQ and GIQ approximations, we express the path coefficients in terms of the grid-point variables as which are all scalars.
In addition to λ we know from Eq. (4) that the control displacements have a series solution. Making the appropriate substitutions we have the shortening vector as where ðd c ; d 1 ; d 2 Þ are the displacement path coefficients. In Eqs. (17) and (51) we have a pair of parametric equilibrium equations in terms of the amplitude a. Approximate equilibrium paths for each control point are computed by evaluating them for a range of values of a.

Results
The theory outlined in Section 3 was implemented in the MATLAB scripting language. In this section we compare results from the current numerical implementation to existing analytical solutions for constant stiffness shells and to finite-element solutions for variable-stiffness shells. In the latter case we make use of the modified Riks algorithm (Crisfield, 1981;Riks, 1979) implemented in ABAQUS 6.12.
In all DQ calculations the weighting matrices have been calculated using Lagrange polynomials as basis functions. The shells surface co-ordinates (x,y) have been discretised by the Chebyshev spacing rule x i ¼ l x ð1 À cos ðϕ i ÞÞ=2 and y j ¼ l y ð1 À cos ðψ j ÞÞ=2 where ϕ i ¼ πðiÀ1Þ=ðN À 1Þ and ψ i ¼ πðjÀ1Þ=ðM À 1Þ.
Laminate stiffness properties have been computed assuming that the plies are two-dimensional orthotropic laminae in plane stress, having material properties given in Table 1. Stein (1959) used the Linstedt-Poincaré method to solve the post-buckling problem asymptotically for an isotropic plate under axial compression with straight edges and simple support. This method was then extended to orthotropic plates by Chandra and Raju (1973). Conveniently, the asymptotic solutions found in these works are solutions to Eqs. (5)- (7). Only a small modification to the existing results is required to derive exact solutions for the path constants; ðλ c ; λ 1 ; λ 2 Þ and ðd c ; d 1 ; d 2 Þ.

Comparison with results for orthotropic plates
Chandra and Raju give a loading parameters β, such that λ ¼ λ c ð1 þ β 2 Þ and λu 0 0 þ v ¼ λu 0 0 þðA 1 u 1 Þβ þ u 2 β 2 ; in which the post-buckling curvature λ 2 has been set equal to λ c a priori. This restriction on the asymptotic solution is made possible by the introduction of an auxiliary amplitude of the linear mode-shape A 1 . The displacement fields u 0 , u 1 and u 2 were found by the method of undetermined coefficients. In the process, values of A 1 ensuring asymptotic convergence of the displacement series were computed. In order to make Chandra and Raju's solutions compatible with the current analysis we remove the restrictions on λ, and set A 1 ¼ 1. Substituting the displacement fields u 0 , u 1 and u 2 into Eqs. (17), we find that the path coefficients are and λ 1 ¼ 0 (i.e. the bifurcation is symmetric), where E x , ν yx and G xy are the effective laminate constants; α is the ratio E y =E x ; h is the laminate thickness; m and n are the number of half sine-waves in the x and y directions. Using Eq. (51) we also have the end-shortening relations as and d 1 ¼ 0. We now use Eqs. (52) and (53) to test the numerical implementation. The boundary conditions of the flat isotropic test plate were Straight edges: U i Sliding edges: N xy ðΦÞ ¼ 0; Simple support: 5a shows the general normalised post-buckling results for a square quasi-isotropic (QI) plate computed using the new formulation. In the pre-buckling regime the plate displays a shortening and a lateral expansion due to Poisson's ratio effects. At the critical point, the plate buckles out-of-plane into a double-sine pattern u 1 . In the post-buckling regime the linear buckling mode is accompanied by a membrane contraction u 2 , resulting in a reduced stiffness of the plate.
Numerical solutions for a angle-ply laminate ½ 7θ 3s were computed for a range of ply angles θ. Fig. 5b shows the buckling loads λ c and the reduction in stiffness k p ¼ λ 2 d c =λ c d 2 . We see that the buckling load increases steadily toward its maximum at a ply angle of 451. This is followed by a reduction towards its minimum at 901. In both the buckling load and the stiffness curves we see a sharp change in behaviour at E561. This is the result of a change in the buckling mode-shape u 1 , in which the axial half wave-number m changes from 1 to 2.
Numerical convergence of the DQ formulation was tested by repeating analyses with different grid densities. It was observed that the buckling load converged relatively 'quickly', with a maximum error of o0:72% with a 9 Â 9 grid. The postbuckling stiffness, however, was much slower to converge, requiring a 19 Â 19 grid to reduce the error to o 1% of the exact value. This is attributed to two effects. Firstly, there is a propagation of error through the analysis. The accuracy of the buckling solution u 1 is dependent on the accuracy of the pre-buckling solution, and the post-buckling solution, u 2 , is dependent on the accuracy of the buckling solution. Secondly, the accuracy of the numerical integration improves as the number of grid-points increases.

Comparison with curved isotropic panels
Next, we compare the numerical model's solutions to analytical results for a long and narrow curved panel with isotropic material properties. Koiter (1956) solved this problem, and derived closed-form solutions for the post-buckling mode shapes and the path coefficients. In the present notation these coefficients are expressed as where S is a non-linear curvature parameter defined as As expected, the path constants for the isotropic shell reduce to those in Eq. (52) with α¼1, as R-1. In addition to the critical load and the curvature, we use the buckling mode-shape u 1 from Koiter (1956) and Eq. (17b) to compute the slope of the path where m is the axial half wave-number. Interestingly, this slope has a dependence on the parity of m. The effect is due to the odd radial displacement term, w=R, which occurs in the circumferential strain component. The occurrence of complete sinusoids has the effect of cancelling out the stretching effects due to w. Odd values of m result in an overall stretching of the surface and an asymmetry of the load-path.
In the numerical analyses the following boundary conditions were used: Periodicity: The panel is assumed to be one section of a complete cylinder separated by stringers; Axial compression: N x ðΦ 1 Sliding edges: N xy ðΦÞ ¼ 0; Free periodic radial support: Note, the shell is free to expand radially at all four boundaries as long as periodicity is enforced. These boundary conditions are valid for finite radii of curvature R. For numerical calculations the magnitude of the shells curvature is limited by illconditioning of the coefficient matrices (R ¼ 1 corresponds to a rigid-body mode). Fig. 6 shows the three modes computed by the present formulation for this structure. In the pre-buckling state the panel shortens, and displays a uniform expansion in the radial direction. It then buckles into a sinusoidal mode-shape w 1 , which is accompanied by compatible membrane stretching fields u 1 and v 1 (given by Eq. (42b)). The post-buckling mode-shape u 2 is comprised of a uniform radial contraction, a uniform shortening, and a smooth sinusoidal component with double the wave-number of u 1 . Fig. 7 shows numerical solutions of the path coefficients for a 300 Â 100 mm panel (ν ¼0.3, h¼1 mm), computed using Eqs. (17) and (50). These results are compared with the exact solutions given by Eq. (54) for a range of curvatures.

Un-deformed panel geometry
Pre-buckling axial displacements, Linear buckling mode-shape, Quadratic buckling mode-shape, We observe that the buckling loads λ c and slopes λ 1 follow the analytical solutions closely up to a relative curvature l y =R ¼ 0:1. At this point Koiter's (1956) assumed mode-shape for u 1 is no longer the critical one, thus we find different trends in the path parameters. In general, we notice a sharp drop in λ 1 (due to the larger wave-number m of u 1 ) and a shallower increase in buckling load with respect to a change in curvature. Fig. 7b shows results for λ 2 and k p for a range of curvatures up to l y =R ¼ 0:1. These curves show a steady decrease in the path curvature up to l y =R ¼ 0:05 (ζ % 0:65) at which point it becomes negative and the shell is unstable. This is reflected by the stiffness k p which shows a steep transition from positive to negative at the same value.
The numerical solutions for both λ c and λ 1 computed using the DQ formulation converged to within an error of o 1% of the exact value with a (9,9) grid (where the mode-shapes agreed with Koiter (1956)). The same accuracy was achieved for the curvature λ 2 with an denser grid of (11,22). It was found that an increased number of grid-points in the y direction was required for the convergence of second-order solutions. This was due to the larger wave-number of the field u 2 and the propagation of the error from u 0 and u 1 .

Comparison with finite-element analysis
Finally, we test the new implementation's ability to predict the post-bifurcation behaviour of variable-stiffness shells. Due to a lack of exact formulae for the path coefficients of these structures, we have relied on the numerical path following methods for validation.
In order to encourage these 'full models' to follow the same equilibrium path as those given by Eqs. (17) and (51) a small proportion of the linear mode-shape u 1 (computed by ABAQUS) was superposed onto the shells nodal co-ordinates prior to the non-linear analysis. The amplitude of this imperfection was set equal to 7 h=100; 000. The sign was chosen to produce the least stable path.
In the finite-element models the S4R shell element was used. A MATLAB script was written to define the mesh, loads, boundary conditions, and to compute the element stiffness properties (on the basis of its centroid location). In all cases the panels under consideration had dimensions 100 Â 300 mm. Each finite-element mesh was chosen such that the linear buckling load λ c converged to four significant figures.
The mesh density was, for the most part, dictated by the degree of variation of the shells stiffness coefficients. The S4R element is formulated with constant material properties. Consequently, the continuous stiffness variation was represented in the finite-element models approximately by a piecewise continuous function. It was found that the densest mesh required, for the given convergence, consisted of 70 elements in the circumferential direction and 210 in the axial direction (14,700 elements and 14,981 nodes).
The boundary conditions of the problem were Straight ends (MPC): U i Axial compression: u 1 Free sides: N y ðΦ j y Þ ¼ 0; Sliding edges: N xy ðΦÞ ¼ 0; Clamped ends and simply supported sides: ∂w=∂xðΦ i The VAT plies are defined by a two-parameter univariate function: θ k ðyÞ ¼ ϕþðφ 1 Àφ 0 Þj2ξ À1jþφ 1 where ξ ¼ y=l y is a normalised circumferential co-ordinate, φ 0 and φ 1 are the tow-angles at the centre and the edges of the plate respectively and ϕ is a constant rotation. These functions are combined with standard laminate notation by including the ply angle as ϕ〈φ 0 jφ 1 〉 (Gürdal and Olmedo, 1993).
In order to test the implementation's ability to compute the pre-buckling and buckling solutions a comparison of numerical results for λ c was made. We also see a strong increase in buckling load as the central fibre angle increases from zero. Graphs have been plotted for a small range of DQ grid densities. In all cases, we see convergence with the FE solution, with a o 1% error given by a (12,24) grid. On the curve for the panel with l y =R ¼ 0:2 a relatively large error occurs when φ ¼ 301 and the grid is (8,16). This is due to a miscalculation of u 1 .
In order to test the post-buckling solutions for variable stiffness shells nine configurations of a panel with a linear ½ 7 0〈φj90〉 4s variable lay-up were chosen (Table 2). Fig. 9 shows converged path coefficients (normalised with respect to λ c ) for panels with different values of l y =R and φ. Table 2 shows particular numerical results for the path coefficients (λ 1 ¼ λ 1 =λ c andλ 2 ¼ λ 2 =λ c ) for the nine test shells. In all cases a (28,56) grid resulted in a convergence to four significant figures. On the basis of these solutions alone we see a strong influence of variable-stiffness on the post-buckling stability. In shells of curvature in the range l y =R ¼ 0:05-0:08 it is possible to find both positive and negative stability for different fibre-paths.  (1,4,7) we see a close agreement even at relatively large values of a. At these ranges of curvature the term w=R Table 2 Shell details used in finite-element simulations and converged path coefficients for each case.
Shell number 7 0〈φj90〉 (degrees) ly=R λ c Â 10 3 (N)λ 2 (mm À 2 )λ 1 Â 10 À 3 (mm À 1 )  unstable shells (3,6,9) a steep drop in load occurs at the point of buckling due to the negative stiffness effects. This behaviour is captured by the asymptotic solution, however, only in the very close neighborhood of the bifurcation. In the unstable shells the Riks results show the stiffness recovering over a relatively short path, which the Koiter model fails to capture. The reason for the disagreement between the two models is that the quadratic expansion of the displacement field v, used in the Koiter model, lacks the degrees of freedom required to capture the correct post-buckled shape u remote from the bifurcation point. The Riks algorithm, on the other hand, enforces equilibrium at every increment and so does not have this problem. In order to compute the correct behaviour using Koiter's method, a higher-order expansion of the displacement field v is required.

Conclusions
An implementation of Koiter's (1945) method for variable-stiffness shell structures was presented. The new formulation was used to model orthotropic plates, isotropic cylindrical panels and variable-stiffness cylindrical panels under axial compression. Numerical results for the displacement fields and path coefficients were in good agreement with the exact formula from Stein (1959), Chandra and Raju (1973), and Koiter (1956), thus validating the use of the pseudo-inverse for implementing the boundary conditions and orthogonality constraints.
Numerical results for the buckling loads, dependent on both the pre-buckling and buckling mode-shapes, were in good agreement with finite-element analysis, thus validating the first two-steps of the implementation. The second-order term could not be compared directly to any FE solutions, they did, however, converge to three decimal places for the test structures with a relatively small number of grid points. The second-order paths computed by the current implementation were compared to arc-length based non-linear solutions. It was found that in all cases the theoretical slope of the path at the bifurcation points was in good agreement with the Riks results. The validity of this solution, however, in the deeper postbuckling range was dependent on the stability of the bifurcation (for the shells tested). Those structures exhibiting stable or neutrally stable bifurcations had better predictions of their post-buckling equilibrium paths in the more distant neighborhood of the bifurcation, while those shells having a negative value of λ 2 showed good agreement at the buckling point but with only a small range of validity. In the full-models, a distinct stiffening of the equilibrium paths is exhibited which cannot be calculated with a second-order model. In panels with supported sides it is known that this stiffening effect always occurs (Gerard and Becker, 1957). On this basis it may be said that in the present analysis the point of minimum stability over the course of the loading is attained.
The strength of the new model is its ability to predict the general stability of the shell at the buckling load with a very small number of linear operations. It is proposed as a potential tool for optimisation studies of variable-stiffness panels where stability is a constraint.
where p and Δ are vectors of discrete forces and displacements at the boundary control points, and A is the area of the shell.
(1) into Eq. (A.1) and collecting terms by powers of u, we write where the functionals P 1 …P 4 are the linear, quadratic, cubic and quartic changes in elastic strain energy for a given change in the displacement fields u. According to Hamilton's principle, the shell is in a state of static equilibrium if the displacement field u satisfies the following variational equation: where W is the space of kinematically admissible displacements, δP is the first variation of the potential P with respect to the arbitrary field η and the notation ðÞ 0 denotes the variational derivative. Evaluating Eq. (A.3) using the Gateaux derivative, we write the equilibrium equation: where τ is an arbitrary scalar, and the subscripts of the functionals refer to the powers of each their arguments in a binomial expansion. For example, P 2 ½u 1 þ u 2 ¼ P 2 ½u 1 þP 11 ½u 1 ; u 2 þP 2 ½u 2 . In order to simplify the notation in the proceeding analysis we introduce the following auxiliary variables: S 1 ½u ¼ P 11 ½u; η; S 2 ½u ¼ P 21 ½u; η and S 3 ½u ¼ P 31 ½u; η which are implicit functions of η. They are given by the following expressions: In the pre-buckling state it is assumed that rotations ∂w 0 =∂x and ∂w 0 =∂y are sufficiently small that the stresses N 0 are approximately equal to A Á e 0 . That is, the higher-order strains E 0 are assumed to be negligible. Substituting the displacements (where v is an arbitrary displacement field,û 0 is the pre-buckling mode and λ is the loading parameter) into Eq. (A.4) yields the following: P 1 ½ηþS 1 ½u 0 þS 1 ½vþS 2 ½u 0 þS 11 ½u 0 ; vþS 2 ½vþS 3 ½u 0 þS 21 ½u 0 ; vþS 12 ½u 0 ; vþS 3 ½v ¼ 0 ðA:6Þ The small-rotations assumption facilitates a considerable simplification of the analysis. At this stage we may assume that S 2 ½u 0 % 0; S 3 ½u 0 % 0 and S 21 ½u 0 ; v % 0 thus reducing the number of terms in (A.6). It should be noted that these approximations limit the scope of the valid initial loading states. Panels having combined loads , membrane constraints at the boundary (Hilburger et al., 2004) or circular perforations  all contradict these assumptions and are, therefore, outside the scope of the current method. Nevertheless, for simple in-plane loads we may find a meaningful variational problem for u 0 by setting v to zero. This gives the following linear variational problem: P 1 ½ηþS 1 ½u 0 ¼ 0; ðA:7Þ Integrating Eq. (A.7) by parts in order to isolate the zeroth-order derivatives of η and noting that η is arbitrary, we find that Eqs. (5) and (9), and either (11) (or (10) and (12)) must hold for equilibrium. Examples of this procedure (which is essentially the applications of Gauss's divergence theorem) may be found in many texts, including Yamaki (1984), van der Heijden (2009) and Reddy (1984).

ðA:9Þ
which implies that S 1 ½u 1 ¼ Àλ c S 11 ½û 0 ; u 1 ðA:10Þ This is a linear eigenvalue problem for the field u 1 and the load λ c . It may be simplified by noting that in S 11 ½û 0 ; u 1 the biggest term by far is Z A ðN T 0 Á E 11 ½u 1 ; ηÞ dA and hence is the only one considered. Applying the divergence theorem, in the same way that it has been to Eq. (A.7), we find that Eqs. (5) and (9), and either (13) (or (10) and (15)) must hold. In general, Eq. (A.10) has no unique solution. In this work we choose the solution which corresponds to the lowest critical load λ c .
Next, we focus on determining the stability of the adjacent state with a finite but small value of the perturbation parameter a. We define three further functionals: R 1 ½v ¼ S 1 ½vþS 11 ½u 0 ; v; R 2 ½v ¼ S 2 ½vþS 12 ½u 0 ; v; and R 3 ½v ¼ S 3 ½v which are implicit functions of η and u 0 . Now, we add a quadratic correction field u 2 to v to give v ¼ au 1 þ a 2 u 2 Substituting this into Eq. (A.8) gives the following expression: R 1 ½u 1 þaðR 1 ½u 2 þR 2 ½u 1 Þ þ Oða 2 Þ ¼ 0 ðA:11Þ The first term is equal to zero at the critical load (see Eq. (A.10)). Ignoring higher-order terms in a we have the following linear problem for u 2 : R 1 ½u 2 ¼ ÀR 2 ½u 1 ðA:12Þ In addition to Eq. (A.12) further constraints on u 2 are required in order to render them unique. This is because at λ ¼ λ c its left-hand side is singular. An efficient choice for this constraint is one which makes the components of the series for v orthogonal in energy-space (Koiter, 1945). The following energy functional F 2 ½u ¼ ÀP 12 ½û 0 ; u ¼ ÀS 2 ½u;û 0 is positive-definite. It, therefore, holds that F 2 ½u þτw Z F 2 ½u with w A W and F 11 ½u; wþτF 2 ½w Z 0 ðA:13Þ We define fields u and w as orthogonal if they satisfy the equation: F 11 ½u; w ¼ 0; ðA:14Þ meaning that the components of energy derived from the two fields may be superposed (i.e. F 2 ½u þ w ¼ F 2 ½uþF 2 ½w).
The proof that Eq. (A.12) has a unique solution as long as u 1 and u 2 are orthogonal is given by van der Heijden (2009). It is shown that the solution for u 2 is independent of the actual choice for the functional used to enforce orthogonality. The only requirements being that it is positive-definite.

A.2. Path parameters
The constants in the series (3) have been determined in the following way. Consider difference of the potential energy of the structure of the in the adjacent and fundamental states. This is P λ ¼ P½λû 0 þvÀP½λû 0 ¼ P 2 ½vþλS 2 ½v;û 0 þP 3 ½vþP 4 ½vþOðλ 2 Þ ð A:15Þ The same expression may be found by calculating a Taylor series for P in terms of λ, about the point λ ¼ λ c (Budiansky, 1974).
According to Koiter (1956), P λ should be minimised with respect to the perturbation parameter a for equilibrium.

Appendix B. Matrix notation
The ○ product has been defined to permit the use of matrix multiplication of sets of DQ variables. The operator has the following simple meaning. Arrays of matrices are multiplied according to the standard matrix dot product. The matrices with the arrays are multiplied in a row-wise manner. For example,