Novel quadrature element formulation for simultaneous local and global buckling analysis of eccentrically stiffened plates

Aerospace Science and Technology 87 (2019) 154–166

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Novel quadrature element formulation for simultaneous local and global buckling analysis of eccentrically stiffened plates Jian Deng a , Xinwei Wang a,∗ , Zhangxian Yuan b , Guangming Zhou a a b

State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China School of Aerospace Engineering, Georgia Institute of Technology, GA 30332-0150, USA

a r t i c l e

i n f o

Article history: Received 14 December 2018 Received in revised form 18 January 2019 Accepted 18 February 2019 Available online 22 February 2019 Keywords: Eccentrically stiffened plate Quadrature plate-stiffener element Local buckling Global buckling High-order finite element method

a b s t r a c t Predicting both buckling load and mode shape of eccentrically stiffened panels correctly is of important and a rather challenging task. In this paper, a novel and efficient quadrature element formulation is developed to fulfill this challenging task. A high-order quadrature plate-stiffener element is proposed by assembling quadrature beam elements to the quadrature plate element via the displacement relations. Thus, the same plate mesh scheme can be used for buckling analysis of stiffened plate with different numbers of stiffeners located at arbitrary positions. Explicit formulations and solution procedures are given. Convergence studies are performed for stiffened plates with different shapes of cross-section of stiffeners and materials. A number of case studies are given. For validations, numerical results are compared with either existing solutions or finite element data. It is demonstrated that high accuracy can be achieved with relatively small number of nodes. Presented formulation is simple, straightforward, and reliable, which can allow a quick and accurate analysis of buckling behavior of eccentrically stiffened plates. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction Eccentrically stiffened plates and shells are common structural elements and widely employed in aeronautical and aerospace structures. Therefore, their static, buckling and dynamic behaviors have been received great attention. A vast body of literature exists and is well documented in [1]. A wide range of analytical, approximate and numerical procedures developed over several decades as well as various engineering applications have been reviewed. Due to the thin-walled in nature, eccentrically stiffened panels are susceptible to buckle [1,2]. Therefore, the buckling behavior is one of the major concerns by researchers and engineers. Due to small difference in buckling loads of the first several modes, however, predicting both buckling load and mode shape of eccentrically stiffened panels correctly is a rather challenging task. Various analytical, approximate and numerical methods, such as Kantorovich method [3], Rayleigh-Ritz method [3–5], the differential quadrature method (DQM) [6], the differential quadrature element method (DQEM) [7,8], and the finite element method (FEM) [9–11], have been developed and/or used to obtain the buckling behavior of stiffened panels thus far.

*

Corresponding author. E-mail address: [email protected] (X. Wang).

https://doi.org/10.1016/j.ast.2019.02.019 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

For stiffened plates and shells having sparsely spaced stiffeners, either global or local instability may be the first buckling mode and thus a general approach which can predict both local and global buckling behavior of eccentrically stiffened panels simultaneously is beneficial for the designers. The conventional FEM is a versatile method [11] and thus widely used at present time [1,9, 10]. However, it is not efficient to employ the FEM during the preliminary design stage, since the dimensions of the stiffened panels and stiffeners are not finalized and to be optimally designed [12]. For laminated composite stiffened panels, the stacking sequence is also to be optimally designed. Therefore, there is still a need for a relatively inexpensive yet accurate method to obtain the buckling behavior of eccentrically stiffened plates and shells. The discrete singular convolution (DSC) method, originated by Wei [13], was used in analyzing static, buckling and dynamic behavior of structural elements by Civalek [14–16] and Lai et al. [17, 18]. The differential quadrature method (DQM) [6,19–21] and the differential quadrature element method (DQEM) [7,8,20,21] were also employed in analyzing the static, buckling and dynamic behavior of beams, plates and shells. It was demonstrated that these strong form methods were simple and could obtain accurate solutions with a minimum number of grid points, especially for structural elements having regular shapes. However, these methods might be still relatively expensive when the eccentrically stiffened

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

plates contain a number of stiffeners, since the number of DSC element and differential quadrature element increases with the increase of number of stiffeners. It was illustrated by Meirovitch and Kwak [22] that the ability of satisfying the differential equation affected more on the convergence rate of the widely used FEM than the one of satisfying the natural boundary conditions. This indicates that a high-order FEM has better convergence rate than the lower order one. The weak form quadrature element method (QEM) [21,23] and the weak form Galerkin finite element method (WG-FEM) [24] may be the two candidates of relatively inexpensive yet accurate numerical methods. Due to high accuracy and convergence rate, weak form highorder elements such as weak form high-order quadrature element and spectral element, are renowned recently [23]. In this paper, the QEM is to be further developed to fulfill the challenging buckling analysis of eccentrically stiffened plates. The QEM was proposed by Striz et al. in 1997 [25]. During the development of the QEM, important contributions were made by Ying and Liu [26], Zhong and Yue [27] and Jin et al. [28]. The major difference of the QEM from the conventional high-order finite element is that it employs non-uniformly distributed nodes. Therefore, the QEM possesses exponential rate of convergence, since a three-dimensional quadrature element is essentially the same as the well-known spectral finite element (SFE) [29]. Up to now, the QEM was successfully used for buckling analysis of thin rectangular plate [21], shear deformable plate [30], planar frameworks [31], and for local buckling analysis of stiffened panels [32]. The main objective of this paper is to present a weak-form high-order quadrature element formulation for simultaneous local and global buckling analysis of eccentrically stiffened plates. The materials of the stiffened plate are isotropic, orthotropic or anisotropic. To raise the computational efficiency, the assemblage of the quadrature beam elements with the quadrature plate element is completed via the displacement relations instead of using coincident nodes. Therefore, the same quadrature element model can be used for buckling analysis of stiffened plates with different numbers of stiffeners located at arbitrary positions. A number of case studies are performed. Finally, some conclusions are drawn based on the results reported in this paper.

Fig. 1. An eccentrically stiffened rectangular plate under uni-axial compression.

Fig. 2. Cross-section of stiffener: (a) I-shaped, (b) T shaped.

For simplicity in presentation, the blade shaped stiffener is denoted by I-shaped stiffener. 2.1. Energy expressions 2.1.1. Energy expressions of rectangular plate Denote u , v , w the middle plane displacement components in x, y and z directions. The strain energy of the entire rectangular plate (the skin) is given by

2. Quadrature element formulations of plate-stiffener system

U= Homogeneous isotropic and anisotropic materials are considered. Since the laminate composite materials can be equivalent to homogeneous anisotropic materials, and thus the formulations are also suitable for stiffened plates made of laminated composite materials. An eccentrically stiffened rectangular plate subjected to uniaxial compression in x direction is schematically shown in Fig. 1. The length, width and uniform thickness of the plate are denoted by a, b and t. A Cartesian coordinate system (x, y , z) is set at the middle plane of the plate and thus −t /2 ≤ z ≤ t /2. The origin of the coordinate system o is located at the plate center. The total number of stiffeners is denoted by N s and all stiffeners are parallel to x axis. The part of shaded area shown in Fig. 1 is commonly used for local buckling analysis [3,32]. In such a case, the unloaded edges are not free but elastically restrained (E). For the stiffened rectangular plate with one stiffener located at y = 0, only one unloaded edge is elastically restrained. Two types of cross-section, called the blade shaped and T shaped cross-section shown in Fig. 2, are considered, where t w and t f denote the thickness of the web and flange, h w is the height of the web, and l f represents the length of the flange, respectively.

155

1 2

a/2 b/2 t /2 {σ } T {ε } dzdydx −a/2 −b/2 −t /2

where the strain vector {ε } is defined as

{ε } =

(1)

⎧ ⎫ ⎨ εx ⎬ ⎩

⎧ ⎪ ⎨

εy = ⎩ γxy ⎭ ⎪

+z

∂u ∂x ∂v ∂y ∂u ∂y

⎧ ⎫ ⎨ w xx ⎬ ⎩

+

⎫ ⎪ ⎬

⎧ ⎪ ⎪ ⎨

∂2 w ∂ x2 ∂2 w +z ∂ y2 ⎪ ⎪ ⎪ ∂v ⎭ ⎩ 2 ∂2 w ∂x ∂ x∂ y

⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭

=

w yy = {ε0 } + z {κ } ⎭ 2w xy

⎧ ⎨ ⎩

ux vy

u y + vx

⎫ ⎬ ⎭

(2)

In Eq. (2), {ε0 } is the middle-plane strain vector and {κ } is the curvature vector. For simplifying in presentation, subscripts x and y on the displacement components are used to represent the differentiation of the displacement with respect to x and y. The stress vector is given by

{σ } =

⎧ ⎫ ⎨ σx ⎬ ⎩

σ y = [Q ] {ε} τxy ⎭

where [ Q ] is the stiffness matrix of the material.

(3)

156

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

Since the skin is relatively thin, thus plane stress is assumed. For orthotropic materials with principal material axes (1,2,3) aligned with the x, y and z axes, elements in matrix [ Q ] are given by

Q 11 = E 11 /(1 − μ12 μ21 ) Q 22 = E 22 /(1 − μ12 μ21 ) Q 13 = Q 31 = Q 23 = Q 32 = 0 Q 12 = Q 21 = μ21 Q 11 = μ12 Q 22 Q 33 = G 12

(4)

where E 11 , E 22 are modulus of elasticity in x and y direction, G 12 is the in-plane shear modulus, and μ12 , μ21 are major and minor Poisson’s ratios, respectively. Integrating Eq. (1) along the thickness direction yields

U=

a/2 b/2

1 2

−a/2 −b/2

{ε0 } {κ }

T

[ A] [B] [B] [D]



 {ε0 } dydx {κ }

(5)

V =

a/2 b/2 t /2

2

wx wy

−a/2 −b/2 −t /2

T

σx τxy τxy σ y



wx wy

 dzdydx

(6)

where σx , σ y are normal stresses in x and y directions, and τxy is the in-plane shear stress. Note that their positive directions for buckling analysis are opposite to the ones for stress and dynamic analysis. 2.1.2. Energy expressions of stiffeners Assume that the stiffened plate has N s stiffeners, which are parallel to the x axis and fully connected to the plate (skin) for simplicity in presentation. Each stiffener is modeled by a beam element. Besides considering the bending about x axis and the torsion, the axial deformation should also be considered for modeling eccentric stiffeners. The strain energy of the nth (n = 1, 2, ..., N s ) stiffener is given by

n

U =

1 2

a/2

⎡



⎣ E n An

−a/2

 n n

+G J

∂θ ∂x

∂ ub (x) ∂x

2

2

 + E n I ny

∂ 2 w (x, yn ) ∂ x2

2

u

v

w

wx

wy

Simply supported (S) (x = ∓a/2) Clamped (C ) (x = ∓a/2) Simply supported (S) ( y = ∓b/2) Free (F ) ( y = ∓b/2) elastically constrained (E) ( y = ∓b/2)

– 0 0 – –

0 0 – – –

0 0 0 – 0

– 0 – – –

– – – – –

The work done by the axial compressive stress stiffener is given by

n

V =

2



a/2

1

σx A

n

−a/2

∂ w (x, yn ) ∂x

2

 + rn2

(7)

where u b and w are the axial and transverse displacements of the neutral axis of the nth stiffener, w (x) = w (x, yn ), θ = w y | y = yn , A n is the cross-sectional area, I ny is the moment of inertia about the neutral axis paralleling to y axis, J n is the torsion constant, E n and G n are the elasticity and shear modulus of the nth stiffener, respectively. It is worth noting that only constant cross-section is considered.

∂θ ∂x

σx on the nth

2  dx

(8)

where rn2 is approximately calculated by [9]





(9)

In Eq. (9), I nz is the moment of inertia about the neutral axis paralleling to z axis, the eccentricity e is the distance between the neutral axis of the stiffener and the x axis, shown in Fig. 1. The second term in Eq. (8) is caused by the warping deformation, known as Wagner effect [34], and can be neglected since it is only a secondary effect. 2.1.3. Essential boundary conditions Since the QEM is to be used for solving the buckling problem of the eccentrically stiffened rectangular plate, only essential boundary conditions are required and given. For edge 1 and edge 3, i.e., x = ∓a/2, simply supported (S) and clamped boundary conditions are considered. For edge 2 and edge 4, i.e., ( y = ∓b/2), simply supported (S) and free (F ) boundary conditions are considered. For the partial plate model used for local buckling analysis [32], edge 2 and edge 4 are rotationally elastically constrained (E), namely, M y = G J w y /2, where G J is the torsion stiffness of the stiffener. A torsion bar element with the torsion stiffness G J /2 is added on each edge of the plate element to model the rotationally constraint. For simplifying in presentation, symbols SFSF denote that edge 1 and edge 3 are simply supported and the other two edges are free, where the edge numbers of the stiffened plate are defined in Fig. 1. The essential boundary conditions considered in this paper are summarized in Table 1. 2.2. Formulation of quadrature beam element Quadrature beam element with Gauss-Lobatto-Legendre (GLL) nodes, denoted by xi (i = 1, 2, ..., N ), is developed. GLL points are roots of following equation,

(1 − ξ 2 )

⎤ ⎦ dx

Boundary conditions

rn2 ≈ I ny + I nz + e 2 A n / A n

where [ A] is in-plane stiffness matrix, [D] is flexural stiffness matrix, and [B] is the coupling stiffness matrix, respectively. For the plate contains only one layer, [ A ] = [ Q ]t, [ B ] = [0] and [ D ] = t 3 [ Q ]/12. If the middle plane of the plate does not coincide with the neutral plane, such as the non-symmetric laminated plate, [ B ] = [0]. If the principal material axes (1,2,3) are not aligned with the x, y and z axes, coordinate transformations are necessary. For general laminated plates, formulas of computing [ A], [B] and [D] are given in [33]. The work done by the in-plane stresses applied on the plate edges is given by

1

Table 1 Essential boundary conditions.

d P N −1 (ξ ) dξ

=0

(10)

where P N −1 (ξ ) is the (N-1)th order Legendre polynomial. Since ξi ∈ [−1, 1], and thus xi = aξi /2 (i = 1, 2, ..., N ). A FORTRAN subroutine and its corresponding MATLAB function are available in [23] for calculations of the abscissas ξi and weights H i of GLL quadrature. According to criterion 5 for selection of displacement functions in developing a finite element, “The displacement functions must be so chosen that at the interelement boundaries, the displacements and their derivatives up to at least one order less than the highest-order derivative in the energy expression must be continuous. In other words, the displacement derivatives with the same

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

order must be finite” [35], different displacement functions can be assumed in formulating the N-node quadrature beam element. Since the stiffener is fully connected to the plate and thus the deflection of the plate at y = yn , i.e., w (x, yn ), is directly used as the deflection of the nth stiffener. The three displacements used in formulating the stiffness matrix are assumed as N 

u b (x) =

li (x)u b (xi ) =

N 

i =1

θ(x) =

N 

li (x)u bi

i =1

li (x)θ(xi ) =

i =1

N 

li (x)θi

(11)

i =1

w (x, yn ) =

N 

ϕi (x) w (xi , yn ) + ψ1 (x) w x (x1 , yn )

+ ψ N (x) w x (xN , yn ) N +2 

h i (x)(δi )n

i =1

The displacements used in formulating the geometric stiffness matrix are assumed as

⎧ N N   ⎪ ⎪ ⎪ w (x, yn ) = li (x) w (xi , yn ) = li (x)( w i )n ⎪ ⎨ i =1

i =1

(12)

N N ⎪   ⎪ ⎪ ⎪ li (x)θ(xi ) = li (x)θi ⎩ θ(x) = i =1

i =1

In Eq. (11) and Eq. (12), Lagrange interpolation function li (x) is defined as N 

li (x) =

x − xk

1

⎧ N ⎪ aE n A n  ⎪ ⎪ H l A li A lm ⎪ (kuim )n = ⎪ ⎪ 2 ⎪ ⎪ l =1 ⎪ ⎪ N ⎨ aE n I ny  w (kim )n = H l B˜ li B˜ lm ⎪ 2 ⎪ ⎪ l = 1 ⎪ ⎪ N ⎪ ⎪ aG n J n  ⎪ θ ⎪ ( k ) = H l A li A lm ⎪ n ⎩ im 2

 ∂ li (x)  A li = ∂ x x=xl ⎧ N % N   ⎪ ⎪ ⎪ (xl − xk ) (xi − xk ) (l = i ) ⎪ ⎪ ⎪ ⎪ k =1 ⎨ k =1 k=l,i k=i = N ⎪  1 ⎪ ⎪ ⎪ (l = i ) ⎪ ⎪ (xl − xk ) ⎪ ⎩



 N  1  = dx x=xi (xi − xk )

dli (x) 

B˜ l( N +2) =

(15)

2

+

1 2

T n

! u"  k

n

! θ"

{θ}nT k

n

ub

{θ}n

n

+

1 2

B˜ li =

!

1

(x1 − xN ) " + 2δl1 (xN − x1 ) " + 2δlN

(x1 − xN ) × B˜ l( N +1) 1

(xN − x1 ) × B˜ l( N +2)

B l1 (xl − x1 )(xl − x N ) + 2 Al1 (2xl − x1 − x N ) (19)

!

1

1

B˜ l1 =

B˜ lN =

It is clear that different w (x, yn ) are used in formulating the stiffness matrix and the geometric stiffness matrix. In other words, inconsistent geometric stiffness matrix is formulated. Numerical experience shows that the difference in buckling stress is negligible whether the consistent or inconsistent geometric stiffness matrix is used in the analysis since N is usually larger than 11. Substituting Eq. (11) into Eq. (7) and integrating the strain energy of the nth (n = 1, 2, ..., N s ) stiffener by GLL quadrature yields

ub

(l, i = 1, 2, ..., N )

(18)

li (x)(x − xi )(x − x N −i +1 )

k =1 k=i

Un =

(17)

k =1 k=l

In Eq. (14), li (x) is Lagrange interpolation function given by Eq. (13) and A ii is computed by

1

(i , m = 1, 2, ..., N + 2)

In Eq. (17), H l is the weight of GLL quadrature corresponding to the abscissa xl , Ali and Alm are the weighting coefficients of the first order derivative of Lagrange interpolation functions li (x) and lm (x) with respect to x at abscissa xl , B˜ li and B˜ lm are the weighting coefficients of the second order derivative of Hermite interpolation functions h i (x) and hm (x) with respect to x at abscissa xl . They can be explicitly calculated by the differential quadrature (DQ) law. For example, A li can be calculated by [21,23]

B˜ l( N +1) =

( i = 1, N )  1 1 l ( x )( x − x ) − A + ψi (x) i N − i + 1 ii x i − x N − i +1 x i − x N − i +1 ( i = 1 , N ) ϕi (x) = (14) 1 ⎪ li (x)(x − x1 )(x − x N ) ⎪ ( x − x )( x − x ) ⎪ 1 N i i ⎩ ( i = 2, 3, · · · , N − 1) x i − x N − i +1

⎧ ⎪ ⎪ ⎪ ⎨

A ii =

n

and B˜ li can be explicitly computed by [21,23]

and Hermite interpolation function h i (x) is defined by [21,23]

ψi (x) =

ub

(13)

xi − xk

k =1 k=i

$

= u b1 , u b2 , ..., u bN nT , {δ}n =  w 1 , w 2 , ..., w N , ( w x )1 , w )n and ( w x )N nT and {θ}n = θ1 , θ2 , ..., θ N nT , elements (kuim )n , (kim (kθim )n in axial, lateral and torsion stiffness matrices [ku ]n , [k w ]n and [kθ ]n can be explicitly given as where

l =1

i =1

=

#

157

B lN (xl − x1 )(xl − x N ) + 2 AlN (2xl − x1 − x N )



[B l1 (xl − x N ) + 2 Al1 ] − A 11 +

1

(20)

(x1 − xN )

[B lN (xl − x1 ) + 2 AlN ] − A N N +

(21)

1

(xN − x1 ) (22)

!

1

(xi − x1 )(xi − xN )

B li (xl − x1 )(xl − x N )

+ 2 Ali (2xl − x1 − xN ) + 2δli

"

( i = 2, 3, · · · , N − 1)

(23)

where δli is Kronecker symbol, A li is given by Eq. (18) and B li can be obtained by

B li =

 N  ∂ 2li (x)  = A lk A ki ∂ x2 x=xl

(24)

k =1

! " {δ}nT k w n {δ}n (16)

Weighting coefficients A lm and B˜ lm can be computed in a similar way. Simply replace subscript i in Eq. (18) to Eq. (24) by subscript m.

158

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

By the same token, substituting Eq. (12) into Eq. (8) and integrating the work done of the nth (n = 1, 2, ..., N s ) stiffener using GLL quadrature yields

Vn =

1 2

! " ! " 1 { w }nT g w n { w }n + {θ}nT g warp n {θ}n

(25)

2

where { w }n =  w 1 , w 2 , ..., w N nT and {θ}n = θ1 , θ2 , ..., θ N nT , elewarp w ments ( g im )n and ( g im )n in geometric stiffness matrices [ g w ]n warp and [ g ]n can be also explicitly computed by

⎧ N ⎪ aσx A n  ⎪ w ⎪ g im ) = H l A li A lm ( n ⎪ ⎨ 2 l =1

N ⎪ aσx A n rn2  ⎪ warp ⎪ ⎪ ( g ) = H l A li A lm ⎩ im n 2

(i , m = 1, 2, ..., N )

(26)

U=

=

M N ab  

8 1 2

j =1 i =1

v (x, y ) =

li (x)

M 

l j ( y )u ( xi , y j ) =

N 

li (x)

M 

i =1

j =1

i =1

j =1

N 

M 

N 

M 

li (x)

i =1

w 1 (x, y ) =

N 

l j ( y ) v ( xi , y j ) =

j =1

li (x)

i =1

M 

li (x)

i =1

l j ( y ) w ( xi , y j ) =

j =1

N  i =1

l j ( y )u i j

(27)

l j ( y)v i j

(28)

M 

l j ( y) w i j

j =1

(29) w 2 (x, y ) =

N +2 

h i (x)

i =1

w 3 (x, y ) =

N  i =1

M 

l j ( y) w i j

(30)

˜ ij h j ( y) w

(31)

j =1

li (x)

M +2  j =1

where Lagrange interpolation function and Hermite interpolation function are given by Eq. (13) and Eq. (14), respectively. In



{ε0 (xi , y j )} {κ (xi , y j )}



(32)

where {q} T = u 11 , u 21 , ....u N M , v 11 , v 21 , .... v N M , w 11 , w 21 , .... w N M , ..., ( w y ) N M , which contains 3M × N + 2( M + N ) degrees of freedom (DOFs), H i and H j are the weights of GLL quadrature, xi = aξi /2 (i = 1, 2, ..., N ), and y j = bη j /2 ( j = 1, 2, ..., M ), here η j are also GLL points. In Eq. (32), the strain matrix at an integration point (x j , y i ) can be explicitly given by using the DQ law as

{ε0 (xi , y j )} =

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

N  M 

A ilx δ jk ulk

l =1 k =1

M N  

y

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬

δil A jk v lk ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ l = 1 k = 1 ⎪ ⎪ ⎪ ⎪ N  M M N  ⎪ ⎪   ⎪ ⎪ y ⎪ ⎪ x ⎪ ⎪ δ A u + A δ v ⎪ ⎪ il lk jk lk il jk ⎩ ⎭ l =1 k =1

(33)

and

⎧ ⎫ N +2  M  ⎪ ⎪ ⎪ ⎪ x ⎪ ˜ ¯ lk ⎪ B il δ jk w ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ l =1 k =1 ⎪ ⎪ ⎪ ⎪ N + 2 M ⎨ ⎬ y ¯ lk {κ (xi , y j )} = δil B˜ jk w ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ l =1 k =1 ⎪ ⎪ ⎪ ⎪ M N ⎪ ⎪   ⎪ ⎪ y ⎪ ⎪ x ⎪ ⎪ ¯ A A w ⎪ ⎪ lk il jk ⎩ ⎭ l =1 k =1

j =1

li (x)

[ A] [B ] [B] [D]

(i = 1, 2, ..., N , j = 1, 2, ..., M )

Let the numbers of the node in x and y directions be denoted by N and M. An N × M-node quadrature rectangular plate element with GLL nodes is formulated. According to criteria 5 for selection of displacement functions [35], five different displacement functions are assumed for formulating the rectangular plate element. The five displacement functions are N 

T

l =1 k =1

2.3. Formulation of quadrature rectangular plate element

u (x, y ) =

{ε0 (xi , y j )} {κ (xi , y j )}

{q}T [k] {q}

l =1

In Eq. (26), Ali and Alm are the weighting coefficients of the first order derivative of Lagrange interpolation functions li (x) and lm (x) with respect to x at abscissa xl and can be explicitly calculated by Eq. (18).

H j Hi

Eq. (30), w i j contains the nodal deflection w i j (i = 1, 2, . . . , N, j = 1, 2, . . . , M) as well as the first order derivative with respect to x at nodes on edges x = ∓a/2, i.e., ( w x )1 j and ( w x ) N j ˜ i j contains the nodal deflection w i j ( j = 1, 2, . . . , M). In Eq. (31), w (i = 1, 2, . . . , N, j = 1, 2, . . . , M) and the first order derivative with respect to y at nodes on edges y = ∓b/2, i.e., ( w y )i1 and ( w y )iM (i = 1, 2, . . . , N). The assumed deflection function w 1 (x, y ) is used to compute w xy in Eq. (5) as well as w x , w y in Eq. (6), w 2 (x, y ) is used to compute w xx , and w 3 (x, y ) is used to compute w y y . In this way, a C 1 compatible plate element can be easily formulated. Substituting the assumed displacement functions into Eq. (5) and integrating the strain energy of the rectangular plate by GLL quadrature yields

(i = 1, 2, ..., N , j = 1, 2, ..., M )

(34)

¯ lk , contains the deflecwhere δ jk and δil are Kronecker symbols, w tion at all nodes and the first order derivative with respect to x at nodes on edges x = ∓a/2 as well as the first order derivative with respect to y at nodes on edges y = ∓b/2, i.e., w lk , ( w x )1k , ( w x )Nk , ( w y )l1 and ( w y )lM , (l = 1, 2, ..., N , k = 1, 2, ..., M ), A ilx and

A jk are computed by Eq. (18), B˜ ilx and B˜ jk are computed by Eq. (19) to Eq. (23), the superscript x on the weighting coefficients means that the differentiation is with respect to x, and the superscript y denotes that the differentiation is with respect to y. The quadrature plate element contains 3N × M + 2( N + M ) degrees of freedom (DOFs), three DOFs at all inner nodes, four DOFs at all boundary nodes, and five DOFs at four corner nodes, respectively. It is seen that using three deflection functions circumvents the difficulty in formulation of a C 1 compatible quadrature plate element without having the DOF of w xy at four corner points. The geometric stiffness matrix [g] of the plate element is also obtained by using GLL quadrature and the DQ law. The displacement given by Eq. (29), i.e., w 1 (x, y ), is used and thus inconsistent geometric stiffness matrix is formulated. Substituting Eq. (29) into Eq. (6) and performing the numerical integration result y

V =

=

y

M N ab  

8 1 2

H j Hi

j =1 i =1

{q}T [g] {q}

w x ( xi , y j ) w y ( xi , y j )

T

σxt τxy t τxy t σ y t



w x ( xi , y j ) w y ( xi , y j )



(35)

where stresses are assumed uniformly distributed across the thickness for simplicity in presentation.

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166



w x ( xi , y j ) w y ( xi , y j )

 =

⎧ N M ⎫  ⎪ ⎪ ⎪ x ⎪ ⎪ A il δ jk w lk ⎪ ⎪ ⎪ ⎨ ⎬ ⎪ ⎪ ⎪ ⎪ ⎩

l =1 k =1 M N   l =1 k =1

(u i )n = u (xi , yn ) =

N 

159

ll (xi )

lm ( yn )ulm =

m =1

l =1

⎪ ⎪ y ⎪ δil A jk w lk ⎪ ⎭

M 

M 

lm ( yn )u im

m =1

(i = 1, 2, ..., N )

(i = 1, 2, ..., N , j = 1, 2, ..., M )

(36)

y

where A ilx and A jk are computed by Eq. (18). It is not difficult to see that the non-zero terms of the geometric stiffness matrix are only related to the DOFs of w lk (l = 1, 2, ..., N , k = 1, 2, ..., M ).

(42)

where ll (xi ) = δil has been applied since xi is the GLL point, lm ( yn ) can be computed by Eq. (13). In short, one has

{u }n = [ T n4 ] {u }

(43)

¯ } can be By the same token, the relation between {δ}n and { w obtained by using Eq. (30) and written as

2.4. Quadrature element formulation of the plate-stiffener system

¯} {δ}n = [ T n1 ] { w The nodes of stiffeners are usually not coincided with the ones of the plate element, and thus the relations between the DOFs of the stiffeners and the ones of the plate element should be found to formulate the matrix equation of the plate-stiffener system. Since the stiffener is fully connected with the plate, thus one has



b

x

u ( xi )

n

= {u (xi , yn )} + e { w x (xi , yn )} = {u }n + e [ A ] {δ}n (37)

where u is the middle plane displacement of the plate in x direction, w x is the first order derivative of the plate deflection with respect to x, and e is the eccentricity of the stiffener to the x axis, [ A x ] is the weighting coefficient matrix of the first order derivative of the deflection w with respect to x, respectively. Substituting Eq. (37) into Eq. (16) results

Un =

  {δ}nT [k w ]n + e 2 [ A x ]T [ku ]n [ A x ] {δ}n 2 2   & u ' 1 1 T x + {u }n e [k ]n [ A ] {δ}n + {δ}nT e [ A x ]T [ku ]n {u }n 1

{u }nT [ku ]n {u }n +

+

2 1 2

1

2

{θ}nT [kθ ]n {θ}n

(38)

¯ } can be obtained by using The relation between { w }n and { w Eq. (29) and written as

¯} { w }n = [ T n2 ] { w

Un =

2

¯} {θ}n = [ T n3 ] { w

1 2

1 2

{δ}nT [kˆ w ]n {δ}n +

{δ}nT [k w −u ]n {u }n +

1 2

1 2

{u }nT [ku − w ]n {δ}n

{θ}nT [kθ ]n {θ}n

For assemblage, the matrix equation of the plate-stiffener system is also written as a partitioned form, namely,

⎡ ˜ uu [k ] ⎣ [k˜ vu ] [k˜ wu ]

[k˜ uv ] [k˜ v v ] [k˜ w v ]

⎫ ⎫ ⎡ ⎤⎧ ⎤⎧ [kuu ] [kuv ] [ku w ] ⎨ {u } ⎬ [0] [0] [0] ⎨ {u } ⎬ ⎣ [k vu ] [k v v ] [k v w ] ⎦ { v } = σx ⎣ [0] [0] [0] ⎦ { v } ⎩ ¯ ⎭ ⎩ ¯ ⎭ [k wu ] [k w v ] [k w w ] {w} [0] [0] [ g w w ] {w} ⎡

(40)

{u } =

⎫ ⎫ ⎡ ⎤⎧ ⎤⎧ [k˜ u w ] ⎨ {u } ⎬ [0] [0] [0] ⎨ {u } ⎬ [k˜ v w ] ⎦ { v } = σx ⎣ [0] [0] [0] ⎦ { v } ⎩ ¯ ⎭ ⎩ ¯ ⎭ {w} [0] [0] [ g˜ w w ] {w} [k˜ w w ] (47)

where each sub-matrix in the stiffness matrix is obtained by

[k˜ uu ] = [kuu ] +

Ns  [ T n4 ]T [ku ]n [ T n4 ],

⎫ ⎧ u ⎪ ⎪ ⎪ ⎬ ⎨ 11 ⎪ u 21

... ⎪ ⎪ ⎪ ⎪ ⎭ ⎩ uNM ⎧ ⎫ ⎨ {u } ⎬ {q} = { v } ⎩ ¯ ⎭ {w}

,

{v } =

⎫ ⎧ v ⎪ ⎪ ⎪ ⎬ ⎨ 11 ⎪ v 21

... ⎪ ⎪ ⎪ ⎪ ⎭ ⎩ vNM

,

¯}= {w

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩

[k˜ u w ] = [ku w ] +

w 11 w 21

⎫ ⎪ ⎪ ⎬

... ⎪ ⎪ ⎭ ( w y )N M

Ns  [ T n4 ]T [ku − w ]n [ T n1 ],

[k˜ v w ] = [k v w ],

n =1

(39)

To illustrate the process of formulating the matrix equation of the plate-stiffener system, the matrix equation of the rectangular plate element is re-written as follows:

where

(46)

n =1

{u }nT [ku ]n {u }n +

+

(45)

¯ } can be obtained by using And the relation between {θ}n and { w Eq. (31) and written as

or in short

1

(44)

[k˜ wu ] = [k wu ] +

Ns  [ T n1 ]T [k w −u ]n [ T n4 ],

[k˜ uv ] = [kuv ],

n =1

[k˜ v v ] = [k v v ], [k˜ w v ] = [k w v ], N s    [ T n1 ]T [kˆ w ]n [ T n1 ] + [ T n3 ]T [kθ ]n [ T n3 ] . [k˜ w w ] = [k w w ] +

[k˜ vu ] = [k vu ],

n =1

The non-zero sub-matrix in geometric stiffness matrix is calculated as

[ g˜ w w ] = [ g w w ] +

 [ T n2 ]T [ g w ]n [ T n2 ] + [ T n3 ]T [ g warp ]n [ T n3 ]

Ns   n =1

,

(41)

The relations between the nodal DOFs of the beam element and the ones of plate element are summarized as follows. Using Eq. (27), elements in {u }n can be expressed in terms of elements in {u } as follows:

The assemblage of one rectangular plate element with N s beam elements is regarded as the N × M-node quadrature plate-stiffener element. It is seen that the matrix size of the quadrature platestiffener element, i.e., Eq. (47), is the same as the one of the quadrature rectangular plate element and independent of the number of the stiffeners. Thus the matrix size for buckling analysis of the eccentrically stiffened rectangular plate is much smaller as compared to the one of the conventional finite element method and the one of the DQEM [7,8], since only one N × M-node quadrature plate-stiffener element is used in the analysis.

160

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

2.5. Solution procedures Although Eq. (47) together with the applied essential boundary conditions can be directly solved by a generalized eigen-solver, however, it is computationally inefficient. Therefore, a different way is used herein. Rewrite Eq. (47) as



[k˜ in ]

[k˜ in− w ] in− w T ˜ [k ] [k˜ w w ]



{qin } ¯} {w

(



= σx

[0] [0] [0] [ g˜ w w ]



{qin } ¯} {w

 (48)

)

where {q in } T = {u } T , { v } T is the vector containing nodal DOFs of the in-plane displacements. Eliminating {q in } yields

¯ } = σx [ g˜ w w ]{ w ¯} [k¯ ]{ w ˜ww

(49)

where [k¯ ] = [k ] − [k ] [k ] [k ]. Note that essential boundary conditions for the in-plane displacements u and v, i.e., u (0, ∓b/2) = 0 and v (∓a/2, 0) = 0, have been applied to get the inverse of the in-plane stiffness matrix, although essential in-plane boundary conditions do not needed for free unloaded boundaries and loaded edges. Thus [k˜ in ] is changed to [k¯ in ]. Other essential boundary conditions for the in-plane displacements u and v can be also used to remove the three rigid body displacements. Applying the essential boundary conditions for the out-of-plane displacement w, Eq. (49) can be rewritten as



[k¯ ee ] [k¯ re ]

[k¯ er ] [k¯ rr ]



{we} {wr }

˜ in− w T



 = σx

¯ in −1

˜ in− w

[0] [0] [0] [ g˜ w w ]



{we} {wr }

 (50)

where { w e } contains non-zero derivative DOFs at boundary nodes to be eliminated, and { w r } contains non-zero displacement DOFs as well as non-zero derivative DOFs at unloaded boundary nodes to be retained. Note that due to considering the warping effect in the formulation of the geometric stiffness matrix, derivative DOFs at unloaded boundary nodes are non-zero and should be included in { w r }. After eliminating { w e }, Eq. (50) can be further rewritten as

[k∗ ]{ w r } = σx [ g˜ w w ]{ w r }

(51)

where [k∗ ] = [k¯ rr ] − [k¯ re ][k¯ ee ]−1 [k¯ er ]. Solving Eq. (51) by a generalized eigen-solver yields the eigenvalues. The lowest eigen-value is the critical buckling load of the eccentrically stiffened rectangular plate. Since the matrix size of Eq. (51) is much smaller than the one of Eq. (47), thus solving Eq. (51) is computationally more efficient than solving Eq. (47). 3. Numerical results and discussion A FORTRAN program is developed. To verify the formulations and program, the buckling stress of eccentrically stiffened rectangular plates subjected to uni-axial compression, schematically shown in Fig. 1, is analyzed. For comparison purposes, the dimensions of all stiffeners are assumed the same. A number of cases with different cross-section shapes and numbers of stiffeners are analyzed by the proposed quadrature plate-stiffener element. It is worth noting that only one N × M-node quadrature plate-stiffener element is used in the analysis and thus the computational effort should be relatively small as compared with the one of the existing DQEM [7,8] and FEM [9,10]. Case 1: Consider first the buckling of an SSSS stiffened rectangular plate with one I-shaped stiffener located at y = 0. The plate dimensions are a = 222.30 mm, b = 111.15 mm, and t = 1 mm. The width (t w ) and height (h w ) of the I-shaped stiffener are 3.0 mm and 7.41 mm, respectively. Isotropic material with elasticity modulus 10,000.0 MPa and Poisson’s ratio μ = 0.3 is considered.

Fig. 3. The critical buckling stress of the SSSS stiffened plate with one eccentric I-shaped stiffener at y = 0.

Therefore, γn = ( E n I ny )/(bD 11 ) ≈ 10.0 and δn = A n /(bt ) ≈ 0.2. Both eccentric and concentric stiffeners are studied. The eccentric stiffener refers to that the lower surface of the stiffener is connected with the upper surface of the plate, while the concentric stiffener refers to that the stiffener is placed with the neutral axis coincides with the one of the plate. The isotropic stiffened rectangular plate is subjected to uni-axial compression and the compressive stress σx is uniformly distributed over the cross-sections of the stiffener and plate at x = ±a/2. A convergence study is performed first. The number of element nodes varies from 9 to 21. Finite element results are included for verifications. The variations of critical buckling stress with N or M (N = M) are shown in Fig. 3. In the figure, notation FEM (beam) denotes that the stiffener is modeled by beam elements and notation FEM (plate) represents that the stiffener is modeled by plate elements. From Fig. 3, it is clearly seen that the rate of convergence of the proposed quadrature plate-stiffener element is high and one 9 × 9 quadrature plate-stiffener element can yield convergent buckling stress. The QEM data are slightly higher than the FEM results. Perhaps the reason is that the QEM does not take the shear deformation of the blade shaped stiffener into consideration. Besides, the developed quadrature beam element is not a full three-dimensional beam element. Note that the data obtained by FEM with beam or plate elements to model the stiffeners are also slightly different. Table 2 summarizes the first two buckling stresses of SSSS stiffened rectangular plates with one I-shaped stiffener. The QEM data are slightly closer to the ones of FEM (beam). If the stiffener is concentrated, the first two buckling modes are global buckling modes and the third one is the local buckling mode. It is seen that the global buckling stress is much lower than the one of eccentrically stiffened rectangular plate. It is interesting to see that the local buckling stress is almost the same for both eccentrically and concentrically stiffened rectangular plates. The first two buckling modes of eccentrically stiffened rectangular plate are shown in Fig. 4 and Fig. 5. They agree quite well with the mode shapes of the FEM. To obtain the local buckling stress, a half plate model is also solved by the QEM. The boundary conditions are SSSE or SESS since the unloaded edges of the plate do not have stiffeners and are simply supported. A half of the GJ is used for the rotationally elastically constrained boundary [32]. The critical buckling stress obtained using one 17 × 17 quadrature plate-stiffener element is 14.82 MPa, which is close to the one of FEM (plate). The corresponding mode shape is shown in Fig. 5c. The buckling mode also has five half waves along x direction and agrees with the one obtained by using a full plate model. It is seen that the SSSE (or SESS) partial plate model may yield conservative local buckling stress for the case considered. From data listed in Table 2 and Figs. 3–5, the formulations of the QEM and written program are validated. It is also worth noting that the warping of the stiffener does not have any effect on the critical buckling stress for this case.

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

161

Table 2 Comparisons of buckling stress of SSSS stiffened plate with one I-shaped stiffener. Mode No (Type)

QEM (eccentric)

QEM (concentric)

FEM (beam)

FEM (plate)

QEM (half plate model)

1 (global) 2 (local)

13.20 MPa 15.37 MPa

5.871 MPa 15.33 MPa (3)∗

13.07 MPa 15.02 MPa

13.02 MPa 14.95 MPa

– 14.82 MPa (1)∗

(.)∗ : mode number.

Fig. 4. The global buckling mode shape of the SSSS stiffened plate with one eccentric I-shaped stiffener at y = 0.

Case 2: Consider next the buckling of an SFSF stiffened rectangular plate with four I-shaped stiffeners equally spaced at y = ±b/2 and y = ±b/6. The plate dimensions are a = 700.0 mm, b = 840.0 mm, and t = 1 mm. The width (t w ) and height (h w ) of the I-shaped stiffener are 3.0 mm and 20.0 mm, respectively. For comparisons, isotropic material with elasticity modulus E n = 72, 000 MPa and Poisson’s ratio μ = 0.33 is considered. Therefore, γn = ( E n I ny )/(bD 11 ) = 25.46 and δn = A n /(bt ) ≈ 0.0714. Both eccentric and concentric I-shaped stiffeners are studied. The stiffened rectangular plate is subjected to uni-axial compression and the compressive stress σx is uniformly distributed over the crosssections of the I-shaped stiffener and plate at x = ±a/2. Table 3 summarizes the critical buckling stress of the SFSF stiffened rectangular plate with four equally spaced I-shaped stiffeners. The existing results obtained by Ritz and Kantorovich methods are included for comparisons. Note that a partial plate model, i.e., the shaded area shown in Fig. 1, was used in using Ritz and Kantorovich methods and thus the boundary conditions are SESE. The QEM data with the full plate model are closer to the Ritz results. Again whether the stiffeners are eccentrically or concentrically located affects the buckling stress only slightly, since the first buckling mode is a local buckling one. The SESE partial plate is also analyzed by using the QEM and FEM, the results of QEM and FEM are almost the same but much smaller than the ones obtained by the QEM with full plate models. As demonstrated in Case 1, the SESE partial plate model may yield conservative local buckling stress for the case considered. The first buckling mode of the eccentrically stiffened rectangular plate with four I-shaped stiffeners, obtained by the QEM, is shown in Fig. 6. Fig. 6a shows the mode shape with a full plate model and Fig. 6b is the mode shape with a partial plate model. It is seen that the two mode shapes agree each other.

To demonstrate the effect of aspect ratio on the local buckling stress, the SFSF stiffened rectangular plate with four equally spaced I-shaped stiffener and different aspect ratios are analyzed by the newly developed method. Full plate model is used. The results are listed in Table 4. It is seen that the aspect ratio affects slightly the magnitude of the local buckling stress. On the other hand, the aspect ratio has an appreciable effect on the buckling mode shape of the stiffened plate. The larger the aspect ratio is, the more number of half waves along x direction. Case 3: Consider the buckling of an SFSF stiffened rectangular plate with four T shaped stiffeners, which are equally spaced at y = ±b/2 and y = ±b/6. The plate dimensions are the same as the ones in Case 2. The thickness of the web and flange is t w = t f = 3 mm, the height of the web is h w = 20 mm, and the length of the flange is l f = 20 mm. Isotropic material with elasticity modulus E n = 72, 000 MPa and Poisson’s ratio μ = 0.33 is considered. Both eccentric and concentric stiffeners are studied. The stiffened rectangular plate is subjected to uni-axial compression and the compressive stress σx is uniformly distributed over the cross-sections of the T shaped stiffeners and plate at x = ±a/2. To investigate if the number of stiffeners and boundary conditions affect the rate of convergence of the proposed element, a convergence study is performed. The number of elemental nodes varies from 9 to 21. Finite element results as well as the data obtained by Ritz and Kantorovich methods are also included for verifications. The variations of critical buckling stress with N or M (N = M) are shown in Fig. 7. Again notation FEM (beam) means that the stiffener is modeled by beam elements. From Fig. 7, it is clearly seen that the rate of convergence of the proposed QEM is also high and one 13 × 13 quadrature plate-stiffener element can yield convergent buckling stress. Com-

162

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

Fig. 5. The local buckling mode shape of the SSSS stiffened plate with one eccentric I-shaped stiffener at y = 0. Table 3 Comparisons of local buckling stress of SFSF stiffened plate with four I-shaped stiffeners. QEM (eccentric)

QEM (concentric)

Ritz method (SESE) [3]

Kantorovich method (SESE) [3]

QEM (SESE) [32]

FEM (SESE) [32]

5.44 MPa (3)∗

5.43 MPa (3)

5.47 MPa (3)

5.23 MPa (3)

4.99 MPa (3)

5.00 MPa (3)

(.)∗ : Number of half waves along x direction.

parisons with Fig. 3 reveal that the rate of convergence is a little lower and thus the number of stiffeners and boundary conditions may affect the rate of convergence slightly. Numerical experience shows that the warping of the stiffener has only a little effect on the critical buckling stress for the problem considered. Table 5 lists the critical buckling stresses of SFSF stiffened rectangular plates with four equally spaced T shaped stiffeners. The existing results obtained by Ritz and Kantorovich methods using partial plate model and finite element result with full plate model are included for comparisons. The QEM buckling stress with full plate model is closer to the one of the FEM. Again whether the stiffeners are eccentrically or concentrically located affects the buckling stress only slightly, since the first buckling mode is local buckling one. The SESE partial plate is also analyzed by using the QEM, the result of partial plate model is much lower than the one with full plate model. This indicates that the SESE partial plate model may yield much conservative buckling stress for the stiffened rectangular plate with T shaped stiffeners, too. The buckling mode of eccentrically stiffened rectangular plate is shown in Fig. 8. It is seen that the mode shape of the QEM

with partial plate model is different from the one of FEM with a full plate model. Both QEM and Kantorovich method with a partial plate model yield a mode shape having only 3 half waves along x direction, while the FEM and QEM with the full plate model and Ritz method with a partial plate model yield mode shapes containing 4 half waves. Perhaps this is the reason why the buckling stress obtained by QEM and Kantorovich method using a partial plate model is much lower than the one using a full plate model. The partial plate model is fully equivalent to the full mode only if the stiffened plate contains infinite repeated partial plate model in the y direction. But in reality the stiffened plate can only be divided into limited numbers (3 for this case) of the partial plate model. As a result, the edge effects at y = ±b/2 are neglected in the partial plate model. Therefore, care should be taken if a partial plate model is to be used to get the local buckling behavior, since the mode shape may be different from the real one. It should be mentioned that the second QEM mode shape of the SESE plate, not shown in this paper to save the space, contains 4 half waves along x direction and the corresponding buckling stress is 5.49

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

163

Fig. 6. Local buckling mode of SFSF eccentrically stiffened plate with four I-shaped stiffeners. Table 4 Effect of aspect ratio on the local buckling stress of SFSF stiffened plate with four I-shaped stiffeners. a/b

0.5

0.75

5/6

1.0

2.0

Buckling stress

5.404 MPa (2)∗

5.635 MPa (3)

5.438 MPa (3)

5.445 MPa (4)

5.470 MPa (8)

(.)∗ : Number of half waves along x direction. Table 5 Comparisons of local buckling stress of SFSF stiffened plate with four T shaped stiffeners. QEM (eccentric)

QEM (concentric)

Ritz method (SESE) [3]

Kantorovich method (SESE) [3]

QEM (SESE) [32]

FEM (eccentric)

5.824 MPa (4)∗

5.815 MPa (4)

5.91 MPa (4)

5.65 MPa (3)

5.40 MPa (3)

5.78 MPa (4)

(.)*: Number of half waves along x direction.

MPa, much closer to 5.78 MPa obtained by FEM with a full plate model. It is seen that the superiority and advantage of the newly developed QEM against the commonly used numerical techniques, including weak form Ritz method and strong form Kantorovich methods, is that the method has good flexibility in the buckling analysis of plate-stiffener system with an arbitrary pattern. Due to its weak form nature, applying various boundary conditions is easy. In addition, the computational effort is similar to the Ritz and Kantorovich methods since only one plate element is used in the analysis. Both Ritz and Kantorovich methods are used to analyze the local buckling based on the partial plate model. However, the partial plate model is an approximate one when the number of stiffeners is small. Besides, the stiffeners should be equally spaced to employ the partial plate model. Present QEM does not have such limitations. Case 4: Consider the buckling of an SSSS eccentrically stiffened rectangular plate with three T shaped stiffeners, which are equally spaced at y = ±b/4 and y = 0. The plate dimensions are a = 1800.0 mm, b = 3600.0 mm, and t = 21 mm. Two types of T cross-section shapes, called P50 and P51 [10] for simplicity in presentation, are considered. For P50, the thickness and height of the web are t w = 20 mm and h w = 50 mm, the thickness and length of the flange are t f = 30 mm and l f = 200 mm. For P51, the thickness and height of the web are t w = 12 mm and h w = 84 mm, the thickness and length of the flange are t f =

Fig. 7. The buckling stress of the SFSF stiffened plate with four eccentric T shaped stiffeners.

15 mm and l f = 100 mm. Isotropic material with elasticity modulus E n = 205, 800 MPa and Poisson’s ratio μ = 0.3 is considered. The stiffened rectangular plate is subjected to uni-axial compression and the compressive stress σx is uniformly distributed over the cross-sections of the T shaped stiffeners and plate at x = ±a/2. Table 6 lists the critical buckling stresses of SSSS eccentrically stiffened rectangular plates with three equally spaced T shaped stiffeners. Finite element results with full plate model are included for comparisons. The QEM buckling stress agrees more closely to the one of the FEM using beam elements to model the stiffeners. The QEM buckling stress with considering warping effect is slightly lower than the ones without considering warping effect. The SESE partial plate is also analyzed by using the QEM. The buckling stress of the stiffened plate with P51 stiffeners is 421.76

164

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

Fig. 8. The local mode shape of the SFSF stiffened plate with four eccentric T shaped stiffeners. Table 6 Comparisons of buckling stress of SSSS eccentrically stiffened rectangular plate with three equally spaced T shaped stiffeners. Type

QEM (with warping effect)

QEM (without warping effect)

FEM [10]

FEM (beam)

FEM (plate)

P50 P51

491.24 MPa 410.80 MPa

495.71 MPa 411.98 MPa

494.0 MPa 410.0 MPa

491.4 MPa 411.0 MPa

487.2 MPa 408.0 MPa

MPa, slightly larger than the finite element one (411.0 MPa) based on a full plate model. Although the buckling mode-shape of the partial plate model is the same as the one of full plate model, the SESE partial plate model yields slightly higher buckling stress than the one of full plate model, quite different from previous cases. In other words, employing a partial plate model may not always yield conservative buckling stress. Therefore, it is desirable to get the buckling stress using a full plate model during design stage. From data listed in Table 6, one may conclude that the cross section area of the T stiffener has an appreciable effect on the buckling load of the stiffened plate. The larger the area of the stiffener is, the higher the buckling load. This is obvious that more loads can be carried by the stiffeners with larger cross section area if the dimensions of the plate remain the same. Besides, the area of the web and flange of the T stiffener might also have different effect on the buckling load. The research is beyond the scope of this paper and thus is not involved in present investigation. Case 5: Consider last the buckling of eccentrically stiffened laminated composite plates with either I-shaped or T shaped stiffeners. The stiffeners are equally spaced along y direction. Three types of eccentrically stiffened laminated composite plates, denoted by U90, U0 and F45, are investigated. The U90 and U0 stiffened laminated plates are made of carbon/epoxy unidirectional laminas and

the F45 stiffened laminated plate is made of carbon/epoxy fabrics. The material properties of carbon/epoxy unidirectional lamina are E 11 = 113, 000 MPa, E 22 = 9, 000 MPa, G 12 = 3, 820 MPa and major Poisson’s ratio μ12 = 0.302. The material properties of carbon/epoxy fabric are G 12 = 2, 350 MPa, E 11 = E 22 = 52, 000 MPa, and Poisson’s ratio μ12 = 0.048. The lay-up of the skin plate in U90 stiffened laminated composite plate is [02 902 ]s and the one of the flange and/or web is [(45/ − 45)4 06 ]s . The lay-up of the skin plate in U0 stiffened laminated plates is [0/45/ − 45/0]s and the one of the flange and/or web is [0/45/ − 45/0]s3 . The lay-up of the skin plate in F45 stiffened laminated plates is [45/ − 45]s and the one of the flange and/or web is [45/ − 45]s3 . For the U90 stiffened laminated composite plate, t = 1.2 mm, t w = t f = 4.2 mm, l f = 20 mm and h w = 17.9 mm. For the U0 stiffened laminated composite plate, t = 1.2 mm, t w = t f = 3.6 mm, l f = 20 mm and h w = 18.2 mm. For the F45 stiffened laminated composite plate, t = 1.31 mm, t w = t f = 3.93 mm, l f = 20 mm and h w = 18.035 mm. The plate dimensions are a = 700.0 mm, b = 840.0 mm for the stiffened plate with four stiffeners and b = 560.0 mm for the stiffened plate with five stiffeners. The eccentrically stiffened laminated composite plate is subjected to uni-axial compression and the compressive stress σx is

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

165

Table 7 Comparisons of local buckling stress of SFSF eccentrically stiffened laminated plates. Stiffened plate ID

Stiffener ID

QEM

Ritz method (SESE) [3]

Kantorovich method (SESE) [3]

ABAQUS

U90

4I 4T

3.65 MPa (2)∗ 3.73 MPa (3)

3.63 MPa (2) 3.77 MPa (3)

3.50 MPa (2) 3.72 MPa (3)

3.60 MPa (2) 3.72 MPa (3)

F45

4I 4T 5T

5.49 MPa (4) 5.62 MPa (4) 22.73 MPa (7)

5.91 MPa (3) 6.02 MPa (4) 24.32 MPa (7)

5.48 MPa (4) 5.59 MPa (4) 22.58 MPa (7)

5.23 MPa (4) 5.54 MPa (4) 22.06 MPa (7)

U0

4I 5I 4T 5T

3.99 MPa (2) 16.45 MPa (5) 4.20 MPa (3) 16.73 MPa (5)

4.13 MPa (2) 17.09 MPa (5) 4.43 MPa (3) 17.55 MPa (5)

3.94 MPa (2) 16.30 MPa (5) 4.24 MPa (3) 16.75 MPa (5)

4.00 MPa (2) 16.45 MPa (5) 4.19 MPa (3) 16.61 MPa (5)

(.)∗ : Number of half waves along x direction. Table 8 Local buckling stress of CFCF eccentrically stiffened laminated composite plates. Stiffened plate ID

Stiffener ID

U90

4I 4T

F45

U0

QEM

Stress ratio

α

ABAQUS

QEM

ABAQUS

4.27 MPa (2)∗ 4.44 MPa (2)

1.17 1.19

1.16 1.21

4.17 MPa (2) 4.51 MPa (2)

4I 4T 5T

5.79 MPa (3) 5.92 MPa (3) 23.01 MPa (7)

1.05 1.05 1.01

1.04 1.04 1.00

5.44 MPa (3) 5.78 MPa (3) 22.13 MPa (7)

4I 5I 4T 5T

4.53 MPa (2) 17.27 MPa (5) 4.76 MPa (2) 17.44 MPa (5)

1.14 1.05 1.13 1.04

1.12 1.05 1.11 1.03

4.49 MPa (2) 17.32 MPa (5) 4.64 MPa (2) 17.06 MPa (5)

(.)∗ : Number of half waves along x direction.

Fig. 9. Local buckling stress of the SFSF stiffened laminated plate with T shaped stiffeners.

uniformly distributed over the cross-sections of the I-shaped or T shaped stiffeners and plate at x = ±a/2. Two combinations of boundary condition, i.e., SFSF and CFCF, are considered. A convergence study is performed first for the F45 eccentrically stiffened laminated composite plates with four and five T shaped stiffeners. The combination of boundary condition is SFSF. The number of element nodes varies from 11 to 21. The existing results obtained by Ritz and Kantorovich methods using partial plate model are included for comparisons. The variations of critical buckling stress with N or M (N = M) are shown in Fig. 9. It is clearly seen that the rate of convergence of the proposed quadrature plate-stiffener element is also high and one 15 × 15 plate-stiffener element can yield convergent buckling stress. The number of stiffeners and material properties affect the rate of convergence of the QEM slightly. The QEM data are slightly closer to the ones obtained by the Kantorovich method for the case considered. Since the buckling mode contains seven half waves along x direction for the plate with five T shaped stiffeners, thus slightly higher number of nodes (N = M) plate-stiffener element is required to get the convergence results. Table 7 summarizes the critical buckling stresses of SFSF eccentrically stiffened laminated composite plates with four or five

equally spaced stiffeners. Eight cases are investigated. The existing results obtained by Ritz and Kantorovich methods using partial plate model and recalculated finite element data using ABAQUS with the full plate model are included for comparisons. It is seen that the QEM buckling stresses with full plate model are closer to the ones obtained by the Kantorovich method and ABAQUS. The data obtained by Ritz method are slightly larger than the ones obtained by other three methods. The predicted number of the half waves along x direction by QEM is also the same as the one reported in [3] and obtained by ABAQUS with the full plate model. Therefore, the developed QEM is further validated and can be used for buckling analysis of eccentrically stiffened laminated composite plates. In experiment, the loaded edges are usually clamped. Thus the buckling behavior of CFCF eccentrically stiffened laminated composite plates with four or five equally spaced stiffeners is analyzed by the developed technique. Results are listed in Table 8. For verifications, the results obtained by ABAQUS are also included. In the table, the stress ratio α is defined as

α = σcrC F C F /σcrS F S F

(52)

where σcrC F C F and σcrS F S F are critical buckling stress of CFCF and SFSF eccentrically stiffened laminated composite plates, respectively. In other words, the buckling load of the CFCF plate equals to the product of the ratio α and the known buckling load of the SFSF plate. It is seen that the QEM buckling stresses and the number of half waves along x direction agree with the ones obtained by the ABAQUS. The stress ratio, often used to get the buckling stress of SFSF eccentrically stiffened laminated composite plate from the test one of CFCF eccentrically stiffened laminated composite plate, is not a constant but varies case-by-cases. The stress ratios obtained by QEM and ABAQUS are quite close and vary in the range between 1.00 and 1.21. Besides the number of half waves along x direction of the CFCF eccentrically stiffened laminated composite

166

J. Deng et al. / Aerospace Science and Technology 87 (2019) 154–166

plate may be also different and usually smaller than the one of SFSF eccentrically stiffened laminated composite plate. Data listed in Table 8 are new and may be served as a reference for developing new numerical methods. 4. Conclusions In this paper, a weak-form high-order quadrature element formulation is presented for simultaneous local and global buckling analysis of eccentrically stiffened plates made of isotropic and laminated composite materials. The quadrature beam elements are connected with the quadrature plate element via the displacement relationship instead of using coincident nodes, and thus the stiffeners can be located anywhere other than along the nodal lines of the plate element. Besides, the resulting matrix size of the N × M-node quadrature plate-stiffener system is small and the same as the one of an N × M-node quadrature plate element, thus presented formulation could be also used for optimization routines in preliminary design of eccentrically stiffened plates. Convergence studies and comparisons have been made to show the rate of convergence and to validate the applicability of the presented formulation. It is found that the rate of convergence of the QEM is high, one 15 × 15-node quadrature plate-stiffener element can yield convergent solutions for eccentrically stiffened plates with different numbers and cross-sections of stiffeners. A number of case studies are presented. Results show good agreement with existing data and finite element solutions. Based on the reported results, one may conclude that presented formulation is simple, straightforward, and reliable, which can efficiently compute the local and global buckling stress of eccentrically stiffened plates made of isotropic or laminated composite materials. Conflict of interest statement The authors have no conflict of interest to declare. Acknowledgements Jian Deng, Xinwei Wang and Guangming Zhou acknowledge the support from the Priority Academic Program Development of Jiangsu Higher Education Institutions. References [1] Bedair, Analysis and limit state design of stiffened plates and shells: a world view, Appl. Mech. Rev. 62 (2009) 020801. [2] P.M. Vuong, N.D. Duc, Nonlinear response and buckling analysis of eccentrically stiffened FGM toroidal shell segments in thermal environment, Aerosp. Sci. Technol. 79 (2018) 383–398. [3] C. Bisagni, R. Vescovini, Analytical formulation for local buckling and postbuckling analysis of stiffened laminated panels, Thin-Walled Struct. 47 (3) (2009) 318–334. [4] S.P. Timoshenko, J.M. Gere, Theory of Elastic Stability, second edition, McGraw– Hill Book Company, New York, 1961. [5] D.G. Stamatelos, G.N. Labeas, K.I. Tserpes, Analytical calculation of local buckling and post-buckling behavior of isotropic and orthotropic stiffened panels, Thin-Walled Struct. 49 (2011) 422–430. [6] Z. Siddiqi, A. Kukreti, Analysis of eccentrically stiffened plates with mixed boundary conditions using differential quadrature method, Appl. Math. Model. 22 (1998) 251–275.

[7] C. Dai, Y. Wang, Differential quadrature element analysis for stability of stiffened plates, J. Nanjing Univ. Aeronaut. Astronaut. 39 (5) (2007) 642–645. [8] L. Jiang, Y. Wang, X. Wang, Buckling analysis of stiffened circular cylindrical panels using differential quadrature element method, Thin-Walled Struct. 46 (4) (2008) 390–398. [9] M. Guo, I.E. Harik, Stability of eccentrically stiffened plates, Thin-Walled Struct. 14 (1992) 1–20. [10] Hughes, B. Ghosh, Y. Chen, Improved prediction of simultaneous local and overall buckling of stiffened panels, Thin-Walled Struct. 42 (6) (2004) 827–856. [11] Z. Xia, F. Ju, K. Sasaki, A general finite element analysis method for balloon expandable stents based on repeated unit cell (RUC) model, Finite Elem. Anal. Des. 43 (2007) 649–658. [12] X. Wang, J.S. Hansena, D.C.D. Oguamanam, Layout optimization of stiffeners in stiffened composite plates with thermal residual stresses, Finite Elem. Anal. Des. 40 (2004) 1233–1257. [13] G.W. Wei, A new algorithm for solving some mechanical problems, Comput. Methods Appl. Mech. Eng. 190 (2001) 2017–2030. [14] Ö. Civalek, Three-dimensional vibration, buckling and bending analyses of thick rectangular plates based on discrete singular convolution method, Int. J. Mech. Sci. 49 (2007) 752–765. [15] Ö. Civalek, Buckling analysis of symmetric laminated composite plates by using discrete singular convolution, Trends Appl. Sci. Res. 2 (6) (2007) 460–471. [16] Ö. Civalek, Free vibration and buckling analyses of composite plates with straight-sided quadrilateral domain based on DSC approach, Finite Elem. Anal. Des. 43 (13) (2007) 1013–1022. [17] S.K. Lai, Y. Xiang, DSC analysis for buckling and vibration of rectangular plates with elastically restrained edges and linearly varying in-plane loading, Int. J. Struct. Stab. Dyn. 9 (3) (2009) 511–531. [18] S.K. Lai, L.H. Zhang, Thermal effect on vibration and buckling analysis of thin isotropic/orthotropic rectangular plates with crack defects, Eng. Struct. 177 (2018) 444–458. [19] Ö. Civalek, Application of differential quadrature (DQ) and harmonic differential quadrature (HDQ) for buckling analysis of thin isotropic plates and elastic columns, Eng. Struct. 26 (2004) 171–186. [20] F. Tornabene, N. Fantuzzi, F. Ubertini, E. Viola, Strong formulation finite element method based on differential quadrature: a survey, Appl. Mech. Rev. 67 (2015) 020801. [21] X. Wang, Differential Quadrature and Differential Quadrature Based Element Methods: Theory and Applications, Butterworth-Heinemann, Oxford, 2015. [22] L. Meirovitch, M.K. Kwak, Convergence of the classical Rayleigh–Ritz method and the finite element method, AIAA J. 28 (8) (1990) 1509–1516. [23] X. Wang, Z. Yuan, C. Jin, Weak form quadrature element method and its applications in science and engineering: a state-of-the-art review, Appl. Mech. Rev. 69 (2017) 030801. [24] L. Mu, J. Wang, G.W. Wei, X. Ye, S. Zhao, Weak Galerkin methods for second order elliptic interface problems, J. Comput. Phys. 250 (2013) 106–125. [25] A.G. Striz, W.L. Chen, C.W. Bert, Free vibration of plates by the high accuracy quadrature element method, J. Sound Vib. 202 (5) (1997) 689–702. [26] Y. Xing, B. Liu, High-accuracy differential quadrature finite element method and its application to free vibrations of thin plate with curvilinear domain, Int. J. Numer. Methods Eng. 80 (2009) 1718–1742. [27] H.Z. Zhong, Z.G. Yue, Analysis of thin plates by the weak form quadrature element method, Sci. China, Phys. Mech. Astron. 57 (8) (2012) 1–10. [28] C. Jin, X. Wang, L. Ge, Novel weak form quadrature element method with expanded chebyshev nodes, Appl. Math. Lett. 34 (2014) 51–59. [29] X. Wang, Z. Yuan, A novel weak form three-dimensional quadrature element solution for vibrations of elastic solids with different boundary conditions, Finite Elem. Anal. Des. 141 (2018) 70–83. [30] H. Zhong, C. Pan, H. Yu, Buckling analysis of shear deformable plates using the quadrature element method, Appl. Math. Model. 35 (10) (2011) 5059–5074. [31] H. Zhong, R. Zhang, H. Yu, Buckling analysis of planar frameworks using the quadrature element method, Int. J. Struct. Stab. Dyn. 11 (2) (2011) 363–378. [32] X. Wang, Z. Yuan, An efficient method for local buckling analysis of stiffened panels, Trans. Nanjing Univ. Aeronaut. Astronaut. 36 (1) (2019) 16–24. [33] J.R. Vinson, T.W. Chou, Composite Materials and Their Use in Structures, Applied Science Publishers LTD, London, 1975. [34] W.F. Chen, T. Atsuta, Theory of Beam-Columns, Volume 2-Space Behavior and Design, McGraw–Hill Book Co., New York, USA, 1977. [35] T.Y. Yang, Finite Element Structural Analysis, Prentice-Hall, Englewood Cliffs, New York, 1986.