An improved Moving Kriging-based meshfree method for static, dynamic and buckling analyses of functionally graded isotropic and sandwich plates

Engineering Analysis with Boundary Elements 64 (2016) 122–136

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

An improved Moving Kriging-based meshfree method for static, dynamic and buckling analyses of functionally graded isotropic and sandwich plates Chien H. Thai a,b, Vuong N.V. Do b, H. Nguyen-Xuan c,d,n a

Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Vietnam Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Vietnam c Department of Physical Therapy, Graduate Institute of Rehabilitation Science, China Medical University, Taichung 40402, Taiwan d Center for Interdisciplinary Research in Technology, Ho Chi Minh City University of Technology, Vietnam b

art ic l e i nf o

a b s t r a c t

Article history: Received 18 June 2015 Received in revised form 1 December 2015 Accepted 8 December 2015

A meshfree method with a modified distribution function of Moving Kriging (MK) interpolation is investigated. This method is then combined with a high order shear deformation theory (HSDT) for static, dynamic and buckling analyses of functionally graded material (FGM) isotropic and sandwich plates. A meshfree method uses the normalized form of MK interpolation under a new quartic polynomial correlation to build the basis shape functions in high order approximations. The Galerkin weak form is used to separate the system equations which is numerically solved by meshfree method. A rotation-free technique extracted from isogeometric analysis is introduced to eliminate the degrees of freedom of slopes. Then, the method retains a highly computational effect with a lower number of degrees of freedom. In addition, the requirement of shear correction factors is ignored and the traction free is at the top and bottom surfaces of FGM plates. Various thickness ratios, boundary conditions and material properties are studied to validate the numerical analyses of the rectangular and circular plates. The numerical results show that the present theory is more stable and well accurate prediction as compared to three-dimensional (3D) elasticity solution and other meshfree methods in the literature. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Sandwich plate structures Functionally graded material (FGM) Meshfree methods Modified moving Kriging

1. Introduction Functionally graded materials (FGMs) were first introduced by a group of Japanese scientists in the last two decades [1,2]. Until now, FGMs have been developed and widely utilized in various engineering fields. FGMs have played an important role and made preferable many disciplines such as shipbuilding, automotive engineering, nuclear power planning and aero-space, medical devices, etc. In this research, the concerned FG material is used to study since the gradual change of its properties can be tailored to many applications and working environments. In the typical FGM plates, the material is continuously changed in microstructure from ceramic in top layer and metal in bottom one. It is known that the ceramic with a low conductivity can well resist to the effects of thermal stress due to the high temperature environment as well as the surface corrosion of components [3–5], meanwhile n Corresponding author at: Center for Interdisciplinary Research in Technology, Ho Chi Minh City University of Technology, Vietnam. E-mail addresses: [email protected] (C.H. Thai), [email protected] (V.N.V. Do), [email protected], [email protected] (H. Nguyen-Xuan).

http://dx.doi.org/10.1016/j.enganabound.2015.12.003 0955-7997/& 2015 Elsevier Ltd. All rights reserved.

the metal is well capable of imposing high loading intensity to the structures. Research studies for the thermal stress analysis of FGMs were reported in [6,7]. Other applications of FGMs can be found in electronic fields [8,9]. In addition, the transient waves by exciting impact loads in plates are very promising in terms of nondestructive evaluation and material characterization [10–12], etc. To use these materials effectively, static, dynamic and buckling analyses are really necessary and therefore the research task has never been ended. Besides various applications of FGM plates to engineering problems, many plate theories have been proposed such as the classical plate theory (CPT) [13,14] with Kirchhoff–Love assumptions for thin plate. The first order shear deformation (FSDT) [15–17] which included shear effects was then proposed for Reissner–Mindlin plates. However, with the assumption of linear in-plane displacement, the results reported in [18] showed that the zero-shear strain/ stress condition has been not satisfied on the plate surfaces and it can lead to inaccurate results. It was also noted that the presence of shear correction factors (SCFs) [19] aims to adjust the unrealistic shear strain component of energy equation which occurs in plate and shell issues. However, such adjustment can lead to complicated computation and inaccuracy results due to the linear distribution of

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

shear components through the plate thickness. To overcome these drawbacks, numerous researches of higher-order shear deformable theory (HSDT), which focus on the justification of distribution function f ðzÞ, have been addressed [20,21]. It is found that better numerical results can be achieved as pointed out in many reports such as third-order shear deformation theory (TSDT) [21,22], fifthorder shear deformation theory (FiSDT) [23,24], trigonometric shear deformation theory [25–28], exponential shear deformation theory (ESDT) [29–31] and so on. More importantly, the requirement of HSDT theory is C1-continuous, which is not able to implement in the standard finite element method. There are some improvements for consecutive condition by using the C° continuous elements or classical Hermite interpolations which involve in transverse and rotations [32,33]. However, the increasing unknown variables are inevitable [32] and arise high computational costs. In this paper, a modified Moving Kriging interpolation function [34,35] is implemented into the meshfree method and the present method does not require any additional variables. In the framework of numerical methods, many meshless methods [36–54] have been shown as efficient computational tools for engineering problems. They were named as element-free Galerkin method (EFG) [36], the reproducing kernel particle method (RKPM) [37], the smoothed particle hydrodynamics (SPH) [38,39], the meshless local Petrov–Galerkin (MLPG) [40,41], etc. Among them, the element-free Galerkin method which employed the moving least squares (MLS) approximation [42,43] was very widely used. However, the major difficulty is that the moving least square function (MLS) does not satisfy the Kronecker delta condition [44]. Hence, the special treatment of the essential boundary condition needs to be provided such as penalty method [48], Lagrange multipliers [49], Transformation method [50]. More importantly, an improved meshfree method with radial point interpolation function [51–54] used the smoothed technique as found in Refs. [55–60] can resolve the adverse condition of Kronecker delta. This will be a great motivation for our forthcoming study. On the other hands, an alternative approach of a combination of Galerkin weak form and Moving Kriging (MK) interpolation function [61] allows us to overcome Kronecker's delta property. Such a method with the essential boundary conditions is easily enforced in a similar way as the conventional finite element method (FEM). However, the MK distribution interpolation based on the original Gaussian exponential function depends on the correlation parameter θ [62]. The evaluation of MK shape function is accurately determined by choosing an arbitrarily correlation factor on a bound of 0.1C500 and yet any publications to search an optimal value θ. This causes the instability in solution and even undetermined solution in several cases. In this work, we investigate a modified distribution of MK shape function and numerically show that the proposed method is insensitive to the correlation parameter θ. By a minor justification, the present formulation is significantly improved with high accuracy. In addition, the in-plane displacement is written in a compact form of higher order displacement field and constant transverse deflection through the plate thickness. Numerical study points out that the present method with the normalization of MK distribution function via quartic polynomial function can produce more stable solutions than the published approach [62]. The power and Mori–Tanaka models are utilized to homogenize the material properties. Numerical examples of isotropic and sandwich FGM plates are shown to validate high performance of the proposed method.

123

2. The higher order shear deformation theory for FGM plates 2.1. Problem formulation A rectangular plate is illustrated with dimension of a, b and plate thickness h as shown in Fig. 1. Three functionally graded material (FGM) plates typically consist of isotropic plates (type 1), sandwich plates using with FGM core and isotropic skins (type 2) and vice versa (type 3). 2.1.1. Isotropic FGM plates (type 1) The FGM plate is composed of two different material phases which are ceramic and metal arranging at the top and bottom of the plate, respectively. The effective moduli of two constituents are homogenized by the rule of mixture or the Mori–Tanaka models which are used to evaluate the effective elastic properties of the grade composite. The volume fraction of two distinct phases is described as [32]
V c ðzÞ ¼

   1 z n þ ; z A  h=2; h=2 ; 2 h

Vm ¼ 1Vc

ð1Þ

where m and c are defined to the metal and ceramic. The power index n is representative to the volume fraction varying through the thickness. The rule of mixture [32] used to estimate the effective moduli of material is defined as P e ¼ P c V c ðzÞ þ P m V m ðzÞ

ð2Þ

where P c ; P m represent the material constituents of the ceramic and the metal, respectively. Young’s modulus (E), Poisson’s ratio (ν), and density (ρ) can follow the rule defined by Eq. (2). However, the rule of mixture does not reflect the interactions among the two materials [63,64]. Meanwhile, the Mori–Tanaka model [65] is assumed to calculate their interactions through the effective bulk and shear modulus given by Ke  Km Vc ¼ ; Km K c  K m 1 þV m K Kþc 4=3 μ m

m

μe  μm Vc ¼ μc  μm 1 þ V m μμc þμfm1

ð3Þ

m

where f1 ¼

μm ð9K m þ 8μm Þ 6ðK m þ 2μm Þ

ð4Þ

The effective Young's modulus Ee and Poisson's ratio respectively now written as Ee ¼

9K e μe ; 3K e þ μe

νe ¼

3K e 2μe 2ð3K e þ μe Þ

νe are ð5Þ

The comparison of the effective Young's modulus corresponding with various power index values of Al/ZrO2 FGM plate between the mixture and Mori–Tanaka schemes is illustrated in Fig. 2. For the homogeneous constituents, the identical values are given by two models. However, when the material is inhomogeneous, the effective moduli of the Mori–Tanaka scheme is lower than that of the rule of mixture one.

Fig. 1. A typical FGM plate.

124

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136




0.5 n=10

0.4

Non−dimentional thickness z/h

n

  h ; h A z3 ; ; top skin 2

Vm ¼ 1Vc

0.3

ð7Þ

n=5

0.2

The thickness index of each layer of bottom-core-top (hb  hc  ht ) is indicated as different ratios “1-1-1”, “2-1-1” and so on. For example, “1-1-1” means that the thickness of 3 layers is identical.

n=3

0.1 0

n=1

−0.1

2.2. The higher order shear deformation theories

−0.2

n=0.5

−0.3

n=0.1

n=0.3

n=0

−0.4 −0.5

z4  z z4  z3

V c ðzÞ ¼

1

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

3

Effective modulus Eeff / Em

Fig. 2. The effective modulus of Al/ZrO2 FGM plate computed by the rule of mixture (in solid line) and the Mori–Tanaka (in dash dot line).

Let us consider a domain Ω in ℝ2 provided in the mid-plane of the plate. The displacements u0, v0, w are indicated in three x, y, z directions while the rotations βx and βy are identified in two planes of x–z and y–z (for the-y and x axes), respectively. General speaking, the higher-order shear deformation theories [28,31,66,67] are defined as uðx; y; zÞ ¼ u0 ðx; yÞ  z∂w ∂x þ f ðzÞβ x ðx; yÞ

vðx; y; zÞ ¼ v0 ðx; yÞ  z∂w ∂y þ f ðzÞβ y ðx; yÞ ;




wðx; y; zÞ ¼ wðx; yÞ

 h h rzr 2 2

ð8Þ

The displacement fields in Eq. (8) can be written in a compact form as u ¼ u0 þ zu1 þ f ðzÞu2 where 8 9 > = < u0 > u0 ¼ v 0 ; > ; : > w

Fig. 3. The sandwich plate with FGM core and isotropic skins.

metal ceramic

z4 z3

ceramic ceramic metal

ht hc

z2 z1

hb

ð9Þ

9 8 > = <  w;x > u1 ¼  w;y ; > > ; : 0

u2 ¼

8 9 > = < βx > > :

βy 0

ð10Þ

> ;

where w;x and w;y are the slopes in x, y directions, respectively. Following the classical theories (i.e, FSDT, CPT), the transverse distribution of the shear strain/stress is not satisfied the traction free on the top and bottom surfaces of plate. To improve these shortcomings, a shape function f ðzÞ through the thickness is taken into account in the higher order displacement field as given in Eq. (8). Several distribution functions are reported in Table 1. The following equation of strains and displacements relationship are described by n oT εp ¼ εxx εyy γ xy ¼ ε0 þ zε1 þ f ðzÞε2

Fig. 4. The sandwich plate with FGM skins and homogeneous core.

n

2.1.2. Sandwich plate with FGM core and isotropic skins (type 2) A sandwich plate (type 2) of FGM core layer and isotropic skins layers is drawn in Fig. 3. The volume fraction of the ceramic and metal phase in the core is defined [32] as
 1 zc n þ ; Vm ¼ 1Vc ð6Þ V c ðzÞ ¼ 2 hc where zc A ½z2 ; z3  and the thickness of core, hc ¼ z3  z2 . 2.1.3. Sandwich plates with isotropic core and FGM skins (type 3) The third case is a sandwich plate with two FGM skin layers and one isotropic core layer. In the bottom layer, the material components are changed from a metal surface (z ¼z1) to a ceramic surface (z¼ z2) and vice versa, in the top skin layer, a ceramic surface (z¼ z3) to a metal surface (z¼z4), as shown in Fig. 4. The ceramic and metal phases are expressed by a volume fraction [76] as
   z  z1 n h ; h A  ; z2 ; bottom skin V c ðzÞ ¼ 2 z2  z1 V c ðzÞ ¼ 1; h A ½z2 ; z3 ; core

γ ¼ γ xz γ yz

oT

0

¼ f ðzÞεs

ð11Þ

where 8 > > <

∂u0 ∂x ∂v0 ∂y 0¼ > > : ∂v0 þ ∂u0 ∂x ∂y

ε

9 8 9 8 2 9 ∂β x > > > > >  ∂∂xw2 > ∂x > > > > > ( ) > > > > > = = = < ∂β < βx 2 y ∂ w  ∂y2 ; εs ¼ ; ε2 ¼ ; ε1 ¼ ∂y β > > > > > y > > > > > ; > > ∂β βx > > ; > :  2 ∂2 w > ; : ∂xy þ ∂∂y ∂x∂y ð12Þ

0

in which f ðzÞ is the derivative function. Table 1 Several shape functions and its derivatives. 0

Model

f ðzÞ

Reddy [32] (TSDT 1984)

z  43z3 =h

Karama [29] (ESDT 2003)

ze  2ðz=hÞ

Nguyen [23] (FiSDT 2013)

7 2 3 2 5 8z  h2 z þ h4 z

Thai [26] (ITSDT 2014)

f ðzÞ 2

2

 h tan  1 2hz  z

2

1 4z2 =h  2 1  42 z2 e  2ðz=hÞ h 7 6 2 10 4 8  h2 z þ h4 z



2 

2 1  2hz = 1 þ 2hz

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

The static problem model can be expressed in a weak form under a distributed load as 8 9T 2 38 9 Z Z Z > < ε0 > = A B E > < ε0 > = 7 s δ ε1 6 δwq0 dΩ 4 B D F 5 ε1 dΩ þ δεs T D εs dΩ ¼ > Ω > Ω Ω :ε > ; :ε > ; E F H 2 2 ð13Þ Z Z Dsij ¼ The material 2 1 Ee 6 ν e Q¼ 4 1  νe 2 0

h=2  h=2 h=2

 h=2

2

ð1; z; z2 ; f ðzÞ; zf ðzÞ; f ðzÞÞQ ij dz;

0 2 f ðzÞ Gij dz

0



1

ð15Þ

Regarding free vibration problem, a weak form can be obtained from the dynamic equation by 8 9T 2 38 9 Z Z > < ϵ0 > = A B E > < ϵ0 > = 6 7 s δ ϵ1 4 B D F 5 ϵ1 dΩ þ δϵs T D ϵs dΩ > Ω > Ω :ϵ > ; :ϵ > ; E F H 2 2 8 9T 2 9 38 Z > < u0 > = I1 I2 I4 > < u0 > = 6 7 ð16Þ þ δ u1 4 I 2 I 3 I 5 5 u 1 dΩ ¼ 0 > > Ω > :u > ; I I5 I6 : u 2 ; 2 4 where Z

h=2  h=2

ϕI ðxÞ ¼

m X

pj ðxÞAjI þ

n X

j

r k ðxÞBkI





ρðzÞ 1; z; z2 ; f ðzÞ; zf ðzÞ; f 2 ðzÞ dz

The matrices A and B can be derived by the following equations A ¼ ðPT R  1 PÞ  1 PT R  1 ;

B ¼ R  1 ðI  PAÞ

ð17Þ

ð18Þ

3. A modified Moving Kriging interpolation function for FGM plate formulation 3.1. Introduction to MK shape functions The MK interpolation is used to build the shape functions and its derivatives through a set of nodal points. The estimation of the arbitrary distribution point uðxi Þ can be defined in the subdomain Ωx and Ωx A Ω. The interpolation domain is based on all nodal values, where xi ði A ½1; nÞ, n is a set of nodes in the support domain Ωx . The MK distribution interpolation can be applied to uh ðxÞ for every x in the support domain Ωx and is expressed as T

u ðxÞ ¼ ½p ðxÞA þ r ðxÞBu

and rðxÞ vector is defined by rðxÞ ¼ fRðx1 ; xÞ Rðx2 ; xÞ…Rðxn ; xÞgT

uh ðxÞ ¼

n X 1

ϕI ðxÞuI

ð27Þ

where r ij ¼ ‖xi  xj ‖, and θ 4 0 is the parameter which is strongly dependent coordinates in the support domain. However, the quality of MK depends heavily on the correlation parameter θ [62] which often causes the instability in the numerical model and its optimal value is still questionable. In this paper, we propose a polynomial function, which uses the quartic polynomial function, normalized as
2
3
4 θ 2 θ 3 θ 4 Rðxi ; xj Þ ¼ 1  6  r ij þ 8  r ij  3  r ij ð28Þ rs rs rs where r s is the scale factor used to normalize the distance. Here r s is taken to be the maximum distance between a pair of nodes in the support domain Ωx . The covariance matrix R A ℝnn is expressed as follows 2 3 1 Rðx1 ; x2 Þ ::: Rðx1 ; xn Þ 6 7 1 ⋯::: Rðx2 ; xn Þ 7 6 Rðx2 ; x1 Þ 7 ð29Þ R½Rðxi ; xj Þ ¼ 6 6 7 ⋮ ⋮ ⋱ ⋮ 4 5 1 Rðxn ; x1 Þ Rðxn ; x2 Þ ⋯ The successive derivatives of the modified MK shape functions can be readily illustrated as

ϕI;i ðxÞ ¼

m X

pj;i ðxÞAjI þ

n X

j

ϕI;ii ðxÞ ¼

r k;i ðxÞBkI

k

m X j

ð20Þ

ð26Þ

where Rðxi ; xÞ is the covariance basis function regarding an arbitrary pair of n nodes ðxi ; xÞ. There have many ways to choose this correlation function in which Gaussian function was widely used. The interpolation function based on the Gaussian form [62] is

ð19Þ

or Eq. (19) can be described in a compact form as

ð23Þ

For example, the basis functions for a two-dimensional problem can be expressed by a linear and quadratic form as follows, respectively,

pT ðxÞ ¼ 1 x y ; ðm ¼ 3Þ and n o ð25Þ pT ðxÞ ¼ 1 x y x2 xy y2 ; ðm ¼ 6Þ

Rðxi ; xj Þ ¼ e  θrij

where N0x , N0y and N 0xy are pre-buckling loads in x, y and x–y directions, respectively.

T

ð22Þ

In Eq. (22), P is the matrix of size n  m, which contains the values of all polynomial basis functions at n nodes in Ωx , 2 3 p1 ðx1 Þ p2 ðx1 Þ ::: pm ðx1 Þ 6 7 ⋮ ⋱ ⋮ 7 6 ⋮ 7 ð24Þ P¼6 6 p1 ðx2 Þ p2 ðx2 Þ … pm ðx2 Þ 7 4 5 p1 ðxn Þ p2 ðxn Þ ⋯ pm ðxn Þ

2

For buckling analysis, a weak form can be expressed as 8 9T 2 38 9 Z > Z < ϵ0 > = A B E > < ϵ0 > = 6 s δ ϵ1 4 B D F 7 5 ϵ1 dΩ þ δϵs T D ϵs dΩ > Ω > Ω :ϵ > ; :ϵ > ; E F H 2 2 2 0 3( ( ) ) T Z N x N 0xy w;x w;x 4 5 þh δ dΩ ¼ 0 w;y w;y N 0xy N 0y Ω

h

ð21Þ

k

pðxÞ ¼ fp1 ðxÞ p2 ðxÞ ::: pm ðxÞgT ð14Þ

matrices are given as 3 νe 0  1 Ee 7 1 0 5; G ¼ 2ð1 þ νe Þ 0 0 ð1  νe Þ=2

ðI 1 ; I 2 ; I 3 ; I 4 ; I 5 ; I 6 Þ ¼

where ϕI ðxÞ is the MK function and can be rewritten as

where vector p(x) can be defined with m basis polynomial terms while I is the unit matrix

where Aij ; Bij ; Dij ; Eij ; F ij ; H ij ¼

125

pj;ii ðxÞAjI þ

n X

r k;ii ðxÞBkI

ð30Þ

k

where ϕI;i ðxÞ and ϕI;ii ðxÞ are the first- and second-order spatial derivatives, respectively.

126

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

In meshfree methods, a set of nodes located in the support domain is used to construct shape functions. The rectangular or circle support domain can be applied, but the circle domain is often used and its radius is defined by dm ¼ αdc

ð31Þ

where dc is a characteristic length related to the nodal space, which is taken to be equal to the nodal spacing with regular mesh nodes. α stands for the scale factor. In practice, the radius of each influence domain is chosen so that it has a sufficient number of nodes to construct shape functions. 3.2. MK approximations applying for a higher-order plate formulation Meshfree with MK approximations can be written for both the description of the geometry (or the physical point) and the displacement field uh of the plate as n

n

xh ¼ ∑ ϕI ðxÞxI ; I

uh ðxÞ ¼ ∑ ϕI ðxÞqI



ð32Þ

I



where xT ¼ x y is the physical coordinate vector. In Eq. (32), ϕI ðxÞ is the modified MK basic function, xI is the h iT is the vector of nodal values and qI ¼ u0I v0I wI β xI βyI nodal degrees of freedom associated with the arbitrary nodal I. Substituting Eq. (32) into Eq. (12), the relationship between inplane/shear strain and displacement can be expressed as

ϵ0 ¼

n X

n X

Bm I qI ; ϵ1 ¼

I¼1

Bb1 I qI ; ϵ2 ¼

I¼1

3

ϕI;x 0 0 0 0 ϕI;y 0 0 0 7 7; 5 ϕI;y ϕI;x 0 0 0

6 60 Bm I ¼4

0

2

0

" BsI ¼

0 0

BsI qI

ð33Þ

I¼1

0 6 b1 0 6 BI ¼ 4 0

3

ϕI;x 0 6 b2 ϕI;y 7 0 0 0 0 6 7; BI ¼ 4 5 0 0 0 ϕI;y ϕI;x 0

n X

Bb2 I qI ; ϵs ¼

I¼1

in which 2

2

n X

 ϕI;xx

0

0

0

 ϕI;yy

0

0

 2ϕI;xy

0

0 0

0 0

0

#

ϕI 0 ; ϕI

ð34Þ

Similarly, substituting Eq. (32) into Eq. (9), the displacement fields u0 , u1 and u2 can be expressed as follows u0 ¼

n X

N0I qI ; u1 ¼

I¼1

where 2 ϕI 60 0 NI ¼ 4 0 2 0 6 N2I ¼ 4 0 0 (

n X

N1I qI and u2 ¼

I¼1

0

0

0

ϕI

0

0

0

ϕI

0

0

0

ϕI

0

0

0

0

0

0

n X

N2I qI

0

3

07 5

0 3 0 7 ϕI 5 0

ð35Þ

2

N1I

0 6 0 ¼4 0

0

−ϕI;x

0

0

−ϕI;y

0

0

0

0

0

7 0 5 and 0 ð36Þ

ð37Þ

where ω; λcr A ℝ are the natural frequency and the critical buckling value, respectively. The global stiffness matrix K can be derived by 28 3 9T 8 m 9 m B > >B > > > > > > > > > 6> 7 2 3 > A B E > > > > > > > 7 > > > Z 6 = = < < 6> 7 6 6 7 sT s s 7 b1 b1 þ B D B 7d Ω ð42Þ K¼ 6 B 4B D F 5 B > > > > 7 Ω6 > > > > > > 6> 7 E F H > > > > > > > > 4> 5 > > ; ; : b2 > : b2 > B B and the load vector f is computed by Z f¼ ΦT q0 dΩ

ð43Þ

where h Φ¼ 0

ð44Þ

Ω

0

ϕI 0 0

The global mass 8 0 9T 2 Z > = I1 6I M¼ 4 2 N1 Ω> : 2> ; I 4 N

i

matrix M is expressed as 38 0 9 I2 I4 > = I3 I5 7 5 N 1 dΩ > > I 5 I 6 : N2 ;

ð45Þ

The geometric stiffness matrix is given by 2 3 0 0 Z

g T N x Nxy 5B g d Ω B 4 0 Kg ¼ h Nxy N 0y Ω

ð46Þ

4. Results and discussions In order to certify the effectiveness of the present theory, several examples of FGM isotropic and sandwich plates are presented. For comparison, a quadratic polynomial basic functions (m ¼6) is used to construct MK interpolation functions. A background mesh with 4  4 Gauss points are used. Four defined types of distributed functions through the plate thickness: cubic [32], quintic [23] polynomials, exponent [29] and trigonometric [26] are illustrated. The material parameters are given in Table 2. Two types of boundary conditions are

ð47Þ

 Clamped (C) u0 ¼ v0 ¼ w ¼ β x ¼ β y ¼ w;n ¼ 0

ð48Þ

The Dirichlet boundary conditions of u0 ; v0 ; w; βx and β y can be easily and directly imposed as in the standard FEM. For the normal slope w;n , it can be directly enforced by assigning zero value of the transverse displacement at nodes that are adjacent to Table 2 Material properties.

# ð38Þ

Substituting Eq. (33), (35) and (37) into Eqs. (13), (16) and (18), the static, free vibration and buckling problems can be written by Kq ¼ f

ð41Þ þ

 Simply supported (S) u0 ¼ w ¼ βx ¼ 0 at y ¼ 0; b v0 ¼ w ¼ β y ¼ 0 at x ¼ 0; a

3

I¼1

ϕI;x 0 0 ϕI;y 0 0

 K  λcr Kg q ¼ f0g

I¼1

The derivations of the displacements are also given by ) n X w;x BgI qI ¼ w;y

where " 0 0 BgI ¼ 0 0

ð40Þ

3

7 0 7; 5 0

0

 K  ω2 M q ¼ f0g

ð39Þ

E (GPa) ν ρ (kg/m3)

Al

SiC

ZrO2-1

ZrO2-2

Al2O3

70 0.3 2707

427 0.17 –

200 0.3 5700

151 0.3 3000

380 0.3 3800

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

Uniform nodal structure;

127

Modified nodal structure for the clamped boundary condition.

Fig. 5. (a) Uniform nodal structure; (b) modified nodal structure for the clamped boundary condition.

the boundary points as same as isogeometric analysis [23]. This is referred to the rotation-free approach in the literature. Its application to FGM plates can be found in [45-47]. This means that in this research we need to define a special distribution of nodes such that the structure of this contribution is similar to that of control points in isogeometric analysis, as shown in Fig. 5. It is clear that using this scheme is more facilitated than several other techniques [48–50] which lead to increasing the unknown variables in the solution formulation. 4.1. Effect of the correlation parameter on the solution accuracy For a sake of illustration, a rectangular plate using a set of 23  23  nodes

 is subjected to a sinusoidal load q0 ¼ q0 sin πax sin πby . The boundary condition is simply supported (SSSS) and Al/SiC material properties are given in Table 2. We need to verify that the method is insensitive to the correlation parameter. In meshless methods, the parameter α of the influence domain can be fixed at 2.4 (of course this value should be chosen to be from 2.0 to 4.0 for the linear analysis [68]). More importantly, the correlation factor θ is an important one to decide the accuracy of solution as well as the quality of shape functions. The fact that the Gaussian function [62] embed in MK distributions is proven to be completely unstable under a small θ value and makes undetermined for all size meshes with a bound in 1C500. The solution becomes relatively stable as gradually increasing θ from 500 to 10,000. However, the accuracy of numerical solutions is not often expected. This paper therefore proposes a new approach based on a normalized MK function. The suggested quartic polynomial interpolation handles easily numerical computations. Table 3 shows that a non-dimensional central deflection using a quartic polynomial interpolation fits well the exact solution [69] for each value of θ from 1C10; 000. As expected, θ does not depend on background mesh sizes or nodal coordinates. As a result, the correlation parameter θ could be taken equal to 1. It is seen that the quartic function can yield the excellent convergence and ensure for the stable numerical computation. 4.2. Static analysis 4.2.1. Isotropic FGM plates Let us consider a FGM plate with a simply supported (SSSS)



 boundary condition and a sinusoidal load q0 ¼ q0 sin πax sin πby . The material property of Al/Al2O3 plate (type 1) is used. The non-dimensional transverse displacement and stresses are

Table 3 Normalized displacement w of a simply supported square plate. Distribution function

θ

1 Gaussian Quartic spline

3D Exact solution [69] 100

500

5000

10,000 a

#N/A #N/A #N/A 1.8655 (  0.6) 1.8768 1.8768 1.8768 1.8768 (0.005)

1.8656 1.8768

1.8767

a The percent of error compared to the exact solution is given within parentheses.

expressed as

 3 10h Ec a b h a b ; ; 0 ; ; ; z and w σ ¼ σ x x 2 2 aq 0 2 2 a4 q 0 a h τ yz ¼ τyz ; 0; z aq 0 2

w¼

ð49Þ

The normalized displacement and axial stress at the plate center are studied. The ratios a/h (side/thickness) of 4, 10 and 100 and various values n of 1, 4, 10, respectively, are used and listed in Table 4. The numerical results are good agreement with several solutions reported in the literatures such as the analytical solutions based on the generalized shear deformation theory (GSDT) by Zenkour et al. [67], Carrera's Unified Formulation (CUF) by Carrera et al. [70] and Brischetto [71], the quasi-3D meshfree solution with the sinusoidal and higher shear deformation theories by Neves et al. [72,73]. As seen in Table 4, although the present method does not take into account the stretching effect (i.e., εz ¼ 0), the displacement and axial stresses obtained are relatively reasonable in comparison with those of several other methods. Note that the small discrepancy of solutions as compared with quasi-3D values is decreased when the plate is thinner because the stretching effect reduces. The numerical results derived from the present method match well with other published ones. In addition, the stress distribution is continuous through the plate thickness and the traction-free on boundary surfaces of plate is ensured. As observed in Fig. 6 with a/ h¼4 and n¼ 1, it is found that the axial stress distribution of the present method well matches for several higher order deformation theories while shear stress distribution is slightly deviated, yet acceptable.

128

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

The second example is a rectangular Al/ZrO2-1 plate subjected to a uniformly distributed load. The Mori–Tanaka technique [65] is used to mix the ceramic and metal materials for this model (type 1). The non-dimensional central deflection is normalized by

Various boundary conditions such as simply supported (SSSS), clamped (CCCC) or supported–free (SFSF), respectively, are used to validate the present method. The results in Table 5 can be compared with the analytical results from the higher-order shear and normal deformable plate theory (HOSNDPT) by Gihooley et al. [74] and the isogeometric analysis based on TSDT by Tran et al. [75]. It can be explained that the stiffness of plates is increased when the boundary forms are from SFSF to CCCC. The displacement amplitudes are also varied and gradually reduced accordingly. Additionally, the higher exponent value n leads to the effect of decreasing the volume fraction V c of Zirconia through the plate thickness, and therefore the magnitude of the deflection is gradually increased. A comparison between the present solution and other ones for different values n is shown in Fig. 7. It is clear that the present method is a strong competitor to other Refs. [74,75].

3

wc ¼

100Em h wc 12ð1  ν2 Þq0 a4

Table 4 The non-dimensional displacement and axial stress for a simply supported (SSSS) Al/Al2O3 square plate under a sinusoidal load. a/h ¼4

εz

Model

a/h ¼ 10

a/h ¼ 100

w

 σ x z ¼ 3h w



 h

σx z ¼ 3

w

 σ x z ¼ h3

0 0 a0 a0 a0 0 a0 0 0 0 0

– 0.7289 0.7171 0.7171 0.6997 0.7308 0.702 0.7299 0.7264 0.7285 0.7268

– 0.7856 0.6221 0.6221 0.5925 0.5806 0.5911 0.5875 0.5844 0.5859 0.5858

0.5889 0.589 0.5875 0.5875 0.5845 0.5913 0.5868 0.5897 0.5893 0.5896 0.5893

1.4894 2.0068 1.5064 1.5064 1.4945 1.4874 1.4917 1.5044 1.5029 1.5036 1.5035

– 0.5625 0.5625 0.5625 0.5624 0.5648 0.5647 0.5624 0.5624 0.5624 0.5624

– 20.149 14.969 14.969 14.969 14.944 14.945 15.0925 15.0923 15.0924 15.0924

GSDT [67] CUF [70] CUF [70] CUF [71] SSDT [72] HSDT [73] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0 0 a0 a0 a0 0 a0 0 0 0 0

– 1.1673 1.1585 1.1585 1.1178 1.1553 1.1108 1.1621 1.1635 1.1648 1.1642

– 0.5986 0.4877 0.4877 0.4404 0.4338 0.433 0.4503 0.4421 0.4456 0.4442

0.8651 0.8828 0.8821 0.8821 0.875 0.877 0.87 0.8827 0.8831 0.8832 0.8832

1.1783 1.5874 1.1971 1.1971 1.1783 1.1592 1.1588 1.1931 1.1896 1.1911 1.1905

– 0.8286 0.8286 0.8286 0.8286 0.8241 0.824 0.8284 0.8284 0.8284 0.8284

– 16.047 11.923 11.923 11.932 11.737 11.737 12.0369 12.0366 12.0367 12.0366

10 GSDT [67] CUF [70] CUF [70] CUF [71] SSDT [72] HSDT [73] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0 0 a0 a0 a0 0 a0 0 0 0 0

– 1.3925 1.3745 1.3745 1.349 1.376 1.3334 1.3939 1.3891 1.3936 1.3900

– 0.4345 0.1478 0.3695 0.3227 0.3112 0.3097 0.3295 0.3217 0.3250 0.3236

1.0089 1.009 1.0072 1.0072 0.875 0.9952 0.9888 1.0104 1.0099 1.0105 1.0100

0.8775 1.1807 0.8965 0.8965 1.1783 0.8468 0.8462 0.8878 0.8844 0.8858 0.8852

– 0.9361 0.9361 0.9361 0.8286 0.9228 0.9227 0.9361 0.9361 0.9361 0.9361

– 11.989 8.9077 8.9077 11.932 8.6011 8.601 8.9826 8.9823 8.9824 8.9824

4

Thickness z/h

GSDT [67] CUF [70] CUF [70] CUF [71] SSDT [72] HSDT [73] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

1

Table 5 The non-dimensional displacement of an Al/ZrO2-1 plate (a/h ¼ 5) under a uniformly distributed traction and different boundary conditions. BCs

0.5063 0.7572 0.8735 0.9780 1.0640 1.1502 1.4465

0.5069 0.7582 0.8745 0.9792 1.0654 1.1520 1.4484

0.5064 0.7574 0.8737 0.9783 1.0643 1.1506 1.4469

SSSS

Ceramic 0.5 1 2 4 8 Metal

0.1671 0.2505 0.2905 0.328 – – 0.4775

0.1716 0.2554 0.2955 0.3334 0.3655 0.3958 0.4903

0.1718 0.2557 0.2959 0.3338 0.3660 0.3963 0.4909

0.1713 0.2550 0.2950 0.3331 0.3652 0.3963 0.4893

0.1716 0.2554 0.2956 0.3336 0.3659 0.3958 0.4903

0.1713 0.2551 0.2952 0.3332 0.3654 0.3951 0.4895

CCCC Ceramic 0.5 1 2 Metal

0.0731 0.1073 0.1253 0.1444 0.2088

0.0734 0.1077 0.1256 0.1447 0.2098

0.0731 0.1072 0.1250 0.1439 0.2089

0.0723 0.1060 0.1237 0.1426 0.2065

0.0727 0.1067 0.1245 0.1434 0.2078

0.0724 0.1062 0.1239 0.1428 0.2068

0.25

0.3

0.35

0.3

0.2

0.2 0.1

−0.1 −0.2

−0.3

−0.3

−0.4

−0.4 −0.2

0

0.2

Axialstress

0.4

0.6

0.8

σ¯x (a/ 2,b/ 2,z )

1

1.2

MK−TSDT MK−ESDT MK−FiSDT MK−ITSDT

0

−0.2

−0.4

ITSDT

0.5074 0.7588 0.8752 0.9797 1.0659 1.1532 1.4497

0.3

−0.5 −0.6

ESDT

0.5088 0.7607 0.8776 0.9830 1.0701 1.1577 1.4537

0.4

−0.1

FiSDT

0.5019 0.7543 0.8708 0.9744 – – 1.4345

0.5

MK−TSDT MK−ESDT MK−FiSDT MK−ITSDT

TSDT Ceramic 0.5 1 2 4 8 Metal

0.4

0

Model

SFSF

0.5

0.1

n

HOSNDPT [74] TSDT [75] MK

Thickness z/h

n

−0.5

0

0.05

0.1

0.15 Shearstress

0.2

τ¯yz (a/ 2, 0,z )

Fig. 6. The distribution axial σ x and shear stress τyz through the plate thickness.

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

The deformation shapes of the plate with different boundary conditions are depicted in Fig. 8.

4.2.2. Sandwich plates with FGM core Let us consider a rectangular sandwich plate (type 2). The Al/ Al2O3 material properties are described for top and bottom skins of

the plate in which Young's moduli for aluminum, alumina are Em ¼70 GPa and Ec ¼ 380 GPa, respectively. The thickness of each skin layer is equal to 0.1h, where h is the plate thickness. The plate



 is subjected to a sinusoidal load q0 ¼ q0 sin πax sin πby . The normalized displacement and the shear stress are given as

 3 10h Ec a b h b h ð50Þ w¼ 4 w ; ; 0 and τxz ¼ τxz 0; ; 2 2 aq0 2 6 a q0 Considering the different (a/h¼ 4, 10) ratios for each relevant exponential value n, the displacement and shear stress at the certain position in the plate thickness are studied. The obtained results from the present solutions are compared with the analytical solution based on CUF by Bischetto [71] and Carrera et al. [70], the meshfree solution based on SSDT and HSDT by Neves et al. [72,73]. It is indicated from Table 6 that the transverse displacement and shear stress are reasonable with the reference solutions. The error corresponding to the displacement evaluation is not decreased when increasing the thickness ratio. Therefore, the improvement of them by the use of layerwise theory will be more viable.

SFSF

SSSS

129

CCCC

4.3. Free vibration analysis

Fig. 7. The normalized vertical displacement of an Al/ZrO2-1 FGM plate with different values n and boundary conditions.

4.3.1. Isotropic FGM plates Let us consider an Al/ZrO2-1 square plate with SSSS boundary condition. Material characteristics are homogenized by the Mori–Tanaka scheme. The good comparison of the first

Table 6 The non-dimensional vertical displacement and transverse shear stress of a SSSS FGM rectangular sandwich plate (type 2) under a sinusoidal load. n

Model

εz

a/h ¼ 4

a/h ¼10

w

 τ xz z ¼ 6h

w

 τ xz z ¼ 6h

1

CUF [71] CUF [70] CUF [70] SSDT [72] SSDT [72] HSDT [73] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

a0 0 a0 0 a0 0 a0 0 0 0 0

0.7628 0.7735 0.7628 0.7744 0.7416 0.7746 0.7417 0.7740 0.7703 0.7725 0.7707

0.2613 0.2596 0.2604 0.2703 0.2742 0.2706 0.2745 0.2916 0.3003 0.2978 0.2985

0.6324 0.6337 0.6324 0.6356 0.6305 0.6357 0.6305 0.6344 0.6338 0.6342 0.6339

0.2605 0.2593 0.2594 0.2718 0.2788 0.2720 0.2789 0.2710 0.2835 0.2795 0.2818

4

CUF [71] CUF [70] CUF [70] SSDT [72] SSDT [72] HSDT [73] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

a0 0 a0 0 a0 0 a0 0 0 0 0

1.0934 1.0977 1.093 1.0847 1.0391 1.0826 1.0371 1.0837 1.0909 1.0899 1.0912

0.2429 0.24 0.24 0.2699 0.2723 0.2671 0.2696 0.2590 0.2812 0.2734 0.2790

0.8321 0.8308 0.8307 0.8276 0.8202 0.8272 0.8199 0.8289 0.8302 0.8300 0.8303

0.2431 0.2398 0.2398 0.2726 0.2778 0.2695 0.2747 0.2474 0.2740 0.2645 0.2717

10

CUF [71] CUF [70] CUF [70] SSDT [72] SSDT [72] HSDT [73] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

a0 0 a0 0 a0 0 a0 0 0 0 0

1.2232 1.224 1.2172 1.2212 1.178 1.2183 1.1752 1.2180 1.2303 1.2279 1.2292

0.215 0.1935 0.1932 0.1998 0.2016 0.1996 0.1995 0.2019 0.2207 0.2141 0.2183

0.8753 0.8743 0.874 0.8718 0.865 0.8712 0.8645 0.8740 0.8763 0.8758 0.8761

0.2174 0.1944 0.1944 0.2021 0.2059 0.2018 0.2034 0.1955 0.2180 0.2098 0.2152

130

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

given in Table 8 corresponding with a/h ¼ 5, 10, 20. Fig. 9 plots the first six mode shapes of the Al/ZrO2-1 simply supported square plate. The normalized frequencies obtained exhibit the high accuracy comparing to the exact solution when the plate thickness reduces. The present results agree well with those indicated in the literature.

normalized frequency between the modified MK model with 3D analytical solution [63], the quasi-3D meshfree solution of SSDT [72] and HSDT [73] as well as the analytical of HOSNDPT [64] is confirmed in Table 7. It is also found that the FiSDT using in meshfree method gives the best results as compared to the 3D exact solution. In addition, the first 10 normalized frequencies are

0.5 0.5

0.4

0

0

0.2 0

−0.5

−0.5 −1

−0.2

−1

−0.4 0

−1.5 −2 1

0.8

0.6

0.4

0.2

0.4

0.2

0 0

0.6

0.8

1

−1.5 0 0.2

0.2 0.4

0.4

0.6

0.8

1

0

0.2

0.4

0.8

0.6

0.6

1

0.8 1 0

0.2

0.4

0.6

0.8

1

Fig. 8. Deformation shapes of Al/ZrO2-1 FGM plates for different boundary conditions: (a) SFSF; (b) SSSS; (c) CCCC. Table 7 The first normalized frequency ω ¼ ωh

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ρm =Em of a SSSS Al/ZrO2-1 square plate with a/h ¼5.

εz

Model

3D [63] HOSNDPT [64] SSDT [72] SSDT [72] HSDT [73] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

n

a0 0 a0 0 a0 0 0 0 0

0

0.5

1

2

3

5

10

– – – – 0.2459 0.2469 0.2460 0.2464 0.2461 0.2463

– – – – 0.2219 0.2228 0.2221 0.2225 0.2223 0.2224

0.2192 0.2152 0.2184 0.2193 0.2184 0.2193 0.2184 0.2188 0.2185 0.2187

0.2197 0.2153 0.2189 0.2198 0.2191 0.22 0.2189 0.2192 0.2190 0.2191

0.2211 0.2172 0.2202 0.2212 0.2206 0.2215 0.2203 0.2205 0.2203 0.2204

0.2225 0.2194 0.2215 0.2225 0.222 0.223 0.2215 0.2218 0.2216 0.2218

– – – – 0.2219 0.2229 0.2210 0.2215 0.2212 0.2214

Table 8 The normalized frequencies ω of a SSSS Al/ZrO2-1 square plate with various ratios a/h. a/h

Model

Modes 1

(2,3)

4

5

6

7

8

9

10

3D [63] HOSNDPT [64] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0.2192 0.2152 0.2184 0.2188 0.2185 0.2187

– 0.4114 0.4107 0.4107 0.4107 0.4107

– 0.4761 0.4788 0.4803 0.4794 0.4801

– 0.4761 0.4788 0.4803 0.4794 0.4801

– 0.582 0.5825 0.5825 0.5825 0.5825

– 0.6914 0.6959 0.6993 0.6973 0.6988

– 0.8192 0.8162 0.8193 0.8182 0.8193

– 0.8217 0.8164 0.8194 0.8184 0.8194

– 0.8242 0.8193 0.8210 0.8193 0.8204

10

3D [63] HOSNDPT [64] SSDT [72] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0.0596 0.0584 0.0596 0.0596 0.0595 0.0596 0.0596 0.0596

– 0.1410 0.1426 0.1426 0.1421 0.1422 0.1421 0.1422

– 0.2058 0.2058 0.2059 0.2054 0.2054 0.2054 0.2054

– 0.2058 0.2058 0.2059 0.2054 0.2054 0.2054 0.2054

– 0.2164 0.2193 0.2193 0.2189 0.2192 0.2190 0.2192

– 0.2646 0.2676 0.2676 0.2639 0.2644 0.2641 0.2643

– 0.2677 0.2676 0.2676 0.2640 0.2645 0.2642 0.2644

– 0.2913 0.291 0.2912 0.2914 0.2914 0.2914 0.2914

– 0.3264 0.3363 0.3364 0.3342 0.3349 0.3345 0.3348

20

3D [63] HOSNDPT [64] SSDT [72] HSDT [73] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0.0153 0.0149 0.0153 0.0153 0.0153 0.0153 0.0153 0.0153

– 0.0377 0.0377 0.0377 0.0376 0.0377 0.0376 0.0376

– 0.0593 0.0596 0.0596 0.0597 0.0597 0.0597 0.0597

– 0.0747 0.0739 0.0739 0.0731 0.0732 0.0731 0.0731

– 0.0747 0.0739 0.0739 0.0731 0.0732 0.0731 0.0732

– 0.0769 0.095 0.095 0.0948 0.0949 0.0948 0.0949

– 0.0912 0.095 0.095 0.0948 0.0949 0.0948 0.0949

– 0.0913 0.1029 0.103 0.1027 0.1027 0.1027 0.1027

– 0.1029 0.1029 0.103 0.1027 0.1027 0.1027 0.1027

5

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

0.4

0.4 0 −0.5

0.2

0.2

0

0

−0.2 −1 1

−0.4 0

0.8 0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1

−0.2

0.2

0.4

0.6

0.8

1 1

0.8

0.6

0.4

0.2

−0.4 1

0

0.3

0.3

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

0.15 0.1

0.2

0.2

131

0.05

0.1

0.1

0

0 −0.1

0

−0.2

−0.1

−0.05 −0.1 −0.15

−0.3 −0.2 1

−0.4 1 0.5 0 0

0.2

0.4

0.6

0.8

1

0.5 0 0

1 0.6 0.8 0.2 0.4

−0.2 1 0.5 0 0

0.2

0.4

0.6

0.8

1

Fig. 9. The first six mode shapes for a SSSS square plate (n¼ 1).

Table 9 The flexural vibration frequencies ωi ði ¼ 1C5Þ of a SSSS sandwich FGM plate 2-1-2 (a/h ¼ 10).

Table 10 The natural frequency ω of a SSSS sandwich plate with other theories (a/h ¼ 10). n

n

Model

ω1

ω2

ω3

ω4

ω5

1

3D elasticity [76] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

1.3018 1.3002 1.3005 1.3008 1.3008

3.1588 3.1447 3.1463 3.1477 3.1479

3.1588 3.1447 3.1463 3.1477 3.1479

4.9166 4.9021 4.9058 4.9090 4.9093

6.0405 5.9469 5.9574 5.9523 5.9570

10

3D elasticity [76] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0.9404 0.9432 0.9438 0.9442 0.9441

2.2862 2.2973 2.3006 2.3028 2.3028

2.2862 2.2973 2.3006 2.3028 2.3028

3.5647 3.6034 3.6113 3.6164 3.6162

4.3844 4.3867 4.3982 4.4057 4.4054

4.3.2. Sandwich plate with FGM skins and isotropic core The (type 3) rectangular sandwich plate with material components in which ceramic material is of core while the metal for the top and bottom skins are used. The thickness of the sandwich plate is described with the ratios of bottom, core and top skins, e.g., hb hc  ht ¼ 2-1-2. The material properties of alumina and aluminum (Al/Al2O3) are given as Ec ¼380E0 , ρc ¼3800ρ0 , Em ¼ 70E0 and ρm ¼2707ρ0 , where ρ0 ¼1 kg/m3 and Ec ¼1 GPa. The non-dimensional frequency of the sandwich plate is given by

ω¼

rffiffiffiffiffiffi

ωa2 ρ0 h

Method

1-0-1

2-1-2

2-1-1

1-1-1

2-2-1

1-2-2

0.5

TSDT [79] SSDT [79] 3D elasticity [76] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

1.4442 1.4443 1.4461 1.4443 1.4450 1.4447 1.4450

1.4841 1.4842 1.4861 1.4842 1.4848 1.4844 1.4847

1.5125 1.5126 1.5084 1.5065 1.5070 1.5067 1.5069

1.5192 1.5193 1.5213 1.5193 1.5198 1.5195 1.5197

1.5520 1.5520 1.5493 1.5472 1.5476 1.5473 1.5476

1.5745 1.5745 1.5766 1.5745 1.5748 1.5746 1.5748

1

TSDT [79] SSDT [79] 3D elasticity [76] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

1.2432 1.2433 1.2447 1.2433 1.2440 1.2437 1.2440

1.3001 1.3002 1.3018 1.3002 1.3008 1.3005 1.3008

1.3489 1.3489 1.3351 1.3335 1.3339 1.3337 1.3339

1.3533 1.3534 1.3552 1.3534 1.3538 1.3536 1.3538

1.4079 1.4079 1.3976 1.3958 1.3961 1.3959 1.3961

1.4393 1.4393 1.4413 1.4394 1.4396 1.4394 1.4395

5

TSDT [79] SSDT [79] 3D elasticity [76] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0.9460 0.9463 0.9448 0.9461 0.9473 0.9468 0.9473

0.9818 0.9820 0.9810 0.9820 0.9829 0.9825 0.9829

1.0743 1.0744 1.0294 1.0308 1.0314 1.0311 1.0314

1.0447 1.0448 1.0453 1.0449 1.0454 1.0452 1.0454

1.1473 1.1474 1.1098 1.1092 1.1096 1.1094 1.1096

1.1740 1.1740 1.1757 1.1741 1.1743 1.1742 1.1743

10

TSDT [79] SSDT [79] 3D elasticity [76] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

0.9284 0.9288 0.9273 0.9285 0.9298 0.9293 0.9299

0.9430 0.9433 0.9418 0.9432 0.9441 0.9438 0.9442

1.0386 1.0455 0.9893 0.9923 0.9930 0.9927 0.9929

0.9955 0.9952 0.9952 0.9957 0.9964 0.9961 0.9964

1.1053 1.0415 1.0610 1.0613 1.0618 1.0616 1.0618

1.1231 1.1346 1.1247 1.1233 1.1236 1.1235 1.1235

Frequency

E0

Regarding different index values (n ¼1, 10), the first five frequencies which are tabulated in Table 9 are computed from the modified MK function and HSDT theories under the SSSS boundary condition. The present results are well compared with the

elasticity solution [76]. The decreasing natural frequencies corresponding to higher index value n are addressed. Table 10 presents the first non-dimensional frequency with six differential thickness ratios for four types of volume fraction indices

132

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1

Geometry and

−0.5

0

0.5

1

distribution of nodes of a circular plate.

Fig. 10. (a) Geometry and (b) distribution of nodes of a circular plate.

Table 11 The buckling load parameter of a clamped thick circular Al/ZrO2-2 plate. n

Method

h/R 0.1

0.2

0.25

0.3

0

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

14.089 14.089 14.086 14.107 14.095 14.104

12.574 12.575 12.522 12.571 12.543 12.566

11.638 11.639 11.594 11.661 11.623 11.653

10.67 10.67 10.643 10.726 10.679 10.717

0.5

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

19.411 19.413 19.419 19.446 19.431 19.443

17.311 17.31 17.250 17.317 17.278 17.308

16.013 16.012 15.965 16.054 16.003 16.043

14.672 14.672 14.647 14.760 14.696 14.737

2

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

23.074 23.075 23.098 23.132 23.113 23.128

20.803 20.805 20.735 20.817 20.772 20.809

19.377 19.378 19.318 19.428 19.367 19.416

17.882 17.881 17.845 17.984 17.909 17.921

10

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

27.133 27.131 27.141 27.176 27.156 27.172

24.423 24.422 24.324 24.407 24.359 24.398

22.725 22.725 22.637 22.749 22.685 22.737

20.948 20.949 20.889 21.030 20.950 21.025

4.4. Buckling analysis

are homogenized by the mixture formulation as Eq. (51) by
 1 z n P e ¼ P c V c ðzÞ þ P m V m ðzÞ where  V m ðzÞ ¼ ; V c ¼ 1  V m; 2 h ð51Þ

4.4.1. Isotropic FGM plates A clamped Al/ZrO2-2 circular plate is enforced by a radial compression load p0 as shown in Fig. 10a. A mesh of nodes distribution is plotted in Fig. 10b. The material properties of Al/ZrO2-2

Table 11 presents the numerical results of critical loading value 3 by expression of pcr ¼ pcr R2 =Dm , where Dm ¼ Em h =12ð1  ν2m Þ. The problem is studied with exponential values n and relevant different ratios of h/R. It is seen that the obtained results agree well with those of TSDT [77] and (UTSDT) [78]. The effects of the increasing

n. From these results, it is seen that the solutions obtained match well with those of TSDT [79], SSDT [79] and 3D elasticity [76].

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

133

1 0 0.5

−0.2 −0.4

0

−0.6 −0.5

−0.8 −1 1 0.5 0 −0.5

−0.5

−1 −1

0.5

0

−1 1

1

0.5

1 0.5

0

0

−0.5 −1

1

−0.5 −1

0.4 0.2

0.5

0 0 −0.2 −0.5

−0.4 1

−1 1

1

0.5 0.5

1 0.5

0

0.5

0

0 −0.5

0

−0.5 −1

−0.5

−0.5 −1

−1

−1

0.5

0.6 0.4 0.2

0 0 −0.2 −0.4 1 0.5

1 0.5

0 −0.5 −1

−0.5

−0.5 1 0.5

1 0.5

0

0

0

−0.5

−1

−1

−0.5 −1

Fig. 11. The buckling modes of circular plate with h/R¼ 0.1 and n ¼2.

of n and the smallest thickness of the plate can lead to the highest critical buckling loading. The exponential value n gives less dependent on the behavior of the circular plate when its value is higher than 10. As shown in Fig. 11, the first six buckling modes are presented for the circular plate with h/R¼ 0.1 and n ¼2. Another study for a FGM circular plate is performed by the static analysis of the vertical displacements for the proposed model and the others [77,78]. Table 12 confirms again the good performance of the present method. 4.4.2. Sandwich plate with FGM skins and isotropic core A final example is a SSSS rectangular FGM sandwich plate with a/h¼ 10. The material components of the sandwich plate are metal for top and bottom skins while ceramic is a core layer with Young’s

modulus Em ¼70 GPa and Ec ¼380 GPa, respectively. The buckling loading parameter under the uniaxial and biaxial is defined by 2 P cr ¼ P cr a3 , where E0 ¼1 GPa. 100h E0

The problem is calculated with various thickness ratios of the FGM plate in which six types hb  hc ht are illustrated in Tables 13 and 14. It is observed that the present solutions of the uniaxial and biaxial buckling loading parameters are compared to those of TSDT [32], SSDT [79] and HSDT [72,73]. A good agreement is confirmed.

5. Conclusions A generalized higher-order shear deformation theory in combination with a modified MK meshfree method has been

134

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

Table 12 3

Ec h 64 of a clamped thick circular Al/ZrO2-2 plate. The non-dimensional deflection w ¼ w12ð1ν 2Þ  q R4 c

n

0

Method

h/R 0.1

0.2

0.25

0.3

0

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

2.2546 2.2553 2.2785 2.2748 2.2766 2.2744

2.5423 2.5478 2.5585 2.5475 2.5528 2.5462

2.7554 2.7660 2.7562 2.7540 2.7549 2.7510

3.0134 3.0317 2.9922 2.9967 2.9791 3.0121

0.5

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

1.6364 1.6369 1.6523 1.6497 1.6510 1.6494

1.8470 1.8509 1.8565 1.8487 1.8525 1.8477

2.0027 2.0105 2.0009 2.0175 2.0123 2.0132

2.1913 2.2046 2.1731 2.1710 2.1638 2.1904

2

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

1.3763 1.3766 1.3891 1.3868 1.3879 1.3866

1.5355 1.5386 1.5451 1.5385 1.5417 1.5378

1.6535 1.6595 1.6549 1.6590 1.6544 1.6573

1.7964 1.8068 1.7859 1.7899 1.7779 1.7952

5

TSDT [77] UTSDT [78] MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

1.2483 1.2486 1.2605 1.2585 1.2595 1.2583

1.3908 1.3935 1.4005 1.3947 1.3974 1.3940

1.4963 1.5017 1.4989 1.4905 1.5003 1.5015

1.6242 1.6334 1.6162 1.6112 1.6093 1.6133

Table 13 The non-dimensional critical buckling parameters (uniaxial) of a SSSS sandwich plate with FGM skins and isotropic core. n

Theory

P cr 1-0-1

2-1-2

2-1-1

1-1-1

2-2-1

1-2-1

0

TSDT [32] SSDT [79] HSDT [72] εz a 0 HSDT [73] εz ¼ 0 MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

13.0050 13.0061 12.9529 13.0051 12.9859 12.9986 12.9906 12.9969

13.0050 13.0061 12.9529 13.0051 12.9859 12.9986 12.9906 12.9969

13.0050 13.0061 12.9529 13.0051 12.9859 12.9986 12.9906 12.9969

13.0050 13.0061 12.9529 13.0051 12.9859 12.9986 12.9906 12.9969

13.0050 13.0061 12.9529 13.0051 12.9859 12.9986 12.9906 12.9969

13.0050 13.0061 12.9529 13.0051 12.9859 12.9986 12.9906 12.9969

1

TSDT [32] SSDT [79] HSDT [72] εz a 0 HSDT [73] εz ¼ 0 MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

5.1671 5.1685 5.0614 5.0785 5.1612 5.1670 5.1643 5.1669

5.8401 5.8412 5.7114 5.7302 5.8333 5.8387 5.8361 5.8384

6.1939 6.1946 6.0547 6.0736 6.1867 6.1909 6.1886 6.1906

6.4647 6.4654 6.3150 6.3356 6.4571 6.4612 6.4590 6.4608

6.9494 6.9498 6.7841 6.8055 6.9411 6.9446 6.9425 6.9443

7.5066 7.5063 7.3200 7.3437 7.4973 7.4989 7.4974 7.4984

5

TSDT [32] SSDT [79] HSDT [72] εz a 0 HSDT [73] εz ¼ 0 MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

2.6582 2.6601 2.6365 2.6468 2.6554 2.6620 2.6594 2.6621

3.0426 3.0441 3.0079 3.0187 3.0397 3.0450 3.0429 3.0450

3.4035 3.4045 3.3626 3.3720 3.4005 3.4043 3.4026 3.4041

3.5796 3.5806 3.5301 3.5415 3.5762 3.5802 3.5785 3.5800

4.1121 4.1129 4.0507 4.0616 4.1083 4.1116 4.1116 4.1115

4.7347 4.7349 4.6470 4.6606 4.7298 4.7312 4.7303 4.7307

10

TSDT [32] SSDT [79] HSDT [72] εz a 0 HSDT [73] εz ¼ 0 MK-TSDT MK-FiSDT MK-ESDT MK-ITSDT

2.4873 2.4893 2.4722 2.4822 2.4843 2.4918 2.4888 2.4921

2.7463 2.7484 2.7205 2.7308 2.7438 2.7496 2.7474 2.7497

3.0919 3.1344 3.0607 3.0694 3.0892 3.0932 3.0915 3.0930

3.1947 3.1946 3.1576 3.1684 3.1918 3.1963 3.1945 3.1961

3.7075 3.1457 3.6617 3.6715 3.7042 3.7079 3.7063 3.7078

4.2799 4.3818 4.2055 4.2179 4.2756 4.2777 4.2766 4.2772

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

proposed for analysis of isotropic and sandwich FGM plates. Various distributed functions across the plate thickness as cubic, quantic polynomials, exponent and trigonometric were investigated. The shear correction factors are negligible and the zero tangential traction boundary condition is satisfied at the top and bottom surfaces of the FGM plates. The essential boundary conditions are enforced directly and simply through a rotation-free technique in a same way as isogeometric analysis without any additional variables. For all examples, the obtained results are in very good agreement compared with the other theories. Especially, the results match well with the 3D solutions and the closed solution forms of SSDT, HSDT were available in the literatures. The suggested method can be well applied to FGM plates which the high reliability is found for both thin and thick cases. It is also observed that the proposed method can produce the high accuracy of solutions for FGM sandwich plates. The important point is that the normalized quartic shape function is integrated into the moving-Kriging distribution interpolation. As a result, the obtained solutions are stable and more accurate while using Gaussian exponential correlation function leads to a large dependence on the correlation parameter θ. Finally, we believe that the proposed method opens more promising for practical applications.

Acknowledgements This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 107.02-2014.24.

References [1] Yamanouchi M, Koizumi M, Hirai T, Shiota I. Proceedings of first international symposium on functionally gradient materials. Sendai, Japan; 1990. [2] Koizumi M. Functionally gradient materials the concept of FGM. Ceram Trans 1993;34:3–10. [3] Hasselman DPH, Youngblood GE. Enhanced thermal stress resistance of structural ceramics with thermal conductivity gradient. J Am Ceram Soc 1978;61:49–53. [4] Liu GR, Han X, Lam KY. Stress waves in functionally gradient materials and its use for material characterization. Compos Part B 1999;30:383–94. [5] Han X, Liu GR, Lam KY, Ohyoshi T. A quadratic layer element for analyzing stress waves in FGMs and its application in material characterization. J Sound Vib 2000;236:307–21. [6] Liu GR, Tani J, Ohyoshi T. Lamb waves in a functionally gradient material plates and its transient response. Theory Trans Jpn Soc Mech Eng 1991;57:603–11. [7] Tani J, Liu GR. SH surface-waves in functionally gradient piezoelectric plates. JSME Int J Ser A Mech Mater Eng 1993;36:152–5. [8] Liu GR, Tani J. Surface-waves in functionally gradient piezoelectric plates. J Vib Acoust Trans ASME 1994;116:440–8. [9] Liu GR, Dai KY, Han X, Ohyoshi T. Dispersion of waves and characteristic wave surfaces in functionally graded piezoelectric plates. J Sound Vib 2003;268:131–47. [10] Han X, Liu GR, Lam KY. Transient waves in plates of functionally graded materials. Int J Numer Methods Eng 2001;52:851–65. [11] Han X, Liu GR. Computational inverse technique for material characterization of functionally graded materials. AIAA J 2003;41:288–95. [12] Han X, Liu GR, Xi ZC, Lam KY. Characteristics of waves in a functionally graded cylinder. Int J Numer Methods Eng 2002;53:653–76. [13] Dawe DJ, LAM SSE, Azizian ZG. Nonlinear finite strip analysis of rectangular laminates under end shortening, using classical plate theory. Int J Numer Methods Eng 1992;35:1087–110. [14] York CB, Williams FW. Buckling analysis of skew plate assemblies: classical plate theory results incorporating largrangian multiplier. Comput Struct 1995;56:625–35. [15] Rolfes R, Rohwer K. Improved transverse shear stresses in composite finite elements based on first order shear deformation theory. Int J Numer Methods Eng 1997;40:51–60. [16] Jonnalagadda KD, Blandfords GE, Tauchert TR. Piezothermo-elastic composite plate analysis using the first order shear deformation theory. Comput Struct 1994;51:79–89. [17] Carrera E. C° Reissner–Mindlin multilayered plate element including Zig-Zag and interlaminar stress continuity. Int J Numer Methods Eng 1996;39:1797– 820.

135

[18] Auricchio F, Sacco E. Refined first-order shear deformation theory models for composite laminates. J Appl Mech 2003;70:381–90. [19] Hosseini-Hashemi S, Rokni Damavandi Taher H, Akhavan H, Omidi M. Free vibration of functionally graded rectangular plates using first-order shear deformation plate theory. Appl Math Model 2010;34:1276–91. [20] Pandya BN, Kant T. Higher order shear deformable theories for flexure of sandwich plates-finite element evaluations. Int J Solids Struct 1988;24:1267– 86. [21] Reddy JN. A simple higher-order theory for laminated composite plates. J Appl Mech 1984;51:745–52. [22] Shi G. A new simple third-order shear deformation theory of plates. Int J Solids Struct 2007;44:4399–417. [23] Nguyen-Xuan H, Thai CH, Nguyen-Thoi T. Isogeometric finite element analysis of composite sandwich plates using a higher order shear deformation theory. Compos Part B Eng 2013;55:558–74. [24] Nguyen Tuan N, Hui D, Lee J, Nguyen-Xuan H. An efficient computational approach for size-dependent analysis of functionally greaded nanoplates. Comput Methods Appl Mech Eng 2015;297:191–218. [25] 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. [26] Thai CH, Ferreira AJM, Bordas SPA, Rabczuk T, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. Eur J Mech A/Solids 2014;43:89–108. [27] Thai Chien H, Kulasegaram S, Tran Loc V, Nguyen-Xuan H. Generalized shear deformation theory for functionally graded isotropic and sandwich plates based on isogeometric approach. Comp and Struct. 2014;141:94–112. [28] Soldatos KP. A transverse shear deformation theory for homogenous monoclinic plates. Acta Mech 1992;94:195–220. [29] Karama M, Afaq KS, Mistou S. Mechanical behavior of laminated composite beam by new multi-layered laminated composite structures model with transverse shear stress continuity. Int J Solids Struct 2003;40:1525–46. [30] Arya H, Shimpi RP, Naik NK. A zigzag model for laminated composite beams. Compos Struct 2002;56:21–4. [31] Aydogdu M. A new shear deformation theory for laminated composite plates. Compos Struct 2009;89:94–101. [32] Reddy JN. Analysis of functionally graded plates. Int J Numer Methods Eng 2000;47:663–84. [33] Thai HC, Tran VL, Tran TD, Nguyen-Thoi T, Nguyen-Xuan H. Analysis of laminated composite plates using higher-order shear deformation plate theory and node-based smoothed discrete shear gap method. Appl Math Model 2012;36:5657–77. [34] Nguyen VP, Rabczuk T, Bordas S, Duflot M. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79:763–813. [35] Belytschko T, Lu YY, Gu L. Element free Galerkin method. Int J Numer Methods Eng 1994;37:229–56. [36] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng 1999;45:601–20. [37] Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. Int J Numer Methods Fluids 1995;20:1081–106. [38] Randles PW, Libersky LD. Recent improvements in SPH modeling of hypervelocity impact. Int J Impact Eng 1997;20:525–32. [39] Johnson GR, Beissel SR. Normalized smoothing functions for SPH impact computations. Int J Numer Methods Eng 1996;39:2725–41. [40] Atluri SN, Shen S. The meshless local Petrov–Galerkin (MLPG) method: a simple and less-costly alternative to the finite element and boundary element methods. Comput Model Eng Sci 2002;3:11–51. [41] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [42] Shepard D. A two-dimensional function for irregularly spaced points, in: Proceedings of the 23rd ACM national conference; 1968. pp. 517–524. [43] Lancaster P, Salkauskas K. Surfaces generated by moving least squares methods. Math Comput 1981;37:141–58. [44] Fernandez-Mendez S, Huerta A. Imposing essential boundary conditions in mesh-free methods. Comput Methods Appl Mech Eng 2004;193:1257–75. [45] Tran Loc V, Thai Chien H, Nguyen-Xuan H. An isogeometric finite element formulation for thermal buckling analysis of functionally graded plates. Finite Elem Anal Des 2013;73:65–76. [46] Nguyen-Xuan H, Tran Loc V, Thai Chien H, Kulasegaram S, Bordas SPA. Isogeometric analysis of functionally graded plates using a refined plate theory. Compos Part B 2014;64:222–34. [47] Tran Loc V, Ly Hung Anh, Lee J, Abdel-Wahab M, Nguyen-Xuan H. Vibration analysis of cracked FGM plates using higher-order shear deformation theory and extended isogeometric approach. Int J Mech Sci 2015;96:65–78. [48] Cho JY, Song YM, Choi YH. Boundary locking induced by penalty enforcement of essential boundary conditions in mesh-free methods. Comput Methods Appl Mech Eng 2008;131:1167–83. [49] Krysl P, Belytschko T. Analysis of thin plates using element-free Galerkin methods. Comput Mech 1995;17:26–35. [50] Liew MK, Wang J, Ng TY, Tan MJ. Free vibration and buckling analyses of sheardeformable plates based on FSDT meshfree method. J Sound Vib 2004;276:997–1017. [51] Liu GR, Wang JG. A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 2002;54:1623–48.

136

C.H. Thai et al. / Engineering Analysis with Boundary Elements 64 (2016) 122–136

[52] Liu GR, Yan L, Wang JG, Gu YT. Point interpolation method based on local residual formulation using radial basis functions. Struct Eng Mech 2002;14:713–32. [53] Dai KY, Liu GR, Lim KM, Han X, Du SY. A meshfree radial point interpolation method for analysis of functionally graded material (FGM) plates. Comput Mech 2004;34:213–23. [54] Dai KY, Liu GR, Lim KM, Gu YT. Comparison between the radial point interpolation and the kriging interpolation used in meshfree methods. Comput Mech 2003;32:60–70. [55] Liu GR, Zhang GY, Dai KY, Wang YY, Zhong ZH, Li GY, Han X. A linearly conforming point interpolation method (LC-PIM) for 2D solid mechanics problems. Int J Comput Methods 2005;2:645–65. [56] Liu GR. A generalized gradient smoothing technique and the smoothed bilinear form for galerkin formulation of a wide class of computational methods. Int J Comput Methods 2008;5:199–236. [57] Zhao X, Liu GR, Dai KY, Zhong ZH, Li GY, Han X. Free-vibration analysis of shells via a linearly conforming radial point interpolation method (LC-RPIM). Finite Elem Anal Des 2009;45:917–24. [58] Chen L, Nguyen-Xuan H, Nguyen-Thoi T, Wu SC. Assessment of smoothed point interpolation methods for elastic mechanics. Int J Numer Methods Biomed Eng 2010;26(12):1635–55. [59] Cui XY, Liu GR, Li GY. A cell-based smoothed radial point interpolation method (CS-RPIM) for static and free vibration of solids. Eng Anal Bound Elem 2010;34:144–57. [60] Liu GR, Jiang Y, Chen L, Zhang GY, Zhang YW. A singular cell-based smoothed radial point interpolation method for fracture problems. Comput Struct 2011;89:1378–96. [61] Gu L. Moving Kriging interpolation and element free Galerkin method. Int J Numer Methods Eng 2003;56:1–11. [62] Bui TQ, Nguyen TN, Nguyen-Dang H. A moving Kriging interpolation-based meshless method for numerical simulation of Kirchhoff plate problems. Int J Numer Methods Eng 2009;77:1371–95. [63] Vel SS, Batra RC. Three-dimensional exact solution for the vibration of functionally graded rectangular plates. J Sound Vib 2004;272:703–30. [64] Qian LF, Batra RC, Chen LM. Static and dynamic deformations of thick functionally graded elastic plate by using higher-order shear and normal deformable plate theory and meshless local Petrov–Galerkin method. Compos Part B: Eng 2004;35:685–97.

[65] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall 1973;21:571–4. [66] Touratier M. An efficient standard plate theory. Int J Eng Sci 1991;29:745–52. [67] Zenkour AM. Generalized shear deformation theory for bending analysis of functionally graded plates. Appl Math Model 2006;30:67–84. [68] Dolbow J, Belytschko T. An introduction to programming the meshless element-free Galerkin method. Arch Comput Methods Eng 1998;5:207–41. [69] Vel SS, Batra RC. Exact solution for thermoelastic deformations of functionally graded thick rectangular plates. AIAA J 2004;40:1021–33. [70] Carrera E, Brischetto S, Cinefra M, Soave M. Effects of thickness stretching in functionally graded plates and shells. Compos Part B: Eng 2011;42:123–33. [71] Brischetto S. Classical and mixed advanced models for sandwich plates embedding functionally graded cores. J Mech Mater Struct 2009;4:13–33. [72] Neves AMA, Ferreira AJM, Carrera E, Roque CMC, Cinefra M, Jorge RMN. A quasi-3D sinusoidal shear deformation theory for the static and free vibration analysis of functionally graded plates. Compos Part B: Eng 2012;43:711–25. [73] Neves AMA, Ferreira AJM, Carrera E, Cinefra M, Roque CMC, Jorge RMN. 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: Eng 2013;44:657–74. [74] Gilhooley DF, Batra RC, Xiao JR, McCarthy MA, Gillespie JW. Analysis of thick functionally graded plates by using higher-order shear and normal deformable plate theory and MLPG method with radial basis functions. Comput Struct 2007;80:539–52. [75] Tran VL, Ferreira AJM, Nguyen-Xuan H. Isogeometric analysis of functionally graded plates using higher-order shear deformation theory. Compos Part B: Eng 2013;51:368–83. [76] Li Q, Lu V, Kou K. Three-dimensional vibration analysis of functionally graded material sandwich plates. J Sound Vib 2008;311:498–515. [77] Ma LS, Wang TJ. Relationship between axisymmetric bending and buckling solutions of FGM circular plates based on third-order plate theory and classical plate theory. Int J Solids Struct 2004;41:85–101. [78] Saidi AR, Rasouli A, Sahraee S. Axisymmetric bending and buckling analysis of thick functionally graded circular plates using unconstrained third-order shear deformation plate theory. Compos Struct 2009;89:110–9. [79] Zenkour AM. A comprehensive analysis of functionally graded sandwich plates: Part 2 buckling and free vibration. Int J Solids Struct 2005;42:5243–58.