Three-dimensional non-linear elasticity-based 3D cubic B-spline finite element shear buckling analysis of rectangular orthotropic FGM plates surrounded by elastic foundations

Accepted Manuscript Three-dimensional non-linear elasticity-based 3D cubic B-spline finite element shear buckling analysis of rectangular orthotropic FGM plates surrounded by elastic foundations M. Shariyat, K. Asemi PII: DOI: Reference:

S1359-8368(13)00540-4 http://dx.doi.org/10.1016/j.compositesb.2013.09.027 JCOMB 2652

To appear in:

Composites: Part B

Received Date: Revised Date: Accepted Date:

10 June 2013 5 September 2013 7 September 2013

Please cite this article as: Shariyat, M., Asemi, K., Three-dimensional non-linear elasticity-based 3D cubic B-spline finite element shear buckling analysis of rectangular orthotropic FGM plates surrounded by elastic foundations, Composites: Part B (2013), doi: http://dx.doi.org/10.1016/j.compositesb.2013.09.027

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

Three-dimensional non-linear elasticity-based 3D cubic Bspline finite element shear buckling analysis of rectangular orthotropic FGM plates surrounded by elastic foundations M. Shariyat†a,b, K. Asemic a

Faculty of Mechanical Engineering, K.N. Toosi University of Technology, Tehran 19991-43344, Iran. b Center of Excellence in Smart Structures and Dynamical Systems c Department of Mechanical Engineering, Amirkabir University of Technology, Tehran, Iran.

Abstract In the present paper, shear buckling analysis of the orthotropic heterogeneous FGM plates is investigated for the first time. Moreover, influence of the Winkler-type elastic foundation is considered. The material properties are assumed to have in-plane orthotropy and transverse heterogeneity. The most accurate approach, i.e., the three-dimensional elasticity is employed instead of using the approximate plate theories. In contrast to all of the available displacement-based buckling analyses that have employed C0-continuous commercial finite element codes or semi-analytical methods, present formulations are C2continuous due to using the proposed 3D cubic B-spline element. Results are derived based on principle of minimum potential energy and a non-linear finite element procedure utilizing a Galerkin-type 3D cubic B-spline solution algorithm. Buckling loads are detected based on a generalized geometric stiffness concept. In this regard, effects of both the prebuckling and buckling states are considered. To present a better imagination and more detailed discussions, details of the buckling mode shapes and the foundation interaction are discussed for plates with simply-supported edges. In addition, the more practical free and clamped edge conditions are also considered. Keywords: A. Plates; B. Anisotropy; B. Buckling; C. Finite element analysis; Orthotropic functionally graded cubic B-spline elements.

1 Introduction While using two-dimensional woven or particle-reinforced composite constructions may lead to high specific strengths, especially when the principal directions of the material properties are coincident with the load transfer path, using transversely-graded distributions of the material properties eliminates the jumps in the transverse distribution of the stress components due to discrete and sometimes sever changes in the material properties of the adjacent layers. Using benefits of both of the mentioned material distributions (in-plane orthotropy and transverse gradation) may lead to an optimal three-dimensional material properties distribution and consequently, an optimized design. Plate-type components have been used extensively in skin panels, underframes, and main structures of the ground, marine, space, and even aerospace vehicles in addition to various civil, industrial, and common components. Shear in-plane loads are sometimes among the common load cases. When the large unstiffened regions of the plate-type structures are vulnerable to sever shear loads, these structures are likely to buckle. Recently, various problems in solid mechanics have been studied wherein the elastic coefficients are no longer constant but they are functions of the position. Indeed, the idea of non-homogeneity in the elastic coefficients is not at all hypothetical, but more realistic [1]. For example, elastic properties in soil may †

Corresponding author, Professor, E-mail addresses: [email protected] and [email protected] Tel.: +98 9122727199; Fax: +98 21 88674748, zip code: 19991-43344.

1

vary considerably with position. Moreover, some phases of the functionally graded materials may have orthotropic characters. For example, in graded composite materials, graded regions are treated as series of perfectly bonded composite layers, each layer being assigned slightly different properties. Continuously varying the volume fraction of the fibers in the transverse direction may also lead to orthotropic functionally graded or heterogeneous orthotropic materials. The heterogeneity of the materials may also stem from effects of the humidity, manufacturing process, production techniques, surface and thermal polishing processes, effect of radiation, and operating under radiation and elevated temperatures [2,3]. Sometimes, both the gradation and orthotropy in the material properties may be designed deliberately. For example, using the electroslag remelting technique, a functionally graded steel may be manufactured by using even only one material [4]. Forming the resulted FGM component through a rolling process, induces directional (orthotropic) material properties as well. Some researches have been accomplished on bending and thermoelastic analyses of the inhomogeneous composite plates. In all of the mentioned researches, the principal directions of the material properties that denote the orthotropy directions, are adopted parallel/perpendicular to the edges of the rectangular [5] or circular [6] plate, panel [7], or cylindrical [8] or conical shells [3]. Some researchers have studied behaviors of rectangular plates fabricated from orthotropic functionally graded materials (FGMs) under mechanical loads, using the 3D elasticity theory. Pan [9] derived an exact solution for a simply supported rectangular FG anisotropic laminated plate using the pseudo-Stroh formalism extending Pagano’s solution to the FG plates. Wen et al. [10] presented a 3D analysis for isotropic and orthotropic FG plates with simply supported edges under static and dynamic loads. Exponent laws were utilized for description of variations of the volume fractions of the constituent phases of the FG plates. A comprehensive literature survey on works published so far on the FGM plates may be found in a paper by Jha et al. [11]. Almost in all analyses performed for the orthotropic FGM plates, exponential laws were used. Numerous researches may be found in literature on buckling and postbuckling analyses of composite/functionally graded plates with (Duc and Tung [12]) or without (Abrate [13], Zhao et al. [14], Oyekoya et al. [15], Lee et al. [16], and Thai and Choi [17]) elastic foundations. Shariyat [18-20] studied mechanical and thermo-electro-mechanical buckling of laminated composite plates with elastic or viscoelastic layers, using high-order or global-local plate theories. Shariyat [21] and Alipour and Shariyat [22,23] investigated buckling of complex FGM plates with or without elastic foundations. Majority of the available studies on buckling of flat plates undergoing mechanical loads have been accomplished under axial or biaxial compressive loads. Nevertheless, some few researchers have investigated buckling of plates under shear loads. Budiansky [24] studied shear buckling of the isotropic clamped rectangular plates. Smith et al. [25] presented a Rayleigh–Ritz method for local buckling analysis of the rectangular unilaterally restrained plates under pure shear, using the classical plate theory. Azikov [26] investigated buckling of composite plates subjected to combined in-plane compression and shear. Loughlan [27] examined effects of the bending-twist coupling on the shear buckling behavior of the laminated composite plates, using a finite strip procedure and the classical plate theory. Kim et al. [28] carried out a buckling and postbuckling analysis for composite plates under pure shear. Han et al. [29] studied postbuckling behavior of laminated composite plates under the combination of in-plane shear, compression, and lateral loading, using a Lagrangian finite element formulation and the first-order sheardeformation plate theory. Lopatin and Korbut [30] studied buckling of clamped orthotropic laminated plates under shear loads, based on the classical plate theory. Papadopoulos and Kassapoglou [31] determined buckling loads of rectangular composite plates consisting of concentric rectangular layups under shear, using the classical plate theory, Ritz energy minimization, and a Navier-type solution. A semi-analytical approach to the buckling analysis of generally supported laminated plates subjected to various combinations of in-plane shear, compression, and tension loads was presented by Shufrin et al. for homogeneous isotropic[32] and orthotropic [33] plates, based on the principle of virtual work, classical plate theory, and a multi-term extended Kantorovich semi-analytical solution. Wu et al. [34] employed the extended spline collocation method (SCM) in the form of a cubic table to analyze shear buckling of thin isotropic homogeneous rectangular plates based on the classical plate theory. The linear buckling loads were determined through an eigenvalue analysis. Explicit expressions were developed by Tarján et 2

al. [35] for buckling analysis of long rectangular plates with constrained edges under shear and linearly varying edge loads, using Rayleigh-Ritz method. Wu and Nie [36] studied shear buckling of the steelconcrete composite slabs considering influence of the interface shear stiffness. Buckling and fracture collapse mechanisms of an elastic rectangular thin-plate with a crack under shear loading were analyzed by Brighenti and Carpinteri [37], for different boundary conditions. Chirica and Beznea [38] analyzed influence of the initial delaminations on buckling behavior of composite plates under shear and compression loadings numerically and experimentally. An analytical solution was developed by Yang et al. [39] for shear buckling of thin anisotropic rectangular plates having arbitrary boundary conditions along their edges or at their corners. Very rare buckling analyses have been presented so far based on the general three-dimensional theory of elasticity [40] or based on a three-dimensional elasticity with a uniform prebuckling stress assumption [41] which is equivalent to the membrane assumption used in the plate theories. Uymaz and Aydogdu [42] studied shear buckling of the rectangular functionally graded plates. The analysis was performed based on the linear strain energy concept and Ritz method. Asemi et al. [43] employed a threedimensional finite element elasticity approach to investigate buckling of heterogeneous functionally graded plates under biaxial compression, shear, tension-compression, and shear-compression load conditions. The foregoing brief review reveals that shear buckling analysis of the orthotropic functionally graded plates has not been performed so far, especially for plates resting on or surrounded by elastic foundations. In the present paper, nonlinear three-dimensional shear buckling of general orthotropic functionally graded plates surrounded by elastic foundations is investigated based on the elasticity approach instead of employing the approximate plate theories. The material properties are assumed to have in-plane orthotropy and transverse heterogeneity. The governing equations are derived using principle of minimum total potential energy and solved based on a Galerkin-type cubic B-spline finite element procedure. Therefore, in contrast to the finite element and semi-analytical procedures used for plates and shells, the displacement field is C2-continuous rather than C0-contiuous. Thus, in addition to satisfying the stress continuity conditions, accuracy of the results is significantly higher. Buckling loads are detected based on a geometric stiffness concept. In this regard, effects of both the prebuckling and buckling states are considered. To present a better imagination and more detailed discussions, details of the mode shapes and the foundation effects are discussed. Finally, a comprehensive parametric study is performed to extract more practical conclusions. 2 Description of the geometry and distribution of the material properties of the plate The orthotropic FGM plate under investigation is illustrated in Fig. 1. Length, width, and thickness of the plate are denoted by a, b, and h, respectively.

Fig. 1 Geometric parameters and coordinate system of the considered orthotropic functionally graded plate that is surrounded by an elastic foundation. 3

It is assumed that the material properties of the plate are in-plane orthotropic and transversely-graded. The rectangular Cartesian coordinate system of the plate that is denoted by xyz is coincident with the principal directions of the material properties of the pale. Hence, the relation between the stress ( σ ) and strain ( ε ) vectors of an arbitrary point of the plate may be presented as: σ = D ε,

σT = σ x

σy

ε = εx

εy

T

σ z τ xy τ xz τ yz , εz

γ xy

γ xz

(1)

γ yz

and D is the elasticity coefficients matrix. As mentioned before, the material properties are assumed to have in-plane orthotropy and transverse heterogeneity. Variations of the elastic coefficients in the transverse direction may be assumed to obey an exponential distribution. Therefore [5,9,10,44],

D ( z ) = Ce

⎛z⎞ n⎜ ⎟ ⎝h⎠

(2)

where, C is the matrix of the elastic coefficients of the bottom layer of the plate and n is the so-called material properties exponent. The plate is supported by a Winkler-type elastic foundation. Since effect of the distributed elastic foundation located above the plate is to some extent similar to that located below the plate, it is assumed that the elastic foundation surrounds the plate (exists above and below the plate). Intensity of each of the distributed stiffnesses of the foundations (located above/below the plate) is chosen as k/2 per unit area (so that the total intensity becomes equal to k in the plate theories). Furthermore, edges of the considered plate are constrained. Buckling may occur due to uniform shear tractions imposed on some or all of the edges of the plate (Fig. 1). 3 The proposed 3D cubic B-spline element Apart from their huge capabilities, the available well-known commercial finite element analysis codes such as ANSYS, NASTRAN, and ABAQUS have several shortcomings such as employing displacementbased formulations that guarantee the displacement continuity conditions only. Therefore, their strain and stress results are not continuous. Moreover, formation of sharp corners within the deformed plate at the mutual nodes of the elements is inevitable, especially when beam or shell elements are employed. One solution to this problem is using Hermitian elements instead of the conventional Lagrangian elements [45-47]. The main disadvantage of this procedure is that the resulting computational run times are not economic. Using B-spline elements may be considered as an alternative solution [48]. Wu [34] presented the only available spline interpolation scheme for buckling analysis of the isotropic homogeneous plates in the finite difference form, using the classical plate theory and neglecting the in-plane displacements of the plate. It is customary to derive the three-dimensional shape functions of the hexahedral elements based on the concept of the separation of variables. To use the standard spline interpolation in each of the x, y, and z directions of the element shown in Fig. 2, intervals of variations of the x, y, and z coordinates of the element are divided into m, n, and s equal sections, respectively. The displacement components associated with the x, y, and z directions are denoted by u, v, and w, respectively. Therefore, one-dimensional variations of the displacement component u of the element along the x axis may be expressed as: m+1

u( x) = ∑αi Bi ( x) = Bα

(3)

i =−1

where Bi (x) represents the local hill-shaped cubic B-spline (basis spline) functions (Fig. 2), α i are the unknown displacement parameters at the control points (knots) of the spline curve and m is the number of the sections. Each of the local splines is a piece-wise defined third-order polynomial which extends over

4

four consecutive sections only. In other words, each local cubic B-spline Bi ( x) has no-zero values over four consecutive sections with the section-knot xi as the center, and is defined as:

⎧ 0, ⎪ 3 ⎪ ( x − xi−2 ) , ⎪ 3 2 3 2 1 ⎪h + 3h ( x − xi −1 ) + 3h ( x − xi−1 ) − 3( x − xi −2 ) , Bi ( x) = 3 ⎨ 6h ⎪ h3 + 3h2 ( xi+1 − x) + 3h ( xi +1 − x)2 − 3( xi+1 − x)3 , ⎪ 3 ⎪( xi+2 − x) , ⎪ 0, ⎩

x < xi −2 , xi −2 ≤ x ≤ xi−1, xi−1 ≤ x ≤ xi , xi ≤ x ≤ xi+1,

(4)

xi+1 ≤ x ≤ xi +2 , xi+2 < x,

where h is the distance between the sections [49].

Fig. 2 In-plane variations of the cubic base spline functions and the nodal degrees of freedom of the full three-dimensional B-spline element used in the present research.

This localized nature of the splines results in a stiffness matrix with a narrow band when applied in a finite element model. Because the local splines extend over more than one section, the real displacement u(x) at a section corresponding to the knot xi may be related to nodal parameters of the three adjacent knots as:

1 (5) (αi−1 + 4αi + αi+1 ) 6 Because Bi −2 ( xi ), Bm+2 ( xi ) = 0 . Therefore, m+3 local splines are required to one-dimensionally subdivide the side into m sections. Since α i are not nodal variables (have not physical meaning), Ganapathi et al. u( xi ) =

[49] proposed using the displacement component and its derivative, i.e., u and u,x instead of the spline parameters α 0 ,α m at the ends of the element and α −1 ,α m +1 outside the element; where the comma symbol stands for the partial derivative. A somewhat similar method was employed by Yuen and Van Erp [50] which used linear combinations of the basis splines near the boundary to replace α −1 ,α 0 and α m ,α m −1 that 5

fall out of the plate (for the first and final elements in the considered coordinate direction, respectively) by the real displacements and first derivatives at x=0 and x=a. Therefore α will be:

α = U, x (0) U (0) α1 ..... αm−1 U (a) U, x (a)

(6)

and U is the nodal value of the u displacement component. In this formulation, the localized nature of the spline functions is retained by transforming only every other nodal parameter, rather than all, into real displacements. By choosing the element mesh such that the required rigid support coincides with an even indexed knot, the modelling of the support follows exactly the same rules as in the finite element method [50]. These concepts may be extended to three-dimensional elements through employing the idea of separation of variables, e.g.,

 u( x, y, z) = ∑∑ ∑αijk Bi ( x)Bj (y)Bk (z) = S( x, y, z)α = S( x, y, z)U(e)

(7)

 ⎧u ⎫ ⎡S ⎪ ⎪ ⎢ δ = ⎨ v ⎬ = ⎢0 ⎪w⎪ ⎢0 ⎩ ⎭ ⎣

(8)

m+1 n+1 s +1

i =−1 j =−1 k =−1

 where S and U ( e ) are the shape functions matrix and the nodal quantities vector, respectively. Eq. (7) may be used to interpret the displacement vector of an arbitrary point of the element in terms of the nodal values:

0  S 0

0⎤ ⎥ 0⎥ δ(e) = Nδ(e)  S⎥⎦

where N and δ( e ) are the total shape functions matrix and vector of the nodal displacement quantities of the element, respectively:

⎧ U(e) ⎫ ⎪ ⎪ δ(e) = ⎨ V(e) ⎬ ⎪W(e) ⎪ ⎩ ⎭

(9)

This results in a displacement field in which the degrees of freedoms at nodes with even indices represent real displacements. In this case, the nodal vectors appearing in Eq. (9) contains the zero to mixed secondorder derivatives (Fig. 2) U( e )T = U U , x U , y U , z U , xy U , xz U , yz ,

V ( e )T = V V, x V, y V, z V, xy V, xz V, yz ,

(10)

W(e)T = W W, x W, y W, z W, xy W, xz W, yz Therefore, the presented 3D cubic B-spline element guarantees C2-continuous displacement results and in contrast to the traditional elements of the well-known commercial finite element analysis codes, continuity conditions of the stresses, strains, and slopes are satisfied through using the present element. 4 The finite element form of the governing equations of the nonlinear buckling problem Generally, the strain vector contains linear and non-linear terms: ε = ε l + ε nl = L δ + ε nl = B δ ( e ) + ε nl where εT = ε x ε y ε z γ xy γ xz γ yz ,

and

6

(11) (12)

⎡ u, x ⎤ ⎢ v ⎥ ⎢ ,y ⎥ ⎢ w, z ⎥ (e) εl = ⎢ ⎥= Bδ u + v ⎢ , y ,x ⎥ ⎢ w, x + u, z ⎥ ⎢ ⎥ ⎣⎢ v, z + w, y ⎦⎥

(13)

⎧ 12 (u,2x + v,2x + w,2x ) ⎫ ⎪ ⎪ ⎪ 12 (u,2y + v,2y + w,2y ) ⎪ ⎪ ⎪ ⎪ 12 (u,2z + v,2z + w,2z ) ⎪ nl ε =⎨ ⎬ ⎪u, x u, y + v, x v, y + w, x w, y ⎪ ⎪ ⎪ ⎪ u, x u, z + v, x v, z + w, x w, z ⎪ ⎪u u + v v + w w ⎪ , y ,z ⎭ ⎩ , y ,z , y ,z

(14)

L is a matrix containing partial differentiation operators and B is the matrix of the derivatives of the shape functions: 0⎤ ⎡∂ x 0 ⎢0 ∂ 0⎥ y ⎢ ⎥ ⎢0 0 ∂z ⎥ L=⎢ B = LN (15) ⎥, ⎢∂ y ∂ x 0 ⎥ ⎢∂ z 0 ∂ x ⎥ ⎢ ⎥ ⎢⎣ 0 ∂ z ∂ y ⎥⎦ Governing equations of the elements may be obtained based on the principle of minimum total potential energy: δ Π = δU − δ W = 0 (16) where U and W are the strain energy and work of the externally applied loads, respectively: δ U = ∫ δ εT σd Ω + ∫ kwδ wdA A Ω (17) δ W = δΠ Ext . = ∫ τ xyδ ut d Γ Γ

Ω , A , and Γ are the volume, upper and bottom areas, and boundary surface of the element, respectively and ut is the tangential displacement component. In the prebuckling state, the displacements may be considered to be small. Therefore, based on Eqs. (11,17), one may write:

δU =

∫

Ω

δ ε T σ d Ω + ∫ kwδ wdA = δ δ ( e ) T A

(∫

Ω

)

B T D B d Ω + ∫ k S T R dA δ ( e ) = δ δ ( e ) T K δ ( e ) A

(18)

where D is the elasticity coefficients matrix and S and R are matrices containing 1 and 0 elements for relating transverse displacement of the nodes of the upper and bottom surfaces of the plate to the total nodal displacement vector of the plate. Since in both of the prebuckling and buckling states, Eq. (16) holds, one may use the following relation for the buckling onset: δ (δ Π ) = δ 2 Π = 0 (19) Therefore, based on Eqs. (17-19): δ δ ( e ) T K δ δ ( e ) + δ 2 Π Ext . = 0 (20)

7

Based on Eq. (11), the general strain energy of the structure contains linear as well as nonlinear terms of the strain-displacement expressions: T 1 1 U = ∫ ε T σd Ω + ∫ kw2δ wdA = ⎡⎢ ∫ (ε l + ε nl ) σd Ω + ∫ kw2δ wdA⎤⎥ (21) A A ⎦ 2 Ω 2⎣ Ω In the prebuckling state, no transverse displacements or large deformations exist. Therefore, only the linear strain components appear. On the other hand, from Eq. (16) or directly from the principle of minimum total potential energy: δ (U − W ) = 0 (22) In other words, work of the externally applied shear tractions is equal to the work of the induced stresses. Therefore, the following expression appears only after buckling:

)

(

U = Π Ext . = 12 ∫ (εnl ) σd Ω T

(23)

Ω

Since in the onset of the buckling occurrence:

τxy =τCr

(24)

at the boundaries, Eq. (23) may be rewritten in the following expanded form, based on Eq. (14):

U = Π Ext . = 12 ∫ τ Cr HxT Hy d Ω

(25)

Ω

where

H = ⎡⎣u, x T x

v, x

w, x ⎤⎦ = δ

(e) T

⎡ u,y ⎤ ⎢ ⎥ N , Hy = ⎢ v,y ⎥ = N,yδ( e ) ⎢⎣ w,y ⎥⎦ T ,x

or:

Π Ext . =

1 2

∫

Ω

τ Cr HxT H y d Ω = 12 δ ( e ) T

(∫ τ Ω

Hence, based on Eqs. (20) and (27):

δ 2 Π = δ δ ( e )T K δ δ ( e ) + δ δ ( e )T

(∫ τ Ω

Cr

Cr

)

N T, x N ,y d Ω δ ( e )

)

N *, xT N *,y d Ω δ δ ( e ) = 0

or in a compact form: δ δ ( e ) T (K + K G )δ δ ( e ) = 0 KG is the so-called geometric stiffness matrix. Noting that δ δ ( e ) ≠ 0 , the following governing equation of instability is obtained from Eq. (29): (K + K G )δ δ ( e ) = 0 Existence of nonzero solution requires that: K + KG = 0 The load amplification factor may be defined as:



τ Cr = λCrτ xy

(27) (28) (29)

(30) (31) (32)



where, τ xy is the imposed shear traction. Therefore, based on Eqs. (18,28,32):  T *T * ∫ B D B d Ω + λCrτ xy ∫ N, x N,y d Ω = 0 Ω

(26)

Ω

or in a compact form: G ( λ Cr ) = 0

(33) (34)

8

5 Condensing the final governing system of equations Before solving the final system of equations, the boundary conditions have to be incorporated in the whole finite element system of equations of the plate. The final form of the resulting system may be written as:  (35) Kδδ = 0 It is not recommended to solve the final equations as they appear. Otherwise, huge numerical computations should be accomplished. The degrees of freedom may be categorized as those that can be  determined directly δ1 and the quantities that contain the second derivatives of the displacement  components and cannot be determined directly δ2 :  ˆ U ˆ ˆ ˆ ˆ V ˆ ˆ ˆ ˆ W ˆ ˆ ˆ δ1 = U U U V V V W W W ,x ,y ,z ,x ,y ,z ,x ,y ,z (36)  ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ δ2 = U U U V V V W W W , xy , xz , yz , xy , xz , yz , xy , xz , yz

( )

( )

ˆ ” denotes a vector contains all U displacement parameters of all nodes of the plate. Expressions like “ U Some of the equations may be merged into the each other in order that the final system contains only the  zero and first derivative degrees of freedom ( δ1 ). In this regard, Eq. (35) may be rewritten as:  ⎡ K 11 K 12 ⎤ ⎧⎪ δ δ1 ⎫⎪ ⎧ 0 ⎫ (37) ⎢ ⎥⎨  ⎬ = ⎨ ⎬ ⎣ K 21 K 22 ⎦ ⎩⎪δ δ 2 ⎭⎪ ⎩ 0 ⎭ Therefore, from the lower part of the resulting system of equations, one may infer:   (38) δ δ 2 = − K −221 K 21δ δ 1 Substituting Eq. (38) into the upper part of the system of equations appeared in Eq. (37), leads to:  K11 − K12 K −221K 21 δ δ1 = 0 (39)

(

)

or in a compact form:  Kδ δ = 0 ⇒ K = 0

(40)

1

The resulting condensed system (40) includes influence of all nodal quantities but its dimensions are considerably smaller.

6 Results and discussions 6.1 Definition of the base information of the treated problems In the present section, results of shear buckling of orthotropic functionally graded plates surrounded by elastic foundations are reported. When using plate theories without transverse flexibility [51-53], i.e., with constant w component, these results are equivalent to those of plates resting on elastic foundations. In addition to the elastic foundation, edges of the plate may also be supported. In this regard, three types of edge conditions are considered. The employed displacement-type edge conditions and their relevant mathematical interpretations are: i) All edges are simply supported but in-plane movable: w (0, y , z ) = w ( a , y , z ) = w ( x,0, z ) = w ( x, b, z ) = 0 (41) ii) Edges parallel to the y-axis (x=0, a) are clamped and the remaining edges are free: u, v, w (0, y , z ) = u, v, w ( a , y , z ) = 0 (42) iii) Edges parallel to the x-axis (y=0, b) are clamped and the remaining edges are free: u, v , w ( x,0, z ) = u, v , w ( x , b, z ) = 0 (43) Therefore, in the two last edge conditions, the shear tractions are imposed only on two opposite edges of the plate. Although this case is the most practical one, it has not been considered by any researcher 9

before. The rectangular plates have material properties with in-plane orthotropy and transverse gradation. The adopted elasticity coefficients matrix of the materials is similar to the one used by Phan [9]:

0 0 0 ⎤ ⎡7.3802 2.3121 1.8682 ⎢ 2.3121 173.406 2.3121 0 0 0 ⎥⎥ ⎢ ⎢1.8682 2.3121 7.380 0 0 0 ⎥ n⎛⎜⎝ hz ⎞⎟⎠ D (z) = ⎢ ⎥ e GPa 0 0 0 3.445 0 0 ⎢ ⎥ ⎢ 0 0 0 0 3.445 0 ⎥ ⎢ ⎥ 0 0 0 0 1.378⎦ ⎣ 0

(44)

The side to thickness ratio is adopted as: b/h=10. Two aspect ratios a/b=1,2 and several exponents of the transverse material properties variations and elastic foundation coefficients are employed in the parametric study. Based on Reddy [54], a/h or b/h<20 denotes a thick orthotropic plate. Therefore, in contrast to the present three-dimensional elasticity approach, results of the plate theories may encounter accuracy problem for the considered conditions. In deriving present results a mesh composed of 20, 20, and 8 elements in the x, y, and z directions respectively (with overall 3200 elements) is employed for the FGM plates. On the other hand, for the orthotropic plates (without transverse material properties gradations) used in the verification example, a mesh containing 20*20*4 elements is utilized. 6.2 Presenting and verification of the results for an orthotropic plate Convergence of the finite element results of the present research is investigated through comparing results of successive refinements of the elements size. The transverse gradation of the material properties considered in the present paper has an exponential form. Therefore, for the considered boundary conditions and the material properties distribution, only the results available on shear buckling of isotropic homogenous plates can be used for verification purposes. In this case, the comparative studies may not lead to a strong verification. For this reason, it is preferred to compare present results with those obtained using a well-known commercial finite element analysis code (ANSYS). ANSYS code uses Lagrangian (C0-continuous) elements, whereas the elements employed in the present research are 3D cubic B-spline (C2-continuous) ones. Therefore, in addition to the quite different natures of formulations of the elements, using present elements is expected to lead to more accurate results. In this regard, homogeneous orthotropic (n=0) plates without elastic foundation are considered. Using the elasticity coefficients matrix appeared in Eq. (44) is equivalent to using the following material properties for the orthotropic plates: E1 = E3 = 6.89 GPa, E2 = 172.2 GPa, G12 = G13 = 3.44 GPa, G23 = 1.37GPa,ν 12 = 0.01,ν 13 = ν 23 = 0.25 Results of the critical shear buckling loads are obtained for two aspect ratios (a/b=1,2) and two types of edge conditions (the first two types of the introduced edge conditions) and compared with those extracted from ANSYS software in Table 1. Table1 A comparison between present results for critical shear buckling loads (in GPa) of the homogenous orthotropic rectangular plates and results of ANSYS code for different boundary conditions and aspect ratios (b/h=10, h=0.1m). Edge condition Aspect ratio a/b=1 a/b=2 Present ANSYS Present ANSYS Simply supported 0.89157 0.91135 0.77305 0.78606 Clamped edges at x=0,a 0.53281 0.54821 0.11241 0.11557

10

(b) 1

1

0.8

0.8

0.6

0.6

y (m )

y (m )

(a)

0.4

0.4

0.2

0.2

0 0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

0.6

0.8

1

x (m)

x (m)

(c) (d) Fig. 3 A comparison between the first two buckling mode shapes predicted by ANSYS [(a) and (b)] and present approach [(c) and (d)], for the homogenous orthotropic square plate (a/b=1).

Results presented in Table 1 reveal that the relative differences between present results and results of ANSYS code are less than 3%. On the other hand, the relevant first two buckling mode shapes predicted based on the present approach and ANSYS are compared in Figs. 3 and 4, for a/b=1 and 2, respectively. Comparing these results reveals that the mode shapes predicted by the present approach are accurate as well. Therefore, both sufficiency of the chosen number of elements and accuracy of the developed formulations are verified. As may be expected, the buckling waves are oblique (they are not parallel to the edges of the plate) and the buckling shape becomes more complicated in higher modes. In isotropic homogeneous plates, the shear buckling waves constitute a 45 degrees angle with respect to sides of the rectangular plate; because origin of the shear buckling is not the shear itself, but the compressive normal stress that is established by the shear stress. As Mohr circle confirms, the maximum compressive stress direction (principal direction) makes a 45 degrees angle with the maximum shear direction (along the dges). But the situation is more complicated in the orthotropic plates; since no principal direction may be 11

defined for the orthotropic plate, due to the existence of the extension-shear coupling. As it will be discussed in the next sections, the first and second buckling modes do not necessarily correspond to formation of one and two buckling waves, respectively. In some circumstances, e.g., non-square plates or plates with orthotropic material properties or elastic foundations, the first buckling mode may be associated with multiple or even local buckling waves.

(b) 1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

(a)

0.4

0.4

0.2

0.2

0 0

0.5

1 x (m)

1.5

2

0 0

0.5

1 x (m)

1.5

(c) (d) Fig. 4 A comparison between the first two buckling mode shapes predicted by ANSYS [(a) and (b)] and present approach [(c) and (d)], for the homogenous orthotropic rectangular plate (a/b=2).

6.3 A parametric study for simply supported orthotropic FGM plates with elastic foundations

In the present section, influences of the material properties exponent, aspect ratio, and foundation stiffness on shear buckling loads of a simply supported orthotropic FGM plate are evaluated. Results of the base case of an orthotropic plate (n=0) without elastic foundation have been reported in the foregoing section. Influences of gradation of the material properties in the transverse direction (n=0.1) on the buckling mode shapes may be evaluated by comparing the first two buckling mode shapes illustrated in Figs. 5 and 6 with those of Figs. 3 and 4, respectively. Figs. 5 and 6 reveal that the transverse gradation of the material properties affects the buckling mode shapes and consequently, the buckling mechanism.

12

2

1

0.8

0.8

0.6

0.6

y (m)

y (m)

1

0.4 0.2 0 0

0.4 0.2

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x (m)

0.6

0.8

1

0.6

0.8

1

x (m)

(b) 1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

(a)

0.4 0.2 0 0

0.4 0.2

0.2

0.4

0.6

0.8

0 0

1

x (m)

(c)

0.2

0.4 x (m)

(d)

Figure 5 First four buckling mode shapes of a simply supported square orthotropic FGM plate (b/h=10, n=0.1) without elastic foundation.

Using an elastic foundation postpones formation of the buckling waves and consequently, the buckling occurrence. So that when the foundation is stiff enough, buckling is accompanied with significant kinking, crimpling, and wrinkling. For a plate without an elastic foundation, occurrence of the global long buckling waves is more likely whereas in presence of a distributed elastic foundation, occurrence of relatively short and narrow local buckling waves is more likely. This fact may be easily deduced by comparing Figs. 7 and 8 that are associated with plates with an elastic foundation with k=1e10 (N/m3) respectively with Figs. 5 and 6. In other words, some of the first modes that appear in plates without elastic foundations may disappear and only modes corresponding to higher strain energies and buckling loads may appear.

13

1

0.8

0.8

0.6

0.6

y (m)

y (m)

1

0.4

0.4

0.2 0 0

0.2

0.5

1 x (m)

1.5

(a)

0.5

1 x (m)

1.5

2

(b) 1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

0 0

2

0.4 0.2 0 0

0.4 0.2

0.5

1 x (m)

1.5

0 0

2

0.5

1 x (m)

1.5

2

(c) (d) Figure 6 First four buckling mode shapes of a simply supported rectangular orthotropic FGM plate

1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

(a/b=2, b/h=10, n=0.1) without elastic foundation.

0.4 0.2 0 0

0.2

0.2

0.4

0.6

0.8

0 0

1

x (m)

(a)

0.4

0.2

0.4

0.6 x (m)

(b)

14

0.8

1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

1

0.4 0.2 0 0

0.4 0.2

0.2

0.4

0.6

0.8

0 0

1

0.2

x (m)

0.4

0.6

0.8

1

x (m)

(c)

(d)

Figure 7 Influence of the elastic foundation on the crimpling and buckling wave patterns of the first four

1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

buckling mode shapes of a square simply supported plate (b/h=10, n=0.1, k=1e10 (N/m3)).

0.4 0.2 0 0

0.2

0.5

1 x (m)

1.5

0 0

2

(a)

0.5

1 x (m)

1.5

2

(b) 1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

0.4

0.4 0.2 0 0

0.4 0.2

0.5

1 x (m)

1.5

0 0

2

0.5

1 x (m)

1.5

2

(c) (d) Figure 8 Influence of the elastic foundation on the crimpling and buckling wave patterns of the first four buckling mode shapes of a rectangular simply supported plate [a/b=2, b/h=10, n=0.1, k=1e10 (N/m3)]. Influences of the aspect ratio, gradation exponent of the material properties, and stiffness of the elastic foundation on the numerical values of the shear buckling loads (tractions) may be inferred from the results given in Table 1. In this regard, two aspect ratios (a/b=1, 2), three material properties gradation exponents (n=0.1, 0.3, 0.5), and four elastic foundation stiffness values (k=0,1e9, 5e9, 1e10 N/m3) are considered and the results are reported for the first four buckling modes of the plate. For plates with 15

elastic foundation, stiffness of the elastic foundation is so chosen that it is comparable with stiffness of the plate. For plates with or without elastic foundations, buckling waves may be formed in lower loads when area of the plate increases, i.e., when higher aspect ratios are employed. When stiffness of the plate in the transverse direction becomes in an order of magnitude equal or less than the stiffness of the elastic foundation, the plate with greater aspect ratio experiences more difficulties in moving away the foundation in the transverse direction. Hence, the buckling loads become higher when greater aspect ratios are employed. Increasing exponent of the gradation of the material properties in the transverse direction leads to higher rigidities for the plate and consequently, higher buckling loads. On the other hand, as expected, the buckling loads increase as the coefficient of the elastic foundation increases. In contrast to the shear buckling loads of the isotropic homogeneous plates, and due to existence of the shear-extension-bending coupling in the orthotropic FGM plate, the differences between values of the buckling loads correspond to the successive buckling modes are not significant (Table 1). The shearextension coupling is due to the in-plane orthotropy of the material properties and the extension-bending coupling is established due to the asymmetric material properties distribution in the transverse direction. Therefore, the shear tractions invoke both the extension and bending modes simultaneously. A fact that makes the situation significantly complicated. Moreover, differences between results of the numerical buckling loads of the plates with a/b=1 and a/b=2 aspect ratios are smaller for higher stiffnesses of the elastic foundation. Therefore, effect of the elastic foundation is the dominant parameter in higher stiffnesses of the elastic foundation. Table 1 Numerical values of the first four shear buckling loads for simply supported orthotropic FGM plates with various aspect ratios, exponents of the material properties gradation, and coefficients of the elastic foundation (b/h=10). Aspect Buckling load ratio (GPa)

k=0 (N/m3)

k=1e9 (N/m3)

k=5e9 (N/m3)

k=1e10 (N/m3)

n=0.1

n=0.3 n=0.5

n=0.1

n=0.3 n=0.5

n=0.1

n=0.3

n=0.5

n=0.1

n=0.3 n=0.5

0.926

1.025

1.136

1.079

1.1813 1.294

1.368

1.486

1.617

1.555

1.685

1.829

0.928

1.027

1.137

1.089

1.1931 1.308

1.368

1.487

1.618

1.556

1.685

1.829

1.241

1.374

1.523

1.329

1.4641 1.615

1.522

1.664

1.823

1.67

1.819

1.985

1.253

1.387

1.538

1.335

1.4711 1.623

1.524

1.666

1.826

1.671

1.820

1.987

0.831

0.919

1.019

1.016

1.109

1.213

1.355

1.466

1.59

1.583

1.707

1.846

0.8431

0.933

1.034

1.021

1.114

1.219

1.355

1.467

1.59

1.583

1.708

1.846

0.971

1.074

1.191

1.127

1.233

1.352

1.441

1.564

1.701

1.656

1.790

1.940

0.973

1.077

1.193

1.138

1.247

1.369

1.442

1.565

1.701

1.656

1.791

1.941

a/b=1

a/b=2

16

6.4 Considering more complicated boundary conditions

In majority of the engineering applications, the shear stresses are not uniformly distributed on all edge of the plate. For example, if two opposite edges of the plate are clamped and the remaining edges are subjected to in-plane shear, the shear tractions distributed on the clamped edges are surely not uniform. Influence of the stiffness of the elastic foundation on the buckling mode shapes of the square orthotropic FGM (n=0.1) plate may be traced through comparing Figs. 9-11 that correspond to successive increasing the stiffness of the elastic foundation form k=0 to 1e9 (N/m3) to 1e10 (N/m3), respectively for plates whose clamped edges are parallel to y axis (x=0,a). Comparing these figures with each other confirms the idea that as the stiffness of the foundation increases, the buckling waves becomes local, shorter and narrower. This idea may also be confirmed by comparing Figs. 12 and 13 that are plotted for rectangular plates (a/b=2) without and with high stiffness [k=1e10 (N/m3)] elastic foundation, respectively. Numerical values of the shear buckling loads listed in Table 2 confirm the aforementioned conclusions. Moreover, in contrast to the previous cases, differences between results of the a/b=1 and a/b=2 aspect ratios are significant. However, since in higher stiffnesses of the elastic foundation effect of the elastic

1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

foundation is the dominant, these differences diminish as the stiffness of the foundation increases.

0.4 0.2 0 0

0.4 0.2

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x (m)

(a)

0.8

1

0.6

0.8

1

(b) 1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

0.6 x (m)

0.4 0.2 0 0

0.4 0.2

0.2

0.4

0.6

0.8

0 0

1

x (m)

0.2

0.4 x (m)

(c) (d) Figure 9 First four buckling mode shapes of a square orthotropic FGM plate with clamped edges at x=0,a (b/h=10 , n=0.1, k=0). 17

1

0.8

0.8

0.6

0.6

y (m)

y (m)

1

0.4 0.2 0 0

0.4 0.2

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x (m)

(a)

0.8

1

0.6

0.8

1

(b) 1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

0.6 x (m)

0.4 0.2 0 0

0.4 0.2

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x (m)

x (m)

(c) (d) Figure 10 First four buckling mode shapes of a square orthotropic FGM plate with clamped edges at x=0,a [b/h=10 , n=0.1, k=1e9 (N/m3)]

1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

The same trends may be observed for a square orthotropic FGM plate whose opposite edges at y=0,b are clamped while the remaining opposite edges are subjected to uniform in-plane shear tractions. These results are reported in Table 3. Since material properties are different in the principal directions of the material properties oriented along x and y axes (i.e., influence of the orthotropic nature of the material properties), these results are not identical to those reported for plates whose x=0,a edges are clamped.

0.4 0.2 0 0

0.2

0.2

0.4

0.6

0.8

0 0

1

x (m)

(a)

0.4

0.2

0.4

0.6 x (m)

(b)

18

0.8

1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

1

0.4 0.2

0.4 0.2

0 0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

x (m)

0.6

0.8

1

x (m)

(c) (d) Figure 11 First four buckling mode shapes of a square orthotropic FGM plate with clamped edges at x=0,a [b/h=10, n=0.1, k=1e10 (N/m3)].

Table 2 Numerical values of the first four shear buckling loads of orthotropic FGM plates with clamped x=0,a edges and various aspect ratios, exponents of the material properties gradation, and coefficients of the elastic foundation (b/h=10). Aspect Buckling load ratio (GPa)

a/b=1

a/b=2

k=0 (N/m3)

k=1e9 (N/m3)

n=0.1

n=0.3 n=0.5

λ1

0.554

0.612

0.677

λ2

0.656

0.725

λ3

0.827

λ4

n=0.1

k=5e9 (N/m3)

k=1e10 (N/m3)

n=0.3 n=0.5

n=0.1

n=0.3

n=0.5

n=0.1

n=0.3 n=0.5

0.85

0.925

1.006

1.072

1.168

1.274

1.181

1.287

1.403

0.802

0.871

0.952

1.042

1.074

1.17

1.277

1.181

1.287

1.404

0.915

1.013

0.975

1.068

1.171

1.132

1.238

1.357

1.219

1.333

1.459

0.883

0.977

1.082

0.976

1.072

1.179

1.134

1.241

1.361

1.219

1.333

1.46

λ1

0.133

0.147 0.1635

0.633

0.683

0.738

0.921

0.995

1.075

1.073

1.159

1.253

λ2

0.242

0.267

0.295

0.648

0.698

0.751

0.921

0.995

1.077

1.073

1.159

1.253

λ3

0.415

0.459

0.508

0.683

0.738

0.801

0.958

1.038

1.125

1.102

1.192

1.292

λ4

0.439

0.486

0.537

0.710

0.773

0.84

0.959

1.038

1.126

1.102

1.192

1.292

19

1

0.8

0.8

0.6

0.6

y (m)

y (m)

1

0.4 0.2 0 0

0.2

0.5

1 x (m)

1.5

0 0

2

(a)

0.5

1 x (m)

1.5

2

0.5

1 x (m)

1.5

2

(b) 1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

0.4

0.4 0.2 0 0

0.4 0.2

0.5

1 x (m)

1.5

0 0

2

1

1

0.8

0.8

0.6

0.6

y (m)

y (m)

(c) (d) Figure 12 First four buckling mode shapes of a rectangular orthotropic FGM plate with clamped edges at x=0,a (a/b=2, b/h=10, n=0.1, k=0).

0.4 0.2 0 0

(a)

0.4 0.2

0.5

1 x (m)

1.5

0 0

2

(b)

20

0.5

1 x (m)

1.5

2

1

0.8

0.8

0.6

0.6

y (m)

y (m)

1

0.4 0.2

0.4 0.2

0 0

0.5

1 x (m)

1.5

0 0

2

0.5

1 x (m)

1.5

2

(c) (d) Figure 13 First four buckling mode shapes of a rectangular orthotropic FGM plate with clamped edges at x=0,a [a/b=2, b/h=10, n=0.1, k=1e10 (N/m3)]. Table 3 Numerical values of the first four shear buckling loads of square orthotropic FGM plates with clamped y=0,b edges and various aspect ratios, exponents of the material properties gradation, and coefficients of the elastic foundation (b/h=10). Buckling load (GPa)

k=0 (N/m3) n=0.1

k=1e9 (N/m3)

n=0.3 n=0.5

k=5e9 (N/m3)

k=1e10 (N/m3)

n=0.1

n=0.3 n=0.5

n=0.1

n=0.3

n=0.5

n=0.1

n=0.3 n=0.5

λ1

0.244

0.269

0.297

0.245

0.27 0.299

0.248

0.274

0.274

0.252

0.279

0.308

λ2

0.244

0.269

0.297

0.245

0.27 0.299

0.248

0.274

0.274

0.252

0.279

0.308

λ3

0.283

0.313

0.346

0.285

0.314 0.348

0.291

0.321

0.321

0.297

0.328

0.362

λ4

0.283

0.313

0.346

0.285

0.314 0.348

0.291

0.321

0.321

0.297

0.328

0.362

7 Conclusions Non-linear shear buckling analysis of orthotropic heterogeneous FGM plates surrounded by Winklertype elastic foundations is investigated in the present paper, for the first time. In this regard, a new 3D cubic B-spline element is introduced and employed to develop the finite element form of the instability equations. The material properties are assumed to have in-plane orthotropy and transverse heterogeneity. Some of novelties of the present research are: - Presenting a nonlinear shear buckling analysis based on a general geometric stiffness matrix concept. - Performing the analysis for the orthotropic FGM plates. Therefore, it is the first time that influence of the shear-extension-bending coupling on the buckling load is investigated. - Incorporating simultaneous effects of the elastic foundation and various edge conditions. - The most accurate three-dimensional elasticity theory is employed instead of the approximate plate theories. - In contrast to all of the available displacement-based buckling analyses that have employed either the C0-continuous commercial finite element codes or semi-analytical methods, present formulations are C2-continuous. 21

- Presenting a 3D cubic B-spline elements. - Effects of both prebuckling and buckling states are considered. Comparing present results with those obtained from the ANSYS finite element analysis code, verifies accuracy of the present predictions. Moreover, the following practical conclusions are extracted based on the obtained results: - In contrast to the isotropic homogeneous plates, the oblique buckling waves do not constitute a 45 degrees angle. - Furthermore, for an orthotropic FGM plate with uniform in-plane shear tractions, differences of the successive buckling loads are remarkably lower than those of plates with two opposite clamped edges. The later edge condition resembles the more encountered edge conditions in the engineering applications. - Due to using a 3D elasticity approach instead of the plate theories and considering the transverse flexibility of the plate, influences of the elastic foundations located at top and bottom surfaces of the plate are different. - For the considered conditions, buckling loads are smaller for higher aspect ratios. - Buckling loads are smaller for edge conditions that permit more movability. - In contrast to the isotropic plates, the buckling waves do not constitute a 45 degrees angle with the edges. Furthermore, the critical buckling load does not necessary correspond to one buckling wave. - Due to the existence of the shear-extension-bending coupling in the orthotropic FGM plate, the transverse gradation of the material properties affects both the buckling mode shapes and the buckling loads of the mentioned plate. - The surrounding elastic foundation postpones the buckling occurrence and leads to local, narrower, and shorter buckling waves. - Presence of the elastic foundation leads to significant kinking, crimpling, and wrinkling at the buckling instant.

References

[1] Kar A, Kanoria M. Generalized thermoelastic functionally graded orthotropic hollow sphere under thermal shock with three-phase-lag effect. Europ J Mech A/Solids 2009; 28: 757–767. [2] Sofiyev AH, Karaca Z. The vibration and buckling of laminated non-homogeneous orthotropic conical shells subjected to external pressure. Europ J Mech A/Solids 2009; 28: 317–328. [3] Sofiyev AH, Kuruoglu N. Buckling analysis of nonhomogeneous orthotropic thin-walled truncated conical shells in large deformation. Thin-Walled Struct 2013; 62: 131–141. [4] Aghazadeh Mohandesi J, Shahosseini MH, Parastar Namin R. Tensile behavior of functionally graded steels produced by electroslag remelting. Metallurgical Mater Trans A 2006; 37A: 21252132. [5] Ootao Y, Tanigawa Y. Three-dimensional solution for transient thermal stresses of an orthotropic functionally graded rectangular plate. Compos Struct 2007, 80: 10–20. [6] Peng X-L, Li X-F. Elastic analysis of rotating functionally graded polar orthotropic disks. Int J Mech Sci 201; 60: 84–91. [7] Wang X, Sudak LJ. Three-dimensional analysis of multi-layered functionally graded anisotropic cylindrical panel under thermomechanical loading. Mech Mater 2008; 40: 235–254. [8] Pelletier JL, Vel SS. An exact solution for the steady-state thermoelastic response of functionally graded orthotropic cylindrical shells. Int J Solids Struct 2006; 43: 1131–1158. [9] Pan E. Exact solution for functionally graded anisotropic elastic composite laminates. J Compos Mater 2003; 37:1903–20. [10] Wen PH, Sladek J, Sladek V. Three-dimensional analysis of functionally graded plates. Int J Numer Meth Eng 2011; 87: 923–942. 22

[11] Jha DK, Kant T, Singh RK. A critical review of recent research on functionally graded plates. Compos Struct 2013; 96: 833–849. [12] Duc ND, Tung HV. Mechanical and thermal postbuckling of higher order shear deformable functionally graded plates on elastic foundations. Compos Struct 2011; 93: 2874–2881. [13] Abrate S. Free vibration, buckling and static deflections of functionally graded plates. Compos Sci Tech 2006; 66: 2383-2394. [14] Zhao X, Lee YY, Liew KM. Mechanical and thermal buckling analysis of functionally graded plates. Compos Struct 2009; 90: 161-171. [15] Oyekoya OO, Mba DU, AM El-Zafrany. Buckling and vibration analysis of functionally graded composite structures using the finite element method. Compos Struct 89 (2009) 134-142. [16] Lee YY, Zhao X, Reddy JN. Postbuckling analysis of functionally graded plates subject to compressive and thermal loads. Comput Meth Appl Mech Eng 2010; 199: 1645–1653. [17] Thai HT, Choi DH. An efficient and simple refined theory for buckling analysis of functionally graded plates. Appl Math Model 2012; 36: 1008-1022. [18] Shariyat M. Dynamic buckling of imperfect laminated plates with piezoelectric sensors and actuators subjected to thermo-electro-mechanical loadings, considering the temperature-dependency of the material properties. Compos Struct 2009; 88: 228-239. [19] Shariyat M. A double-superposition global-local theory for vibration and dynamic buckling analyses of viscoelastic composite/sandwich plates: a complex modulus approach. Arch Appl Mech 2011; 81: 1253-1268. [20] Shariyat M. A nonlinear double-superposition global-local theory for dynamic buckling of imperfect viscoelastic composite/sandwich plates: A hierarchical constitutive model. Compos Struct 2011; 93: 1890-1899. [21] Shariyat M. Vibration and dynamic buckling control of imperfect hybrid FGM plates with temperature-dependent material properties subjected to thermo-electro-mechanical loading conditions. Compos Struct 2009; 88: 240-252. [22] Alipour MM, Shariyat M. Semi-analytical buckling analysis of heterogeneous variable thickness viscoelastic circular plates on elastic foundations. Mech Res Commun 2011; 38: 594-601. [23] Alipour MM, Shariyat M. Semianalytical solution for buckling analysis of variable thickness twodirectional functionally graded circular plates with nonuniform elastic foundations. J Eng Mech 2013; 139: 664-676. [24] Budiansky B, Connor RW. Buckling stresses of clamped rectangular flat plates in shear. NACA Tech. Note No. 1559, 1948. [25] Smith ST, Bradford MA, Oehlers DJ. Elastic buckling of unilaterally constrained rectangular plates in pure shear. Eng Struct 1999; 21: 443–453. [26] Azikov NS, Vasiliev VV, Paterekas AD. Buckling of composite plates under compression and shear. Mech Compos Mater 1990; 2: 351-3. [27] Loughlan J. The shear buckling behaviour of thin composite plates with particular reference to the effects of bend-twist coupling. Int J Mech Sci 2001; 43: 771-792. [28] Kim KD,. Lomboy GR, Han SC. A co-rotational 8-node assumed strain shell element for postbuckling analysis of laminated composite plates and shells. Comput Mech 2003; 30: 330–342. [29] Han S-C, Lee S-Y, Rus G. Postbuckling analysis of laminated composite plates subjected to the combination of in-plane shear, compression and lateral loading. Int J Solids Struct 2006; 43: 57135735. [30] Lopatin AV, Korbut YB. Buckling of clamped orthotropic plate in shear. Compos Struct 2006; 76: 94-98. [31] Papadopoulos L, Kassapoglou C. Shear buckling of rectangular composite plates composed of concentric layups. Composites: Part A 2007; 38: 1425–1430. [32] Shufrin I, Eisenberger M. Shear buckling of thin plates with constant in-plane stresses. Int J Struct Stabil Dyn 2007; 7: 179–192. 23

[33] Shufrin I, Rabinovitch O, Eisenberger M. Buckling of laminated plates with general boundary conditions under combined compression, tension, and shear—A semi-analytical solution. ThinWalled Struct 2008; 46: 925–938. [34] Wu L-Y, Wu C-H, Huang H-H. Shear buckling of thin plates using the spline collocation method. Int J Struct Stabil Dyn 2008; 8: 645–664. [35] Tarján G, Sapkás Á, Kollár LP. Stability analysis of long composite plates with restrained edges subjected to shear and linearly varying loads. J Reinf Plast Compos 2010; 29: 1386-1398. [36] Wu L, Nie J. Elastic buckling of steel-concrete composite slabs subjected to pure shear load. Dongnan Daxue Xuebao/J Southeast Univ 2011; 41: 622-629. [37] Brighenti R, Carpinteri A. Buckling and fracture behaviour of cracked thin plates under shear loading. Mater Des 2011; 32: 1347-1355. [38] Chirica I, Beznea E-F. Buckling behavior of the multiple delaminated composite plates under shear and axial compression. Comput Mater Sci 2012; 64: 173-178. [39] Yang D, Huang Y, Li G. Shear buckling analysis of anisotropic rectangular plates. Yingyong Lixue Xuebao/Chinese J Appl Mech 2012; 29: 220-224. [40] Shariyat M, Eslami MR. Dynamic buckling and postbuckling of imperfect orthotropic cylindrical shells under mechanical and thermal loads, based on the three dimensional theory of elasticity. Trans ASME, J Appl Mech 1999; 66: 476-484. [41] Gu H, Chattopadhyay A. Three-dimensional elasticity solution for buckling of composite laminates. Compos Struct 2000; 50: 29-35. [42] Uymaz B, Aydogdu M. Three dimensional shear buckling of FG plates with various boundary conditions. Compos Struct 2013; 96: 670–682. [43] Asemi K, Shariyat M, Salehi M, Ashrafi H. A full compatible three-dimensional elasticity element for buckling analysis of FGM rectangular plates subjected to various combinations of biaxial normal and shear loads. Accepted in Finite Elements in Analysis and Design. [44] Fernando Ramirez, Paul R. Heyligera, Ernian Panb. Static analysis of functionally graded elastic anisotropic plates using a discrete layer approach. Composites: Part B 37 (2006) 10–20. [45] Shariyat M, Eslami MR. Isoparametric finite-element thermoelastoplastic creep analysis of shells of revolution. Int J Press Vess Pip1996; 68(3): 249-259. [46] Shariyat M. A nonlinear Hermitian transfinite element method for transient behavior analysis of hollow functionally graded cylinders with temperature-dependent materials under thermomechanical loads. Int J Press Vess Piping 2009; 86: 280-289. [47] Zienkiewicz OC, Taylor RL. The finite element method: its basis and fundamentals. 6th edition, Butterworth-Heinemann, 2005. [48] Hollig K. Finite Element Methods with B-Splines. Society for Industrial Applied Mathematics Philadelphia, 2003. [49] Ganapathi M, Patel BP, Saravanan J, Touratier M. Shear flexible curved spline beam element for static analysis. Finite Elem Anal Des 1999; 32: 181-202. [50] Yuen SW, Van Erp GM. Transient analysis of thin-walled structures using macro spline finite elements. Eng Struct 1999; 21: 255–266. [51] Shariyat M, Alipour MM. Semi-analytical consistent zigzag-elasticity formulations with implicit layerwise shear correction factors for dynamic stress analysis of sandwich circular plates with FGM layers. Composites: Part B 2013; 49: 43-64. [52] Kiani Y, Eslami MR. An exact solution for thermal buckling of annular FGM plates on an elastic medium. Compos Part B 2013; 45: 101–110. [53] Neves AMA, Ferreira AJM, Carrera E, Cinefra M, Roque CMC, Jorge RMN, Soares CMM. Static, free vibration and buckling analysis of isotropic and sandwich functionally graded plates using a quasi-3D higher-order shear deformation theory and a meshless technique. Compos Part B 2013; 44: 657–674.

24

[54] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. CRC Press, 2nd Edition, 2004.

25