Buckling analysis of functionally graded microplates incorporating nonlocal strain gradient theory and surface energy effects

Abstract

The purpose of this work is to present a buckling analysis of functionally graded (FG) microplates by combining nonlocal strain gradient theory (NSGT) with the higher-order shear deformation plate theories (HSDPTs). The microplate is assumed to be composed of a combination of ceramic and metal materials, and the material properties are assumed to vary continuously in the thickness direction based on a simple power law. The equilibrium equations and the boundary conditions are derived using the principle of minimum potential energy. Analytical solutions are determined for the critical buckling loads of the rectangular FG microplates with different boundary conditions. The results obtained have been verified by comparison with existing findings in the literature. Furthermore, some numerical illustrations are provided to investigate the effects of nonlocal parameters, the material length scale parameter, shear deformation, aspect ratios, the power-law index, and the surface energy on the buckling response of the rectangular FG microplates.

Introduction

Currently, scientists and engineers acknowledge the significance of using innovative materials to meet both economic and environmental goals. Functionally graded materials (FGMs) are advanced engineering materials that are specifically designed to exhibit certain performance characteristics or functionalities. These materials have a spatial gradation in their structure and/or composition that allows their properties to be customized according to specific requirements. This is achieved through the provision of comprehensively graded compositions, microstructures, and material properties. FGMs were initially introduced in Japan in 1984 to be developed and applied as thermal barrier materials. ¹ With their subsequent widespread application in many different fields, FGMs have attracted the attention of many scientists, particularly materials and structural scientists. Numerous works have been published on the mechanical behavior of FG structures over the last 30 years.^2–4

With the rapid ongoing development of science and technology, micro/nanostructures, which include micro/nanobeams, micro/nanorods, micro/nanowires, micro/nanoplates, and micro/nanoshells, are being applied increasingly widely in various microelectronic devices, including micro/nano-electromechanical systems (MEMS/NEMS). Accurate comprehension of the mechanical performances of these small-sized structures is highly important because they are the fundamental elements in MEMS and NEMS devices. In contrast to macro-sized structures, several experimental studies have concluded that the mechanical behavior of these small-sized structures is strongly size-dependent.^5–7 Classical elasticity theory, because of its lack of length scale parameters, can no longer be used to describe the mechanical behavior of these small-sized structures. To overcome this issue, several nonclassical elasticity theories that integrate the required length scale parameters have been proposed to describe the mechanical responses of these small structures more fully, including nonlocal elasticity theory (NET),^8,9 strain gradient theory (SGT),^10,11 and modified couple stress theory (MCST). ¹² With the integration of additional nonlocal parameters, NET can capture the softening-stiffness effect in nanostructures, while the hardening-stiffness effect can be observed through SGT and MCST because of the introduction of the material length scale parameters. These nonclassical elasticity theories have been used very effectively to describe the mechanical behavior of micro/nanostructures.^13–17

Recently, by integrating NET^8,9 and SGT^10,11 into a generalized nonclassical elasticity theory, Lim et al. ¹⁸ proposed the nonlocal strain gradient theory (NSGT). By integrating the nonlocal parameters and the material length scale parameters into a single elasticity theory, both the softening-stiffness effect and the hardening-stiffness effect can be captured through NSGT. ¹⁸ Many works related to analysis of the mechanical behavior of micro/nanostructures using this theory have been published in the literature. The size-dependent nonlinear bending and free vibration behavior of Euler-Bernoulli and Timoshenko FG beams were examined by Li and Hu. ¹⁹ Using Euler-Bernoulli beam theory and NSGT, Şimşek ²⁰ investigated the size-dependent nonlinear vibration characteristics of an FG nanobeam by considering the position of the neutral surface. Dang et al. ²¹ also used NSGT and Euler-Bernoulli beam theory to study the size-dependent nonlinear vibration behavior of a nanobeam when subjected to electrostatic loading. The size-dependent nonlinear vibration and stability behaviors of an FG porous electrostatically actuated microbeam were reported by Dang and Do ²² , who also combined Euler-Bernoulli beam theory with NSGT. Anh and Hieu ²³ examined the size-dependent nonlinear vibration characteristics of an FG nanobeam that had been subjected to random excitation.

Because of the high surface area-to-bulk volume ratio in micro/nanostructures, the energy associated with the surface atoms becomes very important when compared with that of the bulk part of these small-sized structures. The surface energy effects can be observed based on the surface elasticity theory (SET) proposed by Gurtin and Murdoch,^24,25 in which a solid is assumed to be composed of a bulk component and surface layers of negligible thickness that are perfectly bonded to the bulk part. Accordingly, constitutive equations for the bulk part and the surface layers have been established separately. Hosseini et al. ²⁶ studied the surface energy effect on the biaxial buckling and free vibration behavior of FG nonlocal nanoplates when resting on an elastic foundation based on Kirchhoff plate theory. They concluded that the damped natural frequency of the nanoplate is increased significantly by increasing the number of higher modes, whereas the buckling force ratio is reduced significantly by increasing the number of higher modes when the surface energy effects are considered. ²⁶ The free linear vibration behavior of Kirchhoff and Mindlin plates was investigated by incorporating the surface energy effect and NSGT by Lu et al. ²⁷ . They found that the surface energy effects have an important influence on the vibration behavior of nanoplates with higher length-to-thickness ratios. ²⁷ The isogeometric analysis (IGA) technique was used by Yin et al. ²⁸ to examine the size-dependent static bending and free vibration problems of nonclassical Bernoulli-Euler beams by considering the surface energy effects. The authors found that the influence of the microstructure on the structure’s natural frequency is more significant than the surface energy effect on the nanoscale. ²⁸ The surface energy effect on the thermoelastic vibration behavior of NEMS was investigated by Kumar and Mukhopadhyay, ²⁹ who combined NET and Moore–Gibson–Thompson thermoelasticity theory. The authors concluded that the effects of the surface residual tension and the surface elastic modulus on the deflection and temperature responses of NEMS are insignificant for short time intervals during vibration, but they are considerable for longer time intervals. ²⁹ Based on MCST and Kirchhoff plate theory, the dynamic instability responses of thin nanoplates with consideration of the surface energy effects were investigated by Nguyen and Phan ³⁰ using the IGA method. The results obtained showed that the influence of the surface energy becomes dominant in the case of ultra-thin nanoplates. ³⁰

It should be noted that because it ignores the effects of shear deformation, Kirchhoff plate theory is only suitable for modeling of thin plates. To overcome this limitation of Kirchhoff plate theory, many alternative shear deformation plate theories that incorporate the influence of the transverse shear deformation effect have been developed, including the Mindlin plate theory, ³¹ in which a shear correction factor is required to correct unrealistic variations in the transverse shear stresses and shear strains through the thickness of the plate. To avoid the requirement for the shear correction factor, several higher-order shear deformation plate theories (HSDPTs) have been proposed, such as the third-order shear deformation plate theory, ³² the sinusoidal shear deformation plate theory, ³³ the hyperbolic shear deformation plate theory, ³⁴ and the exponential shear deformation plate theory. ³⁵ The variation behavior of FG nanoplates when resting on an elastic foundation was studied by Bounouara et al. ³⁶ by using NET and the zeroth-order shear deformation plate theory. Chaht et al. ³⁷ investigated the static bending and buckling characteristics of nonlocal FG nanoplates, including the thickness stretching effect, based on a sinusoidal variation of all displacements through the plate thickness without using a shear correction factor. Using exponential shear deformation plate theory, Khorshidi and Fallah ³⁸ examined the size-dependent buckling behavior of nonlocal FG rectangular nanoplates. More recently, a unified size-dependent plate model was developed by Lu et al. ³⁹ to examine the buckling behavior of nanoplates while including the surface effects. The governing equations for their model were established based on NSGT and the HSDPTs.

Because the governing equations of structures are typically composed of a system of complex partial differential equations (PDEs), finding solutions to these governing equations often presents certain challenges. The finite element method (FEM) ³⁵ and various other numerical methods^28–30 provide effective tools to address these problems. To avoid using classical discretization methods such as the FEM entirely, and to realize extremely fast solutions after network training, deep neural networks (DNNs)^40–42 can be used to approximate the solutions to PDEs. The applicability of DNNs to solution of PDEs is broad and is not restricted to specific boundary conditions. Recently, based on the variational principle, a numerical method called the nonlocal operator method was developed to solve systems of PDEs.^43,44 The nonlocal operator can be regarded as an integral form that is ‘equivalent’ to the differential form in the context of a nonlocal interaction model that can be used to solve for the unknown field.^43,44 In addition to the FEM and numerical methods, analytical methods also have limitations because of their dependence on specific types of boundary conditions. However, they allow explicit solutions to be obtained, which is why analytical methods have been used widely by many researchers.^{18–23,26,27,37–39}

To the best of the authors’ knowledge, the mechanical behavior of FG microplates based on NSGT when used in conjunction with SET has not been studied to date. The purpose of this work is to develop a unified size-dependent plate model based on NSGT and the HSDPTs that can be used to examine the buckling responses of FG microplates while considering the surface energy effect and using an analytical approach for the first time. An FG microplate is considered in which the material properties are assumed to vary continuously in the microplate’s thickness direction based on a simple power law. The principle of minimum potential energy is used to derive both the equilibrium equations and the boundary conditions. Analytical solutions for the buckling problem are developed with various boundary conditions. The accuracy of the analytical results obtained has been verified via comparison with existing results in the literature. Additionally, several numerical illustrations are provided to investigate the influence of the nonlocal parameters, the material length scale parameter, shear deformation, aspect ratios, the power-law index, and the surface energy on the buckling responses of rectangular FG microplates.

Theoretical formulations

Functionally graded microplate

A rectangular FG microplate with length a, width b, and thickness h is considered, as illustrated in Figure 1. The FG microplate is subjected to the axially compressive loads $N_{\mathrm{xx}}^0$ and $N_{\mathrm{yy}}^0$ . The coordinate system is selected as shown in Figure 1, where the origin of the coordinate system is located at the corner of the middle plane of the FG microplate, and the x-, y- and z-axes are directed in the length, width, and thickness directions of the FG microplate, respectively. To enable consideration of the surface effect, the microplate is assumed to be composed of the bulk part and two surface layers located at $z=\pm h/2$ . Furthermore, the bulk part of the microplate is assumed to be composed of a mixture of ceramic and metal constituent materials, in which the material properties, including Young’s modulus $E(z)$ , the mass density $\rho (z)$ , and Poisson’s ratio $\nu (z)$ , are all considered to vary along the thickness direction of the microplate. In this work, the simple power law is used to estimate the material properties of the FG microplates:^37,38

$$E(z)=(E_c-E_m){(\frac{z}{h}+\frac{1}{2})}^k+E_m,$$

(1)

$$\rho (z)=({\rho }_c-{\rho }_m){(\frac{z}{h}+\frac{1}{2})}^k+{\rho }_m,$$

(2)

$$\nu \left(z\right)=({\nu }_c-{\nu }_m){(\frac{z}{h}+\frac{1}{2})}^k+{\nu }_m,$$

(3)

where the subscripts “c” and “m” represent the ceramic and metal phases, respectively; the exponent k ( $0\le k<\infty$ ) governs the variations in the volume fractions of the ceramic and metal phases, and it is known as the power-law index. It can be observed from these equations that $E(z)=E_c,\rho (z)={\rho }_c$ and $\nu (z)={\nu }_c$ if $k=0$ (which means that the bulk part of the microplate is composed entirely of ceramic). In contrast, if $k\to \infty$ , it leads to $E(z)=E_m,\rho (z)={\rho }_m$ and $\nu (z)={\nu }_m$ (which means that the bulk part of the microplate is composed entirely of metal).

Figure 1.

Configuration of a rectangular functionally graded (FG) microplate under in-plane axially compressive loading.

Nonlocal strain gradient theory

According to the NSGT, as proposed by Lim et al., ¹⁸ the total stress ( $t_{\mathrm{ij}}$ ) is a function of the classical stress ( ${\sigma }_{\mathrm{ij}}$ ) and the higher-order stress ( ${\sigma }_{\mathrm{ijm}}^{(1)}$ ), as shown below:

$$t_{\mathrm{ij}}={\sigma }_{\mathrm{ij}}-\nabla {\sigma }_{\mathrm{ijm}}^{(1)},$$

(4)

where $\nabla =\partial /\partial x+\partial /\partial y$ is the differential operator; and the classical and higher-order stresses are then defined as: ¹⁸

$${\sigma }_{\mathrm{ij}}={\int }_V{\alpha }_0(\mathit{x}\prime ,\mathit{x},e_{0}a)C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}\prime dV\prime ,$$

(5)

$${\sigma }_{\mathrm{ijm}}^{(1)}=l^2{\int }_V{\alpha }_1(\mathit{x}\prime ,\mathit{x},e_{1}a)C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl},m}\prime dV\prime ,$$

(6)

where ${\alpha }_0(\mathit{x}\prime ,\mathit{x},e_{0}a)$ and ${\alpha }_1(\mathit{x}\prime ,\mathit{x},e_{1}a)$ are two kernel functions that were defined by Eringen;^8,9 $e_{0}a$ and $e_{1}a$ represent two nonlocal parameters that are used to assess the influence of the nonlocal stress effect; the effect of the strain gradient can be captured by introducing the material length scale parameter l; $C_{\mathrm{ijkl}}$ represents the elastic coefficients; and ${\epsilon }_{\mathrm{kl}}$ and ${\epsilon }_{\mathrm{kl},m}$ are the strain and the strain gradient, respectively. The two kernel functions satisfy the following conditions ¹⁸ :

$$[1-{(e_{0}a)}^2{\nabla }^2]{\sigma }_{\mathrm{ij}}=C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}},$$

(7)

$$[1-{(e_{1}a)}^2{\nabla }^2]{\sigma }_{\mathrm{ijm}}^{(1)}=l^2{C}_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl},m},$$

(8)

in which ${\nabla }^2={\partial }^2/\partial x^2+{\partial }^2/\partial y^2$ denotes the Laplacian operator. The generalized constitutive equation can then be obtained by substituting equations (7) and (8) into equation (4), with results as shown below: ¹⁸

$$\begin{array}{c}[1-{(e_{1}a)}^2{\nabla }^2][1-{(e_{0}a)}^2{\nabla }^2]t_{\mathrm{ij}}\hfill \\ =[1-{(e_{1}a)}^2{\nabla }^2]C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}-l^2[1-{(e_{0}a)}^2{\nabla }^2]C_{\mathrm{ijkl}}{\nabla }^2{\epsilon }_{\mathrm{kl}}.\hfill \end{array}$$

(9)

For simplicity, it is assumed that $e_{0}a=e_{1}a=\mathrm{ea}$ . Therefore, the constitutive equation (9) can be reduced to:

$$[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{ij}}=(1-l^2{\nabla }^2)C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}.$$

(10)

The nonlocal strain gradient constitutive equation (10) proposed by Lim et al. ¹⁸ has been used extensively to investigate the mechanical responses of small-scale structures while considering the combined effects of the nonlocal stress and the strain gradient. It should be noted that by setting l = 0 in equation (10), the constitutive equation based on NET^8,9 can be recovered, while setting ea = 0 in equation (10) recovers the constitutive equation based on SGT.^10,11

Surface elasticity theory

The effect of the cohesive forces acting between the atoms on the surfaces of micro/nanostructures can be captured using SET, which was proposed by Gurtin and Murdoch.^24,25 A surface is described as a flexible membrane with zero mathematical thickness that adheres completely to the underlying bulk material. Based on SET^24,25, the constitutive equations for the two surface layers that are underlying the bulk part of the FG microplate can be expressed as:

$$\begin{array}{c}{\sigma }_{\alpha \beta }^{s\pm }={\tau }^s{\delta }_{\alpha \beta }+({\tau }^s+{\lambda }^s){\epsilon }_{\mathrm{kk}}^{s\pm }{\delta }_{\alpha \beta }+2({\mu }^s-{\tau }^s){\epsilon }_{\alpha \beta }^{s\pm }+{\tau }^s{u}_{\alpha ,\beta }^{s\pm },\hfill \\ {\sigma }_{\alpha z}^{s\pm }={\tau }^s{u}_{z,\alpha }^{s\pm },(\alpha ,\beta =x,y),\hfill \end{array}$$

(11)

where ${\delta }_{\alpha \beta }$ is the Kronecker-delta function, ${\lambda }^s$ and ${\mu }^s$ denote Lamé’s constants for the surface layers, which can be determined via atomic simulations, ${\tau }^s$ is the residual surface stress, and $u_{\alpha ,\beta }^{s\pm }$ represents the displacement components of the surface layers. Therefore, the surface stress components can be expressed as follows:

$${\sigma }_{\mathrm{xx}}^{s\pm }={\tau }^s+({\lambda }^s+2{\mu }^s){\epsilon }_{\mathrm{xx}}^{s\pm }+({\lambda }^s+{\tau }^s){\epsilon }_{\mathrm{yy}}^{s\pm },$$

(12)

$${\sigma }_{\mathrm{yy}}^{s\pm }={\tau }^s+({\lambda }^s+2{\mu }^s){\epsilon }_{\mathrm{yy}}^{s\pm }+({\lambda }^s+{\tau }^s){\epsilon }_{\mathrm{xx}}^{s\pm },$$

(13)

$${\sigma }_{\mathrm{xy}}^{s\pm }=2({\mu }^s-{\tau }^s){\epsilon }_{\mathrm{xy}}^{s\pm }+{\tau }^s\frac{\partial u_x^{s\pm }}{\partial y},$$

(14)

$${\sigma }_{\mathrm{yx}}^{s\pm }=2({\mu }^s-{\tau }^s){\epsilon }_{\mathrm{xy}}^{s\pm }+{\tau }^s\frac{\partial u_y^{s\pm }}{\partial x},$$

(15)

$${\sigma }_{\mathrm{xz}}^{s\pm }={\tau }^s\frac{\partial w}{\partial x},$$

(16)

$${\sigma }_{\mathrm{yz}}^{s\pm }={\tau }^s\frac{\partial w}{\partial y}.$$

(17)

In classical elasticity theory, the normal stress component ( ${\sigma }_{\mathrm{zz}}$ ) of the bulk part is usually neglected because it is very small when compared with the other normal stresses. However, the surface equilibrium conditions in SET are not satisfied by this hypothesis. To address this problem, the normal stress component ( ${\sigma }_{\mathrm{zz}}$ ) is assumed to vary linearly through the thickness direction of the FG microplate, as follows: ³⁹

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{\mathrm{zz}}=\frac{1}{2}[(\frac{\partial {\sigma }_{\mathrm{xz}}^{s+}}{\partial x}+\frac{\partial {\sigma }_{\mathrm{yz}}^{s+}}{\partial y})-(\frac{\partial {\sigma }_{\mathrm{xz}}^{s-}}{\partial x}+\frac{\partial {\sigma }_{\mathrm{yz}}^{s-}}{\partial y})]}\\ +\frac{z}{h}[(\frac{\partial {\sigma }_{\mathrm{xz}}^{s+}}{\partial x}+\frac{\partial {\sigma }_{\mathrm{yz}}^{s+}}{\partial y})+(\frac{\partial {\sigma }_{\mathrm{xz}}^{s-}}{\partial x}+\frac{\partial {\sigma }_{\mathrm{yz}}^{s-}}{\partial y})].\hfill \end{array}$$

(18)

Combining equations (16), (17), and (18) leads to:

$${\sigma }_{\mathrm{zz}}=\frac{2{\tau }^{s}z}{h}(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}).$$

(19)

Inclusion of ${\sigma }_{\mathrm{zz}}$ in the normal stresses for the bulk part of the FG microplate allows the stress components to be rewritten as follows:

$${\sigma }_{\mathrm{xx}}=\frac{E}{1-{\nu }^2}({\epsilon }_{\mathrm{xx}}+\nu {\epsilon }_{\mathrm{yy}})+\frac{2\nu {\tau }^{s}z}{(1-\nu )h}(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}),$$

(20)

$${\sigma }_{\mathrm{yy}}=\frac{E}{1-{\nu }^2}({\epsilon }_{\mathrm{yy}}+\nu {\epsilon }_{\mathrm{xx}})+\frac{2\nu {\tau }^{s}z}{(1-\nu )h}(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}),$$

(21)

$${\sigma }_{\mathrm{xy}}=G{\gamma }_{\mathrm{xy}},{\sigma }_{\mathrm{xz}}=G{\gamma }_{\mathrm{xz}},{\sigma }_{\mathrm{yz}}=G{\gamma }_{\mathrm{yz}},$$

(22)

where $G=E/[2(1+\nu )]$ is the shear modulus.

Higher-order shear deformation plate theory

In this work, a unified higher-order shear deformation plate theory that contains various existing plate theories is applied. According to this theory, the displacement field ( $u_x,u_y,u_z$ ) can be expressed as: ³⁹

$$u_x(x,y,z)=-z\frac{\partial w(x,y)}{\partial x}+\Phi (z)[\frac{\partial w(x,y)}{\partial x}+{\phi }_x(x,y)],$$

(23)

$$u_y(x,y,z)=-z\frac{\partial w(x,y)}{\partial y}+\Phi (z)[\frac{\partial w(x,y)}{\partial y}+{\phi }_y(x,y)],$$

(24)

$$u_z(x,y,z)=w(x,y),$$

(25)

where $w(x,y)$ represents the transverse displacement of any point on the middle plane surface of the FG microplate, and ${\phi }_x(x,y)$ and ${\phi }_y(x,y)$ denote the rotation angles about the x- and y-axes, respectively. The function $\Phi (z)$ represents the shape function that describes the distributions of the transverse shear deformation and stress across the thickness direction of the FG microplate. The displacement fields of the different plate theories can be obtained by selecting the appropriate shape function, as follows: ³⁹

$$+\mathrm{Kirchhoff}\mathrm{plate}\mathrm{theory}\left(\mathrm{KPT}\right):\Phi (z)=0.$$

(26)

$$+\mathrm{Mindlin}\mathrm{plate}\mathrm{theory}\left(\mathrm{MPT}\right):\Phi (z)=z.$$

(27)

$$\begin{array}{c}\multicolumn{1}{c}{+\mathrm{Third}-\mathrm{order}\mathrm{shear}\mathrm{deformation}\mathrm{plate}}\\ \multicolumn{1}{c}{\mathrm{theory}\left(\mathrm{TSDPT}\right):\Phi (z)=z(1-\frac{4z^2}{3h^2}).}\end{array}$$

(28)

$$\begin{array}{c}+\mathrm{Sinusoidal}\mathrm{shear}\mathrm{deformable}\mathrm{plate}\hfill \\ \mathrm{theory}\left(\mathrm{SSDPT}\right):\Phi (z)=\frac{h}{\pi }\sin \left(\frac{\pi z}{h}\right).\hfill \end{array}$$

(29)

$$\begin{array}{c}+\mathrm{Hyperbolic}\mathrm{shear}\mathrm{deformable}\mathrm{plate}\hfill \\ \mathrm{theory}\left(\mathrm{HSDPT}\right):\Phi (z)=\mathrm{hsin}h\left(\frac{z}{h}\right)-\mathrm{zcos}h\left(\frac{1}{2}\right).\hfill \end{array}$$

(30)

$$\begin{array}{c}+\mathrm{Exponential}\mathrm{shear}\mathrm{deformable}\mathrm{plate}\hfill \\ \mathrm{theory}\left(\mathrm{ESDPT}\right):\Phi (z)=z{e}^{-2{\left(\frac{z}{h}\right)}^2}.\hfill \end{array}$$

(31)

Using the displacement field given in equations (23)–(25), the strain components of the FG microplate can be expressed as follows:

$${\epsilon }_{\mathrm{xx}}=\frac{\partial u_x}{\partial x}=(\Phi -z)\frac{{\partial }^{2}w}{\partial x^2}+\Phi \frac{\partial {\phi }_x}{\partial x},$$

(32)

$${\epsilon }_{\mathrm{yy}}=\frac{\partial u_y}{\partial y}=(\Phi -z)\frac{{\partial }^{2}w}{\partial y^2}+\Phi \frac{\partial {\phi }_y}{\partial y},$$

(33)

$${\gamma }_{\mathrm{xy}}=\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}=2(\Phi -z)\frac{{\partial }^{2}w}{\partial x\partial y}+\Phi (\frac{\partial {\phi }_x}{\partial y}+\frac{\partial {\phi }_y}{\partial x}),$$

(34)

$${\gamma }_{\mathrm{xz}}=\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}=\frac{\partial \Phi }{\partial z}(\frac{\partial w}{\partial x}+{\phi }_x),$$

(35)

$${\gamma }_{\mathrm{yz}}=\frac{\partial u_y}{\partial z}+\frac{\partial u_z}{\partial y}=\frac{\partial \Phi }{\partial z}(\frac{\partial w}{\partial y}+{\phi }_y).$$

(36)

The governing equation

Expressions for the stress components of the bulk part in terms of the displacement components can be obtained by substituting the strain components in equations (32)–(36) into the equations for the stress components given in equations (20)–(22); then, by using the relation given in equation (10), the constitutive equations for the bulk part of the FG microplate based on NSGT can be obtained as follows:

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{xx}}=(1-l^2{\nabla }^2)\hfill \\ \left\{\begin{array}{c}[\frac{E(\Phi -z)}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z}{(1-\nu )h}]\frac{{\partial }^{2}w}{\partial x^2}+[\frac{E\nu (\Phi -z)}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z}{(1-\nu )h}]\frac{{\partial }^{2}w}{\partial y^2}\hfill \\ +\frac{E\Phi }{1-{\nu }^2}\frac{\partial {\phi }_x}{\partial x}+\frac{E\nu \Phi }{1-{\nu }^2}\frac{\partial {\phi }_y}{\partial y}\hfill \end{array}\right\},\hfill \end{array}$$

(37)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{yy}}=(1-l^2{\nabla }^2)\hfill \\ \left\{\begin{array}{c}[\frac{E\nu (\Phi -z)}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z}{(1-\nu )h}]\frac{{\partial }^{2}w}{\partial x^2}+[\frac{E(\Phi -z)}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z}{(1-\nu )h}]\frac{{\partial }^{2}w}{\partial y^2}\hfill \\ +\frac{E\nu \Phi }{1-{\nu }^2}\frac{\partial {\phi }_x}{\partial x}+\frac{E\Phi }{1-{\nu }^2}\frac{\partial {\phi }_y}{\partial y}\hfill \end{array}\right\},\hfill \end{array}$$

(38)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{xy}}=G(1-l^2{\nabla }^2)\hfill \\ [2(\Phi -z)\frac{{\partial }^{2}w}{\partial x\partial y}+R(\frac{\partial {\phi }_y}{\partial x}+\frac{\partial {\phi }_x}{\partial y})],\hfill \end{array}$$

(39)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{xz}}=G\frac{d\Phi }{\mathrm{dz}}(1-l^2{\nabla }^2)(\frac{\partial w}{\partial x}+{\phi }_x),$$

(40)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{yz}}=G\frac{d\Phi }{\mathrm{dz}}(1-l^2{\nabla }^2)(\frac{\partial w}{\partial y}+{\phi }_y).$$

(41)

In addition, the constitutive equations for the surface layers based on NSGT can be obtained by inserting equations (32)–(36) into equations (12)–(17) and then using equation (10), as shown below:

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{xx}}^{s\pm }=(1-l^2{\nabla }^2)\hfill \\ \left\{\begin{array}{c}{\tau }^s+({\lambda }^s+2{\mu }^s)\{[\Phi (\pm \frac{h}{2})\mp \frac{h}{2}]\frac{{\partial }^{2}w}{\partial x^2}+\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_x}{\partial x}\}\hfill \\ +({\lambda }^s+{\tau }^s)\{[\Phi (\pm \frac{h}{2})\mp \frac{h}{2}]\frac{{\partial }^{2}w}{\partial y^2}+\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_y}{\partial y}\}\hfill \end{array}\right\},\hfill \end{array}$$

(42)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{yy}}^{s\pm }=(1-l^2{\nabla }^2)\hfill \\ \left\{\begin{array}{c}{\tau }^s+({\lambda }^s+2{\mu }^s)\{[\Phi (\pm \frac{h}{2})\mp \frac{h}{2}]\frac{{\partial }^{2}w}{\partial y^2}+\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_y}{\partial y}\}\hfill \\ +({\lambda }^s+{\tau }^s)\{[\Phi (\pm \frac{h}{2})\mp \frac{h}{2}]\frac{{\partial }^{2}w}{\partial x^2}+\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_x}{\partial x}\}\hfill \end{array}\right\},\hfill \end{array}$$

(43)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{xy}}^{s\pm }=(1-l^2{\nabla }^2)\hfill \\ \left\{\begin{array}{c}(2{\mu }^s-{\tau }^s)[\Phi (\pm \frac{h}{2})\mp \frac{h}{2}]\frac{{\partial }^{2}w}{\partial x\partial y}+({\mu }^s-{\tau }^s)\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_y}{\partial x}\hfill \\ +{\mu }^s\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_x}{\partial y}\hfill \end{array}\right\},\hfill \end{array}$$

(44)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{yx}}^{s\pm }=(1-l^2{\nabla }^2)\hfill \\ \left\{\begin{array}{c}(2{\mu }^s-{\tau }^s)[\Phi (\pm \frac{h}{2})\mp \frac{h}{2}]\frac{{\partial }^{2}w}{\partial x\partial y}+({\mu }^s-{\tau }^s)\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_x}{\partial y}\hfill \\ +{\mu }^s\Phi (\pm \frac{h}{2})\frac{\partial {\phi }_y}{\partial x}\hfill \end{array}\right\},\hfill \end{array}$$

(45)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{xz}}^{s\pm }={\tau }^s(1-l^2{\nabla }^2)\frac{\partial w}{\partial x},$$

(46)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]t_{\mathrm{yz}}^{s\pm }={\tau }^s(1-l^2{\nabla }^2)\frac{\partial w}{\partial y}.$$

(47)

To derive both the equilibrium equations and the boundary conditions, the principle of the minimum total potential energy is used, and this principle can be stated as:

$$\delta U=\delta W,$$

(48)

where $\delta U$ represents the variation in the strain energy of the FG microplate and $\delta W$ is the variation in the work done by the externally applied loads. By incorporating the surface energy effect, the variation in the strain energy of the FG microplate can be determined as follows:

$$\delta U={\int }_A{\int }_{-h/2}^{h/2}{t}_{ij}\delta {\epsilon }_{ij}dzdA+{\int }_{S^{+}}{t}_{ij}^{s+}\delta {\epsilon }_{ij}d{S}^{+}+{\int }_{S^{-}}{t}_{ij}^{s-}\delta {\epsilon }_{ij}d{S}^{-}={\int }_A\left[\begin{array}{l}{M}_{xx}\frac{{\partial }^2\delta w}{\partial x^2}+M_{yy}\frac{{\partial }^2\delta w}{\partial y^2}+2M_{xy}\frac{{\partial }^2\delta w}{\partial x\partial y}+P_{xx}\frac{\partial \delta {\phi }_x}{\partial x}\\ +P_{yy}\frac{\partial \delta {\phi }_y}{\partial y}+P_{xy}\left(\frac{\partial \delta {\phi }_y}{\partial x}+\frac{\partial \delta {\phi }_x}{\partial y}\right)+Q_x\left(\frac{\partial \delta w}{\partial x}+\delta {\phi }_x\right)\\ +Q_y\left(\frac{\partial \delta w}{\partial y}+\delta {\phi }_y\right)+Q_x^s\frac{\partial \delta w}{\partial x}+Q_y^s\frac{\partial \delta w}{\partial y}\end{array}\right]dA,$$

(49)

where the following stress resultants are defined as:

$$M_{\mathrm{xx}}={\int }_{-h/2}^{h/2}{t}_{\mathrm{xx}}(\Phi -z)\mathrm{dz}+[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}](t_{\mathrm{xx}}^{s+}-t_{\mathrm{xx}}^{s-}),$$

(50)

$$M_{\mathrm{yy}}={\int }_{-h/2}^{h/2}{t}_{\mathrm{yy}}(\Phi -z)\mathrm{dz}+[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}](t_{\mathrm{yy}}^{s+}-t_{\mathrm{yy}}^{s-}),$$

(51)

$$M_{\mathrm{xy}}={\int }_{-h/2}^{h/2}{t}_{\mathrm{xy}}(\Phi -z)\mathrm{dz}+\frac{1}{2}[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}](t_{\mathrm{xy}}^{s+}+t_{\mathrm{yx}}^{s+}-t_{\mathrm{xy}}^{s-}-t_{\mathrm{yx}}^{s-}),$$

(52)

$$P_{\mathrm{xx}}={\int }_{-h/2}^{h/2}{t}_{\mathrm{xx}}(\Phi -z)\mathrm{dz}+\Phi \left(\frac{h}{2}\right)(t_{\mathrm{xx}}^{s+}-t_{\mathrm{xx}}^{s-}),$$

(53)

$$P_{\mathrm{yy}}={\int }_{-h/2}^{h/2}{t}_{\mathrm{yy}}(\Phi -z)\mathrm{dz}+\Phi \left(\frac{h}{2}\right)(t_{\mathrm{yy}}^{s+}-t_{\mathrm{yy}}^{s-}),$$

(54)

$$P_{\mathrm{xy}}={\int }_{-h/2}^{h/2}{t}_{\mathrm{xy}}\Phi \mathrm{dz}+\frac{1}{2}\Phi \left(\frac{h}{2}\right)(t_{\mathrm{xy}}^{s+}+t_{\mathrm{yx}}^{s+}-t_{\mathrm{xy}}^{s-}-t_{\mathrm{yx}}^{s-}),$$

(55)

$$Q_x=k_s{\int }_{-h/2}^{h/2}{t}_{\mathrm{xz}}\frac{d\Phi }{\mathrm{dz}}\mathrm{dz},$$

(56)

$$Q_y=k_s{\int }_{-h/2}^{h/2}{t}_{\mathrm{yz}}\frac{d\Phi }{\mathrm{dz}}\mathrm{dz},$$

(57)

$$Q_x^s=t_{\mathrm{xy}}^{s+}+t_{\mathrm{xy}}^{s-},$$

(58)

$$Q_y^s=t_{\mathrm{yz}}^{s+}+t_{\mathrm{yz}}^{s-},$$

(59)

in which $k_s$ represents the shear correction factor, which is selected to be 5/6 for the Mindlin microplate and to be 1 for the rest. The variation of work done by external applied loads is determined as:

$$\delta W={\int }_A(N_{\mathrm{xx}}^0\frac{\partial w}{\partial x}\frac{\partial \delta w}{\partial x}+N_{\mathrm{yy}}^0\frac{\partial w}{\partial y}\frac{\partial \delta w}{\partial y})\mathrm{dA}.$$

(60)

Using equations (49) and (60), equation (48) leads to:

$$\frac{\partial P_{\mathrm{xx}}}{\partial x}+\frac{\partial P_{\mathrm{xy}}}{\partial y}-Q_x=0,$$

(61)

$$\frac{\partial P_{\mathrm{yy}}}{\partial y}+\frac{\partial P_{\mathrm{xy}}}{\partial x}-Q_y=0,$$

(62)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2{M}_{\mathrm{xx}}}{\partial x^2}+\frac{{\partial }^2{M}_{\mathrm{yy}}}{\partial y^2}+2\frac{{\partial }^2{M}_{\mathrm{xy}}}{\partial x\partial y}-\frac{\partial Q_x}{\partial x}-\frac{\partial Q_y}{\partial y}-\frac{\partial Q_x^s}{\partial x}-\frac{\partial Q_y^s}{\partial y}}\\ +N_{\mathrm{xx}}^0\frac{{\partial }^{2}w}{\partial x^2}+N_{\mathrm{yy}}^0\frac{{\partial }^{2}w}{\partial y^2}=0,\hfill \end{array}$$

(63)

and the corresponding boundary conditions:

$$\delta {\phi }_x=0\mathrm{or}{P}_{\mathrm{xx}}{n}_x+P_{\mathrm{xy}}{n}_y=0,$$

(64)

$$\delta {\phi }_y=0\mathrm{or}{P}_{\mathrm{xy}}{n}_x+P_{\mathrm{yy}}{n}_y=0,$$

(65)

$$\frac{\partial \delta w}{\partial x}=0\mathrm{or}{M}_{\mathrm{xx}}{n}_x+M_{\mathrm{xy}}{n}_y=0,$$

(66)

$$\frac{\partial \delta w}{\partial y}=0\mathrm{or}{M}_{\mathrm{xy}}{n}_x+M_{\mathrm{yy}}{n}_y=0,$$

(67)

$$\begin{array}{c}\multicolumn{1}{c}{\delta w=0\mathrm{or}(\frac{\partial M_{\mathrm{xx}}}{\partial x}+\frac{\partial M_{\mathrm{xy}}}{\partial y}+N_{\mathrm{xx}}^0\frac{\partial w}{\partial x}-Q_x-Q_x^s)n_x}\\ +(\frac{\partial M_{\mathrm{yy}}}{\partial y}+\frac{\partial M_{\mathrm{xy}}}{\partial x}+N_{\mathrm{yy}}^0\frac{\partial w}{\partial y}-Q_y-Q_y^s)n_y=0,\hfill \end{array}$$

(68)

where $n_x=\cos \theta$ and $n_y=\sin \theta$ are the components of the unit normal vector $n$ when oriented clockwise at an angle of θ from the positive x-axis.

The expressions given for the stress resultants in terms of the displacement components can be obtained by substituting equations (37)–(47) into equations (50)–(59), as shown below:

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]M_{\mathrm{xx}}=(1-l^2{\nabla }^2)\hfill \\ (A_{11}\frac{{\partial }^{2}w}{\partial x^2}+B_{11}\frac{{\partial }^{2}w}{\partial y^2}+C_{11}\frac{\partial {\phi }_x}{\partial x}+D_{11}\frac{\partial {\phi }_y}{\partial y}),\hfill \end{array}$$

(69)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]M_{\mathrm{yy}}=(1-l^2{\nabla }^2)\hfill \\ (B_{11}\frac{{\partial }^{2}w}{\partial x^2}+A_{11}\frac{{\partial }^{2}w}{\partial y^2}+D_{11}\frac{\partial {\phi }_x}{\partial x}+C_{11}\frac{\partial {\phi }_y}{\partial y}),\hfill \end{array}$$

(70)

$$\begin{array}{c}\multicolumn{1}{c}{[1-{(\mathrm{ea})}^2{\nabla }^2]M_{\mathrm{xy}}=(1-l^2{\nabla }^2)[E_{11}\frac{{\partial }^{2}w}{\partial x\partial y}+F_{11}(\frac{\partial {\phi }_x}{\partial y}+\frac{\partial {\phi }_y}{\partial x})],}\end{array}$$

(71)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]P_{\mathrm{xx}}=(1-l^2{\nabla }^2)\hfill \\ [G_{11}\frac{{\partial }^{2}w}{\partial x^2}+H_{11}\frac{{\partial }^{2}w}{\partial y^2}+I_{11}\frac{\partial {\phi }_x}{\partial x}+J_{11}\frac{\partial {\phi }_y}{\partial y}],\hfill \end{array}$$

(72)

$$\begin{array}{c}[1-{(\mathrm{ea})}^2{\nabla }^2]P_{\mathrm{yy}}=(1-l^2{\nabla }^2)\hfill \\ [H_{11}\frac{{\partial }^{2}w}{\partial x^2}+G_{11}\frac{{\partial }^{2}w}{\partial y^2}+J_{11}\frac{\partial {\phi }_x}{\partial x}+I_{11}\frac{\partial {\phi }_y}{\partial y}],\hfill \end{array}$$

(73)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]P_{\mathrm{xy}}=(1-l^2{\nabla }^2)[K_{11}\frac{{\partial }^{2}w}{\partial x\partial y}+L_{11}(\frac{\partial {\phi }_y}{\partial x}+\frac{\partial {\phi }_x}{\partial y})],$$

(74)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]Q_x=S_{11}(1-l^2{\nabla }^2)[\frac{\partial w}{\partial x}+{\phi }_x],$$

(75)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]Q_y=S_{11}(1-l^2{\nabla }^2)[\frac{\partial w}{\partial y}+{\phi }_y],$$

(76)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]Q_x^s=2{\tau }^s(1-l^2{\nabla }^2)\frac{\partial w}{\partial x},$$

(77)

$$[1-{(\mathrm{ea})}^2{\nabla }^2]Q_y^s=2{\tau }^s(1-l^2{\nabla }^2)\frac{\partial w}{\partial y},$$

(78)

where:

$$\begin{array}{c}\multicolumn{1}{c}{A_{11}={\int }_{-h/2}^{h/2}[\frac{E{(\Phi -z)}^2}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z(\Phi -z)}{(1-\nu )h}]\mathrm{dz}}\\ +2({\lambda }^s+2{\mu }^s){[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}]}^2,\hfill \end{array}$$

(79)

$$\begin{array}{c}\multicolumn{1}{c}{B_{11}={\int }_{-h/2}^{h/2}[\frac{E\nu {(\Phi -z)}^2}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z(\Phi -z)}{(1-\nu )h}]\mathrm{dz}}\\ +2({\lambda }^s+{\tau }^s){[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}]}^2,\hfill \end{array}$$

(80)

$$C_{11}={\int }_{-h/2}^{h/2}\left[\frac{E\Phi (\Phi -z)}{1-{\nu }^2}\right]\mathrm{dz}+2({\lambda }^s+2{\mu }^s)\Phi \left(\frac{h}{2}\right)[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}],$$

(81)

$$D_{11}={\int }_{-h/2}^{h/2}\left[\frac{E\nu \Phi (\Phi -z)}{1-{\nu }^2}\right]\mathrm{dz}+2({\lambda }^s+{\tau }^s)\Phi \left(\frac{h}{2}\right)[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}],$$

(82)

$$E_{11}={\int }_{-h/2}^{h/2}2G(\Phi -z{)}^2\mathrm{dz}+2(2{\mu }^s-{\tau }^s){[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}]}^2,$$

(83)

$$F_{11}={\int }_{-h/2}^{h/2}G\Phi (\Phi -z)\mathrm{dz}+(2{\mu }^s-{\tau }^s)\Phi \left(\frac{h}{2}\right)[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}],$$

(84)

$$\begin{array}{c}\multicolumn{1}{c}{G_{11}={\int }_{-h/2}^{h/2}[\frac{E\Phi (\Phi -z)}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z\Phi }{(1-\nu )h}]\mathrm{dz}}\\ +2({\lambda }^s+2{\mu }^s)\Phi \left(\frac{h}{2}\right)[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}],\hfill \end{array}$$

(85)

$$\begin{array}{c}\multicolumn{1}{c}{H_{11}={\int }_{-h/2}^{h/2}[\frac{E\nu \Phi (\Phi -z)}{1-{\nu }^2}+\frac{2\nu {\tau }^{s}z\Phi }{(1-\nu )h}]\mathrm{dz}}\\ +2({\lambda }^s+{\tau }^s)\Phi \left(\frac{h}{2}\right)[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}],\hfill \end{array}$$

(86)

$$I_{11}={\int }_{-h/2}^{h/2}\left[\frac{E{\Phi }^2}{1-{\nu }^2}\right]\mathrm{dz}+2({\lambda }^s+2{\mu }^s){\Phi }^2\left(\frac{h}{2}\right),$$

(87)

$$J_{11}={\int }_{-h/2}^{h/2}\left[\frac{E\nu {\Phi }^2}{1-{\nu }^2}\right]\mathrm{dz}+2({\lambda }^s+{\tau }^s){\Phi }^2\left(\frac{h}{2}\right),$$

(88)

$$K_{11}={\int }_{-h/2}^{h/2}2G\Phi (\Phi -z)\mathrm{dz}+2(2{\mu }^s-{\tau }^s)\Phi \left(\frac{h}{2}\right)[\Phi \left(\frac{h}{2}\right)-\frac{h}{2}],$$

(89)

$$L_{11}={\int }_{-h/2}^{h/2}G{\Phi }^2\mathrm{dz}+(2{\mu }^s-{\tau }^s){\Phi }^2\left(\frac{h}{2}\right),$$

(90)

$$S_{11}={\int }_{-h/2}^{h/2}{k}_{s}G{\left(\frac{d\Phi }{\mathrm{dz}}\right)}^2\mathrm{dz}.$$

(91)

To obtain the governing equations in terms of the displacements for the FG microplate based on NSGT, the expressions for the stress resultants given in equations (69)–(78) are substituted into Equations (61)–(63), as shown below:

$$\begin{array}{c}{G}_{11}\frac{{\partial }^{3}w}{\partial x^3}+(H_{11}+K_{11})\frac{{\partial }^{3}w}{\partial x\partial y^2}+I_{11}\frac{{\partial }^2{\phi }_x}{\partial x^2}+(J_{11}+L_{11})\frac{{\partial }^2{\phi }_y}{\partial x\partial y}\hfill \\ +L_{11}\frac{{\partial }^2{\phi }_x}{\partial y^2}-S_{11}(\frac{\partial w}{\partial x}+{\phi }_x)=0,\hfill \end{array}$$

(92)

$$\begin{array}{c}{G}_{11}\frac{{\partial }^{3}w}{\partial y^3}+(H_{11}+K_{11})\frac{{\partial }^{3}w}{\partial x^2\partial y}+L_{11}\frac{{\partial }^2{\phi }_y}{\partial x^2}+(J_{11}+L_{11})\frac{{\partial }^2{\phi }_x}{\partial x\partial y}\hfill \\ +I_{11}\frac{{\partial }^2{\phi }_y}{\partial y^2}-S_{11}(\frac{\partial w}{\partial y}+{\phi }_y)=0,\hfill \end{array}$$

(93)

$$\begin{array}{l}(1-l^2{\nabla }^2)\left[\begin{array}{l}{A}_{11}\frac{{\partial }^{4}w}{\partial x^4}+2(B_{11}+E_{11})\frac{{\partial }^{4}w}{\partial x^2\partial y^2}+A_{11}\frac{{\partial }^{4}w}{\partial y^4}\\ +C_{11}\left(\frac{{\partial }^3{\phi }_x}{\partial x^3}+\frac{{\partial }^3{\phi }_y}{\partial y^3}\right)+(D_{11}+2F_{11})\left(\frac{{\partial }^3{\phi }_y}{\partial x^2\partial y}+\frac{{\partial }^3{\phi }_x}{\partial x\partial y^2}\right)\\ -S_{11}\left(\frac{\partial {\phi }_x}{\partial x}+\frac{\partial {\phi }_y}{\partial y}\right)-(2{\tau }^s+S_{11})\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}\right)\end{array}\right]\\ +\left[1-{(ea)}^2{\nabla }^2\right]\left(N_{xx}^0\frac{{\partial }^{2}w}{\partial x^2}+N_{yy}^0\frac{{\partial }^{2}w}{\partial y^2}\right)=0.\end{array}$$

(94)

In this work, the rectangular FG microplate has several boundary conditions to be considered, including fully simply supported (SSSS), fully clamped (CCCC), or clamped at x = 0 and a but simply supported at y = 0 and b (CCSS). The mathematical representations of the boundary conditions are thus given as follows:

For the SSSS FG microplate:

$$\begin{array}{c}w=0,{\phi }_x=0,M_{\mathrm{xx}}=0\mathrm{at}x=0\mathrm{and}x=a,w=0,\hfill \\ {\phi }_y=0,M_{\mathrm{yy}}=0\mathrm{at}y=0\mathrm{and}y=b.\hfill \end{array}$$

(95)

For the CCCC FG microplate:

$$\begin{array}{c}w=0,{\phi }_x={\phi }_y=0\mathrm{at}x=0\mathrm{and}x=a,w=0,\hfill \\ {\phi }_x={\phi }_y=0\mathrm{at}y=0\mathrm{and}y=b.\hfill \end{array}$$

(96)

For the CCSS FG microplate:

$$\begin{array}{c}\mathit{w}=0,{\mathbf{\phi }}_{\mathit{x}}={\mathbf{\phi }}_{\mathit{y}}=0\mathrm{at}\mathit{x}=0\mathrm{and}\mathit{x}=\mathit{a},\hfill \\ \multicolumn{1}{c}{\mathit{w}=0,{\mathbf{\phi }}_{\mathit{y}}=0,{\mathit{M}}_{\mathit{yy}}=0\mathrm{at}\mathit{y}=0\mathrm{and}\mathit{y}=\mathit{b}.}\end{array}$$

(97)

Analytical solutions

In this section, the analytical solutions for the buckling analysis of the FG microplate with several boundary conditions will be determined. For the purposes of the in-plane biaxial buckling analysis, it is assumed that $N_{\mathrm{xx}}^0=N_{\mathrm{yy}}^0=N_{\mathrm{cr}}$ . When using the three types of boundary conditions described above, the displacement components are selected as follows: ³⁹

$${\phi }_x(x,y)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{\Theta }_{\mathrm{xmn}}\frac{d{X}_m(x)}{\mathrm{dx}}{Y}_n(y),$$

(98)

$${\phi }_y(x,y)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{\Theta }_{\mathrm{ymn}}{X}_m(x)\frac{d{Y}_n(y)}{\mathrm{dy}},$$

(99)

$$w(x,y)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{W}_{\mathrm{mn}}{X}_m(x)Y_n(y),$$

(100)

where ${\Theta }_{\mathrm{xmn}},{\Theta }_{\mathrm{ymn}}$ , and $W_{\mathrm{mn}}$ are arbitrary parameters that must be found, and m and n are the numbers of half-waves in the x and y directions, respectively. The mode functions $X_m(x)$ and $Y_n(y)$ are selected based on each of the types of boundary conditions for the FG microplate. ³⁹

For the SSSS FG microplate:

$$X_m(x)=\sin ({\alpha }_{m}x);Y_n(y)=\sin ({\beta }_{n}y),{\alpha }_m=m\pi /a;{\beta }_n=n\pi /b.$$

(101)

For the CCCC FG microplate:

$$\begin{array}{c}\multicolumn{1}{c}{X_m(x)=\sin ({\alpha }_{m}x)-\sinh ({\alpha }_{m}x)}\\ -\frac{\sin ({\alpha }_{m}a)-\sinh ({\alpha }_{m}a)}{\cos ({\alpha }_{m}a)-\cosh ({\alpha }_{m}a)}[\cos ({\alpha }_{m}x)-\cosh ({\alpha }_{m}x)],\hfill \\ \multicolumn{1}{c}{Y_n(y)=\sin ({\beta }_{n}y)-\sinh ({\beta }_{n}y)}\\ -\frac{\sin ({\beta }_{n}b)-\sinh ({\beta }_{n}b)}{\cos ({\beta }_{n}b)-\cosh ({\beta }_{n}b)}[\cos ({\beta }_{n}y)-\cosh ({\beta }_{n}y)],\hfill \end{array}$$

$$\begin{array}{c}{\alpha }_1=4.73/a,{\alpha }_2=7.853/a,{\alpha }_3=10.996/a,\hfill \\ {\alpha }_4=14.137/a\mathrm{and}{\alpha }_m=(m+0.5)\pi /a(m\ge 5),\hfill \\ {\beta }_1=4.73/b,{\beta }_2=7.853/b,{\beta }_3=10.996/b,\hfill \\ {\beta }_4=14.137/b\mathrm{and}{\beta }_n=(n+0.5)\pi /b(n\ge 5).\hfill \end{array}$$

(102)

For the CCSS FG microplate:

$$\begin{array}{c}\multicolumn{1}{c}{X_m(x)=\sin ({\alpha }_{m}x)-\sin h({\alpha }_{m}x)}\\ -\frac{\sin ({\alpha }_{m}a)-\sin h({\alpha }_{m}a)}{\cos ({\alpha }_{m}a)-\cos h({\alpha }_{m}a)}[\cos ({\alpha }_{m}x)-\cos h({\alpha }_{m}x)],\hfill \\ Y_m(y)=\sin ({\beta }_{n}y).\hfill \end{array}$$

(103)

Substitution of the chosen solutions given in equations (98)–(100) into the governing equations (92)–(94), multiplying each equation by its corresponding eigenfunction, and then integrating the resulting equations over the domain leads to the algebraic equation below:

$$\begin{array}{c}\left[\begin{array}{ccc}{m}_{11}\hfill & m_{12}\hfill & m_{13}\hfill \\ \multicolumn{1}{c}{m_{21}}& m_{22}\hfill & m_{23}\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{mn}}{m}_{31}}& {\gamma }_{\mathrm{mn}}{m}_{32}\hfill & {\gamma }_{\mathrm{mn}}{m}_{33}+{\chi }_{\mathrm{mn}}{N}_{\mathrm{cr}}(c_1+c_2)\hfill \end{array}\right]\hfill \\ \left[\begin{array}{c}{\Theta }_{\mathrm{xmn}}\hfill \\ \multicolumn{1}{c}{{\Theta }_{\mathrm{ymn}}}\\ \multicolumn{1}{c}{W_{\mathrm{mn}}}\end{array}\right]=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right],\hfill \end{array}$$

(104)

where:

$$\begin{array}{cc}\multicolumn{1}{c}{m_{11}={\int }_0^a{\int }_0^b}& (I_{11}\frac{d^3{X}_m}{d{x}^3}{Y}_n+L_{11}\frac{d{X}_m}{\mathrm{dx}}\frac{d^2{Y}_n}{d{y}^2}-S_{11}\frac{d{X}_m}{\mathrm{dx}}{Y}_n)\hfill \\ \frac{d{X}_m}{\mathrm{dx}}{Y}_n\mathrm{dxdy},\hfill \end{array}$$

(105)

$$m_{12}={\int }_0^a{\int }_0^b(J_{11}+L_{11})\frac{d{X}_m}{\mathrm{dx}}\frac{d^2{Y}_n}{d{y}^2}\frac{d{X}_m}{\mathrm{dx}}{Y}_n\mathrm{dxdy},$$

(106)

$$\begin{array}{cc}\multicolumn{1}{c}{m_{13}={\int }_0^a{\int }_0^b}& [G_{11}\frac{d^3{X}_m}{d{x}^3}{Y}_n+(H_{11}+K_{11})\frac{d{X}_m}{\mathrm{dx}}\frac{d^2{Y}_n}{d{y}^2}-S_{11}\frac{d{X}_m}{\mathrm{dx}}{Y}_n]\hfill \\ \frac{d{X}_m}{\mathrm{dx}}{Y}_n\mathrm{dxdy},\hfill \end{array}$$

(107)

$$m_{21}={\int }_0^a{\int }_0^b(J_{11}+L_{11})\frac{d^2{X}_m}{d{x}^2}\frac{d{Y}_n}{\mathrm{dy}}{X}_m\frac{d{Y}_n}{\mathrm{dy}}\mathrm{dxdy},$$

(108)

$$\begin{array}{c}\multicolumn{1}{c}{m_{22}={\int }_0^a{\int }_0^b(I_{11}{X}_m\frac{d^3{Y}_n}{d{y}^3}+L_{11}\frac{d^2{X}_m}{d{x}^2}\frac{d{Y}_n}{\mathrm{dy}}-S_{11}{X}_m\frac{d{Y}_n}{\mathrm{dy}})}\\ X_m\frac{d{Y}_n}{\mathrm{dy}}\mathrm{dxdy},\hfill \end{array}$$

(109)

$$\begin{array}{c}\multicolumn{1}{c}{m_{23}={\int }_0^a{\int }_0^b(G_{11}{X}_m\frac{d^3{Y}_n}{d{y}^3}+(H_{11}+K_{11})\frac{d^2{X}_m}{d{x}^2}\frac{d{Y}_n}{\mathrm{dy}}-S_{11}{X}_m\frac{d{Y}_n}{\mathrm{dy}})}\\ X_m\frac{d{Y}_n}{\mathrm{dy}}\mathrm{dxdy},\hfill \end{array}$$

(110)

$$\begin{array}{c}\multicolumn{1}{c}{m_{31}={\int }_0^a{\int }_0^b(C_{11}\frac{d^4{X}_m}{d{x}^4}{Y}_m+(D_{11}+2F_{11})\frac{d^2{X}_m}{d{x}^2}\frac{d^2{Y}_n}{d{y}^2}-S_{11}\frac{d^2{X}_m}{d{x}^2}{Y}_n)}\\ X_m{Y}_n\mathrm{dxdy},\hfill \end{array}$$

(111)

$$\begin{array}{c}\multicolumn{1}{c}{m_{32}={\int }_0^a{\int }_0^b(C_{11}{X}_m\frac{d^4{Y}_m}{d{y}^4}+(D_{11}+2F_{11})\frac{d^2{X}_m}{d{x}^2}\frac{d^2{Y}_n}{d{y}^2}-S_{11}{X}_m\frac{d^2{Y}_n}{d{y}^2})}\\ X_m{Y}_n\mathrm{dxdy},\hfill \end{array}$$

(112)

$$\begin{array}{c}\multicolumn{1}{c}{m_{33}={\int }_0^a{\int }_0^b\left[\begin{array}{c}{A}_{11}\frac{d^4{X}_m}{d{x}^4}{Y}_n+2(B_{11}+E_{11})\frac{d^2{X}_m}{d{x}^2}\frac{d^2{Y}_n}{d{y}^2}+A_{11}{X}_m\frac{d^4{Y}_n}{d{y}^4}\hfill \\ -(2{\tau }^s+S_{11})(X_m\frac{d^2{Y}_n}{d{y}^2}+\frac{d^2{X}_m}{d{x}^2}{Y}_n)\hfill \end{array}\right]}\\ X_m{Y}_n\mathrm{dxdy},\hfill \end{array}$$

(113)

$${\gamma }_{\mathrm{mn}}=1+({\alpha }_m^2+{\beta }_n^2)l^2;{\chi }_{\mathrm{mn}}=1+({\alpha }_m^2+{\beta }_n^2)(\mathrm{ea}{)}^2,$$

(114)

$$c_1={\int }_0^a{\int }_0^b\left(\frac{d^2{X}_m}{d{x}^2}{Y}_n\right)X_m{Y}_n\mathrm{dxdy};c_2={\int }_0^a{\int }_0^b\left(X_m\frac{d^2{Y}_n}{d{y}^2}\right)X_m{Y}_n\mathrm{dxdy}.$$

(115)

The critical buckling force $N_{\mathrm{cr}}$ can then be realized by setting the determinant of the coefficient matrix to zero, as shown below:

$$N_{\mathrm{cr}}=\frac{\begin{array}{c}(m_{11}{m}_{22}{m}_{33}+m_{12}{m}_{23}{m}_{31}+m_{21}{m}_{32}{m}_{13}-m_{13}{m}_{22}{m}_{31}\hfill \\ -m_{11}{m}_{23}{m}_{32}-m_{12}{m}_{21}{m}_{33}){\gamma }_{\mathrm{mn}}\hfill \end{array}}{(m_{12}{m}_{21}-m_{11}{m}_{22})(c_1+c_2){\chi }_{\mathrm{mn}}}.$$

(116)

Numerical results

In this section, several numerical illustrations are provided to verify the reliability of the results obtained and evaluate the influence of the important parameters on the buckling responses of the FG microplates based on NSGT and the HSDPTs. The constituents of the FG microplate are selected as follows: ³⁸ ceramic (Al₂O₃: $E_c=380\mathrm{GPa},{\rho }_c=3800\mathrm{kg}/m^3,{\nu }_c=0.3$ ) and metal (Al: $E_m=70\mathrm{GPa},{\rho }_m=2702\mathrm{kg}/m^3,{\nu }_c=0.3$ ). The surface layers to be considered are composed of silicon and have the following material properties: ³⁹ ${\mu }^s=-2.774N/m,{\lambda }^s=-4.488N/m$ , and ${\tau }^2=0.605N/m$ . Additionally, the dimensionless buckling force is defined as:

$${\overline{N}}_{\mathrm{cr}}=\frac{100}{E_{m}h}{N}_{\mathrm{cr}}.$$

(117)

Comparison study

The first comparison is performed for the nonlocal SSSS FG nanoplate when it is subjected to in-plane biaxial compression without surface energy effects. The dimensionless critical buckling forces ( $N_{\mathrm{cr}}^{*}=\frac{N_{\mathrm{cr}}{a}^2}{E_m{h}^3}$ ) obtained in this work and those realized by Khorshidi and Fallah ³⁸ are compared in Table 1. These dimensionless critical buckling forces were obtained based on ESDPT and NET. Note that the results for the nonlocal FG nanoplate can be determined from the present results by simply letting the material length scale parameter (l) be equal to zero ( $l=0$ ). Close agreement can be observed between the present results and the results of Khorshidi and Fallah. ³⁸ There is a slight difference between the present results and those of Khorshidi and Fallah; ³⁸ however, because the displacement field considered in the present study did not take the effects of the mid-plane displacements into account.

Table 1.

Comparison of dimensionless buckling forces ( $N_{\mathrm{cr}}^{*}$ ) for nonlocal fully simply supported (SSSS) FG nanoplate when subjected to in-plane biaxial compression.

$\sqrt{\mathrm{ea}/a}$	a/b	a/h	$k=0$		$k=0.5$		$k=1$		$k=2$		$k=10$
			Present	Ref. [38]	Present	Ref. [38]	Present	Ref. [38]	Present	Ref. [38]	Present	Ref. [38]
0.0	0.5	5	5.3834	5.3807	3.5428	3.5414	2.7361	2.7350	2.1175	2.1166	1.5378	1.5367
		10	5.9276	5.9255	3.8589	3.8572	2.9702	2.9694	2.3129	2.3118	1.7528	1.7515
		20	6.0802	6.0797	3.9463	3.9454	3.0357	3.0346	2.3674	2.3665	1.8163	1.8151
	1.0	5	8.0225	8.0212	5.3201	5.3191	4.1182	4.1170	3.1736	3.1729	2.2411	2.2400
		10	9.2940	9.2924	6.0645	6.0633	4.6717	4.6709	3.6329	3.6318	2.7274	2.7260
		20	9.6790	9.6772	6.2850	6.2839	4.8353	4.8341	3.7694	3.7686	2.8845	2.8833
0.1	0.5	5	4.7903	4.7898	3.1536	3.1525	2.4374	2.4369	1.8853	1.8841	1.3690	1.3679
		10	5.2756	5.2748	3.4345	3.4336	2.6444	2.6433	2.0588	2.0579	1.5605	1.5592
		20	5.4137	5.4121	3.5133	3.5121	2.7025	2.7013	2.1073	2.1066	1.6169	1.6158
	1.0	5	6.6999	6.6989	4.4431	4.4422	3.4396	3.4383	2.6503	2.6499	1.8715	1.8708
		10	7.7620	7.7605	5.0649	5.0637	3.9017	3.9009	3.0343	3.0331	2.2773	2.2766
		20	8.0828	8.0819	5.2491	5.2480	4.0385	4.0372	3.1482	3.1474	2.4092	2.4080
0.2	0.5	5	3.6039	3.6028	2.3722	2.3713	1.8822	1.8813	1.4183	1.4172	1.0294	1.0289
		10	3.9688	3.9676	2.5836	2.5827	1.9893	1.9882	1.5490	1.5479	1.1739	1.1728
		20	4.0718	4.0709	2.6428	2.6417	2.0327	2.0319	1.5855	1.5846	1.2163	1.2154
	1.0	5	4.4835	4.4822	2.9731	2.9723	2.3017	2.3006	1.7741	1.7730	1.2526	1.2517
		10	5.1938	5.1925	3.3890	3.3881	2.6113	2.6101	2.0304	2.0294	1.5244	1.5233
		20	5.4089	5.4076	3.5129	3.5114	2.7021	2.7012	2.1070	2.1059	1.6121	1.6112

The second comparison is presented for the nonlocal strain gradient square nanoplates when subjected to in-plane biaxial compression and with consideration of the surface effects. Several plate theories and three types of boundary conditions are considered in this comparison. The dimensionless buckling forces ( ${\hat{N}}_{\mathrm{cr}}=\frac{100N_{\mathrm{cr}}}{\mathrm{Eh}}$ ) obtained in this work and those determined by Lu et al. ³⁹ are compared and the results are shown in Tables 2 to 4. It can be seen that the results presented here are the same as those of Lu et al. ³⁹ from the literature. Note that the results for the homogeneous nanoplates can be determined by letting the power-law index k tend to infinity ( $k\to \infty$ ). In this comparison calculation, the value of the power-law index k was set at 10,000.

Table 2.

Comparison of dimensionless buckling forces ( ${\hat{N}}_{\mathrm{cr}}$ ) for nonlocal strain gradient square nanoplate when subjected to in-plane biaxial compression (SSSS, $a/h=10$ ).

$\mathrm{ea}/h$	Plate theory	SSSS
		$l/h=0$		$l/h=1$		$l/h=2$
		Present	Ref. [39]	Present	Ref. [39]	Present	Ref. [39]
0.0	KPT	1.8470	1.8470	2.2116	2.2116	3.3053	3.3053
	MPT	1.8006	1.8006	2.1560	2.1560	3.2223	3.2223
	TSDPT	1.7949	1.7949	2.1492	2.1492	3.2121	3.2121
	SSDPT	1.7946	1.7946	2.1489	2.1489	3.2116	3.2116
	HSDPT	1.7950	1.7950	2.1493	2.1493	3.2122	3.2122
	ESDPT	1.7945	1.7945	2.1487	2.1487	3.2113	3.2113
1.0	KPT	1.5425	1.5425	1.8470	1.8470	2.7604	2.7604
	MPT	1.5038	1.5038	1.8006	1.8006	2.6911	2.6911
	TSDPT	1.4990	1.4990	1.7949	1.7949	2.6826	2.6826
	SSDPT HSDPT ESDPT	1.4988	1.4988	1.7946	1.7946	2.6821	2.6821
		1.4990	1.4990	1.7950	1.7950	2.6827	2.6827
		1.4986	1.4986	1.7945	1.7945	2.6819	2.6819
2.0	KPT	1.0321	1.0321	1.2358	1.2358	1.8470	1.8470
	MPT	1.0062	1.0062	1.2048	1.2048	1.8006	1.8006
	TSDPT	1.0030	1.0030	1.2010	1.2010	1.7949	1.7949
	SSDPT	1.0028	1.0028	1.2008	1.2008	1.7946	1.7946
	HSDPT	1.0030	1.0030	1.2010	1.2010	1.7950	1.7950
	ESDPT	1.0027	1.0027	1.2007	1.2007	1.7945	1.7945

Table 3.

Comparison of dimensionless buckling forces ( ${\hat{N}}_{\mathrm{cr}}$ ) for nonlocal strain gradient square nanoplate when subjected to in-plane biaxial compression (CCSS, $a/h=10$ ).

$\mathrm{ea}/h$	Plate theory	CCSS
		$l/h=0$		$l/h=1$		$l/h=2$
		Present	Ref. [39]	Present	Ref. [39]	Present	Ref. [39]
0.0	KPT	3.0176	3.0176	3.9905	3.9905	6.9093	6.9093
	MPT	2.8413	2.8313	3.7442	3.7442	6.4828	6.4828
	TSDPT	2.8097	2.8097	3.7156	3.7156	6.4334	6.4334
	SSDPT	2.8086	2.8086	3.7141	3.7141	6.4308	6.4308
	HSDPT	2.8098	2.8098	3.7158	3.7158	6.4337	6.4337
	ESDPT	2.8081	2.8081	3.7134	3.7134	6.4296	6.4296
1.0	KPT	2.2819	2.2819	3.0176	3.0176	5.2247	5.2247
	MPT	2.1410	2.1410	2.8313	2.8313	4.9022	4.9022
	TSDPT	2.1247	2.1247	2.8097	2.8097	4.8648	4.8648
	SSDPT	2.1238	2.1238	2.8086	2.8086	4.8629	4.8629
	HSDPT	2.1248	2.1248	2.8098	2.8098	4.8651	4.8651
	ESDPT	2.1234	2.1234	2.8081	2.8081	4.8620	4.8620
2.0	KPT	1.3179	1.3179	1.7428	1.7428	3.0176	3.0176
	MPT	1.2365	1.2365	1.6352	1.6352	2.8313	2.8313
	TSDPT	1.2271	1.2271	1.6228	1.6228	2.8097	2.8097
	SSDPT	1.2266	1.2266	1.6221	1.6221	2.8086	2.8086
	HSDPT	1.2278	1.2272	1.6228	1.6228	2.8098	2.8098
	ESDPT	1.2264	1.2264	1.6218	1.6218	2.8081	2.8081

Table 4.

Comparison of dimensionless buckling forces ( ${\hat{N}}_{\mathrm{cr}}$ ) for nonlocal strain gradient square nanoplate when subjected to in-plane biaxial compression (fully clamped or CCCC, $a/h=10$ ).

$\mathrm{ea}/h$	Plate theory	CCCC
		$l/h=0$		$l/h=1$		$l/h=2$
		Present	Ref. [39]	Present	Ref. [39]	Present	Ref. [39]
0.0	KPT	3.9876	3.9876	5.7719	5.7719	11.1248	11.1248
	MPT	3.6725	3.6725	5.3158	5.3158	10.2458	10.2458
	TSDPT	3.6366	3.6366	5.2638	5.2638	10.1455	10.1455
	SSDPT	3.6347	3.6347	5.2611	5.2611	10.1403	10.1403
	HSDPT	3.6368	3.6368	5.2641	5.2641	10.1461	10.1461
	ESDPT	3.6339	3.6339	5.2599	5.2599	10.1378	10.1378
1.0	KPT	2.7549	2.7549	3.9876	3.9876	7.6858	7.6858
	MPT	2.5372	2.5372	3.6725	3.6725	7.0785	7.0785
	TSDPT	2.5124	2.5124	3.6366	3.6366	7.0092	7.0092
	SSDPT	2.5111	2.5111	3.6347	3.6347	7.0056	7.0056
	HSDPT	2.5126	2.5126	3.6368	3.6368	7.0096	7.0096
	ESDPT	2.5105	2.5105	3.6339	3.6339	7.0039	7.0039
2.0	KPT	1.4294	1.4294	2.0689	2.0689	3.9876	3.9876
	MPT	1.3164	1.3164	1.9054	1.9054	3.6725	3.6725
	TSDPT	1.3035	1.3035	1.8868	1.8868	3.6366	3.6366
	SSDPT	1.3029	1.3029	1.8858	1.8858	3.6347	3.6347
	HSDPT	1.3036	1.3036	1.8869	1.8869	3.6368	3.6368
	ESDPT	1.3025	1.3025	1.8854	1.8854	3.6339	3.6339

KPT: Kirchhoff plate theory; MPT: Mindlin plate theory; TSDPT: Third-order shear deformation plate theory; SSDPT: Sinusoidal shear deformable plate theory; HSDPT: Hyperbolic shear deformable plate theory; ESDPT: Exponential shear deformable plate theory.

4.2. Parametric study

In this subsection, the effects of some important parameters on the buckling behavior of the FG microplates will be examined. For this purpose, the selected FG microplate thickness is $h=10\mu m$ .

The mechanical behavior of micro/nanostructures is strongly influenced by size-dependent effects. Numerous experimental and theoretical studies have been reported that have focused on the development of techniques to quantify the nonlocal parameters and the material length scale parameters across the different materials, along with various nonclassical elasticity theories^9,11,45–49. Lam et al. ¹¹ determined the material length scale parameter experimentally for a homogeneous epoxy microbeam in MCST to be 17.6 μm. However, based on SGT, Kahrobaiyan et al. ⁴⁵ inferred the material length scale parameter for an epoxy microbeam to be 11.01 μm. Using experimental results, Lei et al. ⁴⁶ reported values of the material length scale parameter for a nickel cantilever microbeam of 0.843 μm and 1.553 μm corresponding to MCST and SGT, respectively. The material length scale parameter for a polymer SU-8 microbeam in SGT was reported by Liebold and Müller, ⁴⁷ based on an atomic force microscopy experiment, to be 2.5 μm. Through matching of the dispersion curves for plane waves in NET and the Born-Karman model of lattice dynamics in a phonon analysis, the value of e₀ computed by Eringen ⁹ was $e_0\approx 0.39$ . The nonlocal parameter value for carbon nanotubes was estimated to be in the $0<\mu <2\mathrm{nm}$ range based on theoretical and experimental results reported by Zhang et al. ⁴⁸ and Wang. ⁴⁹ From the data available in the literature, it can be concluded that the values of the nonlocal and material length scale parameters are dependent on both the materials used to form the micro/nanostructures and the nonclassical elasticity theory that is used. To date, no experimental or theoretical data for the nonlocal and material length scale parameters of FG microplates have been reported in the literature. Therefore, in this study, to examine the effects of both the nonlocal parameter (ea) and the material length scale parameter (l) on the buckling response of FG microplates, the values of these parameters are normalized with respect to the thickness (h) of the FG microplate in numerical simulations.

The effects of the nonlocal parameter and the material length scale parameter on the buckling behavior of FG microplates with consideration of the surface layers are illustrated by the results presented in Table 5, and in Figures 2 and 3. The dimensionless buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) values for FG square microplates with respect to different boundary conditions and several plate theories are tabulated in Table 5. Figures 2 and 3 show the variations in the dimensionless buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) of the FG microplate versus the dimensionless nonlocal parameter ( $\mathrm{ea}/h$ ) and the dimensionless material length scale parameter ( $l/h$ ), respectively. The results shown in Figures 2 and 3 were examined using TSDPT. As expected, the results showed that when the nonlocal parameter is equal to the material length scale parameter ( $\mathrm{ea}=l$ ), the dimensionless critical buckling forces based on NSGT are then equal to those obtained based on classical elasticity theory ( $\mathrm{ea}=l=0$ ). The FG microplates exhibit the hardening-stiffness effect when the material length scale parameter is larger than the nonlocal parameter ( $l>\mathrm{ea}$ ); this means that the dimensionless critical buckling forces of FG microplates based on NSGT are always greater than those of classical FG microplates. In contrast, when the nonlocal parameter is greater than the material length scale parameter ( $\mathrm{ea}>l$ ), the dimensionless critical buckling forces of FG microplates based on NSGT are always smaller than those of classical FG microplates; this means that the FG microplates exert the softening-stiffness effect. Additionally, Figures 2 and 3 show that the nonlocal parameter leads to a reduction in the critical buckling forces, while the material length scale parameter leads to an increase in the critical forces for the FG microplates. This observation can be explained as follows: the nonlocal parameter has the effect of reducing the stiffness of the FG microplate, whereas the material length scale parameter increases the stiffness of the FG microplate. To enable better observation of the softening- and hardening-stiffness effects in the FG microplates, Figure 4 shows the variations in the dimensionless critical buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) versus the scale ratio ( $\delta =l/\mathrm{ea}$ ) for various boundary conditions. The results presented in Figure 4 were determined by using TSDPT. It is obvious that the dimensionless critical buckling forces of the classical FG microplates are not dependent on the scale ratio ( $\delta$ ). Furthermore, Figure 4 shows that the dimensionless critical buckling forces of FG microplates based on NSGT are always smaller than those of classical FG microplates if the scale ratio is smaller than unity; in addition, the dimensionless critical buckling forces of FG microplates based on NSGT are always greater than those of classical FG microplates if the scale ratio is higher than unity. When the scale ratio δ is equal to unity (i.e., $\delta =1$ ), the results for NSGT match those obtained from classical theory. It can also be seen from Table 5 and Figures 2 to 4 that the dimensionless critical buckling forces of the FG microplates have the highest values corresponding to the CCCC boundary condition and the lowest values corresponding to the SSSS boundary condition. In addition, the dimensionless critical buckling forces for KPT are always higher than those observed for the other plate theories.

Table 5.

Dimensionless critical buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) values of FG square microplates with consideration of the surface layers ( $a=20h,k=5$ ).

$\mathrm{ea}/h$	Plate theory	SSSS			CCCC			CCSS
		$l/h=0.0$	$l/h=0.5$	$l/h=1.0$	$l/h=0.0$	$l/h=0.5$	$l/h=1.0$	$l/h=0.0$	$l/h=0.5$	$l/h=1.0$
0.0	KPT	1.024	1.026	1.032	2.748	2.760	2.797	1.967	1.973	1.992
	MPT	1.005	1.007	1.013	2.619	2.631	2.666	1.891	1.897	1.915
	TSDPT	1.002	1.004	1.010	2.595	2.606	2.641	1.876	1.882	1.901
	SSDPT	1.002	1.004	1.009	2.593	2.605	2.640	1.876	1.882	1.900
	HSDPT	1.002	1.004	1.010	2.595	2.606	2.641	1.876	1.883	1.901
	ESDPT	1.001	1.003	1.009	2.593	2.604	2.639	1.875	1.881	1.899
0.5	KPT	1.022	1.024	1.030	2.736	2.748	2.785	1.960	1.967	1.986
	MPT	1.003	1.005	1.011	2.607	2.619	2.654	1.885	1.891	1.909
	TSDPT	1.000	1.002	1.008	2.583	2.595	2.629	1.870	1.876	1.894
	SSDPT	1.000	1.002	1.007	2.582	2.593	2.628	1.870	1.876	1.894
	HSDPT	1.000	1.002	1.008	2.583	2.595	2.629	1.870	1.876	1.895
	ESDPT	0.999	1.001	1.007	2.581	2.593	2.627	1.869	1.875	1.893
1.0	KPT	1.016	1.018	1.024	2.700	2.712	2.748	1.942	1.948	1.967
	MPT	0.997	0.999	1.005	2.573	2.584	2.619	1.867	1.873	1.891
	TSDPT	0.994	0.996	1.002	2.549	2.560	2.595	1.853	1.858	1.876
	SSDPT	0.994	0.996	1.002	2.548	2.559	2.593	1.852	1.858	1.876
	HSDPT	0.994	0.996	1.002	2.549	2.560	2.595	1.853	1.859	1.876
	ESDPT	0.994	0.996	1.001	2.547	2.558	2.593	1.851	1.857	1.875

Figure 2.

Variations in the dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ versus the dimensionless nonlocal parameter ea/h ( $a/h=20,l/h=0.2,b=1.5a$ ).

Figure 3.

Variations in the dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ versus the dimensionless material length scale parameter l/h ( $a/h=10,\mathrm{ea}/h=0.3,b=a$ ).

Figure 4.

Variation in the dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ versus the scale ratio δ ( $a/h=20,k=1,b=a$ ).

The influence of the power-law index k on the buckling behavior of the FG microplate is illustrated in Figure 5 and Table 6. Figure 5 shows the variation in the dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ with respect to the power-law index k based on TSDPT and with consideration of the surface layers. The dimensionless critical buckling forces of the FG microplates with consideration of the surface layers and with respect to several differential plate theories are listed in Table 6. Clearly, the critical buckling force decreases as the value of the power-law index increases, especially for smaller values of the power-law index. This reflects the fact that increasing the value of the power-law index k leads to a reduction in the stiffness of the FG microplates, which in turn reduces the critical buckling force. Furthermore, the results in Table 6 show that the dimensionless critical buckling force of the FG microplates decreases as the width to length ratio ( $b/a$ ) increases.

Figure 5.

Variation in the dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ versus the power-law index k ( $a=30h,\mathrm{ea}=0.1h,l=0.2h$ ).

Table 6.

Dimensionless critical buckling forces ( ${\overline{N}}_{\mathrm{cr}}$ ) of FG microplates with consideration of the surface layers ( $a=30h,\mathrm{ea}=0.1h,l=0.2h$ ).

b/a	Plate theory	k = 0			k = 1			k = 10
		SSSS	CCCC	CCSS	SSSS	CCCC	CCSS	SSSS	CCCC	CCSS
0.5	KPT	2.730	8.761	3.239	1.617	1.617	1.918	0.939	3.013	1.114
	MPT	2.688	8.324	3.178	1.592	4.929	1.882	0.920	2.817	1.086
	TSDPT	2.688	8.325	3.178	1.592	4.929	1.882	0.917	2.783	1.081
	SSDPT	2.688	8.325	3.178	1.592	4.929	1.882	0.916	2.782	1.081
	HSDPT	2.688	8.325	3.178	1.592	4.929	1.882	0.917	2.783	1.081
	ESDPT	2.689	8.327	3.178	1.592	4.931	1.882	0.917	2.783	1.081
1	KPT	1.091	2.931	2.097	0.646	1.736	1.242	0.375	1.008	0.721
	MPT	1.084	2.883	2.069	0.642	1.707	1.225	0.372	0.986	0.708
	TSDPT	1.084	2.883	2.069	0.642	1.707	1.225	0.372	0.982	0.706
	SSDPT	1.084	2.883	2.069	0.642	1.707	1.225	0.372	0.982	0.706
	HSDPT	1.084	2.883	2.069	0.642	1.707	1.225	0.372	0.982	0.706
	ESDPT	1.084	2.883	2.069	0.642	1.707	1.225	0.372	0.982	0.706
2	KPT	0.682	2.184	2.124	0.404	1.293	1.257	0.235	0.751	0.731
	MPT	0.679	2.156	2.096	0.402	1.277	1.241	0.233	0.738	0.718
	TSDPT	0.679	2.156	2.096	0.402	1.277	1.241	0.233	0.736	0.716
	SSDPT	0.679	2.156	2.096	0.402	1.277	1.241	0.233	0.736	0.716
	HSDPT	0.679	2.156	2.096	0.402	1.277	1.241	0.233	0.736	0.716
	ESDPT	0.679	2.156	2.096	0.402	1.277	1.241	0.233	0.736	0.716

The effect of the length-to-thickness ratio ( $a/h$ ) on the buckling behavior of the FG microplates is illustrated by the results presented in Figure 6 and Table 7. The results obtained show that the dimensionless critical buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) of the FG microplates decrease when the value of $a/h$ increases. Furthermore, the effect of the thickness on the buckling behavior of the FG microplates was also investigated. Figure 7 represents the variation in the dimensionless critical buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) versus the thickness (h) of the FG microplates for TSDPT. The figure shows that the thickness of the FG microplates has an interesting effect on the buckling behavior of these microplates for two cases of the boundary conditions in particular, comprising the SSSS and CCSS conditions; in contrast, for the CCCC boundary condition, the dimensionless critical buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) increases with increasing FG microplate thickness. However, the critical buckling forces do not change greatly when the magnitude of the FG microplate thickness is changed. The largest relative differences observed in the dimensionless critical buckling forces for the SSSS and CCSS boundary conditions are approximately 0.29% and 0.31%, respectively.

Figure 6.

Variation in the dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ versus the length-to-thickness ratio $a/h$ ( $a=b,\mathrm{ea}=0.1h,l=0.2h,k=1$ ).

Table 7.

Dimensionless critical buckling force ( ${\overline{N}}_{\mathrm{cr}}$ ) of FG square microplates with consideration of the surface layers ( $\mathrm{ea}=0.2h,l=0.1h$ ).

a/h	Boundary condition	$k=0$		$k=10$		$k=20$
		With surface layers	Without surface layers	With surface layers	Without surface layers	With surface layers	Without surface layers
5	SSSS	31.3061	31.3064	9.7330	9.7331	8.2244	8.2246
	CCCC	62.5059	62.5064	17.6793	17.6796	15.3691	15.3694
	CCSS	49.0208	49.0212	14.2508	14.2511	12.2863	12.2866
10	SSSS	9.2347	9.2347	3.0809	3.0808	2.5460	2.5459
	CCCC	22.5819	22.5819	7.2179	7.2179	6.0468	6.0468
	CCSS	16.6334	16.6334	5.3829	5.3830	4.4917	4.4917
20	SSSS	2.4157	2.4155	0.8243	0.8242	0.6762	0.6761
	CCCC	6.3244	6.3243	2.1302	2.1301	1.7549	1.7548
	CCSS	4.5612	4.5611	1.5422	1.5421	1.2689	1.2688
50	SSSS	0.3917	0.3915	0.1347	0.1345	0.1103	0.1101
	CCCC	1.0469	1.0468	0.3589	0.3588	0.2940	0.2939
	CCSS	0.7502	0.7501	0.2574	0.2573	0.2108	0.2107

Figure 7.

Variation in the dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ versus the thickness h ( $a=b,\mathrm{ea}=0.1h,l=0.2h,k=10$ ).

Finally, the effect of the surface layers on the buckling behavior of the FG microplates was also investigated and the results are shown in Table 7. The dimensionless critical buckling force ${\overline{N}}_{\mathrm{cr}}$ was obtained by applying TSDPT for this investigation. The results in Table 7 demonstrate that the effect of the surface layers on the buckling behavior of the FG microplates is dependent on the value of the length-to-thickness ratio (a/h). Specifically, when the ratio is small (e.g., $a/h=5$ ), the dimensionless critical buckling forces of the FG microplates with surface layers are only a little smaller than those of the FG microplates without the surface layers. However, for the higher values of the length-to-thickness ratio (e.g., $a/h=20,50$ ), the dimensionless critical buckling forces of the FG microplates with surface layers are slightly larger than those of the corresponding FG microplates without the surface layers. This result can be explained by the fact that the influence of the surface layers on the buckling behavior of the FG microplates is dependent on the choice of surface materials.

Concluding remarks

A size-dependent buckling analysis of FG microplates is performed in this work by using NSGT with incorporation of surface effects and higher-order shear deformation plate theory. The governing equations for the FG microplates are derived by using the principle of minimum total potential energy. The analytical solutions are used to perform buckling analyses of the FG microplates with three types of boundary conditions, including SSSS, CCCC, and CCSS. Comparisons of the results obtained with the existing results in the literature are performed for different elasticity theories. The effects of some important parameters on the buckling responses of the FG microplates are also considered. From the results obtained in this work, the following conclusions can be drawn:

• The buckling forces obtained for the higher-order shear deformation plate theories are always lower than those obtained for both KPT and MPT because they consider the shear deformation effect.

• The critical buckling forces of the FG microplates can be increased by increasing the material length scale parameter (l) or by reducing the nonlocal parameter (ea). The FG microplates demonstrate the hardening-stiffness effect when the material length scale parameter is higher than the nonlocal parameter but demonstrate the softening-stiffness effect when the opposite is true. The results for the classical FG microplates can be obtained from the nonlocal strain gradient FG microplates by setting $\mathrm{ea}=0$ and $l=0.$

• The power-law index has an important influence on the buckling behavior of the FG microplates. The critical buckling forces decrease as the power-law index (k) increases, especially at smaller values of the power-law index. The critical buckling forces for the full metal and full ceramic microplates can be acquired by setting the value of k to approach infinity and to be equal to zero, respectively.

• The critical buckling forces decrease when the length-to-thickness ratio (a/h) and the width-to-length ratio (b/a) increase. Additionally, the thickness (h) does have an effect on the buckling response of the FG microplates but it is minimal.

• The effect of the surface layers on the buckling behavior of the FG microplates is dependent on the value of the length-to-thickness ratio (a/h).