A three-dimensional surface elastic model for vibration analysis of functionally graded arbitrary straight-sided quadrilateral nanoplates under thermal environment

In this paper, a three-dimensional (3D) size-dependent formulation is developed for the free vibrations of functionally graded quadrilateral nanoplatessubjectedtothermalenvironment.TheplatemodelisconstructedwithintheframeworksoftheGurtin–Murdochsurfaceandthe3Delasticitytheories.Inthisway,theeffectofsurfacefreeenergyandallthecomponentsofstressandstraintensorsareconsideredwithoutanyinitialassumptiononthemasthereisnoneedtoassumethevariationoftransversenormalstressinsidethebulkmaterialinadvance.Thevariationaldifferentialquadratureapproachandthemappingtechniqueareappliedtoderiveadiscretizedweakformofthegoverningequations.Thepresentsolutionmethodbypassesthetransformationanddiscretizationofthehigherorderderivativesappearingintheequationsofthestrongform.Theeffectsofsurfacestress,thermalenvironment,materialgradientindexandgeometricalpropertiesonthesize-dependentvibrationalbehaviorofquadrilateralnanoplatesareinvestigated.Itisobservedthatthethermalloadintensifiestheeffectofsurfacefreeenergyonthenaturalfrequencyofthenanoplates.Thepresentmodelisexactintheextentofthecontinuummodelsandcanbeemployedforstructureswithanythickness–spanratios.


INTRODUCTION
Mechanical behavior of nanoscale structural elements including nanobeams and nanoplates is the focus of attention in many research works due to their several applications in many devices, especially nanoelectromechanical systems (NEMS). Among the theoretical approaches used for the mechanical analysis of nanostructures, continuum models are widely used because their computational cost is less than that of atomistic approaches like molecular dynamics simulations.
A critical issue in the continuum modeling of nanostructures is related to the size-dependent behavior of materials at nanoscale. The size effects become significant when physical systems are scaled down to deep submicron dimensions. The size effects, due to longrange interatomic interactions, can be accounted for by the modified continuum theories such as the nonlocal elasticity, strain gradient and couple stress theories [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. In this respect, Malekzadeh et al. [1] studied the small-scale effect on the free vibration response of orthotropic arbitrary straight-sided quadrilateral nanoplates using the nonlocal elasticity and first-order shear deformation theories. Alibeygi Beni and Malekzadeh [2] analyzed the free vibration of non-prismatic skew nanoplates with different boundary conditions considering the small-scale effect by means of the same theories. Ansari et al. [3] incorporated the interatomic potential into a nonlocal plate theory in order to develop a Young's modulus-independent nonlocal model for investigating the vibrations of graphene sheets. Mohammadi et al. [4] employed the nonlocal elasticity theory to describe the size-dependent shear buckling behavior of an orthotropic rectangular nanoplate in thermal environment. Natural frequencies of the functionally graded (FG) nanoplates were investigated by Khorshidi et al. [5] based on the exponential shear deformation theory and nonlocal continuum mechanics. The couple stress theory was employed by Ansari et al. [6] to develop a size-dependent plate model for the forced vibrations of FG microplates under a harmonic excitation. Ghorbanpour-Arani et al. [7] adopted the sinusoidal shear deformation and Eringen's nonlocal theories to study dynamic buckling of FG nanoplates resting on an orthotropic elastic medium with inclusion of size effects. Farajpour et al. [8] considered the higher order deformations along with the higher and lower order nonlocal effects for the bucking of orthotropic nanoplates in thermal environment using the higher order nonlocal strain gradient theory. Kolahchi [9] carried out a comparative study on the static response, the natural frequencies and the buckling loads of sandwich nanoplates based on nonlocal refined zigzag, sinusoidal shear deformation, Mindlin's and classical plate theories using the differential cubature method. He found that the refined zigzag theory can accurately predict the static and vibration responses of sandwich nanoplates. Based on nonlocal plate theory, Jamalpoor et al. [10] determined closed-form solutions for natural frequencies and buckling loads of double-elastic nanoplate systems under initial external electric and magnetic potentials. Also, using an analytical Hamiltonian-based method, Fan et al. [11] derived exact solutions for size-dependent forced vibration of completely free orthotropic rectangular nanoplates based on Kirchhoff plate theory. Ghorbanpour-Arani et al. [12] adopted a nonlocal refined zigzag theory to investigate the buckling behavior of embedded nanoplates subject to external and magnetic potentials. Thermoelastic vibration of nanoplates was studied by Zenkour [13] using a mixed variational formulation based on nonlocal elasticity and the classical Mindlin's plate theories. Based on the nonlocal strain gradient theory and Navier method, Sharifi et al. [14] studied analytical vibrations of the FG nanoplates integrated by piezoelectric layers. Also, a Navier solution was developed by Zur et al. [15] for free vibration and buckling of sandwich FG material (FGM) nanoplates integrated with a piezomagnetic face sheet based upon the modified higher order sinusoidal shear deformation theory and nonlocal continuum mechanics. Three-dimensional (3D) continuum models were constructed for nonlocal free vibration behavior of a rectangular nanoplate [16] and nonlocal buckling of an embedded composite nanoplate made of quasicrystal coatings [17].
Among various size effects, the influence of surface free energy can be mentioned that dates back to the work of Gibbs [18] who presented the concept of surface stress in solids. This concept can be explained as follows. Since the energy of atoms at and near a surface of a solid body is usually different from that of atoms in the bulk phase, forming surface results in an additional free energy named as the surface free energy; the surface stress is defined based on the variation of the surface free energy with the surface strain [19]. After Gibbs, several researchers investigated the effects of surfaces on the behavior of solids (see, for example, [20][21][22]).
Surface stresses in a large-scale structure can be ignored because of their small values as compared to bulk stresses. However, in nanostructures that have high surface-to-volume ratios, the surface free energy and surface stress play important roles in the mechanical behavior and make the response of structure size dependent. Therefore, in the past two decades, many attempts have been made at capturing the effects of surfaces on the behavior of nanomaterials [23][24][25][26][27][28][29].
Among approaches that belong to the continuum-based group, the model developed by Gurtin and Murdoch [30,31] in 1970s has attracted a lot of attention from the researchers to analyze the mechanical characteristics of nanoscale structures with the consideration of surface effects. The applications of the Gurtin-Murdoch (GM) model to different problems of nanostructures such as vibration [32][33][34][35][36][37], bending [38,39], pull-in instability [40,41] and buckling [32,[42][43][44][45][46] in both linear and nonlinear regimes have been reported in the literature.
Also, some researchers considered the coupling effects of surface energies, nonlocality and thermal loading using Eringen's nonlocal and surface elasticity theories in predicting the buckling behavior [47], isogeometric free vibration response [48] and pull-in instability [49] of circular nanoplates, thermomechanical vibration of double-layer nanosheets [50], vibrations of a nanoplate sinking in an incompressible fluid [51], and statics and dynamics of annular piezoelectric nanoplates [52]. Using the Navier method, Pang and Zhang [53] presented closed-form solutions for natural frequencies of the viscoelastic rectangular nanoplate considering high-order surface stress and small-scale effects. Recently, applying surface stress relations, nonlocal elasticity and nonlinear von Karman theories, Ebrahimi and Hosseini [54] investigated nonlinear dynamic instability of nanoplates embedded on a visco-Pasternak foundation.
Arbitrary shaped nanoplates such as rhombic and trapezoidal ones are broadly employed in modern industries like biomedical devices, NEMS, actuators and sensors [55]. So, to properly design and manufacture nanodevices, a thorough knowledge of the mechanical behaviors of these nanostructured materials is required. Against the large number of research works into the behavior of rectangular nanoplates, there are few studies on the size-dependent modeling of oblique nanoplates [1,2,[55][56][57][58][59][60]. In the majority of available works, the effect of surface free energies has been ignored. Malekzadeh et al. [61] developed a size-dependent nonlinear plate model including both the small-scale and surface effects using the classical plate theory in conjunction with the nonlocal and surface elasticity theories for the free vibrations of skew nanoplates. Karimi and Shahidi conducted buckling [62] and free vibration [63] analyses of skew magnetoelectrothermoelastic nanoplates based on a refined plate model and surface and nonlocal elasticity theories. The governing equations were solved by the Galerkin approach and verified by the Navier solution.
To the authors' knowledge, no research has been carried out on the 3D modeling of oblique nanoplates accounting for the influences of surface stress and thermal loading. In the existing size-dependent continuum models based on the GM surface theory, the bulk substrate is modeled using two-dimensional (2D) theories. As well known, such theories ignore the transverse normal stress in the thickness direction. This shortcoming causes some of the interface boundary conditions between the bulk and the surface layers, and consequently equilibrium equations on the surface not to be satisfied. This issue necessitates that the variation of normal stress inside the bulk material be assumed in advance [32]. Prompted by this, in this work, a variational formulation is proposed in order to study the 3D free vibration behavior of arbitrary straight-sided quadrilateral nanoscale plates made of FGMs and subjected to thermal environment in the presence of the surface effects. The 3D theory of elasticity along with the GM surface theory is employed to establish an exact size-dependent model without any initial assumption. Using this model, in contrast with the size-dependent 2D models, the transverse normal stress is directly determined from the solution of governing equations and does not need to be taken in advance.
Moreover, in the simulations of the non-rectangular structures available in the literature, a mapping method is used to transform the irregular physical domain into a regular computational one. The first-and higher order derivatives of displacement field in the physical domain are calculated in terms of the derivatives in the computational domain. This task needs the calculation of inverse transformation matrices that may be computationally expensive in some cases. In contrast, in the present solution approach called variational differential quadrature (VDQ) [59,[64][65][66], a variational framework is used to carry out the discretization through which the mapping process is conducted by calculating just the first-order derivative of the displacement in the physical domain. Also, by constructing integral operators whose role is the integration of energy functional in a computationally efficient manner, the set of discretized governing equations of weak form can be readily identified and solved. By this, the need for derivation of the strong statement of governing equations through minimizing the energy functional and then the transformation and discretization of the higher order derivatives is removed (against the conventional differential quadrature formulation in which the strong form of equations is discretized). This feature can be of great significance, as the derivation of differential equations through minimizing the energy functional may be a complicated task for the size-dependent continuum models. Thus, the VDQ method is computationally effective, while being of high accuracy and fast rate of convergence. Finally, selected numerical results are given to study the effects of surface stress, thermal environment, material gradient index and geometrical properties on the size-dependent vibrations of quadrilateral nanoplates under different edge supports.

CONSTITU TIVE REL ATIONS
Consider an FG arbitrary straight-sided quadrilateral nanoplate with thickness h as shown in Fig. 1. A Cartesian coordinate system (x, y, z) is located on the bottom surface of the nanoplate to label the displacement components of an arbitrary material point in the x, y and z directions as u, v and w, respectively. It is assumed that the nanoplate is made from a mixture of metal and ceramic. The bottom (s−) and top (s+) surfaces of the nanoplate are taken to be ceramic-rich and metal-rich, respectively. E s , ρ s and τ s show surface Young's modulus, mass density and residual stress, respectively. It is also considered that the material properties of the bulk (elastic material coefficients C ij and mass density ρ) vary in the thickness direction according to the power-law function as follows [64,67]: where the superscripts c and m stand for the ceramic and metallic phases, respectively, and k denotes the material gradient index or power-law exponent. For an isotropic FGM, C ij are related to elastic material properties, Young's modulus (E) and Poisson's ratio (ν) of the nanoplate, or equivalently Lamé constants λ and μ, by the following formulas [64,67]: in which the Poisson's ratio is taken to be constant through the thickness. Further, as the displacement field over the thickness of the nanoplate is considered to be continuous, the Poisson's ratio of the material surface and bulk remains the same [33]. The classical continuum mechanics is not capable of capturing the atomic features of the nanostructures. In this work, the modified continuum elasticity based on the GM surface theory is introduced to incorporate size effects into the classical continuum approach. Following the GM surface elasticity theory [30], it is considered that the nanoplate consists of two phases: bulk and surface. According to [30], the physical properties of the nanoplate in the neighborhood of its surface are taken to be different from those of its interior. Also, the surface phase is modelled as a material surface. Another important assumption is that the surface phase adheres to the bulk phase without slipping.
Herein, the 3D elasticity theory is used for writing the constitutive equations of bulk phase. Hence, no assumption is made on the components of the stress tensor and all of them are directly determined from a numerical solution. Based on the 3D linear small-strain elasticity theory, the components of the strain tensor corresponding to the bulk phase are related to the displacement components as follows: The components of the stress tensor are related to the components of the strain tensor via the material elastic coefficients matrix C. So, for an elastic isotropic material, the thermoelastic constitutive relations are given by where σ M = C Cε, T = T − T 0 is the temperature rise from the reference temperature T 0 at which the plate is stress-free and α denotes the coefficient of thermal expansion. Also, the matrices C and α are of the following form: In the pre-bucking stage corresponding to the thermal loading, the in-plane and out-of-plane behavior of the nanoplate is like that of a rigidly fixed plate in extension and a simply supported one in bending [68]. So, the in-plane and out-of-plane components of normal stress are obtained by Eq. (4) as follows: The buckling phenomenon is a significant mode of failure for thin to moderately thick FG plates [67]. Accordingly, the transverse stress components at the equilibrium state can be ignored [67,69]. Thus, setting σ Th zz = 0 in Eq. (8), one has zz = Substituting Eq. (9) back into Eqs (6) and (7), the absolute values of the pre-buckling stresses for the thermal loading are obtained as follows: σ Th The pre-buckling stresses in Eqs (10) and (11) are equivalent to those of [69] presented in terms of Lamé constants. Due to the interaction between the elastic surface and bulk material, the nanoplate is subjected to in-plane loads in various directions leading to surface stresses. These surface stresses can be evaluated by the constitutive relations of the surface phase as given by Gurtin and Murdoch [30,31] as follows: in which λ s = E s ν (1+ν)(1−2ν) and μ s = E s 2(1+ν) are surface Lamé constants. Moreover, τ s indicates the surface residual stress. Using Eqs (12) and (13), the components of surface stress at the top surface and the bottom surface of the nanoplate are derived as

GEOMETRICAL MAPPING
In order to employ the VDQ method, the 3D computational domain should be of the form of a cube. Thus, the non-rectangular physical domain needs to be mapped into a rectangular one in the computational space. To this end, a four-node nanoplate element is used in the natural ξ -η plane (see Fig. 2). Based on the standard procedure utilized in the finite element approach, the transformation formulations are given by [65] x = 1 4 in which x i and y i (i = 1,…, 4) denote the coordinates of node i in the physical domain. Using Eqs (19) and (20), the derivatives of any function defined in one domain can be obtained in terms of those of the other domain. By means of the chain rule, the first-order derivative of any function in the physical domain can be calculated in terms of the one in the computational domain as follows: Using the transformation matrix [T 11 ], the variational formulation is developed in the computational domain.

VDQ APPROACH 4.1 Discretization of the strain and stress components
Based upon the strain-displacement relations given in Eq. (3) and the transformation rule considered in Eq. (21), the vector of strain components for the bulk phase can be decomposed to where T ξη i j as the element located on the ith row and jth column of For computing the surface and bulk elastic energies, the strains related to the bulk and material surfaces are treated separate from each other. According to Eqs (14)- (18), the vector of the surface strain components corresponding to the surface stress components can be defined as whereĒ in the physical domain is expressed asĒ and in the computational domain using Eq. (21) one has Note that for any arbitrary vector V, the notation (V)| s ± indicates the components of V calculated at the surfaces. From Eqs (14) whereẼ, in the physical domain, and τ s are of the following form: and in the computational domain, regarding Eq. (21), one can arrive at Now, for discretizing the stress and strain components, the generalized differential quadrature (GDQ) technique is applied to construct the differential matrix operators. The displacement functions u(ξ , η, z, t), v(ξ , η, z, t) and w(ξ , η, z, t) are defined in the computational domain with discrete points distribution as [ξ 1 , …, ξ N ], [η 1 , …, η M ] and [z 1 , …, z P ]. The following column vectors of the size NM × 1 are then introduced: in which l = 1, 2, …, P. In fact, u l , v l and w l separately contain the functional values of u(ξ , η, z, t), v(ξ , η, z, t) and w(ξ , η, z, t) at the grid points (ξ i , η j , z l , t) with i = 1, 2, …, N and j = 1, 2, …, M. By means of Eqs (30)- (32), the nodal values of the displacement field over the whole domain can be calculated as where Note that u, v and w are column vectors of the size NMP × 1. By using the Kronecker tensor product ⊗ defined in Appendix and Eq. (34), the approximation of first-order derivatives of the displacement functions with respect to the coordinate directions, e.g. u(ξ , η, z, t), at the grid points over the whole domain can be taken in the following vector form: where the coefficients of u are those differential matrix operators in which D (1) z , D (1) η and D (1) ξ are the weighting coefficient matrices corresponding to the first-order derivative with respect to z, η and ξ , respectively, in the GDQ approach (see Appendix). Also, I z , I η and I ξ stand for identity tensors with P 2 , M 2 and N 2 components, respectively. According to Eqs (22) and (33)-(37), the vector of the bulk strain components is expressed in the discretized form as follows: where with T 11 11 = I z ⊗ diag t 11 11 , T 11 21 = I z ⊗ diag t 11 21 , andt 11 11 = vec t 11 11 ,t 11 12 = vec t 11 12 ,t 11 21 = vec t 11 21 ,t 11 22 = vec t 11 22 , T 11 22 (ξ N , η 1 ) T 11 22 (ξ N , η 2 ) · · · T 11 22 Definition of the operators diag( ) and vec( ) used in Eqs (42)-(48) is presented in Appendix. Likewise, the discretized vector of the surface strain components s ± given by Eqs (24) and (26) is written as in whichĒ denotes the discretized form of the matrixĒ in the computational space given bȳ Moreover, the vector of the surface stress components σ s ± , Eqs (27) and (29), is similarly transformed into the discrete form as in whichẼ s± is the discretized form of the matrixẼ at the top and bottom surfaces of the nanoplate.

Quadratic representation of energy functionals
The elastic strain energy of the bulk phase U b can be considered as the sum of mechanical and thermal works separately denoted by U Mb and U Tb , i.e. U b = U Mb + U Tb , where [65,67,69] or in a matrix notation where ∀ or∀ represents the volume of the elastic body so that ∫∀( )d∀ = ∫ ξ N ξ 1 ∫ η M η 1 ∫ z P z 1 ( )dξ dη dz and In the computational domain, with the use of Eq. (21), relationship (56) becomes with The kinetic energy of the bulk phase is also given by Analogous to the mechanical strain energy of the bulk, the surface strain energy can be written in a quadratic form. So, based on Eqs (24) and (27), the elastic strain energy due to the surface stresses is computed as where the capital letters S + and S − show the area of the top and bottom surfaces, respectively, and ∫ S ( )dS = ∫ ξ N ξ 1 ∫ η M η 1 ( )dξ dη . Also, the kinetic energy due to the surface effects can be computed like the kinetic energy of the bulk as Since U b = U Mb + U Tb , the total strain and kinetic energies of the nanoplate are then calculated as Work done by the surface residual stress is also formulated as Using Eq. (17), the preceding relationship can be rewritten in a matrix format in the computational space as where Now, the discretized forms of the stress and strain tensors as well as the displacement vector U are applied so as to transform the energy functionals into a discrete form. To this end, the Hadamard product "•" is also used that, herein, performs the element by element product of two vectors (see Appendix). By Eqs (38), (54) and (58), the discretized forms of the strain energy functionals of the bulk phase are written as In addition,σ Th xx ,σ Th yy and C signify square matrices of the size 3NMP × 3NMP and 6NMP × 6NMP, respectively, whose expressions and those of G ξ and G η are given bȳ with From Eq. (60), the kinetic energy functional of the bulk phase is also discretized as Based on Eqs (49)-(53), (61) and (62), for the discrete form of the surface free energies functionals, one can write Finally, using Eq. (66), for the discretized form of the functional corresponding to the work done by the surface residual stress, one has where q s+ = t s+ ⊗ I NMP , q s− = t s− ⊗ I NMP , G = 0 0 T ξη

Discretized governing equations and solution
The elastic energies are numerically calculated from the discretized form of energy functionals. This is accomplished by 2D integral operators for the surface free energies, and 3D integral operators are used to evaluate the energies of the bulk material. Consider the mode shape function u(ξ , η, z, t) defined in a cubic domain with the same mesh point distribution as before. Similar to the method used in the approximation of the partial derivatives of u(ξ , η, z, t), the triple integral ∫ ξ N ξ 1 ∫ η M η 1 ∫ z P z 1 u(ξ, η, z)dξ dη dz can be calculated via the Kronecker tensor product as in whichS z ,S η andS ξ are one-dimensional (1D) integral operators in the z, η and ξ directions constructed from the GDQ differential operators and Taylor series expansion (see Appendix). Besides, the double integral ∫ η M η 1 ∫ ξ N ξ 1 u(ξ, η, z, t )dξ dη can be approximated as where the coefficient of u is a 2D operator that integrates a function defined in a 3D computational domain.
To calculate the energies of the bulk phase, the integral of the functional is triple type, and the 3D operators similar to that of Eq. (82) are required, which are introduced as MultiplyingS ub by Eq. (68) andŜ tb by Eqs (69) and (76) results in the strain and kinetic energies of the bulk phase as For surface energies, 2D operators should be used. Since the energies of the top and bottom surfaces are calculated, the operators are defined in such a way that the integration is made only for these surfaces. In other words, the nodal values of the functional at the grid points distributed on the layers z = z 1 and z = z P , are only considered in the integration process. Thus, by Eq. (83), the following augmented operators are introduced: whereS us− andS us+ are separately utilized for the computation of the strain energy as a result of the bottom and top surface stresses. Also,S ts− andS ts+ are used in the evaluation of the resulting kinetic energy. Based on Eqs (78)-(80) and (87)-(89), the surface strain and kinetic energies as well as the work done by the surface residual stress are obtained as So, combining Eqs (85), (86), (90) and (91) based on Eqs (63) and (64), the total strain and kinetic energies are derived as The following operators are also introduced: which lead to Let a be a row vector, and b and c be column vectors, then a (b • c) = b T diag(a)c. By means of this equality and Eqs (96)-(101), one can rewrite Eqs (93) and (94) in the following form: Similarly, using Eqs (81), (92), (102), (103), and the introduced equality, the work done by the surface residual stress becomes Let the configuration of the nanoplate be prescribed at t = t 1 and t = t 2 , applying Hamilton's principle to the energy functionals given by Eqs (104)-(106) leads to Integrating by parts and then vanishing the coefficient of δU T gives This equation is the discretized form of the governing equations including the surface stress effects that can be stated in terms of domain matrices as follows: where To impose the suitable mechanical boundary conditions to the nanoplate taking the surface effects into account, the generalized resultant stresses are defined as [26] Let n x and n y represent the components of the unit normal vector to the boundary of the nanoplate in the x and y directions, respectively. The boundary conditions along the skew edge are written as either u n = n x u + n y v = 0 or either u s = −n y u + n x v = 0 or The boundary conditions on the surfaces z = 0, h are also written as By the substitution of Eqs (3), (4) and (14)- (18) into relationships (113)-(117), one can express the boundary conditions in terms of the displacement variables. For example, the expression σ nn in Eq. (114) becomes σ nn = C 11 n 2 x + C 12 n 2 y T 1 ξη + 2C 66 n x n y T 2 ξη u + C 12 n 2 x + C 22 n 2 y T 2 ξη + 2C 66 n x n y T 1 ξη v Similar process can be implemented for other components of the stress. For the simply supported (S) and clamped (C) boundary conditions, the mechanical boundary conditions at the edge are stated as follows: Clamped: Simply supported: To apply the edge conditions, their discretized form is directly inserted into the domain matrices given by Eqs (110)-(112) like what is generally done in the GDQ approach. In order to treat the free vibration of the nanoplate, any type of bending force should be dropped. So, if ω represents the natural frequency of the nanoplate, by setting F = 0 and employing the mode shape functions of the following general form: Eq. (109) can be recast to an eigenvalue problem on domain from which the natural frequencies and the associated mode shapes are evaluated.

RESULTS AND DISCUSSION
In this section, selected numerical results are presented on the vibrational characteristics of quadrilateral FG nanoplates under thermal environment with various geometries considering surface stress effects. Schematics of the rhombic (skew) and quadrilateral nanoplates with the geometrical parameters considered in the numerical calculations are shown in Figs 3 and 4, respectively.  It is assumed that the temperature of the whole material varies from a reference value (T 0 ) to a final value. The reference temperature is taken as T 0 = 300 K. It is also considered that the nanoplates are subject to different boundary conditions, including SSSS, CCCC, CCSS and CCSC. Also, it is assumed that the FGM consists of Al as metal phase and Si as ceramic phase with the following bulk and surface properties: Moreover, the Poisson's ratio of the FG nanoplate is considered to be ν = 0.27, otherwise specified.
The convergence of presented numerical solution approach is studied in Table 1. In this table, the convergence of dimensionless fundamental frequencyω i = ω i h 2(1 + ν )ρ c /E c of a simply supported skew nanoplate is given for two values of a/h and three values of skew angle θ . One can find that the results tend to converge by increasing the number of discrete points. Table 1 shows that converged results are obtained for N = M ≥ 13 and P ≥ 9 demonstrating a fast convergence speed. Therefore, N = M = 13 and P = 9 are selected for generating the following numerical results.
Also, to validate the present approach, a comparison is made between the present results and those reported by Assadi [33] in Table 2. This table indicates the second natural frequency to first natural frequency ratio of a simply supported homogeneous rectangular nanoplate for different aspect ratios and thicknesses based on the surface elasticity theory. The material properties of the nanoplate are considered as = 68.5 × 10 9 N/m 2 , ρ = 2700 kg/m 3 , ν = 0.35, τ s = 0.910 N/m, E s = 6.090 N/m and ρ s = 5.46 × 10 −7 kg/m 2 . It is observed that the present results are in close agreement with the ones reported in [33]. Table 1 Convergence of dimensionless fundamental frequency of a simply supported FG skew nanoplate considering surface effects when k = 1, a/b = 1, a = 10 nm and T = 200 K. In Figs 5 and 6, the first frequency parameter of FG rhombic nanoplates is plotted versus temperature change. The frequency parameter is defined asω i = ω i a 2 √ ρ c h/D in which D = E c h 3 /(12(1 − ν 2 )). In Fig. 5, the effect of skew angle on the curves can be seen, whereas Fig. 6 shows the effect of material gradient index. Both figures illustrate that the frequency quickly decreases, as the temperature difference goes up. If the temperature difference further increases, the frequency becomes smaller and smaller until it approaches zero at buckling load. The numerical results corresponding to different boundary conditions show that for an imposed temperature change, the natural frequencies of the SSSS and CCCC nanoplate are the smallest and the largest, respectively. Also, the natural frequency of the former approaches zero at a lower temperature difference as compared to the latter. So, it is concluded that the stiffer the edge support, the larger the natural frequency and then the larger the thermal buckling load. Moreover, from Fig. 5 one can find that increasing the skew angle has a decreasing influence on the first frequency parameter of FG rhombic nanoplates. In other words, it is discerned that the skew nanoplates have larger vibration natural frequency as well as buckling load than the rectangular ones.
Furthermore, from Fig. 6 it is seen that an increase in the martial gradient index leads to the increase of the frequency since it causes the volume fraction of the ceramic phase and therefore the stiffness of the nanoplate to rise.
To examine the effect of thermal loading on higher modes of vibrations, the frequency parameters of the FG skew nanoplate corresponding to the lowest three modes are plotted versus the temperature change in Fig. 7. It is noted that with increase of the temperature, the resonant frequency of higher modes reaches zero at larger values of temperature difference as compared to that of fundamental mode. It is expectable because the thermal buckling load related to modes 2, 3 and higher is larger than the fundamental buckling load. Figures 8 and 9 show the variations of the first natural frequency ratio of FG rhombic nanoplates with thickness for various skew angles and selected materials, respectively. The natural frequency ratio is defined as the ratio of calculated frequency considering surface effects to the calculated frequency without considering surface effects. In fact, Figs 8 and 9 illustrate the effect of surface stress on the vibrational behavior of nanoplates. As observed, the natural frequency ratio of all the nanoplates is obtained to be larger than unity. It means that for the selected material properties of the bulk and surface, the surface effect leads to the increase of the natural frequency. It is observed that such effect is very pronounced at small values of plate thickness. However, the surface effects tend to diminish as the thickness increases, and these effects can be neglected for sufficiently thick nanoplates. It can be explained by the fact that when the thickness is small, the surface free energy is considerable as compared to the energy associated with the bulk phase. However, by increasing the thickness, the surface free energy can be ignored in comparison with the energy of bulk part.  In addition, Fig. 8 reveals that the frequency ratio becomes larger when the skew angle increases. In other words, the effects of surfaces on the natural frequency of the rectangular nanoplates are more prominent as compared to the skew ones. Furthermore, according to Fig. 9, at a given value of thickness, the frequency ratios of nanoplates made of Si are smaller than those of nanoplates made of FGM, and the latter are also smaller than those of nanoplates made of Al. As a finding, the effect of surface energies on the natural frequency of a fully ceramic nanoplate is less pronounced than the effect on the frequency of an FGM or a fully metal nanoplate.
Depicted in Fig. 10 is the variation of natural frequency ratio of the FG rhombic nanoplate with thickness corresponding to the first three modes of vibration. It is seen from the figure that, as the mode number increases, the frequency ratio decreases so that the fundamental mode is of the largest frequency ratio. In other words, the effect of surface stress on the frequency of the higher modes of vibrations is less significant than the fundamental mode. This result is confirmed by the observation made on higher modes of buckling behavior of Mindlin rectangular nanoplates in [38].
In Fig. 11, the first natural frequency ratio is plotted against temperature change corresponding to three values of the aspect ratio b/a. The results indicate that the difference between the natural frequencies obtained by taking surface effects into account and those computed by neglecting such effects intensifies as the temperature change goes up. It means that when a nanoplate is subjected to a thermal environment, the influence of surface stress on its vibrational behavior becomes more prominent. It is also observed that the increase in the aspect ratio has an increasing effect on the frequency ratio or on the surface energy effects, especially for large temperature differences. It is due to the fact that for a fixed value of a, an increase of the aspect ratio leads to the increase of the surface area, and therefore, it causes the surface effects to rise.
A plot of natural frequency ratio of the FG rhombic nanoplate versus temperature change associated with the lowest three modes of vibration is shown in Fig. 12. Again, the first mode is observed to have the largest natural frequency ratio. Moreover, numerical results reveal that as the temperature difference goes up to 550 • K, the frequency ratio of modes 1, 2 and 3 increases by ∼71%, ∼25% and ∼14%, respectively. As a finding, the intensification of surface stress effect due to thermal loading on the first vibrational mode is much greater than other vibrational modes.
Presented in Fig. 13 are the variations of first frequency parameter of the FG quadrilateral nanoplate shown in Fig. 4 with asymmetry angle (β) for various cord ratios (b/c). Note that β = 0 gives a symmetric trapezoidal nanoplate. The figure shows that the frequency parameter increases by increasing the asymmetry angle. As a finding, the frequency of an asymmetric quadrilateral nanoplate is higher than that of a symmetric one. Another observation is that at a given asymmetry angle, the nanoplates vibrate at larger frequencies when a larger value of the cord ratio is selected.
Finally, the influence of thermal environment on the frequency ratio-cord ratio curve of a fully simply supported FG asymmetric quadrilateral nanoplate is highlighted in Fig. 14. It is noted that as the cord ratio goes up, the frequency ratio diminishes. Because with the increase of the cord ratio (by keeping the parameters a and b unchanged), the surface area of the nanoplate is reduced, the surface  effects get nullified. Again, it is also found from this figure that an increase in temperature difference makes the surface stress effects more pronounced, especially at small core ratios.
It is worth mentioning that as the present size-dependent model is based on the 3D theory of elasticity and bypasses any assumption on the components of stress and strain tensors, it is more precise than its 2D counterparts. The model provides direct calculation of transverse shear and normal stresses satisfying the boundary conditions on the top and bottom surfaces from field equations. This can be worthy in the case of dynamic or bending and thermal buckling simulations of nanostructures in which the through-thickness stress distribution may be important. In dynamic problems where the time response is of interest, the model can provide the examination     of the effect of surface on the stress distribution at any time in addition to that on the resonant frequency, and in the static problems, the effect of surface on the stress distribution can be identified at critical loads. Solving these problems may be viewed as future perspectives of the current model.

CONCLUSIONS
In the context of the GM model accounting for effects of the surface stresses, the free vibrations of arbitrary straight-sided quadrilateral nanoplates were analyzed in this work based on an efficient numerical variational approach. It was considered that the nanoplates are made from FGMs and are subjected to thermal environment. For the constitutive relations of bulk phase of nanoplate, the 3D elasticity theory was utilized by which no assumption was made on the components of the stress tensor, and all of them were directly determined from a numerical solution. This feature makes the proposed model promising to be employed for structures with any thickness-span ratios. The results of the model would be used as a benchmark to assess the accuracy of the 2D models including the surface effects. A DQ-based variational method named as VDQ was adopted to formulate a weak form solution of the considered problem. In this approach, differential operators using the DQ rule were first introduced so as to discretize the bulk and surface strains.
In the next step, integral operators were constructed to calculate the elastic energies of the nanoplate. Then, the set of discretized governing equations in terms of domain matrices was derived based on Hamilton's principle. The convergence and validity of the calculations were first shown, and then some numerical examples were given to investigate the vibrational responses of FG rhombic and asymmetric quadrilateral nanoplates with different boundary conditions. In these examples, the effects of surface stress, material gradient index, skew angle, asymmetry angle, aspect ratio, cord ratio and temperature change were studied. It was shown that the surface stress considerably affects the frequencies of very thin nanoplates due to the important role of surface free energy. It was also revealed that such effect is dependent on the skew angle and the selected material for the nanoplate so that the fully ceramic skew nanoplates are less affected by the surface energies as compared to the rectangular FGM ones. Another conclusion was that the thermal environment can decrease the natural frequency and make the surface effects more significant, especially when the nanoplate vibrates at the fundamental frequency. Also, investigation of the different modes of vibration showed that the resonant frequency of fundamental mode is the most sensitive to the surface stress effect as compared to that of higher modes. The frequency of the nanoplate with stiffer edge constraints and oblique edge was observed to approach zero at a higher temperature difference. Moreover, the results indicated that increasing the asymmetry angle leads to the increase of the first frequency. In future works, we aim to use the same formulation for addressing other problems including bending and buckling. We also aim to extend the present formulation from linear to nonlinear. The capability of the present 3D model in reflecting the surface effect and all the stress components may lay a foundation for future research on nonlinear dynamics of nanostructures. (A2)

A.2 1D integral operator
Let the one-variable function h(k) be defined in [k 1 , k Q ], (k = ξ , η, z; Q = N, M, P). The integral of g(k) over this interval can be approximated by where h includes the functional values of h(k) andS k is read as the integral operator. Taylor series expansion of the function h(k) about the point k i leads to Combining two preceding equations and using the GDQ approach for the approximation of the partial derivative of a function, one can get