CC BY
0Contents lists available at ScienceDirect
Journal of the Mechanics and Physics of Solids
journal homepage: www.elsevier.com/locate/jmps
Initial post-buckling of variable-stiffness curved panels
S.C. White n, G. Raju, P.M. Weaver
Advanced Composites Centre for Innovation and Science, Department of Aerospace Engineering, Queen's Building, University Walk, United Kingdom
CrossMark
ARTICLE INFO
Article history:
Received 18 December 2013
Received in revised form
1 July 2014
Accepted 5 July 2014
Available online 14 July 2014
Keywords: Shell buckling Post-buckling Koiter's method Variable angle tow Differential quadrature
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. © 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC
BY license (http://creativecommons.org/licenses/by/3.0/).
1. 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, 2005), causing linear design calculations to be non-conservative. There have been many theoretical investigations into these effects. Cox and Pribram (1948) characterised the non-linear load path of isotropic panels qualitatively and found that, in the same manner as complete cylinders (von Kármán and Tsien, 1941), a degree of unloading and shortening occurs immediately after buckling (sometimes termed "double-backing"). This behaviour is in contrast to flat plates, whose load paths are monotonic and roughly piecewise-linear. Koiter (1956) quantified the effect of curvature on a panel's stability by computing the tangent stiffness of the load path in the near vicinity of the buckling point. It was found that a small
* Corresponding author. E-mail address: simon.white@bristol.ac.uk (S.C. White).
http://dx.doi.org/10.1016/jjmps.2014.07.003
0022-5096/© 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/).
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 (1965, 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. (1999, 2001, 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., 2013a, 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).
Fig. 1. Load vs. end-shortening curves for isotropic (v=0.3) curved panels, under compression in the axial (x) direction (with a range of initial curvatures k). Radial displacements, tangential rotations and axial rotations are zero on all edges. Tangential and axial displacements are free on all edges. [Reproduced from Pope (1965)].
End shortening
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.
2. 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 9k(x,y) and that the laminate has no significant coupling between membrane and bending deformation (i.e. B = 0). The panel has length lx, width ly, thickness h and a radius of curvature R. The variation of the reinforcement direction results in a laminate with the following constitutive properties:
N = A(x, y)-e, and M = D(x, y)-K, (1)
where A and D are the membrane and bending stiffness matrices, and
e = [Cx ey Cxy]T, K = [Kx Ky Kxy ]T, N = [Nx Ny Nxy]T,
and M = [Mx My Mxy]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
e = e+E and k = k
e(u) =
du Idv W dx Vdy +1
dy dx)
E(U) = 2
dW2 dW'
2^w dW dx dy
k(u) = -
rd2W d2W 2 d2 W
dx2 dy2 dxdy
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 =R2 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 z-direction are larger than those in the x and y directions.
Fig. 2. The panel geometry with the current co-ordinate system.
2.1. 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:
X = Xc + aXi +a2 Â2 + ■■■ (3)
u(M, a) = u0+au1 +a2U2, (4)
where a is a perturbation parameter. The fields u0, u1 and u2 are the pre-buckling and linear and quadratic buckling modes. The coefficients Xa M and M2 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 M; thus u = u0 = Mu0. At some critical load value, Mc, the shell bifurcates with a corresponding linear mode-shape u1. In the analysis, the stability of the linear mode-shape is determined by computing the second-order field u2 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 u1 and u2), 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 u0, u1 and u2 have been derived in A.1 using the principle of stationary potential energy. Introducing the following notation for the stress-resultants:
Nm = A • e(um), Mm = D • k(um), and N mm = A • E(um),
with m a {0,1,2}, we write
Pre-buckling:
dNxo àNxyo 0 dx dy '
dNyo dNxyo 0 dy dx = ;
d2Mxa ! 2d2MXyo dx2 " dxdy d2Myo + dy2 +"?=o-
Buckling:
dNxi dNXyi_0 dx dy '
dNyi dNxyi _0 dy dx = ;
d2Mxi 2d2MXyi dx2 ' " dxdy d2Myi + dy2 Nyi + R
Post-buckling:
dNx2+ dNxy2 _ dx dy ' dNxii dNxyii
dx dy '
dNy2 + dNxy2 dy dx ' dNyii dNxyii
dy dx '
d2Mx2 ! 2d2Mxy2 dx2 " dxdy d2My2 + dy2 +¥+ Xc!
(5a) (5b) (5c)
(6a) (6b)
d Ídw1 dwA , d ( dw1 dwA
Tx{N*0*r+»^J - XTy[Nx°-3"+»«i») (6c)
{*№£+ »-w) +1^+n-t?))
(7a) (7b)
d / dwi „ dwA d / dwi „ dwA Nyii
=--Nx1—"+Nvv"—"--Nx1—"+Nx1—1 H—y"", (7c)
x1 dx xy1 dy ) dy\ x1 dx x1 dx ) R ' y 7
L(N^N^) ^ {»^»■»WwdA = (8,
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 u1 and u2. It is required to ensure uniqueness of the solution to u2 (Koiter, 1945).
2.2. Boundary conditions
All external forces in the model are applied at the boundary (see Fig. 3). The edges corresponding to x = {0, lx} are termed
ends and given the symbols 01 and 02. Similarly, the edges corresponding to y = {0, ly} are termed sides and have the
symbols 01 and 02.
Bending boundary conditions are identical for each of the displacement fields u0, u1 and u2 (due to the linearity of the changes in curvature k). They are written as
Ends 0kx: Mxm = 0, or ^ = 0, (9a)
and Qxm = 0 or wm = 0, (9b)
Sides 0ly: Mym = 0, or — = 0, (9c)
and Qym = 0 or Wm = 0, (9d)
where m is the state (i.e. 0 pre-buckling, 1 buckling and 2 post-buckling). The fields Q and Qy are the transverse shear stress resultants (for a shallow-shell) and are given by the expressions:
Qx = mx + d_Mxy and Qy = dMy+dMxxi
dx dy dy dx
In our model the following (optional) membrane constraints are also common to u0, u1 and u2, they are
Ends 0: Um =[U'k(y) + u'km V'x(y) + v'xm Wm]T (10a)
Sides 0ly: Um = [Uy(x) + uym vy(x) + vym Wmf (10b)
The displacement functions Ukx and Vxk, and Uy and Vy, are the fixed axial and circumferential displacements at the ends and sides, respectively. The scalars uxm, vlkm, ulym, and vlym 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.
2.2.1. 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 A. Fixed stress boundary conditions are
Ends 0k: Nxc = F'k(y), = S'k(y) (11a)
Sides 0ly: Nye = Fy(x), = Sy(x) (11b)
where Fk and Sy, and Slx and F'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 finite-element formulations. Forces are applied to the edge control points, each of which has two degrees of freedom. We write the
equilibrium equations of these points
Ends 0
Nxody = f'X
dy = p'k
Sides 0
fix rix
'y'- J Nyo dx = ply, J^ Nxyo dx = fy
where the variables fk, pk, fy and p'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.
2.2.2. Buckling and post-buckling membrane boundary conditions
Due to the non-linearity in the strains e, the membrane force boundary conditions are slightly different for each of the displacement fields u0, u1 and u2. At edges where standard stress boundary conditions (11) have been applied, we have the following
Buckling:
Ends 0k- Nxi = 0, Nxyi = 0
Sides 0
Ny1 = 0 Nxy1 = 0
Post-buckling:
Ends 0k: Nx2 = —Nx11, Nxy2 = — Nxy11
Sides 0
Ny2 = — Nyll,
Nxy2 = —Nxy11
(13a) (13b)
(14a) (14b)
Alternatively, at edges where force boundary conditions (12) have been applied, we have Buckling:
Ends 0
Nx1 dy = 0,
Nxy1 dy = 0
Sides 0y :
Ny1 dx = 0,
Nxy1 dx = 0
Post-buckling:
Ends 0
Sides 0l
Nx2 dy =— / Nx11 dy,
Ny2 dx = —
/ Nxy2 dy = — / 00
Nyn dx,
Nxy2 dx = —
Nxy11 dy
Nxy11 dx
(16a) (16b)
Eqs. (15) and (16) must be used in conjunction with Eq. (10). 2.3. Path coefficients
Assuming that solutions for u0, u1 and u2 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 A for a given change in the perturbation parameter a. Expanding A 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:
Ac = —1jA (NÏ • e1 +MÎ ■ k^dA
-2d/a( N1 • £0dA
• A • E1 +e2 • A • e2 + k2 • D • k2+2eï • A • £n+2e2 • A • E1+2AcN0 • dA
(17b) (17c)
D = (N0 ■ Ei) dA
dw1 dw2 dw1 dw2 / dw1 dw2 dw1 dw2\ dx dy dy dx J
dx dx dy dy
The stability of a shell in its initial post-buckling range may be determined from the path parameters. The slope indicates asymmetry of A 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 = A2/Ac), 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.
3. Differential quadrature implementation
In this paper a mixed differential-integral quadrature (GDQ-GIQ) technique has been used to solve the integro-differential 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
df=tx • f
f = [f 1-../N]T and fi = f(xt)
where the matrix Tx 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| 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 x M grid of points (xi}yj) and a two-dimensional function-array f (xi}yj) =fj. As a matter of computational convenience it is useful to hold the discretised function fij 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
f =[f11---f1M f21 ■ ■ -f2M ■■■ fN1 ■■fNM]T
where it now has dimensions (NM x 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 fj and f it is possible to write the partial derivatives of the two-dimensional field f(x,y) in DQform as
—f = Fq • f,
(N - 1)M+1
Fig. 4. A two-dimensional Chebyschev grid plotted on the surface of a curved panel with the grid-point numbering system.
-f = Fqr. f
where the matrices Fq and Fqr both have dimensions (NM x 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, GIQuses 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
g =y f (x) dx,
we write the discrete indefinite-integral of f in vector form as
g =[g1--gN] where gi = g(x), This is computed approximately by the expression:
g ^(Tx)+ • f+1
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
iq = (Ty
where Iq has the same dimensions as Tq. The use of the pseudoinverse operator is required because Tq 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
jx = (IXN}- I} and jy = (4 - ¡1 ) (24)
which have dimensions (1 x N) and (1 x M), respectively. The vectors jx and jy 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 jx and jy are cast into the matrices Jx and Jy, such that
J f (x, y) dx-Jx ■ f with } = 1...M
f (xi, y) dy-Jy ■ f with i = 1...N
where Jx and Jy have dimensions (M x NM) and (N x NM), respectively. Multiplying the vector of function values f by the matrix Jx results in an (M x 1) vector of definite x-direction integrals of f at different co-ordinates yj. Likewise, multiplication of f by Jy yields a vector of y-direction integrals at different co-ordinates xi.
3.1. Variable discretisation
Making use of the field discretisation in Eq. (19), we represent the displacement fields u, v and w at the DQgrid points by the (NM x 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) we obtain the strains and the curvatures at the grid points as
" ex " Fx 0 0 ' u "
ey = 0 Fy I/R v = E ■ U (26a)
exy Fy Fx 0 w
kx 0 0 Fxx
ky = - 0 0 Fyy
kxy 0 0 2Fxy
where I is an (NM x NM) identity matrix, E and K are the strain and curvature co-efficient matrices, and
U = [u v w]T.
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
sym. Dii sym.
A16 " ex "
A26 o ey
A66 exy
D16 kx
D26 o ky
D66 kxy
= A°(£ • U)
= Do(K • U)
where the vectors Aj and Dj contain the stiffness coefficient values at all the grid points and have dimensions (NM x 1). The o matrix product is defined in Appendix B. On the basis of Eq. (27), it is possible to write the stresses and bending moments in terms of the displacement vector directly. Collecting coefficients of U we get
nx "Nx"
ny = Ny • U = N • U
n xy Nxy
mx Mx
my = My • U = M • U
mxy Mxy
where (Nx, Ny, Nxy) and (Mx, My, Mxy) are the membrane and bending stress-resultant coefficient matrices. They each have dimensions (NM x 3NM).
3.1.1. Common boundary conditions
Bending boundary conditions (9) and kinematic membrane boundary conditions (10) are common to each of the vector fields u0, u1 and u2. The bending boundary conditions are expressed as
Ends 0ix:
Sides 0>y :
(Mx)tf • U = 0, or (F% • w = 0
(Fx • Mx + Fy • Mxy)<p • U = 0 or
• w = 0
(My)^ • U = 0, or (Fy )<p • w = 0
(Fy • My + Fx • Mxy)0 • U = 0 or
• w = 0
(29a) (29b) (29c) (29d) (29e) (29f) (29g) (29h)
where the 0 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:
BB -u = 0
where BB is the bending boundary coefficient matrix which has dimensions ((2N+2M) x NM).
Next, by sampling the prescribed membrane boundary displacements Uj, Vk, 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:
C • U = V (31)
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 x NM) and ((2N+2M) x NM).
3.2. Linear algebra problem for îi0
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
Fx o F1" o -
F • N • U o = o
{[Fxx Fyy 2Fxy] • M + R -1 Ny} • Uo = 0
Lxu Lxv Lxw "Uo "
Lyu Lyv Lyw Vo = L • U o = o
Lzu Lzv Lzw Wo
where Eq. (32) is the membrane equilibrium equation, Eq. (33) is the bending equilibrium in the equation, and îi0 is the pre-buckling displacement vector. Combining Eqs. (32) and (33) and collecting the various coefficients of îi0 we may write the following system of linear equations:
where the matrices Lj 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%, Sk, Fy and Sy 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:
BM • îio = Fo (35)
where BM 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
Ends 0
jyT °(NX )<p
U o = fkx and
Sides 0
°(Nxy)<p • Uo = p'X
jxT°(Ny)0 • Uo = ft and
(36b) (36c)
jxT °(Nxy)<p • U o = ply
The integration vectors jx and jy have been used to integrate the stresses along each edge based on the global displacement vectors. Coefficients of U'0 in Eq. (36) are collected and assembled into BM. The forces fX, pX, fy and p'y are contained in the vector Fo.
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 U0, hence
■ L - + - o -
BM F o
_ C _ . V _
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^ = p^ = 0; u1 = 0;f^ = X.
3.3. Linear algebra problem for U1
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 w1. 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
AC(FX .{nxo°FX + nxyo°Fy}+-Fy • {nyooFy+nxyo°FX})-wi = XcGZw • wi, (38)
where G,™ is the so-called geometric stiffness matrix. The unconstrained eigenvalue problem for U1 is, therefore, given by
L • U1 = XcG • U1 (39)
where Ac is the critical buckling load. The array G has only one non-zero sub-matrix G^. 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
C • U1 = 0; BB • U1 = 0; and BM • U1 = 0 (40)
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
L* • U1 = XcG* • U1 (41)
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 u1 (see A.1), the eigenvalue problem only concerns the surface normal displacement components w1. 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
K • w1 = XcG*, • w1 (42a)
= P • w1
L L* -1
L L L
and K = L*„ + [L*u L*
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/RV ~4{ }, 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.
3.4. Linear algebra problem for U2
Converting both sides Eqs. (7) into DQ form, we find that the equilibrium equation for U2 is given by
{£+XcG} • U2 = Q (43)
where Q is a (3NM x 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 U2 unique. Converting into DQform, using the integration vectors and matrices, we may write it as
jy • Jx • {(nxo°Fx • wi)°Fx + (nxyo°Fx • wi)°Fy + (ny0°Fy • wi)°F* + (n^F* • wi)°Fx} • w2 = 0 (44)
which is a single row equation. Collecting coefficients of U2, Eq. (44) simplifies to
S • w2 = 0 and [0 0 S] = S • U2 = 0 (45)
where S is the orthogonality matrix with dimensions (1 x 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
C • U2 = 0, BB • U2 = 0, and BM • U2 = F2 (46)
The problem for the post-buckling displacements U2 is, therefore, given by the overdetermined system of Eqs. (43), (44) and (46). We solve this system of equations using the pseudoinverse operator:
"{£+AcGY + ' Q '
BM F 2
3.5. 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 £[ui] and £11[ui, Uj] at the grid-points as
2!(Fx • wi)o(Fx • wi;
1( Fy • wO o( F • wO (Fx • wi)o(Fy • wi)
(Fx • wi)o(Fx • wj) (Fy • wi)o(Fy • wj) (Fx • wi)o(Fy • wj) + (Fx • wj)o(Fy • wi)
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
D = jy • Jx •[(N • Uo)ToRi] 1.
j • Jx • |_(N • U1)1 o(E • Ui) + (M • U1)1 o(K • U1 )j
Ai = "2Djy • Jx • [(N • U1)Tor,
A = -Dj • Jx • [r! o(AoR1 ) + (N • U2)ToE • U2 + (M • U1)To(K • U1)
+ 2(N • Ui)ToRi2 + 2(N • U2)ToRi +2Ac(N • Uo)ToRz]
(50a) (50b)
(50c) (50d)
which are all scalars.
In addition to X we know from Eq. (4) that the control displacements have a series solution. Making the appropriate substitutions we have the shortening vector as
A(a) = A(a)A 0+Aia+A2a2 = AcA 0 + (Ai 0^i)a+(A2 + 0^2)a2 = dc+di a+d2a2
where (dc; d1; d2) 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.
4. 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
Xi = lx(1 - cos (0*))/2 and yj = ly(1 - cos (^,))/2
where = n(i - 1)/(N-1) and ^ = n(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.
Table 1
Uni-directional carbon-fibre ply properties.
Material E11 (N/mm2) E22 (N/mm2) G12 (N/mm2) v12 t (mm)
IM7 8552 163 x 103 12 x 103 5 x 103 0.3 0.131
4.1. Comparison with results for orthotropic plates
Stein (1959) used the Linstedt-Poincare 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; (Xc,X1; X2) and
(dc, di, ¿2).
Chandra and Raju give a loading parameters ft, such that
X = Xc(1 + ft2) and Xu0 +v = Xu0 + (A1 u1)ft+u2ft2,
in which the post-buckling curvature X2 has been set equal to Xc a priori. This restriction on the asymptotic solution is made possible by the introduction of an auxiliary amplitude of the linear mode-shape A1. The displacement fields u0, u1 and u2 were found by the method of undetermined coefficients. In the process, values of A1 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 X, and set A1 = 1. Substituting the displacement fields u0, u1 and u2 into Eqs. (17), we find that the path coefficients are
Xc =—2 3n h-(Exa2l4n4 + Exal4m4 + {2ExaVyX + 4GXy(a-v2x)}l2xllm2n2),
12l;l3m2(a - v2x)
„ Exn2h(a(4n4 +14 m4)
X2 = --Mh-y—- (52)
2 16l2l3m2
and X1 = 0 (i.e. the bifurcation is symmetric), where Ex, vyx and Gxy are the effective laminate constants; a is the ratio Ey/Ex; 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
. Xclx . (X2lx , lx(mn\2 .
dc= - m d2 = yMxxiy+A-) 1- (53)
and d1 = 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: Ux = % = U]y = Vy = 0;
• Axial compression: u1 = v1 = 0, f^ = -1, u1 = vj = 0, p^ = 0;
• Sliding edges: Nxy(0) = 0;
• Simple support: Mx(^) = My(@ly) = 0, w(@) = 0
Fig. 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 u1. In the post-buckling regime the linear buckling mode is accompanied by a membrane contraction u2, resulting in a reduced stiffness of the plate.
Numerical solutions for a angle-ply laminate [ + ff\3s were computed for a range of ply angles d. Fig. 5b shows the buckling loads Xc and the reduction in stiffness kp = X2dc/Xcd2. We see that the buckling load increases steadily toward its maximum at a ply angle of 45°. This is followed by a reduction towards its minimum at 90°. In both the buckling load and the stiffness curves we see a sharp change in behaviour at e 56°. This is the result of a change in the buckling mode-shape u1, 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 < 0.72% with a 9 x 9 grid. The post-buckling stiffness, however, was much slower to converge, requiring a 19 x 19 grid to reduce the error to < 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 u1 is dependent on the accuracy of the pre-buckling solution, and the post-buckling solution, u2, is dependent on the accuracy of the buckling solution. Secondly, the accuracy of the numerical integration improves as the number of grid-points increases.
Q (Degrees)
-------- Exact —o— DQ (N = M = 9)
—.— DQ (N = M = 19) —n— DQ (N =M = 15)
Fig. 5. (a) Normalised post-buckling equilibrium path with displacement modes superimposed for a square QI flat panel. (b) Buckling load kc and post-buckling stiffness kp — (A2dc/Acd2) vs. fibre orientation ft The plate is 100 mm square, and has the stacking sequence [ + ft]3 s.
4.2. 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
. Eh3n2 Ehly
Xc — ^-—I--% and
c 3(1 - v2)ly 4n2R
. En2 h(1 -9S) 3El3(1 -v2) r54,
A2 ——sy---(54)
where S is a non-linear curvature parameter defined as
S — 2 [(r2 + 1)-2+Z-4(r2 +1 )2+z-4-1]-1
r — 1
with Z — ^P2) -ffiffi (55)
As expected, the path constants for the isotropic shell reduce to those in Eq. (52) with a — 1, as R —i. In addition to the critical load and the curvature, we use the buckling mode-shape u1 from Koiter (1956) and Eq. (17b) to compute the slope of the path
2Elyh(1 -(-1)m) M — n2Rm '
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: Nx(01x) = -Nx(0;) = XF, Ny0) = 0;
• Sliding edges: Nxy(0) = 0;
• Free periodic radial support: Mx(0ix)=My(0iy) = 0, Qx(0x) = Qy(0y) = 0
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 ill-conditioning 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 w1, which is accompanied by compatible membrane stretching fields u1 and v1 (given by Eq. (42b)). The post-buckling mode-shape u2 is comprised of a uniform radial contraction, a uniform shortening, and a smooth sinusoidal component with double the wave-number of u1.
Fig. 7 shows numerical solutions of the path coefficients for a 300 x 100 mm panel (v=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.
Pre-buckling axial displacements, u o
Un-deformed panel geometry
Linear buckling Quadratic buckling
mode-shape, u i mode-shape, U2
Fig. 6. Mode-shapes for a shallow (i.e. Z < 1) isotropic shell under axial compression (lx=300 mm, ly = 100 mm, R=1000 mm, v=0.3). (a) Pre-buckling mode u0 with axial displacements superimposed. (b) Buckling modes u1 and u2 with the displacement magnitudes (vector norms) superimposed.
--------------- Koiter [11] -=-DQ (9,9) -«-DQ (21,21)
-o-DQ (21,41) -.-DQ (21,67)
........Koiter [11] - -DQ (15,30)
—□—DQ (9,18) - -DQ (11,22)
Fig. 7. Numerical results for the path coefficients of a curved isotropic panel (ly=100 mm, h = 1 mm, v=0.3, lx=300 mm). (a) Normalised buckling load Xc = Xc3ly(1 - v2)/Eh3n2 and slope X] = Xi(Xc/Xc) vs. relative curvature ly/R. (b) Normalised path curvature X2 = X2 /Xc and post-buckled stiffness kp vs. relative curvature ly/R.
We observe that the buckling loads Xc and slopes X1 follow the analytical solutions closely up to a relative curvature ly/R = 0.1. At this point Koiter's (1956) assumed mode-shape for u1 is no longer the critical one, thus we find different trends in the path parameters. In general, we notice a sharp drop in X1 (due to the larger wave-number m of u1) and a shallower increase in buckling load with respect to a change in curvature.
Fig. 7b shows results for X2 and kp for a range of curvatures up to ly/R = 0.1. These curves show a steady decrease in the path curvature up to ly/R = 0.05 (Z ~ 0.65) at which point it becomes negative and the shell is unstable. This is reflected by the stiffness kp which shows a steep transition from positive to negative at the same value.
The numerical solutions for both Xc and X1 computed using the DQ formulation converged to within an error of < 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 X2 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 u2 and the propagation of the error from u0 and u.
4.3. 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 (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 + 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 x 300 mm. Each finite-element mesh was chosen such that the linear buckling load Xc 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'x = 0, u1 = 0;
• Axial compression: uj = vj = 0, f2 = -1, uj = vj = 0, p2 = 0;
• Free sides: Ny(0iy) = 0;
• Sliding edges: Nxy(0) = 0;
• Clamped ends and simply supported sides: dw/dx(@'x) = 0, My(0y) = 0, w(@) = 0;
The VAT plies are defined by a two-parameter univariate function:
8k(y) = <P + (P1 - ^0)12£ -11+ P1
where £ = y/ly is a normalised circumferential co-ordinate, p0 and p1 are the tow-angles at the centre and the edges of the plate respectively and p is a constant rotation. These functions are combined with standard laminate notation by including the ply angle as p(p0\p1) (Gurdal 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 Xc was made.
Fig. 8 shows normalised results for panels with a linear [ + 0(p\0)]4 s variable lay-up (see inset Fig. 8) and with a range of different initial curvatures ly/R. The results have been normalised with respect to a flat panel, with p = 0. As expected we see an increase in the buckling load with an increase in curvature (reflecting Eq. (54a). 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 < 1% error given by a (12,24) grid. On the curve for the panel with ly/R = 0.2 a relatively large error occurs when p = 30° and the grid is (8,16). This is due to a miscalculation of Ui.
In order to test the post-buckling solutions for variable stiffness shells nine configurations of a panel with a linear [ + 0(p\90)]4s variable lay-up were chosen (Table 2).
Fig. 9 shows converged path coefficients (normalised with respect to Xc) for panels with different values of ly/R and p. Table 2 shows particular numerical results for the path coefficients (X 1 = X1/Xc and X2 = X2/Xc) for the nine test shells.
20 30 40
^ (Degrees)
- dq (8,16)
- abaqus
-dq (12,24) -(10,20)
Fig. 8. A comparison of buckling eigenvalues for a panel with a [ + 0<p\0)]4s variable lay-up. The buckling loads have been normalised with respect to the buckling load of an equivalent flat (i.e. ly/R = 0), [0]8 panel (Eq. (52).
S.C. White et al. / J. Mech. Phys. Solids 71 (2014) 132-155 149
Table 2
Shell details used in finite-element simulations and converged path coefficients for each case.
Shell number + 0(<^|90) (degrees) ly/R Ac x 103 (N) Â2 (mm-2) Â1 x 10-3 (mm-1)
1 90 0.02 1.160 0.223 1.862
2 90 0.083 1.341 0.024 6.349
3 90 0.1 1.425 - 0.052 6.857
4 75 0.02 1.278 0.199 1.774
5 75 0.067 1.415 0.028 4.937
6 75 0.1 1.62 - 0.156 4.997
7 60 0.02 1.621 0.15 0.870
8 60 0.05 1.714 0.014 0.959
9 60 0.08 1.879 - 0.192 0.016
(Degrees)
^ (Degrees)
Fig. 9. Numerical results for normalised slopes X1 = Xi/Xc and curvatures X2 = X2/Xc of curved panels with a [ + 0(^|90>]4s variable lay-up.
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 ly/R = 0.05-0.08 it is possible to find both positive and negative stability for different fibre-paths.
Fig. 10a-c shows comparison between the Riks and our quadrature solutions. The load factors and control displacements have all been normalised with respect to the numerical values of Xc computed by the FE analysis. In the cases of the most stable shells (1,4,7) we see a close agreement even at relatively large values of a. At these ranges of curvature the term w/R
A 1 = 0.02 D
/ ly/ u ^.......
VS. \\ s \ \ \ ^ X \ N ly/R = 0.083 )
ls/ R = 0.1 <p = 90°
0.99 1.00
1.01 A/dc
1.02 1.03
1.03 1.02 1.01 1
0.99 0.98 0.97
0.96 0.
y R = 0 02
R = 0 w 067
ly/R = 0.1
<p = 75°
1 1.02 1.04 1.06 1.08 A/dc
y/R = 0.02 (4)
— /R =( >.05 -
/A s____ , ©
/ ly/ R =0. 08 = 6(
/ H> )
0.96 0.
1 1.02 1.04 1.06 1.08 A/dc
Fig. 10. Riks (solid) and DQ-Koiter (dashed) equilibrium paths for variable-stiffness cylinders in the post-buckling range. (a) Shells 1-3 (p = 90°). (b) Shells 7-9 (p = 75°). (c) Shells 4-6 (p = 60°).
in the circumferential strain ey is small and the bifurcations are approximately symmetric. Next, shells (2,5,8) are almost neutrally stable and have an asymmetric bifurcation. The sign of a was chosen such that it corresponded to the least stable equilibrium branch. In these shells good agreement between Riks and DQ-Koiter is found for moderate values of a, however, a distinct stiffening in the more remote neighborhood of the bifurcation is observed in the Riks results. Finally, in the
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.
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 post-buckling 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 X2 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.
Acknowledgements
The authors would like to thank the UK Engineering and Physical Sciences Research Council (EPSRC) for funding the work carried out under the ACCIS DTC. Grant No. EP/G036772/1.
Appendix A. Derivation of the governing equations
A.1. Governing equations
Consider the structure shown in Fig. 2. Under external loading the total potential energy is equal to the total strain energy, minus the work done by loads. We write this as
where p and A are vectors of discrete forces and displacements at the boundary control points, and A is the area of the shell. Substituting Eqs. (1) into Eq. (A.1) and collecting terms by powers of u, we write
5. Conclusions
P[U]=Pi +P2 +P3 +P4
Pi[u]=-pT ■ A(u),
?4[U]= iy (ET • A • e) dA.
where the functionals P1 ...P4 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:
SP = P'n = 0 with 8n e W (A.3)
where W is the space of kinematically admissible displacements, SP is the first variation of the potential P with respect to the arbitrary field n and the notation ()' denotes the variational derivative. Evaluating Eq. (A.3) using the Gateaux derivative, we write the equilibrium equation:
d:P[u+tn] = Pi [n] + P11 [u,n] + P21 [u,n] + P31 [u,n] = 0 (A.4)
where t is an arbitrary scalar, and the subscripts of the functionals refer to the powers of each their arguments in a binomial expansion. For example, P2[u1 + u2] = P2[u1]+P11[u1, u2]+P2[u2]. In order to simplify the notation in the proceeding analysis we introduce the following auxiliary variables:
S1[u] = P„[u, n], S2[u] = P21 [u, n] and S3[u] = P31 [u, n]
which are implicit functions of n. They are given by the following expressions:
51 [u; n] = i [e(u)T • A • e(n)+k(u)T • D • k(n)] dA
52 [u; n] = 2 f [e(n)T • A • E(u)+e(u)T • A • E„ (u, n)] dA
53 [u; n] = / [E(u)T • A • E11 (u, n)] dA
En[ u1, u2] =
dw1 dw2 dw1 dw2 / dw1 dw2 dw1 dw2\ dx dx dy dy V dx dy dy dx )
In the pre-buckling state it is assumed that rotations dw0/dx and dw0/dy are sufficiently small that the stresses N0 are approximately equal to A • e0. That is, the higher-order strains E0 are assumed to be negligible. Substituting the displacements
u = u0+v with u0 = Au0 (A. 5)
(where v is an arbitrary displacement field, u0 is the pre-buckling mode and A is the loading parameter) into Eq. (A.4) yields the following:
P1 [n]+S1 [uo] + S1 [v] + S2 [uo] + S11 [uo, v]+S2 [v] + S3 [uo] + S21 [uo, v] + S12 [uo, v]+S3 [v] = 0 (A. 6)
The small-rotations assumption facilitates a considerable simplification of the analysis. At this stage we may assume that
S2[uo]~0, S3[uo]~0 and S21 [uo, 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 (Hilburger et al., 2001), membrane constraints at the boundary (Hilburger et al., 2004) or circular perforations (Hilburger et al., 2001) 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 u0 by setting v to zero. This gives the following linear variational problem:
P1 [n]+S1 [uo] = 0, (A. 7)
Integrating Eq. (A.7) by parts in order to isolate the zeroth-order derivatives of n and noting that n 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).
Assuming that the pre-buckling displacements are known, we may now formulate a problem for v. Simplifying Eq. (A.6) we have
S1 [v] + S11 [uo, v] + S2 [v] + S12 [uo, v] + S3 [v] = 0
(A. 8)
The type of buckling assumed in this analysis is a bifurcation, that is, the structure's stability changes because a compatible adjacent state is met at a critical value of the loading parameter X = Xc. We must, therefore, determine a displacement vector au1 such that Eq. (A.8) is satisfied as a —0. Substituting v = au1 we find
Si [ui]+S11 [uo, ui ]+a(S2 [ui ]+S12 [uo, Ui, n]) + O(a2) = 0 (A.9)
which implies that
Si[ui]=-XcS„[u o, ui] (A. 10)
This is a linear eigenvalue problem for the field u1 and the load Xc. It may be simplified by noting that in S11 [u0, u1] the biggest term by far is
[ (N0 • E„[ui, n]) 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 Xc.
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:
Ri[v]=Si[v]+Sii[u0, v], R2[v] = S2[v]+Si2[u0, v], and R¡[v] = S3[v]
which are implicit functions of n and u0. Now, we add a quadratic correction field u2 to v to give
v = au1 + a2 u2
Substituting this into Eq. (A.8) gives the following expression:
Ri [ui ] + a(Ri [u2] +R2 [ui ]) + O(a2) = 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 u2:
Ri[u2]=-R2 [ui] (A. 12)
In addition to Eq. (A.12) further constraints on u2 are required in order to render them unique. This is because at X= Xc 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
F2[u]= -Pi2[u 0, u]= -S2[u; u 0]
is positive-definite. It, therefore, holds that
F?[u+tw] > F2 [u] with w a W
F11 [u, w] + tF2 [w] > 0 (A. 13)
We define fields u and w as orthogonal if they satisfy the equation:
Fii[u, w] = 0, (A. 14)
meaning that the components of energy derived from the two fields may be superposed (i.e. F2[u+w] =F2[u]+F2[w]).
The proof that Eq. (A.12) has a unique solution as long as u1 and u2 are orthogonal is given by van der Heijden (2009). It is shown that the solution for u2 is independent of the actual choice for the functional used to enforce orthogonality. The only requirements being that it is positive-definite.
Applying Gauss' theorem to Eq. (A.12) we find that Eqs. (7) and (9), and either (14) (or (10) and (16)) must hold for equilibrium. Expanding Eq. (A.13) we also find that Eq. (8) must hold for orthogonality.
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
PX = P[Xu 0 + v] - P[Xu 0] =P2[v]+XS2[v; Ü0]+P3[v]+P4[v] + O(X2) (A. 15)
The same expression may be found by calculating a Taylor series for P in terms of X, about the point X = Xc (Budiansky, 1974). According to Koiter (1956), PX should be minimised with respect to the perturbation parameter a for equilibrium. Substituting v = au1 + a2u2 into Eq. (A.15) gives
PX = a2 (P2 [ui ]+XS2 [ui; u 0]) + a3 (P„ [ui, u2]+P3 [ui ]) + a4 (P4 [ui ]
+P21 [ui, u2]+P2 [u2] + XS2 [u2; u 0]) + O(a5) (A. 16)
Differentiating Eq. (A.16) with respect to a, substituting in the series in Eq. (3) and noting that P11[u1, u2] equals zero, we may write the following:
P2 [U1 ] + ÂcS2 [U1 ; u 0] + a (2 P3 [u ]+h S2 [U1 ; u 0] ) + • -
2a2 (?4[U1] +P21 [U1,U2]+P2[U2] +hcS2 [U2;û0] +1^2[U1 ;û0^ + O(a3) = 0 (A. 17)
Requiring that each of the coefficients of a in this series are equal to zero gives the following expressions for the coefficients in Eq. (3):
hc = - P2[u1] , (A. 18a)
c S2 [U1; it 0]' ^ 7
. 3 P3[u1] rA18b
h = "2 S^^' (A:18b)
, 2P2 [U2] + P4[U1 ] + P21 [U1 ; U2] + hcS2 [U2 ; it0] (A 18c)
2 S2 [U1 ; it 0] 1 ■ c
Expanding the functionals in these expressions gives Eqs. (17a-c).
Appendix B. Matrix notation
The o 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,
[a b]c
a1 c12
a2 c22
References
Batdorf, S.B., 1947, A Simplified Method of Elastic-Stability Analysis for Thin Cylindrical Shells. Technical Report, NACA, Langley Field, VA. Bauld Jr., N.R., Khot, N.S., 1982. A numerical and experimental investigation of the buckling behavior of composite panels. Comput. Struct. 15 (4), 393-403. Bellman, R., Kashef, B.G., Casti, J., 1972. Differential quadrature: a technique for the rapid solution of non-linear partial differential equations. J. Comput. Phys. 10 (1), 40-52.
Blom, A.W., Setoodeh, S., Hol, J.M., Gurdal, Z., 2008. Design of variable-stiffness conical shells for maximum fundamental eigenfrequency. Comput. Struct. 86 (9), 870-878.
Budiansky, B., 1974. Theory of buckling and post-buckling behavior of elastic structures. Adv. Appl. Mech. 14, 1-65.
Casciaro, R., Garcea, G., Attanasio, G., Giordano, F., 1998. Perturbation approach to elastic post-buckling analysis. Comput. Struct. 66 (5), 585-595.
Chandra, R., Raju, B.B., 1973. Postbuckling analysis of rectangular orthotropic plates. Int. J. Mech. Sci. 15 (1), 81-97.
Cox, H.L., Clenshaw, W.J., 1941. Compression Tests on Curved Plates of Thin Sheet Duralumin. A.R.C. R.&M. No. 1849.
Cox, H.L., Pribram, E., 1948. The elements of the buckling of curved plates. J. R. Aeronaut. Sci. 52 (453), 551-565.
Crisfield, M.A., 1981. A fast incremental/iterative solution procedure that handles 'snap-through'. Comput. Struct. 13 (1), 55-62.
Donnell, L.H. 1933. Stability of Thin-Walled Tubes Under Torsion (No. NACA-R-479). California Inst. of Tech., Pasadena (USA).
El Naschie, M.S., 1991. Stress, Stability and Chaos in Structural Engineering: An Energy Approach. McGraw-Hill, New York.
Garcea, G., Madeo, A., Zagari, G., Casciaro, R., 2009. Asymptotic post-buckling FEM analysis using corotational formulation. Int. J. Solids Struct. 46 (2), 377-397.
Gerard, G., Becker, H., 1957. Handbook of Structural Stability, Part III - Buckling of Curved Plates and Shells. N.A.C.A. T.N. 3783.
Gurdal, Z., Tatting, B.F., Wu, C.K., 2008. Variable stiffness composite panels: effects of stiffness variation on the in-plane and buckling response. Compos.
Part A: Appl. Sci. Manuf. 39 (5), 911-922. Gurdal, Z., Olmedo, R., 1993. In-plane response of laminates with spatially varying fiber orientations: variable stiffness concept. AIAA J. 31 (4), 751-758. Hilburger, M.W., Britt, V.O., Nemeth, M.P., 1999. Buckling behaviour of compression-loaded quasi-isotropic curved panels with a cutout. In: 40th AIAA/
ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, April 12-15. Hilburger, M.W., Nemeth, M.P., Starnes, Jr., J.H., 2001. Nonlinear and buckling behavior of curved panels subjected to combined loads. In: 42nd AIAA/ASME/
ASCE/AHS Structures, Structural Dynamics and Materials Conference and Exhibit, 16-19 April, Seattle, WA. Hilburger, M.W., Britt, V.O., Nemeth, M.P., 2001. Buckling behaviour of compression-loaded quasi-isotropic curved panels with a circular cutout. Int. J. Solids Struct. 38, 1495-1522.
Hilburger, M.W., Nemeth, M.P., Riddick, J.C., Thornburgh, R.P., 2004. Effects of elastic edge restraints and initial prestress on the buckling response of compression-loaded composite panels. In: Presented at the 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, April 2004.
Hyer, M.W., Lee, H.H., 1991. The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes. Compos. Struct. 18, 239-261.
Jackson, K.B., Hall, A.H., 1947. Curved Plates in Compression. National Research Council of Canada, Aeronautical Report No. AR-1.
Kheyrkhahan, M., Peek, R., 1999. Postbuckling analysis and imperfection sensitivity of general shells by the finite element method. Int. J. Solids Struct. 36 (18), 2641-2681.
Kim, B.C., Potter, K., Weaver, P.M., 2012. Continuous tow shearing for manufacturing variable angle tow composites. Compos. Part A: Appl. Sci. Manuf. 43 (8), 1347-1356.
Kirchhoff, G., 1850. Uber das Gleichgewicht und die Bewegung einer elastischen Scheibe. J. die reine angew. Math. 40, 51-88.
Koiter, W.T., 1956. Buckling and post-buckling behaviour of a cylindrical panel under axial compression. Versl. Verh. (Report S. 476) 20, 71-84.
Koiter, W.T., 1945. On the Stability of Elastic Equilibrium. Ph.D. Techische Hooge School, Delft.
Lopes, C.S., Camanho, P.P., Gurdal, Z., Tatting, B.F., 2007. Progressive failure analysis of tow-placed, variable-stiffness composite panels. Int. J. Solids Struct. 44 (25), 8493-8516.
Madenci, E., Barut, A., 1994. Pre- and post-buckling of curved, thin, composite panels with cutouts under compression. Int. J. Numer. Methods Eng. 37, 1499-1510.
Panda, S.K., Ramachandra, L.S., 2010. Postbuckling analysis of cross-ply laminated cylindrical shell panels under parabolic mechanical edge loading. Thin-Walled Struct. 48 (8), 660-667.
Penrose, R., 1955. A generalized inverse for matrices. In: Proceedings of Cambridge Philosophical Society, vol. 51, No. 3, pp. 406-413.
Peterson, J.P., 1969. Buckling of Stiffened Cylinders in Axial Compression and Bending: A Review of Test Data. National Aeronautics and Space Administration.
Pope, G.G., 1965. On the Axial Compression of Long, Slightly Curved Panels. British A.R.C., R. & M. No. 3392.
Pope, G.G., 1968. The buckling behaviour in axial compression of slightly-curved panels including the effects of shear deformability. Int. J. Solids Struct. 4, 323-340.
Rahman, T., Ijsselmuiden, S.T., Abdalla, M.M., Jansen, E.L., 2011. Postbuckling analysis of variable stiffness composite plates using a finite element-based perturbation method. Int. J. Struct. Stab. Dyn. 11 (04), 735-753.
Raju, G., Wu, Z., Kim, B.C., Weaver, P.M., 2012. Prebuckling and buckling analysis of variable angle tow plates with general boundary conditions. Compos. Struct. 94 (9), 2961-2970.
Raju, G., Wu, Z., Weaver, P. M., 2013. Postbuckling analysis of variable angle tow plates using differential quadrature method. Compos. Struct. 106, 74-84, http://dx.doi.org/10.1016lj.compstruct.2013.05.010.
Reddy, J.N., 1984. Energy and Variational Methods in Applied Mechanics. Wiley, New York396-403.
Riks, E., 1979. An incremental approach to the solution of snapping and buckling problems. Int. J. Solids Struct. 15 (7), 529-551.
Sherbourne, A.N., Pandey, M.D., 1991. Differential quadrature method in the buckling analysis of beams and composite plates. Comput. Struct. 40 (4), 903-913.
Shu, C., 2000. Differential Quadrature and Its Application in Engineering. Springer ISBN: 1-85233-209-3.
Shu, C., Du, H., 1997. A generalized approach for implementing general boundary conditions in the GDQ free vibration analysis of plates. Int. J. Solids Struct. 34 (7), 837-846.
Stein, M., 1959. Loads and Deformations of Buckled Rectangular Plates. Technical Report No. R-40, NASA.
Thornburgh, R.P., Hilburger, M.W., 2005. Identifying and Characterizing Discrepancies Between Test and Analysis Results of Compression-Loaded Panels. Technical Report, October, NASA.
Thornburgh, R.P., Hilburger, M.W., 2006, A numerical and experimental study of compression-loaded composite panels with cutouts. In: 47th AIAA/ASME/ ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 1 -28, Newport RI, AIAA.
van der Heijden, A.M. (Ed.), 2009. WT Koiter's Elastic Stability of Solids and Structures. Cambridge University Press, Cambridge.
von Karman, T., Tsien, H.-S., 1941. The buckling of thin cylindrical shells under axial compression. J. Spacecr. Rockets 40 (6).
Wang, X., Tan, M., Zhou, Y., 2003. Buckling analyses of anisotropic plates and isotropic skew plates by the new version differential quadrature method. Thin-Walled Struct. 41 (1), 15-29.
Wu, K.C., 2008, Design and analysis of tow-steered composite shells using fiber placement. In: American Society for Composites. In: 23rd Annual Technical Conference, Memphis, TN, September 9-11, 2008.
Wu, Z., Raju, G., Weaver, P.M., 2013a. Postbuckling analysis of variable angle tow composite plates. Int. J. Solids Struct. 50 (10), 1770-1780.
Wu, Z., Weaver, P.M., Raju, G., 2013b. Postbuckling optimization of variable angle tow composite plates. Compos. Struct. 103, 34-42.
Yamaki, N., 1984. Elastic Stability of Circular Cylindrical Shells. North-Holland Publishing, Amsterdam.
Zeng, H., Bert, C.W., 2001. A differential quadrature analysis of vibration for rectangular stiffened plates. J. Sound Vib. 241 (2), 247-252.
Zhang, Y., Matthews, F.L., 1983. Initial buckling of curved panels of generally layered composite materials. Compos. Struct. 1 (1), 3-30.
