Scholarly article on topic 'Buckling analysis of variable angle tow, variable thickness panels with transverse shear effects'

Buckling analysis of variable angle tow, variable thickness panels with transverse shear effects Academic research paper on "Civil engineering"

Share paper
Academic journal
Composite Structures
OECD Field of science
{"Variable stiffness laminates" / Buckling / "Transverse shear" / "Differential quadrature"}

Abstract of research paper on Civil engineering, author of scientific article — R.M.J. Groh, P.M. Weaver

Abstract The manufacture of advanced composite panels with variable fibre angles can lead to laminates with a flat profile on one side and a smooth, curved profile on the other. When modelling these laminates in two-dimensional form the flat plate assumptions may no longer accurately capture the structural behaviour. In this paper the buckling behaviour of laminates with one-dimensional fibre variations and symmetric stacking sequences is investigated. The assumptions of modelling the three-dimensional profile as a flat plate or a cylindrical panel are assessed, taking into account the effects of transverse shear deformation. The governing differential equations are solved in the strong form using the Differential Quadrature method and validated by 2D finite element models. The validity of the two modelling approaches is assessed by comparing the solutions to a 3D finite element model capturing the actual shape of the laminate. It is suggested that the buckling event of these variable angle tow, variable thickness laminates is characterised more accurately by “shell-like” than by “plate-like” behaviour. The idea of investigating the effects of two-dimensional fibre orientations with their associated doubly curved topologies is proposed.

Academic research paper on topic "Buckling analysis of variable angle tow, variable thickness panels with transverse shear effects"

Composite Structures 107 (2014) 482-493

Contents lists available at ScienceDirect

Composite Structures

journal homepage:

Buckling analysis of variable angle tow, variable thickness panels with transverse shear effects q

R.M.J. Groh *, P.M. Weaver

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


The manufacture of advanced composite panels with variable fibre angles can lead to laminates with a flat profile on one side and a smooth, curved profile on the other. When modelling these laminates in two-dimensional form the flat plate assumptions may no longer accurately capture the structural behaviour. In this paper the buckling behaviour of laminates with one-dimensional fibre variations and symmetric stacking sequences is investigated. The assumptions of modelling the three-dimensional profile as a flat plate or a cylindrical panel are assessed, taking into account the effects of transverse shear deformation. The governing differential equations are solved in the strong form using the Differential Quadrature method and validated by 2D finite element models. The validity of the two modelling approaches is assessed by comparing the solutions to a 3D finite element model capturing the actual shape of the laminate. It is suggested that the buckling event of these variable angle tow, variable thickness laminates is characterised more accurately by "shell-like" than by "plate-like" behaviour. The idea of investigating the effects of two-dimensional fibre orientations with their associated doubly curved topologies is proposed.

© 2013 The Authors. Published by Elsevier Ltd. All rights reserved.


Article history:

Available online 29 August 2013


Variable stiffness laminates Buckling Transverse shear Differential quadrature

1. Introduction

The idea of tailoring the structural performance of composite laminates by spatially varying the point-wise fibre orientations over the planform has been explored since the early 1990s. For example, the work by Hyerand Lee [1] and Hyerand Charette [2] showed that such variable angle tow (VAT) laminates can improve the stress concentration around holes by arranging the fibres in the direction of critical load paths. In recent years the use of fibre re-inforced composites in primary aircraft structures has led to increased interest in VAT technology. Numerous works have shown that tailoring the in-plane stiffness of a plate allows pre-buckling stresses to be re-distributed to supported regions and thereby improve the buckling behaviour [3-7]. In this manner VAT technology has been shown to improve the buckling performance of a composite fuselage window section by 12% compared to an equivalent straight-fibre laminate [8] and alleviate the pressure pillowing of fuselage sections [9]. Recent results by Wuetal. [10] show that VAT plates with linear fibre variations can be designed to exhibit smaller stiffness reductions in the post-buckling regime than their straight-fibre counterparts. Besides, the optimum fibre orientations for increasing the

q This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial-No Derivative Works License, which permits non-commercial use, distribution, and reproduction in any medium, provided the original author and source are credited.

* Corresponding author. Tel.: +44 0 117 33 15793. E-mail address: (R.M.J. Groh).

buckling load are similar to those for minimising the transverse displacement in the post-buckling regime [11].

Currently, the major technology for manufacturing VAT laminates is the automated fibre placement (AFP) technique developed since the 1980s. AFP uses a robotic fibre placement head that deposits multiple pre-impregnated tows of ''slit-tape'' allowing cutting, clamping and restarting of every single tow. While the robotic head follows a specific fibre path, tows are heated shortly before deposition and then compacted onto the substrate using a special roller. Due to the high fidelity of current robot technology AFP machines can provide high productivity and handle complex geometries [12]. However, in AFP steering is accomplished by bending the tows of fibres in-plane, which leads to local fibre buckling on the inside radii of the curved tow and thus limits the steering radius of curvature [13]. Furthermore, if individual tows are placed next to each other by shifting the reference path along a specific direction, tow gaps and overlaps are inevitably required to cover the whole surface. To overcome the drawbacks of AFP the continuous tow shearing technique (CTS) was developed which uses shear deformation to steer fibres at the point of application [14]. This technique not only allows much tighter radii of curvature but tow gaps and overlaps are also avoided by tessalating tows on the substrate [14]. One feature of CTS is that in order to maintain the volume fraction of fibre the thickness of a tow inherently increases as it is sheared. The relation between unsheared tow thickness t0 and sheared tow thickness th is

to cos 8

0263-8223/$ - see front matter © 2013 The Authors. Published by Elsevier Ltd. All rights reserved.

where h is the shearing angle of the tow. Consequently, the thickness of a ply may locally increase by a factor of 4 if the fibre tow is sheared through an angle of 75°. As the laminate is cured on a tool plate one side of the laminate remains flat while the other resembles a curved panel as depicted in Fig. 1. The effects of this onesided thickness profile in terms of local three-dimensional stress fields or buckling behaviour is currently unexplored.

Furthermore, since composite laminates are more affected by transverse shear effects than isotropic materials [15,16] these local areas of increased thickness make predictions of the buckling behaviour using Classical Laminate Analysis (CLA) [17] overly conservative. Buckling optimisation algorithms using the Finite Element method (FEM) can become computationally expensive, and do not readily shed physical insight into the fundamental mechanisms in which these thickness variations could influence the buckling behaviour. In fact, there is an inherent need in the aerospace industry for rapid yet accurate analysis tools to aid in preliminary design studies. For this reason a reduced 2D equivalent single-layer formulation for the flexural behaviour of VAT plates incorporating transverse shear effects was developed [18]. The model accurately predicts global bending and buckling behaviours up to thickness to length ratios of 1:10 while restricting the number of unknowns to only four variables to enhance computational efficiency. It was shown that CLA predictions result in discrepancies in excess of 15% for thickness to length ratios of 1:20. Thus, considering the extent of thickness variation possible using CTS, transverse shear effects need to be incorporated if the spatial fibre distribution is to be optimised accurately for buckling and post-buckling behaviour.

In this work, the elastic stability of VAT panels with linear fibre variations manufactured using the CTS technique is analysed using two different approaches. The aim is to compress the three-dimensional shape of the laminate as seen in Fig. 1 onto an equivalent single layer. However, due to the differences in topology between the top and bottom surfaces the shape of the chosen reference

plane is not immediately obvious. First, the VAT panel is assumed to behave as a flat, non-symmetric plate and the aforementioned 2D equivalent single-layer formulation is applied to solve the buckling problem. Second, and the major novelty of this work, is to assume that the smooth thickness variations cause the structure to behave more like a curved panel, in which case the equilibrium equations of cylindrical shells must be used.

This premise creates the opportunity for interesting optimisation studies since non-prismatic fibre variations cause the curved side of the VAT panel to take on more complex, doubly-curved shapes. In this case the need for using the most general set of shell kinematics may reveal different and perhaps enhanced buckling solutions compared to previous surveys based on flat plates [5,6,8,9]. Furthermore, the relatively shallow curvature produced by the thickness variations suggests that the unstable nature of shell post-buckling paths may be contained to a large extent. Thus, the idea of a structure with higher "shell-like" buckling loads in conjunction with stable "plate-like" characteristics in the post-buckling regime is suggested. For the purposes of this paper we are solely interested in solving the buckling problem of VAT panels with prismatic fibre variations that cause cylindrical shapes of the reference surface. Buckling and post-buckling studies of more general VAT panels as discussed above will be the topic of future work.

A large number of studies regarding the solution of von Kármán's non-linear differential equations for the postbuckling behaviour of plates using either energy methods [19-22] or Fourier-series expansions [23,24] were published in the early to mid-part of the 20th century. In the present work, the perturbation technique demonstrated by Stein [25] for rectangular isotropic plates and later developed by Chandra and Raju [26] for symmetrically laminated, orthotropic plates is applied. In this technique the non-linear differential equations are reduced to an infinite set of linear differential equations by expanding the unknown functional fields in a power series of an arbitrary perturbation parameter about the point of buckling. The solution of the first equation of the infinite set provides the pre-buckling stresses in the structure, while the eigenvalues and eigenmodes of the second equation yield the buckling loads and modes. Finally, succeeding equations beyond this point relate to large-deflection postbuckling solutions.

The rest of the paper is laid out as follows. Section 2 briefly describes the theoretical background of the aforementioned 2D equivalent single-layer model and outlines the governing equations needed to model the laminate as a shell. Furthermore, the perturbation technique for solving the ensuing non-linear shell equations is presented in more detail. In Section 3, the governing equations are solved in the strong form for various load and boundary conditions. The analytical results are validated by 2D FEM and 3D FEM techniques and the accuracy of the "plate-like" and "shell-like" modelling approaches discussed over a range of realistic fibre trajectories. Finally, conclusions are drawn in Section 4.

2. Theory

2.1. Modelling of variable thickness cross-section

Linear fibre variations in one direction only (i.e prismatic variation) can conveniently be defined by the notation Q = /(T0\T1), where / denotes the rotation of the fibre path with respect to the x-axis, while T0 and T1 are the fibre angles at the ply centre and a characteristic length d from the centre respectively [3]. To fill the planform the fibre trajectories are then shifted perpendicular to the steering direction /. A 0{0\ -70) VAT layer is drawn schematically in Fig. 1.

Tool Surface

Fig. 1. CTS manufactured 0{0|-70) layer with reference path shifted in the Y-direction showing fibre orientations and thickness variations. The 3D structure can be compressed onto an equivalent single layer resembling either a cylindrical shell (Curve 1) or a non-symmetric flat plate (Curve 2).

The thickness variation produced by CTS means that the VAT laminate now resembles a flat plate on one side and a cylindrical shell of varying radius of curvature on the other. If this 3D structure is compressed onto an equivalent single-layer coincident with the mid-plane of the thickness at each point, the reference layer will resemble a cylindrical shell (Curve 1 in Fig. 1). Alternatively the hypothesis is made that the VAT laminate continues to behave as a flat plate with the reference surface coincident with the mid-plane of the nominal (unsheared) laminate thickness (Curve 2 in Fig. 1). Thus, the equivalent single-layer is assumed to behave as a non-symmetric laminate in areas of increased thickness even if the nominal stacking sequence is symmetric. In this work we do not attempt to solve the coupled stretching-bending problem of the non-symmetric plate but use the more straightforward, albeit less accurate technique, of substituting the reduced stiffness matrices A* = -A - BD-1B and D* = D - BA-1B [27], and reduced flexibility matrix d* = D* into the equations for a symmetric laminate derived in Section 2.2. The applicability of these two modelling approaches is validated by comparing the results to 3D FEM solutions.

2.2. Equivalent single-layer plate model with transverse shear effects

The equivalent single-layer plate model of the 3D structure, as depicted in Fig. 2, including the effects of transverse shear deformation is derived in detail in [18]. The linear variations of in-plane stresses ~k = (rXk ryk rxyk )T through the thickness of each ply with bending moments M = (Mx My Mxy)T as defined in CLA, are retained. For a symmetric laminate the stresses in each ply k are

ak = QkzdM (2)

where Q^ is the transformed stiffness matrix of each ply k and d is the bending flexibility matrix of the complete laminate.

Thus, no formal assumption of the displacement variables u, v, w is made and the curvatures Kx, Ky, Kxy are not explicitly related to the second partial derivatives of w. The transverse shear stresses are derived by integrating the in-plane stresses Eq. (2) in Cauchy's in-plane equilibrium equations in the thickness z-direction. The constants of integration are eliminated by enforcing boundary conditions of shear stress continuity at each ply interface and vanishing transverse shear stresses on the top and bottom surfaces. Consequently, the transverse shear stress profile is defined to be quadratic in each ply a priori such that a shear correction factor, as in Mindlin's FSDT [28], is not required. Thus for a symmetric laminate the transverse shear stresses sk = (sxz syz)T for a particular layer k is given by [29],

Sk = -CBk(x,y,z) dM (3)

Bk(x,y,z)=-Qk(z2 - z2_1) + Qj(z2 - z2_1

/ 9x 9y \

( 0 9 9

\ 9y 9x /

By observing Eq. (4) it can be seen that Bk defines the partial in-plane/out-of-plane coupling stiffness matrix up to the kth ply under consideration. This approach of capturing both the in-plane and transverse shear stresses with just the three moments Mx, My and Mxy improves the computational efficiency compared to many other enhanced plate models that have been published in recent decades. For example Higher-Order Shear Deformation Theories (HSDT) extend the assumption on the deformation field (u, v,w) from linear to a cubic variation in the through-thickness co-ordinate z, thereby removing the need for a shear correction factor and improving the accuracy of the tranverse shear strains. A large number of polynomial [30-32], trigonometric [33-37], hyperbolic [38,39] and exponential [40,41] shear functions have been investigated in recent years. However, these theories assume between 3 (FSDT: w0 and two rotations) to 7 (Third Order SDT: w0 and six rotations) unknown variables to represent the deformation field. To capture the discontinuous, Zig-Zag intra-ply deformations this approach has also been extended to Layer-Wise theories by assuming independent displacement fields for each layer. However, by being formulated on a displacement assumption these theories cannot a priori guarantee interlaminar continuity for the transverse stresses. In this respect Murakami's equivalent single-layer approach [42] based on Reissner's Mixed Variational Theorem assuming two independent fields, one for displacements and the other for transverse stresses, elegantly guarantees transverse continuity and Zig-Zag behaviours a priori. While the equivalent single-layer model in this paper cannot capture Zig-Zag effects nor accurate local 3D stresses for very thick laminates, it is similar in accuracy to HSDT for moderately thick laminates with the added benefit of guaranteeing interlaminar continuity and reduced number of variables.

Equilibrium equations for a VAT panel in pure bending are derived using the Hellinger-Reissner mixed variational principle. By definition the total potential energy of the 3D system is given by,

n = 2 (rT e + ST c + rzez)dx dy dz -(run + TnSus + T^w)ds dz (5)

The variation of the functional n must vanish in such a manner that transverse equilibrium is satisfied. The shear stress resultants Qx and Qy are related to the transverse load intensity p by the conventional transverse equilibrium equation,

9@t++p - 0.

Following Reissner's mixed variational principle [43] Eq. (6) is added to the variation of the functional P using the Lagrange multiplier k(x,y). Replacing the strains in Eq. (5) with constitutive relations we find,

»{2 Z0

lx rly rZn / __ , __ , .

<5^ / (~T Q -1r + sT G-1t) dz dy dx-

J0 J0 Jz0 ^

(r nUn + TnsUs + T„zW)ds dz+

Jz0 Jr rlx r h

i i" <91+f+p) dy * - 0.

Fig. 2. 3D structure condensed to a single layer coincident with the mid-plane (grey).

Note, that in Eq. (7) we have neglected the effects of normal transverse stress rz. This is in contradiction of Koiter's Recommendation

that a refinement of Love's first approximation theory should take into account both transverse shear and normal stresses [44]. However, in the present work we are primarily interested in global buckling effects rather than accurate localised ply-by-ply stresses. Furthermore, as representative of common aerospace panels, the thickness to length ratios of the analysed panels do not exceed 1:35. Under these circumstances the effects of transverse normal stresses on the global elastic behaviour is considered to be negligible compared to the five other stress components (rx, ry, axy, sxz,

After substituting Eq. (2) into Eq. (7) the variations with respect to the unknown fields M and Lagrange multiplier k(x,y) are performed using the calculus of variations. Thus,

SMXy : (di6 + a3i)Mx + (£¡26 + «32)My + (de6 + «33 )Mxy + o^

dMx x dMy x -Mxy y dMx y dMy X dx + «32 dx + «33 dx + «31 dy + «32 dy

, ay @M*y + axx d_Mx. axx d_My , axx @2Mx^+ ayy

_x I ayy d lv'y ayy _i_ «L «y ^"y

@y2 + «32 @y2 + «33 @y2 + «31 @xÔy + «32 dxdy

«xy d2Mxy

33 dxdy

d2w dxdy

dM (dM )-

aM + a

;d_M dy

, d2M dxdy

+ *X (dMx)-

2 -*-» (dMxy)

d2X dy2

: (dMy)

*L , ^ + p )(dk)]dydx

@x@y @y

bM + bx f + , §) dM-

(X - W)dQn-

nx ■

u'x dMx - < n

@X HTy-

ny{ — + u0, I dM,

where nx and ny are the directional cosines of the tangent vector to the boundary curve C, and a and b terms are 3 x 3 matrices that represent the shear coefficients of the present theory. To guarantee that dP = 0 the coefficients of the variations of M, Qn and k must equate to zero in both the double and the line integrals. In this manner, the governing field equations are derived from the double integrals and the pertinent boundary conditions from the line integrals,

SMX: (d„ + an)Mx dMx X dx~' ~12 dx

+ a?3 dM*+ax @2Mx+«s d2M+«»@2M

(£¡12 + ai2)My + (¡¡16 + <Xi3)Mxy

x @ My x @ Mxy

xx d J",

<12 dx2 2

I ^^ <y @ml

II dy + <12 dy

x i <ll @_ml , <ll d_My axy + «12 «„2 + <13 »„2 + <11

+ «x1y3

d2Mx „ d2My + <12

13 dxdy @_w dx_

dMy : (¡12 + a_1)Mx + (£¡22 + a__)My + (£¡26 + a_3)Mxy + a_1 dMr

. d2My

dMx, dx

+ <22 ^r + «_3^zt + «2^ + «__-y + «_3 ^ + <_1

+ + «23

«__ dxdy

+ «2L

x2 d_My

xy + «y2y1

— + «__ S + «53

dxdy '

r + «__

dMy dM,

d_Mx dx2

y + «_3-

<* + 2d_Mxy

Thus, the flexural behaviour of the plate for a specific transverse loading p is described by just four unknown fields Mx, My, Mxy and w that are related by the four differential Eq. (9). Note that the effects of transverse shear deformation are simply incorporated into the present model by modifiying the classic constitutive relations of CLA with higher order derivatives of the three moments. Thus Eqs. 9a, 9c and 9d include differential terms of Mx, My and Mxy multiplied by various a coefficients that can be interpreted as transverse shear flexibilities.

The displacement boundary conditions are derived from the coefficients of the line integral,

dMx : fiuMx + b1_My + b13Mxy + bxn ddxx + b1_ ^ + b13

dMxy by @Mx by dMy by dMxy dx +b11 dy +b12 dy +b13 dy '

Mx x My

dMy : b_M + b__My + b_3Mxy + b_1 -dxx + b_2 ^xL + b_3

@mxl+b_21 dMx+b_2 -m-+b_ -Mi.

dx 1 b21 d^ """b22 dy b23 dy y dy

dMxy : b31Mx + b3_My + b33 Mxy + bx31 @M + b3_ -Mx- + bX33

d-Mi+b31 -Mi+b3_ -M-+b3 -Mi-

-nx — -nv

-dx 1 b31 dy b32 dy b33 dL dy y w

= nxuy + n-ux

(10c) (10d)

dQn : X = w = w

which have to be applied if the force boundary conditions Mx = Mx, ML = ML, Mxy = Mxy or Qn = Qn are not specified.

2.3. Equivalent single-layer plate buckling problem

Since the in-plane stiffness matrix A varies spatially from point to point for a VAT laminate, even constant applied edge loads or displacements result in non-uniform in-plane stress distributions rx(x,L), ry(x,y) and rxy(x,y). It is important to accurately determine the in-plane stress distribution as an input to the buckling problem. By observation of Cauchy's in-plane equilibrium equations,

drx dffxy dQxz _ o

dx + dy + dz =

= nyUy

dffxy dffv dr

the varying stress distributions rx(x,y), ry(x,y) and rxy(x,y) for a VAT plate under pure in-plane loading may require transverse shear stresses to satisfy equilibrium. However, in the present work the effect of transverse shear stresses on the pre-buckling problem is assumed to be negligible such that any non-zero drx/dx and dry/dy are assumed to be fully equilibrated by drxy/dy and drxy/dx respectively. Substituting the constitutive relations for the in-plane load resultants given by,

„1 A du dv . fdu dv\ ... .

Nx = A11 du + A12 dy + + dx) (11a)

.du . dv . fdu dv

NV = A12 @u + A22 dV + A26{ (Ту + dV

du dv du dv Nxy — Ai6 + A26 ^ + A66( + KV)

where u and v are the mid-plane deformations in the defined x- and y-directions, into the in-plane equilibrium equations written in terms of the stress resultants

dx + dy

X dNxy_ 0 and

dNxy 0

dx + dy


we get

A11u,xx + 2A16u ,xy + A66u,yy + A16 v,xx + (A12 + A66) v,xy

+ A26V yy + (An,x + A16v)u,x + (A16,x + A66v)u,v + (A16,x + A66;v)V,x + (A12,x + A26y ) Vj, = 0

A16u,xx + (A12 + A66)u,xv + A26u,yy + A66 v,xx + 2A26 v,xy + A22V yy + (A16,x + A12y)ux + (A66,x + A26v)uy + (A66,x + A26;v)V ,x + (A26,x + A22y )V j, 0

which can be solved for u and v by prescribing the boundary conditions u = U or Nxy = Nxy on x; u = U or Nx = Nx on y; v = V or Ny = Ny on x and v = V or Nxy = Nxy on y. Once u and v are known under the prescribed boundary conditions the in-plane stress resultants N0, N0, N0y are found using Eq. (11). Finally, the transverse pressure p in Eq. (9d) is replaced by the out-of-plane components of the in-plane stress resultants,

d2Mx ~dxr

xv , d My _ M0 @2w , n M0 @2w

dxdy dy


x dx xy dxdy

ю @2w

such that the eigenvalue problem of Eq. (9) can be solved.

2.4. Shallow shell model using first order shear deformation theory

To model the VAT laminate as a moderately thick, 2D single-layer shell as depicted in Fig. 3 we develop a formulation based

Fig. 3. Schematic representation and curvilinear co-ordinate system ({1, {2, f) of a doubly-curved shell.

on the First Order Shear Deformation Theory (FSDT). In this respect the following four assumptions are made [45],

1. Transverse normals to the reference surface of the shell before deformation are inextensible.

2. Normals to the reference surface do not necessarily remain normal after deformation i.e. assumptions on the rotations of the cross-section with respect to the undeformed state, /1 and /2, are relaxed such /1 — - dw/af1 and /1 — - dwjd{2.

3. The transverse normal stress may be neglected in comparison with the stresses acting in the direction parallel to the reference surface.

4. All strains are sufficiently small i.e. e < 1 and Hooke's Law applies.

Under these circumstances the FSDT displacement field in terms of the orthogonal curvilinear co-ordinates ({1, {2,f) shown in Fig. 3 is assumed to be

u(nb ¿2, 0=u0(¿1, ¿2, C)+C/1(f1, ¿2, f)

v(¿1, ¿2, f) — V0(f1, ¿2, f) + f/2(Î1, ¿2, f) W(Î1, ¿2, f) = W0(Î1, ¿2, f)

(14a) (14b) (14c)

where (u0, v0,w0) are the displacements of a point (fi, f2,0) on the mid-surface of the shell, and / and /2 are the rotations of the normal to the reference surface. It is shown in a number of works [45-47] that the linear deformations in orthogonal curvilinear coordinates of a general shell under these conditions are given by,

1 + f/Rr

-fx1 ю2

£2 + f£2 1 + f/R2

1 + f/R1 1 + f/R2

£4 " 1 + 4/R2 , where

£0 = 1 (du0

1 + f/R1

a1 V@Î1

£0 = 1 (dvi

2 аД

v0 da1 a1

TTÏÏir + 1T W0 a2 d^2 R1

u0 da1 a2

0 1 fdw0 , a2 N

£0 = Ж + a2/2 - R2 v ")

0 1 (dw0 a1 \

e° = ал%+a1/1- R1u0)

X°-l (dVl

X — aAÔf1

x1 — — 1a1

d/i4 9Î1

ÔÎ1 '

u0 @a1

~ä~2 Ж2

/2 da1 Ъ W2

/1 da1 " di2


(15c) (15d-e)

(16a) (16b) (16c) (16d)

- V0 da2) (16e-f )

fdu0 V0 dai

®2 — а Д d^2 a1 df1

£2 — — 2 a2

Ц+/1 da2) (16g-h)

a1 df1

Ю2 — — 2 a2

- z2 da^) (16i-j)

/2 da2 07 d^

and the quantaties a1 and a2 are the curvilinear scale factors that represent the length of the vectors tangent to curves of constant f1 and n2. In the present work, fibre angles may only vary in one direction and are always symmetric about the centre of the panel. In this case, only two possible cylindrical configurations, as shown

in Fig. 4, can occur. Note that the global co-ordinate system (X, Y,Z) is used to describe the fibre paths and loads applied on the structure. The local co-ordinate system (x,y,z) is used to describe the shell reference surface and deformation variables. The circumferential co-ordinatey is used in preference to an angular one (0) since it results in simpler equations and allows for easier correlations between a curved and flat panel [48]. Thus for cylindrical shells using the co-ordinate convention described above Eq. (16) reduce greatly in complexity since

¿1 = x, R =1, a1 = U2 = 0 = y/R, R2 = R, a2 = R

In this analysis the radius of curvature of the panel R is calculated from the thickness profile th of the 3D structure. Since this is a trigonometric function of the shearing angle Eq. (1) the radius of curvature may vary around the circumferential co-ordinate y as follows,

R(y)-1/j =

(1 + dte/dyf2 fît»/dy2

to create a generally elliptic reference surface. Returning to the kinematic Eq. (16) we see that for cylinders of elliptic cross-section the only term that explicitly captures the change in radius of curvature da2/df1 = dR(y)/dx must vanish. Thus, it is interesting to point out that for VAT panels with fibre variations in one direction the partial derivative of R in the prismatic direction does not feature in the governing equations. However, in the general case, with varying fibre angles in both X and Y directions, the shape of the VAT panel may take on more complex, doubly-curved shapes. Based on this premise, fibre path optimisation studies regarding maximum buckling loads may reveal entirely different and perhaps enhanced optimal solutions compared with previous flat plate analyses [5,6,8,9]. Of course singly- or doubly-curved shells exhibit more complex and often unstable post-buckling characteristics that have to be accounted for. A general optimisation study investigating these ideas will be the topic of future work.

Based on the observations made above, the linear kinematics of Eq. (16) can be greatly simplified. In addition, these equations are supplemented with the assumptions underlying Donnell's non-linear theory for shallow shells [49]

1. The shell is sufficiently shallow such that L < R and curvature changes are small.

2. Shell deformations u and vare much smaller than the thickness dimension h so that they are linear functions of w only. However, w can be of the same order as h.

3. Derivatives of w are small but their squares and products are of the same order as the linear strains.

Following these assumptions the effect of f/R in the shell kinematics is neglected and the non-linear von Kármán strains applicable for deformations dominated by the normal displacement w are introduced. Thus the final in-plane strain components are given by

Cy \6xy/

( 2( @t)2\

dWoY dy ) dwo dwo , dx dy /

I duo . 1 /dWtA2 \ dx + 2 V dx !

dvo Wo i l(9Wo\ dy + R + 2 dy J

duo dvo dwo dwo > 9y + dx + dx 9y

d/x ■

\ 9/L i 9Íí.

\ 9v + dx /

= Co + zr

Eq. (18) reveal the relative simplicity in which Donnell's equations incorporate the stretching-bending effect induced by the initial curvature of the cylindrical shell compared to the more involved Flügge deep-shell or Sanders equations. Compared to the kinematics of a flat plate only the w/R term in the expression for ey is introduced to capture the change in circumference with changing radius [48]. Finally, the transverse shear strains are given by,

fc) - < C4) - № ,)

where the shear correction factor k = 5/6 has been introduced to account for the parabolic variation of shear strain through the thickness. The governing equilibrium equations are derived using the principle of virtual work

Fig. 4. Two VAT panels with corresponding thickness variations. Note that the co-ordinate system (x,y,z) describing the shell reference surface is not always aligned with the global co-ordinate system (X, Y, Z) used to defining the fibre orientations and applied loads.

dV — J {axSex + Gyèey + aXy SeXy + Xxz èy^ + tyzáyyz)dV " / /(rnèUn + XnsèUs + T„zèw)dzds

where ffn, x„5, inz are the prescribed stresses on the boundary con-

tour. The stress resultants per unit width are defined as usual

(Nx, Ny, Nxy) = / (ax, ay, a^dz (21a)

-t/2 t/2 -t/2

(Mx, My, Mxy)= / z(rx, ay, rxy)dz -t/2

(Qx, Qy) = k (Xxz, Xyz)dz -t/2

The non-linear kinematics of Eqs. (18) and (19) are substituted into Eq. (20) and the calculus of variations applied on the unknown functional fields u0, v0, w0, /x and /y taking into account Eq. (21). Thus the resulting Euler equations that define equilibrium are defined by

dNx dN: dx +

-xy = 0,

xy i dNy i Qy _ o

dNy dy

@ Mx @ M

@x @Qx

@y @Qy

- Qx = o,

+ - N + Al N

dx + dy R + dx[ '

dx dwo

+^ - Q- o

(22a-b) (22c-d)

d (n. @wi + m, @wO)

The Euler Eq. (22) above represent a system of five non-linear partial differential equations for the coupled stretching-bending problem of a cylindrical shell that can be solved for the five unknown fields u, v, w, /x and /y when the boundary conditions in Eq. (23) are specified

u — u or Nxnx + Nxyny — Nxnx + Nxyny v — v or Nxynx + Nyny — Nxynx + Nyny

/x — /x or Mxnx + Mxyny — Mxnx + Mxyny /y — /y or Mxynx + Myny — Mxynx + Myny

w — WW or (Nxnx + Nxyny) — + (Nxynx + Nyny) — + Qxnx + Qyny — (Qx nx + Qyny

(23a) (23b) (23c) (23d)

The set of non-linear differential Eq. (22) are reduced to an infinite set of linear differential equations by expanding the unknown functional fields in a power series of an arbitrary perturbation parameter a about the point of buckling. Since we are dealing with a coupled stretching-bending problem all unknown functional fields have to be expanded in terms of both even and odd powers of the arbitrary perturbation parameter. Likewise it is required to expand the internal stress resultants in terms of the arbitrary parameter. Thus,

~ = (U0 v 0 W0 /x /y)T

F — (Nx Ny

Nxy Mx

My Mxy Qx Qy)'

F,F)= F0,F0


F2,F2 a

= ]T( F(i',F(9) a

where the variables multiplied with the a0 and a1 are the displacements and stress resultants corresponding to the pre-buckling and buckling solutions respectively. The higher order variables are those related to the post-buckling solutions.

For a shell of N anisotropic layers the constitutive law for each layer k is given by

where the in-plane transformed stiffnesses Qij i, j = 1, 2, 6 are defined as in CLA [17]. The transverse stiffness terms for a layer k with angle hk to the principle material coordinates are given by

(ax) (k) Q11 Q12 o o Q16" (k) f€x)

ay Q12 Q22 o o Q 26 £y

xxz = o o Q44 Q45 o yxz

xyz o o Q45 Q 55 o yyz

axy Q16 Q26 o o Q 66. VxyJ

Q44 — Q44 cos2 9k + Q55 sin2 h

Q45 — (Q55 - Q44) cos 8k sin ek

M55 — Q55 cos2 hk + Q44 sin2 h with Q44 — G23 and Q55 — G1:

(26c) (26d)

It is clear from Eq. (26c) that ifthe condition oftransverse orthotro-py is assumed such that G13 = G23, the coupling term Q45 vanishes. Integrating the stresses of Eq. (25) through the thickness in regard of Eqs. (18), (19) and (21) we get,

A44 A45

A45 A55

( @w0 + /x \ V@wO+/y - Vt)

where the laminate stiffness terms A(X, Y), B(X, Y) and D(X, Y) are defined as in CLA but may vary depending upon the location (X, Y) on the shell reference surface. The analysis in this work is restricted to symmetric VAT laminates such that the coupling matrix B vanishes and is disregarded hereon. Writing the stress resultants of Eq. (27) using the perturbed expansion of the displacements given in Eq. (24) the stress resultants in terms of the perturbation parameter a are obtained. In the following, only the first two terms of the infinite series are explicitly shown since in this study we are primarily interested in the pre-buckling and buckling solutions. Thus,

Ai6 ¿26 ¿66.

híw0 dwL , dwo dwl

r dy dy

'A+dwo ü +dw0

dx + dx dy + dx dy

a1 + a2 [..

D11 D12 d16 '

d12 D22 D26 d16 d26 d66

a1 + a2 [..

d/x I "yy

V dy + dx /

0 + /0 V A* + /1.,1 |a1 + a2[..

0 + /J-

where we have assumed that the partial derivatives of dw00=dx and dw00=dy are negligible in the pre-buckling problem. Since the pertur-

duo . dvo

d/x , d/

V dy + dx !

bation parameter a is arbitrary the power series can only vanish if the coefficient of each parameter a(n) is set equal to zero. Substituting the perturbed stress resultant expressions of Eq. (28) into the governing partial differential Eq. (22) and enforcing the requirement of vanishing coefficients we obtain an infinite set of linear equations,

@N0 dN0i = 0 "

dx + dy dx + dy + R

- Q0 = 0

9M0 -My Sx + -y

-Ma i £M° _ Q 0—0

dx + 8y Qy = 0

@q0 @q0_ n0 = 0

ax + 9y R 0

9x + 9l u

-»i 9n1 q1 — 0 9x + 9l + R

di+d--■ - q1 = 0 dJw+îs1 - q1 = 0

9q1 «j ¿(n

9x. + 9l r + dx. \ :

where the first set of Eq. (29) is solved to obtain the pre-buckling solution (u0, v0, w0, /0, /0 ) for a given set of boundary conditions,

whereupon the second set of Eq. (30) is solved to obtain the buckling eigenvalues. Buckling is assumed to be an instantaneous event such that there can be no change in load or displacement on the boundary. Thus non-zero boundary conditions of the pre-buckling problem are now set to zero in the buckling problem.

2.5. Differential quadrature method

Recently, the Differential Quadrature Method (DQM) has been shown to be a fast, accurate and computationally efficient technique for solving the variable coefficient, higher order differential equations for bending [18], buckling [7] and post-buckling [50] of VAT plates. In the latter study the accuracy of the DQM approach was validated for a number of free, simply-supported and clamped plate boundary conditions. In comprehensive surveys by Viola, Tornabene and Fantuzzi the DQM was effectively applied to analyse the free vibration [51,52] and static behaviour [53] of doubly-curved, straight fibre laminates based on a large number of higher-order shear deformation theories.

Differential quadrature is a numerical discretization technique that approximates the partial derivatives of a functional field with respect to a specific spatial variable using a linear weighted sum of all the functional values in the domain. For example, the nth partial derivative at the ith discretization point is,

f i= 1,2,

where xi is the set of discretization points in the x-direction, typically defined by the non-uniform Gauss-Labotto-Chebychev distribution,

1 - cos( NTT p

i = 1,2,..., N

where An are the weighting coefficients of the nth derivative and repeated index j means summation from 1 to N. The weighting coefficients are determined by approximating f(x) using test functions. In this work, Lagrange polynomials defined by Quan and Chang [54] are used,

gk(x) =

(x - xk)M(1)(xk)'

k = 1,2,..., N

where M(x) = ]J(x - j M(1)(x,) = fl (x - xk) j=1 1N


) A(1) =

xi - xk

xi xik=1,k-ijxi xk N1

and A« = V —

" y.- — :

for i — j

>xi - xk

The higher order weighting coefficients are then obtained by direct

matrix multiplication: [A(m)] = [A(1)][A(m-1)], m = 2, 3.....N - 1 [55].

Thus, a differential equation is reduced to a system of algebraic equations, with each equation corresponding to a gridpoint in the domain, such that the unknown functional values at each gridpoint are found by solving the system of equations with the appropriate boundary conditions. Using this technique the in-plane problem for the non-symmetric flat plate of Eq. (13) is written as

I<1u K1v

K2u K2v

KU = 0

and once u and v are known across the plate, the general eigenvalue problem defined by Eq. (9) is solved,

K11 K12 K13 K14 0 0 0 0 1 (Mx\

K21 K22 K23 K24 - x 0 0 0 0 1 My

K31 K32 K33 K34 0 0 0 0 1 Mxy

.K41 K42 K43 K44 0 0 0 N44. 1 \ w /

= (K - xN)M = 0

where k is the buckling eigenvalue, K is a matrix of DQ coefficients, bending flexibilities and shear coefficients, X the vector of unknown functional values at the grid points, and N a matrix of DQ coefficients multiplied by pre-buckling in-plane stress resultants. Of course, the same methodology can be applied for solving the shell buckling problem. In this case five differential equations are discre-tised for both the pre-buckling and buckling problems since the inplane and out-of-plane behaviours are coupled.

The boundary equations can be discretised in a similar fashion and then applied using Shu and Du's general approach [56]. Here the discretised mesh is separated into functional unknowns on the boundary and functional unknowns vi in the interior. By discretising the field equations at the interior points and the boundary conditions at the boundary points we get,

Field :[K>b]- Ub + [K>>]- X> = k[N>>]X> Boundary :[Kbb] Ub + [Kb>]-lJ> = 0

) {[Kii]-[Kib][Kbb]-1 [Kw]- k[N>>]}X> = 0 )(K - kN) X> = 0

and thus the buckling eigenvalue and mode shape is found. It is important to point out that the DQ stiffness matrices are densely populated, and due to the variable coefficient nature of the differential equations, different rows of the matrices [K>>] and [Kbb] can have largely different magnitudes. To condition these matrices each row is multiplied by a unique normalisation factor. This normalisation vector is calculated from the arithmetic mean of the absolute values of all elements in the matrix divided by the maximum absolute value of the elements along each row.

3. Results and model validation

Both the flat plate and curved panel models were used to simulate the buckling behaviour of a square VAT laminate with unit in-plane dimensions (lx = ly = 1) under four different boundary conditions A1, A2, B1 and B2 illustrated in Fig. 5a and b. For case A the

n1 nlss1+ ni,9w°+ n0,9s

transverse edges Y = 0 and Y = 1 are free to move in-plane while for case B they are constrained. Note, that since the fibre orientations may vary along the edges the boundary can deform if it is free to do so (see Fig. 5a). For both cases A and B all four edges of the panel were either pinned (A1 and B1) or clamped (A2 and B2). Finally, in all four cases a uniform end-shortening of DX = 1 mm was applied. As is apparent from Fig. 4, in the local shell co-ordinate system this corresponds to an applied loading in terms of v0 and u0 for fibre variations in the X- and Y-directions respectively.

Orthotropic carbon-epoxy material properties were defined as E11 = 150 GPa, E22 = 8.8 GPa, G12 = 4.8 GPa, G13 = G23 = 4 GPa and v12 = 0.35. Simulations were run on several symmetric 4-ply VAT laminates with nominal span to thickness ratios (L/t) of 100:1 and 50:1. For a 4-ply [90 ± (0|70)]s these aspect ratios would reduce by a factor of 3 in the areas of maximum thickness buildup, thereby accentuating transverse shear effects.

The numerical results were verified using the FEM package ABAQUS by creating both the non-symmetric plate and symmetric shell models using quadratic 2D S8R FSDT elements. In order to represent the VAT fibre and thickness distributions a unique straight fibre composite layup was defined for each S8R element on the reference surface. However, for the non-symmetric plate the element offset was set relative to Curve 2 shown in Fig. 1 in order to model the non-symmetric thickness variation. Furthermore a 3D model of the real structure, meshed using linear C3D8R brick elements with enhanced hourglassing control, was created to serve as the benchmark for all results. For this 3D model the variable thickness outline of the VAT panel was created as a solid part and then uniformly meshed using a single brick element per ply. The variable fibre paths were discretised by assigning rotated material orientations for each element depending on its location in the XYZ-space. Mesh convergence studies were performed for each configuration to find the number of DQ grid points/FE elements (N) required in the X- and Y-directions for the DQM/FE analysis to obtain converged results. Subsequently, meshsizes of Ndq =20 x 20, N2D = 50 x 50 and N3D = 4 x 40 x 40 were used throughout the analyses.

Tables 1 and 2 compare the numerical results obtained by the two analytical models (DQ Plate and DQ Shell) with the three FEM solutions for a [0 ± <4510)]s and a [90 ± (0|45)]s laminate. The buckling mode shape of a [90 ± <0|45)]s laminate under Case A1 boundary conditions obtained by the DQ Plate and DQ Shell models are compared to the 3D FEM solution in Fig. 6. The mode shapes can be seen to be very similar and had same accuracy for all configurations tested.

Tables 1 and 2 show that the results of the DQ plate model correlate closely with the 2D FE plate analysis for all configurations. Since the 2D FE code in ABAQUS directly solves the coupled stretching-bending problem of the non-symmetric plate, this correlation suggests that the reduced stiffness approach employed in the DQ plate model is accurate for these laminates. The DQ shell and 2D FE shell results are similarly well correlated with good overall accuracy compared to the 3D FE benchmark. Since pre-buckling transverse shear effects were ignored in the DQ plate theory but are incorporated in the 2D FE plate model, it can be concluded that transverse shear effects, on prebuckling, are largely negligible for linear fibre variations with local length to thickness ratios of up to L/t = 35:1.

The results of the DQ plate and 2D FE plate models generally correlate well with the 3D FEM results. Most results are within 5% of the full 3D FE solution and for a number of configurations the results are within 1%. On the other hand, Table 1 clearly highlights a striking inaccuracy for both laminates under boundary condition B1. In this case the DQ plate and 2D FE buckling eigenvalues are actually lower than the 3D FEM benchmark. One would actually expect the results of the 3D model to be lower since the effects of transverse normal stresses are taken into account to give more degrees of freedom. One explanation for this is because the flat plate model disregards the added stiffness effects of the initial curvature. Apart from two exceptions the 2D FE shell results are always closer to the 3D FE benchmark than the 2D FE plate results. Even though this trend is not as pronounced for DQ plate compared to DQ shell results, the general trend is for the shell results to be more closely matched to the 3D FE benchmark.

Tables 1 and 2 also shows that this trend varies with the laminate configuration and type of boundary condition. Indeed, Fig. 7a shows that both the DQ shell and 2D FE shell models are increasingly more accurate as the fibre angle T2, and therefore degree of initial curvature, is increased for a pinned [90 ± <0|T2)]s laminate. Here the variations in fibre angle have been constrained to 70° since this represents the manufacturing limit of the CTS technique. At this maximum steering angle of T2 = 70° the predictions of both the 2D FE and DQ plate models are conservative by a considerable margin of 26%. A similar observation is made for increasing angle T1 for a [0 ± <T1|0)]s laminate in Fig. 7b, where the discrepancy is 22%. Comparing this to Fig. 7c and d clamping the edges of the panel reduces the importance of initial curvature on the buckling results such that the 2D FE plate model remains accurate even for high steering angles. In Fig. 7b and d it can be seen that for high steering angles the DQ plate solution diverges from the 2D FE plate

Al: All Edges Pinned All All Edges Clamped

(a) Case A: Transverse edges free to move in-plane



Bl: All Edges Pinned B2: All Edges Clamped

(b) Case B: Transverse edges constrained

Fig. 5. VAT plate illustrating the four different boundary conditions investigated.

Table 1

Comparison of DQM Plate, 2D FEM Plate, DQM Shell, 2D FEM Shell and 3D FEM buckling eigenvalues (mm) for various VAT laminates with pinned edges and boundary conditions A1 and B1.


L/t (nominal,thickest)

DQ Plate

FEM S8R Plate

DQ Shell

FEM S8R Shell


Al: All edges pinned. Transverse edges free

[0 ± <45|0)]s 100,70 0.186

50,35 0.731

[90 ± <0|45)]s 100,70 1.044

50, 35 4.031

Bl: All edges pinned. Transverse edges constrained

[0 ± <45|0)]s 100,70 0.090

50, 35 0.355

[90 ± <0|45)]s 100,70 0.370

50, 35 1.443

0.180 0.710 1.050 4.100

0.090 0.354 0.359 1.409

0.192 0.705 1.042 4.076

0.092 0.334 0.381 1.428

0.180 0.712 1.043 4.090

0.091 0.340 0.388 1.526

0.185 0.747 1.018 3.863

0.092 0.391 0.406 1.744

Table 2

Comparison of DQM Plate, 2D FEM Plate, DQM Shell, 2D FEM Shell and 3D FEM buckling eigenvalues (mm) for various VAT laminates with clamped edges and boundary conditions A2 and B2.


L/t (nominal,thickest)

DQ Plate

FEM S8R Plate

DQ Shell

FEM S8R Shell


A2: All edges clamped. Transverse edges free

[0 ± <45|0)]s 100,70 0.506

50, 35 1.883

[90 ± <0|45)]s 100,70 1.483

50, 35 5.560

B2: All edges clamped. Transverse edges constrained

[0 ± <45|0)]s 100,70 0.256

50, 35 0.969

[90 ± <0|45)]s 100,70 0.785

50, 35 2.951

0.482 1.834 1.492 5.676

0.237 0.906 0.743 2.838

0.481 1.8 1.488 5.679

0.239 0.893 0.735 2.804

0.482 1.839 1.492 5.677

0.238 0.909 0.743 2.846

0.497 1.873 1.534 5.663

0.238 0.915 0.742 2.835

(a) DQM Plate (N=30x30) (b) DQM Shell (N=30x30) (c) C3D20R FEM (N=4x40x40)

Fig. 6. Buckling mode shape of a [90 ± <0|45)]s laminate under Case A1.

result. We believe this is due to the fact that as the bending-stretching coupling effects increase with steering angle, the approximate solution based on the reduced stiffness matrices becomes increasingly inaccurate.

In Fig. 7c it is also observed that the DQ plate results are inaccurate beyond the T2 = 30° datapoint and do not follow the general trend of the other curves. We believe this happens due to numerical ill-conditioning. In the governing differential Eqs. 9a, 9b and 9c some of the moment terms are multiplied by the flexibility terms dij which can be up to 15 orders of magnitude lower than the DQ weighting coefficients of the geometric curvature terms (e.g. 2w/ dx2). Thus elements along the rows of the stiffness matrix K, where each row corresponds to the governing equation at a certain grid-point, can be of largely different magnitudes. Under these conditions the matrix can have a non-zero determinant but the 2-

norm condition number may be large leading to an ill-conditioned matrix. Thus, the use of dij flexibility terms is an inherent drawback of the chosen theoretical framework in terms of solving the problem numerically using the DQ method. This numerical instability is not encountered to the same extent in the FSDT shell model since the A and D matrices are of similar orders of magnitude as the geometric curvature terms. However, the assembled DQ stiffness matrix K is naturally dense, non-symmetric and relatively large due to solving a system with 4 or 5 unknown functional fields (1600 x 1600 or 2000 x 2000 for N = 20 x 20). It was observed that for some boundary conditions and fibre variations the condition number of K for the buckling problem increases with steering angle as shown for two cases in Table 3. Furthermore, even for the straight fibre laminate the condition number is indeed very high. In the equivalent cases of solving the equations using CLA (i.e. only

20 30 40 50

Angle T2 (deg)

Case Bl: Xcr vs. T2 for a [90±<0|T2>] laminate

0 10 20 30 40 50 60 70

Angle T1 (deg)

(b) Case Bl: Xcr vs. Ti for a [0±<Ti|0>]s laminate

3D FEM DQ Plate -0-2D FE Plate DQ Shell 2D FE Shell

0 10 20 30 40 50 60 70

Angle T2 (deg)

(c) Case B2: Xcr vs. T2 for a [90±<0|T2>]S laminate

3D FEM DQ Plate -$-2D FE Plate DQ Shell 2D FE Shell

20 30 40 50

Angle T1 (deg)

(d) Case B2: ACT- vs. Ti for a [0±<Ti|0>] laminate

Fig. 7. Comparison of DQ Plate, DQ Shell, 2D FE Plate, 2D FE Shell and 3D FEM buckling eigenvalue kcr vs. versus fibre angle T1 or T2 for either clamped or pinned edges with transverse edges constrained.

Table 3

Condition number of the DQ stiffness matrix K when solving the buckling problem with increasing steering angle.

[0 ± (1,10) ]s DQ Shell under Case A2

Ti 0 15 30 45 60 70

Cond (K) 1.26E+16 1.70E+17 4.64E+16 2.59E+16 2.83E+16 8.74E+16

[90 ± (01T2)]s DQ Plate under Case B2

T2 0 15 30 45 60 70

Cond (K) 4.68E+08 5.15E+08 9.45E+08 3.19E+09 9.71E+09 1.44E+11

three variables u, v, w) produces 1e3 < cond (K) < 1e4. Thus applying the DQM to solve the buckling problem of variable angle, variable thickness laminates using enhanced shear models may cause large, dense matrices that introduce numerical inaccuracies. To improve the accuracy and stability in the future, pre-conditioning techniques and DQ schemes that produce box-banded or uniformly banded weighting matrices will be tested.

Finally, it is observed in Fig. 7b and c that the DQ shell results are slightly lower than the 2D FE shell results for values of T1 and T2 greater than 60°. One explanation for this could be ill-conditioning effects as discussed above, since for T1, T2 > 60° both the DQ plate and shell results are seen to diverge. However, for the panel dimensions assumed here high shearing angles create local radius of curvature to length ratios (R/L) close to 1.5:1. Under these circumstances the assumption of a shallow shell is unjustified

and a non-shallow shell representation such as the Sanders-Koiter non-linear theory [57] would be expected to be more accurate.

4. Conclusions

The Continuous Tow Shearing technique recently developed to improve the manufacture of laminates with variable angle tows, produces laminates with a flat profile on one side and a smooth curved profile on the other. The buckling behaviour of such laminates with prismatic fibre variations and symmetric stacking sequences was investigated based on the premise that the three-dimensional profile can either be modelled as a flat, non-symmetric plate or a cylindrical, symmetric panel. In both approaches the effects of transverse shear deformation were taken into account. The governing differential equations were solved in the strong

form using the Differential Quadrature method and shown to be in good agreement with 2D finite element models based on the same assumptions. It was shown that the analytical and finite element models based on the cylindrical shell formulation correlate closest when compared to the 3D benchmark. Thus it is suggested that the buckling event of these variable angle tow, variable thickness laminates is characterised more accurately by "shell-like" than by "plate-like" behaviour.

However, further analytical and experimental evidence in terms of the post-buckling path of the 3D structure is needed in order to clarify this concept in detail. Unstable load drops in the post-buckling regime would point towards "shell-like" behaviour and validate the findings of the present work. Nonetheless, the CTS technique creates the possibility of designing smooth, curved panels which add another dimension to the optimisation problem of laminates with variable angle tows. Finally, by extending the analysis to laminates with non-prismatic fibre variations the idea of increasing the buckling loads by inducing more complex, doubly curved profiles is proposed.


The authors would also like to sincerely thank the EPSRC for funding the ACCIS Doctoral Training Centre (EP/G036772/1) at Bristol.


[1] Hyer MW, Lee HH. The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes. Compos Struct 1991;18:239-61.

[2] Hyer MW, Charette RF. The use of curvilinear fiber format in composite structure design. AIAAJ 1991;29(6):1011-5.

[3] Gurdal Z, Olmedo R. In-plane response of laminates with spatially varying fiber orientations: variable stiffness concept. AIAAJ 1993;31(4):751-8.

[4] Gurdal Z, Tatting BF, Wu CK. Variable stiffness composite panels: effects of stiffness variation on the in-plane and buckling response. Compos: Part A 2008;39:911-22.

[5] Setoodeh S, Abdalla MM, Ijsselmuiden ST, Gurdal Z. Design of variable stiffness composite panels for maximum buckling load. Compos Struct 2008;87:109-17.

[6] Wu Z, Weaver PM, Gangadharan R, Kim BC. Buckling analysis and optimisation of variable angle tow composite plates. Thin-Walled Struct 2012;60:163-72.

[7] Raju G, Wu Z, Kim BC, Weaver PM. Prebuckling and buckling analysis of variable angle tow plates with general boundary conditions. Compos Struct 2012;94(9):2961-70.

[8] van den Brink WM, Vankan WJ, Maas, R. Buckling optimized variable stiffness laminates for a composite fuselage window section. In: Proceedings of the 28th international congress of the aeronautical sciences, Brisbane (Australia); 2012.

[9] Alhajahmad A, Abdalla MM, Gurdal Z. Optimal design of tow-placed fuselage panels for maximum strength with buckling considerations. J Aircraft 2010;47(3):775-82.

[10] Wu Z, Gangadharan R, Weaver PM. Postbuckling analysis of variable angle tow composite plates. Int J Solids Struct 2013;50:1770-80.

[11] Wu Z, Weaver PM, Gangadharan R. Postbuckling optimisation of variable angle tow composite plates. Compos Struct 2013;103:34-42.

[12] Blom AW. Structural performance of fiberplaced variable-stiffness composite conical and cylindrical shells. 2nd revised ed., PhD Thesis. University of Delft, Holland; 2010.

[13] Beakou A, Cano M, Le Cam J-B, Verney V. Modelling slit tape buckling during automated prepreg manufacturing: A local approach. Compos Struct 2011;93:2628-35.

[14] Kim BC, Potter K, Weaver PM. Continuous tow shearing for manufacturing variable angle tow composites. Compos: Part A 2012;43(8):1347-56.

[15] Carrera E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. Arch Comput Meth Eng 2002;9(2):87-140.

[16] Reddy JN. Energy and variational methods in applied mechanics. John Wiley & Sons; 1984. p. 354-88.

[17] Jones RM. Mechanics of composite materials. 2nd revised ed. London (UK): Taylor & Francis Ltd.; 1998.

[18] Groh RMJ, Weaver PM, White S, Raju G, Wu Z. A 2D equivalent single-layer formulation for the effect of transverse shear on laminated plates with curvilinear fibres. Compos Struct 2013;100:464-78.

[19] Von Kärmän T, Sechler EE, Donnell LH. The strength of thin plates in compression. ASME Trans 1932;54(2):53-7.

[20] Marguerre K, Trefftz E. Über die Tragfähigkeit eines längsbelasteten Plattenstreifens anch Überschreiten der Beullast. ZfaMM 1937;17(2):85-100.

[21] Kromm A, Marguerre K. Behaviour of a plate-strip under shear and compressive stresses beyond the buckling limit. NACA TM 870; 1938.

[22] Koiter WT. De meedragende breedte bij groote overschrijding der knikspanning voor verschillende inklemming der plaatranden. Rep. S. 287. Nationaal lucht- en ruimtevaartlaboratorium, Amsterdam; 1943.

[23] Levy S. Bending of rectangular plates with large deflections. NACA Rep. 737, 1942.

[24] Hu PC, Lundquist EE, Batdorf SB. Effect of small deviations from flatness on effective width and buckling of plates in compression. NACA TN 1124,1946.

[25] Stein M. Loads and deformations of buckled rectangular plates. NASA TR R-40, 1959.

[26] Chandra R, Raju BB. Postbuckling analysis of rectangular orthotropic plates. Int JMech Sci 1973;15:81-97.

[27] Ashton JE. Approximate solutions for unsymmetrically laminated plates. J Compos Mater 1969;3:189.

[28] Mindlin RD. Influence of rotary inertia and shear on flexural motion of isotropic elastic plates. ASME J Appl Mech 1951;18:31-8.

[29] Cosentino E, Weaver PM. An enhanced single-layer variational formulation for the effect of transverse shear on laminated orthotropic plates. Eur J Mech A/ Solids 2010;29:567-90.

[30] Ambartsumyan SA. On theory of bending plates. Isz Otd Tech Nauk AN SSSR 1958;5:6977.

[31] Reissner E. On transverse bending of plates, including the effect of transverse shear deformation. Int J Solids Struct 1975;11:56973.

[32] Reddy JN. A refined shear deformation theory for the analysis of laminated plates. NASA contractor report 3955; 1986.

[33] Levy M. Memoire sur la theorie des plaques elastique planes. J Math Pures Appl 1877;30:219306.

[34] Stein M. Nonlinear theory for plates and shells including the effect of transverse sharing. AIAAJ 1986;24:153744.

[35] Touratier M. An efficient standard plate theory. Int J Eng Sci 1991;29:90116.

[36] Karama M, Abou Harb B, Mistou S, Caperaa S. Bending, buckling and free vibration of laminated composite with a transverse shear stress continuity model. Compos Part B Eng 1998;29B:22334.

[37] Ferreira AJM, Roque CMC, Jorge RMN. Analysis of composite plates by trigonometric shear deformation theory and multiquadrics. Comput Struct 2005;83:222537.

[38] Soldatos KP. A transverse shear deformation theory for homogeneous monoclinic plates. Acta Mech 1992;94:195220.

[39] Neves AMA, Ferreira AJM, Carrera E, Cinefra M, Roque CMC, Jorge RMN, et al. Free vibration analysis of functionally graded shells by a higher-order shear deformation theory and radial basis functions collocation, accounting for through-the-thickness deformations. Eur J Mech A Solid 2013;37:2434.

[40] Karama M, Afaq KS, Mistou S. Mechanical behaviour of laminated composite beam by the new multi-layered laminated composite structures model with transverse shear stress continuity. Int J Solids Struct 2003;40:152546.

[41] Mantari JL, Oktem AS, Guedes Soares C. Static and dynamic analysis of laminated composite and sandwich plates and shells by using a new higherorder shear deformation theory. Compos Struct 2011;94:3749.

[42] Murakami H. Laminated composite plate theory with improved in-plane responses. ASME J Appl Mech 1986;53:661666.

[43] Reissner E. The effect of transverse shear deformation on the bending of elastic plates. J Appl Mech 1945;12:A69-77.

[44] Koiter WT. A consistent first approximation in the general theory of thin elastic shells. In: Proceedings of the first symposium on the theory of thin elastic shells, Amsterdam; 1960. p. 12-23.

[45] Reddy JN. Mechanics of laminated composite plates and shells. 2nd ed. Boca Raton (Florida, USA): CRC Press LLC; 2004.

[46] Reddy JN. Exact solutions of moderately thick laminated shells. J Eng Mech 1984;110(5):794809.

[47] Pirrera A, Weaver PM. Geometrically nonlinear first-order shear deformation theory for general anisotropic shells. AIAAJ 2009;47(3):767-82.

[48] Donnell LH. Stability of thin-walled tubes under torsion, vol. 479, National Advisory Committee for Aeronautics; 1935.

[49] Yamaki N. Elastic stability of circular cylindrical shells. The Netherlands: Elsevier Science Publishers B.V. Amsterdam; 1984.

[50] Raju G, Wu Z, Weaver PM. Postbuckling analysis of variable angle tow plates using differential quadrature method. Compos Struct 2013;106:74-84.

[51] Viola E, Tornabene F, Fantuzzi N. General higher-order shear deformation theories for the free vibration analysis of completely doubly-curved laminated shells and panels. Compos Struct 2013;95:639666.

[52] Tornabene F, Viola E, Fantuzzi N. General higher-order equivalent single layer theory for free vibrations of doubly-curved laminated composite shells and panels. Compos Struct 2013;104:94117.

[53] Viola E, Tornabene F, Fantuzzi N. Static analysis of completely doubly-curved laminated shells and panels using general higher-order shear deformation theories. Compos Struct 2013;101:59-93.

[54] Quan JR, Chang CT. New insights in solving distributed system equations by the quadrature method. 1. Analysis. Comput Chem Eng 1989;13(7):779-88.

[55] Shu C. Differential quadrature and its application in engineering. Springer Verlag; 2000. p. 37,208.

[56] Shu C, Du H. Implementation of clamped and simply supported boundary conditions in the GDQ free vibration analysis of beams and plates. Int J Solids Struct 1997;23:819-35.

[57] Sanders Jr JL. Nonlinear theories for thin shells. QAppl Mech 1963;21:21-36.