A Strong Form Meshless Method for the Solution of FGM Plates

Abstract

Laminated composite structures suffer from failure because of concentrations of gradient fields on interfaces due to discontinuity of material properties. The rapid development of material science enables designers to replace classical laminated plate elements in aerospace structures with more advanced ones made of functionally graded materials (FGM), which are microscopic composite materials with continuous variation of material coefficients according to the contents of their micro-constituents. Utilization of FGM eliminates the inconvenience of laminated structures but gives rise to substantial changes in structural design This paper deals with the presentation of a strong formulation meshless method for the solution of FGM composite plates. Recall that the fourth-order derivatives of deflections are involved in the governing equations for plate structures. However, the high-order derivatives of field variables in partial differential equations (PDE) lead to increasing inaccuracy of approximations. For that reason, the decomposition of the high-order governing equations into the second-order PDE is proposed. For the spatial approximation of field variables, the meshless moving least square (MLS) approximation technique is employed. The reliability (numerical stability, convergence, and accuracy) as well as computational efficiency of the developed method is illustrated by several numerical investigations of the response of FGM plates with the transversal gradation of material coefficients under stationary and/or transient mechanical and thermal loadings.

1. Introduction

With the development of advanced composite materials, their application within aerospace engineering has become more widespread, and laminated composite plates are used as fundamental structural elements [1,2,3,4,5,6,7,8]. However, the discontinuities of material coefficients on the interfaces between two or more layers with different material properties can lead to concentrations of gradient fields and also to catastrophic failure of aerospace structures [9,10,11]. To overcome the problem of the delamination of laminated composite plates the functionally graded materials (FGM) are looking like promising alternatives of laminated structures. Due to the great potential in several branches of engineering applications, the FGMs are in the focus of researchers, which can be confirmed by a huge number of papers devoted to analyses with using various plate bending theories for FGM plates under static and/or dynamic loadings [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]. FGM structure can be characterized as a composite fabricated by mixing two or more different material micro-constituents [28,29]. In continuum models, the effective material properties, such as elastic modulus, Poisson’s ratio, mass density, etc., of the functionally graded composites are determined by the volume fraction distribution of the dispersed phases. In the literature several approaches are available for FGM modeling, such as the Mori–Tanaka scheme [30,31], composite cylindrical assemblage model [32,33], the simplified strength of materials method [34,35], and one of the most frequently used approach the rule of mixture where the material coefficients of multiphase materials are related directly to the volume fractions and individual coefficients of the constituents [36,37].

Despite FGMs being microscopically highly non-homogeneous structures, they exhibit continuous properties changing smoothly with respect to the spatial coordinates in macro-structural elements. The continuum models of FGM have been widely used in many engineering applications, for instance, the static behavior of FG rectangular plates was studied by Reddy [16], and also an axisymmetric formulation for circular and annular FG plate bending was observed [38] within the third-order shear deformation plate theory. Xu and Zhou [39] studied rectangular FG plates with variable plate thickness and exponential gradation of elastic modulus in the lateral direction as a 3D problem.

In the last three decades, the meshless formulations have become an attractive alternative to classical mesh-based discretization methods (such as FEM, BEM, etc.), because of various advantages [40] resulting from the fact that only nodes are used for approximation instead of elements and no elements are needed for background integration like in the element-free Galerkin (EFG) method [41,42]. Additional advantages of the meshless formulations arise in numerical solutions of boundary value problems for partial differential equations (PDE) with variable coefficients like in the case of functionally graded materials (FGM), since the complexity is not increased as compared with the case of homogeneous media. The maximum order of spatial derivatives in weak formulations is lower by one than the order of governing PDE.

The main drawback of meshless methods is the worse computational efficiency than that in mesh-based methods because the shape functions are not expressed as elementary function. This handicap is solved in the present paper by decreasing the number of evaluations of shape functions by using the strong formulation, which can be shown to be a limit case of local weak formulations with approaching the size of the local subdomain to zero. The decrease in the accuracy of approximation of high-order derivatives in plate bending problems is overcome by decomposition of governing equations into governing equations with lower-order derivatives. This is performed at the price of an increase in nodal unknowns. The moving least square (MLS) approximation scheme [43] with using the so-called central approximation node (CAN) concept [44,45] is employed for approximation of field variables. In regards to the numerical evaluation of the derivatives of field variables, we shall use the standard as well as the modified differentiation [46,47] with using not higher than first-order derivatives of shape functions even for higher-order derivatives of field variables. Numerical tests for accuracy and computational efficiency of these approaches are presented in this paper. Moreover, the modified evaluation of the shape functions and their derivatives will be employed [46,47] in order to obey correctly the requirement of completeness of the used set of shape functions. The accuracy and convergence study is accomplished on boundary value problems for bending of thin rectangular FGM plate, where the exact solutions are available and can be used as the benchmark solutions. The attention is paid also to numerical parametric study with respect to various parameters of gradations of the material coefficients of the FGM plates. Several numerical examples are presented to investigate the static and transient response of elastic FGM plates, and thermomechanical response of the FGM plates with transversal and/or in-plane gradation of material properties under stationary and/or transient mechanical and thermal loadings. The main goal of this paper is to present a well-developed computational method (with attributes mentioned in the previous paragraph) for a broad class of plate bending problems. The functional dependence of material coefficients (such as Young modulus, mass density, specific heat capacity, thermal expansion, heat conduction coefficient) gives rise to coupling effects among the field variables including the interaction between the in-plane deformation and bending modes. The unified formulation covers three different kinds of deformation assumptions known as the classical Kirchhoff–Love theory for thin plate bending, the first-order shear deformation theory and the third-order shear deformation theory. The governing equations and possible boundary conditions are developed from the variation principles. In Section 2.1, the formulation is derived for FGM elastic plates under static and/or dynamic loadings. The coupled thermoelastic problems are considered in Section 2.2. Section 3 is devoted to numerical implementation of meshless discretization with using strong formulation and the MLS approximation for spatial dependence of field variables. Finally, in Section 4, we present illustrative examples in order to demonstrate the reliability of the computational method, its universality and discuss the coupling phenomena within three various plate bending theories.

2. Governing Equations

2.1. Governing Equation for Elastic FGM Plates

In this chapter, we shall derive the unified formulation for all three plate bending theories, such as the Kirchhoff–Love theory (KLT), the first-order shear deformation plate theory (FSDPT), and the third-order shear deformation plate theory (TSDPT), in which the displacement field $v_i(x,x_{3,}t)$ can be expressed in terms of the in-plane displacements $u_{\alpha }(x,t)$, transversal deflections $w(\mathbf{x},t)$ and rotations of the normal to the mid-surface ${\phi }_{\alpha }(x,t)$ by

$$v_i(x,x_3,t)={\delta }_{i\alpha }\left\{u_{\alpha }(x,t)+\left[c_1\varphi (x_3)-x_3\right]w_{,\alpha }(x,t)+c_1\varphi (x_3){\phi }_{\alpha }(x,t)\right\}+{\delta }_{i3}w(x,t)$$

(1)

where $i=1,2,3$, $\alpha =1,2$, $\varphi (x_3):=x_3-c_2\psi (x_3)$, $\psi (x_3):=\frac{4}{3}{\left(\frac{x_3}{h}\right)}^2{x}_3$, $x_3\in \left[-h/2,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}h/2\right]$, $x\in \Omega$, with $\Omega$ being the plate area in the mid-surface (Figure 1).

The key factors $c_1$ and $c_2$ in Equation (1) were introduced for the proper selection among various plate bending theories within the unified formulation, and are specified as:

(i)

for KLT: $c_1=0$, $c_2=0$,

(ii)

for FSDPT: $c_1=1$, $c_2=0$,

(iii)

for TSDPT: $c_1=1$, $c_2=1$.

According to the Hooke law, the 3D stresses in the linear elastic and isotropic continuum are given as

$${\sigma }_{ij}(x,x_3,t)=\frac{E}{1-{\nu }^2}\frac{1-\nu }{H}\left[H{e}_{ij}(x,x_3,t)+\nu {\delta }_{ij}{e}_{kk}(x,x_3,t)\right]$$

(2)

where $E,\hspace{0.17em}\nu$ stand for the elastic modulus and Poisson’s ratio, respectively. The total strains are expressed as $e_{ij}=\left(v_{i,j}+v_{j,i}\right)/2$, while

$$H=1-\chi \nu , \chi =\left\{\begin{array}{c}1,\hspace{0.17em}\hspace{0.17em}\mathrm{for} plane\hspace{0.17em}\hspace{0.17em}stress\hspace{0.17em}\hspace{0.17em}problems\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ 2,\hspace{0.17em}\hspace{0.17em}for 3D\hspace{0.17em}\hspace{0.17em}problems or\hspace{0.33em}plane\hspace{0.17em}\hspace{0.17em}strain\hspace{0.17em}\end{array}\right.$$

(3)

By substitution of Equation (1) into Equation (2), the stresses inside an elastic plate considered within three plate bending theories are obtained

$${\sigma }_{\alpha \beta }(x,x_3,t)=\frac{E}{1-{\nu }^2}\frac{1-\nu }{H}\left[{\tau }_{\alpha \beta }^{(u)}(x,t)+c_1\varphi (x_3){\tau }_{\alpha \beta }^{(\phi )}(x,t)+\left(c_1\varphi (x_3)-x_3\right){\tau }_{\alpha \beta }^{(w)}(x,t)\right]$$

$${\sigma }_{\alpha 3}(x,x_3,t)=\frac{E}{1+\nu }\frac{c_1}{2}{\varphi }^{\prime }(x_3)\left[w_{,\alpha }(x,t)+{\phi }_{\alpha }(x,t)\right], {\varphi }^{\prime }(x_3)={\partial }_3\varphi (x_3)$$

(4)

$${\sigma }_{33}(x,x_3,t)=\frac{E\nu }{1-{\nu }^2}\frac{1-\nu }{H}\left[u_{\gamma ,\gamma }(x,t)+c_1\varphi (x_3){\phi }_{\gamma ,\gamma }(x,t)+\left(c_1\varphi (x_3)-x_3\right)w_{,\gamma \gamma }(x,t)\right]$$

in which the strain contributions associated with in-plane displacements, rotations, and transversal displacements (deflections) are

$${\tau }_{\alpha \beta }^{(u)}(x,t):=H{\epsilon }_{\alpha \beta }(x,t)+\nu {\delta }_{\alpha \beta }{\epsilon }_{\gamma \gamma }(x,t),{\tau }_{\alpha \beta }^{(\phi )}(x,t):=H{\eta }_{\alpha \beta }(x,t)+\nu {\delta }_{\alpha \beta }{\eta }_{\gamma \gamma }(x,t)$$

$${\tau }_{\alpha \beta }^{(w)}(x,t):=H{w}_{,\alpha \beta }(x,t)+\nu {\delta }_{\alpha \beta }{w}_{,\gamma \gamma }(x,t)$$

(5)

with ${\epsilon }_{\alpha \beta }=\left(u_{\alpha ,\beta }+u_{\beta ,\alpha }\right)/2$,${\eta }_{\alpha \beta }=\left({\phi }_{\alpha ,\beta }+{\phi }_{\beta ,\alpha }\right)/2$.

Within this paper, we utilize the Einstein summation rule with respect to repeated subscripts, and the comma preceding a subscript denotes the derivative with respect to the corresponding coordinate. Thus,

$$w{(\mathbf{x})}_{,\gamma \gamma }=w_{,11}(\mathbf{x})+w_{,22}(\mathbf{x})=\frac{{\partial }^{2}w(\mathbf{x})}{\partial x_1\partial x_1}+\frac{{\partial }^{2}w(\mathbf{x})}{\partial x_2\partial x_2}$$

(6)

Making use of the rule of mixture for micro-composites with two constituents, one can obtain material properties for FGM plates with power law gradation in the transversal direction as

$$E(x,x_3)=E_0{E}_H(x)E_V(x_3), E_V(x_3)=1+\zeta {\left(\frac{1}{2}\pm \frac{x_3}{h}\right)}^p$$

$$\rho (x,x_3)={\rho }_0{\rho }_H(x){\rho }_V(x_3), {\rho }_V(x_3)=1+\epsilon {\left(\frac{1}{2}\pm \frac{x_3}{h}\right)}^g$$

(7)

where the Young modulus $E$, and mass density $\rho$ are factorized using the subscript 0, H, V for the reference values, factors responsible for in-plane (horizontal) gradation and transversal gradation, respectively. The levels of gradation are determined by $\zeta$ and $\epsilon$, while $p$ and $g$ are exponents of the power law gradations. In addition, the thickness of the FGM plate is allowed to be variable, $h(x)$.

For the derivation of the equations of motion together with the boundary and initial conditions for mechanical fields, Hamilton’s principle is utilized

$$\delta U-\delta W_e-\delta K=0$$

(8)

where the variations of the elastic deformation energy, the kinetic energy, and the energy of external transversal forces are expressed as

$$\delta U={\displaystyle \underset{0}{\overset{T}{\int }}{\displaystyle \underset{\Omega }{\int }\left({\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{ij}(x,x_3,t)\delta e_{ij}(x,x_3,t)d{x}_3}\right)d\Omega }dt}$$

$$\delta K={\displaystyle \underset{0}{\overset{T}{\int }}{\displaystyle \underset{\Omega }{\int }\left({\displaystyle \underset{-h/2}{\overset{h/2}{\int }}\rho {\dot{v}}_i\delta {\dot{v}}_{i}d{x}_3}\right)d\Omega }}\hspace{0.17em}dt$$

(9)

$$\delta W_e={\displaystyle \underset{0}{\overset{T}{\int }}{\displaystyle \underset{\Omega }{\int }{t}_3(x,t)\delta w}(x,t)d\Omega dt}$$

in which T is the considered time period.

The time derivative of field variables is marked by a dot over the variable. Performing the integrations in a transversal direction, we obtain

$$\begin{array}{l}\delta U={\displaystyle \underset{0}{\overset{T}{\int }}{\displaystyle \underset{\Omega }{\int }\left\{T_{\alpha \beta }(x,t)\delta u_{\alpha ,\beta }(x,t)-M_{\alpha \beta }^{(w)}(x,t)\delta w_{,\alpha \beta }(x,t)+M_{\alpha \beta }^{(\phi )}(x,t)\delta {\phi }_{\alpha ,\beta }(x,t)+\right.\hspace{0.17em}}}\\ \left.+T_{3\beta }^{(w\phi )}(x,t)\left[\delta w_{,\beta }(x,t)+\delta {\phi }_{\beta }(x,t)\right]\right\}d\Omega dt\end{array}$$

(10)

$$\begin{array}{cc}\hfill & \delta K=\underset{0}{\overset{T}{\int }}\underset{\Omega }{\int }\left\{\left[{\tilde{I}}^{\left(uu\right)}\left(\mathbf{x}\right){\ddot{u}}_{\beta }(\mathbf{x},t)+{\tilde{I}}^{\left(u\phi \right)}\left(\mathbf{x}\right){\ddot{\phi }}_{\beta }(\mathbf{x},t)+{\tilde{I}}^{\left(uw\right)}\left(\mathbf{x}\right){\ddot{w}}_{,\beta }(\mathbf{x},t)\right]\delta u_{\beta }(\mathbf{x},t)+\right.\hfill \\ \hfill & +\left[{\tilde{I}}^{\left(\phi u\right)}\left(\mathbf{x}\right){\ddot{u}}_{\beta }(\mathbf{x},t)+{\tilde{I}}^{\left(\phi \phi \right)}\left(\mathbf{x}\right){\ddot{\phi }}_{\beta }(\mathbf{x},t)+{\tilde{I}}^{\left(\phi w\right)}\left(\mathbf{x}\right){\ddot{w}}_{,\beta }(\mathbf{x},t)\right]\delta {\phi }_{\beta }(\mathbf{x},t)+\hfill \\ \hfill & +\left[{\tilde{I}}^{\left(wu\right)}\left(\mathbf{x}\right){\ddot{u}}_{\beta }(\mathbf{x},t)+{\tilde{I}}^{\left(w\phi \right)}\left(\mathbf{x}\right){\ddot{\phi }}_{\beta }(\mathbf{x},t)+{\tilde{I}}_2^{\left(ww\right)}\left(\mathbf{x}\right){\ddot{w}}_{,\beta }(\mathbf{x},t)\right]\delta w_{,\beta }(\mathbf{x},t)+\hfill \\ \hfill & \left.+{\tilde{I}}_1^{\left(ww\right)}\left(\mathbf{x}\right)\ddot{w}(\mathbf{x},t)\delta w(\mathbf{x},t)\right\}d\Omega dt\hfill \end{array}$$

(11)

where the following integral quantities for in-plane stresses

$$T_{\alpha \beta }(x,t):={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{\alpha \beta }(x,x_3,t)d{x}_3}$$

(12)

stress couples

$$M_{\alpha \beta }^{(\phi )}(x,t):=c_1{\displaystyle \underset{-h/2}{\overset{h/2}{\int }}\varphi (x_3){\sigma }_{\alpha \beta }(x,x_3,t)d{x}_3}, M_{\alpha \beta }^{(w)}(x,t):={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{x}_3{\sigma }_{\alpha \beta }(x,x_3,t)d{x}_3}-M_{\alpha \beta }^{(\phi )}(x,t)$$

(13)

and shear stresses

$$T_{3\beta }^{(w\phi )}(x,t):=c_1{\displaystyle \underset{-h/2}{\overset{h/2}{\int }}\left\{\left[\left(1-c_2\right)\kappa +c_2\right]-c_2{\psi }^{\prime }(x_3)\right\}{\sigma }_{3\beta }(x,x_3,t)d{x}_3}$$

(14)

represent the considered fields averaged through the plate thickness.

The mass inertia terms are defined as

$$\left\{{\tilde{I}}^{\left(uu\right)},{\tilde{I}}^{\left(u\phi \right)},{\tilde{I}}^{\left(uw\right)}\right\}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}\left\{1,\hspace{0.33em}{c}_1\varphi (x_3),\hspace{0.33em}\left(c_1\varphi (x_3)-x_3\right)\right\}}\rho (x_3)d{x}_3$$

$$\left\{{\tilde{I}}^{\left(\phi u\right)},{\tilde{I}}^{\left(\phi \phi \right)},{\tilde{I}}^{\left(\phi w\right)}\right\}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}\left\{c_1\varphi (x_3),\hspace{0.33em}{c}_1{\varphi }^2(x_3),\hspace{0.33em}{c}_1\left({\varphi }^2(x_3)-x_3\varphi (x_3)\right)\right\}}\rho (x_3)d{x}_3$$

(15)

$$\begin{array}{l}\left\{{\tilde{I}}^{\left(wu\right)},{\tilde{I}}^{\left(w\phi \right)},{\tilde{I}}_1^{\left(ww\right)},{\tilde{I}}_2^{\left(ww\right)}\right\}=\\ \\ ={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}\left\{\left(c_1\varphi (x_3)-x_3\right),\hspace{0.33em}{c}_1\left({\varphi }^2(x_3)-x_3\varphi (x_3)\right),\hspace{0.33em}1,\hspace{0.33em}{\left(c_1\varphi (x_3)-x_3\right)}^2\right\}}\rho (x_3)d{x}_3\end{array}$$

The shear correction factor, $\kappa$, represents the Reissner modification of the shear stresses in the case of the FSDPT. During the derivation of Equation (11), we were keeping in mind that the variations $\delta u_{\beta }$, $\delta {\phi }_{\beta }$, $\delta w$ are vanishing at $t=0$ and $t=T$, hence $\delta w_{,\beta }(x,t=0)=0\hspace{0.33em}$ and $\delta w_{,\beta }(x,t=T)=0$ as well. The explicit expressions for the averaged fields defined by Equations (12)–(14) and mass inertia terms defined by Equation (15) can be found in [48].

Certain rearrangements [48,49] in $\delta U-\delta W$ can lead to expression without derivatives of variations of field variables. Then, the set of governing equations at $x\in \Omega$, $t\in [0,\hspace{0.17em}T]$ and restrictions on the boundary edge $\partial \Omega$ at $t\in [0,\hspace{0.17em}T]$ are obtained:

governing equations

$$T_{\alpha \beta ,\beta }(x,t)+{\tilde{I}}^{(uu)}(x){\ddot{u}}_{\alpha }(x,t)+{\tilde{I}}^{(u\phi )}(x){\ddot{\phi }}_{\alpha }(x,t)+{\tilde{I}}^{(uw)}(x){\ddot{w}}_{,\alpha }(x,t)=0$$

$$\begin{array}{l}{M}_{\alpha \beta ,\alpha \beta }^{(w)}(x,t)+T_{3\beta ,\beta }^{(w\phi )}(x,t)+\\ +{\tilde{I}}_1^{(ww)}(x)\ddot{w}(x,t)-{\left[{\tilde{I}}^{(wu)}(x){\ddot{u}}_{\beta }(x,t)+{\tilde{I}}^{(w\phi )}(x){\ddot{\phi }}_{\beta }(x,t)+{\tilde{I}}_2^{(ww)}(x){\ddot{w}}_{,\beta }(x,t)\right]}_{,\beta }=-t_3(x,t)\end{array}$$

(16)

$$M_{\alpha \beta ,\beta }^{(\phi )}(x,t)-T_{3\alpha }^{(w\phi )}(x,t)+{\tilde{I}}^{(\phi u)}(x){\ddot{u}}_{\alpha }(x,t)+{\tilde{I}}^{(\phi \phi )}(x){\ddot{\phi }}_{\alpha }(x,t)+{\tilde{I}}^{(\phi w)}(x){\ddot{w}}_{,\alpha }(x,t)=0$$

boundary restrictions

$$n_{\beta }(x)T_{\alpha \beta }(x,t)\delta u_{\alpha }(x,t)=0$$

$$n_{\alpha }(x)n_{\beta }(x)M_{\alpha \beta }^{(w)}(x,t)\delta \left(\frac{\partial w}{\partial n}(x,t)\right)=0$$

(17)

$$n_{\beta }(x)M_{\alpha \beta }^{(\phi )}(x,t)\delta {\phi }_{\alpha }(x,t)=0$$

$$\left\{V(x,t)-n_{\beta }(x)\left[{\tilde{I}}^{(wu)}(x){\ddot{u}}_{\beta }(x,t)+{\tilde{I}}^{(w\phi )}(x){\ddot{\phi }}_{\beta }(x,t)+{\tilde{I}}_2^{(ww)}(x){\ddot{w}}_{,\beta }(x,t)\right]\right\}\delta w(x,t)=0$$

where

$$V(x,t):=n_{\alpha }(x)\left(M_{\alpha \beta ,\beta }^{(w)}(x,t)+T_{3\alpha }^{(w\phi )}(x,t)\right)+\frac{\partial }{\partial t}{T}^{(w)}(x,t)-{\displaystyle \sum _c\delta (x-x^c)\left[\left[T^{(w)}(x^c,t)\right]\right]}$$

is the generalized shear force, $T^{(w)}(x,t):=t_{\alpha }(x)n_{\beta }(x)M_{\alpha \beta }^{(w)}(x,t)$ is the twisting moment defined on the plate edge $\partial \Omega$. Furthermore, the jump at a corner on the oriented boundary edge $\partial \Omega$ is defined as

$$\left[\left[A\left({\mathbf{x}}^c\right)\right]\right]:=A(x^c-0)-A(x^c+0)$$

It is appropriate to introduce the dimensionless fields

$$u_{\beta }^{\ast }(x,t):=\frac{u_{\beta }(x,t)}{h_0}, {\phi }_{\beta }^{\ast }(x,t):={\phi }_{\beta }(x,t), w^{\ast }(x,t):=\frac{w(x,t)}{h_0}$$

(18)

in which the dimensionless variables are defined as

$$x_{\beta }^{\ast }:=\frac{x_{\beta }}{L}, x_3^{\ast }:=\frac{x_3}{h_0}=h^{\ast }(x)z, h(x)=h_0{h}^{\ast }(x), t^{\ast }:=\frac{t}{T}$$

(19)

where $L$ is the characteristic length for the in-plane dimensions and $h_0$ is the thickness of the plate.

In what follows, for simplicity, we shall write $g^{\ast }(x,t)$ instead of $g^{\ast }(x^{\ast },t^{\ast })$.

Substituting Equations (18) and (19) into Equation (16), one can obtain the governing equation in sense of dimensionless formulation

$$T_{\alpha \beta ,\beta }^{\ast }(x,t)+I^{(uu)}(x){\ddot{u}}_{\alpha }^{\ast }(x,t)+I^{(u\phi )}(x){\ddot{\phi }}_{\alpha }^{\ast }(x,t)+I^{(uw)}(x){\ddot{w}}_{,\alpha }^{\ast }(x,t)=0$$

$$\begin{array}{l}{M}_{\alpha \beta ,\alpha \beta }^{\ast (w)}(x,t)+{\left(\frac{L}{h_0}\right)}^2{T}_{3\beta ,\beta }^{\ast (w\phi )}(x,t)+I_1^{(ww)}(x){\ddot{w}}^{\ast }(x,t)-\\ -{\left[I^{(wu)}(x){\ddot{u}}_{\beta }^{\ast }(x,t)+I^{(w\phi )}(x){\ddot{\phi }}_{\beta }^{\ast }(x,t)+I_2^{(ww)}(x){\ddot{w}}_{,\beta }^{\ast }(x,t)\right]}_{,\beta }=-t_3^{\ast }(x,t)\end{array}$$

(20)

$$M_{\alpha \beta ,\beta }^{\ast (\phi )}(x,t)-{\left(\frac{L}{h_0}\right)}^2{T}_{3\alpha }^{\ast (w\phi )}(x,t)+I^{(\phi u)}(x){\ddot{u}}_{\alpha }^{\ast }(x,t)+I^{(\phi \phi )}(x){\ddot{\phi }}_{\alpha }^{\ast }(x,t)+I^{(\phi w)}(x){\ddot{w}}_{,\alpha }^{\ast }(x,t)=0$$

The dimensionless transversal mechanical loading is defined as $t_3^{\ast }(x,t):=\left(L^4/D_0{h}_0\right)t_3(x,t)$, where $D_0=\frac{E_0{(h_0)}^3}{12(1-{\nu }^2)}$ is the bending stiffness. For elastodynamic plate bending problems, the explicit expressions for the semi-integral fields (12)–(14) and the mass inertia coefficients (15) are given in [48]. The boundary restrictions are obtained from (17) by replacing all fields and mass inertia coefficients by their dimensionless counterparts.

In Equation (20) the fourth-order derivatives of transversal deflections and third-order derivatives of in-plane displacements and rotations are included, which can lead to a reduction in the accuracy of approximations of high-order derivatives of field variables. To overcome this shortcoming, we employ the decomposition technique on the derived set of the partial differential equations (PDE) by introducing additional field variables

$$m^{\ast }(x,t):={\nabla }^2{w}^{\ast }(x,t), s_{\alpha }^{\ast }(x,t):={\nabla }^2{u}_{\alpha }^{\ast }(x,t), f_{\alpha }^{\ast }(x,t):={\nabla }^2{\phi }_{\alpha }^{\ast }(x,t)$$

(21)

Then, the final set of governing equations for primary field variables $\left\{w^{\ast }(x,t),\hspace{0.17em}{m}^{\ast }(x,t),\hspace{0.17em}{u}_{\beta }^{\ast }(x,t),\hspace{0.17em}{s}_{\beta }^{\ast }(x,t),\hspace{0.17em}{\phi }_{\beta }^{\ast }(x,t),\hspace{0.17em}{f}_{\beta }^{\ast }(x,t)\right\}$ is given by Equations (20) and (21), which include derivatives not higher than second-order [47].

2.2. Governing Equation for Thermoelastic FGM Plates

Similarly, like in the previous chapter, we shall derive the governing equation for the FGM plates under thermal loadings. Within the linear thermoelastic and isotropic media, the 3D stresses are given as

$${\sigma }_{ij}(x,x_3,t)=\frac{E}{1-{\nu }^2}\frac{1-\nu }{H}\left[H{e}_{ij}(x,x_3,t)+\nu {\delta }_{ij}{e}_{kk}(x,x_3,t)-\alpha (1+\nu )\left(\theta (x,x_3,t)-{\theta }_0\right){\delta }_{ij}\right]$$

(22)

where $E,\hspace{0.17em}\nu$, $\alpha$ denote the elastic modulus, Poisson’s ratio, and the linear thermal expansion coefficient, respectively.

Furthermore, $\theta$, ${\theta }_0$ and $e_{ij}$ stand for the temperature, reference temperature, and total strains, respectively.

By substitution of Equation (1) into Equation (22), the stresses inside an elastic plate considered within three plate bending theories are obtained as

$${\sigma }_{\alpha \beta }(x,x_3,t)=\frac{E}{1-{\nu }^2}\frac{1-\nu }{H}\left[{\tau }_{\alpha \beta }^{(u)}(x,t)+c_1\varphi (x_3){\tau }_{\alpha \beta }^{(\phi )}(x,t)+\left(c_1\varphi (x_3)-x_3\right){\tau }_{\alpha \beta }^{(w)}(x,t)-\alpha {\delta }_{\alpha \beta }{\tau }^{(\theta )}(x,x_3,t)\right]$$

$${\sigma }_{\alpha 3}(x,x_3,t)=\frac{E}{1+\nu }\frac{c_1}{2}{\varphi }^{\prime }(x_3)\left[w_{,\alpha }(x,t)+{\phi }_{\alpha }(x,t)\right]$$

(23)

$${\sigma }_{33}(x,x_3,t)=\frac{E\nu }{1-{\nu }^2}\frac{1-\nu }{H}\left[u_{\gamma ,\gamma }(x,t)+c_1\varphi (x_3){\phi }_{\gamma ,\gamma }(x,t)+\left(c_1\varphi (x_3)-x_3\right)w_{,\gamma \gamma }(x,t)-\frac{\alpha }{\nu }{\tau }^{(\theta )}(x,x_3,t)\right]$$

in which the strain contributions associated with in-plane displacements, rotations, and transversal deflections ${\tau }_{\alpha \beta }^{(u)}(x,t)$, ${\tau }_{\alpha \beta }^{(\phi )}(x,t)$, ${\tau }_{\alpha \beta }^{(w)}(x,t)$ are defined by Equation (5) and the temperature contribution is

$${\tau }^{(\theta )}(x,x_3,t):=(1+\nu )\left(\theta (x,x_3,t)-{\theta }_0\right)$$

(24)

The power law gradation in the transversal direction for material coefficients such as elastic modulus, thermal expansion coefficient, heat conduction coefficients, specific heat capacity, and mass density is defined by the following relations:

$$E(x,x_3)=E_0{E}_H(x)E_V(x_3), E_V(x_3)=1+\zeta {\left(\frac{1}{2}\pm \frac{x_3}{h}\right)}^p$$

$$\alpha (x,x_3)={\alpha }_0{\alpha }_H(x){\alpha }_V(x_3), {\alpha }_V(x_3)=1+\xi {\left(\frac{1}{2}\pm \frac{x_3}{h}\right)}^r$$

$$k(x,x_3)=k_0{k}_H(x)k_V(x_3), k_V(x_3)=1+\omega {\left(\frac{1}{2}\pm \frac{x_3}{h}\right)}^s,$$

$$c(x,x_3)=c_0{c}_H(x)c_V(x_3), c_V(x_3)=1+\chi {\left(\frac{1}{2}\pm \frac{x_3}{h}\right)}^q,$$

$$\rho (x,x_3)={\rho }_0{\rho }_H(x){\rho }_V(x_3), {\rho }_V(x_3)=1+\epsilon {\left(\frac{1}{2}\pm \frac{x_3}{h}\right)}^g$$

(25)

The parameters $\zeta ,\hspace{0.17em}\xi ,\omega ,\chi \epsilon$ stand for the levels of power law gradations in transversal direction, while $p,r,s,q,g$ are exponents of these gradations.

In order to develop the formulation consistent with elastic deformation assumptions (with the known $x_3$-dependence), we approximate the temperature field by the second-order power series expansion with respect to the $x_3$-coordinate as

$$\theta (x,x_3,t)\approx {\theta }_0+{\vartheta }_0(x,t)+z{\vartheta }_1(x,t)+z^2{\vartheta }_2(x,t), z=\frac{x_3}{h}\in [-0.5,\hspace{0.17em}\hspace{0.17em}0.5]$$

(26)

where ${\vartheta }_s(x,t)$ are a new field variables with $s=0,1,2$. Then, one can write

$${\tau }^{(\theta )}(x,x_3,t)={\displaystyle \sum _{s=0}^2{z}^s{\tau }^{({\vartheta }_s)}(x,t)}, {\tau }^{({\vartheta }_s)}(x,t):=(1+\nu ){\vartheta }_s(x,t)$$

(27)

The governing equations and the boundary restrictions for mechanical fields are formally given by Equations (16) and (17), respectively. Recall that in case of thermoelasticity, the semi-integral fields (12)–(14) include also thermal contributions [49] in contrast to case of elastodynamics. In the governing Equation (16), we have only five PDEs for eight independent field variables, which means we need three more equations to complete the set of governing equations for FGM plate bending problems in classical thermoelasticity.

The heat conduction equation within coupled classical thermoelasticity is given by

$${\left(k{\theta }_{,j}\right)}_{,j}-\rho c\dot{\theta }-\frac{E\alpha }{H}{\theta }_0{\dot{v}}_{j,j}=0$$

(28)

in which $k$, $\rho$, and $c$ is the heat conduction coefficient, mass density, and the specific heat capacity, respectively.

In view of the previously introduced three field variables ${\vartheta }_s(x,t)$, the heat conduction equation is still dependent on $x_3$ as

$$\begin{array}{l}{\left[k_H(x)k_V(x_3){\left({\displaystyle \sum _{s=0}^2{z}^s{\vartheta }_s(x,t)}\right)}_{,j}\right]}_{,j}-\frac{{\rho }_0{c}_0}{k_0}{\rho }_H(x){\rho }_V(x_3)c_H(x)c_V(x_3)\left[{\displaystyle \sum _{s=0}^2{z}^s{\dot{\vartheta }}_s(x,t)}\right]-\\ -\frac{E_0{\alpha }_0}{k_{0}H}{E}_H(x)E_V(x_3){\alpha }_H(x){\alpha }_V(x_3){\theta }_0{\dot{v}}_{\beta ,\beta }(x,x_3,t)=0\end{array}$$

(29)

In order to accomplish the 2D formulation, we propose to integrate Equation (29) through the thickness of the plate with getting the heat conduction equation in an averaged sense

$$\begin{array}{l}{\displaystyle \sum _{s=1}^2{C}^{(\theta {\vartheta }_s)}(x){\vartheta }_s^{\ast }(x,t)}+{\displaystyle \sum _{s=0}^2{G}_{\beta }^{(\theta {\vartheta }_s)}(x){\vartheta }_{s,\beta }^{\ast }(x,t)}+{\displaystyle \sum _{s=0}^2{G}^{(\theta {\vartheta }_s)}(x){\vartheta }_{s,\beta \beta }^{\ast }(x,t)}+\\ +{\displaystyle \sum _{s=0}^2{I}^{(\theta {\vartheta }_s)}(x){\dot{\vartheta }}_s^{\ast }(x,t)}+C^{(\theta u)}(x){\dot{u}}_{\beta ,\beta }^{\ast }(x,t)+C^{(\theta w)}(x){\dot{w}}_{,\beta \beta }^{\ast }(x,t)+C^{(\theta \phi )}(x){\dot{\phi }}_{\beta ,\beta }^{\ast }(x,t)=0\end{array}$$

(30)

For the details on coefficients $C^{(\cdot \cdot )}$, $G^{(\cdot \cdot )}$, and $I^{(\cdot \cdot )}$ see [50].

The dimensionless fields

$$u_{\beta }^{\ast }(x,t):=\frac{u_{\beta }(x,t)}{h_0}, {\phi }_{\beta }^{\ast }(x,t):={\phi }_{\beta }(x,t), w^{\ast }(x,t):=\frac{w(x,t)}{h_0}, {\vartheta }_s^{\ast }(x,t):=\frac{{\vartheta }_s(x,t)}{{\theta }_0}$$

(31)

are introduced by using the dimensionless variables as

$$x_{\beta }^{\ast }:=\frac{x_{\beta }}{L},x_3^{\ast }:=\frac{x_3}{h_0}=h^{\ast }(x)z, h(x)=h_0{h}^{\ast }(x),t^{\ast }:=\frac{t}{T}$$

(32)

Furthermore, it is written $g^{\ast }(x,t)$ instead of $g^{\ast }(x^{\ast },t^{\ast })$ for the reason of simplicity.

The other two equations result from the 3D boundary conditions considered on the top and bottom surfaces of the plate. These BCs can be given by an arbitrary combination of the following boundary conditions

(i)

Dirichlet type:

$${\vartheta }_0^{\ast }(x,t)\pm \frac{1}{2}{\vartheta }_1^{\ast }(x,t)+\frac{1}{4}{\vartheta }_2^{\ast }(x,t)=\frac{\overline{\theta }(x,\pm h/2,t)}{{\theta }_0}-1$$

(33)

(ii)

Neumann type:

$$\pm k_0{k}_H(x)k_V(\pm 1/2)\left[{\vartheta }_1^{\ast }(x,t)\pm {\vartheta }_2^{\ast }(x,t)\right]=\overline{q}(x,\pm h/2,t)/{\theta }_0$$

(34)

(iii)

Robin type:

$$A\left[1+{\vartheta }_0^{\ast }(x,t)\pm \frac{1}{2}{\vartheta }_1^{\ast }(x,t)+\frac{1}{4}{\vartheta }_2^{\ast }(x,t)\right]\pm B{k}_0{k}_H(x)k_V(\pm h/2)\left[{\vartheta }_1^{\ast }(x,t)\pm {\vartheta }_2^{\ast }(x,t)\right]=0$$

(35)

The $\overline{q}(x,\pm h/2,t)$ and $\overline{\theta }(x,\pm h/2,t)$ are the prescribed values of the heat flux and temperature on the top and bottom surfaces of the plate, respectively. Taking into account that, $\theta (x,x_3=0,t)={\theta }_0+{\vartheta }_0(x,t)$, the boundary conditions on $\partial \Omega$ can be given as

(i)

Dirichlet type:

$${\left.{\vartheta }_0^{\ast }(x,t)\right|}_{\partial \Omega }=\overline{\theta }(x,x_3=0,t)/{\theta }_0-1$$

(36)

(ii)

Neumann type:

$${\left.-k_0{k}_H(x)k_V(0)n_{\beta }(x){\vartheta }_{0,\beta }^{\ast }(x,t)\right|}_{\partial \Omega }=\left(L/{\theta }_0\right)\overline{q}(x,0,t)$$

(37)

(iii)

Robin type:

$${\left.A{\vartheta }_0^{\ast }(x,t)\right|}_{\partial \Omega }+{\left.B{k}_0{k}_H(x)k_V(0)n_{\beta }(x){\vartheta }_{0,\beta }^{\ast }(x,t)\right|}_{\partial \Omega }=0$$

(38)

with $\overline{\theta }(x,x_3=0,t)$ and $\overline{q}(x,0,t)$ being the prescribed values of the temperature and heat flux on the boundary edge of the plate.

Substituting Equations (31) and (32) into Equation (16), one can obtain the dimensionless form of the governing equations which are formally the same as (20), but the semi-integral fields include also the thermal contributions [49].

Thus, the complete set of governing equations for field variables describing the process of plate bending in classical thermoelasticity is given by Equations (20) and (30) plus two equations properly selected from (33)–(35) according to prescribed thermal boundary conditions on the top and bottom of the plate. The possible boundary conditions result form (17) and (36)–(38).

To overcome the shortcoming of decreasing the accuracy of approximation due to high-order derivatives in governing equations we utilize the decomposition formulation introduced in Section 2.1.

Then, the final set of governing equations for primary field variables $\left\{w^{\ast }(x,t),\hspace{0.17em}{m}^{\ast }(x,t),\hspace{0.17em}{u}_{\beta }^{\ast }(x,t),\hspace{0.17em}{s}_{\beta }^{\ast }(x,t),\hspace{0.17em}{\phi }_{\beta }^{\ast }(x,t),\hspace{0.17em}{f}_{\beta }^{\ast }(x,t),\hspace{0.17em}{\vartheta }_0^{\ast }(x,t),\hspace{0.17em}{\vartheta }_1^{\ast }(x,t),\hspace{0.17em}{\vartheta }_2^{\ast }(x,t)\right\}$ is given by Equations (21), (30) and (35), which include derivatives not higher than the second-order [47].

3. Meshless Approximations of Field Variables by Moving Least Square Approximation

In this paper, we shall employ the central approximation node (CAN) concept of the standard MLS approximation. By considering $x^q$ as the CAN for the approximation at a point $x$ the number of nodes implied into the approximation at $x$ is reduced a-priori from $N$ to $N^q$, where $N$ is the number of all nodes and $N^q$ is the number of nodes supporting the approximation at $x^q$, i.e., the number of nodes in the set ${ℳ}^q={\{\forall {\mathbf{x}}^a;w^a({\mathbf{x}}^q)>0\}}_{a=1}^N$, where $w^a(x)$ is the weight function associated with the node $x^a$ and taken at the field point $x$. In this paper, we employ the Gaussian weights. The MLS-CAN approximation is given as

$$u(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}^{\overline{a}}{\varphi }^{(q,a)}(x)}, \overline{a}=n(q,a)$$

(39)

where $\overline{a}$ is the global number of the $a$-th node from the $N^q$ nodal points, ${\widehat{u}}^a$ is a nodal unknown different from the nodal value $u(x^a)$, and ${\varphi }^a(x)$ is the shape function associated with the nodal point $x^a$. Generally, the nearest node to the field point $x$ is selected as the CAN node.

In this paper, we shall utilize two types for differentiation: (a) the standard approach (D0), and (b) the modified approach ($D1$) for the evaluation of derivatives of field variables.

Within the standard approach, the derivative of $u(x)$ can be approximated by differentiating Equation (39), i.e.,

$$u_{,i}(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}^{\overline{a}}{\varphi }_{,i}^{(q,a)}(x)}, u_{,ij}(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}^{\overline{a}}{\varphi }_{,ij}^{(q,a)}(x)}, u_{,ijk}(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}^{\overline{a}}{\varphi }_{,ijk}^{(q,a)}(x)}$$

(40)

For the details on the shape functions and their derivatives see [47].

In the modified approach, the derivatives are approximated by using the MLS shape functions ${\varphi }^{\overline{a}}(x)$ and nodal values ${\widehat{u}}^{\overline{h}}$ and the first-order derivatives of the MLS shape functions. So,

$$u_{,i}(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}_i^{\overline{a}}{\varphi }^{(q,a)}(x)}, u_{,ij}(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}_{ij}^{\overline{a}}{\varphi }^{(q,a)}(x)}, u_{,ijk}(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}_{ijk}^{\overline{a}}{\varphi }^{(q,a)}(x)}.$$

(41)

From (40₁), we have

$$u_{,i}(x^c)\approx {\displaystyle \sum _{h=1}^{N^c}{\widehat{u}}^{\overline{h}}{\varphi }_{,i}^{(c,h)}(x^c)}={\displaystyle \sum _{h=1}^{N^c}{f}_i^{ch}{\widehat{u}}^{\overline{h}}}, \mathrm{with} \overline{h}=n(c,h), f_i^{ch}={\varphi }_{,i}^{(c,h)}(x^c)$$

(42)

while from (41₁), we have

$$u_{,i}(x^c)\approx {\displaystyle \sum _{a=1}^{N^c}{\widehat{u}}_i^{\overline{a}}{\varphi }^{(c,a)}(x^c)}={\displaystyle \sum _{h=1}^{N^c}{e}^{ca}{\widehat{u}}_i^{\overline{a}}}, \mathrm{with} \overline{a}=n(c,a), e^{ca}={\varphi }^{(c,a)}(x^c)$$

(43)

Extending the definitions of matrices $f_i^{ch}$ and $e^{ca}$ to all nodes as

$$E^{cd}:=\left\{\begin{array}{c}{e}^{ca}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}d=\overline{a}\hfill \\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em} d\ne \overline{a}\hspace{0.17em}\hfill \end{array}\right., F_i^{cg}:=\left\{\begin{array}{c}{f}_i^{ca}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}g=\overline{a}\hfill \\ 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em} g\ne \overline{a}\hfill \end{array}\right., \mathrm{with} \overline{a}=n(c,a), a\in \{1,\hspace{0.17em}2,\hspace{0.17em}\dots \hspace{0.17em},\hspace{0.17em}{N}^c\}$$

(44)

The Equations (42) and (43) are rewritten as

$$u_{,i}(x^c)={\displaystyle \sum _{d=1}^N{E}^{cd}{\widehat{u}}_i^d}={\displaystyle \sum _{g=1}^N{F}_i^{cg}{\widehat{u}}^g}$$

hence,

$${\widehat{u}}_i^d={\displaystyle \sum _{g=1}^N{\displaystyle \sum _{c=1}^N{(E^{-1})}^{dc}{F}_i^{cg}{\widehat{u}}^g}={\displaystyle \sum _{g=1}^N{G}_i^{dg}{\widehat{u}}^g}}, G_i^{dg}:={\displaystyle \sum _{c=1}^N{(E^{-1})}^{dc}{F}_i^{cg}}$$

(45)

Differentiating (41₁), one can obtain

$$u_{,ij}(x)\approx {\displaystyle \sum _{a=1}^{N^q}{\widehat{u}}_i^{\overline{a}}{\varphi }_{,j}^{(q,a)}(x)}\Rightarrow u_{,ij}(x^c)\approx {\displaystyle \sum _{a=1}^{N^c}{\widehat{u}}_i^{\overline{a}}{\varphi }_{,j}^{(c,a)}(x^c)}={\displaystyle \sum _{g=1}^N{F}_j^{cg}{\widehat{u}}_i^g}$$

(46)

From (41₂), we have

$$u_{,ij}(x^c)\approx {\displaystyle \sum _{a=1}^{N^c}{\widehat{u}}_{ij}^{\overline{a}}{\varphi }^{(c,a)}(x^c)}={\displaystyle \sum _{d=1}^N{E}^{cd}{\widehat{u}}_{ij}^d}$$

(47)

and by comparison to (46) and (47), one can receive

$${\widehat{u}}_{ij}^d={\displaystyle \sum _{g=1}^N{\displaystyle \sum _{c=1}^N{(E^{-1})}^{dc}{F}_j^{cg}{\widehat{u}}_i^g}={\displaystyle \sum _{g=1}^N{G}_j^{dg}{\widehat{u}}_i^g}={\displaystyle \sum _{g,h=1}^N{G}_j^{dg}{G}_i^{gh}{\widehat{u}}^h}}$$

(48)

where we employed the expression for ${\widehat{u}}_i^g$ by Equation (45).

In a similar way,

$${\widehat{u}}_{ijk}^d={\displaystyle \sum _{g,b,h=1}^N{G}_k^{dg}{G}_j^{gb}{G}_i^{bh}{\widehat{u}}^h}$$

(49)

Substituting Equations (45), (48) and (49), into Equation (41), one can obtain the relationship for modified evaluation of derivatives of primary field variables. In this D1 approach, only the first-order derivatives of the MLS shape function ${\varphi }^{\overline{a}}(x)$ and nodal values ${\widehat{u}}^h$ are utilized for the evaluation of higher-order derivatives of the approximated function. The evaluation of $G_i^{dg}$, due to the necessity of inversion of $E^{cd}$ can lead to the lower computational efficiency of the modified differentiation approach as compared with the standard differentiation approach. This is the cost, what we have to pay for better accuracy given by the modified differentiation approach.

4. Numerical Examples

In this chapter, we shall present the results for FGM elastic and/or thermoelastic square plates with the power law gradation of material coefficients. The main geometrical parameter is the length of all sides of the plate $L=1$ m. The length to thickness ratio is $L/h_0=50$. The sides of the FG plates are considered to be clamped. The material coefficients are graded according to Equation (6) and/or Equation (25), while Poisson’s ratio is considered to be constant $\nu =0.3$.

The nodes utilized for the discretization of the analyzed domain are distributed uniformly. The distance between neighbour nodes is marked by $\delta$. The other parameters in the MLS-CAN approximation technique are the radius of the interpolation domain ${\rho }^a=3.001\hspace{0.17em}\delta$ and shape function parameter $c^a=\delta$. In all numerical experiments the cubic polynomial basis $m=10$ is used.

4.1. Verification of Presented Numerical Method

In the literature, we can easily find the analytical solution for clamped elastic square plates with constant bending stiffness [51]. For the study of the accuracy and convergence of the proposed computational method, the above-mentioned analytical solution will be utilized as benchmark solution. The error norm characterizing the accuracy of numerical solutions of BVP is given by

$$error=\frac{{\left({\displaystyle \sum _{a=1}^N{\left[w({\mathbf{x}}^a)-w^{ex}({\mathbf{x}}^a)\right]}^2}\right)}^{1/2}}{{\left({\displaystyle \sum _{a=1}^N{\left[w^{ex}({\mathbf{x}}^a)\right]}^2}\right)}^{1/2}}\hspace{0.33em}100(\%)$$

(50)

where N represents the number of all nodes within the discretized domain, $w({\mathbf{x}}^a)$ is the computed value of the field variable by the present method, while $w^{ex}({\mathbf{x}}^a)$ is the exact value of the field variable at the nodal point ${\mathbf{x}}^a$.

First of all, to verify the proposed computational method, we shall present the comparison of results by the proposed method and the analytical solution by [51] for elastic plate bending problems. The homogeneous, as well as FGM square plates, are considered. The plates in all computations are subjected to uniform, stationary transversal loading $q*=1$. The power law transversal gradation of elastic modulus follows Equation (6), while elastic modulus at the bottom surface of the plate is $E_{(0)}^{}=120\times {10}^9\mathrm{Pa}$, the level of power law gradation $\zeta =0$ for the plate with constant bending stiffness, while $\zeta =1$ and $\zeta =3$ for the FGM plate with linear $p=1\hspace{0.33em}$ and/or non-linear $p=2$ transversal gradation, respectively. In Figure 2. it is seen that the results for the deflections by the present computational method are in excellent agreement with the analytical results [51].

A study on the accuracy, convergence, and computational efficiency of the present method was carried out too. The results of such investigation are shown in Figure 3. Good convergence rates are achieved by both employed differentiation techniques with increasing the number of nodes, i.e., with decreasing parameter δ. In regards to the computational efficiency, from Figure 3b it is clearly seen that in the case of the modified differentiation (D1) the computational time is not increased significantly, as compared to standard differentiation (D0).

As regards the comparison of the computational efficiency of the strong vs. weak formulations using the same number of nodes, it depends also on the amount of discretized equations [46,52]. The total computational time is composed of the time needed for creation of the system matrix ($t_{sm}$) and the time needed for solution of the system of discretized equations ($t_{sol}$). Both the $t_{sm}$ and $t_{sol}$ are increasing with increasing the amount of nodes, and $t_{sm}$ is approximately 10 times smaller in strong formulation than in case of weak formulation, while $t_{sol}$ is independent on the kind of discretization. Thus, the computational efficiency of the strong formulation is ever better than that of the weak formulation, but the differences are diminishing with increasing the number of nodes when $t_{sol}$ is becoming dominant.

4.2. Study of $w-u_i$ Coupling in FGM Plates under Static Mechanical Loading

Now, we shall consider the FGM plate under static mechanical loading $q*=1$, $t=0$(the acceleration terms of Equation (20) are equal to zero).

First, we shall investigate the FGM plate with the transversal power law gradation of elastic modulus, focusing on the effect of $w-u_i$ coupling. According to [49], it is evident that the coefficients $C^{(uw)}$, $C^{(u\phi )}$, $C^{(wu)}$, $C^{(w\phi )}$, responsible for coupling between in-plane displacements $u_{\alpha }$ and bending measures (transversal deflections $w$ and rotations ${\phi }_{\alpha }$), are proportional to the parameter $\zeta$ and vanish in homogeneous plates, where $\zeta =0$. In order to assess the effect of material gradation on the $w-u_i$ coupling, we introduce the multiplication factor $c$ in the above coupling coefficients by replacing $C^{(··)}$ with $c{C}^{(··)}$. Then, the value $c=1$ corresponds to the real physical FGM plate, while $c=0$ corresponds to the artificial situation when the $w-u_i$ coupling in the FGM plate is switched off. Furthermore, the plane stress ($\chi =1$) and the 3D stress assumption are considered ($\chi =2$) due to the questionability of plane stress formulation in FG plates with the transversal gradation of elastic modulus. Figure 4 and Figure 5 present the three-dimensional map of variation of the deflection and in-plane deformation fields with in-plane coordinates within homogeneous as well as FGM plates. Figure 6a shows the effect of power law gradation of elastic modulus on the reduction of deflections compared to the deflections of homogeneous, as well as the influence of two factors $c$ and $\chi$ on deflections. Figure 6b illustrates that the in-plane displacements are induced only if the coupling factor is not switched off ($c=1$), and the effect of $\chi$-factor on in-plane displacements is shown too. The reduction of the maximum deflection is around 20% lower in the case of plane stress formulation ($\chi =1$) than in the case of 3D stress formulation ($\chi =2$).

In Figure 7a we can see that the reduction of the maximum value of deflection with respect to the $w_{ref}^{*}=w^{*}(\zeta =0,\hspace{0.17em}\chi =1)$ is increasing with increasing the level of gradation under keeping the volume contents of the constituents to be constant ($p=const$). In Figure 7b we can see that the increasing value of the exponent $p$ is leading to the decreasing reduction in the maximum deflections. This can be explained by the fact that the volume content of the constituent with higher elastic modulus is decreasing with increasing the exponent $p$.

4.3. Response of FGM Plates under Transient Mechanical Loading

In this chapter, we shall present numerical experiments for the considered FGM plates with clamped edges, and transversally graded elastic modulus and mass density, subjected to Heaviside impact load. The initial conditions of IBVP are vanishing. Various combinations of the power law gradation exponents (p, g) and levels of gradation ($\epsilon ,\zeta$) are considered. We consider also damping of the plate vibrations in view of the Rayleigh proportional damping model $\left[C\right]=\alpha \left[M\right]+\beta \left[K\right]$ [53], with choosing the parameters $\alpha =0.001$ and $\beta =0.0001$. To a numerical solution of the second-order ordinary differential equations (ODE) we apply the Newmark method [54] which is the implicit one-step time integration method.

Firstly, we shall present the time variation of central deflection and in-plane deformations for various combinations of gradations of the elastic modulus and mass density (Figure 8):

a.

Without gradation of material coefficients ($\epsilon =0,\zeta =0$).

b.

With transversal gradation of mass density only ($\epsilon =1,\zeta =0$).

c.

With transversal gradation of elastic modulus only ($\epsilon =0,\zeta =1$).

d.

With transversal gradation of elastic modulus and mass density ($\epsilon =1,\zeta =1$).

From Figure 8, one can see that both the amplitudes and frequencies of transversal deflection and the in-plane displacements are affected by the transversal gradation of the elastic modulus and/or mass density. The effect of the pure gradation of mass density on amplitudes is negligible as compared with the effect of the gradation of elastic modulus. However, the influence of each of the gradation coefficients $\left(\zeta ,\hspace{0.17em}\hspace{0.17em}\epsilon \right)$ on frequencies of vibrations is almost the same in magnitude but opposite in sign.

In what follows, we pay attention to the investigation of the role of higher values of the parameters $\zeta$ and $\epsilon$ and/or the gradation exponents $p$ and $g$ on the changes within the dynamic response of FGM plates with the transversal gradation of elastic modulus and mass density. From Figure 9 one can observe a significant influence of the level of gradation of each of the material coefficients on the frequencies of the vibrations of the central deflections and in-plane displacements. On the other hand, the amplitudes of vibrations are affected significantly only by the level of gradation of the elastic modulus. Since the frequency of oscillations is decreasing with the increasing value of ε and increases with the increasing value of ζ, one can tune the frequency rather sensitively by choosing the levels of gradations.

Figure 10 illustrates the study of the influence of the gradation exponents ($p,g$) on the dynamic response of the FGM plate. From Figure 10a we can see that the higher value of exponent p leads to a lower reduction in deflection, while the exponent g affected the frequency of vibrations remarkably, but the amplitudes are almost unaffected.

4.4. Response of FGM Plates under Transient Thermal Loading

The main focus of this chapter is on the investigation of the transient response of thermoelastic FGM plates subjected to thermal loadings defined as: ${\overline{\theta }}^{*}(0,x_2,x_3,t)=0$, ${\overline{\theta }}^{*}(1,x_2,x_3,t)=20H(t)$, ${\overline{q}}^{*}(x_1,0,x_3,t)=0$, ${\overline{q}}^{*}(x_1,1,x_3,t)=0$, ${\overline{q}}^{*}(\mathbf{x},\pm h/2,t)=0$, ${\overline{\theta }}^{*}(\mathbf{x},x_3,t=0)=0$, where $H(t)$ is standing for the Heaviside’s unit step function.

We shall present the results for temperature, in-plane deformation, and deflection fields for the case of pure transversal gradation of material coefficients. In Figure 11, it is illustrated that the transversal gradation of elastic modulus and/or thermal expansion coefficient has no influence on the time variation of temperature neither in the case of linear nor non-linear power law gradation. However, the power law gradation of material coefficients associated with heat conduction Equation (31), i.e., gradation parameters $\omega ,\hspace{0.17em}\chi ,\hspace{0.17em}\epsilon$(the gradation levels for heat conduction coefficients, specific heat capacity, and mass density) results in remarkable changes in the time evolution of temperature fields. It is seen that the linear and/or non-linear case of gradation of the heat conduction coefficient accelerates the time evolution of temperature fields, while the gradation of the mass density and/or specific heat capacity decelerates the time evolution of temperature fields.

The time variation of in-plane displacements is presented in Figure 12. In contrast to the temperature fields, the in-plane displacements are clearly affected by the power law gradation of the thermal expansion coefficient ($\alpha$); however, the impact of the gradation of elastic modulus ($E$) is negligible. The effect of the power law gradation of the heat conduction coefficients, specific heat capacity, and mass density on the time variation of in-plane displacements is similar to that on temperature fields.

From Figure 13. it is obvious that the power law gradation of elastic modulus and thermal expansion coefficient plays an important role in the deflection response of FGM thermoelastic plates subjected to transient thermal loadings. It can be seen that the transversal gradation of other material properties associated with the heat conduction equation has no influence on the deflection response of the FGM plate.

5. Conclusions

Three main conclusions can be drawn from the research presented in this paper:

(i)

The proposed unified formulation for plate bending problems allows:

• A unique treatment of plate bending problems with simple switching among three plate bending theories (classical Kirchhoff–Love theory, first-order shear deformation theory, third-order shear deformation theory), differing in various deformation assumptions;
• Physically correct derivation of governing equations and possible boundary conditions using variation principles;
• To study the response of linear elastic plates with functionally graded material properties to static as well as dynamic mechanical and thermal loadings;
• Comparison of results by three various theories using the same mathematical treatment.

(ii)

The developed advanced meshless method is characterized by:

• Efficient solution of systems of partial differential equations with variable coefficients (due to functional gradation of material coefficients) with the same demands as in the case of homogeneous media;
• Decomposition of the original fourth-order PDE to the system of the second-order PDE;
• Enhanced accuracy of approximation of higher-order derivatives;
• Improvement of computational efficiency by using the strong formulation;
• Overcoming shortcomings of the standard finite element method with preserving its universality.

(iii)

A study of coupling effects by numerical simulations revealed:

• The reliability (convergence and accuracy) and computational efficiency of developed numerical techniques;
• The functional dependence of material coefficients gives rise to coupling effects among the field variables including the interaction between the in-plane deformation and bending modes;
• To assess the influence of functional gradation parameters on the response of FGM plates to external loadings.

The proposed treatment of FGM plates could become a useful tool for designers because of the unified formulation including the assumptions of three various plate bending theories and owing to the wide, universal advanced numerical method for multi-physical plate bending problems.