Closed-form approximate solution for linear buckling of Mindlin plates with SRSR-boundary conditions

Journal Pre-proofs Closed-form approximate solution for linear buckling of Mindlin plates with SRSR-boundary conditions M. Beerhorst, S. Thirusala Suresh Babu PII: DOI: Reference:

S0263-8223(19)31270-X https://doi.org/10.1016/j.compstruct.2020.112037 COST 112037

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

9 April 2019 16 December 2019 3 February 2020

Please cite this article as: Beerhorst, M., Suresh Babu, S.T., Closed-form approximate solution for linear buckling of Mindlin plates with SRSR-boundary conditions, Composite Structures (2020), doi: https://doi.org/10.1016/ j.compstruct.2020.112037

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

Closed-form approximate solution for linear buckling of Mindlin plates with SRSR-boundary conditions

M. Beerhorsta,∗, S. Thirusala Suresh Babu b a German Aerospace Center (DLR) Lilienthalplatz 7, D-38108 Braunschweig, Germany b Nanyang Technological University (Graduate Student) 50 Nanyang Avenue, Singapore 639798

Abstract

The present work deals with the buckling analysis of rectangular Mindlin plates consisting of laminated composites with symmetrical, balanced lay-up or isotropic materials. Along the longitudinal edges the plate is rotationally restrained by springs. The transversal edges are simply supported. In agreement with common notation, the boundary conditions are abbreviated as follows: simply supported (S), rotationally restrained (R), and fully clamped (C). As loading situation axial compression is considered. Aiming at high computational eciency the problem is solved by the Rayleigh-Ritz-method. Since well suited shape functions for deection and rotations with very few variables are used, closed-form approximate solutions for the buckling load can be obtained. For verication exact transcendental solutions and/or high delity nite element analyses are employed. Additionally, results are compared to those of existing closed-form approximate solutions. Apart from the special case of simply supported longitudinal edges where all methods yield exact or nearly exact results, the present method shows clear advantages: 1. Due to the type of shape functions it is able deal with unsymmetrical boundary conditions. 2. For the case of both longitudinal edges fully clamped where all closed-form approximate solutions show the greatest deviations the present method is signicantly more accurate. Keywords:

buckling, plate, Mindlin, FSDT, closed-form, rotationally restrained 1. Introduction

In many elds of lightweight design, structures of composite materials nd an ever increasing use. In mechanical models of a real structure, such as ship hulls or aircraft fuselages, the overall structure is often decomposed into substructures, which can be idealized as rods, beams, plates, or shells. Over the decades various analysis methods ranging from global approaches, over the nite strip method (FSM) to the nite element method (FEM) have been developed ∗ Corresponding author

Email addresses:

[email protected] (M.

Beerhorst),

[email protected]

(S. Thirusala Suresh Babu)

Preprint submitted to Elsevier

December 16, 2019

to analyze the stability behavior of such structures. For certain structural congurations exact transcendental solutions have been derived which are often used for method verication. When applying the latter three in design stages of highly iterative character like predesign and optimization, suciently accurate results can often only be achieved at the cost of high computational eort. Because of this, the focus of the present paper is on the development of a sufciently accurate closed-from solution which delivers results in fractions of a second. Thin plates are normally analyzed using the classical laminate plate theory (CLT) which is based on the assumptions of Kirchho's plate theory. An extensive number of works was published on buckling of isotropic (e. g. [1, 2]) and anisotropic plates (e. g. [3]). For thicker plates or plates with relatively low transverse shear stinesses transverse shear deformation has to be taken into account. The simplest form of doing this is to use rst order shear deformation theory (FSDT, [4]). Besides that, various higher order shear deformation theories (HSDT) have been developed [5]. Exact transcendental solutions of the buckling problem can be obtained with so-called Levy-type approaches. For this purpose, an approach must be found for one coordinate direction that exactly fullls the boundary conditions. In the other direction, complex exponential functions, for example, are then used as shape functions. The solutions are exact within the numerical accuracy limits, but can only be determined iteratively due to the transcendental nature of the solution. The success of the iteration heavily depends on the selection of an appropriate starting point because the resulting equations are highly nonlinear and can return complex results during the iteration. Exact solutions for buckling of orthotropic plates according to CLT with longitudinal edges free, simply supported or clamped can be found in [6]. In [7] the case of rotionally restrained and free longitudinal edges is considered. Isotropic plates with linearly varying compressive loads also known as in-plane bending were examined in [8]. FSDT was used in [9] to obtain solutions for isotropic plates under axial compressive load with simply supported transverse edges and various boundary conditions at the longitudinal edges. Further solutions for biaxial loading, two opposing simply supported edges and the other edges either free, hinged or fully clamped were published in [10]. Orthotropic plates and dierent plate theories (CLT, FSDT, HSDT) were considered in [11]. Besides the exact solutions many approximation methods exist which rely on increasing the number of variables in the shape functions to increase the accuracy. Widely used are multi-purpose methods like FEM ([5, 12, 13]) and FSM ([1416]). General Ritz type solution approaches using shape functions with a large number of variables, for example series expansions, were used in [1719]. By far computationally most ecient are so called closed-form solutions which are of explicit nature and directly deliver the buckling load depending on the number of buckling half waves without solving systems of equations. In the special case of a plate with all sides simply supported even the exact solution can be obtained in a closed-form manner (see e. g. [5]). For thin plates with various boundary conditions a huge number of closed-form solutions has been published (e. g. [2026]). The number of works taking into account transverse shear deformation is signicantly lower. FSDT was applied in [27], and [28], 2

kb

y

a

0

Nxx b

k0

z

x

Figure 1: Structural conguration showing the plate with longitudinal edges elastically restrained against rotation under axial compression

Figure 2: Idealization of proled beams or stiened plates by segmentation

HSDT in [29]. Due to the symmetric nature of the approaches only symmetric buckling mode shapes and thus symmetric boundary conditions were taken into account. In the present paper an approach similar to [28] is used but with polynomial shape functions for the deection. These have already shown to deliver results of very high accuracy when dealing with thin plates ([30], [26]). Another advantage is the ability to describe unsymmetrical boundary conditions. To account for transverse shear eects, rotations are also described by polynomial shape functions. As structural conguration a plate of length a, width b and thickness t under 0 axial compression Nxx is considered (see Fig. 1). While the transversal edges are simply supported, the longitudinal edges are elastically restrained against rotation. By assuming dierent spring stinesses k0 at y = 0 and kb at y = b unsymmetric boundary conditions can be accounted for. Possible applications of this model for the local buckling analysis of proled beams or stiened plates are shown in Fig. 2. The structures are segmented and each segment is treated as a plate. This process is also referred to as discrete plate analysis. Restraining eects from adjacent segments are taken into account by springs.

3

2. Governing equations 2.1. Classical formulations, constitutive relations

The displacements u, v , and w in x-, y -, and z -direction of an arbitrary point of the plate are: (1) (2) (3)

u(x, y, z) = u0 (x, y) + zψx (x, y), v(x, y, z) = v0 (x, y) + zψy (x, y), w(x, y) = w0 (x, y).

Herein, the subscript 0 is related to the middle plane of the plate. The rotations of the cross-sections about the x-, y -axis, respectively, are ψy and ψx . A change in the displacement eld causes the following strains: εxx = εyy = εzz = γxy = γxz = γyz =

∂u ∂x ∂v ∂y ∂w ∂z ∂u ∂v + ∂y ∂x ∂u ∂w + ∂z ∂x ∂v ∂w + ∂z ∂y

∂u0 ∂ψx +z , ∂x ∂x ∂v0 ∂ψy = +z , ∂y ∂y

(4)

=

(5) (6)

= 0, ∂u0 ∂v0 + +z ∂y ∂x ∂w0 = + ψx , ∂x ∂w0 = + ψy . ∂y =



∂ψx ∂ψy + ∂y ∂x



,

(7) (8) (9)

Multiplying these with a matrix containing reduced stinesses Qi,j (i, j = 1, 2, 6) and shear stinesses Q44 and Q55 yields the stresses in a single layer:    σ11  Q11         Q12 σ22  τ12 =  0       0  τ23      τ13 0

Q12 Q22 0 0 0

0 0 Q66 0 0

0 0 0 KQ44 0

with the shear correction factor K and E11 , 1 − ν12 ν21 = G23 ,

E22 , 1 − ν12 ν21 = G13 ,

  0 ε11         0    ε22   0  γ12 ,   0  γ23        γ13 KQ55

(10)

ν12 E22 , 1 − ν12 ν21 = G12 . (11)

Q11 =

Q22 =

Q12 =

Q44

Q55

Q66

The resulting forces and moments are determined by integration over the laminate thickness:  0   Nx  σxx              Ny0   Z t/2  σyy   0 Nxy = τxy dz,    −t/2    τyz  Q          y    τxz Qx

   0  Mx  Z t/2 σxx  M0 = σyy zdz.  0y  −t/2  τ  Mxy xy

4

(12)

Substituting the stresses (10) into the above equation and using the relation Aij =

Z

t/2

Qij dz,

Z

Dop =

−t/2

t/2

(13)

Qop z 2 dz, −t/2

yields the material law of a symmetric orthotropic Mindlin plate:  0 Nx      0    N  y    0      N   xy    0  Mx =  My0        0   Mxy         Qy       Qx

A11 A12   0    0   0   0   0 0 

A12 A22 0 0 0 0 0 0

0 0 A66 0 0 0 0 0

0 0 0 D11 D12 0 0 0

0 0 0 D12 D22 0 0 0

0 0 0 0 0 D66 0 0

0 0 0 0 0 0 KA44 0

  ∂u0   0 ∂x     ∂v 0     0    ∂y    ∂v0  ∂u0    + 0  ∂y ∂x      ∂ψx   0  ∂x .  ∂ψy  0   ∂y       ∂ψx + ∂ψy  0     ∂y ∂x       ∂w0 0  + ψ   y  ∂y     ∂w 0 KA55 + ψx ∂x

(14) The bar over the stinesses Q in (13) indicates that they have been transformed into the global (x, y, z)-coordinate system. 2.2. Buckling analysis

To obtain the buckling load the total elastic energy of the system will be minimized by using the Rayleigh-Ritz method. The energy Π consists of three parts: (15) Π = Πi + Πe + Πs . These are the energy of the plate,   2 2 2 ∂ψx ∂ψy ∂ψx ψx ψy ∂ψy + 2D12 + D66 D11 + D22 + ∂x ∂x ∂y ∂y ∂y ∂x 0 0  2  2  ∂w0 ∂w0 + KA44 dydx, (16) + ψy + KA55 + ψx ∂y ∂x

1 Πi = 2

Z

a

Z

b





the external load, Πe = −

0 Nxx 2

Z

a 0

Z

b 0



∂w0 ∂x

2

dydx,

(17)

and the springs, Πs =

1 2

Z

a 0


k0 ψy2 (y = 0) + kb ψy2 (y = b) dx.

(18)

For the unknown deection w0 and rotations ψx and ψy a general variable separable approach is chosen: w0 = W w1 (x)w2 (y),

ψx = Xψx1 (x)ψx2 (y),

5

ψy = Y ψy1 (x)ψy2 (y).

(19)

This makes it possible to separate the integrals in the energy expressions (16) (18). For better overview, abbreviations I1 =

Z

I4 =

Z

I7 =

Z

I10 =

a 0



dψx1 dx

2

dx,

I2 =

Z

I5 =

Z

I8 =

Z

a 2 ψx1 dx, 0 a

w12 dx,

Z0 a 0

a 0 a 0

dψx1 ψy1 dx, dx dψy1 ψx1 dx, dx

a

w1 ψy1 dx, 0

I3 =

Z

a 2 ψy1 dx, 0

2 dψy1 I6 = dx, dx 0 2 Z a dw1 I9 = dx, dx 0 Z

a



dw1 ψx1 dx, dx

(20)

and J1 =

Z

b 2 ψx2 dy, 0

2 dψx2 dy, dy 0 2 Z b dw2 dy, J7 = dy 0 Z b J10 = w2 ψx2 dy, J4 =

J2 =

b

Z



J5 = J8 =

Z Z Z

b 0 b 0 b 0

dψy2 ψx2 dy, dy

J3 =

dψx2 ψy2 dy, dy

J6 =

dw2 ψy2 dy, dy

J9 =

2 C0,rot = ψy2 (y = 0),

Z Z Z

b 0



dψy2 dy

2

dy,

b 2 ψy2 dy, 0 b

w22 dy, 0

2 Cb,rot = ψy2 (y = b), (21)

0

are introduced. The total energy (15) then becomes: 1 1 1 Π = D11 X 2 I1 J1 + D12 XY I2 J2 + D22 Y 2 I3 J3 + D66 X 2 I4 J4 2 2 2 1 1 + D66 XY I5 J5 + D66 Y 2 I6 J6 + KA44 W 2 I7 J7 + KA44 W Y I8 J8 2 2 1 1 2 2 + KA44 Y I3 J6 + KA55 W I9 J9 + KA55 W XI10 J10 2 2 1 1 0 1 2 + KA55 X I4 J1 − Nxx W 2 I9 J9 + Y 2 I3 (k0 C0,rot + kb Cb,rot ) . 2 2 2

(22)

At the minimum point all rst derivatives of the energy must vanish: ∂Π ∂Π ∂Π = 0, = 0, = 0, ∂W ∂X ∂Y

(23)

which yields:
∂Π 0 = KA44 I7 J7 + KA55 I9 J9 − Nxx (24) I9 J9 W ∂W + (KA55 I10 J10 ) X + (KA44 I8 J8 ) Y = 0, (25) ∂Π = (KA55 I10 J10 ) W + (D11 I1 J1 + D66 I4 J4 + KA55 I4 J1 ) X ∂X + (D12 I2 J2 + D66 I5 J6 ) Y = 0, (26) ∂Π = (KA44 I8 J8 ) W + (D12 I2 J2 + D66 I5 J5 ) X ∂Y + (D22 I3 J3 + D66 I6 J6 + KA44 I3 J6 + (k0 C0,rot + kb Cb,rot )) Y = 0.

(27)

6

By collecting coecients of W , X , and Y as: λ11 = KA44 I7 J7 + KA55 I9 J9 , λ11 = I9 J9 , λ22 = D11 I1 J1 + D66 I4 J4 + KA55 I4 J1 , λ33 = D22 I3 J3 + D66 I6 J6 + KA44 I3 J6 + (k0 C0,rot + kb Cb,rot ) , λ12 = KA55 I10 J10 , λ13 = KA44 I8 J8 , λ23 = D12 I2 J2 + D66 I5 J5 ,

(28)

the minimization requirement in vector-matrix-notation becomes: 0 λ11 − Nxx λ11  λ12 λ13

    λ13 W  0 λ23  X = 0 .     0 Y λ33

λ12 λ22 λ32



(29)

The above homogeneous system of equations only has non-trivial solutions if the determinant of the coecient matrix is zero. From this requirement the buckling load is calculated: 0 = Nxx

λ11 (λ22 λ33 − λ223 ) − λ12 (λ12 λ33 − λ13 λ23 ) + λ13 (λ12 λ23 − λ13 λ22 ) . λ11 (λ22 λ33 − λ223 )

In case of Kirchho-type plates the expression simplies to 0 Nxx =

λ22 + λ33 + 2λ23 . λ11

(30) (31)

To calculate the critical buckling load the critical number of buckling halfwaves

mcrit which minimizes the buckling load needs to be determined. Since m can only take integer values an easy to implement algorithm is to start with m = 1 and subsequently increase it by one. If the buckling load calculated for m = n+1 is larger then the one calculated for m = n, then m = n is the critical number 0 of buckling halfwaves and Nxx (m = n) is the critical buckling load. An even more ecient way to obtain mcrit is to drop the assumption that m is an integer. This technique is known as 'relaxation' in the context of optimization problems with integer variables. If m is a continuous variable, mcrit 0 can be determined by demanding that the rst derivative of Nxx with respect to m vanishes: 0 ∂Nxx = 0. (32) ∂m Since the maximum order of m in the numerator and denominator is six and

only even exponents are present:

0 Nxx =

f (m6 , m4 , m2 , m0 ) , f (m6 , m4 , m2 )

(33)

the following substitution is made: ms = m2 .

7

(34)

Now, the orders of ms are 0 = Nxx

f (m3s , m2s , m1s , m0s ) . f (m3s , m2s , m1s )

(35)

For better overview the coecients of the polynomials in ms are denoted as Ni in the numerator and Di in the denominator: P3 Ni mis 0 = P3i=0 . Nxx i i=1 Di ms

(36)

The expressions for the coecients are given in the appendix. Applying the derivation with respect to ms yields  0 ∂Nxx = (D2 N3 − D3 N2 )m4s + (2D1 N3 − 2D3 N1 )m3s ∂ms + (D1 N2 − D2 N1 − 3D3 N0 )m2s − 2D2 N0 ms − D1 N0

.

 . . . . (37)

Since the maximum order of ms in the numerator is four the zeros of this expression can be calculated algebraically. In all observed cases all four solutions ms,crit,i (i = 1, . . . , 4) were real, but only one was positive. This one is kept because only real solutions for mcrit are of interest. Then, mcrit is calculated as p mcrit = ms,crit , (38) again omitting the negative solution. Inserting this value into the buckling load expression (30) gives a lower bound for the buckling load. To this value the buckling load converges if the aspect ratio ab approaches innity. From the continuous real valued mcrit the two possible solutions for the integer valued mcrit are determined by: mlcrit = max (1, oor(mcrit )) , mucrit = max (1, ceil(mcrit )) .

(39) (40)

Then the critical buckling load can be directly obtained by inserting the above values into the expression of the buckling load (30):
0 0 0 Nxx,crit = min Nxx (mlcrit ), Nxx (mucrit ) .

(41)

The halfwave number which leads to the lower buckling load is the critical one. In this case no iteration over m is necessary. To present the results in a generic manner the non-dimensional buckling coecient kx =

0 Nxx,crit b2 √ π 2 D11 D22

(42)

is used. 2.3. Boundary conditions

The boundary conditions of the considered structural conguration (see Fig. 1) are: 8

at y = [0, b] :

At x = [0, a] :

(43) (44) (45)

w0 = [0, 0], ψy = [0, 0], 0 = [0, 0], Mxx

w0 = [0, 0], ψx = [0, 0], 0 Myy = [k0 ψy , −kb ψy ].

(46) (47) (48)

2.4. Shape functions

When employing the Ritz-method it is mandatory that the shape functions satisfy the kinematic boundary conditions involving deection w0 and rotations ψx and ψy . Additional fulllment of the natural boundary conditions regarding the bending moments improves the accuracy of the solution. While the general form of the shape functions was already given in (19), particular shape functions in longitudinal and transversal direction yet need to be determined. Since the transversal edges are simply supported trigonometric shape functions are used in longitudinal direction ((49), (51), (53)). They identically satisfy boundary conditions (43) - (45). w1 (x) = sin

 mπx  a

,  y 2

y + W2 b b

(49)  y 3

 y 4

.

(50)

 mπx  mπ , cos a a  y 2  y 3  y 4 y + W3 + W4 . ψx2 (y) = + W2 b b b b

(51)

w2 (y) =

+ W3

b

+ W4

b

ψx1 (x) =

ψy1 (x) = sin ψy2 (y) =

 mπx  a

,  2 y

1 + W2 b b

b

(52) (53)

+ W3

3  y 2 b

b

+ W4

4  y 3 b

b

.

(54)

In transversal direction the shape functions need to be exible enough to simulate dierent degrees of rotational restraint ranging from the limit cases of simply supported to fully clamped. This is achieved by using a fourth order polynomial in y for deection function w2 ((50), [30]) and rotation function ψx2 (52) and a third order polynomial for rotation function ψy2 (54). The unknown coecients W2 to W4 can be obtained from boundary conditions (46) - (48) as k0 b , 2D22 −12D22 2 − 5b (k0 + 3/5kb ) D22 − kb b2 k0 , W3 = D22 (kb b + 6D22 )

(55)

12D22 2 + 4b (k0 + kb ) D22 + kb b2 k0 . 2D22 (kb b + 6D22 )

(57)

W2 =

W4 =

9

(56)

Now that the shape functions are fully specied, the integrals (20) and (21) can be evaluated: 1 m4 π 4 , 2 a3 I3 = I7 = I8 ,

1 m2 π 2 , 2 a I4 = I5 = I6 = I9 = I10 .

I1 =

J1 =

b 5



W2 2 +



I2 = −

5 10W4 5 W3 + + 3 7 2



I3 =

5 W2 + W 3 2 + 7



1 , 2

I4 = −I2 ,

(58)

  5 5 5 5 , W4 + 2 W 3 + W 4 2 + W 4 + 4 9 3 3

70W2 2 + (210W3 + 294W4 + 105) W2 + 126W3 2 + (315W4 + 210) W3 + 180W4 2 + 315W4 , 105b 20W2 2 + (60W3 + 80W4 ) W2 + 60W3 2 + 180W3 W4 + 144W4 2 J3 = , 5b3 140W2 2 + (315W3 + 336W4 + 210) W2 + 189W3 2 + (420W4 + 210) W3 + 240W4 2 + 210W4 + 105 J4 = 105b J 4 = J5 = J6 = J7 = J8 . (59) J1 = J9 = J10 , J2 =

C0,rot =

1 , b2

Cb,rot =



W2 W3 W4 1 +2 +3 +4 b b b b

2

.

(60)

3. Results

In this section the results of the developed analytical approximation procedure are compared to results obtained by other methods. Regarding the boundary conditions the transversal edges are simply supported. At the longitudinal edges three dierent sets of boundary conditions are considered: both simply supported (SS), one fully clamped and one simply supported (CS), or both fully clamped (CC). Verication baseline are the exact solutions given in [10] for isotropic plates with various boundary conditions and high delity nite element analyses (FEA) for all materials and boundary conditions. As element type the eight-noded layered element Shell281 [31] is chosen. The width of the plate is kept constant and the length is varied to achieve the desired aspect ratio. A mesh division of 10 elements in transversal direction and the corresponding division in longitudinal direction to obtain nearly square shaped elements proved to be suciently high during convergence studies. The results are also compared to existing closed-form approximate solutions ([27], [28]). Since both use trigonometric shape functions they will be referred to as 'trig1' and 'trig2', respectively. The newly developed approach uses trigonometric and polynomial shape functions and will be referred to as 'poly'. For SS boundary conditions the 'trig2' closed-form solution [28] is identical to the exact solution. Three types of materials are considered: isotropic (Aluminum 7075-T6), orthotropic (HexPly M21 UD IMA), and woven fabric (HexPly M21 Woven AS4C). Material properties have been taken from the material database of ESAComp [32]. The lay-up for orthotropic plates is ((0/90)n )s and for woven fabric ((±45)n )s . Subscript n denotes the number of repetitions and s the symmetry. In case of orthotropic material the layer thickness is 0.125 mm, in case of woven fabric a ±45 layer is 0.250 mm thick. Table 1 lists the material properties of a single layer. The width of the plate is xed to b = 100 mm 10

5.0 Trig1

Poly

Exact, Trig2

b/t = 200 ↓ 10 ↓

4.5

kx

FEM

4.0

3.5

3.0 0.5

5↓ 1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 3: Buckling coecient vs. aspect ratio of isotropic plates with SS boundary conditions. and the apect ratio is varied between 0.5 and 5. Dierent ratios of the relative thickness b/t are examined to show the inuence on the nondimensional buckling load. The values are given in Tab. 2 together with total plate thickness t and the corresponding number of layer repetitions n.

Table 1: Material properties Iso Ortho Woven

E11 /GPa

71.00 162.00 63.65

E22 /GPa

71.00 10.00 63.65

G12 /GPa

27.31 5.20 4.20

G13 /GPa

27.31 4.45 4.00

G23 /GPa

27.31 4.45 4.00

ν12

0.30 0.30 0.05

Table 2: Relative thicknesses b/t, thicknesses t and number of layer repetitions n b/t t/mm n

200 0.5 1

100 1 2

50 2 4

20 5 10

10 10 20

5 20 40

Figures 3 to 5 depict the variation of nondimensional buckling coecient with respect to the plate aspect ratio for simply supported, rectangular plates made up of isotropic, orthotropic and woven fabric layers, respectively. All buckling curves show the expected garlands shape with peaks at the positions where the halfwave numbers change. With increasing aspect ratio the amplitudes of the peaks decrease and the results converge to the minima of the respective curves. Transverse shear eects are negligible for thin plates with b/t = 200, but become clearly relevant for b/t ≤ 10. The results of the dierent methods are identical or very close together. Numerical values for the assessment of the accuracies of the dierent methods are given in Tab. 3. It lists the root mean squared relative errors (RMSRE) with respect to the exact solutions. The RMSRE have been calculated for 46 evenly distributed a/b-values in the interval from 0.5 to 5.0. As mentioned above, the 'trig2' approach delivers the exact solution. Therefore it is used as reference and 11

4.0 Trig1

Poly

Exact, Trig2

FEM

3.5

b/t =

kx

3.0

2.5

200 ↓

2.0

10 ↓ 5↓

1.5

1.0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 4: Buckling coecient vs. aspect ratio of orthotropic plates with SS boundary conditions.

4.0 Trig1

Poly

Exact, Trig2

FEM

3.5

b/t =

kx

3.0

200 ↓ 10 ↓

2.5

2.0

5↓

1.5

1.0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 5: Buckling coecient vs. aspect ratio of woven fabric plates with SS boundary conditions.

12

7.0 Poly

Exact

FEM

6.5

b/t = 200 ↓

6.0

kx

5.5

10 ↓

5.0 4.5

5↓

4.0 3.5 3.0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 6: Buckling coecient vs. aspect ratio of isotropic plates with CS boundary conditions.

no deviations are listed in Tab. 3 for this approach. For isotropic material 'trig1' also delivers the exact solution and the error in case of 'poly' is less than 0.0%. That is why these results are omitted in the table. In case of other materials the accuracy of 'trig1' and 'poly' is also very high. Table 3: SS boundary conditions: RMSRE of buckling coecients for dierent approaches

Mat

b/t

Trig1 [27]

Poly

Ortho

200 100 50 20 10 5

0.0% 0.0% 0.1% 0.3% 0.6% 1.7%

0.1% 0.1% 0.1% 0.1% 0.0% 0.0%

Woven

200 100 50 20 10 5

0.0% 0.0% 0.0% 0.2% 0.5% 1.2%

0.1% 0.1% 0.1% 0.1% 0.0% 0.0%

Figures 6 to 8 show buckling curves for plates with CS boundary conditions made up of isotropic, orthotropic and woven fabric layers, respectively. Since the 'trig' approaches are not able to represent unsymmetrical boundary conditions only 'poly' results are compared to the reference results. The closed-form approach accurately describes the inuence of aspect ratio and relative thickness. This can also be seen by the RMSRE in Tab. 4. Considering all materials, the maximum RMSRE are [0.5%, 3.9%, 4.8%] for b/t = [200, 10, 5]. Results for plates made up of isotropic, orthotropic and woven fabric layers, respectively, with longitudinal edges fully clamped are shown in Fig. 9 to 11. 13

Table 4: CS boundary conditions: RMSRE of buckling coecients for poly approach

Mat

Poly 0.5% 0.5% 0.6% 1.0% 2.3% 4.8% 0.1% 0.1% 0.3% 1.5% 3.9% 3.1% 0.1% 0.1% 0.3% 1.4% 3.7% 4.0%

b/t

200 100 50 Iso 20 10 5 200 100 Ortho 50 20 10 5 200 100 Woven 50 20 10 5

5.0 Poly

FEM

b/t = 200 ↓

4.5 4.0

kx

3.5

10 ↓

3.0 2.5 2.0

5↓

1.5 1.0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 7: Buckling coecient vs. aspect ratio of orthotropic plates with CS boundary conditions.

14

5.0 Poly

FEM

b/t = 200 ↓

4.5 4.0

kx

3.5

10 ↓

3.0 2.5 2.0

5↓

1.5 1.0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 8: Buckling coecient vs. aspect ratio of woven fabric plates with CS boundary conditions. 9.0 Trig1

Trig2

Poly

Exact

FEM

b/t = 200 ↓

8.0

10 ↓

kx

7.0

6.0

5↓

5.0

4.0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 9: Buckling coecient vs. aspect ratio of isotropic plates with CC boundary conditions. For thin plates, results are nearly identical to reference results in case of the 'poly' approach. The 'trig' approaches both show already larger deviations for thin plates, especially in case of the isotropic plate. Since their results are nearly identical for b/t = 200, the curves lie on top of each other. For all b/t-values the improvement brought by the newly developed method is signicant. Table 5 lists the RMSRE. The maximum RMSRE are [0.2%, 4.1%, 6.4%] for b/t = [200, 10, 5] in case of the 'poly' approach compared to [4.5%, 7.8%, 10.7%] in case of the best trig approach. 4. Concluding remarks

Polynomial shape functions which partly share variables in deection and rotation functions have been used to achieve closed-form approximate solutions for buckling loads and numbers of buckling halfwaves of orthotropic Mindlin type plates.

15

Table 5: CC boundary conditions: RMSRE of buckling coecients for dierent approaches

Mat

b/t

200 100 50 Iso 20 10 5 200 100 Ortho 50 20 10 5 200 100 Woven 50 20 10 5

Trig1 Trig2 Poly [27] [28] 4.5% 4.5% 0.2% 4.5% 4.5% 0.2% 4.7% 4.7% 0.4% 5.8% 5.5% 1.2% 8.6% 7.8% 3.3% 12.4% 10.7% 6.4% 2.3% 2.3% 0.1% 2.4% 2.4% 0.2% 2.8% 2.7% 0.5% 5.1% 4.4% 1.9% 7.7% 6.8% 4.1% 3.6% 3.0% 2.2% 2.3% 2.3% 0.1% 2.5% 2.4% 0.2% 2.9% 2.7% 0.4% 5.2% 4.4% 1.8% 8.1% 6.9% 4.0% 4.9% 4.7% 3.2%

7.0 Trig1

Trig2

Poly

FEM

b/t = 200 ↓

6.0

kx

5.0

4.0

10 ↓ 3.0

2.0

1.0 0.5

5↓ 1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 10: Buckling coecient vs. aspect ratio of orthotropic plates with CC boundary conditions.

16

Trig1

6.0

Trig2

Poly

FEM

b/t = 200 ↓

5.0

4.0

kx

10 ↓

3.0

5↓

2.0

1.0 0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

a/b

Figure 11: Buckling coecient vs. aspect ratio of woven fabric plates with CC boundary conditions. With these closed-form solutions results can be generated with negligible computational eort. This is a major advantage compared to methods like FSM or FEM for application in analysis processes with highly iterative character, for example predesign and optimization. Due to the analytical nature exact sensitivities can be calculated if needed for optimization. Implementation into analysis tools is also very easy. For a range of boundary conditions, materials and relative plate thicknesses the accuracy of the newly developed method has been benchmarked against previously developed closed-form solutions. Results from detailed FE models and exact transcendental solutions have been taken as references. Compared to the other closed-form solutions the present method has clear advantages: 1. Due to the type of shape functions it is able to deal with unsymmetrical boundary conditions. 2. For the boundary condition with both longitudinal edges fully clamped where all closed-form approximate solutions show the greatest deviations with respect to reference results the present method is signicantly more accurate. For relative plate thicknesses b/t = [200, 10, 5] the maximum root mean squared relative errors are [0.2%, 4.1%, 6.4%] compared to [4.5%, 7.8%, 10.7%].

Appendix Coecients N0 to N3 and D1 to D3 from (36): 2

2

2

N0 = A44 A55 I 3 I 4 I 7 J1 J6 J7 − A44 A55 I 4 I 8 J1 J82 2

2

2

+ A44 D66 I 3 I 4 I 7 J4 J6 J7 − A44 D66 I 4 I 8 J4 J82 + A44 A55 C0,rot k0 I 3 I 4 I 7 J1 J7 + A44 A55 Cb,rot kb I 3 I 4 I 7 J1 J7 + A44 A55 D22 I 3 I 4 I 7 J1 J3 J7 + A44 C0,rot D66 k0 I 3 I 4 I 7 J4 J7 + A44 Cb,rot D66 kb I 3 I 4 I 7 J4 J7 + A44 D22 D66 I 3 I 4 I 7 J3 J4 J7 ,

17

(61a)

2

2

2

2

N1 = A44 D11 I 1 I 3 I 7 J1 J6 J7 − A44 D11 I 1 I 8 J1 J82 + A44 A55 I 3 I 4 I 9 J1 J6 J9 2

2

2 − A44 A55 I 3 I 10 J6 J10 + 2A44 A55 D12 I 2 I 8 I 10 J2 J8 J10 + A44 A55 D66 I 3 I 4 I 9 J4 J6 J9 + A44 A55 D66 I 4 I 6 I 7 J1 J6 J7 + 2A44 A55 D66 I 5 I 8 I 10 J5 J8 J10 + A44 C0,rot D11 k0 I 1 I 3 I 7 J1 J7 + A44 Cb,rot D11 kb I 1 I 3 I 7 J1 J7 + A44 D11 D22 I 1 I 3 I 7 J1 J3 J7 2

2 2 − A44 D12 I 2 I 7 J22 J7 − 2A44 D12 D66 I 2 I 5 I 7 J2 J5 J7 + A44 D66 I 4 I 6 I 7 J4 J6 J7 2

2

2

2

2 2 − A44 D66 I 5 I 7 J52 J7 + A55 C0,rot k0 I 3 I 4 I 9 J1 J9 − A55 C0,rot k0 I 3 I 10 J10

+

2 2 2 2 2 A55 Cb,rot kb I 3 I 4 I 9 J1 J9 − A55 Cb,rot kb I 3 I 10 J10 + A55 D22 I 3 I 4 I 9 J1 J3 J9 2 2 2 A55 D22 I 3 I 10 J3 J10 + A55 C0,rot D66 k0 I 3 I 4 I 9 J4 J9

− + A55 Cb,rot D66 kb I 3 I 4 I 9 J4 J9 + A55 D22 D66 I 3 I 4 I 9 J3 J4 J9 ,

(61b) N2 = A44 A55 D11 I 1 I 3 I 9 J1 J6 J9 + A44 D11 D66 I 1 I 6 I 7 J1 J6 J7 2

2

2

2 + A55 D66 I 4 I 6 I 9 J1 J6 J9 − A55 D66 I 6 I 10 J6 J10 + A55 C0,rot D11 k0 I 1 I 3 I 9 J1 J9 + A55 Cb,rot D11 kb I 1 I 3 I 9 J1 J9

(61c)

2

2 + A55 D11 D22 I 1 I 3 I 9 J1 J3 J9 − A55 D12 I 2 I 9 J22 J9

2

2 2 − 2A55 D12 D66 I 2 I 5 I 9 J2 J5 J9 + A55 D66 I 4 I 6 I 9 J4 J6 J9 − A55 D66 I 5 I 9 J52 J9 , N3 = A55 D11 D66 I 1 I 6 I 9 J1 J6 J9 , (61d)

D1 = A44 A55 I 3 I 4 I 9 J1 J6 J9 + A44 D66 I 3 I 4 I 9 J4 J6 J9 + A55 C0,rot k0 I 3 I 4 I 9 J1 J9 + A55 Cb,rot kb I 3 I 4 I 9 J1 J9 + A55 D22 I 3 I 4 I 9 J1 J3 J9 + C0,rot D66 k0 I 3 I 4 I 9 J4 J9 + Cb,rot D66 kb I 3 I 4 I 9 J4 J9 + D22 D66 I 3 I 4 I 9 J3 J4 J9 ,

(62a)

D2 = A44 D11 I 1 I 3 I 9 J1 J6 J9 + A55 D66 I 4 I 6 I 9 J1 J6 J9 + C0,rot D11 k0 I 1 I 3 I 9 J1 J9 + Cb,rot D11 kb I 1 I 3 I 9 J1 J9

(62b)

2

2 + D11 D22 I 1 I 3 I 9 J1 J3 J9 − D12 I 2 I 9 J22 J9

2

2 2 − 2D12 D66 I 2 I 5 I 9 J2 J5 J9 + D66 I 4 I 6 I 9 J4 J6 J9 − D66 I 5 I 9 J52 J9 , D3 = I 9 J9 D11 D66 I 1 I 6 J1 J6 , (62c)

with Aii = KAii ,

i = 4, 5.

(63)

The I i can be obtained from the Ii in (58) by omitting all occurrences of m.

Data availability The raw/processed data required to reproduce these ndings cannot be shared at this time due to technical or time limitations.

18

References [1] F. Bleich, Buckling strength of metal structures, McGraw-Hill, New York, 1952 (1952). [2] S. P. Timoshenko, J. Gere, Theory of elastic stability, McGraw-Hill, New York, 1961 (1961). [3] G. J. Turvey, I. H. Marshall, Buckling and postbuckling of composite plates, Chapman and Hall, London, 1995 (1995). [4] R. D. Mindlin, Inuence of rotatory inertia and shear on exural motion of isotropic, elastic plates, Journal of Applied Mechanics 18 (1951) 3138 (1951). [5] J. Reddy, Mechanics of laminated composite plates and shells, 2nd Edition, CRC PRESS, Boca Raton et al., 2004 (2004). [6] E. J. Brunelle, G. A. Oyibo, Generic buckling curves for specially orthotropic rectangular plates, AIAA Journal 21 (8) (1982) 11501156 (1982). [7] L. C. Bank, J. Yin, Buckling of orthotropic plates with free and rotationally restrained unloaded edges, Thin-Walled Structures 24 (1996) 8396 (1996). [8] J.-H. Kang, A. W. Leissa, Exact solutions for the buckling of rectangular plates having linearly varying in-plane loading on two opposite simply supported edges, International Journal of Solids and Structures 42 (14) (2005) 42204238 (Jul. 2005). doi:10.1016/j.ijsolstr.2004.12.011. URL http://www.sciencedirect.com/science/article/pii/S0020768304006705 [9] E. J. Brunelle, Buckling of transversely isotropic Mindlin plates, AIAA Journal 9 (6) (1971) 10181022 (Jun. 1971). doi:10.2514/3.6326. URL https://doi.org/10.2514/3.6326 [10] S. Hosseini-Hashemi, K. Khorshidi, M. Amabili, Exact solution for linear buckling of rectangular Mindlin plates, Journal of Sound and Vibration 315 (1) (2008) 318342 (Aug. 2008). doi:10.1016/j.jsv.2008.01.059. URL http://www.sciencedirect.com/science/article/pii/S0022460X08001119 [11] A. Khdeir, L. Librescu, Analysis of symmetric cross-ply laminated elastic plates using a higher-order theory: Part IIBuckling and free vibration, Composite Structures 9 (4) (1988) 259277 (Jan. 1988). doi:10.1016/02638223(88)90048-7. URL http://www.sciencedirect.com/science/article/pii/0263822388900487 [12] G. Rao, J. Venkataramana, K. Raju, Stability of moderately thick rectangular plates using a high precision triangular nite element, Computers & Structures 5 (4) (1975) 257259 (Nov. 1975). doi:10.1016/00457949(75)90028-0. URL http://www.sciencedirect.com/science/article/pii/0045794975900280 [13] K. Rohwer, Improved transverse shear stinesses for layered nite elements, Tech. Rep. DFVLR-FB-88-32, DFVLR (1988).

19

[14] G. Zou, P. Qiao, Higher-Order nite strip method for postbuckling analysis of imperfect composite plates, Journal of Engineering Mechanics 128 (9) (2002) 10081015 (2002). [15] H. C. Bui, J. Rondal, Buckling analysis of thin-walled sections by semianalytical MindlinReissner nite strips: A treatment of drilling rotation problem, Thin-Walled Structures 46 (6) (2008) 646652 (Jun. 2008). doi:10.1016/j.tws.2007.12.003. URL http://www.sciencedirect.com/science/article/pii/S0263823107003023 [16] H. R. Ovesy, E. Zia-Dehkordi, S. A. M. Ghannadpour, High Accuracy Postbuckling Analysis of Moderately Thick Composite Plates Using an Exact Finite Strip, Computers and Structures 174 (C) (2016) 104112 (Oct. 2016). doi:10.1016/j.compstruc.2015.09.014. URL https://doi.org/10.1016/j.compstruc.2015.09.014 [17] R. Vescovini, L. Dozio, M. D'Ottavio, O. Polit, On the application of the ritz method to free vibration and buckling analysis of highly anisotropic plates, Composite Structures 192 (2018) 460  474 (2018). doi:https://doi.org/10.1016/j.compstruct.2018.03.017. URL http://www.sciencedirect.com/science/article/pii/S0263822318305233 [18] J. Yang, D. Chen, S. Kitipornchai, Buckling and free vibration analyses of functionally graded graphene reinforced porous nanocomposite plates based on chebyshev-ritz method, Composite Structures 193 (2018) 281  294 (2018). doi:https://doi.org/10.1016/j.compstruct.2018.03.090. URL http://www.sciencedirect.com/science/article/pii/S0263822318307864 [19] V. Oliveri, A. Milazzo, A rayleigh-ritz approach for postbuckling analysis of variable angle tow composite stiened panels, Computers & Structures 196 (2018) 263  276 (2018). doi:https://doi.org/10.1016/j.compstruc.2017.10.009. URL http://www.sciencedirect.com/science/article/pii/S004579491730768X [20] L. P. Kollar, G. S. Springer, Mechanics of composite structures, Cambridge University Press, Cambridge, 2003 (2003). [21] P. Qiao, L. Shan, Explicit local buckling analysis and design of berreinforced plastic composite structural shapes, Composite Structures 70 (4) (2005) 468483 (2005). [22] C. Mittelstedt, Local buckling of wide-ange thin-walled anisotropic composite beams, Archive of Applied Mechanics 77 (2007) 439452 (Jul. 2007). doi:10.1007/s00419-006-0102-0. URL https://doi.org/10.1007/s00419-006-0102-0 [23] G. Tarján, L. P. Kollár, Buckling of axially loaded composite plates with restrained edges, Journal of Reinforced Plastics and Composites 29 (23) (2010) 35213529 (2010). doi:10.1177/0731684410381153. URL https://doi.org/10.1177/0731684410381153 [24] M. Beerhorst, M. Seibel, C. Mittelstedt, Fast analytical method describing the postbuckling behavior of long, symmetric, balanced laminated composite plates under biaxial compression and shear, Composite Structures 94 (6) (2012) 20012009 (2012). 20

[25] L. P. Kollár, Buckling of rectangular composite plates with restrained edges subjected to axial loads, Journal of Reinforced Plastics and Composites 33 (23) (2014) 21742182 (2014). doi:10.1177/0731684414555880. URL https://doi.org/10.1177/0731684414555880 [26] M. Beerhorst, M. Seibel, C. Mittelstedt, Ecient approximate solutions for postbuckling of axially compressed orthotropic plates with rotationally restrained and free longitudinal edges, Archive of Applied Mechanics 88 (3) (2018) 461476 (Mar. 2018). doi:10.1007/s00419-017-1319-9. URL https://doi.org/10.1007/s00419-017-1319-9 [27] H. W. March, Eects of shear deformation in the core of a at rectangular sandwich panel, Tech. Rep. 1583, Forest Products Laboratory (1948). [28] T. Kuehn, H. Pasternak, C. Mittelstedt, Local buckling of shear-deformable laminated composite beams with arbitrary cross-sections using discrete plate analysis, Composite Structures 113 (2014) 236248 (Jul. 2014). [29] J. Herrmann, T. Kühn, T. Müllenstedt, S. Mittelstedt, C. Mittelstedt, Closed-Form Approximate Solutions for the Local Buckling Behavior of Composite Laminated Beams Based on Third-Order Shear Deformation Theory, in: H. Altenbach, F. Jablonski, W. H. Müller, K. Naumenko, P. Schneider (Eds.), Advances in Mechanics of Materials and Structural Analysis: In Honor of Reinhold Kienzler, Advanced Structured Materials, Springer International Publishing, Cham, 2018, pp. 175205 (2018). doi:10.1007/978-3-319-70563-7_8. URL https://doi.org/10.1007/978-3-319-70563-7_8 [30] J. Rhodes, J. Harvey, Plates in uniaxial compression with various support conditions at the unloaded boundaries, International Journal of Mechanical Sciences 13 (9) (1971) 787802 (Sep. 1971). doi:10.1016/00207403(71)90047-6. URL http://www.sciencedirect.com/science/article/pii/0020740371900476 [31] ANSYS Inc., Release 15.0 Documentation for ANSYS, Cannonsburg, Pennsylvania, U.S.A., 2013 (2013). [32] Altair Engineering, Inc., ESAComp 4.7 Documentation, Troy, Michigan, U.S.A., 2018 (2018).

21