Two-step procedure for fast post-buckling analysis of composite stiffened panels

Computers and Structures 128 (2013) 38–47

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Two-step procedure for fast post-buckling analysis of composite stiffened panels Riccardo Vescovini, Chiara Bisagni ⇑ Dipartimento di Ingegneria Aerospaziale, Politecnico di Milano, Via La Masa 34, 20156 Milano, Italy

a r t i c l e

i n f o

Article history: Received 26 July 2012 Accepted 11 June 2013

Keywords: Analytical formulation Stiffened panels Composite materials Post-buckling

a b s t r a c t The paper presents an analytical formulation for the post-buckling analysis of composite aeronautical panels with omega stiffeners loaded in compression and shear. The formulation relies on an energy principle and the method of Ritz. In the first step, the panel is an assembly of plate elements, and the buckling analysis is performed. In the second step, the panel is an elastically restrained skin, and the post-buckling behaviour is studied. The comparisons with finite element analyses and experimental results from the literature reveal the ability of the formulation to assess the post-buckling response. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The European consortium MAAXIMUS (More Affordable Aircraft structure through eXtended, Integrated, & Mature nUmerical Sizing) [1], consisting of 58 industries and universities, works to demonstrate the fast development of a highly-optimized composite airframe. The project is divided into a physical platform, focused on the development of composite technologies for low weight aircraft, and a virtual platform aiming to reduce the time required for the identification of the best structural solutions. This paper presents part of the activity performed by Politecnico di Milano during MAAXIMUS in the context of the virtual platform. The work focuses on the development of a fast computational method, which can be used to assess the post-buckling response of composite stiffened panels during the early design stages. The fast tool aims to improve the computational efficiency of the current design procedure, which mainly relies on expensive finite element analyses. Not only the tool can be used for a fast optimization of composite fuselages, but also to account for more degrees of freedom compared to today standard practice. Indeed, the amount of design solutions to be explored can be increased, and consequently the efficiency of the final design can be improved. For this reason, the availability of a fast method represents a crucial aspect to move from a more costly and empirical optimization to a faster structural design optimization loop. Among the fast design methods in the literature, analytical and semi-analytical formulations represent an attractive strategy thanks to their computational efficiency. ⇑ Corresponding author. Tel.: +39 0223998390. E-mail address: [email protected] (C. Bisagni). 0045-7949/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruc.2013.06.002

The simplest approach to study a stiffened panel is to consider the portion of the skin between the stiffeners, assuming the edges as simply supported or clamped. For these cases, closed-form solutions are available both for isotropic [2,3] and composite materials [4–6]. A more refined level of approximation consists in representing the stiffener as an elastic restraint, such as in the analytical formulations of Rhodes and Harvey [7] and van der Neut [8] for isotropic panels and by Mittelstedt and Beerhorst [9] for composite panels subjected to compression loadings. The post-buckling of isotropic stiffened panels is studied with a semi-analytical procedure where the stiffener is modeled as a beam by Brubak and Hellesland [10,11] and Mijukovic et al. [12]. Two of the few analytical works focusing on the post-buckling of elastically restrained composite panels are those ones presented by Chai et al. [13] and by Boay and Wah [14]. The methods of analysis are based on the semi-energy approach in combination with von Kármán large deflection equations. The post-buckling of panels elastically restrained by torsion bars and torsion springs was studied by Bisagni and Vescovini for the analysis of composite stiffened panels, developing a semianalytical procedure [15] and closed-form solutions [16]. An even more refined modeling strategy consists in describing the stiffened panel as an assembly of plate elements. The main advantage is the ability of the formulation to capture local stiffener instabilities, providing at the same time an accurate representation of the skin/stiffener interaction. The drawback of the plate assembly approach is represented by a higher number of degrees of freedom, which turns in a reduction of the computational effectiveness of the analytical procedure. A plate assembly approach is presented by Byklum and Amdahl [3], who developed a computational model for local post-buckling analysis of stiffened isotropic panels. The model considers J-stiff-

39

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

ened panels where the web is modeled as a two dimensional plate element. Buermann et al. [17] extended the approach to the case of curved isotropic stiffened panels, including in the model also frames and skin doublers by means of a beam representation. Both the formulations of Refs. [3,17] are limited to isotropic materials. Despite the relatively large amount of analytical and semi-analytical methods, quite a few of them allow the study of flat and curved composite panels under combined loading conditions. The present work discusses the development of a fast tool for the analysis of the post-buckling response of stiffened panels. The tool offers the advantage of considering a wide set of configurations, including composite and isotropic materials, flat and curved panels, as well as loading conditions of compression and shear. The formulation tries to combine the advantages of the plate assembly representation with the ones of the elastic restraint approach. The use of the fast tool is discussed for the assessment of the buckling and post-buckling response of omega stiffened panels. The analysis procedure is illustrated by describing the input data phase. The results are provided in terms of force–displacement curves, out of plane displacements and maximum stress failure index. The accuracy of the results is demonstrated by comparison with finite element analyses and with experimental data taken from the literature.

2. Description of the formulation The analytical formulation is developed for the buckling and post-buckling analysis of composite flat and curved panels stiffened with omega stringers. The loading conditions are compression and shear, that can be applied separately or in combination. The formulation is divided into two steps, as shown in Fig. 1, and the problem is formulated with an energy approach together with the method of Ritz. In the first step the buckling analysis is performed referring to a structural model based on a plate discretization of the stiffened panel. The approach guarantees a refined description of the panel, and both the skin and the stiffener are explicitly represented as plate elements. The number of degrees of freedom required by the plate representation can be easy handled in the context of a linear analysis, where the buckling load is obtained from the solution of an eigenvalue problem. When assessing the post-buckling response, the governing equations are nonlinear, and it is important to guarantee the computational efficiency through a proper reduction of the number of unknowns. For this reason, the second step is based on a simpler structural model, where the stiffened panel is modeled by considering the portion of skin between two stiffeners. The elasticity of the restraint provided by the stiffeners at the skin longitudinal edges is accounted for by means of torsion springs along the longitudinal edges, whose stiffness is determined from

Fig. 1. Overview of the two-step post-buckling procedure.

the buckling analysis. The choice of the elastically restrained skin model is justified by the fact that the first structural instability of aeronautical panels usually regards the skin. For this reason, the initial post-buckling response can be studied by considering the nonlinear behaviour of the skin, while the response of the stiffener can be reasonably assumed linear. The laminates composing the stiffened panel have an arbitrary number of layers. They are modeled as thin plate elements referring to Kirchhoff hypothesis together with classical lamination theory [18,19]. The effect of transverse shear deformation is not taken into account, and the formulation can be applied to study the behaviour of panels for which the buckling halfwave length is larger compared to the panel thickness. The constitutive equation of each laminate is:

8 9 2 Nx > A11 > > > > > 6 > > > > N A12 y > > > > 6 > = 6 6 A16 xy ¼6 6B > M x > > > 6 > > 6 11 > > > > My > > 4 B12 > > > > : ; M xy B16

A12

A16

B11

B12

A22

A26

B12

B22

A26

A66

B16

B26

B12

B16

D11

D12

B22

B26

D12

D22

B26

B66

D16

D26

9 38 B16 > x > > > > > y > > > B26 7 > > 7> > > > 7> < = 7 c B66 xy 7 7 > w;xx > D16 7> > > > > 7> > > w;yy > > D26 5> > > > : ; 2w;xy D66

ð1Þ

where Nx,Ny,Nxy are the membrane forces per unit length, Mx,My,Mxy the bending and twisting moments per unit length, x,y,cxy the membrane strains, w the out of plane displacement and the derivatives of w represent the curvatures. The comma followed by the coordinate denotes the differentiation with respect to that coordinate. The coefficients Aik and Dik define the in plane and out of plane stiffness of the laminate, respectively, and the coefficients Bik determine the coupling between in plane and out-of-plane behaviour. Each laminate composing the panel can be made of an arbitrary number of layers under the assumptions that:  the laminate is symmetric with respect to the midplane, i.e., Bik = 0  the coupling between extension and shear is null, i.e., A16 = A26 = 0 The formulation accounts for the terms D16 and D26 as typical aeronautical panels can exhibit coupling between bending and twisting. 3. Linear buckling analysis In the first step, the buckling analysis is performed referring to a plate representation of the panel cross section. The approach analyzes a multistiffened panel by studying a representative unit composed of two stiffeners, one bay and two half-bays. An example is reported in Fig. 2(a), where the plate representation of an omega stiffened panel is shown with eleven plate elements. Five elements model the skin, including the portion of skin under the stiffeners, and six elements model the stiffeners. A generic plate element

Fig. 2. Plate assembly representation of the stiffened panel: (a) subdivision in plate elements, (b) generic element dimensions and reference system.

40

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

composing the panel cross-section is reported in Fig. 2(b), where a denotes the length common to all the elements, while b is the width of the element. A Cartesian coordinate system is taken on the element midsurface on lower right corner, with the x-axis directed along the longitudinal direction, the y-axis along the transverse direction, and the z-axis along the normal direction, according to the right-hand rule. Similar plate element representations can be used for the analysis of panels with different stiffener geometries, such as blade, J and T stiffeners, even if the present work focuses on omega stiffened panels. The formulation is developed starting from the total potential energy P of the structure, which is:

P¼

Np X

Pi þ U p

ð2Þ

i

where Np is the number of elements that compose the representative unit, and Up is a penalty term which is added to the functional to enforce compatibility of the rotations along the edges of adjacent plates. The potential energy of each plate element is obtained as the sum of three contributions:

Pi ¼ U b þ U m þ V load

ð3Þ

where Ub is the strain energy for bending, Um is the membrane energy and Vload is the potential of the external loads. The membrane energy term is introduced only for curved elements due to the coupling between in plane and out of plane behaviour. It is neglected in case of flat plates, as the in plane and out of plane response are uncoupled. Referring to Kirchhoff thin plate theory and classical lamination theory, the bending energy of the i-th element is expressed as:

Z  1 D11 w2;xx þ 2D12 w;xx w;yy þ D22 w2;yy þ 4D66 w2;xy 2 A  þ 4D16 w;xx w;xy þ 4D26 w;yy w;xy dA

Ub ¼

ð4Þ

where A is the area of the plate element, w and Dik are the out of plane displacement and the bending stiffness of the plate element. The membrane energy is written in terms of Airy stress function F, and is defined as:

Um ¼

1 2

Z   a11 F 2;yy þ a22 F 2;xx þ 2a12 F ;xx F ;yy þ a66 F 2;xy dA

ð5Þ

A

where aik is the in plane compliance of the laminate. The potential of the external loads is obtained as the sum of the contributions of compression and shear loads:

V load ¼

1 kNx 2

Z A

w2;x dA þ kN xy

Z

w;x w;y dA

ð6Þ

A

where k is the multiplier of the pre-buckling condition. Without loss of generality, the axial resultant Nx is determined by assuming that the plate elements undergo the same shortening. Indeed, in a multi-stringer aeronautical structure, the axial strain compatibility condition is usually enforced by a rib or a frame. As the shear load carrying capacity of a stiffened panel is mainly demanded to the skin, the shear force Nxy is introduced only on the skin elements. The same approach is found in [20]. The essential boundary conditions of the problem are expressed as:

8 on x ¼ 0; a > : Du ¼ 0 along the edges of adjacent elements ij

ð7Þ

where Duij is the difference of rotation between the two elements i and j expressed in function of the out of plane displacement.

The first two conditions express null out of plane displacement along the four edges of the plate element, and are due to the local buckling hypothesis. The third condition regards the continuity of rotations between adjacent plate elements. The shape functions used to describe the out of plane displacement over these elements is a double trigonometric series of sine terms:

w¼

M;N X mpx npy qmn sin sin a b m;n

ð8Þ

where qmn are the unknown amplitudes, while M and N are the number of shape functions along the longitudinal and transverse direction. It can be observed that the first two essential conditions of Eq. (7) are identically satisfied, whereas the third one is imposed with a penalty approach as described in the next paragraph. The unknown amplitudes of the generic i-th element are collected into the column vector qi, and the minimum potential energy principle is applied:

@ Pi @ðU b þ U m Þ @V load ¼ þ ¼0 @qi @qi @qi

ð9Þ

The trigonometric expansions of the problem unknowns are then substituted into Eq. (9) and, for each element composing the section, the analytical integration of the terms is carried out. The two terms of Eq. (9) are re-written as:

@ðU b þ U m Þ ~ ¼ Ki qi @qi

@V load ¼ kGi qi @qi

ð10Þ

~ i and Gi are the stiffness and the loading stiffness matrices where K of i-th element. The degrees of freedom of the structure are obtained by collecting the unknowns qi into a global vector of unknowns q, defined as:

q¼

n

qT1

qT2

   qTNp

oT

ð11Þ

The stiffness and the loading stiffness matrices of the single plate ~ i and Gi are assembled according to the definition of elements K ~ and G. Eq. (11) to obtain K The last step of the procedure regards the enforcement of the compatibility condition of the rotations along the longitudinal edges of adjacent plate elements. The condition is imposed by introducing penalty terms, which are torsion springs along the common edges of adjacent plates. Their expression is:

U p;ij ¼

1 kt 2

Z 0

a

Du2ij dx

ð12Þ

where kt is the stiffness of the spring. An additional term, K, is obtained from the differentiation of Eq. (12) with respect to the unknown amplitudes q. The term is added to the stiffness matrix, obtaining:

~ þK K¼K

ð13Þ

The assembled buckling equations are finally written as:

ðK þ kGÞq ¼ 0

ð14Þ

where K and G are the stiffness and the loading stiffness matrices, while q is the vector collecting the unknown amplitudes describing the out of plane deflections of the plate elements. The eigenvalue problem of Eq. (14) is solved numerically and, to improve the efficiency of the procedure, it is decomposed into a set of smaller subproblems accounting for the symmetry and antisymmetry of the buckling modes.

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

4. Post-buckling analysis The structural model used in the second step is reported in Fig. 3, and is representative of the skin of the stiffened panel. The length and the width are denoted with a and b, respectively. The transverse edges are simply supported, while the stiffeners are accounted for by introducing a torsion spring of stiffness kt along the longitudinal edges. The stiffeners are also characterized by an axial stiffness, which is denoted as EA. It is observed that, for a closedsection stiffener, the axial stiffness EA accounts not only for the contribution of the stiffener itself, but also for the portion of skin under the stiffener. To represent the presence of the surrounding structure, whose effect is to support the in plane motion of the panel, the longitudinal edges are free to move along the transverse direction, but are forced to remain straight. A Cartesian coordinate system is taken over the skin midsurface with the x-axis directed along the longitudinal direction, and the yaxis along the transverse direction. It is observed that the elastically restrained model, compared to the plate assembly one, offers the advantage of limiting the computational time. Indeed, it allows a reduction of the number of the degrees of freedom, and does not require to iterate at each step of the incremental procedure to determine the fraction of load carried by each plate element [3,17]. A further complication in the use of the plate assembly model would be represented by the need for finding a different solution of the compatibility equation for each different expression of the out of plane displacement. 4.1. Evaluation of the elastic restraint The stiffness of the elastic restraint kt is determined by imposing that the buckling load of the elastically restrained panel is equal to the buckling load of the plate assembly model. The approach assumes that the stiffeners restraint the rotation of the skin edges, therefore representing an intermediate condition between the simply-supported edges and the clamped edges [3], and carry axial loads according to their axial stiffness EA. The buckling eigenvalue problem is derived referring to the minimum potential energy principle. In this case, the total potential energy is expressed as:

P ¼ U b þ U m þ V load þ U k

ð15Þ

where the term Uk is added to the terms previously considered in the linear analysis. It represents the energy stored in the torsion spring, and is given by:

1 U k ¼ kt 2

Z 0

a

 

w2;y  dx y¼0

1 þ kt 2

Z

a 0

 

w2;y 

dx y¼b

41

After applying the minimum potential energy principle, the problem is written in the form of Eq. (14), where the matrix K is now function of the unknown stiffness kt. The panel buckling load is obtained as the sum of two contributions. The first one is the load acting on the skin, which is obtained by integrating the skin buckling force per unit length along the panel width. The second term is the load acting on the stiffeners, and can be derived assuming that both the skin and the stiffeners undergo the same displacement. The total load Ptot is:

Ptot ¼

Z 0

b

   Nx  

dy þ x¼0

EA DU a

ð17Þ

where DU is the buckling shortening. An iterative procedure is implemented, where the value of kt is progressively increased starting from kt = 0, until the difference between the buckling loads of the elastically restrained panel and the plate assembly model is below a certain tolerance value. The obtained stiffness kt is then used to perform the post-buckling analysis. 4.2. Derivation of the governing equations The formulation for the post-buckling analysis is organized as shown in the scheme of Fig. 4. The out of plane displacement is described with an arbitrary number of trigonometric shape functions. The bending energy, the strain energy stored in the elastic restraint, as well as the potential of the external loads are directly calculated as function of the out of plane displacement amplitudes. The membrane energy is written in terms of Airy stress function, whose expression is obtained imposing the exact solution of the nonlinear compatibility equation. It is worth observing that the membrane energy is introduced not only in the case of curved panels, but also for flat panels, as the in plane and out of plane behaviour are coupled in the nonlinear post-buckling range. The nonlinear governing equations are derived imposing the stationarity of the total potential energy and referring to the method of Ritz. The numerical solution is performed with an arc-length procedure which allows to capture eventual snap-back and snapthrough responses. The nonlinear compatibility equation for a laminated panel with no shear/extension coupling and characterized by a radius of curvature R is:

a11 F ;yyyy þ ð2a12 þ a66 ÞF ;xxyy þ a22 F xxxx ¼ w2;xy  w;xx w;yy þ 2w0;xy w;xy  w0;xx w;yy  w0;yy w;xx

ð16Þ 

ðw þ w0 Þ;xx R

ð18Þ

where w0 are the initial imperfections. They are introduced to account for small deviations from the nominal configuration. Furthermore, the presence of small imperfections facilitate the numerical solution of the problem, avoiding problems related to the convergence. The fulfillment of the compatibility requirement of Eq. (18) allows to impose the equilibrium of the solution referring to the minimum potential energy principle, which is written in terms of out of plane displacement and Airy stress function. In particular, the total potential energy is obtained as:

P ¼ Um þ Ub þ Uk þ V c þ V s

Fig. 3. Elastically restrained skin.

ð19Þ

where the terms Um and Ub are the membrane and the bending energy, whose expression is reported in Eqs. (4) and (5). It is worth observing that the strain energy expression of Eq. (19) is simplified, as the contribution due to the stiffeners axial stiffness is neglected.

42

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

Fig. 4. Formulation for the nonlinear analysis.

Indeed, it was observed by Brubak and Hellesland [10] that the introduction of this term determines a significant increase of the computational time, up to a factor of about 30, while playing a minor role in case of skin local buckling conditions. The same approach is found in the formulation of Paik and Lee [21]. In any case, the stringer axial stiffness EA is accounted for to compute the total load carried by the panel according to Eq. (17). The last two terms of Eq. (19) are the potential of the applied loads of compression and shear. They are written as:

 1 V c ¼ Nx b a11 F ;yy þ a12 F ;xx  w2;x  w;x w0;x dx 2 0 Z   V s ¼ Nxy w;x w;y þ w;x w0;y þ w;y w0;x dA Z

a

w0 ¼

m;n

2



r;s;t;u¼1 brstu  4 a11 nb þ ð2a12

 qrs qtu þ qrs q0tu þ qtu q0rs  2  2  4 i þ a66 Þ ma nb þ a22 ma

ð24Þ

where the array of scalars brstu is defined as:

brstu ¼ rstu  r 2 u2 brstu ¼ rstu þ r 2 u2



if r þ t ¼ m and s þ u ¼ n



if jr  tj ¼ m and js  uj ¼ n if jr  tj ¼ m and s þ u ¼ n

ð25Þ

if r þ t ¼ m and js  uj ¼ n

and:

ð21Þ

~f mn ¼ b ~mn q þ q0  mn mn

ð26Þ

with:

mpx npy sin sin a b

ð22Þ

where q0mn are the known amplitudes of w0. The trigonometric representation for the Airy stress function is derived from the substitution of the out of plane displacement into Eq. (18). The expression of the Airy stress function that identically satisfies the compatibility equation is found as [3,17]:

1 F ¼  Nx y2 þ Nxy xy þ 2 þ

4a2 b

h

ð20Þ

The variational principle expressed by the functional of Eq. (19) is solved applying the method of Ritz. The out of plane displacement is represented using the trigonometric expansion of Eq. (8). In this case, the essential boundary conditions regard the null displacement along the four edges of the panel, and coincide with the first two conditions of Eq. (7). The use of the shape functions of Eq. (8) ensures the fulfillment of the two essential conditions. The same shape functions of Eq. (8) are used to describe the initial imperfections:

q0mn

PM;N;M;N

1



A

M;N X

fmn ¼

2M;2N X

fmn cos

m;n

M;N X ~f mn sin mpx sin npy a b m;n

mpx npy cos a b ð23Þ

where the first two terms describe a state of uniform axial and shear force. The third term is responsible for the force redistribution due to the large deflections, and the fourth term describes the distribution of the in plane forces related to the coupling between in plane and out of plane behaviour due to the panel curvature. The Airy stress function amplitudes fmn and ~f mn do not introduce further degrees of freedom in the problem, as they are expressed as function of the amplitudes qmn from the solution of the compatibility equation. In particular, after substitution of Eqs. 8, 22 and 23 into Eq. (18), it is obtained:

 2 1 ~mn ¼ m h i b a Rp2 n4 a11 þ mn2 ð2a12 þ a66 Þ þ m4 a22 b

ab

ð27Þ

a

The nonlinear governing equations are derived applying the method of Ritz. The shape functions of Eqs. 8, 22 and 23 are substituted into the expression of Eq. (19), and the first variation of the total potential energy is set to zero:

@P ¼0 @q

ð28Þ

which are M  N equations, where q is the column vector collecting the unknown amplitudes qmn. 4.3. Numerical solution The nonlinear Eqs. (28) can be solved with a Newton–Raphson incrementation technique. However, the post-buckling response of the panels is sometimes characterized by phenomena such as snap-backs or snap-throughs, making necessary the use of more appropriate solution techniques. The strategy here implemented refers to the work of Steen [22], also implemented by Byklum [3] and Buermann et al. [17] in the context of semi-analytical formulations for the post-buckling of isotropic panels. The method is an arc-length procedure, implemented with a perturbation approach. Following the approach of Ref. [22], the rate parameter g, which is a pseudo-time, is introduced. The loading conditions are defined through the introduction of the load parameter K, which is a scalar that becomes part of the solution. The governing equations are then written as function of the amplitudes q and the load parameter K, and are derived minimizing the total energy P of Eq. (19) in the rate form:

43

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

@2P @ 2 P @q @2P @K ¼ þ ¼0 2 @ g@q @q @ g @q@ K @ g

ð29Þ

which is re-written in a more compact notation as:

Kq;g þ GK;g ¼ 0

ð30Þ

The introduction of the load parameter K requires one further equation to solve the set of Eqs. (30). The constraint equation is then introduced, and the resulting system of nonlinear equation results:

(

Kq;g þ GK;g ¼ 0

K2;g þ qT q ¼ 1

ð31Þ

The nonlinear Eqs. (31) are finally solved with a perturbation approach [22], where the unknowns q and K are represented with a first order Taylor expansion. 5. Implementation of the fast tool The fast tool is implemented in a computer code working in Matlab [23] environment. It results suited to perform quick analyses, to study different configurations and to realize parametric studies. Thanks to the user-friendly graphical interface, the program allows for a simple introduction of the input data. The length and the width of the panel, as well as the dimensions of the stiffener, are firstly introduced as inputs. Then, the program requires the material properties in the form of elastic moduli and ply strengths. The values can be set manually or can be loaded from a database. Furthermore, the stacking sequences of the skin and the stiffener are introduced. The fast tool offers the possibility of choosing between two solution procedures. The first one is the linear buckling analysis, and relies on the plate assembly model, which is the first step of the procedure. The second possibility is the post-buckling analysis, and refers to the complete two-step procedure. The input is completed by setting up the analysis parameters. The program requires to define the pre-buckling load and the number of shape functions used in the plate assembly model. In case of post-buckling analysis, further inputs are expected. They include the number of shape functions for the post-buckling procedure, the total applied load, the number of buckling eigenvalues used as initial imperfections and their amplitudes. Different outputs can be requested on the basis of the chosen type of analysis. In case of buckling analysis, the user can require a summary of the results, including the buckling load, the corresponding mode and the linear stiffness. In case of post-buckling analysis, the available outputs are the force–displacement curve, the deformed shape, and the contour of three different failure criteria (maximum stress, maximum strain, Tsai-Hill). 6. Analysis of an omega stiffened panel The fast tool is applied to study a curved omega stiffened panel subjected to loading conditions of combined compression and shear. The load history is divided in two phases: 1. pure shear load per unit length of 62 N/mm 2. compression load with an axial displacement of 2.10 mm, and shear load per unit length kept constant at 62 N/mm

one bay, and two half bays. The height of the stiffener is 24 mm, the crown top is 27 mm long and the web angle measured from the normal to the skin is equal to 10°. The skin consists of quasiisotropic laminate of 12 plies with a stacking sequence of [45°/ 45°/0°/90°/ 45°/45°]s, for a total thickness of 1.5 mm. The stringers are composed of a laminate of 9 plies with stacking sequence of [45°/0°/ 45°/0°/90°/0°/ 45°/0°/45°] and a total thickness of 1.125 mm. The results computed with the fast tool are compared with analyses performed using the commercial finite element code Abaqus [24]. The finite element model of the panel is reported in Fig. 5. The skin and the stiffeners are modeled with four-noded S4R shell elements, with a mesh of about 5  5 mm chosen on the basis of a preliminary convergence analysis. The loaded edges are simply supported, whereas a constraint equation is applied to the longitudinal edges to enforce periodic conditions in terms of out of plane displacements and rotations around the longitudinal axis. Two analyses are performed, eigenvalue and nonlinear static. The input data are introduced by means of the graphical user interface. It is interesting to highlight the choice of the analysis parameters, which are the number of shape functions and the initial imperfections. They can be defined case by case on the basis of a convergence analysis. In the first step, the out of plane deflections are represented using 10  10 functions for the skin, and 10  7 for the stiffener elements, which are a good compromise between computational efficiency and accuracy of the results. Considering that the plate assembly model exploits the symmetries and the anti-symmetries of the buckling modes, the total number of degrees of freedom for the buckling analysis is equal to 235. In the second step, the skin deflection is represented with a number of 12  12 functions. The resulting problem is given by 145 equations, consisting in 144 equilibrium equations and 1 constraint equation requested by the arc-length procedure. Compared to the first step, the number of shape functions is increased. This choice is motivated by the need to calculate the stresses, which are related to the second derivatives of the displacements, and consequently characterized by a slower convergence rate. Regarding the introduction of the initial imperfections, it is important to avoid the convergence of the solution to the unbuckled configuration. The initial imperfections can be determined through a Fourier analysis of experimental data, if available. Otherwise it is possible to use a linear superposition of buckling modes, which is the strategy here adopted. In the present study, the initial imperfections are taken equal to the first buckling mode, with an amplitude to thickness ratio of 5%. The outputs are then defined. In particular, the results of the first step analysis are required in terms of buckling load and buckling mode, while the results of the post-buckling analysis are defined as force displacement curve, deformed surface at different load levels and contour of the maximum stress index. The buckling load obtained from the first step of the procedure is characterized by a compressive force per unit length equal to 59.2 N/mm, and a shear force per unit length of 62 N/mm. The buckling mode shape reveals 6 skewed halfwaves on the skin, whereas the plate elements composing the omega stiffener do

Table 1 IM7/8825 ply mechanical properties. Material properties

The panel is made from carbon/epoxy unidirectional IM7/8552, whose properties are summarized in Table 1. It has dimensions 300  536 mm, and a radius of curvature equal to 2075 mm. A sketch is reported in Fig. 5. The panel, representative of a larger structure with several stringers, is composed of two stiffeners,

Young’s modulus E11 [MPa] Young’s modulus E22 [MPa] Shear modulus G12 [MPa] Poisson’s ratio m12 Ply thickness [mm]

150000 9000 5300 0.32 0.125

44

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

Fig. 5. Curved omega stiffened panel: (a) finite element model, (b) stiffener cross section.

2.5 2 1.5 1 0.5 0 0

Analytical ABAQUS 0.5

1

1.5

2

2.5

3

Fig. 6. Force–displacement curve during the compressive loading phase.

not undergo significant out of plane deflections. The stiffness of the elastic restraint, which is computed at the beginning of the postbuckling analysis procedure, is equal to 2200 N. The force–displacement curve is reported for the compressive loading phase in Fig. 6, where the results are normalized with

respect to buckling conditions under compression load. A drop of stiffness is observed when u/ucr is below the unity, as the applied shear load has the effect of reducing the buckling load. The drop of stiffness after buckling, expressed as ratio between pre- and post-buckling stiffness, is approximately 0.80. A slight difference is observed from the comparison with the numerical result, which is equal to approximately 0.76. In this case, the analytical model is stiffer than the numerical one, and higher loads are obtained for the same value of applied displacement. At the maximum value of imposed displacement, the analytical load is 3% higher than the numerical one. The out of plane displacements in the post-buckling field are illustrated in Fig. 7. The deformed configuration is characterized by six halfwaves, whose skewness is due to the combined effect of the bending/twisting coupling and the shear pre-loading. The comparison with the finite element results reveals the ability of the formulation to represent a relatively complex out of plane deflected pattern, as the one displayed in presence of shear. The fast tool predicts a range of displacements in the range between 3.1 mm and 3.1 mm, while the finite element out of plane displacement are in the range between 3.8 mm and 4.2 mm. The formulation is able to provide also an estimate of the maximum stress distribution on the skin in order to identify the most critical areas which may exhibit damage onset.

Fig. 7. Contour of the out of plane displacement: (a) analytical, (b) Abaqus.

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

45

Fig. 8. Contour of the maximum stress failure index: (a) analytical, (b) Abaqus.

Table 2 T300/934 ply mechanical properties. Material properties Young’s modulus E11 [MPa] Young’s modulus E22 [MPa] Shear modulus G12 [MPa] Poisson’s ratio m12 Ply thickness [mm]

124000 11400 5500 0.3 0.145

The stress field is plotted at the last step of the solution. In this example, the maximum stress criterion is computed for each ply of the laminate. The contour is reported for the ply undergoing the highest value of failure index, which is automatically identified by the fast tool. In this case, the procedure identifies the second ply of the skin, which is oriented at  45°, as the most stressed one. The contour is reported in Fig. 8, with a maximum value equal to 0.50 in proximity of the longitudinal edges. The numerical results of Fig. 8 are characterized by a maximum stress failure index of 0.65. The results obtained with the fast tool can be used as a preliminary estimate of the stress level exhibited by the structure. The

level of approximation introduced by the analytical approach is mainly related to the development by means of a mixed formulation, where the unknowns are the out of plane displacement and the Airy stress function. Consequently, the formulation does not handle the in plane displacements, which have a significant influence on the shear response of the structure. It is observed that the loading condition of compression and shear represents the most complex scenario for the present method, and much smaller differences are obtained in the pure compression loading case. In any case, the fast tool allows for a correct assessment of the deformed configuration shape, so providing useful information regarding the areas undergoing the largest out of plane displacement, as well as the areas and the plies subjected to the highest stress level. Moreover, the gain in terms of CPU time is significant, as the analysis required 43 s versus the 250 s required for the finite element approach, using a Core2 Duo 3.00 GHz with 2 GB of RAM. It is important to note that the difference of time between the two approaches would be much greater in the context of preliminary design optimization, where thousands of analyses are performed. It is also observed that the present method does not require geometric modeling and finite element meshes have not to be created.

Fig. 9. Force–displacement curve of the panel reported in Ref. [26].

46

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

Fig. 10. Experimental [25] and analytical configurations in the post-buckling range: (a) and (b) three halfwave configuration, (c) and (d) mode jumping, (e) and (f) five halfwave configuration.

7. Panel with mode jumping The ability of the formulation to capture mode jumpings is assessed by the comparison with an experimental test published in literature [26,25], where an omega stiffened composite panel exhibits a mode transition in the post-buckling field. The panel is stiffened by two omega stiffeners. The distance from the end tabs is 425 mm, whereas the portion of skin between the stiffeners is equal to 186 mm. The configuration is characterized by the presence of two external bays 40 mm long, whose longitudinal edges are free. The height of the omega stiffener is 25 mm, the length of the crown top is 20 mm, while the portion of the skin under the stiffener is 35 mm. Both the skin and the stiffener are layered with the quasi-isotropic lay-up [45°/ 45°/0°/90°]s.

The panel is made from T300/934 unidirectional prepreg, whose elastic properties are summarized in Table 2, but, following the study of Falzon and Stevens [26], the nominal material properties are reduced by 15.66% in order to avoid an overestimation of the linear pre-buckling stiffness. The fast tool is applied to perform the post-buckling analysis of the panel reported in literature using the two-step approach. It is worth noting that the two external bays do not provide any contribution to the total load carried by the structure in the post-buckling field. Indeed, the external longitudinal edges are free, and consequently buckle at relatively small load levels. In any case, it is expected that these two portions of the structure provide stiffness in the linear pre-buckling field. In the first step of the procedure, the trigonometric expansion is arrested to 10  10 terms for the deflection of the skin, and to 10  7 terms for the deflection of the stiffener elements. In the second step, the number of terms describing the out of plane response of the elastically restrained panel is increased to 12  12. The initial imperfections are introduced as a linear superposition of the first four buckling modes, with a maximum amplitude to thickness ratio equal to 5%. The force–displacement response is reported in Fig. 9, where the analytical results are compared to the experimental ones. The fast tool predicts a buckling load of 9 kN, which is very close to the experimental one of 10 kN. The first buckling mainly involves the skin, and the stiffeners do not experience any out of plane deflections in the initial post-buckling field. For this reason, a relatively small reduction of stiffness can be observed in the force–displacement curve just after the first buckling. The analytical and numerical curves are almost coincident in the initial post-buckling range, until an imposed displacement of approximately 1.5 mm. In the deep post-buckling field, the differences between experimental and analytical results become a little larger and the analytical and the experimental loads differ by 14% when the imposed displacement is 3 mm. This discrepancy is due to the inability of the analytical model to account for the nonlinear response of the stiffeners, including local stiffener instabilities, which becomes relevant at high load levels. From the analytical force–displacement curve it can be observed the mode jumping response of the panel. Indeed while for low levels the load carried by the panel increases monotonically, at 57 kN the structure starts to unload, and progressively turns to another buckled configuration. A similar mode-jumping response is observed in the experimental results, even tough it can not be easily detected from the experimental curve. For this reason, it is useful to present the results in terms of deformed configurations. The more relevant phases are shown in Fig. 10, where the experimental MoirF pattern is compared to the analytical contour of the out of plane deflections. In the initial post-buckling range the buckled configuration is characterized by three halfwaves, as shown in Fig. 10(a) and (b). When the load is increased, the midbuckle enlarges and begins to split into two smaller buckles, as revealed by Fig. 10(c) and (d). The transition terminates with a five halfwaves configuration shown in Fig. 10(e) and (f). Referring to Fig. 10, it can be seen that the experimental mode jumping is observed at 66 kN, while the analytical one at 57 kN. In any case, the result is satisfactory if it is considered that the mode jumping is a highly nonlinear phenomenon that can be hardly reproduced with a high level of accuracy even making use of more sophisticated tools, such as the finite element method. Furthermore, the analytical model relies on several assumptions, including the simplified representation of the stiffener and the simply supported loaded edges and does not consider the small imperfections regarding the loading and boundary conditions that always characterize a physical test. It is worth to note the ability of the fast tool

R. Vescovini, C. Bisagni / Computers and Structures 128 (2013) 38–47

to detect the presence of a mode-jumping and to correctly predict the pattern before and after it. 8. Conclusions An analytical formulation has been developed for the postbuckling analysis of composite stiffened panels based on a two step procedure, which is developed referring to an energy principle and the method of Ritz. In the first step the buckling load is determined by modeling the panel as an assembly of plate elements, while in the second step the post-buckling analysis is performed by assuming the panel as an elastically restrained skin, considering so the post-buckling response with a relatively small number of degrees of freedom. The approach offers the advantage of allowing the study of composite panels subjected to combined loading conditions of compression and shear. The comparison with finite element analyses reveals that the fast tool can be confidently used to predict the drop of stiffness due to the buckling of the skin and to trace the force–displacement curve. Agreement between analytical and numerical predictions is observed in terms of post-buckled pattern, as the number of halfwaves and their skew angle are correctly predicted. Furthermore, the tool can be used to obtain an approximate representation of the stress index distribution of the skin, and to identify the ply undergoing the highest stress level. The implementation of an arc-length solution procedure allows to capture mode jumping phenomena, in case they are present, as demonstrated by comparison with experimental results taken from the literature. The reduced amount of time required to perform the analysis suggests the use of the tool within a preliminary design optimization loop, where buckling and post-buckling requirements are introduced as constraints or compose the objective function. The tool can also be inserted as a module for the post-buckling analysis in the context of a multidisciplinary design procedure, or eventually coupled to a finite element code to obtain a hybrid design procedure. Acknowledgements The research leading to these results has been partially funded by the European Commission Seventh Framework Programme FP7/ 2007–2013 under Grant agreement No. 213371, MAAXIMUS (www.maaximus.eu). References [1] Maaximus Project website. . [2] Levy S. Bending of rectangular plates with large deflections. TM 737, NACA, 1942.

47

[3] Byklum E, Amdahl J. A simplified method for elastic large deflection analysis of plates and stiffened panels due to local buckling. Thin-Walled Struct 2002;40(11):925–53. [4] Harris GZ. Buckling and post-buckling of orthotropic laminated plate. In: 16th AIAA/ASME/SAE Structures, Structural Dynamics, and Materials Conference, AIAA-1975-813, May 27–29 1975, Denver, CO. [5] Romeo G, Frulla G. Analytical/experimental behavior of anisotropic rectangular panels under linearly varying combined loads. AIAA J 2001;39(5):932–41. [6] Diaconu CG, Weaver PM. Postbuckling of long unsymmetrically laminated composite plates under axial compression. Int J Solids Struct 2006;43(22– 23):6978–97. [7] Rhodes J, Harvey JM. Plates in uniaxial compression with various support conditions at the unloaded boundaries. Int J Mech Sci 1971;13(9):787–802. [8] van der Neut A. The stiffness in compression of imperfect elastically restrained plate strips at various in-plane boundary conditions. Report LR-245, Delft, 1977. [9] Mittelstedt C, Beerhorst M. Closed-form buckling analysis of compressively loaded composite plates braced by omega-stringers. Compos Struct 2009;88(3):424–35. [10] Brubak L, Hellesland J. Semi-analytical postbuckling and strength analysis of arbitrarily stiffened plates in local and global bending. Thin-Walled Struct 2007;45(6):620–33. [11] Brubak L, Hellesland J. Semi-analytical postbuckling analysis of stiffened imperfect plates with a free or stiffened edge. Comput Struct 2011;89(17– 18):1574–85. [12] Mijukovic O, Coric B, Pavlovic MN. Transverse-stiffener requirements for the post-buckling behaviour of a plate in shear. Thin-Walled Struct 1999;34(1):43–63. [13] Chai GB, Banks WM, Rhodes J. The instability behaviour of laminated panels with elastically rotationally restrained edges. Compos Struct 1991;19(1):41–65. [14] Boay CG, Wah KP. Postbuckling formulation for symmetric laminated panels with various edge support conditions. Mech Adv Mater Struct 2001;8(1):15–28. [15] Bisagni C, Vescovini R. Analytical formulation for local buckling and postbuckling analysis of stiffened laminated panels. Thin-Walled Struct 2009;47(3):318–34. [16] Vescovini R, Bisagni C. Single-mode solution for post-buckling analysis of composite panels with elastic restraints loaded in compression. Compos: Part B 2012;43(3):1247–58. [17] Buermann P, Rolfes R, Tessmer J, Schagerl M. A semi-analytical model for local post-buckling analysis of stringer- and frame-stiffened cylindrical panels. Thin-Walled Struct 2006;44(1):102–14. [18] Hyer MW. Stress analysis of fiber-reinforced composite materials. New York: McGraw-Hill; 1998. [19] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. Boca Raton: CRC Press; 2004. [20] Loughlan J. The buckling performance of composite stiffened panel structures subjected to combined in-plane compression and shear loading. Compos Struct 1994;29(2):197–212. [21] Paik JK, Lee MS. A semi-analytical method for the elastic–plastic large deflection analysis of stiffened panels under combined biaxial compression/ tension, biaxial in-plane bending, edge shear, and lateral pressure loads. ThinWalled Struct 2005;43(3):375–410. [22] Steen E. Application of the perturbation method to plate buckling problems. Technical Report 98-1, University of Oslo, Department of Mathematics, Mechanics Division, 1998. [23] . Matlab. User’s Guide 2011. Mathworks, Inc.; 2011. [24] . ABAQUS, version 6.11. User’s Manual. Providence, RI, USA: SIMULIA World Headquarters; 2011. [25] Falzon BG. The behaviour of damage tolerant hat-stiffened composite panels loaded in uniaxial compression. Compos Part A: Appl Sci Manuf 2001;32(9):1255–62. [26] Falzon BG, Steven GP. Buckling mode transition in hat-stiffened composite panels loaded in uniaxial compression. Compos Struct 1997;37(2):253–67.