ELSEVIER

Contents lists available at ScienceDirect

Composite Structures

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

Buckling analysis of stiffened variable angle tow panels

Broderick H. Coburn *, Zhangming Wu, Paul M. Weaver

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

ARTICLE INFO ABSTRACT

Variable angle tow (VAT) laminates have previously shown enhanced buckling performance compared to conventional straight fibre laminates. In this study, an analytical method is developed for the buckling analysis of a novel blade stiffened VAT panel to allow this potential to be more fully exploited. The pre-buckling and buckling analysis, performed on a representative section of a blade stiffened VAT panel, are based on a generalised Rayleigh-Ritz procedure. The buckling analysis includes a first order shear deformation theory by introducing additional shape functions for transverse shear and is therefore applicable to structures with thick skins relative to characteristic length. Modelling of the stiffener is achieved with two approaches; idealisation as a beam attached to the skin's midplane and as a rigidly attached plate. Comparing results with finite element analysis (Abaqus) for selected case studies, local buckling errors for the beam model and plate model were found to be less than 3% and 2% respectively, whilst the beam model error for global buckling was between 3% and 10%. The analytical model provides an accurate alternative to the computationally expensive finite element analysis and is therefore suitable for future work on the design and optimisation of stiffened VAT panels.

© 2013 Elsevier Ltd. All rights reserved.

CrossMark

Article history:

Available online 31 December 2013

Keywords: Buckling Rayleigh-Ritz Stiffened panel Transverse shear Variable angle tow

1. Introduction

In recent decades the ability to tailor the stiffness and strength of composite structures has seen their increased use in aerospace applications [1]. Traditional tailoring is achieved by treating the fibre orientation of each ply as a variable and optimising the stacking sequence for laminate performance. Recent advancements of automated tape/fibre laying technologies has led to the possibility of having variable angle tow (VAT) laminates where the fibre orientation can change over the plane of a ply. This results in laminates with varying in-plane and out-of-plane stiffnesses in the xy-plane providing designers with additional degrees of freedom and tailorability. Previous research on VAT plates has shown significant improvement in the stress distribution around holes [2-4] and buckling and post-buckling performance [5-7]. Increases in buckling load are primarily attributed to a redistribution of load to boundaries where the structure is constrained in out-of-plane displacement.

The majority of work to date on the design and optimisation of VAT laminates has utilised finite element analysis (FEA). Although relatively accurate, the fine meshes required to capture the mode shapes makes FEA of VAT laminates computationally expensive [8,9].

Recently, the differential quadrature method [10] and Rayleigh-Ritz energy method [7] have been shown to be viable alternatives

* Corresponding author. Tel.: +44 7869681735. E-mail address: broderick.coburn@bristol.ac.uk (B.H. Coburn).

0263-8223/$ - see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruct.2013.12.029

to FEA which are accurate, robust and computationally efficient and hence suitable for optimisation studies.

Although there has been a significant amount of research in the area of VAT laminates, in most studies the application is limited to simple geometries, i.e. plates and shells, with general boundary conditions. A large gap exists between the current understanding and analysis techniques of VAT and that required to apply a very promising technology to practical structural configurations.

A potential application, exploiting the enhanced buckling performance, is to use a VAT laminate as the skin of a stiffened panel. Here, the VAT skin would redistribute in-plane loads to the stiffen-ers which are then required to act as 'panel breakers' forcing a nodal line. Expected gains are an increased buckling performance allowing the design of lighter structures.

Stiffened panels are commonly used on aircraft as primary structures such as wing covers and fuselage panels [11]. Stiffened panels typically consist of a plate braced by longitudinal stiffeners and are an efficient configuration for carrying compressive loads, particularly when buckling is a design driver as is the case for aircraft wing covers [11]. A stiffened panel can fail via a variety of mechanisms including skin-stiffener debonding [12], material strength failure and buckling. Buckling failure predominantly occurs in one of two modes as shown in Fig. 1; local, where the stiffeners act as 'panel breakers' forcing the skin to buckle locally between the stiffeners and global, where both plate and stiffeners buckle out-of-plane. Confining the buckling mode to be local is preferential to global as it, in general, leads to lighter designs and greater post-buckling stiffness. The local mode's higher

Fig. 1. Stiffened panel local and global buckling mode shapes (yz-plane). (a) Full panel. (b). Representative section.

post-buckling stiffness is due to the unbuckled stiffeners carrying load in the post-buckling regime.

Buckling of stiffened panels has received considerable attention dating as far back as 1921 by Timoshenko [13] who used the Ritz-method to analyse isotropic longitudinally and transversely stiffened plates subject to compression, shear and bending. Continued interest in this field has seen many publications in the past century with current research focussing heavily on composite stiffened panels [14-30]. Local buckling analysis methods can be split into three categories based on the consideration of the stiffener. The first method treats the stiffener as a simple support which facilitates fast closed-form solutions to be obtained but assumes null torsional restraint and hence underestimates the buckling load [14]. The second method models the torsional restraint by replacing the stiffener blade with an equivalent torsional spring or beam attached to the skin's midplane [15-20]. This method is often sufficiently simple to obtain accurate closed-form solutions, however is strictly only valid for an unloaded stiffener and assumes no stiff-ener blade buckling or warping. Correction factors reducing the effective restraint in the case of an axially applied load to the stiff-ener have been proposed and provide improved solutions when load is carried by the stiffener [16,17]. The third method models both the skin and stiffener as plates [21,23,24] allowing local buckling modes of the stiffener and the interaction between the skin and stiffener to be captured. This higher fidelity approach has an increased computational cost but provides a more robust solution than the elastic restraint method.

The analysis of global buckling modes can be approached by either replacing the stiffeners with an equivalent smeared layer or by treating the stiffeners as discrete elements. The smeared approach [25,27] replaces the stiffeners with a distributed A-, B- and D-matrix over the entire panel. This approach is only valid for wide panels consisting of several closely spaced stiffeners. Alternatively, treating the stiffeners as discrete elements enables local interaction effects between the skin and stiffener to be captured and has no restrictions on stiffener spacing [23,25-28]. This is commonly achieved by replacing the stiffener with a beam element attached to the skin's midplane.

In all approaches for both local and global buckling the Ray-leigh-Ritz method is used extensively due to the ease of including the stiffener into the formulation and applying boundary conditions [13,19-26].

Several analysis and optimisation packages have been proposed to solve the linear buckling problem for stiffened panels. PANDA2 [31] uses simple models for the prebuckling, buckling and post-buckling of composite stiffened panels to obtain optimal solutions under a variety of loading conditions. Local and general buckling loads are calculated with the use of either closed-form expressions or discretised models of panel cross sections. The VICONOPT program [32] uses an exact finite strip theory for prismatic strips

where the boundary conditions at the ends are enforced using either polynomial or trigonometric shape functions.

Most research on the buckling of stiffened panels idealises both the skin and stiffener as thin plates neglecting transverse shear effects. However, stiffened panels in aerospace applications often have thick skins, relative to characteristic length, where transverse shear effects must be considered. For isotropic materials, a ratio of thickness to characteristic length of 1/10 results in approximately a 5% error by neglecting transverse shear [33]. For composite materials the ratio of in-plane to transverse shear moduli can be a factor of ten or more than isotropic materials and deformations due to transverse shear become important at even lower ratios of thickness to characteristic length [33,34].

Consideration of shear deformation was first proposed by Timoshenko [35] for a one-dimensional beam and later extended to plates by Reissner [36] and Mindlin [37]. The first order shear deformation theory (FSDT) used by Timoshenko, Reissner and Mindlin requires the use of a shear correction factor to approximate the distribution of transverse shear strain through the thickness. A FSDT has successfully been incorporated into the buckling analysis of sandwich structures and thick plates using Rayleigh-Ritz energy methods by Libove and Bartdorf [38], Dawe and Roufa-eil [39] and Ko and Jackson [40] by assuming shape functions for the shear strain in the xz- and yz-direction.

The FSDT provides accurate solutions to moderately thick plates and is therefore useful for practical cases of stiffened panels used in aerospace applications, however it is strictly only applicable to iso-tropic materials and may have significant error for very high thickness to width ratios. Recently the effect of transverse shear deformations on VAT plates was investigated by Groh et al. [33] by extending the equivalent single layer approach of Weaver and Cosentino [41]. VAT plates were found to be more affected by transverse shear than corresponding homogeneous quasi-isotropic laminates.

To the best of the authors' knowledge, the buckling of a VAT laminate at a structural level remains unexplored. The current contribution extends the work of Wu et al. [7] to develop an analytical model with a generalised Rayleigh-Ritz approach to solve the pre-buckling and buckling problem of a blade stiffened VAT panel including transverse shear deformations. Prebuckling analysis is first required to determine the varying stress field in the skin and constant stress in the stiffener. Restricting the model to prismatic sections enables the skin and stiffener to be treated in isolation for prebuckling. The buckling analysis is performed with two approaches for modelling the stiffener; a beam stiffener model and a plate stiffener model. In both cases the effect of transverse shear in both the skin and stiffener and the axial loading applied to the stiffener is explored and quantified. It should be noted, that although the method presented in this paper is referred to as analytical it is, strictly speaking, a semi-analytical method; analytical in formulation but requiring numerical integration. The numerical routine for the pre-buckling and buckling analysis was implemented in MATLAB R2012a.

The paper is structured as follows. Section 2 provides an overview of the analytical method including the assumptions and boundary conditions. Section 3 introduces the VAT orientation distribution. Section 4 details the prebuckling analysis of the VAT skin and straight fibre blade stiffener. Section 5 details the buckling analysis of the stiffened panel with two approaches for capturing the stiffener behaviour; a beam stiffener model and a plate stiffener model. Section 6 outlines the FEA model developed and used for comparison and validation of the analytical method. Section 7 presents results and a discussion of the analytical model and FEA for different stiffened panel configurations showing the applicability of the model to realistic stiffened panel configurations with thick sections and finally the paper is concluded in Section 8.

2. Analysis overview

During flight an aircraft wing is subject to bending resulting in the upper wing cover experiencing compressive loading which can be approximated by uniform end-shortening. In reality, a linear increase in compressive strain from the tip of the stiffener to the skin's outer surface would be present, however as the distance of the stiffened skin from the wing box global neutral axis is significantly larger than the depth of the stiffened panel this variation is considered negligible. Wing covers in general are supported by spars in the longitudinal direction and ribs in the transverse direction as shown in Fig. 2. The restraint on the stiffened panel by the spars and ribs is complex, however for simplicity they are both assumed to provide a simply supported boundary condition (preventing out-of-plane displacement of the skin) resulting in conservative results. In this study, the spars are additionally assumed to prevent any translation in the y-direction hence inducing biaxial loading in the panel. Only the skin is assumed to be connected to the ribs which provide a simply supported boundary condition, the stiffener is free to rotate about the skin's midplane at the location of the ribs.

Wing covers supported between spars and adjacent ribs are generally wide and contain several equally spaced stiffeners. When local buckling occurs, under compressive loading, repeating displacement modes occur between stiffener elements thus allowing the entire stiffened panel to be modelled by a representative section containing a single stiffener element and half a stiffener bay either side [42] as shown in Figs. 1 and 2.

A symmetry condition is required to be enforced along the skin's longitudinal edges to model the repeating mode shape, this is achieved by setting

Sw/Sy — 0

Cyz = 0

where w is the out-of-plane displacement of the skin in the z-direction and yyz the transverse shear of the skin in the yz-direction. This new boundary condition is henceforth referred to as the symmetric boundary condition. It should be noted that the use of this representative section with the symmetry condition is not valid for in-plane shear loading cases or skin laminates with extension-bending (B-matrix), extension-shear (A16,A26) or bending-twisting (D16,D26) coupling.

When global buckling occurs the representative section no longer represents a repeating unit as the boundary condition from the spar creates a shallow curvature in the y-direction, Ky (Fig. 1a). Global buckling behaviour is, however, dominated by x-direction curvature, Kx, due to the energy required to bend the stiffener, and the shallow Ky has minimal influence on the buckling load. Hence, the representative section used for the local buckling is also

valid for global buckling. Despite not representing global Ky, local Ky behaviour between stiffener elements is still captured in global buckling modes.

The use of a representative section significantly reduces the problem complexity whilst maintaining sufficient detail to allow results to be applicable to full, multi-stiffener, panels. Reduction of the problem complexity additionally allows a deeper understanding and physical insight to be obtained for buckling of stiffened VAT panels. Despite the suitability and advantages of the representative section with the symmetric boundary condition it seldom appears in literature.

Boundary conditions for the representative section are summarised in Fig. 2 and can be split into prebuckling boundary conditions and buckling boundary conditions. Further details are provided for the prebuckling and buckling boundary conditions in Sections 4 and 5 respectively.

3. VAT laminates

In this study, the blade stiffener laminate is constrained to straight fibres only and VAT laminates are only considered for the skin with the fibre orientation variation in the y-direction. The non-linear fibre orientation for each ply is expressed using the Lagrange polynomials method proposed by Wu et al. [7] in the form,

h(y>-£ ^ni^j)

where yj and yn are the y-coordinates of reference points and the coefficient of each term, Tn, is the fibre angle at the specific reference point, yn. For simple interpretation of results all cases used for model validation only consider a linear variation of fibre angle and hence the fibre orientation reduces to

h(y) = Tc + 2(Ti - To)jj (3)

where T0 and T1 are the fibre orientations of the skin at the location of the stiffener and symmetric boundary condition respectively and b is the width of the representative section (distance between stiffeners) as shown in Fig. 3. The fibre orientation of a VAT ply is designated by (To|Ti).

4. Prebuckling analysis

Herein, the fibre orientation is limited to variations in the y-direction only and the stiffened panel is therefore prismatic. For the loading case of end-shortening no coupling or interaction effects exists between the skin and the stiffener and they can be

Fig. 2. Representative section in stiffened panel analysis. Coordinate systems shown are the local skin or global (xyz) and local stiffener (x'y'z'). Boundary conditions and loadings with a * are only used for determining the prebuckling stress field and are removed for the buckling analysis.

Fig. 3. VAT linear fibre angle variation for <0°| ±45°). T0 is the fibre variation along x = 0 (centre line underneath stiffener) and T1 along y = ±b/2 (panel longitudinal edges).

treated separately in the prebuckling analysis. The approach for determining the skin's stress distribution used herein is the same as in Wu et al. [43] where the in-plane equilibrium equations are expressed using Airy's stress function and the condition of compatibility and prescribed displacement boundary conditions are satisfied through minimisation of the total complementary energy. A brief overview of the procedure is now provided, full details can be found in [43].

For the prebuckling analysis the skin's loaded transverse edges are subject to uniform end-shortening and the boundary conditions are

u(x = a/2) = -Ax/2 u(x = -a/2) = Ax/2

where a is the length of the panel in the x-direction (distance between rib bays), u is the skin in-plane displacement in the x-direc-tion and Ax the end-shortening applied to both the stiffener and skin. Similarly, the same uniform end-shortening is applied to the stiffener's loaded transverse edges which have the boundary conditions

u'(x' = a/2) = -Ax/2 u'(x' = -a/2) = Ax/2

where u0 is the stiffener in-plane displacement in the x0-direction. The skin's longitudinal edges, y = ±b/2, are constrained in y-direction translation, inducing biaxial compression and the boundary condition applied is

v (y = ±b/2) = 0

where v is the skin in-plane displacement in the y-direction. The stiffener's longitudinal edges along the y0-direction are free to expand.

The VAT skin considered is confined to laminates with null B-matrix terms. The total complementary energy of the VAT skin can be expressed using Airy's stress function [44] as

Ps = 1

_ d2 U d2U 2a12 ôy2+a22

( d2f\2 _ d2U d2U I-I - 2a-=---

' dy2 dxdy

d2U d2U

b/2 ôF u dxdy

a/2 d2u @xu v d2u

a/2 - dxdy

d2u d2u dydx

dx2 dxdy

b/2 d2 U wu d2U

b/2 - dxdy

a/2 d2U dU v d2U

a/2 - dxdy

x=-a/2

y=-b/2

where aij are terms of the skin a = A-1 matrix and U is Airy's stress function. The in-plane displacements, u and v, along the panel boundaries in Eq. (7) are defined as per the in-plane loading and boundary conditions provided in Eqs. (4)-(6).

The total complementary energy is expressed in normalised coordinates, n = 2x/a and g = 2y/b, and the stress function, U, is assumed to have the following form:

U(n, g) = Uo (n, g) + Ui(i, g)

where U0 satisfies the stress distribution along the boundaries and U1 the stress distribution in the interior region. The form of U0 is assumed to be

Uo(n, g)=/l(f)+/2(g)

= £ CtWk(g)

= Z dtwd(n)

where K is the number of terms in the series, ck and dk are coefficients of the stress function along the boundaries and wck and wkd are admissible functions, here Legendre polynomials are chosen for the admissible functions. Legendre polynomials were chosen because they enable localised behaviour, due to VAT, to be captured accurately with less terms relative to periodic trigonometric functions [45]. The interior region's stress function is expressed in the form:

P-1 Q-1

Ui(«, g) = £E /pqXp (n)Yq(g)

p=0 q=0

where P and Q are the number of terms and Xp and Yp are admissible functions in the x and y-directions respectively. To satisfy the stress free condition for U1 Legendre polynomials are used with circulation functions [7,46,47] for Xp and Yp,

Xp(n) = (1 - n2) Lp(n) Yq (g) = (1 - g2)2Lq(g)

U// U/c U/d

Uc Ucc Ucd

U/d uTd Udd

c = Px0

where Li is the ith term of the Legendre polynomial. Substituting Eqs. (8)-(14) into Eq. (7), evaluating the integrals with numerical integration and minimising the total complementary energy with respect to the unknown coefficients the following set of linear equations, expressed in matrix form, are obtained,

where for example, U/c is the factor of the unknown coefficient c that arises from minimising the total complementary energy with respect to / and Px0 is the constant that arises when minimising the total complementary energy with respect to c. Further details regarding the terms in Eq. (15) can be found in Wu et al. [43]. Solving for the coefficients the prebuckled stress distribution of the VAT skin subject to uniform end-shortening is obtained.

As the stiffener is constrained to contain only straight fibre laminates the prebuckled stress distribution is constant. The free edge ensures no biaxial stress state is induced in the stiffener and the prebuckling stress resultant in the x-direction, Nx>st., is simply obtained using the stiffener laminate equivalent Young's modulus,

Ex,stAstAx

where Exst: is the equivalent Young's modulus of the stiffener laminate in the x-direction given by Ex>st. = 1/(a11>st.fst ), Ast: the

cross-sectional area of the stiffener in the yz-plane and h the height of the stiffener blade.

5. Buckling analysis

The Rayleigh-Ritz energy method is used to solve the buckling problem for the stiffened panel using the stress distribution obtained in the prebuckling analysis. The skin is modelled as a thick plate and the stiffener is considered with two approaches; a beam model and a plate model.

The boundary conditions applied to the skin are identical for both the beam model and plate model. In both cases the prebuckling boundary conditions on u, uU and v are removed for the buckling analysis. During the buckling analysis the skin's loaded transverse edges, x = ±a/2, are simply supported and constrained to have null transverse shear in the yz-direction such that

w(x = ±a/2) = yyz{x — ±a/2) = 0

where yyz is the transverse shear of the skin in the yz-direction. The skin's longitudinal edges, y = ±b/2, are subject to the symmetric boundary condition where rotation along the y-direction and transverse shear in the yz-direction are null, hence we have

dW (y = ±b/2)=yyz(y = ±b/2) = 0

The total potential energy of the skin is the sum of the bending, transverse shear (xz- and yz-direction) and in-plane potential energy,

PTPE.skin — (Pbend. + Ptrans. shear + Pin-plane)

Due to the thickness, tsk:, to width, b, ratio of practical stiffened panels potentially being as large as 1/20 transverse shear effects can be significant [33]. A FSDT is included in the analysis by using a reduced bending energy term and a transverse shear energy term in the total potential energy [38-40,48]. Eq. (19) expanded for the case of a VAT laminate with a FSDT is:

1 fa/2 fb/2 PTPE.skin — 2

rb/2 4d (d2w дУЛ2 (d2w дуЛ

J-b/2 Dl\äx2+D22[w-ЦТ)

(@2w_ @уД дуЛ

у дх2 дх) ^ду2 ду )

@2W ду^ дУуЛ

+D66{2дХдУ-^х)

+2D дУУхЛ(2д_ъ«_дуЛ

+ 16 ^ дх2 дх) дхду ду дх)

J -a/2 . k ia/2

@2w дуп ду2 ду

a/2 tb/2 ^ ,

@х@у @у

@2W дУх;

дУУп дх

-b/2 fb/2

х у у

Су; У^)]^

АхАу

-a/2 -b/2

у 2 х х2 у

д2Ф @w@w

х у х у

АхАу (20)

where w is the deflected shape of the skin, yxz and yyz the shear strain in the xz- andyz-directions respectively, Dj the plate bending stiffness matrix terms which vary over the skin, k the Timoshenko shear factor which is taken as 5/6 for rectangular sections [39,48], Gxz and Gyz the transverse shear stiffness in the xz- and yz-directions respectively and k the loading factor. The total potential energy is then expressed in normalised coordinates. The unknown function w, is represented by a series expansion containing Legendre polynomials and circulation functions to enforce boundary conditions,

M-1N-1 W — / 7 Amn m—0 n—0

— £EAmn [(П2 - 1)Lm]

/(g2 - 1)^ + 1

where M and N are the number of terms in the x- and y-directions respectively and Amn the unknown coefficients. The constant of integration of the indefinite integral is zero. This series expansion ensures null out-of-plane displacement at the transverse edges (П = ±1) and null rotation, dw/dg = 0, along the longitudinal edges (g = ±1).

The unknown functions for the transverse shear, yxz and yyz, are similarly given by:

E-1 F-1 Cxz = £ YBefLeLf

e—0 f— 0

У у; — £ £Cgh [(n2 - 1)Lg ][(g2 - 1)Lh)]

g— 0 h— 0

where for yxz, E and F are the number of terms in the x- and y-directions respectively and Bef the unknown coefficients, and for yyz, G and H are the number of terms in the x- and y-directions respectively and Cgh the unknown coefficients.

(18) 5.1. Beam stiffener model

The total potential energy of the panel when modelling the stiffener as beam is given by the sum of the skin and stiffener contributions,

PTPE — PTPE,skin + PTPE,beam-stiffener

where the total potential energy of the beam stiffener is expanded to

PTPE,beam-stiffener = № ,bend. "" PAG, trans, shear "" PGJ,tors. "" PEA,in-plane)stiffener

The beam stiffener's energy terms in Eq. (25) represent an equivalent beam that lies on the skin's midplane and are the bending, transverse shear (xz-direction), torsional due to twisting and potential due to in-plane loads. The skin and the stiffener are assumed to be rigidly attached and no slipping is allowed between the two components. Hence, the boundary conditions applied to the skin along this attachment line translate to the stiffener when modelled as a beam. The beam is simply supported at its ends, x = ±a/2, allowing the stiffener to rotate about the skin's midplane for global buckling modes.

The bending potential energy of a Timoshenko beam is given by

1 fa/2 PEI.bend — о /

2 J-a/2

Ex.st.ht

@2w @ух

у—0

where Ex>st..ist. is the flexural rigidity of the beam about the y-axis and yxz st the transverse shear deformation of the beam in the xz-direction. The displacement of the stiffener is constrained to be equal to the displacement of the skin's midplane, however, the normal to the midplane rotation, /xzst., and transverse shear displacement, yxz st are free. In reality, the transverse shear displacement at the interface of the skin and the stiffener must be equal but for a FSDT an average over the depth of the section is considered and this is not required to be equal for the skin and stiffener.

Determination of the stiffener's equivalent Ex>st./st. requires estimation of local neutral axis location. The neutral axis lies somewhere between the midplane of the skin and the mid-height of the stiffener depending on the relative in-plane and flexural stiffness of the skin and stiffener [49]. For most practical cases the location of the stiffener neutral axis is only slightly above the skin's

д2Ф Î6w\2 д2Ф/âw\2

midplane and the assumption that the neutral axis lies on the skin's midplane is valid [25]. Thus,

Ex.st. ^st — Ex.

The transverse shear strain energy of the beam is given by

AG.trans . shear :

1 iaß

2 J-a/2

fcAst. GxZ

tT^st.]

where Gxz,st. is the transverse shear stiffness of the stiffener in the xz-direction. The value of GXz,st. is in the xz-direction for the global coordinate system, however when considering the local coordinate system of the stiffener as a laminate it is the in-plane shear stiffness Gx'y st . where x! and y' are local coordinates of the stiffener laminate (Fig. 2). The layup of the stiffener web significantly influences the stiffener's ability to resist global transverse shear deformation, Gxzst. can vary from 5 GPa ([0°]n) up to 50 GPa ([±45°]n) for a typical aerospace grade prepreg. Additionally, the section over which this shear acts is very stubby, transverse shear deformations for composite plates is considered important for ratios larger than 1/20, for the case of blade stiffeners the ratio can be larger than 5/1 and transverse shear effects in the xz-direction can significantly reduce global buckling loads.

The torsional restraint of the stiffener is taken into account by treating the stiffener as a De Saint Venant torsion bar [48] and determining the energy due to the beam rotation,

1 ra/2 /_, fd fdw

PGJ,tors. =~ / GJ ^ ^7- Cy

dx\ dy 'yz

where GJ is the effective torsional restraint. For thin blade stiffeners (tst /h < 1/10) GJ is given by

a-g tst.h

GJ — Gxz,st. 3

For thick blade stiffeners GJ is given by Nemeth [50] as

GJ — Geq . ,stJeq . ,st. (31)

Geq . — \J Gxz,st. Gxy,st.

j _ I "xz,st tst.h Jeq. — o

Jxy ,st.

p—1,2,3,.

^96 tst.

p5 h \G,

1 -(-1)p

tanM ^(G^

2tst . \Gxz,st .

where xy- and xz-directions are in the global coordinate system. Ne-meth's formulation for determining GJ takes into account shear deformation in the xy-direction of the stiffener and is used, unless otherwise stated, for all cases presented in this study.

The above expressions for GJ are strictly only valid when no axial load is applied to the stiffener, to account for an axial stiffener load a reduction factor is applied to GJ [15-17],

GJ red . — GJ X r

^crit. ,sk . ,s . s

^crit .,;

where ecrit st s :S: and ecrit sk: s : s: are the buckling strains of the stiffener blade and skin respectively assuming both are simply supported along the skin-stiffener attachment line and r is the torsional reduc-

tion factor. The stiffener buckling strain is calculated with the following closed-form solution [51]

Ex.st. tst

nWavesp2D1 a2

where nwaves is the number of half waves in the longitudinal direction and D11st. and D66 st. are in the stiffener's local coordinate system (x'y'z'). Determination of ecrit>sk. s.s. is not possible with a closed-form solution and a full prebuckling and buckling analysis on the unstiffened VAT skin is required. However, all stiffness matrices computed for the unstiffened VAT skin are reused in the full stiffened panel analysis and the computational costs due to this additional step are minor.

The potential energy due to the in-plane loading is simply given

EA,in-plane .

k r E A

2 L/2 Ex st Ast. läx

where Ex;st Ast is the stiffener axial stiffness.

The introduction of a beam into the model only requires one additional shape function, compared to an unstiffened panel, for the stiffener transverse shear in the xz-direction, cxzst , which is given in series expansion by

— YPtLt

where T is the number of terms and Dt the unknown coefficients. Substituting Eqs. (19)-(23)and (25)-(38) and the prebuckling solutions, U and Nx,st, into Eq. (24), evaluating the integrals with numerical integration and minimising with respect to all coefficients of the four shape functions, Amn, Bef, Cgh, Dt, a set of linear equations is obtained which are expressed in matrix form,

[(Ksk.+ Kst.)+k(Lsk.+ Lst.)]A = 0 (39)

where Ksk contains bending and transverse shear stiffness matrices of the VAT panel, Kst the bending, transverse shear and torsional stiffness matrices of the stiffener, Lsk and Lst the stability matrices of the skin and stiffener due to the in-plane stress fields respectively and the vector A the coefficients of all the shape functions,

A —[Amn Bef Cgh

The critical buckling load is given by the lowest eigenvalue of Eq. (39) which is then used with the prebuckling solutions to determining the buckling load.

5.2. Stiffener plate model

The stiffener plate model is only applicable to local buckling modes, with the attachment line connecting the skin's midplane and stiffener's midplane constrained in out-of-plane displacement. Modelling the stiffener as a plate the total potential energy is expressed as the sum of the contributions from the skin and stiffener with an additional penalty term,

PTPE = PTPE,skin + PTPE,plate-stiffener + PTPE, penalty (41)

The penalty term is required to enforce compatibility of rotation at the line of attachment between the skin and stiffener plate elements. The total potential energy of the stiffener, similarly to the skin, is the sum of the potential energy terms for the bending, transverse shear and in-plane loads,

PTPE,plate-stiffener = (Pbend . + Ptrans .shear + Pin-plane)stiffener (42)

As the stiffener is now being modelled as a plate additional boundary conditions are required. The loaded stiffener's edges, x' = ±a/2, are simply supported in the stiffener's local coordinate

crit. ,st. ,s . s

system, see Fig. 2, and constrained to have null transverse shear in the y z -direction. This boundary condition simulates the stiffened panel being connected to another stiffened panel section in the x-direction and is given by

wst.(x' = ±a/2) = yyV >st(x' = ±a/2) = 0 (43)

where wst. is the stiffener plate out-of-plane displacement in the z!-direction and yy,z, st the stiffener transverse shear in the y'z'-direction. Along the skin-stiffener attachment line, y = y' = 0, both the skin and stiffener are constrained in out-of-plane displacement

w(y = 0) = 0, Wst.(y' = 0) = 0

To force a node in the skin at y — 0 the constant 1 in the deflected shape function, Eq. (21), is replaced with a 0.

The plate stiffener total potential energy, Eq. (42), is expanded

@Cx'z' ,st.

TPE,plate-stiffener :

1 ia/2 ih D (@2Wst.

2 J-a/2 Jo Dl1st (v^x^ ^"âx^)

i a2wst. _ dyy,z, ,SL \2 . ^ dy02 dy' J

(@Wst. @yx-z-,stA /@2Wst. @yy-z-,st. A

^ dx'2 dx' J ^ dy02 dy' J

Ô2Wst.

__@yx'z' ,st.

dx' dy' dy'

d2wst dyx z ,st

dyyz' ,st. Y

' dx' J

( d2 Wst. dyx'z' ,st.A

^ dx'2 dx' )

d2Wst.

+ 2D26,st

dyx'z' ,st.

dyy'z' ,st.A

( d2 Wst. dyy'z' ,st. A

^ dy'2 dy' )

d2Wst.

dyx z ,st

y z ,st

,a/2 rh

1 ra/2 fh r /

1 /-a/2 70 M^'^

a/2 Jo

kN,- ~ ra/2

dx'dy'

Gy'z',st.yjz',st.) ]dx'dy'

"x',st. + 2

-a/2 0

where yxz st and yy,z, st are the shear strain in the x'z'- and y'z'-directions respectively, Dj,st. the stiffener plate bending stiffness matrix terms (constant) and Gxz,,,st. and Gyz,st. the transverse shear stiffness in the x'z'- and y'z'-directions respectively, all in the stiffener's local coordinate system. Similarly to the analysis procedure for the skin, the total potential energy is expressed in normalised coordinates, t! — 2x'/a and g' = y'/h, and the unknown functions wst., yxzst and yy,z, st are expanded into the following series, satisfying the required boundary conditions,

i-1 J-1

Wst. = ££Xj[(t'2 - 1)Li][g'Lj] i=0 j—0

'x'z',st. — £y^bpzLPLz D—0 z=0

y z ,st

— EEsuv[(t2 - 1)Lu][g'Lv]

where I - J, ! - Z and U - V are the number of terms in the x' -and y'-directions for the Wst., yxzst and yy,z,st series' respectively and % btz and suv are the unknown coefficients for the Wst , yx z ,st and yy,z, st series' respectively.

The penalty term in Eq. (41) is analogous to a torsional spring located along the skin-stiffener attachment line [24,25] and is given by

penalty —

penalty

Sw ¥ - yy

_ ( SWsst.

, V sy

y'z ,st .

dx (49)

where kpenaity is the spring torsional stiffness. The value of kpenaity was determined by increasing the order of kpenaity until convergence was achieved, in this study 1. 0 x 106 N was found to be sufficient to force the compatibility condition.

The solution procedure hereafter is the same as for the beam stiffener model. Substituting Eqs. (19)-(23)and (42)-(49) and the prebuckling solutions, U and Nxst., into Eq. (41), evaluating the integrals with numerical integration and minimising with respect to all six coefficients of the shape functions, Amn, Bef, C

j , lJUZ, Luv:

mn; Bef; Cgh,

a set of linear equations in the matrix form of Eq. (39) are obtained, where Kst now contains the bending and transverse shear stiffness matrices of the stiffener modelled as a plate and Lst. is the stiffener stability matrix due to the in-plane stress field. The vector A contains the coefficients of all the shape functions,

Bef Cgh Xij btz su

The critical buckling load is given by the lowest eigenvalue of Eq. (39).

6. Finite element analysis

FEA for the prebuckling and buckling analysis was performed with Abaqus. A script was developed to allow specification of fibre orientation for individual elements as per Eq. (2), simulating a VAT laminate. The S4R shell element was chosen for the skin and the compatible S4 shell element for the stiffener. The S4 element was required for the stiffener due to S4R elements experiencing hour-glassing when subject to in-plane bending as is the case for stiffener global buckling.

Uniform end-shortening to the skin and stiffener transverse edges was applied as stress perturbation only, thereby allowing the stiffener to rotate about the skin's midplane during the buckling analysis. The boundary conditions applied as stress perturbation only (prebuckling) were

u(x = a/2) = u' (x' = a/2) = -Dx/2 u(x = -a/2) = u' (x' = -a/2) = Dx/2 v (y = ±b/2) = 0

For the buckling analysis all stress perturbation (prebuckling) boundary conditions were relaxed. Simply supported boundary conditions were applied on the loaded skin transverse edges and additionally the rotation of the normal to the midplane in the yz-direction, /yz, was set to zero. Boundary conditions for the transverse shear profiles were unable to be directly assigned in Abaqus, instead the rotation of the normal to the midplane in conjunction with the out-of-plane displacement along edges were used to achieve the same effect. The boundary conditions along the skin transverse edges were

w(x — ±a/2) — /yz(x — ±a/2) — 0

The skin's longitudinal edges, y — ±b/2, are subject to the symmetric boundary condition where /yz is set to zero

r-1Z-1

Фуг(У = ±b/2) = 0

7. Results and discussion

is equivalent to setting dw/dy to zero. The loaded stiffener's edges, X = ±a/2, are simply supported in the stiffener's local coordinate system and constrained to have null rotation of the normal to the midplane in the y'z'-direction, /y,z,st and the boundary conditions are

wst (X = ±a/2) = фуг, st (X = ±a/2) = 0

The prebuckling stiffness and buckling loads were determined by summing nodal forces along the skin and stiffener shell edges at the location of the applied displacement.

Table 1

Carbon fibre unidirectional prepreg mechanical properties used for case studies.

Eii (GPa) E22 (GPa) V12 (-) G12 (GPa) G13 = G23 (GPa)

161 11.38 0.32 5.17 3.98

Fig. 4. Normalised in-plane stress resultant distribution along skin and stiffener edges for VAT <0°| ±45°) skin and 0° stiffener case. Analytical results are obtained using K = 5 and K = 10 for the number of terms in the stress function, Ф0, along the boundary. Due to the sections being prismatic there is no variation in the central region and results are independent of the number of terms of P and Q. The following stress resultants not shown on the plot are: Nxy,sk. = Ny,st. = Nxy,st. = 0.

The stiffened panel dimensions used for model validation were based on Nagendra et al. [29,30]; width (distance between stiffener bays) 200 mm, length 750 mm and stiffener height 60 mm. Material properties of the carbon fibre unidirectional prepreg selected for the study are provided in Table 1. A FEA mesh density of 320 x 60 x 20 elements in the x-, y- and z-directions respectively was selected to achieve converged results. Four layups were considered for the skin, ±45°, quasi-isotropic (QI), VAT (0°|±45°) and VAT (0°| ±90°) and three layups for the stiffener ±45°, QI and 0°. Skin and stiffener thicknesses were fixed by first designing a ±45° skin and QI stiffener panel aiming for local and global buckling to occur at Nxsmeared и 1. 0 kN/mm where Nx smeared is the total load taken by the skin and stiffener per unit panel width. The chosen thickness for the skin and stiffener was 8 mm.

The layups used for the VAT and ±45° laminates were specially orthotropic [±6, T^l^s consisting of eight 1 mm thick layers, thus eliminating extension-shear (Ai6,A26), bending-stretching (Bj) and bending-twisting couplings (Di6,D26). For QI laminates equivalent smeared properties were used to remove stacking sequence dependence.

The analytical model initially used 5 terms for prebuckling shape functions (P - Q - K), 7 terms for skin and stiffener out-of-plane deflected shape functions (M - N, I - J) and 5 terms for skin and stiffener transverse shear shape functions (E - F, G - H, T,Y - Z, U - V). However, for cases where convergence was not achieved the number of terms was increased.

The speed of the current semi-analytical method is implementation and platform specific, however, for all cases in this paper the computation of the results with the analytical method implemented in MATLAB R2012a was found to require less computational time than Abaqus FEA. The number of degrees of freedom used in the analytical model for the buckling analysis ranged from as little as 49(M - N = 7) for straight fibre cases neglecting transverse shear with the beam stiffener model up to 1026(M - N, I - J = 15 and E - F, G - H,Y - Z, U - V = 12) compared to the FEA having 156,006 degrees of freedom.

7.1. Prebuckling

The analytical prebuckling stiffness for all configurations was found to be in close alignment with FEA with less than 0.3% error in all cases. Due to all cases being prismatic, the in-plane forces per unit length in the y- and xy-direction, Ny and Nxy respectively,

For the cases neglecting transverse shear effects setting /yz to zero

Table 2

Comparison of local buckling loads between analytical models and FEA.

Layup FEA Analy.: Beam with GJ Analy.: Beam with GJred Analy.: Full plate model

local local r Nx,local local

Skin Stiff. (kN/mm) (kN/mm) Error (-) (kN/mm) Error (kN/mm) Error (%)

±45° ±45° 0.96 0.97 +0.99% 0.95 0.96 +0.63% 0.96 -0.06

±45° QI 1.14 1.16 +1.80% 0.70 1.14 +0.43% 1.14 -0.01

±45° 0° 1.53a a a a a a 1.53a +0.03a

QI ±45° 2.34 2.38 +2.05% 0.91 2.37 +1.49% 2.33 -0.23

QI QI 2.53 2.70 +6.77% 0.47 2.55 +1.16% 2.52 -0.07

QI 0° 2.49a a a a a a 2.51a,b +0.92a,b

<0°|±45°) ± 45° 1.88 1.94 +3.03% 0.96 1.93 +2.20% 1.90b +0.70b

<0°|±45°) QI 1.91 2.01 +4.99% 0.75 1.94 +1.65% 1.92 +0.39

(0°| ±45°) 0° 1.98a a a a a a 1.99a,b +0.31ab

<0°| ±90°) ± 45° 2.69 2.76b +2.57%b 0.89b 2.73b +1.15%b 2.68b -0.43b

<0°| ±90°) QI 2.92 3.13b +7.15 %b 0.34b 2.96b +1.45%b 2.92b -0.01b

<0°| ±90°) 0° 2.10a a a a a a 2.14a,c +1.70ac

Error: <8% <3% <2%

a Local buckling of the stiffener blade occurred prior to skin buckling, this buckling mode is not predicted with the beam model. b Number of terms increased for convergence, PQK = 10, MNIJ = 10 and EFGHYZUVT = 8. c Number of terms increased for convergence, PQK = 12, MNIJ = 15 and EFGHYZUVT = 12.

are constant over the planform and only a variation in Nx occurs along the y-direction. Fig. 4 shows the Nx distribution over the skin planform for the VAT (0°| ± 45°) skin and 0° stiffener case for the analytical model and FEA. The analytical model error for constant in-plane stress resultants was less than 0.3%. Whilst the error in the maximum Nxskin for the VAT laminates was less than 3% and 1% for 5 and 10 terms in the prebuckling shape functions respectively.

7.2. Local buckling

Local buckling results for different configurations are provided in Table 2 for the beam stiffener model (with and without the GJ reduction factor), the plate stiffener model and FEA. The first local buckling mode shape for the VAT (0°| ±45°) skin and QI stiffener panel is provided in Fig. 5 for FEA and analytical models. The first local buckling mode for all 0° stiffener cases was due to stiffener blade buckling and is not captured with the beam stiffener model. The skin and stiffener transverse shear profiles obtained with the plate stiffener model and FEA are provided in Figs. 6 and 7 respectively for the VAT (0°| ± 45°) skin and QI stiffener panel first local buckling mode. Qualitatively it can be seen that the analytical model accurately captures the FSDT profiles compared to the FEA. The transverse shear profile in the y'z'-direction for the stiff-ener was largely dominated by numerical noise for both the analytical model and FEA.

The beam model using the thick laminate and unreduced, due to axial loading, GJ (Eq. (31)) has less than <8% error in local buckling for all cases, the highest error being for cases where the stiffener blade's buckling eigenvalue, as per Eq. (36), was close to the skin's critical eigenvalue as indicated by values of r less than 1 in Table 2. In these cases, the consideration of the stiffener blade having a tendency to buckle out-of-plane in the local coordinate system is not considered. Including the reduction factor and using GJred: (Eq. (35)) in the analyses reduced the local buckling error for all cases below 3%, those with the lowest r and corresponding highest error being most affected. The analytical model, considering both the skin and stiffener as plates, predicted all local buckling loads to within 2% error and was able to capture buckling modes originating from stiffener blade buckling. The plate stiffener model has an increased computational cost compared to the beam stiffener model due to the additional shape functions and larger stiffness matrix. However, the plate stiffener model is expected to be accurate for a wider range of configurations and different plate boundary conditions when compared to the beam stiffener model and can allow VAT laminates in both the skin and stiffener to be explored.

To gain an understanding of the magnitude of the stiffener tor-sional restraint, skin transverse shear and stiffener transverse shear effects on the local buckling load, analyses were performed

Fig. 6. Transverse shear profile for the first local buckling mode of the VAT <0° | ± 45°) skin and QI stiffener panel in normalised coordinates. (a) Analytical plate stiffener model yxz. (b) FEA yxz. (c) Analytical plate stiffener model yyz. (d) FEA yyz.

with FEA, removing in turn each of the aforementioned effects with results provided in Table 3.

The effect of removing the torsional restraint, simulated by replacing the stiffener with a simply supported line, decreased the local buckling load by up to 18%. All analytical models, being identical for null torsional restraint, captured this behaviour within 1% of FEA results. This indicates that the current model is able to accurately capture the buckling of simply supported thick VAT plates where shear deformations are significant and any error in the full model with the stiffener is not due to FSDT in the skin.

Removing the FSDT for the skin or stiffener was achieved with FEA by using a combination of S3 (thick shell) and STRI3 (thin shell) elements in Abaqus. Removing transverse shear of the skin

Fig. 7. Stiffener x'z'-direction transverse shear, yx,z.sti, profile for the first local buckling mode of the VAT <0°| ±45°) skin and QI stiffener panel in normalised coordinates. (a) Analytical plate stiffener model. (b) FEA

only increased the buckling load by up to 4% and was captured accurately by the beam stiffener model, with GJred., and plate stiffener model within 2% and 1% respectively. However, the beam stiffener model error using the unreduced GJ was up to 7%. Removing the FSDT from the stiffener resulted in an increase in buckling load by up to 6% for the FEA. Removing transverse shear from the beam stiffener model was achieved by using the thin laminate expression for GJ provided in Eq. (30). The beam stiffener model using the reduced GJ was able to match FEA within 4% when stiffener transverse shear was neglected, whilst the plate stiffener model within 1%.

All three effects investigated (torsional restraint, skin transverse shear and stiffener transverse shear) introduce significant error if disregarded and should be captured for practical cases of stiffened panels. The beam stiffener model using GJred. has less error in all cases compared to the beam stiffener model using the unreduced GJ, indicating that neglecting stiffener loading when considering torsional restraint can be a significant source of error.

The majority of the local buckling error for the beam stiffener model is expected to be due to difficulty in accurately capturing the stiffener torsional restraint. More specifically, determining a value for GJ which is valid for axially loaded stiffeners with transverse shear deformations. Nemeth's GJ formulation (Eq. (35)) for thick sections used in conjunction with a linear reduction factor due to axial loading provides accurate results for the cases

considered, however, it is expected that the validity of this simple method may not hold for a wider range of cases. The stiffener plate model is expected to be applicable for a wide range of cases, albeit with increased computational costs.

7.3. Global buckling

Global buckling analytical results are obtained using the beam stiffener model which captures out-of-plane (global z-direction) stiffener displacement by replacing the stiffener with a one-dimensional beam attached on the skin's midplane. During prebuckling analysis both the skin and stiffener are subject to the uniform end-shortening, however for the buckling analysis the stiffener is free rotate about the skin's midplane. Global buckling results compared to FEA are provided in Table 4 and the first global mode shape for the VAT (0°| ± 45°) skin and QI stiffener case illustrated in Fig. 8.

For cases with matching stiffener laminates changing the skin laminate had little effect on the global buckling load, Nx,smeared, which differed by less than 9% between the skin laminates. However, the difference in buckling eigenvalue (end-shortening) differed between cases of fixed stiffener laminate by up to 90%. This

Table 4

Comparison of global buckling results between the analytical beam model and FEA.

Layup FEA (kN/mm) Analytical (kN/mm)

Skin Stiff. kglo. Nx,glo. Nx,glo. Error (%)

±45° ±45° 1.54 0.92 0.95 +3.21

±45° QI 3.46 2.54 2.67 +5.07

±45° 0° 3.72 3.92 4.13 +5.37

QI ±45° 1.19 0.94 0.97 +3.20

QI QI 2.91 2.70 2.87 +6.34

QI 0° 3.31 4.13 4.46 +8.00

(0° | ±45°) ±45° 0.80 0.97 1.01 +3.48

(0° | ±45°) QI 1.95 2.64 2.81a +6.54a

(0° | ±45°) 0° 2.31 3.87 4.23a +9.46a

(0° | ±90°) ±45° 1.28 0.95 0.98a +3.36a

(o° | ±90°) QI 3.15 2.75 2.94a +6.91a

(0° | ±90°) 0° 3.53 4.21 4.67a +8.46a

Error <10%

a Number of terms increased for convergence, PQK = 10, MN = 10 and EFGHT = 8.

Table 3

Effect of stiffener torsional restraint, skin transverse shear and stiffener transverse shear on FEA local buckling results. The value in the parenthesis is the error compared to the full FEA including all effects.

Layup Nx,local (kN/mm)

Skin Stiff. Full model No torsional restraint No skin transverse shear No stiffener transverse shear

±45° ±45° 0.96 0.90 (-6.23%) 0.98 (±2.74%) 0.98 (±2.16%)

±45° QI 1.14 1.10 (-2.83%) 1.17 (±2.66%) 1.14 (±1.02%)

±45° 0° 1.53a - 1.57 (+2.48%)a 1.53 (+0.14%)a

QI ±45° 2.34 2.07 (-11.58%) 2.43 (±3.93%) 2.40 (±2.70%)

QI QI 2.53 2.42 (-3.97%) 2.61 (±3.29%) 2.58 (±2.17%)

QI 0° 2.49a - 2.52 (+1.27%)a 2.57 (+3.14%)a

(0° | ±45°) ± 45° 1.88 1.56 (-17.02%) 1.91 (±1.26%) 1.99 (±5.38%)

(0° | ±45°) QI 1.91 1.74 (-8.93%) 1.93 (±1.18%) 1.96 (±2.80%)

(0° | ±45°) 0° 1.98a - 2.00 (+1.21%)a 2.00 (+1.15%)a

(0° | ±90°) ± 45° 2.69 2.41 (-10.46%) 2.77 (±2.84%) 2.80 (±3.82%)

(0° | ±90°) QI 2.92 2.86 (-1.84%) 2.99 (±2.49%) 2.98 (±1.98%)

(0° | ±90°) 0° 2.10a - 2.11 (+0.51%)a 2.19 (+4.45%)a

Error range from full FEA (-18% ! -1%) (0% ! ±5%) (0% ! ±6%)

Analytical error - beam stiffener (GJ) <1% <7% <7%

Analytical error - beam stiffener (GJred.) <1% <2% <4%

Analytical error - plate stiffener <1% <1% <1%

a Local buckling of the stiffener blade occurred prior to skin buckling. Analytical model errors are compared to FEA when neglecting effects.

behaviour indicates that whilst the stiffener alone primarily dictates the buckling load it is the overall in-plane stiffness of the panel (skin and stiffener) which must be considered when determining the buckling strain.

The inclusion of transverse shear deformation in the xz-direc-tion (global coordinate system) for the stiffener significantly affected results. The transverse shear modulus (xz-direction) of the stiffener in Eq. (28) is the local in-plane shear modulus (x'y'-direc-tion) of the stiffener laminate and hence is dependant on stacking sequence. For the material system used this shear modulus can range from 5 GPa for a 0° laminate up to 42 GPa for a ±45° laminate. Maximising Ex>st.Jst. for global buckling, for fixed geometry, is achieved by increasing the proportion of 0° plies in the laminate and hence transverse shear deformations can become significant. For the case of 0° stiffeners the transverse shear deformation reduces the buckling load by up to 50%.

The assumption of the stiffener's neutral axis being located at the skin's midplane is expected to be the main source of error in global buckling results. In reality, it can be located anywhere between the skin's midplane and mid-height of the stiffener. For fixed dimensions and skin laminate the neutral axis shifts away from the skin's midplane with increasing Exst. resulting in an increased error in the global buckling load for the current model. This is evidenced by switching, for any fixed skin laminate, from a ±45° to QI to 0° stiffener laminate, i.e. increasing Exst., in Table 4. Improved estimations of Exst./st. may be obtained by extending the approach of Seide [49] who considered the relative stiffness and flexural rigidity of the skin and stiffener to determine an effective Exst./st. for isotropic stiffened panels. It should be noted that the 0° stiffeners which have the largest error are not practical designs due to the premature buckling of the stiffener blade.

8. Conclusion

In this paper, an analytical model has been presented for the prebuckling and buckling analysis of novel blade stiffened VAT panels utilising the Rayleigh-Ritz energy method. The buckling analysis includes a FSDT and is applicable to practical stiffened panel configurations containing thick sections. The boundary conditions considered were simply supported for the skin transverse edges and symmetric, with null y-translation, for the skin longitudinal edges. In all cases, the stiffened panel is subject to uniform end-shortening.

Prebuckling was performed for the skin and stiffener in isolation, minimising the total complementary energy for the skin and using a simple Young's modulus relationship for the stiffener and was found to be in good agreement with FEA. The buckling analysis was performed by minimising the total potential energy of the stiffened panel. The stiffener was modelled with two approaches; as a beam element and a plate. The beam stiffener model is applicable to both local and global modes whilst the plate stiff-ener model is only valid for local modes.

Results for the analytical models and FEA (Abaqus) were obtained for selected cases of straight fibre and VAT laminates. For practical configurations the capturing of transverse shear and torsional restraint effects where found to be necessary in obtaining accurate results. The beam stiffener model local buckling error was less than 3% compared to FEA whilst the plate stiffener model error was less than 2% compared to FEA. The plate model enabled local modes originating from stiffener blade buckling to be captured and allows the possibility of VAT blade stiffener laminates to be explored. However, the plate stiffener model has addtional shape functions and consequently an increased computational cost compared to the beam model. Global buckling analysis was performed using the beam stiffener model with error for the analytical model ranging from 3% to 10% compared to FEA.

The developed analytical model provides an accurate alternative to the computationally expensive FEA and is therefore suitable for design and optimisation of stiffened VAT panels. Future work will be undertaken to include the stiffener foot in the analysis, improve the estimation of the location of the stiffener local neutral axis, trial a wider range of cases and boundary conditions and perform optimisation studies on the buckling of stiffened VAT panels.

Acknowledgements

The authors gratefully acknowledge the support of the EPSRC under its ACCIS Doctoral Training Centre Grant, EP/G036772/1.

References

[1] Kassapoglou C. Design and analysis of composite structures: with applications to aerospace structures. Wiley; 2010.

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

[3] 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(3):239-61.

[4] Lopes CS, Gurdal Z, Camanho PP. Tailoring for strength of composite steered-fibre panels with cutouts. Compos Part A: Appl Sci Manuf 2010;41(12):1760-7.

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

[6] Weaver PM, Potter KD, Hazra K, Saverymuthapulle MAR, Hawthorne MT. Buckling of variable angle tow plates from concept to experiment. In: AIAA/ ASME/ASCE/AHS/ASC 50th structures, structural dynamics, and materials conference; 2009.

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

[8] Setoodeh S, Abdalla MM, IJsselmuiden ST, Gurdal Z. Design of variable-stiffness composite panels for maximum buckling load. Compos Struct 2009;87(1):109-17.

[9] IJsselmuiden ST, Abdalla MM, Gurdal Z. Optimization of variable-stiffness panels for maximum buckling load using lamination parameters. AIAA J 2010;48(1):134-43.

[10] 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.

[11] Niu MCY. Airframe structural design. Conmilit Press LTD.; 1988.

[12] Yap JWH, Scott ML, Thomson RS, Hachenberg D. The analysis of skin-to-stiffener debonding in composite aerospace structures. Compos Struct 2002;57(1-4):425-35.

[13] Timoshenko S. Uber die stabilitat at versteifter platten. Eisenbau 1921;12:147-63.

[14] Weaver PM. Approximate analysis for buckling of compression loaded long rectangular plates with flexural/twist anisotropy. Proc Roy Soc A: Math Phys Eng Sci 2006;462(2065):59-73.

[15] Bleich F. Buckling strength of metal structures. McGraw-Hill; 1952.

[16] Qiao P, Davalos JF, Wang J. Local buckling of composite FRP shapes by discrete plate analysis. J Struct Eng 2001;127(3):245-55.

[17] Kollar LP. Local buckling of fiber reinforced plastic composite structural members with open and closed cross sections. J Struct Eng 2003;129(11):1503-13.

[18] Paik JK, Thayamballi AK. Buckling strength of steel plating with elastically restrained edges. Thin-Walled Struct 2000;37(1):27-55.

[19] Bisagni C, Vescovini R. Analytical formulation for local buckling and post-buckling analysis of stiffened laminated panels. Thin-Walled Struct 2009;47(3):318-34.

[20] Mittelstedt C, Beerhorst M. Closed-form buckling analysis of compressively loaded composite plates braced by omega-stringers. Compos Struct 2009;88(3):424-35.

[21] Williams JK, Stein M. Buckling behavior and structural efficiency of open-section stiffened composite compression panels. AIAA J 1976;14(11): 1618-26.

[22] Falzon BG, Steven GP. Buckling mode transition in hat-stiffened composite panels loaded in uniaxial compression. Compos Struct 1997;37(96):253-67.

[23] Mittelstedt C. Closed-form analysis of the buckling loads of uniaxially loaded blade-stringer-stiffened composite plates considering periodic boundary conditions. Thin-Walled Struct 2007;45(4):371-82.

[24] Vescovini R, Bisagni C. Buckling analysis and optimization of stiffened composite flat and curved panels. AIAA J 2012;50(4):904-15.

[25] Vescovini R. Analytical formulation for buckling and post-buckling analysis and optimization of composite stiffened panels. Ph.D. thesis. Politecnico di Milano, Italy.

[26] Seide P, Stein M. Compressive buckling of simply supported plates with longitudinal stiffeners. NACATN-1825; 1949.

[27] Block DL, Card MF, Mikulas MM. Buckling of eccentrically stiffened orthotropic cylinders. NASA TN D-2960; 1965.

[28] Steen E. Elastic buckling and postbuckling of eccentrically stiffened plates. Int J Solids Struct 1989;25(7):751-68.

[29] Nagendra S, Gurdal Z, Haftka RT, Starnes JH. Buckling and failure characteristics of compression-loaded stiffened composite panels with a hole. Compos Struct 1994;28(1):1-17.

[30] Nagendra S, Jestin D, Gurdal Z, Haftka RT, Watson LT. Improved genetic algorithm for the design of stiffened composite panels. Comput Struct 1996;58(3):543-55.

[31] Bushnell D. PANDA2 - program for minimum weight design of stiffened, composite, locally buckled panels. Comput Struct 1987;25(4):469-605.

[32] Butler R, Williams FW. Optimum design features of VICONOPT, an exact buckling program for prismatic assemblies of anisotropic plates. AIAA J 1990:1068.

[33] 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.

[34] Reddy JN. Energy and variational methods in applied mechanics. John Wiley and Sons; 1984.

[35] Timoshenko S. Theory of elasticity. McGraw-Hill; 1934.

[36] Reissner E. On the theory of bending of elastic plates. J Math Phys 1944;23:184-91.

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

[38] Libove C, Batdorf SB. A general small-deflection theory for flat sandwich plates. NACA TN-1528; 1948.

[39] Dawe DJ, Roufaeil OL. Rayleigh-Ritz vibration analysis of Mindlin plates. J Sound Vib 1980;69(3):345-59.

[40] Ko WL, Jackson RH. Combined compressive and shear buckling analysis of hypersonic aircraft structural sandwich panels. NASA TM-4290; 1991.

[41] 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(4):567-90.

[42] Herencia JE, Weaver PM, Friswell MI. Optimization of long anisotropic laminated fiber composite panels with T-shaped stiffeners. AIAA J 2007;45(10):2497-509.

[43] Wu Z, Raju G, Weaver PM. Buckling of VAT plates using energy methods. In: AIAA/ASME/ASCE/AHS/ASC 53rd structures, structural dynamics, and materials conference; 2012.

[44] Washizu K. Variational methods in elasticity and plasticity. Pergamon Press; 1975.

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

[46] Jaunky N, Knight NF, Ambur DR. Buckling of arbitrary quadrilateral anisotropic plates. AIAA J 1995;33(5):938-44.

[47] Smith ST, Bradford MA, Oehlers DJ. Numerical convergence of simple and orthogonal polynomials for the unilateral plate buckling problem using the Rayleigh-Ritz method. Int J Numer Methods Eng 1999;44(11):1685-707.

[48] Timoshenko S, Gere JM. Theory of elastic stability. McGraw-Hill; 1961.

[49] Seide P. The effect of longitudinal stiffeners located on one side of a plate on the compressive buckling stress of the plate-stiffener combination. NACA TN-2873; 1953.

[50] Nemeth MP. A treatise on equivalent-plate stiffnesses for stiffened laminated-composite plates and plate-like lattices. NASA TP-2011-216882; 2011.

[51] Kollâr LP, Springer GS. Mechanics of composite structures. Cambridge University Press; 2003.