A simple FSDT-based meshfree method for analysis of functionally graded plates

Engineering Analysis with Boundary Elements 79 (2017) 1–12

Contents lists available at ScienceDirect

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

A simple FSDT-based meshfree method for analysis of functionally graded plates

MARK

⁎

Tan-Van Vua, , Ngoc-Hung Nguyenb, Amir Khosravifardc, M.R. Hematiyanc, Satoyuki Tanakad, Tinh Quoc Buie,f,⁎⁎ a

Faculty of Civil Engineering, Ho Chi Minh University of Architecture, 196 Pasteur, District 3, Ho Chi Minh City 70000, Vietnam Faculty of Civil Engineering, Thu Dau Mot University, 6 Tran Van On, Binh Duong Province 70000, Vietnam c Department of Mechanical Engineering, Shiraz University, Shiraz 71936, Iran d Graduate School of Engineering, Hiroshima University, Japan e Institute for Research and Development, Duy Tan University, Da Nang City, Vietnam f Tokyo Institute of Technology, Department of Civil and Environmental Engineering, 2-12-1-W8-22, Ookayama, Meguro-ku, Tokyo 152-8552, Japan b

A R T I C L E I N F O

A BS T RAC T

Keywords: Functionally graded material Meshfree methods First order shear deformation theory Moving Kriging interpolation Plate

Modeling of mechanical behavior of plates has been accomplished in the past decades, with different numerical strategies including the finite element and meshfree methods, and with a range of plate theories including the first-order shear deformation theory (FSDT). In this paper, we propose an efficient numerical meshfree approach to analyze static bending and free vibration of functionally graded (FG) plates. The kinematics of plates is based on a novel simple FSDT, termed as S-FSDT, which is an effective four-variable refined plate theory. The S-FSDT requires C1-continuity that is satisfied with the basis functions based on moving Kriging interpolation. Some major features of the approach can be summarized: (a) it is less computationally expensive due to having fewer unknowns; (b) it is naturally free from shear-locking; (c) it captures the physics of sheardeformation effect present in the conventional FSDT; (d) the essential boundary conditions can straightforwardly be treated, the same as the FEM; and (e) it can deal with both thin and thick plates. All these features will be demonstrated through numerical examples, which are to confirm the accuracy and effectiveness of the proposed method. Additionally, a discussion on other possible choices of correlation functions used in the model is given.

1. Introduction Functionally graded (FG) materials are heterogeneous composite materials in which the material properties are spatially varied in a predetermined direction. Due to the superior advantages of mechanical behaviors, FG plates are widely used in numerous applications in advanced engineering fields including automotive, aerospace, civil, infrastructure, nuclear power plants, and many others [1]. It is widely agreed that the most popular methods for modeling the deformation of a plate are based on two theories, the classical plate theory (CPT) and the first-order shear deformable theory (FSDT). Since in the CPT, the transverse shear deformation is neglected, only thin plates can be treated by this theory. In the contrary, the deflection and rotations are independently defined in the FSDT so that the shear-deformation effect can be taken into account to accurately capture the physical behavior of thick plates. In the theory of the FSDT, only C°-continuity is required

⁎

for formulating interpolation functions of dependent variables [2–6]. The FSDT coupled with the finite element method (FEM) to analyze mechanical behaviors of FG plates has been successfully utilized and numerous works are available in the open literature, e.g., see [2–7]. However, it is well-known that the FSDT based FEM approaches using low-order standard shape functions without special treatments inherently produce inaccurate shear strains when dealing with thin plates, due to the shear-locking phenomenon. This phenomenon is mainly due to the inability of simulating Kirchhoff-type constraints in thin plates. To circumvent this awkward and difficult shortcoming, different approaches with additional expense and implementation complexity have been proposed including the MITC family of elements [8], the assumed strain method [9], field consistent approach [10], and isogeometric finite element [11]. Recently, Thai and Choi [12] reformulated the FSDT into a more simple but effective form, known as simple first-order shear deformation theory (S-FSDT), which is accomplished by decomposing of

Corresponding author at: Ho Chi Minh City University of Architecture, Vietnam. Corresponding author at: Institute for Research and Development, Duy Tan University, Da Nang City, Vietnam & Tokyo Institute of Technology, Japan. E-mail addresses: [email protected] (T.-V. Vu), [email protected], [email protected] (T.Q. Bui).

⁎⁎

http://dx.doi.org/10.1016/j.enganabound.2017.03.002 Received 29 November 2016; Received in revised form 1 February 2017; Accepted 18 March 2017 0955-7997/ © 2017 Elsevier Ltd. All rights reserved.

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

the transverse displacements into bending and shear parts. In this formulation, the variable defining the rotation degree of freedom is eliminated in terms of the bending transverse displacements. As a result, the weak form contains only four variables rather than the usual five, and hence the computational expenses are reduced. Although the FEM is being extensively used in various analyses of FG plates, some inherent drawbacks from the FEM however exits, for instance, the cumbersome tasks of mesh generation; a high order, i.e., C1 consistency, a time-consuming procedure for element connectivity; and re-meshing in moving boundary problems. Therefore, a family of meshfree methods based on a set of scattered nodes rather than mesh of elements has been introduced in literature [13], which aim at overcoming the aforementioned drawbacks, especially in mesh distortion and remeshing problems. Interestingly, by introducing basis functions with natural C1 continuity, the numerical meshfree methods will be free from shear-locking, and consequently are capable of accurately handling the physical phenomenon with straightforward numerical implementation. Recently, Chen et al. [14] presented the meshless local natural neighbor interpolation method, for free vibration analysis of moderately thick FG plates, based on the FSDT. Using the three-node triangular FEM shape functions as the test function, the order of integrands involved in domain integrals are reduced in the proposed method. This approach is known as a special case of the generalized meshless local Petrov–Galerkin method. It has been applied to the analysis of FG viscoelastic materials [15] and also nonlinear steady and transient heat transfer of two-dimensional structures [16]. The main objective of this paper is to develop an efficient but accurate numerical approach based on meshfree moving Kriging (MK) method [17] incorporated with the new effective plate theory [12], the S-FSDT, to analyze mechanical behaviors of composite FG plates. One major advantage of MK method is the convenient way of imposing the essential boundary conditions since the MK shape functions automatically possess the Kronecker's delta property where the standard meshfree methods own many difficulties. This is a great benefit to treat the essential boundary conditions in a similar manner as the FEM, consequently its implementation is straightforward. Here in this work, we are particularly interested in applying our proposed approach, the S-FSDT based MK meshless method, to free vibration and static deflection analysis of FG plates. The present method not only avoids the shear-locking phenomenon, handles the essential boundary conditions effectively, is less computationally expensive, but also accurately deals with both thin and thick plates. The rest of the paper is structured as follows. A brief review of the FG plates is given in Section 2, followed by the formulation of the SFSDT in Section 3. Meshfree formulation for bending and free

Fig. 1. Geometry notation and coordinates of an FG plate.

vibration analyses of plates is presented in Section 4. Numerical validations and discussions are presented in Section 5. Other possible choices of correlation functions are addressed in Section 6. Some conclusions drawn from the study are given in Section 7. 2. Functionally graded plates Consider an FG plate as shown in Fig. 1, which is made from a mixture of two material phases, for example, a metal and a ceramic with thickness h . The bottom and top faces of the plate are assumed to be fully metallic and ceramic, respectively [1–5]. In this study, the Poisson's ratio ν is assumed to be constant for simplicity, whereas the Young's modulus E and the density ρ are assumed to vary continuously through the thickness. the effective material properties are computed by a power law distribution with Voigt's rule of mixtures. The Young's modulus E (z ) and mass density ρ (z ) are hence given by

E (z ) = Em + (Ec − Em ) Vc

(1)

ρ (z ) = ρm + (ρc − ρm ) Vc

(2)

where the subscripts m and c represent the metallic and ceramic constituents, respectively; Vc = (0.5 + z / h )n is the volume fraction of the ceramic; and n is the gradient index, which governs the volume fraction gradation. Fig. 2 shows the variation in the ceramic volume Vc with respect to the thickness ratio z / h for different values of the index n [1]. 3. Formulation of the simple first-order shear deformation theory (S-FSDT) In this section, we briefly present the S-FSDT, for a rigorous explanation of the theory one can refer to [12,15]. In the standard

Fig. 2. Variation of ceramic volume fraction Vc versus the thickness ratio z /h for different values of the index n .

2

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

FSDT the three-dimensional displacement field (u , v, w ) can be expressed in terms of five unknown variables as follows:

u (x, y, z ) = u 0 (x, y) + zϕx (x, y)

(3)

v (x, y, z ) = v0 (x, y) + zϕy (x, y)

(4)

w (x, y, z ) = w0 (x, y)

(5)

(6)

Table 1 Material properties of the FG plates used for the analysis. Properties

E (GPa ) υ ρ (kg /m3)

70 0.3 2707

Alumina ( Al2 O3 )

Zirconia (ZrO2 )

380 0.3 3800

200 0.3 5700

w (x, y, z ) = wb (x, y) + ws (x, y)

(8)

⎧ ⎫ ∂u 0 ∂ 2w − z 2b ⎪ ⎪ ∂x ∂x ⎪ ⎧ εx ⎫ ⎪ ∂v0 ∂ 2w ⎪ − z 2b ⎪ εy ⎪ ⎪ ∂x ∂y ⎪ ⎪γ ⎪ ⎪ ⎪ ⎪ 2w ∂ u ∂ v ∂ ⎨ xy ⎬ = ⎨ 0 + 0 − 2z b ⎬ ∂ y ∂ x ∂ x ∂ y ⎪ γxz ⎪ ⎪ ⎪ ⎪γ ⎪ ⎪ ⎪ ∂ws ⎩ yz ⎭ ⎪ ⎪ ∂x ⎪ ⎪ ∂ws ⎪ ⎪ ∂y ⎭ ⎩

Ceramic

Aluminum ( Al )

(7)

In contrast to the standard FSDT, the displacement fields in Eqs. (6)–(8) for the S-FSDT contain only four unknowns, u 0 (x, y),v0 (x, y),wb (x, y),ws (x, y). Because the rotations are obtained from partial derivatives of the bending component wb (x, y), discretization of the S-FSDT is inherently free from the shear-locking. The three physically relevant boundary conditions for the S-FSDT are the clamped, simply supported and free condition, as in the Kirchhoff-Love or classical plate theory. In contrast, the FSDT has five physically relevant boundary conditions, hard clamped, soft clamped, hard simply supported, soft simply supported and free. Interested readers are referred to Arnold and Falk [18] for an in-depth discussion on this subject. Although the boundary conditions in the FSDT are more descriptive than the S-FSDT in practices, we have not found this issue to be a problem in our work as we have demonstrated in the numerical examples that good agreements with the FSDT are achieved. By making the usual small strain assumptions, the strain–displacement relations can be expressed as

where u 0 (x, y),v0 (x, y), and w0 (x, y) are the displacements at the midplane of a plate in the x, y, z directions, respectively; ϕx (x, y), ϕy (x, y) are the rotations of a transverse normal along the y, x axes. To derive the simple FSDT (S-FSDT), the following assumptions are made: (a) the transverse displacement w0 (x, y) is divided into bending wb (x, y) and shear ws (x, y) components, i.e., w0 (x, y) = wb (x, y) + ws (x, y); (b) the rotations in the FSDT is expressed in terms of the bending component only ϕx = −∂wb (x, y)/∂x , ϕy = −∂wb (x, y)/∂y . Finally, Eqs. (3–5) can be rewritten as follows:

u (x, y, z ) = u 0 (x, y) − z∂wb (x, y)/∂x

v (x, y, z ) = v0 (x, y) − z∂wb (x, y)/∂y

(9)

Table 2 Normalized deflection of an Al /ZrO2 square plate with different length-thickness ratio a /h and gradient indices n. Boundary condition

a /h

Method

n=0

n = 0.5

n=1

n=2

SSSS

5

S-FSDT based IGA [20] FSDT based IGA [20] FSDT based kp-Ritz [21] FSDT based ES-DSG3 [22] Present

0.1717 0.1717 0.1722 0.1703 0.1723

0.2324 0.2324 0.2403 0.2232 0.2331

0.2719 0.2719 0.2811 0.2522 0.2723

0.3115 0.3115 0.3221 0.2827 0.3116

20

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.1440 0.1440 0.1450

0.1972 0.1972 0.1982

0.2310 0.2310 0.2318

0.2628 0.2628 0.2634

100

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.1423 0.1423 0.1432

0.1949 0.1949 0.1960

0.2284 0.2284 0.2292

0.2597 0.2597 0.2603

5

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.3164 0.3175 0.3189

0.4299 0.4314 0.4322

0.5032 0.5049 0.5043

0.5752 0.5772 0.5751

20

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.2800 0.2803 0.2834

0.3835 0.3838 0.3869

0.4493 0.4497 0.4517

0.5109 0.5114 0.5125

100

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.2777 0.2777 0.2777

0.3805 0.3805 0.3805

0.4458 0.4458 0.4458

0.5068 0.5069 0.5069

5

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.5083 0.5089 0.5053

0.6918 0.6926 0.6874

0.8099 0.8108 0.8042

0.9247 0.9258 0.9177

20

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.4614 0.4615 0.4595

0.6319 0.6321 0.6289

0.7404 0.7406 0.7363

0.8420 0.8422 0.8368

100

S-FSDT based IGA [20] FSDT based IGA [20] Present

0.4584 0.4584 0.4566

0.6281 0.6281 0.6251

0.7360 0.7360 0.7319

0.8367 0.8367 0.8316

SFSS

SFSF

3

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

Table 3 Normalized deflection of Al /Al2 O3 thin plate under different gradient indices and boundary conditions. Boundary condition

Method

n=0

n = 0.5

n=1

n=2

n=5

n = 10

SFSF

S-FSDT [20] FSDT [20] FSDT-MK Present

1.4302 1.4302 1.4250 1.4245

2.2062 2.2062 2.1951 2.1944

2.8692 2.8693 2.8490 2.8481

3.6770 3.6770 3.6432 3.6419

4.3483 4.3483 4.3118 4.3100

4.7740 4.7740 4.7442 4.7421

SSSS

S-FSDT [20] FSDT [20] FSDT-MK Present

0.4438 0.4438 0.4473 0.4469

0.6846 0.6847 0.6878 0.6873

0.8904 0.8904 0.8904 0.8899

1.1411 1.1411 1.1356 1.1348

1.3494 1.3494 1.3454 1.3443

1.4816 1.4816 1.4843 1.4830

SCSC

S-FSDT [20] FSDT [20] FSDT-MK Present

0.2096 0.2097 0.2106 0.2105

0.3232 0.3234 0.3240 0.3239

0.4204 0.4205 0.4198 0.4196

0.5387 0.5389 0.5357 0.5355

0.6372 0.6375 0.6346 0.6342

0.6996 0.7000 0.6997 0.6992

CCCC

S-FSDT [20] FSDT [20] FSDT-MK Present

0.1384 0.1384 0.1384 0.1384

0.2135 0.2135 0.2126 0.2126

0.2776 0.2776 0.2750 0.2750

0.3557 0.3558 0.3503 0.3503

0.4208 0.4209 0.4152 0.4152

0.4621 0.4622 0.4586 0.4586

(a) SFSF

(b) SSSS

(c) SCSC

(d) CCCC

Fig. 3. Deflection of an Al /Al2 O3 thin plate with different boundary conditions and gradient index n = 10 : (a) SFSF, (b) SSSS, (c) SCSC, (d) CCCC.

or in matrix form

ε=

ε0 0

{} { } +

− zκ γ

⎧ ∂u 0 ∂x ⎪ ⎪ ∂v0 ε0 = ⎨ ∂y ⎪ ⎪ ∂u 0 + ⎩ ∂y

(10)

with

⎫ ⎪ ⎪ ⎬, ⎪ ∂v0 ⎪ ∂x ⎭

⎧ ∂ 2wb ⎫ ⎪ ∂x 2 ⎪ ⎪ ∂ 2w ⎪ κ = ⎨ 2b ⎬, ⎪ ∂y ⎪ ⎪ 2 ∂ 2wb ⎪ ⎩ ∂x ∂y ⎭

⎧ ∂ws ⎫ ⎪ ∂x ⎪ γ = ⎨ ∂ws ⎬ ⎪ ∂y ⎪ ⎩ ⎭ (11)

The constitutive relations are derived from Hooke's law by the following equation: 4

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al. q

uh (x) =

∑ ϕI (x ) uI

(16)

I =1

where ϕI (x ) are the MK shape functions defined by m

ϕI (x) =

q

∑ pj (x) AjI

+

j =1

∑ rk (x) BkI (17)

k =1

and matrices A and B are calculated by

A = (PT R−1P)−1PT R−1

B=

R−1(I

(18) (19)

− PA)

where I is the unit matrix, and the vector p (x) in Eq. (15) is the polynomial with m basis functions

pT (x ) = [ p1 (x), p2 (x), p3 (x)...., pm (x)]

The matrix P (q × m ) is the collected values of the polynomial basis functions

Fig. 4. The computational time for solving the bending deflection of a SSSS FG plate subjected to a uniform load using the developed method.

σ = Dm (z )(ε0 − zκ ),

τ = Ds (z ) γ

⎡ p1 (x1) p2 (x1) ⎢ p (x ) p2 (x2) P=⎢ 1 2 ⎢ ⋮ ⋮ ⎢ p (x ) p (x ) ⎣ 1 q 2 q

(12)

with

σ = [ σx σy τxy ]T ,

τ = [ τxz τyz ]T

(20)

(13)

⋯ pm (x1) ⎤ ⎥ ⋯ pm (x2) ⎥ ⋱ ⋮ ⎥ ⋯ pm (x q)⎥⎦

(21)

and the term r (x) in Eq. (15) is given by

and

r T (x) = [R (x1, x), R (x2, x), ....R (x q , x)]

⎡1 υ ⎤ 0 E (z ) ⎢ ⎥ , D (z ) = kE (z ) ⎡ 1 0 ⎤ Dm (z ) = 1 0 υ s 2(1 + υ) ⎢⎣ 0 1 ⎥⎦ 1 − υ2 ⎢⎣ 0 0 (1 − υ)/2 ⎥⎦

where R (xi , xj) is the correlation function between pairs of the nodes at xi and xj , it is the covariance of the value u (x): R (xi , xj) = cov[u (xi), u (xj)] and R (xi , x) = cov[u (xi), u (x)]. In this study, the widely used Gaussian function is adopted as the correlation function, and a correlation parameter θ > 0 is introduced to best fit the model [17,19]. In this work, unless stated otherwise, we simply use θ = 3 throughout the paper.

(14)

where k is the shear correction factor, and k = 5/6 is taken in this work. 4. Meshless formulation for bending and free vibration of plates

(22)

2

R (xi , xj) = e−θrij 4.1. The moving Kriging shape function

x2

y2

where r ij = xi − xj . The quadratic basic = [1 x y xy ] is employed for the numerical analysis. The correlation matrix R[R (xi , xj)]q × q is given by

The MK shape functions interpolation technique and their derivatives are briefly introduced in this subsection. For thorough description of the method and its mathematical properties one may refer to [17,19]. Consider the distribution function u (x) in a sub-domain Ωx , so that Ωx ⊆ Ω . Assuming that its values can be interpolated based on the nodal values at xi (i ∈ [1, q]), wherein q is the total number of the nodes in Ωx , the approximated function uh (x) can be expressed as

uh (x) = [pT (x) A + r T (x) B] u (x)

(23)

pT (x)

⎡ 1 R (x1, x2) ⎢ R ( x , x ) 1 2 1 R [R (xi , xj )] = ⎢ ⎢ ⋮ ⋮ ⎢ R (x , x ) R (x , x ) q 1 q 2 ⎣

⋯ R (x1, x q) ⎤ ⎥ ⋯ R (x2, x q)⎥ ⎥ ⋱ ⋮ ⎥ ⋯ 1 ⎦

(24)

For thin plate problems, not only the first-order derivatives of shape functions are required, but also the second-order derivatives are needed to be computed, these derivatives are obtained by direct

(15)

or

Fig. 5. Convergence curves of the central deflection for different values of the scaling factor α, in a SSSS FG with a /h = 100 , and gradient index n = 1.

5

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

differentiation of Eq. (17), as follows: m

ϕI . i (x) =

∑ pj,i (x) AjI

+

j

∑ rk,i (x) BkI

⎡I I ⎤ m = ⎢ 0 1 ⎥, ⎣ I1 I2 ⎦

q

∑ pj,ii (x) AjI

+

∑ rk,ii (x) BkI

j

(26)

k

⎧ uh ⎫ ⎪ ⎪ u1 = ⎨ v h ⎬ = ⎪ w h + wh ⎪ ⎩ b s ⎭

∫−h/2 ρ (z)(1, z, z2) dz,

⎡u ⎤ u = ⎢ u1 ⎥ ⎣ 2⎦

(37)

⎧ ∂w h /∂x ⎫ ⎪ b ⎪ u2 = ⎨ ∂w h /∂y ⎬ = ⎪ b ⎪ ⎩ 0 ⎭

q

∑ N1I uI , 1

q

∑ N2I uI I

(38)

and

(27)

where dc defines the characteristic length relative to the nodal spacing near the point of interest, and α is a scaling factor. It should be noted that the shape function ϕI (x) at node xj possesses the delta function property.

⎧1 for I = j ϕI (xj) = δIj = ⎨ ⎩ 0 for I ≠ j

h /2

(I0, I1, I2 ) =

with

In meshfree approaches [19], the influence domain is usually a circle or sphere, defined by a radius and centered at the point of interest. This domain is used to determine the scattered nodes which are used for interpolation. The following expression is taken to compute the size of the support domain

d m = αdc

(36)

where

(25)

k

m

ϕI , ii (x) =

∫Ω δ εT DεdΩ + ∫Ω δ γT DsγdΩ = ∫Ω δ uT mu¨ dΩ

q

N1I

⎡ ϕI 0 0 0 ⎤ ⎢ ⎥ = ⎢ 0 ϕI 0 0 ⎥ , ⎢⎣ 0 0 ϕI ϕI ⎥⎦

N2I

⎡ 0 0 ϕI , x 0 ⎤ ⎢ ⎥ = ⎢ 0 0 ϕI , y 0 ⎥ ⎢⎣ 0 0 0 0 ⎥⎦

(39)

By substituting Eqs. (31) and (34) into Eqs. (33) and (36) the formulations of the static, free vibration are rewritten in the following form:

(28)

(K − ω 2 M) u = 0

Ku = F,

An important point with regard to the effect of the correlation parameter on the quality of the shape functions could be highlighted. However, it is not related to any physical aspect of the problem. Until now, there exist no exact theoretical rules to find out the optimal value of the correlation parameter for all problems. In practice, its reasonable value could be determined by numerical experiments to ensure the accuracy and consistency of the solution based on the properties of each type of problems. For more information, curious readers may refer to [19]. Nonetheless, other possible choices of the correlation functions will be discussed in Section 6 of this paper.

(40)

where the global stiffness matrix K is computed through

K=

⎧

m ⎫T ⎡

∫Ω ⎨⎩ BBb ⎬⎭

Bm ⎫ Dm B ⎤ ⎧ ⎢ ⎥ ⎨ ⎬ dΩ + ⎣ B D b ⎦ ⎩ Bb ⎭

∫Ω (Bs)T DsBsdΩ

(41)

The load vector F is computed as follows:

F=

∫Ω

f NdΩ

where

NI = [0 0 ϕI ϕI ]T

(42)

The global mass matrix M is expressed as:

⎧

T 1⎫

⎧

1⎫

∫Ω ⎨⎩ NN2 ⎬⎭ m ⎨⎩ NN2 ⎬⎭ dΩ

4.2. Discrete governing equations

M=

In the parametric domain in terms of meshfree method, the generalized displacements in the middle plane of the plate are approximated by Eq. (15), in which

In this paper, a background cell with 16 Gaussian points is used for the numerical integration. It should be pointed out that the second order derivatives of the shape functions should be evaluated in the weak formulation of the S-FSDT. This is in contrast to the original FSDT which requires only first order derivatives. This might lead to a loss of accuracy in the FEM, however, since the meshfree methods are generally capable of evaluating higher order derivatives with better accuracy, it does not pose a difficulty, especially if the support domain size is selected appropriately. This point is assessed by the numerical examples, where accurate results are obtained by the proposed meshfree method.

T

uh = [ uh v h wbh wsh ]

(29)

and

uI = [ uI vI wbI wsI ]T

(30)

By substituting Eq. (15) into Eq. (11), one can obtain: q

ε0 =

q

∑ BmI uI ,

κ=

I

q

∑ BbI uI ,

γ=

I

∑ BsI uI

(31)

I

(43)

5. Numerical results and discussions

with

BmI

⎡ ϕI , x 0 0 0 ⎤ ⎡ 0 0 ϕI , xx 0 ⎤ ⎡ 0 0 0 ϕI , x ⎤ ⎢ ⎥ ⎢ ⎥ ⎥ = ⎢ 0 ϕI , y 0 0 ⎥ , BbI = ⎢ 0 0 ϕI , yy 0 ⎥B,mI = ⎢ ⎣⎢ 0 0 0 ϕI , y ⎥⎦ ⎢ϕ ⎥ ⎢0 0 ϕ ⎥ 0 0 0 ϕ I , xy ⎣ I ,y I ,x ⎦ ⎣ ⎦

To illustrate the accuracy of the present approach described in the previous sections, the static bending and free vibration of homogeneous and FG plates with square shape are considered and analyzed. Two types of FG plates of Al / ZrO2 and Al / Al2 O3 are used throughout this work. The material properties of these FG plates are listed in Table 1. The Young's modulus and mass density are evaluated using Voigt's rule of mixtures, see Eqs. (1–2). For convenience, the boundaries of the FG plate are denoted as follows: completely free (F), simply supported (S) or fully clamped (C) edges. Both thin and thick plates are considered through the specified thickness-span aspect ratios. Unless stated otherwise a regular set of 13×13 scattered nodes is used throughout the analyses.

(32)

For the static problem, the weak form can be expressed as follows:

∫Ω δ εT DεdΩ + ∫Ω δ γT DsγdΩ = ∫Ω δ (wb + ws ) fdΩ

(33)

where f is the transverse loading per unit area and

⎡ε ⎤ ε = ⎢ 0 ⎥, ⎣κ⎦

⎡ m B⎤ D = ⎢D ⎥, ⎣ B Db ⎦

Ds =

h /2

∫−h/2 Ds (z) dz

(34)

with

Dm =

h /2

∫−h/2 Dm (z) dz

h /2

B=−

∫−h/2 z Dm (z) dz Db

5.1. Static bending analysis

h /2

=

∫−h/2 z2Dm (z) dz

The first example deals with an Al / ZrO2 square plate subjected to a uniform load P = 1N / mm2 . Different boundary conditions such as SSSS, SFSS, SFSF, various values of the length to thickness ratios a/h =5, 20

(35) For the free vibration, the weak form can be expressed: 6

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

Table 4 First normalized natural frequency of an Al /Al2 O3 square plate for a /h = 2, 10, 20 .

a /h

Method

2

Analytical solution

n=0

n = 0.5

n=1

n=4

n = 10

2D-HOT [20] S-HSDT [20] FSDT-IGA [20]

0.9400 0.9297 0.9265 0.9270

0.8232 0.8110 0.8060 0.8070

0.7476 0.7356 0.7330 0.7350

0.5997 0.5924 0.6111 0.6136

0.5460 0.5412 0.5640 0.5652

2D-HOT [20] S-HSDT [20] FSDT-IGA [20]

0.0578 0.0577 0.0577 0.0575

0.0492 0.0490 0.0490 0.0489

0.0443 0.0442 0.0442 0.0442

0.0381 0.0381 0.0382 0.0383

0.0364 0.0364 0.0366 0.0366

2D-HOT [20] FSDT-Ritz [20] FSDT – IGA [20]

0.0148 0.0146 0.0148 0.0148

0.0125 0.0124 0.0125 0.0125

0.0113 0.0112 0.0113 0.0111

0.0098 0.0097 0.0098 0.0098

0.0094 0.0093 0.0094 0.0094

Present 10

Analytical solution

Present 20

Analytical solution

Present

well with those obtained by other mentioned methods. For all considered boundary conditions and length to thickness ratios, close agreement is found between the present method, S-FSDT-based IGA and FSDT-based IGA methods. It is noted that, the deflection magnitude decreases when the boundary condition changes from SFSF, SFSS to SSSS, since the structural stiffness increases. The deflection magnitude decreases as the gradient index n reduces and also as the length to thickness ratio is increased. Increasing the volume fraction index implies that the property of plates gets more and more close to that of the metallic constituent, and consequently the stiffness of the resulting plates gradually reduces. In the next example, an Al / Al2 O3 thin plate with a length to thickness ration of a / h = 100 under a uniform load P = 1N / mm2 is considered to examine the effect of the boundary conditions and gradient index on the central deflection of the FG plate. The central deflection of this example is 10w E h3

also normalized by w = c 4c . Table 3 presents the normalized central Pa deflections obtained respectively by using the S-FSDT based IGA [20], standard FSDT based IGA [20], and the present S-FSDT based meshfree method. In addition, we also present in Table 3 numerical result of the normalized central deflection calculated by the meshfree moving Kriging method but using the classical FSDT instead, which is named as FSDTMK. Good agreements of the normalized deflections for different gradient indices and boundary conditions of FG plates among considered approaches is obtained. All these results confirm the accuracy of the developed S-FSDT based meshfree method. Similar to the previous Al / ZrO2 plate, it can be seen in this example that the magnitude of the deflection gradually increases when the boundary condition changes from CCCC and SCSC to SSSS and SFSF.

Fig. 6. Convergence curves of the first normalized natural frequency for various scaling factors α , in a SSSS FG with a /h = 10 and the gradient index n = 1.

and 100, and different values of the gradient index n=0, 0.5, 1 and 2, are considered. The normalized central deflection of the plate

w =

100wc Em h3 2 ) Pa4 12(1 − υm

is used to study the accuracy of the obtained results.

Table 2 reports the normalized central deflections, where the obtained results are compared with S-FSDT based isogeometric analysis (IGA) [20], FSDT based IGA [20] and FSDT based kp-Ritz [21], and FSDT based ES-DSG3 [22] method. Clearly, the numerical results of the normalized central deflections obtained by the present method agree

Table 5 The first five normalized natural frequencies of an Al /ZrO2 thick plate with various ratios a /h and gradient index n = 1.

a /h

Model

Mode 1

Mode 2

Mode 3

Mode 4

Mode 5

5

TSDT based Meshless [21] HSDT based Meshless [24] 3D-SSDT based Meshless [25] 3D-HSDT based Meshless [26] Present

0.2188 0.2152 0.2193 0.2192 0.2149

– 0.4114 – – 0.4284

– 0.4114 – – 0.4284

– 0.4761 – – 0.4773

– 0.4761 – – 0.4773

10

TSDT based Meshless [21] HSDT based Meshless [24] 3D-SSDT based Meshless [25] 3D-HSDT based Meshless [26] Present

0.0593 0.0584 0.0596 0.0596 0.0599

0.1431 0.1410 0.1426 0.1426 0.1435

0.1431 0.1410 0.1426 0.1426 0.1435

0.2035 0.2058 0.2058 0.2059 0.2153

0.2035 0.2058 0.2058 0.2059 0.2153

20

TSDT based Meshless [21] HSDT based Meshless [24] 3D-SSDT based Meshless [25] 3D-HSDT based Meshless [26] Present

0.0149 0.0149 0.0153 0.0153 0.0153

0.0378 0.0377 0.0377 0.0377 0.0378

0.0378 0.0377 0.0377 0.0377 0.0378

0.0595 0.0593 0.0596 0.0596 0.0586

0.0749 0.0747 0.0739 0.0739 0.0721

7

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

Table 6 The first five normalized natural frequencies of an Al /Al2 O3 thin plate with various boundary conditions and gradient indices. Method

Mode 1

Mode 2

Mode 3

Mode 4

Mode 5

0

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

56.5660 56.5526 56.5584 56.4791 56.6084

94.7565 94.6609 94.7388 94.7141 93.2207

215.6539 215.3651 215.5711 215.6299 210.6905

228.7079 228.5948 228.5829 – 229.8464

274.4598 274.1872 274.2876 – 270.4348

0.5

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

47.8972 47.8872 47.8913 47.7452 47.9540

80.2349 80.1627 80.2210 80.1576 79.0067

182.6041 182.3880 182.5386 182.4411 178.6432

193.6578 193.5819 193.5617 – 194.8400

232.3979 232.2005 232.2649 – 228.5910

1

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

43.1596 43.1511 43.1544 43.0872 43.2374

72.2984 72.2369 72.2861 72.2001 71.2856

164.5401 164.3570 164.4815 164.3911 161.2857

174.5012 174.4412 174.4179 – 175.8363

209.4085 209.2445 209.2924 – 206.5470

2

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

39.2395 39.2316 39.2347 39.1666 39.3354

65.7314 65.6747 65.7197 65.6400 64.9015

149.5922 149.4233 149.5365 149.0583 146.9400

158.6496 158.5933 158.5722 – 160.1133

190.3849 190.2330 190.2767 – 188.3231

0

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

115.9253 115.8932 115.8926 115.8695 115.7376

289.7845 289.6145 289.5806 289.7708 290.6618

289.7845 289.6145 289.5806 – 290.6618

463.5955 463.1163 463.0741 463.4781 465.9584

579.5359 579.2636 578.7215 – 580.9486

0.5

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

98.1595 98.1354 98.1343 98.0136 98.1726

245.3740 245.2532 245.2169 245.3251 247.0682

245.3740 245.2532 245.2169 – 247.0682

392.5471 392.1950 392.1448 392.4425 397.6142

490.7185 490.6291 490.0963 – 495.2124

1

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

88.4501 88.4296 88.4280 88.3093 88.4983

221.1011 221.0019 220.9643 221.0643 223.1878

221.1011 221.0019 220.9643 – 223.1878

353.7127 353.4173 353.3613 353.6252 359.3324

442.1697 442.1505 441.6348 – 448.5040

2

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

80.4160 80.3972 80.3953 80.3517 80.6550

201.0155 200.9234 200.8879 200.8793 204.0120

201.0155 200.9234 200.8879 – 204.0120

321.5761 321.3032 321.2475 321.4069 330.2422

401.9929 401.9628 401.5008 – 411.6235

0

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

170.0240 169.8890 169.9230 170.0196 169.9640

321.4690 321.1655 321.1937 321.4069 320.1753

407.1180 406.5798 406.5707 – 413.5334

555.3781 554.3949 554.5021 555.2809 555.7793

600.2165 599.7747 599.3170 – 594.9363

0.5

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

143.9675 143.8672 143.8904 143.8179 143.1032

272.2028 271.9835 271.9916 272.1090 271.0817

344.7257 344.3519 344.3090 – 342.0457

470.2635 469.5518 469.5928 470.0770 471.1181

508.2296 508.0145 507.5430 – 505.9031

1

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Exact [27] Present

129.7269 129.6422 129.6605 129.6496 129.9227

245.2758 245.0938 245.0927 245.1310 246.0544

310.6242 310.3216 310.2664 – 311.5795

423.7400 423.1484 423.1599 423.6904 429.6741

457.9482 457.8229 457.3585 – 456.0445

2

CPT-neu based IGA [20] FSDT based IGA [20]

117.9435 117.8652

222.9939 222.8253

282.4052 282.1232

385.2402 384.6922

n (a) SFSF

(b) SSSS

(c) SCSC

8

416.3375 416.2095 (continued on next page)

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

Table 6 (continued) Method

Mode 1

Mode 2

Mode 3

Mode 4

Mode 5

S-FSDT based IGA [20] Exact [27] Present

117.8818 117.8104 117.9340

222.8238 222.8111 223.2098

282.0750 – 279.5867

384.7018 385.0672 389.7883

415.7952 – 412.6522

0

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Present

211.3372 211.1203 211.1468 211.6949

431.0061 430.3532 430.3633 437.2917

431.0061 430.3532 430.3633 437.2917

635.4464 633.9937 634.1625 642.5621

772.7523 771.9821 770.8950 796.3914

0.5

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Present

178.9493 178.7886 178.8047 178.2686

364.9528 364.4940 364.4639 369.0948

364.9528 364.4940 364.4639 369.0948

538.0607 537.0131 537.0816 541.1439

654.3228 654.0679 652.9193 675.9409

1

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Present

161.2484 161.1129 161.1242 161.0227

328.8502 328.4760 328.4308 328.6780

328.8502 328.4760 328.4308 328.6780

484.8293 483.9606 483.9866 488.7393

589.5860 589.5289 588.3962 591.5320

2

CPT-neu based IGA [20] FSDT based IGA [20] S-FSDT based IGA [20] Present

146.6016 146.4765 146.4868 146.9611

298.9753 298.6270 298.5884 297.6900

298.9753 298.6270 298.5884 297.6900

440.7781 439.9730 439.9988 441.7803

536.0119 535.9244 534.9293 531.8659

n

(d) CCCC

normalized natural frequencies ω* = ωh ρc / Ec obtained by different methods for various values of gradient indices are tabulated in Table 4. In this numerical example, analytical solutions are available, see [20], and are used as reference herein. In Table 4, a remarkable agreement between the results of different approaches for the natural frequency is found for all length to thickness ratios and gradient indices. It can be concluded that the magnitude of natural frequency increases as the value of the length to thickness ratio decreases. Also, the magnitude of the first normalized natural frequency decreases as the value of the gradient index decreases. A further study on the effect of the scaling factor on the accuracy of the first natural frequency of the FG plate is also performed. Fig. 6 depicts the fundamental natural frequency of the plate, obtained with various values of α and with different sets of nodes. In the figure a wide range of the scaling factor, from 2 to 6, is considered. It is observed that as the number of nodes increases, the obtained results converge to the reference solution given by the FSDT approach [20]. Here again, it is concluded that when 2.5 ≤ α ≤ 3.5 accurate results can be obtained. The best result is obtained with α = 3.0 . It is observed that for α = 2.0 or α ≥ 4.0 , the accuracy of the results decreases. In the next example of natural frequency analysis, an Al / ZrO2 square thick plate with simply supported boundary condition, the gradient index n = 1 and various different length to thickness ratios a / h is considered. In Table 5, the first five normalized natural frequencies ω* = ωh ρm / Em calculated by the present method are compared with those from the TSDT based meshless global collocation of Ferreira et al. [23], HSDT meshless based Petrov-Galerkin of Qian et al. [22], and the quasi-3D-SSDT,-HSDT based meshless of Neves et al. [23,24]. It is also found that the present method offers acceptable results as compared with the reference ones. To further investigate the effects of the boundary conditions on the natural frequency, an Al / Al2 O3 square thin plate with a length to thickness ratio of a / h = 100 under different boundary conditions and gradient indices is analyzed. The first five normalized natural frequencies ω* = ωπ 2 (a2 / h ) ρm / Em obtained with different methods are listed in Table 6. Here again, the results obtained by the present S-FSDT based meshfree method are in a remarkable agreement with the results derived from the FSDT-based IGA [20], CPT-neu-based IGA [20], and

The deflection magnitude increases as the gradient index n increases. Some deflected shapes of transverse displacement are depicted in Fig. 3 for the various boundary conditions and gradient index n = 10 . The effect of the scaling factor on the computational cost (CPUtime) of the proposed method is investigated herein. Five regular sets of 5 × 5,7 × 7, 9 × 9, 11 × 11 and 13 × 13 nodes are used for the analysis of a SSSS FG plate (n = 1). The density of the scattered nodes is characterized by the distance of two neighboring nodes, i.e. Δh . The computational time (CPU-time) required for the direct full-matrix solver is evaluated. All the tasks are performed using a personal computer with a Core-i5 CPU @ 2.5 GHz and 4 GB of RAM. Fig. 4 shows the measured CPU-time versus the nodal spacing Δh for several values of the scaling factors in a logarithmic scale. It is clear that the CPU-time increases as the scaling factors increases. It is easy to understand that the increase in the scaling factor leads to the increase of the number of scattered nodes in the domain of influence, and consequently increases the computational time. As a result, the scaling factor has a strong impact on the computational time. For comparison of the computational times of the conventional FSDT-MK and the new S-FSDT-MK, we have evaluated the time required for generation of the stiffness matrix of the two approaches using the same setting and computer resource. It is concluded that while the FSDT-MK requires 25.26 s to generate the stiffness matrix for a set of 13 × 13 nodes, the SFSDT-MK effectively performs this task in 16.77 s. The effect of the scaling factor α on the accuracy of the obtained normalized deflection is also investigated. Fig. 5 shows that the normalized deflection at center of the FG plate varies with the scaling factor when different sets of nodes are used. In this case, the value of the correlation parameter θ is set to 3.0. Fig. 5 depicts the convergence behavior of the central deflection for various values of α . The results of this figure imply that the central deflection converges to the reference solution obtained by FSDT approach [20] when 2.0 < α < 3.5. For α > 4.0 , the accuracy of the results deteriorates. 5.2. Free vibration analysis A fully simply supported Al / Al2 O3 square thick plate with different length to thickness ratios a / h is analyzed. The results of the first 9

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

Fig. 7. The first six mode shapes of a simply supported Al /Al2 O3 thin plate: (a) mode 1, (b) mode 2, (c) mode 3, (d) mode 4, (e) mode 5, and (f) mode 6.

10

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

Fig. 8. The first six mode shapes of a fully clamped Al /Al2 O3 thin plate: (a) mode 1, (b) mode 2, (c) mode 3, (d) mode 4, (e) mode 5, and (f) mode 6.

a gradient index of n = 1 is taken.

analytical solutions [27]. It is clear that, as the boundary conditions change from SFSF to SSSS, SCSC, CCCC, the values of the natural frequency gradually increase. Furthermore, values of the natural frequency increase as the gradient index decreases. The first six Eigen modes of the thin plate with simply supported and fully clamped boundary conditions are visualized in Figs. 7 and 8 respectively, where

6. Other possible choices of the correlation functions In Eq. (22), the widely used Gaussian function adopted as the correlation function R (xi , xj) has been taken for the analysis. This 11

Engineering Analysis with Boundary Elements 79 (2017) 1–12

T.-V. Vu et al.

fully thank all the anonymous reviewers for their constructive comments and suggestions to improve the quality of the manuscript.

Gaussian function inherently is dependent upon the so-called correlation parameter θ , which however is not related to any physical aspects of the problem. As far as we know there are no analytical expressions or theoretical rules to find out the optimal value of the correlation parameter for all problems. In practice, its reasonable value could be determined through numerical examinations to ensure the accuracy and consistency of the solution based on the properties of each type of problems. In order to get rid of this correlation parameter, other possible correlation functions are suggested to use in future works. As inspired by the non-local continuum damage mechanics concept [28], some nonlocal weight functions, which depend only on the distance between the source point and the target point, can be used. In this way, the Gaussian function is used as correlation function but the correlation parameter θ is taken to be independent of the nodal distribution. Instead, the correlation parameter is replaced by another factor, which relates to the distribution characteristics of the model as follows:

R (xi , xj) = e

⎛ ⎞ −⎜⎜ 1 ⎟⎟ r ij2 ⎝ 2lc2 ⎠

References [1] Bui QT, Do VT, Ton THL, Doan HD, Tanaka S, Pham TD, Nguyen-Van T-A, Yu TT, Hirose S. On the high temperature mechanical behaviors analysis of heated functionally graded plates using FEM and a new third-order shear deformation plate theory. Compos Part B: Eng 2016;92:218–41. [2] Reddy JN. Analysis of functionally graded plates. Int J Numer Methods Eng 2000;684:663–84. [3] Castellazzi G, Gentilini C, Krysl P, Elishakoff I. Static analysis of functionally graded plates using a nodal integrated finite element approach. Compos Struct 2013;103:197–200. [4] Singha MK, Prakash T, Ganapathi M. Finite element analysis of functionally graded plates under transverse load. Finite Elem Anal Des 2011;47(4):453–60. [5] Reddy JN, Chin C. Thermo mechanic alanalysis of functionally graded cylinders and plates. J Therm Stress 1998;212(6):593–626. [6] Liu P, Bui QT, Zhu D, Yu TT, Wang JW, Yin SH, Hirose S. Buckling failure analysis of cracked functionally graded plates by a stabilized discrete shear gap extended 3-node triangular plate element. Compos Part B: Eng 2015;77:179–93. [7] Zhu P, Zhang LW, Liew KM. Geometrically nonlinear thermo mechanical analysis of moderately thick functionally graded plates using a local Petrov–Galerkin approach. Comput Mech 2004;35(4):685–97. [8] Bathe KJ, Brezzi F, Cho SW. The MITC7 and MITC9 plate bending elements. Comput Struct 1989;32(3–4):797–814. [9] Bathe KJ, Dvorkin EN. A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation. Int J Numer Methods Eng 1983;21(2):367–83. [10] Somashekar BR, Prathap G, Babu CRA. Field-consistent four-noded laminated anisotropic plate/shell element. Comput Struct 1987;25(3):345–53. [11] Yu TT, Yin SH, Bui QT, Hirose S. A simple FSDT-based isogeometric analysis for geometrically nonlinear analysis of functionally graded plates. Finite Elem Anal Des 2015;96:1–10. [12] Thai HT, Choi DH. A simple first-order shear deformation theory for the bending and free vibration analysis of functionally graded plates. Comput Struct 2013;101:332–40. [13] Swaminathan K, Naveenkumar DT, Zenkour AM, E. Carrera. Stress, vibration and buckling analyses of FGM plates-A state-of-the-art review. Comput Struct 2015;120:10–31. [14] Chen SS, Xu CJ, Tong GS, Wei X. Free vibration of moderately thick functionally graded plates by a meshless local natural neighbour interpolation method. Eng Anal Bound Elem 2015;61:114–26. [15] Xu CJ, Tong GS. A meshless local natural neighbour interpolation method to modeling of functionally graded viscoelastic materials. Eng Anal Bound Elem 2015;52:92–8. [16] Li QH, Chen SS, Luo XM. Using meshless local natural neighbour interpolation method to solve two-dimensional nonlinear problems. Int J Appl Mech 2016;8(5):1650069. [17] Gu L. Moving kriging interpolation and element free Galerkin method. Int J Num Methods Eng 2003;56:1–11. [18] Arnold D, Falk R. The boundary layer for the Reissner–Mindlin plate model. SIAM J Math Anal 1990;21(2):281–312. [19] Bui QT, Nguyen NT, Nguyen-Dang H. A moving kriging interpolation-based meshless method for numerical simulation of Kirchhoff plate problems. Int J Numer Meth Eng 2009;77:1371–95. [20] Yin SH, Jack SH, Yu TT, Bui QT, Boras PAS. Isogeometric locking-free plate element: a simple first order shear deformation theory for functionally graded plates. Comp Struct 2014;118:121–38. [21] Lee YY, Zhao X, Liew KM. Thermo elastic analysis of functionally graded plates using the element-free kp-Ritz method. Smart Mater Struct 2009;18(3):35007. [22] Nguyen-Xuan H, Tran LV, Nguyen-Thoi T, Vu-Do HC. Analysis of functionally graded plates using an edge-based smoothed finite element method. Compos Struct 2011;93(11):3019–39. [23] Ferreira AJM, Batra RC, Roque CMC, Qian LF, Jorge RMN. Natural frequencies of functionally graded plates by a meshless method. Compos Struct 2006;75:593–600. [24] Qian LF, BatraR C, 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. Composites: Part B 2004;35:685–97. [25] 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 PartB: Eng 2012;43:711–25. [26] 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. PartB: Eng 2013;44:657–74. [27] Baferani AH, Saidi AR, Jomehzadeh E. An exact solution for free vibration of thin functionally graded rectangular plates. In: Proceedings Inst. Mech. Eng. Part C J. Mech. Eng. Sci. Vol. 225(C3); 2011, p. 526–3. [28] Jirasek M. Nonlocal damage mechanics. Rev Eur Genie Civ 2007;11(7–8):993–1021. [29] Bui QT, Khosravifard A, Zhang Ch, Hematiyan MR, Golub MV. Dynamic analysis of sandwich beams with functionally graded core using a truly meshfree radial point interpolation method. Eng Struct 2013;47:90–104. [30] Hematiyan MR, Khosravifard A, Bui QT. Efficient evaluation of weakly/strongly singular domain integrals in the BEM using a singular nodal integration method. Eng Anal Bound Elem 2013;37:691–8. [31] Khosravifard A, Hematiyan MR, Bui QT, Do VT. Accurate and efficient analysis of stationary and propagating crack problems by meshless methods. Theor Appl Fract Mech 2017;87:21–34. http://dx.doi.org/10.1016/j.tafmec.2016.10.004. [32] Racz D, Bui QT. Novel adaptive meshfree integration techniques in meshless methods. Int J Numer Methods Eng 2012;90:1414–34.

(44)

where lc is the internal length factor of the model, which can be taken as the average distance between nodes in the model. Another possible choice of the correlation function for the present method is the truncated quartic polynomial function [28], which is explicitly expressed as follows:

⎧⎛ ⎪ ⎜1 − R (xi , xj) = ⎨ ⎝ ⎪ ⎩0

rij2 ⎞

2

⎟ lc2 ⎠

if 0 ≤ rij < lc if lc ≤ rij

(45)

It is obvious that the new correlation functions do not depend on the correlation parameter. 7. Conclusions and outlook In this paper, a simple yet effective first order shear deformation plate theory integrated with meshfree moving Kriging method has been developed for the analysis of free vibration and static deflection of functionally graded ceramic-metal plates. By using the simplifying assumptions to the traditional first order shear deformation plate theory, the number of unknown variables and governing equations are reduced by one, hence making it simpler and more effective to implement. The present approach can be applied to thick as well as thin plates without using special techniques or treatments, e.g., shear-locking effect. The essential boundary conditions are enforced directly in the same way as the conventional finite element method. The numerical results of deflections and natural frequencies calculated by the present S-FSDT based meshfree formulation are compared with analytical, finite element, other meshless, and isogeometric methods, and very good agreements are obtained. The developed S-FSDT meshless method can be considered as an alternative approach for modeling mechanical behavior of plates. Additionally, we have discussed various possible choices of the correlation functions, free from the correlation parameter. The present formulation could be further enhanced by taking into account the effective numerical integrations, which have recently developed by the authors, e.g., either the Cartesian transformation method (CTM) [29–31] or the adaptive one [32]. The motivation of employing those new numerical integrations is to avoid the background cells, consequently forming truly meshfree methods. Acknowledgements The financial support by Ho Chi Minh City University of Architecture, Vietnam is gratefully acknowledged. The authors grate12