Finite element formulation of a generalized higher order shear deformation theory for advanced composite plates

Composite Structures 96 (2013) 545–553

Contents lists available at SciVerse ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Finite element formulation of a generalized higher order shear deformation theory for advanced composite plates J.L. Mantari, C. Guedes Soares ⇑ Centre for Marine Technology and Engineering (CENTEC), Instituto Superior Técnico, Technical University of Lisbon, Av. Rovisco Pais, 1049-001 Lisbon, Portugal

a r t i c l e

i n f o

Article history: Available online 14 August 2012 Keywords: Functionally graded plates Finite element Shear deformation Static analysis

a b s t r a c t This paper presents a generalized higher order shear deformation theory (HSDT) and its finite element formulation for the bending analysis of advanced composite plates such as functionally graded plates (FGPs). New shear strain shape functions are presented. The generalized HSDT accounts for non-linear and constant variation of in-plane and transverse displacement respectively through the plate thickness, complies with plate surface boundary conditions and do not require shear correction factors. The generalized finite element code is base on a continuous isoparametric Lagrangian finite element with seven degrees of freedom per node. Numerical results for different side-to-thickness ratio, aspect ratios, volume fraction, and simply supported boundary conditions are compared. Results show that new nonpolynomial HSDTs solved by the proposed generalized finite element technique are more accurate than, for example, the well-known trigonometric HSDT, having the same DOFs. It is concluded that some nonpolynomial shear strain shape functions are more effective than the polynomials counterparts. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Functionally graded materials (FGMs) are advanced composite materials that have been proposed, developed and successfully used in industrial applications since 1984 [1]. FGMs, initially used in the aerospace as thermal barrier material, nowadays have extensive application in nuclear reactors, chemical plants, and other applications such as biomechanical, optical, automotive, electronic, mechanical, civil and shipbuilding industries. FGMs are both macroscopically and microscopically heterogeneous advanced composites which are for example made from a mixture of ceramics and metals with continuous composition gradation from pure ceramic on one surface to full metal on the other. Such gradation leads to a smooth change in the material profile as well as the effective physical properties. In addition, FGMs have interesting thermomechanical properties that can alleviate or eliminate the deformation of structural components [2]. Its counterpart, classical composite materials, suffers from discontinuity of material properties at the interface of the layers and constituents of the composite. Therefore the stress fields in these regions create interface problems and thermal stress concentrations under high temperature environments. Furthermore, large plastic deformation of the interface may trigger the initiation and propagation of cracks in the material [3]. These problems can be decreased by gradually chang-

⇑ Corresponding author. Tel.: +351 218 417607; fax: +351 218 474015. E-mail address: [email protected] (C. Guedes Soares). 0263-8223/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruct.2012.08.004

ing the volume fraction of constituent materials and tailoring the material for the desired application, as in FGMs. Literature survey shows that many papers dealing with static and dynamic behaviour of functionally graded materials (FGMs) have been published recently. An interesting literature review of this work may be found in the paper of Birman and Byrd [4]. More recent work on 3D exact and closed-form analytical solutions and an updated review can be found in Carrera et al. [5] and Mantari and Guedes Soares [6]. Thus the following analysis addresses mainly the finite element formulation. It can be said that not much work has been done for the finite element (FE) analyses of advanced composites structures such as functionally graded plates (FGPs) as for classical composites. Often, the studies on FGPs are largely devoted to thermal stress analysis [7–9] and fracture analysis FG plates and shells [10,11]. In addition, some FE models have been already proposed for the study of FG plates and shells [7,8,12,13]. Recently, Chinosi and Della Croce [14] used a mixed interpolated finite element to study cylindrical shells made of functionally graded materials. Della Croce and Venini [15] previously developed a finite element for functionally graded Reissner–Mindlin plates. Carrera et al. [16] employed the concept of virtual displacements to obtain finite element solutions of functionally graded plates subjected to transverse mechanical loadings. Later Reissner’s Mixed Variational Theorem (RMVT) and CUF were adapted to use in multilayered structures embedding FG layers in [17,18]. Talha and Singh [19] investigated the free vibration and static analysis of functionally graded plates using

546

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553

the finite element method by employing a quasi-3D higher order shear deformation theory. Carrera et al. [5] studied the effects of thickness stretching in FG plates and shells. The importance of the transverse normal strain effects in mechanical prediction of stresses for FG plates was pointed out. In fact, this work is an extension of other several FGM papers published by using Carrera’s Unified Formulation (CUF) [16–18,20,21]. Other refined higher-order models reproduced by CUF for example adopting the Mixed Interpolation of Tensorial Components (MITCs) finite element scheme for the approximation of FGPs was given by Cinefra et al. [22–24]. In the context of generalized or unified formulations in classical composites, it is important to notice the work done by Soldatos [25], Carrera (CUF) [26] and Demasi (GUF) [27,28]. As mentioned, Carrera et al. [16], extended the Carrera’s Unified Formulation to advanced composites plates. Besides the powerful CUF there exists a generalized formulation proposed by Zenkour [29], which were extended to cover the stretching effect in Zenkour [30]. Also the shear deformation theory proposed by Matusanga [31–33] based on polynomial shape strain functions, needs to be mentioned in this group. The generalized shear deformation theory of FGP presented by Zenkour [29,30] and the one proposed here are similar to the one formulated by Soldatos [25] for laminated composites. Normally non-polynomial shear strain shape functions, such as trigonometric, trigonometric hyperbolic, exponential, etc., can be used in this type of generalized formulation. However, new shear strain shape functions are introduced in this paper. Recently, Neves et al. [34,35] and Ferreira et al. [36] developed a quasi-3D hybrid shear deformation theory for the static and free vibration analysis of functionally graded plates by using collocation with radial basis functions. Their formulation can be seen as a generalization of the original CUF, by introducing different nonpolynomial displacement fields for in-plane displacements and polynomial displacement field for the out-of-plane displacement. Mantari and Guedes Soares [6,37] presented bending results of FGP by using a new non-polynomial HSDT. In [37] the stretching effect was included, while in Mantari and Guedes Soares [38] different non-polynomial shear strain shape functions were considered in this context. Therefore, shape strain functions need to be further explored because they are more effective than polynomial ones [5,39]. In the present paper, a set of new HSDTs not including the stretching effect is discussed in a generalized way by using only non-polynomial HSDTs. The present generalized theory allows the inclusion of an odd shear strain shape function, f(z). Then, adequate distribution of the inplane and transverse shear strains and the tangential stress-free boundary conditions is guaranteed, and thus a shear correction factor is not required. The generalized finite element code is based on a continuous isoparametric Lagrangian finite element with seven degrees of freedom per node. Numerical results for different side-to-thickness ratio, aspect ratios, volume fraction, and simply supported boundary conditions are compared. Results show that some of the new HSDT solved by the proposed generalized finite element technique are more accurate than, for example, the well-known trigonometric HSDT, having the same DOFs. It is concluded that some non-polynomial shear strain shape functions are richer than polynomials counterparts. 2. Generalized displacement field A FGP of uniform thickness h is shown in Fig. 1. The rectangular Cartesian coordinate system x, y, z, has the plane z = 0, coinciding with the mid-surface of the plate. The material is inhomogeneous and the material properties vary through the thickness with a simple power-law distribution, which is given below:

Fig. 1. Geometry of a functionally graded plate.

PðzÞ ¼ ðPt  Pb ÞgðzÞ þ Pb ;  gðzÞ ¼

z 1 þ h 2

ð1aÞ

n ;

ð1bÞ

where P denotes the effective material property, Pt and Pb denote the property of the top and bottom faces of the panel, respectively, and n is the power-law exponent that specifies the material variation profile through the thickness. The effective material properties of the shell, including Young’s modulus, E and shear modulus, G vary according to Eq. (1a), and Poisson ratio, m is assumed to be constant. It is important to note that n is a parameter that dictates the material variation profile through the plate thickness and takes values greater than zero. Recently developed HSDTs [40–43], in which the displacement of the middle surface is expanded as an odd non-polynomial function of the thickness coordinate and the transverse displacement taken to be constant through the thickness is considered in a generalized displacement field, as it is given below:

@w  þ f ðzÞh1 ; @x @w v ðx; y; zÞ ¼ v ðx; yÞ þ z½y h2   þ f ðzÞh2 ; @y  y; zÞ ¼ w; wðx;

 ðx; y; zÞ ¼ uðx; yÞ þ z½y h1  u

ð2a-cÞ

where u(x, y), v(x, y), w(x, y), h1(x, y) and h2(x, y) are the five unknown displacement functions of middle surface of the panel,  whilst y ¼ f 0 2h . As example, for the well-known trigonometric plate theory (TPT) [44], the shape strain function is   f ðzÞ ¼ ph sin ph z and y ¼ 0. Inthe case of the tangential trigonomet2 z ric HSDT [42], f ðzÞ ¼ tan 5h and y ¼  sec 5hð0:1Þ), for other shear strain shape functions, see Table 1. In the derivation of the necessary equations, small strains are assumed (i.e., displacements and rotations are small, and obey Hooke’s law). The linear strain expressions derived from the displacement model of Eqs. (2a-c), valid for thin, moderately thick and thick plate under consideration are as follows:

exx ¼ e0xx þ ze1xx þ f ðzÞe2xx ; eyy ¼ e0yy þ ze1yy þ f ðzÞe2yy ; exy ¼ e0xy þ ze1xy þ f ðzÞe2xy ; exz ¼ e0xz þ f 0 ðzÞe3xz ; eyz ¼ e0yz þ f 0 ðzÞe3yz ;

ð3a-eÞ

547

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553 Table 1 Shear strain shape functions for the generalized HSDT. HSDT

y

f(z)

Reddy [8] Touratier [44] and Zenkour [29], Soldatos [25] Karama et al. [39]

f(z) = z

3

f ðzÞ ¼ zm2ðz=hÞ þ y z

Mantari et al. [40]

f ðzÞ ¼ sinðphzÞ em cosð h Þ þ y z   z f ðzÞ ¼ sinh hz em cosh ðhÞ þ y z

@h1 @ w  2; @x @x @h @2w 2 ¼ y  2; @y @y

e2xx ¼

e0yy

e1yy

e2yy

@ v @u ; þ @x @y

e1xy ¼ y

@h2 @h1 @2w ; þ y 2 @x@y @x @y

e2xy ¼

@h2 @h1 þ ; @x @y

e0xz ¼ y h1 e3xz ¼ h1 ; e0yz ¼ y h2 e3yz ¼ h2 : ð4a-mÞ The linear constitutive relations are given below:

8 9 > < rxx > = > :

2

Q 11

ryy ¼ 6 4 Q 12 > sxy ;b Q 16    sxz Q 55 ¼ syz s Q 45

Q 12 Q 22 Q 26 Q 54 Q 44

9 38 Q 16 > < exx > = 7 Q 26 5 eyy ; > :c > ; Q 66 xy b ) (

cxz cyx

EðzÞ ; 1  m2

ð5a-bÞ

; s

Q 12 ¼ m Q 11 ;

G ¼ Q 44 ¼ Q 55 ¼ Q 66 ¼

EðzÞ : 2ð1 þ mÞ

y ¼ mhp y ¼

2 1 2

cosh ð12Þþm sinh h

ðÞ

em cosh ð2Þ 1

Bb ¼ ½Bm  þ z½Bb  þ f ðzÞ½Bb ;

ð9a-bÞ

Bs ¼ ½Bs  þ f 0 ðzÞ½Bs ;

where matrices of Bb and Bs (see Eqs. (9a-b)), for the example given in Section 4, are given in Appendix A. The potential energy Pe of the plate with a mid-surface area X and volume V, with a load vector P corresponding to the sevendegrees-of-freedom of a point on the mid-plane, can be represented as, Y

in which, r = {rxx, ryy, sxy, sxz, syz}T and e = {exx, eyy, cxy, cxz, cyz}T are the stresses and the strain vectors with respect to the plate coordinate system. The Qij expressions in terms of engineering constants are given below:

Q 11 ¼ Q 22 ¼

2

where [B] is the strain-displacement matrix in the Cartesian coordinate system. The [B] matrix can be divided in two parts, one which contains the bending terms and other containing the shear terms, as follows:

@h1 ; @x @h2 ¼ ; @y

e1xx ¼ y

mh

m

pz

where

e0xy ¼

y ¼ m sec2 pffiffiffim y ¼ 1ln

2

Present

@u ; @x @v ¼ ; @y

y = 0

Mantari et al. [41]

e0xx ¼

y = cosh(1/2)

2

f ðzÞ ¼ ze2ðz=hÞ f(z) = tan mz + yz

Mantari et al. [42]

2

y = 4/(3h2) y = 0

  f ðzÞ ¼ ph sin ph z     f ðzÞ ¼ h sinh hz  z cosh 12

¼ U s  W ext ;

e

Z Z ZZZ 1 1 1 feb gT ½Q bðzÞ feb gdxdydz eT r dV þ cT s dV ¼ 2 V 2 V 2 ZZZ 1 fcs gT ½Q sðzÞ fcs gdxdydz; þ 2 Z W ext ¼ wo qdxdy: Us ¼

ð10a-cÞ

X

Considering the Eqs. (7), the potential energy of an element can be rewritten as follows:

1 1 2 2 1 T ¼ fdg fK e gfdg  fdgT fPe g; 2

Pe ¼ fdgT fK be gfdg þ fdgT fK se gfdg  fdgT fPe g; ð6Þ

The modulus E, G and the elastic coefficients Qij vary through the thickness according to Eqs. (1a) and (1b).

ð11Þ

where

K be ¼ K mm þ K mb þ K bb þ K mb þ K bb þ K b b ; 3. Finite element formulation

and

In the present work, a four-nodded quadrilateral C0 continuous isoparametric element with seven-degrees-of-freedom per node is employed. The generalized displacements included in the present theory can be expressed as follows:

K se ¼ K ss þ K ss þ K s s ;

k X d¼ N i di ;

ð7Þ

where d = {uo, vo, wo, ow/ox, ow/oy, h1, h2}T, di is the displacement vector corresponding to node i, Ni is the shape function associated with the node i and k is the number of nodes per element, which is four in the present study. Considering the Eq. (7), the strain vector {e} can be expressed in terms of d containing nodal degrees of freedom as,

feb g ¼ ½Bb fdg; fcs g ¼ ½Bs fdg;

ð8a-bÞ

ð13Þ

Matrices of Kbe and Kse (see Eqs. (12) and (13)) are given in the following equations

K mm ¼ BTm f ð1ÞBm ;

i¼1

ð12Þ

K mb ¼ BTm f ð2ÞBb þ BTb f ð2ÞBm

K mb ¼ BTm f ð4ÞBb þ BTb f ð4ÞBm ; K b b ¼

; K bb ¼ BTb f ð3ÞBb ;

K bb ¼ BTb f ð5ÞBb þ BTb f ð5ÞBb

BTb f ð6ÞBb ; ð14a - fÞ

K ss ¼ BTs gð1ÞBs ;

K ss ¼ BTs gð2ÞBs þ BTs gð2ÞBs ;

K s s ¼ Bs gð3ÞBs ; f ð1; 2; 3; 4; 5; 6Þ ¼

ð15a-cÞ Z

hk

hk1

 2 Q bðzÞ 1; z; z2 ; fðzÞ ; zfðzÞ ; fðzÞ dz;

ð16Þ

548

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553

gð1; 2; 3Þ ¼

Z

hk

hk1

 0 02 Q sðzÞ 1; fðzÞ ; fðzÞ dz;

ð17Þ

where the expressions Q bðzÞ and Q sðzÞ can be obtained by using Eqs. (1a-b) and (5a-b).

fPge ¼

Z Z

½Nw0 T qdX:

ð18Þ

The equilibrium equations can be obtained by minimizing Pe with respect to {d} as,

½K e fdg ¼ fPe g;

ð19Þ

where [Ke] is the element stiffness matrix and {Pe} is the nodal load vector. In what follows, the problem under consideration is solved for the simply supported boundary conditions (SS3) and they are given at all four edges as follows:

uðx; 0Þ ¼ uðx; bÞ ¼ v ð0; yÞ ¼ v ða; yÞ ¼ 0; wðx; 0Þ ¼ wðx; bÞ ¼ wð0; yÞ ¼ wða; yÞ ¼ 0; @w=@xðx; 0Þ ¼ @w=@xðx; bÞ ¼ @w=@yð0; yÞ ¼ @w=@yða; yÞ ¼ 0; hcx ðx; 0Þ ¼ hcx ðx; bÞ ¼ hcy ð0; yÞ ¼ hcy ða; yÞ ¼ 0: ð20Þ Integrations in Eqs. (10a-c) are carried out numerically by Gauss quadrature integration rule. A selective integration technique is adopted for the calculations of the shear stiffness matrix in order to avoid numerical disturbances such as shear locking, which may appear in a full integration scheme. The stiffness matrix of all the elements is calculated and assembled together to form the overall stiffness matrix [K] of the panel and then the static problem is solved. 4. Numerical results and discussion In this paper a set of non-polynomial HSDTs and finite element formulation for functionally graded plates is presented. However, polynomial HSDTs can be also created. Different HSDTs are obtained or implemented by using the present generalized formulation; see the shear strain shape functions in Table 1. In this paper, only non-polynomial HSDTs are discussed. Notice that Touratier’s HSDT [44] (see also Zenkour [29]) and one developed by Karama et al. [39] can be derived as special cases of the HSDTs formulated by Mantari et al. [40,41] (see Table 1). Moreover, Soldatos’ HSDT can be reproduced as a special case of the new HSDT presented in this paper by doing m = 0 (see Table 1).

The bending problem of this section has been studied by several authors, and their results were used for comparison purposes. Wu and Li [46] solved this example problem by using a RMVT-based third order shear deformation theory (TSDT). Zenkour [29] used a generalized HSDT and the well-known trigonometric plate theory (TPT) [44]. Wu et al. [45] used RMVT-based meshless collocation and element-free Galerkin methods for the quasi-3D analysis of multilayered composite and FG plates. Recently, Mantari et al. [47] also provided results. Results from the presented generalized HSDT, which consider several shear strain shape functions taken from several HSDTs [40–43] are presented. These HSDTs are m parameter dependent, see Table 1. In the computation of each of the new HSDTs, two referential values of m for each HSDT are taken, as can be noticed in Tables 2 and 3. As mentioned above, the HSDTs developed in [40,41,43] and the present one gives as special cases the HSDTs developed by Touratier [44], Karama et al. [39], and Soldatos [25], respectively. Therefore, it can be said that a representative set of non-polynomial shear strain shape functions are considered in this study. Then, the results obtained by using the generalized HSDT (with different shear strain shape functions) solved by finite element technique are compared with those mentioned above, see Table 2. Table 2 presents results of in-plane displacements, normal stress, and shear stresses at the specified positions of the FGP for n = {1, 2, 4, 8}, h = 1 m, and a/h = 10. In general, the results of displacements, normal stress and in-plane shear stress are in good agreement with all the theories compared in this section and particularly with the referential quasi-3D results provided by Wu et al.  , produces a [45], see Table 2. Results of in-plane displacement, u relative error compared with the referential solution which slightly increases as n increases; the maximum relative error reported is about 2.3%, see Present ([40]; m = 0.5) in Table 2. In the case of ver the relative error increases as n increases tical central deflection, w, and then decreases; the maximum relative error reported is about  xx , as in the case of in-plane 0.7%. For the normal stress results, r

0.5

0.25

4.1. Bending analysis of FGM Al/Al2O3 plates The results are presented for the simply supported plate under bi-sinusoidal transverse loads of intensity q. The static analysis was conducted using aluminium (bottom, Al) and alumina (top, Al2O3). The following material properties are used for computing the numerical results [45]:

Et ¼ 380 GPa;v t ¼ 0:3; Eb ¼ 70 GPa;v b ¼ 0:3:

_ z

0

-0.25

ð21Þ

The following non-dimensional quantities are used:

    3 3 100h Et b 10h Et a b  ¼ u 0; ; z ; wðzÞ w ; ;z ; 4 4 2 2 2 qa qa   h a b h  xy ðzÞ ¼ rxy ð0; 0; zÞ; ; ;z ; r r xx ðzÞ ¼ rxx qa 2 2 qa   h b z 0; ; z ; z ¼ ; r xz ðzÞ ¼ rxz qa 2 h

-0.5

 ðzÞ ¼ u

0

0.2

n=0.02 n=0.1

0.4

_ g (z)

n=0.2 n=0.5

0.6

n=1.0 n=2.0

0.8

1

n=10 n=300

ð22a-fÞ where q denotes the transverse load. Fig. 2 presents the through thickness distribution of the volume fraction gðzÞ for different values of n.

Fig. 2. Variation of volume fraction gðzÞ through panel thickness for different values of the power-law index, n.

549

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553 Table 2 Displacement and stresses of FGP under bi-sinusoidal load of a thick (a/h = 10) FGP. n

Theories

 ðh=4Þ u

 ð0Þ w

r xxðh=3Þ

bar rxyðh=3Þ

r Axzðh=6Þ

r Bxzðh=6Þ

1

Present (m = 7, N = 21) Present (m = 0, N = 21) Present ([42]; m = 1/(5 h), N = 21) Present ([42]; m = p/(2 h), N = 21) Present ([40]; m = 0.5, N = 21) Present ([40]; m = 0, N = 21) Present ([41]; m = 2.85, N = 21) Present ([41]; m = e, N = 21) HSDT ([41]; m = 2.85) HSDT [47] Generalized TSDT [29] RMVT-Based TSDT [46] RMVT-Based Collocation [45] RMVT-Based Galerkin [45]

0.6324 0.6328 0.6328 0.6332 0.6324 0.6327 0.6325 0.6326 0.6406 0.6398 0.6626 0.6414 0.6436 0.6436

0.5849 0.5851 0.5851 0.5850 0.5848 0.5850 0.5850 0.5850 0.5888 0.5880 0.5889 0.5890 0.5876 0.5876

1.4857 1.4851 1.4851 1.4852 1.4861 1.4852 1.4854 1.4854 1.4891 1.4888 1.4894 1.4898 1.5062 1.5061

0.6059 0.6058 0.6058 0.6058 0.6060 0.6058 0.6059 0.6058 0.6110 0.6109 0.6110 0.6111 0.6112 0.6112

0.2572 0.2599 0.2599 0.2580 0.2542 0.2598 0.2589 0.2590 0.2624 0.2566 0.2622 0.2506 0.2509 0.2511

0.2084 0.2605 0.2596 0.2898 0.1915 0.2441 0.2277 0.2293 0.2624 0.2566 0.2622 0.2506 0.2509 0.2511

2

Present (m = 7, N = 21) Present (m = 0, N = 21) Present ([42]; m = 1/(5 h), N = 21) Present ([42]; m = p/(2 h), N = 21) Present ([40]; m = 0.5, N = 21) Present ([40]; m = 0, N = 21) Present ([41]; m = 2.85, N = 21) Present ([41]; m = e, N = 21) HSDT ([41]; m = 2.85) HSDT [47] Generalized TSDT [29] RMVT-Based TSDT [46] RMVT-Based Collocation [45] RMVT-Based Galerkin [45]

0.8853 0.8861 0.8861 0.8868 0.8852 0.8858 0.8856 0.8856 0.8971 0.8957 0.9281 0.8984 0.9015 0.9013

0.7522 0.7523 0.7523 0.7521 0.7520 0.7523 0.7523 0.7523 0.7572 0.7564 0.7573 0.7573 0.7572 0.7571

1.3915 1.3912 1.3912 1.3916 1.3919 1.3911 1.3912 1.3912 1.3948 1.3940 1.3954 1.3960 1.4129 1.4133

0.5394 0.5394 0.5394 0.5395 0.5395 0.5394 0.5394 0.5394 0.5440 0.5438 0.5441 0.5442 0.5437 0.5436

0.2734 0.2744 0.2745 0.2707 0.2704 0.2751 0.2748 0.2748 0.2778 0.2741 0.2763 0.2491 0.2495 0.2495

0.2366 0.2872 0.2864 0.3138 0.2196 0.2718 0.2560 0.2575 0.2778 0.2741 0.2763 0.2491 0.2495 0.2495

4

Present (m = 7, N = 21) Present (m = 0, N = 21) Present ([42]; m = 1/(5 h), N = 21) Present ([42]; m = p/(2 h), N = 21) Present ([40]; m = 0.5, N = 21) Present ([40]; m = 0, N = 21) Present ([41]; m = 2.85, N = 21) Present ([41]; m = e, N = 21) HSDT ([41]; m = 2.85) HSDT [47] Generalized TSDT [29] RMVT-Based TSDT [46] RMVT-Based Collocation [45] RMVT-Based Galerkin [45]

1.0335 1.0352 1.0352 1.0365 1.0331 1.0346 1.0340 1.0341 1.0481 1.0457 1.0941 1.0502 1.0548 1.0541

0.8763 0.8760 0.8759 0.8753 0.8761 0.8762 0.8763 0.8763 0.8820 0.8814 0.8819 0.8815 0.8826 0.8823

1.1734 1.1739 1.1740 1.1750 1.1737 1.1736 1.1734 1.1734 1.1772 1.1755 1.1783 1.1794 1.1935 1.1941

0.5616 0.5617 0.5617 0.5619 0.5616 0.5616 0.5616 0.5616 0.5665 0.5662 0.5667 0.5669 0.5674 0.5671

0.2610 0.2581 0.2577 0.2516 0.2588 0.2601 0.2611 0.2610 0.2614 0.2623 0.2580 0.2360 0.2360 0.2362

0.2470 0.2875 0.2855 0.3062 0.2323 0.2758 0.2632 0.2644 0.2614 0.2623 0.2580 0.2360 0.2360 0.2362

8

Present (m = 7, N = 21) Present (m = 0, N = 21) Present ([42]; m = 1/(5 h), N = 21) Present ([42]; m = p/(2 h), N = 21) Present ([40]; m = 0.5, N = 21) Present ([40]; m = 0, N = 21) Present ([41]; m = 2.85, N = 21) Present ([41]; m = e, N = 21) HSDT ([41]; m = 2.85) HSDT [47] Generalized TSDT [29] RMVT-Based TSDT [46] RMVT-Based Collocation [45] RMVT-Based Galerkin [45]

1.0588 1.0609 1.0609 1.0625 1.0584 1.0602 1.0595 1.0595 1.0737 1.0709 1.1340 1.0763 1.0840 1.0830

0.9686 0.9686 0.9686 0.9678 0.9682 0.9688 0.9688 0.9688 0.9749 0.9737 0.9750 0.9747 0.9727 0.9739

0.9424 0.9432 0.9431 0.9444 0.9425 0.9427 0.9425 0.9425 0.9451 0.9431 0.9466 0.9477 0.9568 0.9622

0.5803 0.5804 0.5804 0.5806 0.5804 0.5803 0.5803 0.5803 0.5853 0.5850 0.5856 0.5858 0.5886 0.5883

0.2134 0.2120 0.2120 0.2069 0.2112 0.2134 0.2139 0.2139 0.2145 0.2140 0.2121 0.2263 0.2251 0.2261

0.2006 0.2331 0.2326 0.2479 0.1886 0.2238 0.2137 0.2146 0.2145 0.2140 0.2121 0.2263 0.2251 0.2261

stresses, the relative error increases as n increases; the maximum relative error reported is about 2%. In the case of in-plane shear  xy , the maximum error stresses results vertical central deflection, r reported is 1.4%. Because a representative difference in results are reported for shear stresses (due to the fact that different shear strain shape functions are used), as also shown in Mantari et al. [47], then, discussion  xz , will be detailed in what follows. of the shear stresses results, r It is important to stress out that the present generalized formulations just include five unknowns, but the generalized finite element formulation assumes 7DOFs because it is a non-conform Axz ing type. With this assumption (see Eq. (7)), the shear stresses r are obtained but the surface boundary conditions are not exactly satisfied, although considerably well-modelled (see Fig. 6).

However, the free surface boundary condition can be complied  Bxz , with in the post-process computation of the shear stresses, r

 @w FEM Þ by assuming that @ðx;yÞ  @ðw (see Fig. 7). The original shear @ðx;yÞ FEM

stresses results and the ones with such additional assumption  Bxz ) are presented in Figs. 6 and 7 and Tables 2 and 3. The  Axz ; r (r superscript ‘‘A’’ means original results (these FEM shear stresses results, keeping the 7DOFs in the entire computation, which do not exactly comply the free surface boundary conditions, as seen in Fig. 6), and the superscript ‘‘B’’ means that shear stresses results consider the above mentioned additional assumption included in the post-process computation in order to comply exactly with the free surface boundary conditions, as shown in Fig. 7.

550

Table 3 Vertical deflection and stresses of FGP under bi-sinusoidal load of a thick (a/h = {4, 10}) FGP for n = {0, 0.5, 1, 1.5, 2}. Theory

n=0

n = 0.5

n=1

n=2

 w

r xx

 Axz

r

 Bxz

r

 w

r xx

 Axz

r

 Bxz

r

 w

r xx

 Axz

r

 Bxz

r

 w

r xx

r Axz

r Bxz

4

Present (m = 7, N = 21) Present (m = 0, N = 21) Present ([42]; m = 1/(5 h), N = 21) Present ([42]; m = p/(2 h), N = 21) Present ([40]; m = 0.5, N = 21) Present ([40]; m = 0, N = 21) Present ([41]; m = 2.85, N = 21) Present ([41]; m = e, N = 21) HSDT ([42]; m = 1/(5 h)) HSDT [47] HSDT ([41]; m = 2.85) HSDT [8]

0.0034 0.0034 0.0034 0.0034 0.0034 0.0034 0.0034 0.0034 0.0035 0.0034 0.0035 0.0035

0.8314 0.8332 0.8332 0.8327 0.8302 0.8330 0.8324 0.8325 0.8365 0.8479 0.8424 0.8365

0.2481 0.2357 0.2360 0.2269 0.2520 0.2402 0.2442 0.2439 0.2375 0.2708 0.2532 0.2377

0.2053 0.2329 0.2324 0.2496 0.1981 0.2241 0.2156 0.2164 – – – –

0.0035 0.0035 0.0035 0.0035 0.0035 0.0035 0.0035 0.0035 0.0035 0.0035 0.0035 0.0035

0.9723 0.9747 0.9748 0.9743 0.9708 0.9744 0.9736 0.9737 0.9828 0.9959 0.9888 0.9789

0.2508 0.2386 0.2389 0.2299 0.2546 0.2429 0.2469 0.2466 0.2406 0.2744 0.2564 0.2408

0.2058 0.2348 0.2344 0.2526 0.1983 0.2255 0.2166 0.2174 – – – –

0.0036 0.0036 0.0036 0.0036 0.0036 0.0036 0.0036 0.0036 0.0036 0.0036 0.0036 0.0036

1.0489 1.0514 1.0514 1.0507 1.0473 1.0510 1.0503 1.0504 1.0555 1.0712 1.0637 1.0556

0.2481 0.2357 0.2360 0.2269 0.2520 0.2402 0.2442 0.2439 0.2375 0.2708 0.2532 0.2377

0.2053 0.2329 0.2324 0.2496 0.1981 0.2241 0.2156 0.2164 – – – –

0.0037 0.0037 0.0037 0.0037 0.0037 0.0037 0.0037 0.0037 0.0037 0.0037 0.0037 0.0037

1.1290 1.1312 1.1312 1.1300 1.1273 1.1311 1.1304 1.1305 1.1337 1.1518 1.1431 1.1338

0.2417 0.2282 0.2285 0.2188 0.2460 0.2330 0.2374 0.2370 0.2287 0.2628 0.2447 0.2289

0.2063 0.2307 0.2303 0.2450 0.2001 0.2231 0.2156 0.2163 – – – –

10

Present (m = 7, N = 21) Present (m = 0, N = 21) Present ([42]; m = 1/(5 h), N = 21) Present ([42]; m = p/(2 h), N = 21) Present ([40]; m = 0.5, N = 21) Present ([40]; m = 0, N = 21) Present ([41]; m = 2.85, N = 21) Present ([41]; m = e, N = 21) HSDT ([42]; m = 1/(5 h)) HSDT [47] HSDT ([41]; m = 2.85) HSDT [8]

0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027

1.9875 1.9881 1.9881 1.9879 1.9870 1.9881 1.9879 1.9879 1.9943 1.9990 1.9967 1.9943

0.2505 0.2373 0.2376 0.2281 0.2547 0.2420 0.2463 0.2460 0.2383 0.2730 0.2546 0.2386

0.2108 0.2378 0.2373 0.2539 0.2038 0.2293 0.2210 0.2217 – – – –

0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0027 0.0028 0.0028 0.0028 0.0028

2.3207 2.3216 2.3218 2.3214 2.3202 2.3215 2.3212 2.3213 2.3530 2.3399 2.3371 2.3289

0.2529 0.2402 0.2411 0.2311 0.2572 0.2448 0.2490 0.2486 0.2414 0.2762 0.2578 0.2417

0.2109 0.2397 0.2413 0.2570 0.2039 0.2307 0.2219 0.2227 – – – –

0.0028 0.0028 0.0028 0.0028 0.0028 0.0028 0.0028 0.0028 0.0028 0.0028 0.0028 0.0028

2.4979 2.4987 2.4987 2.4984 2.4972 2.4986 2.4984 2.4984 2.5065 2.5129 2.5098 2.5065

0.2505 0.2373 0.2376 0.2281 0.2547 0.2420 0.2463 0.2460 0.2383 0.2730 0.2546 0.2386

0.2108 0.2378 0.2373 0.2539 0.2038 0.2293 0.2210 0.2217 – – – –

0.0029 0.0029 0.0029 0.0029 0.0029 0.0029 0.0029 0.0029 0.0029 0.0029 0.0029 0.0029

2.6748 2.6756 2.6756 2.6750 2.6742 2.6756 2.6753 2.6754 2.6830 2.6905 2.6869 2.6831

0.2442 0.2299 0.2302 0.2200 0.2488 0.2350 0.2396 0.2392 0.2295 0.2651 0.2462 0.2298

0.2119 0.2356 0.2352 0.2492 0.2058 0.2282 0.2210 0.2217 – – – –

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553

a/h

551

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553

 Axz , it can be said that Regarding to the shear stresses results, r the HSDT that produces more accurate results is the Present ([40]; m = 0.5) in Table 2; followed by the Present ([43]; m = 7), Present ([42]; m = p/(2h)) and Present ([41]; m = 2.85); the maximum error reported when n = 1 is 1.2% which increases as n increases, except in the case of n = 8, in which the results are under predicted. In the case of n = 8, the ones that produce closer results is the Present ([41]; m = 2.85), followed by the Present ([43]; m = 7), see Table 2. Another important aspect to notice is that the present FEM results are close to their closed-form solution, as shown for example, the Present HSDT ([40]; m = 0.5) and HSDT [47], Present ([40]; m = 0) and trigonometric plate theory (TPT) [29,44], or the Present ([41]; m = 2.85) and HSDT ([41]; m = 2.85).  Bxz , it can be said that With respect to the shear stresses results, r the HSDT that produce more accurate results are the Present ([40]; m = 0) and the Present ([41]; m = 2.85) in correspondence to both  Axz ) and their closed-form solutions. In order their counterparts (r to have a better or general idea, in what follows, a thicker FGP is studied. Figs. 3–5 show the through thickness distributions of displace ), normal stresses (r  xx ), and in-plane shear stresses ment (u  xy ) at the specific positions for the FGP with power-law in(sigma dex, n = 4, h = 1 m, and a/h = 5, respectively. The results are compared with the quasi-3D results [45]. In general, the results are in excellent agreement with the quasi-3D results. Figs. 6 and 7 show the through thickness distribution of the transverse shear stresses at the specific positions for the FGP with power-law index, n = 4, h = 1 m, and a/h = 5. Fig. 6 shows the distri Axz , while Fig. 7 the distribution of r  Bxz . From Fig. 6, a good bution of r approximation to Quasi-3D solution can be notice, but not all the HSDTs comply with top and bottom free surface boundary conditions. This may be attributed to the assumption in the description of the 7DOFs in the finite element technique, see Eq. (7). Fig. 7 shows such compliance but the results are not in good agreement with the quasi-3D solutions except in the case of the Present HSDT ([41]; m = 2.85) and Present HSDT ([40]; m = 0). Therefore, these HSDTs need to be further explored in future investigations.

0.5 0.4 0.3 0.2 0.1

_ z

0 -0.1 -0.2 Present (m=0) Present ([41], m=2.85) Present ([40], m=0) Present ([42], m=1/(5h)) Quasi-3D

-0.3 -0.4 -0.5 -1.2

-0.6

0

0.6

1.2

_

1.8

2.4

σxx  xx (a/2, b/2, z), through Fig. 4. Distribution of non-dimensionalized normal stress, r the thickness of a thick (a/h = 5) FGP (n = 4).

0.5 0.4 0.3 0.2 0.1

_ z

0 -0.1

0.5

-0.2

Present (m=0) Present ([41]), m=2.85) Present ([40], m=0) Present ([42], m=1/(5h)) Quasi-3D

0.4 0.3

-0.3 -0.4 -0.5 -1.2

0.2 0.1

_ z

Present (m=0) Present ([41], m=2.85) Present ([40], m=0) Present ([42], m=1/(5h)) Quasi-3D -1

-0.8

-0.6

-0.4

_ σxy

-0.2

0

0.2

0.4

 xy (0, 0, z), Fig. 5. Distribution of non-dimensionalized in-plane shear stress, r through the thickness of a thick (a/h = 5) FGP (n = 4).

0 -0.1

4.2. Bending analysis of FGM Al/ZrO2 plates

-0.2

The results are presented for simply supported plates under bisinusoidal transverse loads of intensity q. The static analysis was conducted using Aluminium (bottom, Al) and Zirconia (top, ZrO2). The following material properties are used for computing the numerical results [8].

-0.3 -0.4 -0.5 -2

-1

0

1

_ u

2

3

4

 (0, b/2, z), Fig. 3. Distribution of non-dimensionalized in-plane displacement, u through the thickness of a thick (a/h = 5) FGP (n = 4).

Et ¼ 151 GPa;tt ¼ 0:3; Eb ¼ 70 GPa;tb ¼ 0:3: The following non-dimensional quantities are used:

ð23Þ

552

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553

0.5 0.4 0.3 0.2 0.1

_ z

0 -0.1 -0.2 Present (m=0) Present ([41], m=2.85) Present ([40], m=0) Present ([42], m=1/(5h)) Quasi-3D

-0.3 -0.4 -0.5

0

0.08

_ A σ

0.16

0.24

0.32

xz

 AXZ (0, b/2, z), through the Fig. 6. Distribution of non-dimensionalized shear stress, r thickness of a thick (a/h = 5) FGP (n = 4).

0.5 0.4 0.3 0.2 0.1

_ z

0 -0.1

5. Conclusions

-0.2 Present (m=0) Present ([41], m=2.85) Present ([40], m=0) Present ([42], m=1/(5h)) Quasi-3D

-0.3 -0.4 -0.5 0

0.08

0.16

_ B σ

0.24

0.32

xz

 BXZ (0, b/2, z), through the Fig. 7. Distribution of non-dimensionalized shear stress, r thickness of a thick (a/h = 5) FGP (n = 4).

  D ð2; 2Þ a b ; ; 0 ; w qa4 2 2   h b ¼ rxz 2 0; ; 0 ; qa 2

 ¼ w

r xz

computed as a function of the side-to-thickness-ratio, a/ h = {4, 10}, for different values of power law index, n = {0, 0.5, 1, 2}. In this table, it can be noticed that the central deflection results are well approximated by the finite element method in all the  xx , are also well precases. The non-dimensional axial stresses, r dicted by the FEM and the analytical solution when compared with the well-known HSDT proposed by Reddy [8], with the exception of the closed-form HSDT results provided by Mantari et al. [47],  xz which are slightly higher. The non-dimensional shear stress, r (0, b/2, 0) results, in general, present representative relative difference in results due to the fact that several shear strain shape functions were used, as mentioned in the previous example problem. Particularly, there exist differences in the results of shear stres Axz , between the HSDT solved in closed-form and by using the ses, r finite element technique proposed in this paper, see Present ([40]; m = 0.5) and HSDT [47] in Table 3. However, not all the HSDTs present such representative difference in results, for example the differences between Present HSDT ([42]; m = 1/(5h)) and HSDT ([42]; m = 1/(5h)) are lesser than the ones between Present HSDT ([41]; m = 2.85) and HSDT ([41]; m = 2.85), see the italicized cells in Table 3. Therefore, it can be concluded that the shear strain shape function plays an important role in the prediction of the shear stresses, and in order to obtain an appropriate shear strain shape function, an exhaustive comparison analysis with other example problems need to be performed. So far it is possible to notice that there exist other nonpolynomial HSDTs which produce more accurate results than the  Bxz , again well-known trigonometric HSDT [29]. In the case of, r  Axz and its closed-form the HSDT which produce close results to r solution is the Present ([42]; m = 1/(5h)), which is also comparable with the results presented by Reddy [8]. In summary it will interesting to perform further computations with different new shear strain shape functions, also include the transverse expansion of the plate, and then provide general conclu Axz and r  Bxz . For more details of this aspect, sion regarding to both r readers may consult to the interesting work developed by Carrera et al. [5], in which the importance of the stretching effect was stressed out.

r xx ¼ rxx

  h a b h ; ; ; qa2 2 2 2 ð24a-cÞ

where q denotes the transverse load. Results from the presented generalized HSDT, which consider several shear strain shapes functions taken from several HSDTs mentioned above are presented and compared with analytical closed-form solution provided by Reddy [8], Mantari et al. [47], and the results of closed-form solutions of the HSDTs [41,42] adapted to FGPs in this paper. Table 3 presents results for non-dimensionalized vertical  normal stresses, r  xx , and shear stresses, r  xz , displacement, w,

A generalized higher order shear deformation theory (HSDT) and its non-conforming finite element solution for the bending analysis of functionally graded plates (FGPs) are presented. The generalized HSDT accounts for non-linear and constant variation of in-plane and transverse displacement respectively through the plate thickness, complies with plate surface boundary conditions and do not require shear correction factors. The generalized finite element code is based on a continuous isoparametric Lagrangian finite element with seven degrees of freedom per node. Numerical results for different side-to-thickness ratio, volume fraction, and simply supported boundary conditions are compared. Results show that new non-polynomial HSDTs solved by the proposed generalized finite element technique are more accurate than, for example, the wellknown trigonometric HSDT, having the same DOFs. It is concluded that some non-polynomial shear strain shape functions are more effective than polynomials counterparts. Benchmark results for the displacement and stresses of functionally graded plates are obtained, which can be used for the evaluation of other HSDTs solved in closed-form or by other numerical method such as meshless. Acknowledgment The first author has been financed by the Portuguese Foundation of Science and Technology under the contract number SFRH/ BD/66847/2009.

J.L. Mantari, C. Guedes Soares / Composite Structures 96 (2013) 545–553

Appendix A. Definition of matrices given in Eq. (9a-b)

2

@N @x

6 0 Bm ¼ 6 4

@N @y

0

0 0 0 0 0 ...

7 0 0 0 0 0 . . . 7; 5 0 0 0 0 0 ...

@N @y @N @x

2

0 0 0  @N @x 6 6 Bb ¼ 4 0 0 0 0 0 0 0  @N @y 2

3

3

0

C @N @x

 @N @y

0

 @N @x

C @N @y

7 C @N . . . 7; @y 5 C @N ... @x

0 ...

3

@N @x

0 0 0 0 0

6 0 0 0 0 0 0 Bb ¼ 6 4 0 0 0 0 0 @N @y

@N @y @N @x

0 ...

7 . . . 7; 5 ...

ðA1a-cÞ

and

" Bs ¼

Bs ¼

0 0 0 0



@N @x @N @y

N

0

CN

0 ...

0

N

0

CN . . .

0 0 0 0 0 N

0 ...

0 0 0 0 0 0

N ...

# ;

;

ðA2a-bÞ

where

C ¼ P1x ¼ P 1y ¼ f 0 ðzÞjh : 2

ðA3Þ

References [1] Koizumi M. The concept of FGM, ceramic transactions. Funct Grad Mater 1993;34:3–10. [2] Wetherhold RC, Wang SS. The use of functionally graded materials to eliminate or control thermal deformation. Compos Sci Technol 1996;56:1099–104. [3] Vel SS, Batra RC. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. J Sound Vib 2004;272:703–30. [4] Birman V, Byrd LW. Modeling and analysis of functionally graded materials and structures. ASME Appl Mech Rev 2007;60:195–216. [5] Carrera E, Brischetto S, Cinefra M, Soave M. Effects of thickness stretching in functionally graded plates and shells. Composites: Part B 2011;42:123–33. [6] Mantari JL, Guedes Soares C. Bending analysis of thick exponentially graded plates using a new trigonometric higher order shear deformation theory. Compos Struct 2012;94(6):1991–2000. [7] Reddy JN, Chin CD. Thermomechanical analysis of functionally graded cylinders and plates. J Thermal Stress 1998;21:593–626. [8] Reddy JN. Analysis of functionally graded plates. Int J Numer Methods Eng 2000;47:663–84. [9] Tanaka K, Tanaka Y, Watanabe H, Poterasu VF, Sugano Y. An improved solution to thermoelastic materials designed in functionally graded materials: scheme to reduce thermal stresses. Comput Methods Appl Mech Eng 1993;106:377–89. [10] Erdogan F. Fracture mechanics of functionally graded materials. Compos Eng 1995;5:770–3. [11] Jin JH, Batra RC. Some basic fracture mechanics concepts of functionally graded materials. J Mech Phys Solids 1996;44:1221–35. [12] He XQ, Liew KM, Ng TY, Sivashanker S. Active control of FGM plates subjected to a temperature gradient: modeling via finite element method based on FSDT. Int J Numer Methods Eng 2001;52:1253–71. [13] Liew KM, He XQ, Kitipornchai S. Finite element method for the feedback control of FGM shells in the frequency domain via piezoelectric sensors and actuators. Comput Methods Appl Mech Eng 2004;193:257–73. [14] Chinosi C, Della Croce L. Mixed-interpolated finite elements for functionally graded cylindrical shells. Compos Struct 2010;92:2314–20. [15] Della Croce L, Venini P. Finite elements for functionally graded ReissnerMindlin plates. Comput Methods Appl Mech Eng 2004;193:705–25. [16] Carrera E, Brischetto S, Robaldo A. A variable kinematic model for the analysis of Functionally Graded Materials plates. AIAA J 2008;46:194–203. [17] Brischetto S, Carrera E. Advanced mixed theories for bending analysis of functionallygradedplates. Comput Struct 2008;88:1474–83.

553

[18] Carrera E, Brischetto S, Cinefra M, Soave M. Refined and advanced models for multilayered plates and shells embedding functionally graded material layers. Mech Adv Mater Struct 2010;19:604–22. [19] Talha M, Singh BN. Static response and free vibration analysis of fgm plates using higher order shear deformation theory. Appl Math Model 2010;34:3991–4011. [20] Brischetto S, Leetsch R, Carrera E, Wallmersrerger T, Kröplin B. Thermomachanical bending of functionally graded plates. J. Thermal Stress 2008;31(3):286–308. [21] Brischetto S. Classical and mixed advanced models for sandwich plates embedding functionally graded cores. J Mech Mater Struct 2009;4:13–33. [22] Cinefra M, Carrera E, Brischetto S, Belouettar S. Thermo-mechanical analysis of functionally graded shells. J. Thermal Stress 2010;33:942–63. [23] Cinefra M, Belouettar S, Soave M, Carrera E. Variable kinematic models applied to free vibration analysis of functionally graded material shells. Eur J Mech A Solids 2010;29:1078–87. [24] Cinefra M, Soave M. Accurate vibration analysis of multilayered plates made of functionally graded materials. Mech Adv Mater Struct 2011;18(1):3–13. [25] Soldatos KP. A transverse shear deformation theory for homogeneous monoclinic plates. Acta Mech 1992;94:195–220. [26] Carrera E. Theories and finite elements for multilayered plates and shells: a unified compact formulation with numerical assessment and benchmarks. Arch Comput Meth Eng 2003;10:215–96. [27] Demasi L. 13 Hierarchy plate theories for thick and thin composite plates: The generalized unified formulation. Compos Struct 2008;84:256–70. [28] Demasi L. 16 Mixed plate theories based on the generalized unified formulation. Part I: Governing equations. Compos Struct 2009;87:1–11. [29] Zenkour AM. Generalized shear deformation theory for bending analysis of functionally graded plates. Appl Math Model 2006;30:67–84. [30] Zenkour AM. Benchmark trigonometric and 3-D elasticity solutions for an exponentially graded thick rectangular plate. Appl Math Model 2007;77:197–214. [31] Matsunaga H. Free vibration and stability of functionally graded plates according to 2-D higher-order deformation theory. Compos Struct 2008;82:256–70. [32] Matsunaga H. Stress analysis of functionally graded plates subjected to thermal and mechanical loadings. Compos Struct 2009;87:344–57. [33] Matsunaga H. Thermal buckling of functionally graded plates according to a 2D higher-order deformation theory. Compos Struct 2009;90:76–86. [34] Neves AMA, Ferreira AJM, Carrera E, Roque CMC, Cinefra M, Jorge RMN, et al. Bending of FGM plates by a sinusoidal plate formulation and collocation with radial basis functions. Mech Res Commun 2011;38:368–71. [35] Neves AMA, Ferreira AJM, Carrera E, Roque CMC, Cinefra M, Jorge RMN, et al. A quasi-3D sinusoidal shear deformation theory for the static and free vibration analysis of functionally graded plates. Composites: Part B 2012;43:711–25. [36] Ferreira AJM, Carrera E, Cinefra M, Roque CMC, Polit O. Analysis of laminated shells by a sinusoidal shear deformation theory and radial basis functions collocation, accounting for through-the-thickness deformations. Composites: Part B 2011;42:276–84. [37] Mantari JL, Guedes Soares C. A novel higher-order shear deformation theory with stretching effect for functionally graded plates. Composites: Part B (2012). http://dx.doi.org/10.1016/j.compositesb. 2012.05.036. [38] Mantari JL, Guedes Soares C. Generalized hybrid quasi-3D shear deformation theory for the static analysis of advanced composite plates. Compos Struct 2012;94(8):2561–75. [39] Karama M, Afaq KS, Mistou S. Mechanical behavior of laminated composite beam by the new multilayered laminated Composite Structures model with transverse shear stress continuity. Int J Solid Struct 2003;40(6):1525–46. [40] Mantari JL, Oktem AS, Guedes Soares C. A new higher order shear deformation theory for sandwich and composite laminated plates. Composites: Part B 2012;43:1489–99. [41] Mantari JL, Oktem AS, Guedes Soares C. Static and dynamic analysis of laminated composite and sandwich plates and shells by using a new higherorder shear deformation theory. Compos Struct 2011;94:37–49. [42] Mantari JL, Oktem AS, Guedes Soares C. A new trigonometric shear deformation theory for isotropic, laminated composite and sandwich plates. Int J Solids Struct 2012;49:43–53. [43] Mantari JL, Guedes Soares C. Analysis of isotropic and multilayered plates and shells by using a generalized higher-order shear deformation theory. Compos Struct 2012;94(8):2640–56. [44] Touratier M. An efficient standard plate theory. Int J Eng Sci 1991;29(8):901–16. [45] Wu CP, Chiu KH, Wang YM. RMVT-based meshless collocation and elementfree Galerkin methods for the quasi-3D analysis of multilayered composite and FGM plates. Compos Struct 2011;93:923–43. [46] Wu CP, Li HY. An RMVT-based third-order shear deformation theory of multilayered functionally graded material plates. Compos Struct 2010;92:2591–605. [47] Mantari JL, Oktem AS, Guedes Soares C. Bending response of functionally graded plates by using a new higher order shear deformation theory. Compos Struct 2012;94:714–23.