Some approaches to buckling analysis of flexurally anisotropic composite plates, subjected to combined in-plane loading

Thin-Walled Structures 98 (2016) 375–383

Contents lists available at ScienceDirect

Thin-Walled Structures journal homepage: www.elsevier.com/locate/tws

Some approaches to buckling analysis of flexurally anisotropic composite plates, subjected to combined in-plane loading Sergey Selyugin AIRBUS Operations GmbH, Kreetslag 10, 21129 Hamburg, Germany

art ic l e i nf o

a b s t r a c t

Article history: Received 2 February 2015 Received in revised form 9 October 2015 Accepted 12 October 2015

Buckling of flexurally anisotropic composite plates under combined loading is considered. The purpose of the paper is to develop approaches for predicting buckling loads for the plates, basing on buckling analysis results for separate loadings. Theoretical aspects of the problem are explored. The derivatives of the buckling interaction curves are analytically derived. The finite predicting formulae are proposed. The Galerkin-type numerical method for combined loading buckling analysis is developed. The method is based on buckling modes of separate loadings. Numerical examples are discussed. Convergence of the method is shown. Influence of the flexural anisotropy on the buckling behavior is demonstrated. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Buckling Plate Flexural anisotropy Combined loading

1. Introduction Composite materials are widely used now in many applications like aerospace, marine, civil structures, etc. In the majority of cases the structural components made of the materials are thin-walled structures like plates, shells, panels, etc. For the components, in particular, for the plates the buckling phenomenon is one of the main driver in designing them. Our object of interest in the current paper is the plates. Talking about plates, their loading is, as a rule, a combination of simple loadings like compression or tension, shear, in-plane bending, etc. First we say several words on a role of bending–twisting coupling in pure shear buckling of composite plates. As it is indicated in, for example, [1,2], the shear buckling eigenvalues of composite plates may be “markedly different” for positive and negative shear loading directions. This fact must be taken into account in analysis of combined loading of composite plates. Many features inherent in plate buckling under combined loading were profoundly explored for the case of isotropic materials. Usually the buckling behavior of a plate is described by the so-called buckling interaction curve in the plane of two load quantities corresponding to every separate loading participating in the combined loading (in case of combination of three loadings there must be a surface in the corresponding loading space). The features of the curves are comprehensively described in many publications and textbooks, for references we indicate the books of Timoshenko and Gere [3], Kollar [4] (Section 2), Alfutov [5], E-mail address: [email protected] http://dx.doi.org/10.1016/j.tws.2015.10.008 0263-8231/& 2015 Elsevier Ltd. All rights reserved.

Perelmuter and Slivker [6]. Among the features one may indicate the Dunkerley formula (an inequality linking buckling eigenvalues of combined and separate loadings) and the Föppl–Papkovich theorem (called as the Papkovich theorem in some books) on convexity of the buckling interaction curve. The convexity of the curve in case of beams under combined loading was also shown in the study of Schäfer [7]. It is worth mentioning that one may find an elegant derivation of the Dunkerley formula for plate buckling in the Section 8 of the book of Collatz [8]. Composite plates retain some features inherent in the isotropic case, but the anisotropy of the former plates leads to new features and peculiarities. The features appear for combined loading cases also. The buckling studies of anisotropic plates under combined loading were started in 1930s in a series of papers devoted to plywood structures. At the period the anisotropic material was widely used in the airplane structures. The proper references to corresponding papers of British, German and Soviet researches may be found in the above mentioned books. In particular, the approach accounting for flexural plate anisotropy and leading to a finite shear–compression buckling interaction formula was developed by Balabukh (indicated in the book of Lekhnitskii [9]). His known formula is, in fact, a parabola, coming through points of precise values in buckling interaction plane, the points corresponding to pure compression and pure shear of both loading directions. Upon introduction of composite materials into practice in 1970ies the composite plate buckling under combined loading was widely explored. The most comprehensive results (including simple interaction formulae) are available for long plates made of specially

376

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

orthotropic materials, i.e. without consideration of flexural anisotropy (see the review in [10]). In particular, in [10] the interaction formulae are derived using regression analysis. The formulae lead to prediction accuracy for combined loading about 5%. The general flexurally anisotropic composite plates are not so widely explored; one may indicate the papers of Nemeth [11–14]. The papers contain a lot of numerical analysis results for such plates with various anisotropy parameters. In [11,12] a parametric numerical study of long flexurally anisotropic symmetrically laminated plates under combined loading is performed. The loading considered is a combination of tension (compression), shear and in-plane bending. Some generic design charts are presented. The charts may be used for interaction curve prediction with precision of several percents. In [13] a parametric study of long anisotropic plates under combined loading is presented. An important conclusion is made that flexural anisotropy is even more important in case of combined loading than for separate loadings. In [14] a family of design curves is created accounting for flexural anisotropy of a composite plate restrained against inplane lateral and shear deformations. The family is based on numerous direct buckling analysis results for combined loading. Among other papers devoted to direct buckling analysis results for composite plate under combined loading one may indicate [15–21]. In the book of Kassapoglou [15] the combined plate shearcompression is analyzed in case of absence of the flexural anisotropy. The interaction is qualitatively considered using approximate solution with two sine*sine terms and following the Galerkin method. The paper [16] is devoted to creation of the semi-analytical Kantorovich-type approach to direct analysis of composite plates under combined loading. Presented examples demonstrate convergence of the approach. The buckling results for separate loadings are not used in the method. The paper [17] is devoted to buckling FEM analysis of elements used in construction. The paper [18] introduces a discrete model to describe the buckling of stiffened panel beam under a complex loading environment. Numerical results are presented. In [19] a series of tests is detailed, carried out to investigate the behavior of a number of optimized fiber composite plates of differing geometry. The plates are loaded by shear and in-plane bending. The purpose of the study is to assess the suitability of analytical and finite element analysis techniques for predicting structural behavior in the considered conditions. The papers [20,21] take the transverse shear effect into account also. The current paper deals with flexurally anisotropic composite plates subjected to combined loading. The purpose of the paper is to create approaches for predicting buckling loads for flexurally anisotropic composite plates under combined loading, basing on buckling analysis results for separate loadings. The prediction ability is important for the design process as it allows a designer accounting for variations in internal structural load distribution (following from structural changes) without buckling re-analysis. The prediction should be made by some simple approximate formulae, as well as by a method allowing reaching any pre-determined solution accuracy. In the paper the Section 2 is devoted to the theoretical consideration of the problem and analysis of regularities. In the Section 3 the proposed formulae and the developed numerical method are described. The Section 4 is devoted to the numerical results obtained and their consideration. Everywhere in the examples the shear loading is one of the loadings combined with some others. We consider the shear loading direction giving the

lowest buckling eigenvalue, keeping in mind the comments at the beginning of the current paper. The conclusions are made in the Section 5.

2. Theoretical consideration In the paper we consider a rectangular flat thin composite plate. The plate is made of composite material and has anisotropic characteristics, in particular, demonstrates flexurally anisotropic behavior. The plate is loaded by some combined in-plane loading. The sign convention for the separate loadings is shown in the Fig. 1. The coordinate axes x, y are parallel to lateral sides of the plate. The origin is located at the left bottom corner of the plate. In general, the plate is loaded by compression (tension) in two directions and shear. The flows of tension/compression Nx, Ny and the shear flow Nxy may vary along edges and are in equilibrium. The Classical Laminated Plate Theory (CLPT) is used for description of plate deflections at buckling (see the books of Gibson [22] and Reddy [23]). The plate is made of fiber-reinforced material, with the lay-up being symmetric. The boundary conditions considered are the simple support and/or clamping. Following [22,23], we consider the buckling equilibrium of the infinitely small element of the plate. The equation of the force balance in the direction perpendicular to x-y plane is written as:

∂ 2My ∂ 2Mxy ∂Q y ∂ 2Mx ∂Q x + +2 − − =0 2 2 ∂x∂y ∂x ∂y ∂x ∂y

(1)

where Mx, My, Mxy, Q x, Q y respectively are the x and y bending moments, the twisting moment, the buckling-induced forces acting normally to x and y edges of the element, with the moments and forces being calculated for unit element edge length. The Classical Laminated Plate Theory (CLPT), based on linear elasticity relations for every layer and the Kirchhoff hypotheses, is used in (1) for obtaining the deflections w (see [22,9,23]). The equation is written as follows: D11

∂ 4w ∂x 4

+ 2 (D12 + 2D 66 ) − Nx

∂ 2w ∂x

2

− Ny

∂ 4w ∂x 2∂y 2

∂ 2w ∂y

2

+ D 22

− 2Nxy

∂ 4w ∂y 4

+ 4D16

∂ 4w ∂x 3∂y

+ 4D 26

∂ 2w =0 ∂x∂y

∂ 4w ∂x∂y 3

(2)

where Dij , i = 1, 2, 6; j = 1, 2, 6, are the elements of the bending stiffness matrix D, coupling the bending/twisting moments and various second derivatives of the deflection w with respect to x and y. The coupling is written as:

Fig. 1. Rectangular plate. Sign convention.

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

⎡ Mx ⎤ ⎡ D D ⎥ ⎢ 11 12 ⎢ ⎢ My ⎥ = ⎢ D12 D22 ⎢ M ⎥ ⎢⎣ D ⎣ xy ⎦ 16 D26

⎡ ∂ 2w ⎤ ⎥ ⎢ − ∂x2 ⎥ ⎢ D16 ⎤ ⎢ ⎥ ⎥ ∂ 2w ⎥ D26 ⎥ ⎢ − 2 ⎢ ∂y ⎥ D66 ⎥⎦ ⎢ ⎥ 2 ⎢− 2 ∂ w ⎥ ⎢⎣ ∂x∂y ⎥⎦

The buckling Eq. (2) with proper boundary conditions is the socalled eigenvalue problem. Considering Nx, Ny, Nxy as composed of (0) , multiplied by a parasome initial force flow values Nx(0), N y(0), Nxy meter λ, we see that the eigenvalue problem has the classical form

L (w ) + λT (0) (w ) = 0 (3)

¯, Below we will denote the left-hand side of (3) as a vector M and the vector at the right-hand side (which is multiplied to D matrix) as k¯ . Also we use the notation

∂ 2w ∂ 2w ∂ 2w kx = − ; ky = − ; k xy = − 2 2 ∂x∂y ∂x ∂y

(4)

Then (3) is rewritten in the matrix-vector form as:

¯ = Dk¯ M

(5)

For Q x, Q y as functions of external loads the following relation is valid (see [22]):

∂Q y ∂Q x + = Nx k x + Ny k y + 2Nxy k xy ∂x ∂y

(7)

Eq. (2) directly follows from (1),(3), and (6). To complete the boundary-value problem, one must attach two boundary conditions at every edge, considered together with (2). The zero-displacement conditions are:

w=0

for

x = 0, a;

y = 0, b

(8)

For the simple support case the following conditions are valid:

Mx = 0

y = 0, b

(9)

My = 0

x = 0, a

(10)

where the moments correspond to (3). For the plate, clamped at its boundary, the conditions, to be added to (8), are:

∂w =0 ∂x

y = 0, b

∂w =0 ∂y

x = 0, a

(11)

(12)

It is also possible to have some edges clamped and other edges simply supported. The free edge conditions are not considered in the paper. Eq. (2) may be rewritten in the form:

L (w ) + T (w ) = 0

(13)

T (0) (w ) = − Nx(0)

2 ∂ 2w ∂ 2w (0) ∂ w − N y(0) 2 − 2Nxy 2 ∂x∂y ∂x ∂y

⎛

Π˜ =

⎛

∂ 4w ∂x

4

+ 4D16

+ 2 (D12 + 2D66 ) ∂ 4w

∂x 3∂y

+ 4D26

∂ 4w ∂ 4w + D22 4 2 2 ∂x ∂y ∂y

∂ 2w ∂ 2w ∂ 2w T (w ) = − Nx 2 − Ny 2 − 2Nxy ∂x∂y ∂x ∂y

⎞2

⎝

+2D16

+ D12

⎛ ∂ 2w ⎞2 1 ∂ 2w ∂ 2w ⎜ 2⎟ D + 22 2 ⎝ ∂y ⎠ ∂x2 ∂y2

⎛ ∂ 2w ⎞2⎞ ∂ 2w ∂ 2w ∂ 2w ∂ 2w ⎜ ⎟⎟ D 2 2 D + + 26 66 ⎝ ∂x∂y ⎠ ⎟⎠ ∂x2 ∂x∂y ∂y2 ∂x∂y

(18)

The work, performed by the external forces, is

W= −

1 2

⎡

⎞2

⎛

∬ dS ⎢⎢ Nx ⎝ ∂∂wx ⎠ ⎜

⎟

⎣

⎤ ⎛ ∂w ⎞2 ∂w ∂w ⎥ + Ny ⎜ ⎟ + 2Nxy ∂x ∂y ⎥⎦ ⎝ ∂y ⎠

(19)

The total energy U of the system is written as:

(20)

U = Π˜ − W

Due to energy balance at the buckling state w the energy variation δU is equal to zero:

(21)

δU = δΠ˜ − δW = 0

where δ is the variation symbol. Making variations in (18) and (19) and taking into account the boundary conditions of all considered types and in-plane equilibrium used in (7), we obtain the equilibrium Eq. (2). Hence, the stationarity conditions of (20) with respect to w, under simple support or clamped boundary conditions at every edge, lead to the buckling Eq. (2). The statement is a content of the kinematic variational principle describing the plate buckling. The value of U at the stationarity point is equal to zero. The result is easily obtained, using integration by parts and the boundary conditions. Following the book of Vashizu [24], the kinematic variational principle may be rewritten in the form of stationarity of the following ratio (under the same boundary conditions as before):

⎡ Π˜ ⎤ δ ⎢ (0) ⎥ = 0 ⎣W ⎦

(22)

where

⎛

⎞2

∬ dS ⎢⎢ 21 Nx(0)⎝ ∂∂wx ⎠ ⎜

⎣

⎟

+

2 ⎤ 1 (0)⎛ ∂w ⎞ (0) ∂w ∂w ⎥ Ny ⎜ ⎟ + Nxy 2 ∂x ∂y ⎥⎦ ⎝ ∂y ⎠

(23)

The eigenvalue, corresponding to the solution w of the considered buckling problem, may be calculated as the Rayleigh ratio:

∂ 4w ∂x∂y 3

2

∬ dS ⎜⎜ 21 D11⎜⎝ ∂∂xw2 ⎟⎠

W (0) = −

L (w ) = D11

(17)

Resolving the eigenvalue problem, we should determine the value(s) of λ and the corresponding deflection field(s) w. As a rule, the minimal eigenvalue is used in further structural analysis. We indicate now the kinematic variational principles corresponding to (2) and proper boundary conditions. The total strain potential energy Π˜ of the structure at the buckling state is

⎡

where

(16)

where

(6)

The relation (6) with the use of the in-plane equilibrium equations is written as:

∂Q y ∂Q x ∂ ⎛ ∂w ∂w ⎞ ∂ ⎛ ∂w ∂w ⎞ + = − + Nxy + Ny ⎜ Nx ⎟− ⎜ Nxy ⎟ ∂x ∂y ∂x ⎝ ∂x ∂y ⎠ ∂y ⎝ ∂x ∂y ⎠

377

(14)

(15)

λ=

Π˜ W (0)

(24)

It is easy to see that, making variation of (24), we obtain (21) in numerator.

378

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

Talking further on combined loading, we mean that separate loadings like tension or compression in one direction (x or y), shear, in-plane bending are applied simultaneously. The curve corresponding to buckling under such loading (as a curve in plane of corresponding force flows, for example, Nx and Nxy ) we will call the buckling interaction curve or just buckling curve. It is known from the theory of elastic stability (see [3–5]) that in case of structures made of isotropic materials the buckling interaction curve is convex. The statement is known as the Papkovich theorem (called sometimes as the Föppl–Papkovich theorem). There are several proofs of the Papkovich theorem (originally proved by Papkovich in a rather complicated way). Observing the recent one of them indicated in the book of Perelmuter and Slivker [6, vol. 1 Section C.1 ], one may conclude that the proof is valid for buckling of anisotropic plates also. The proof is based on positivity of the second variation of the structural potential energy in case of stable structure and linearity of the variation with respect to structural stress distribution. Considering the behavior of the buckling interaction curve, it is possible to calculate the derivative(s) of the curve, using the known implicit function theorem [25] and (22), (24). Using the notations of the Fig. 1 we have:

shear−compression

∂w ∂w ∂y

∫ dS ∂x dNx = − dNxy 1/2 ∫ dS

( )

(25)

∂w ∂w dS ∂x ∂y

∫ dNxmax = − 2 (y − b / 2) dNxy 1/2 ∫ dS b

∂w 2 ∂x

( )

(26)

where it is supposed that the in-plane bending at the sides parallel to y axis is created by the distributed normal force flow

Nx = Nxmax

Nxy01 + 2Nxy02 +

f 0/ Nxy01Nxy02 Nx0

−

2 f1/ Nxy01Nxy 02

(

Nx0 Nxy01 − Nxy02

2 (y − b/2) b

(27)

The derivatives for other combined cases may be calculated in a similar way. The above approach like (25), (26) may be used for calculation of partial derivatives of the buckling interaction surface in case of combination of more than two loadings.

)

⎛ f / Nxy01Nxy02 ⎞ ⎟⎟ p = − q ⎜⎜ Nxy01 + Nxy02 + 0 Nx0 ⎝ ⎠

(29)

(30)

f0/ , f1/

where are the derivatives of the buckling interaction curve Nx (Nxy ) at points Nx = 0 (pure positive-directed shear) and Nxy = 0 (pure compression), respectively. Below we call (28)–(30) as the 4th order polynomial also. The formulae (28)–(30) give precise values of Nx , Nxy and derivatives of Nx (Nxy ) at points of the coordinate axes. Indeed, calculating derivatives of Nx (Nxy ) at points of the coordinate axes, we obtain:

dNx dNxy

Nxy= 0

⎛ 1 1 p⎞ = − Nx0 ⎜ + + ⎟ = f1/ Nxy02 q⎠ ⎝ Nxy01

Nx = Nx0

dNx dNxy

Nxy = Nxy01

⎧ ⎪ ( Nxy01 − Nxy02 ) ⎡ ⎪ 2 ⎤⎫ / Nxy01) − p ( Nxy01) + q⎥ ⎬ = Nx 0 ⎨ = f0 ( ⎢ ⎪ ⎣ ⎦⎪ ⎩ Nxy01Nxy02 ⎭

Nx = 0

∂w 2 ∂x

where the integrals are calculated over the plate surface, the plate is compressed in x direction, w is the deflection,

shear–in-plane bending

2 Nxy 01Nxy02

q= −

(31)

(32)

Combining (31), (32) and omitting cumbersome algebraic transformations, we obtain (29), (30). For comparison with (28) we indicate the Balabukh formula:

⎛ Nxy ⎞ ⎛ Nxy ⎞ − 1⎟ ⎜ − 1⎟ Nx/Nx0 = − ⎜ N N ⎝ xy01 ⎠ ⎝ xy02 ⎠

(33)

The formula gives precise values of Nx , Nxy at points of the coordinate axes. In case of shear, combined with in-plane bending, the formula of “shifted ellipse” is proposed. It is written as follows: 2 ( Nxmax/Nxmax 0 ) =

⎛ Nxy ⎞ ⎛ Nxy ⎞ −⎜ − 1⎟ ⎜ − 1⎟ ⎝ Nxy01 ⎠ ⎝ Nxy02 ⎠

(34)

max where Nxmax is described in (27) and Nx0 corresponds to buckling under separate in-plane bending. The center of the ellipse (34) is shifted from the origin upwards or downwards because of difference between values of Nxy01 and Nxy02 . Observing the generic charts of [12], one may expect that (34) may give a reasonable quick estimation of the interaction curve for a long plate in case of flexural anisotropy parameters γ , δ (see [12]) being lower than 0.5.

3. Proposed analysis approaches 3.2. Special Galerkin-type approach 3.1. Finite interpolation formulae In the sub-section a possibility to use finite formulas for predicting buckling at combined loadings is considered for the cases of shear–compression and shear-in-plane bending. In case of shear-compression a generalization of the Balabukh formula (see [9]) is proposed. It is written as (hereinafter Nx, Nxy are buckling values for combined loading)

Nx/Nx0 =

⎞ ⎛ Nxy ⎞⎡ 2 1 ⎛ Nxy ⎤ ⎜ − 1⎟ ⎜ − 1⎟ ⎣ ( Nxy ) − p ( Nxy ) + q⎦ q ⎝ Nxy01 ⎠ ⎝ Nxy02 ⎠

(28)

where the subscript 0 means the eigenvalue for one loading, 01 corresponds to the shear eigenvalue in one (considered) direction, 02 corresponds to the shear eigenvalue in case of loading in opposite direction (it has the “-“ sign). In fact, (28) gives Nx as a fourth order polynomial of Nxy . In (28) the quantities p and q are determined as follows:

It is proposed to analyze the buckling behavior under combined loading using the developed and described below Galerkin-type approach (the generic description of the original Galerkin approach one may find in the book of Kantorovich and Krylov [26]). The proposed approach uses eigen buckling modes for separate loadings for performing the corresponding combined loading buckling analysis. It is obvious that the modes satisfy both essential and natural boundary conditions (8)–(10) or (8), (11) and (12). It is supposed that the combined set of buckling modes for separate loadings is a complete set of functions for the boundary problem (2), (8)–(12) of the combined loading. It is also important to note that eigen modes for both loading directions of every separate loading should be used. In general, because of flexural anisotropy, the buckling results for separate loadings are different for direct and opposite loading directions. Observing the variational principle (21) it is clear that the Galerkin approach seeks an energy extremum of (20) using functions

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

satisfying the essential and the natural boundary conditions (see [26]). In our case, following the formalism of the Galerkin approach and denoting the above-indicated eigen functions as wi, i = 1, .. . , N , the eigen mode for combined loading is written as

β=

D12 + 2D66 D11D22 D16

γ= 4

3 D11 D22

4

3 D11D22

N

w=

(35)

i=1

where αi, i = 1, .. . , N are the unknown factors. Substituting (35) in proper relations of the above Section, we obtain the following system of equations for the factors:

⎡ ⎛

N

⎞

⎛

⎠

⎝ i=1

N

⎞⎤

∫S dS ⎢⎢ L ⎜⎜ ∑ αi wi ⎟⎟ + λT (0) ⎜⎜ ∑ αi wi ⎟⎟ ⎥⎥ w j=0 ⎣ ⎝ i=1

⎠⎦

j = 1, …, N (36)

or N

N

⎛

∑ αi ⎝ ∫ ⎜

i=1

=0

S

D16

δ=

∑ αi wi

⎞ ⎛ L (wi ) wj dS⎟ + λ ∑ αi ⎜ ⎠ ⎝ i=1

∫S T (0) (wi ) w j dS⎞⎠

⎟

, j = 1, …, N

(37)

where S is the plate area. The system of relations (37) is a system of homogeneous linear equations for the factors αi, i = 1, .. . , N . The determinant of the system must be equal to zero. Resolving the non-linear algebraic equation of N-th order for the determinant and taking the first positive root λ of the equation, we obtain the eigenvalue for the combined loading. The integrals in (37) are calculated with account of eigen-function properties of particular wi, i = 1, .. . , N for corresponding separate loadings.

4. Examples and numerical results In the Section the results for composite plates under combined loadings are discussed. The plates are made of T300-carbon/ 5208epoxy tape material (see also http://composite.about.com/library/ data/blc-t300-5208.htm) of thickness 0.125 mm. The material is produced by Hexcel Company. The properties of the material are shown in Table 1. The plates have symmetric [[þ 45°,  45°]4]S lay-up (the outermost layer is þ45°, the 0° direction corresponds to the x axis of Fig. 1, the angle is positive in case of rotation from x to y), the total thickness is equal to 2 mm. The plate width is 200 mm. Two aspect ratio a/b values are considered, a/b = 2 and a/b = 3. The plates are simply supported. The results hereinafter in the Section 4 are presented in non-dimensional form, with the stress flow value of 50 N/ mm ( 50,000 N/m) being used for non-dimensioning. The components of the bending stiffness matrix Dij are D11 ¼D22 ¼37.8 Nm, D12 ¼28.2 Nm, D16 ¼D26 ¼5.36 Nm, D66 ¼31.1 Nm, where the elements of the bending stiffness D matrix are determined in (3) (see also [22], Section 7.3). The orthotropy parameter β and the flexural anisotropy parameters γ , δ mentioned in the papers [10–12] are respectively equal to 2.39, 0.142, 0.142 in our case. The formulae for the parameters are:

379

(38)

For long plates the value up to 3.0 is indicated in [10] as the largest value for the buckling interaction formulae proposed there. Below we will use the buckling interaction formulae of [10], intended for use with orthotropic laminates, for comparison. We repeat here that the used in the Section 4 (positive) direction of shear loading corresponds to the lowest shear buckling eigenvalue. It should be mentioned prior to discussion of further results that the precision of 1% in buckling eigenvalue determination is acceptable for majority of applications, for example, aerospace applications. The results for the above-described approaches hereinafter in the Section 4 are compared (for justification) with high-accuracy numerical solutions based on FEM. 4.1. Combined shear-compression. The plate of a/b = 3 is considered as loaded by compression and shear of maximal values Nx = − 50 N/mm ( 50,000 N/m), Nxy = 50 N/mm ( 50,000 N/m). Direct buckling analysis is performed using high-accuracy numerical approach. It is revealed that at the point with Nx/Nxy = 0.8 there is a duplicated eigenvalue equal to 1.14. Left and right derivatives of the buckling curve Nx (Nxy ) at the point are equal to  1.19 and  1.05 (13.3% difference), respectively. In Fig. 2 the buckling interaction curves determined in accordance with the above interpolation formulae are presented, compared to high-accuracy numerical solution. Also at the Figure there are results for the formula (8-1) of the paper of Weaver and Nemeth [10] (with the value of β = 2.39). It is seen at the Fig. 2 that the 4th order polynomial gives results closer to the high-accuracy solution, as compared to ones according to the Balabukh formula, underestimating the buckling

Table 1 T300/5208 Carbon/Epoxy Unidirectional Prepreg. E11, GPa

E22, GPa

G12, GPa

ν12

181

10.3

7.17

0.28

Fig. 2. Numerical results for finite formulae (in comparison with the high-accuracy solution): critical shear stress flow vs critical compression stress flow, a/b¼ 3.

380

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

Table 2 Results for special Galerkin-type approach, combined shear-compression, a/b¼ 3. Number of buckling modes used

Maximal difference in first buckling eigenvalues for combined loading, %

One mode (first or second “positive”) shear þ one mode (first or second) of compression Two modes (first and second) of “positive” shear þ two modes (first and second) of compression Three modes (first, second and third) of “positive” shear þ three modes (first, second and third) of compression

2. 1.7 0.5

values. At the same time, neglecting flexural anisotropy according to the formula of Weaver and Nemeth we come to some difference in obtained results, as compared to the high-accuracy solution, with the results being higher than ones of the solution. This indicates that in case of a/b¼ 3.0 the formula overestimates the buckling eigenvalue. The result is clear from physical standpoint as the account of flexural anisotropy may decrease or increase the buckling values [13] and we deliberately use the shear loading direction leading to the decrease. The various results for the special Galerkin-type approach are shown in Table 2. It is seen from the Table that three “positive” modes of shear (hereinafter “positive” means under loading in positive direction) plus three modes of compression give the results which may be used in practice. The results of the row “One first or second mode shearþ one first or second mode compression” of the Table 2 have been obtained in the following way. For Nx/Nxy ≤ 0.8 (shear-dominating buckling) the first mode in pure “positive” shear and the second mode in compression have been used (the ratio of 0.8 corresponds to the multiple eigenvalue mentioned above). For Nx/Nxy ≥ 0.8 (compression-dominating buckling) the second mode in pure “positive” shear and the first mode in compression have been used. The corresponding non-dimensional results are shown in Fig. 3. It is seen at the Figure that the results rather good enough and largest deviation from the high-accuracy numerical solution is at the center of the plot. The corresponding first and second shear buckling modes, as well as the first and the second compression buckling modes are shown in Figs. 4–7. Due to flexural anisotropy the patterns (as compared to specially orthotropic case) have extra inclination to y axis of Fig. 1.

Fig. 4. First positive shear buckling mode, a/b ¼3.

Fig. 5. Second positive shear buckling mode, a/b¼ 3.

Fig. 6. First compression buckling mode, a/b¼ 3.

Fig. 3. Critical shear stress flow vs critical compression stress flow, a/b¼3.

Fig. 7. Second compression buckling mode, a/b ¼3.

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

Fig. 8. Comparison of Nxy (Nxmax ) curves for a/b¼ 2.

4.2. Combined shear-in-plane bending The plates of a/b = 2 and a/b = 3 are considered as loaded by shear ( 50,000 N/m) and linearly varying in-plane bending (27) with maximal intensity of tension-compression stress flow Nxmax = 500 N/mm ( 500,000 N/m), with the plates being compressed at the upper corners in the in-plane bending case. First we discuss the plate with moderate aspect ratio a/b = 2. Fig. 8 shows the curves Nxy (Nxmax ) (the curves are shown in non-dimensional form) for several approaches. A comparison between the “shifted ellipse” approach and highaccuracy numerical solution leads to a conclusion that the “shifted ellipse” approach gives a rather close approximation of the results for the considered case of the plate under the combined loading (within 2% maximal difference in buckling eigenvalue). The results for special Galerkin-type approach in case of use of two modes (one “positive” shear mode and one “positive” in-plane bending mode) show a big difference between the curves. The choice of the modes (first, second or even third) is made in a way to minimize the difference with the high-accuracy numerical solution. The results for special Galerkin-type approach in case of use of six modes (three “positive” shear modes and three “positive” inplane bending modes) demonstrate up to 1.3% maximal difference with high-accuracy numerical results, with the only point above 1% difference being the second one point (counting from top) at the proper plot. Also the plate with a/b = 3 has been analyzed. Fig. 9 shows the comparison of two curves corresponding to the “shifted ellipse” approach and the high-accuracy numerical solution. It is seen the “shifted ellipse” approach gives conservative results with about 5-7% difference to the high-accuracy numerical solution. Also at the Figure the results for the formula (8–2) of the paper of Weaver and Nemeth [10] are presented. In the formula the effect of flexural anisotropy is neglected. Comparison of the results of the Figure leads to a conclusion that neglecting flexural anisotropy for the plate one obtains over-estimation of buckling values (because of deliberate choice of shear loading direction leading to decrease of the buckling eigenvalue). Also the conclusion of [13] on higher influence of flexural anisotropy in case of

381

Fig. 9. Comparison of Nxy (Nxmax ) curves for a/b¼ 3.

combined loading (as compared to separate loadings) is confirmed. The results for the special Galerkin-type approach are shown in Table 3. It is seen at the Table that even after using twenty modes (ten for shear and ten for in-plane bending) with positive direction of separate loadings it is not possible to obtain 1% closeness to the high-accuracy numerical solution. The closeness may be reached by using 16 modes corresponding to different signs of loadings, namely four initial “positive” and four initial “negative” of both shear and in-plane bending. This results confirms the conclusion that for composite plates with flexural anisotropy it is necessary to use the buckling modes of both (direct and opposite) loading directions. The conclusion demonstrates a considerable difference in buckling analysis and results for composite plates as compared to isotropic plates. Comparing the results of Narita and Leissa [27] and Narita [28] with the results of the Table 3, one may say that instead of solving the determinant 100*100 indicated as acceptable for the Ritz method at the paper, the 16*16 is good enough for the currently considered plate subjected to combined loading and analyzed according to the Galerkin-type approach. 4.3. Combined shear-in-plane bending – longitudinal compression (tension) A possibility to analyze with the developed Galerkin-type approach the loading combined of more than two separate loadings is checked. The plate with a/b¼ 3 is considered. It is loaded by shear with Nxy = 50 N/mm ( 50,000 N/m), in-plane bending with Nxmax = 500 N/mm ( 500,000 N/m) and longitudinal compression (tension) with Nx = ± 50 N/mm (7 50,000 N/m). For the case of compression as a part of the loading the precise eigenvalue is equal to 0.478. The Galerkin-type analysis with four “positive” and four “negative” eigen modes of shear buckling, four “positive” and four “negative” eigen modes of in-plane bending buckling and four eigen modes of compression buckling gives the result with 0.2% accuracy. For the case of tension as a part of the loading the precise eigenvalue is equal to 0.839. The Galerkin-type analysis with four “positive” and four “negative” eigen modes of shear buckling, four

382

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

Table 3 The results for the special Galerkin-type approach, a/b ¼3. Number of modes used

Maximal difference in first buckling eigenvalues for combined loading, %

Four modes (first “positive” and first “negative” of both shear and in-plane bending) Eight modes (two initial “positive” and two initial “negative” of both shear and in-plane bending) Twelve modes (three initial “positive” and three initial “negative” of both shear and in-plane bending) Twenty modes: ten “positive” modes of shear and ten “positive” modes of in-plane bending Sixteen modes (four initial “positive” and four initial “negative” of both shear and in-plane bending)

 20% 46%

   

o 2%  1.6% o 1%

predict with high accuracy the buckling results for combined loading, basing on buckling results for separate loadings; the number of combined separate loadings checked for the proposed Galerkin-type procedure is up to three and may be extended above three; the interaction curve prediction of the proposed Galerkin-type procedure is based on use of, in general, several initial buckling modes for every separate loading; the buckling modes of every separate loading (used for the proposed Galerkin-type procedure) correspond to both direct and opposite loading directions; the convergence rate of the proposed Galerkin-type approach is higher than in case of the Ritz approach based on standard sine functions.

Acknowledgment The author is thankful to colleagues for fruitful discussions. The author is very much indebted to the reviewers for suggesting improvements to the material.

Fig. 10. Shear stress flow vs compression stress flow, a/b¼ 3., Nxmax/Nxy = 10.

“positive” and four “negative” eigen modes of in-plane bending buckling and four eigen modes of compression buckling gives the result with 0.4% accuracy. It is important to emphasize that for the case of tension involved in the combined loading we use the compression buckling modes. Fig. 10 shows a cross-section of the buckling interaction surface in case of Nxmax/Nxy = 10. The plot presents the curve Nxy vs Nx calculated with high-accuracy numerical approach. The proposed Galerkin-type approach gives the results (not shown at the plot) with maximal deviation of 0.3% as compared to the high-accuracy numerical solution. Following the way used in this sub-Section, the proposed Galerkin-type approach of buckling analysis under combined loading may be extended to an arbitrary number of separate loadings.

5. Conclusions Basing on the results of the study, the following conclusions are made:

 in the considered examples the proposed finite formulae give a



conservative estimation of buckling eigenvalue corresponding to combined loading, contrary to the results for the specially orthotropic material approach, over-estimating the buckling eigenvalues; using the proposed Galerkin-type approach it is possible to

References [1] J. Loughlan, The influence of bend–twist coupling on the shear buckling response of thin laminated composite plates, Thin-Walled Struct. 34 (1999) 97–114. [2] J. Loughlan, The shear buckling behaviour of thin composite plates with particular reference to the effects of bend-twist coupling, Int. J. Mech. Sci. 43 (2001) 771–792. [3] S.P. Timoshenko, J.M. Gere, Theory of Elastic Stability, McGraw-Hill, London, 1961. [4] L. Kollar, Structural Stability in Engineering Practice, E & FN Spon, London, 2003. [5] N.A. Alfutov, Stability of Elastic Structures, Springer, New York, 2000. [6] A.V. Perelmuter, V. Slivker, Handbook of Mechanical Stability in Engineering, in Three volumes, World Scientific Publishing Company, Singapore, 2013. [7] H. Schäfer, Beitrag zur berechnung des kleinsten Eigenwertes eindimensionaler Eigenwertprobleme. (in German) [D] 36 S. Diss. Hannover, Braunschweig, Vervielfältigungsanstalt E. Hunold,1937. [8] L. Collatz, Eigenwertaufgaben mit technischen Anwendungen. 2.Durchgesehene Auflage, Geest & Portig K-G, Leipzig, 1963. [9] S.G. Lekhnitskii, Anisotropic Plates, Gordon-Breach, New York, 1968. [10] P.M. Weaver, M.P. Nemeth, Improved design formulae for buckling of orthotropic plates under combined loading, AIAA J. 46 (9) (2008) 2391–2396. [11] NemethMP. Buckling behavior of long symmetrically laminated plates subjected to combined loadings. NASA-TP-3195, 1992. [12] M.P. Nemeth, Buckling of symmetrically laminated plates with compression, shear, and in-plane bending, AIAA J. 30 (12) (1992) 2959–2965. [13] NemethMP. Buckling behavior of long anisotropic plates subjected to combined loads. NASA TP-3568, 1995. [14] M.P. Nemeth, Buckling of long compression-loaded anisotropic plates restrained against inplane lateral and shear deformations, Thin-Walled Struct. 42 (2004) 639–685. [15] C. Kassapoglou, Design and Analysis of Composite Structures, with Applications to Aerospace Structures, Wiley, USA, 2010. [16] I. Shufrin, O. Rabinovitch, M. Eisenberger, Buckling of laminated plates with general boundary conditions under combined compression, tension, and shear —a semi-analytical solution, Thin-Walled Struct. 46 (2008) 925–938.

S. Selyugin / Thin-Walled Structures 98 (2016) 375–383

[17] Q. Liu, P. Qiao, X. Guo, Buckling analysis of restrained orthotropic plates under combined in-plane shear and axial loads and its application to web local buckling, Compos. Struct. 111 (2014) 540–552. [18] H.A. Kim, C.A. Featherstone, J. Ussel, P.A. Williams, Introducing a discrete modelling technique for buckling of panels under combined loading, Struct. Multidiscip. Optim. 36 (1) (2008) 3–13. [19] C.A. Featherston, A. Watson, Buckling of optimised flat composite plates under shear and in-plane bending, Compos. Sci. Technol. 65 (2005) 839–853. [20] N. Rastogi, Buckling of laminated composite plates subjected to combined inplane bi-axial and shear loads. In: Proceedings of the 50th AIAA/ASME/ASCE/ AHS/ASC Structures, Structural Dynamics and Materials Conference, Palm Springs, California, 4–7 May 2009. [21] N.G.R. Iyengar, A. Chakraborty, Study of interaction curves for composite laminate subjected to in-plane uniaxial and shear loadings, Compos. Struct. 64 (2004) 307–315.

383

[22] R.F. Gibson, Principles of Composite Material Mechanics, McGraw-Hill, Inc, New York, 1994. [23] J.N. Reddy, Mechanics of Laminated Composite Plates and Shells. Theory and analysis, 2nd edition, CRC Press, USA, 2004. [24] K. Vashizu, Variational Methods in Elasticity and Plasticity, Pergamon Press, Oxford, 1982. [25] G.A. Korn, T.M. Korn, Mathematical Handbook for Scientists and Engineers, Dover Publ., USA, 2000. [26] L.V. Kantorovich, V.I. Krylov, Approximate Methods of Higher Analysis, Noordhoff, Gröningen, 1964. [27] Y. Narita, A.W. Leissa, Buckling studies for simply supported symmetrically laminated rectangular plates, Int. J. Mech. Sci. 32 (11) (1990) 909–924. [28] Y. Narita, Series and Ritz-type buckling analyses, in: G.J. Turvey, I.H. Marchall (Eds.), Buckling and Postbuckling of Composite Plates, Chapman & Hall, London, 1995, pp. 33–57.