Effect of Thickness Stretching on Bending and Free Vibration Behaviors of Functionally Graded Graphene Reinforced Composite Plates

Abstract

The focus of this paper is the effect of thickness stretching on the static and dynamic behaviors of functionally graded graphene reinforced composite (FG-GRC) plates. The bending and free vibration behaviors of FG-GRC plates under simply supported conditions are studied based on two plate theories, with or without taking into account the thickness stretching effect, respectively, and the effect of thickness stretching on FG-GRC plates is analyzed by comparing the calculated results of the two types of plate theories. The properties of composite materials are estimated by the modified Halpin-Tsai model and rule of mixture, Hamilton’s principle is used to construct its governing equation, and the Navier solution method is used to find the closed solution. The numerical results show that the effect of thickness stretching depends mainly on the transverse anisotropy of the FG-GRC plates, and the FG-GRC plates are most significantly affected by the thickness stretching when the graphene nanoplatelets (GPLs) are asymmetrically distributed, and the effect of thickness stretching tends to increase as the total number of layers and the weight fraction of GPLs increase.

1. Introduction

The concept of functionally graded materials has been proposed by Japanese materials scientists since 1984 [1], among which traditional metal-ceramic functionally graded materials have been widely used in aerospace engineering. Owing to their excellent thermal insulation and high strength, functionally graded materials can be used to make high-performance parts in spacecraft [2,3]. This paper deals with a new class of functionally gradient materials, which enhance the mechanical properties by adding reinforcement to the matrix materials.

Graphene nanoplatelets (GPLs) have attracted a significant amount of attention from researchers since their discovery by Novoselov et al. in 2004 [4]. Owing to the high Young’s modulus and high strength of GPLs, we can obtain composites with excellent mechanical properties by incorporating GPLs into the matrix materials. Since GPLs have a larger specific surface area than carbon nanotubes (about three times greater than carbon nanotubes), the mechanical properties of GPLs-reinforced composites are considerably improved over the previously studied carbon nanotubes reinforced composites [5]. The excellent performance of GPLs-reinforced composites gives them a very broad range of application scenarios in aerospace, civil engineering, biomedical, etc. Therefore, it is important to study the mechanical behaviors of functionally graded graphene reinforced composites (FG-GRC).

Some researchers have chosen classical plate theory [6] (Kirchhoff theory) for the analysis of FG-GRC. Yang et al. [7] investigated the critical buckling load of spinning FG-GRC cylindrical shells by using classical plate theory. Gao et al. [8] studied the FG-GRC porous plates using classical plate theory and solved the nonlinear free vibration of the plates using the differential quadrature method. Li et al. [9] investigated the nonlinear vibration of FG-GRC sandwich porous plates based on classical plate theory.

Classical plate theory neglects the effects of shear and normal deformation and is only applicable to thin FG-GRC sheet/shell, so more researchers have adopted first-order shear deformation plate theory (FSDPT) [10] and third-order shear deformation plate theory (TSDPT) [11]. Song et al. [12] analyzed the free and forced vibrations of FG-GRC plates employing FSDPT. Garcia-Macias et al. [13] calculated the bending and free vibration of FG-GRC plates based on FSDPT using the finite element method, and also calculated the mechanical behaviors of carbon nanotube reinforced plates for comparison purposes. Mohammod et al. [14] performed a free vibration analysis of FG-GRC curved nanobeams on Pasternak foundation based on the displacement field described by FSDPT. Guo et al. calculated nonlinear bending [15] and free vibration [16] of an FG-GRC quadrilateral plate based on FSDPT, using the element-free IMLS-Ritz method. Reddy et al. [17] used FSDPT and finite element method to analyze the free vibration characteristics of FG-GRC plates under various boundary conditions.

Due to the assumption of a constant transverse shear strain in the transverse section, the FSDPT must be shear-corrected, and this correction factor is not only related to the materials and geometric parameters, but also to the loads as well as to the boundary conditions. Since the shear correction factor is difficult to determine, TSDPT becomes the favored choice of many researchers. The advantage of TSDPT is that it does not need to consider the shear correction factor, and the calculation result of TSDPT is also closer to the 3D solution than FSDPT. Shen et al. [18] established geometric nonlinear equations of motion with von Kármán strain form based on higher-order shear deformation plate theory to study the significant amplitude vibration in FG-GRC plates under an elastic foundation and thermal environment. Kiani [19] used Reddy-based TSDPT to obtain the total strain energy and performed thermal post-buckling analysis on FG-GRC plates. Gholami et al. [20] investigated the geometrically nonlinear harmonically excited vibration of FG-GRC rectangular plates under different edge conditions using TSDPT. Li et al. [21] studied the bending, buckling and free vibration of FG-GRC plates using isogeometric analysis based on FSDPT and TSDPT.

In addition to FSDPT and TSDPT, there are some articles that use other plate theories. Thai et al. [22] investigated the bending, buckling and free vibration behaviors of FG-GRC plates based on refined shear deformation plate theory. Gholami and Ansari [23] performed a large-deflection geometric nonlinear analysis of FG-GRC plates based on sinusoidal shear deformation plate theory. Based on sinusoidal shear deformation plate theory, Arefi et al. [24] reported the effect of small-scale effects on the natural frequencies of FG-GRC nanoplates on elastic foundations.

All the plate theories mentioned in the above review do not take into account the effect of thickness stretching, which can lead to large errors when we calculate moderately thick and thick plates. To compensate for the disadvantages of those theories, many researchers have proposed plate theories that take into account the effect of thickness stretching. For example, Zenkour [25] proposed a quasi-3D higher-order shear deformation plate theory based on sinusoidal functions of in-plane and transverse displacements, while the quasi-3D theory [26] is based on the parabolic function of transverse displacement and the cubic function of in-plane displacement. Both of the above two plate theories contain six unknown quantities in the displacement field. There are also three high-order shear deformation plate theories with nine unknowns proposed by Neves [27,28,29].

Senthilnathan et al. [30] simplified the TSDPT by dividing the transverse displacement into bending and shear components so that the number of unknowns was reduced by one compared to the TSDPT. They applied this assumption to the analysis of the isotropic plates and proposed a refined plate theory (RPT). This hypothesis of dividing the displacement into a bending part and a shear part was first proposed by Huffington [31]. Shimp [32] published a detailed analysis of homogeneous plates based on RPT in 2002 and proposed two extended theories based on RPT. Mechab [33] first extended RPT to the analysis of functionally graded materials. Since then, a large number of documents [34,35,36,37,38,39,40] have analyzed and studied the mechanical behaviors of functionally graded materials based on RPT, which showed the good applicability of RPT in the mechanical analysis of functionally graded materials.

The calculation result of RPT is very close to TSDPT [30], and because it reduces one degree of freedom, the calculation process is greatly simplified. However, RPT does not consider the effect of thickness stretching, so its shortcomings are the same as TSDPT. In order to make up for the shortcomings of RPT, Thai [41] proposes an improved refined plate theory (IRPT) that takes into account the effect of thickness stretching. IRPT corrects the disadvantage of RPT in the calculation of moderately thick and thick plates by adding a thickness stretching effect term in the displacement field. The calculation results [41] in functionally graded materials show that the IRPT that takes into account the effect of thickness stretching is more accurate than the RPT that does not consider the effect of thickness stretching. Compared with the previously mentioned plate theories that consider the effect of thickness stretching, IRPT contains the least number of unknowns [41].

The importance of the effect of thickness stretching for the analysis of the mechanical behaviors of metal-ceramic functionally graded plates was emphasized by Carrera et al. [42]. Usually, the effect of thickness stretching can be disregarded for thin isotropic plates, but due to the specificity of the structure of the functionally graded materials, ignoring the effect of thickness stretching in the analysis of mechanical behaviors may lead to inaccurate calculation results [42]. For moderately thick and thick plates, neglecting the effect of thickness stretching produces non-negligible errors, which is evident in both multilayered composite plates and metal-ceramic functionally graded plates.

To the authors’ knowledge, few researchers have further studied the effect of thickness stretching on the mechanical behaviors of FG-GRC plates in detail. In the present paper, attention is focused on the effect of thickness stretching on the static and dynamic behaviors of FG-GRC plates. The novelty of this paper is reflected in the fact that the importance of the thickness stretching effect is investigated by comparing the numerical results based on two plate theories, with or without taking into account the thickness stretching effect, respectively. In this paper, RPT and IRPT are used to study the bending and free vibration behaviors of FG-GRC plates. The properties of composite materials are estimated by the modified Halpin-Tsai model and rule of mixture. Hamilton’s principle is used to construct its governing equation, and the Navier solution method is used to find the closed solution. The effects of thickness stretching, distribution modes of GPLs, and some parameters, such as the total number of layers $N_L$, and weight fraction $f_{\mathrm{G}}$ of GPLs on frequencies and deflections of the FG-GRC rectangular simply supported plates are discussed in detail.

2. Effective Material Properties

Consider an FG-GRC plate with length a, width b, and thickness h (Figure 1) using the cartesian coordinate system x, y, z, with plane z = 0, coinciding with the midplane of the plate. $N_L$ is the total number of layers of FG-GRC plates, and it is worth noting that this paper is only for FG-GRC plates with an even number of layers. GPLs are used as fillers and four distribution patterns proposed by Song et al. [12,43] are used in this paper. In Figure 2, taking the plate with $N_L=6$ as an example, a schematic diagram of the four distribution modes of GPLs is given [12,43].

In the UD mode, the GPLs content is uniformly distributed. In the FG-O mode, the GPLs content decrease from the middle surface to the surface of the plate. In the FG-X mode, the distribution is opposite to in the FG-O mode, and the GPLs content decreases from the surface to the middle surface of the plate. In the FG-A mode, the GPLs content increases from the upper surface to the lower surface.

The effective Young’s modulus of the kth layer of the FG-GRC plates, $E_c^{\left(k\right)}$, can be calculated based on the Halpin-Tsai micromechanics model [44] and Voigt-Reuss model [45], as follows:

$$E_c^{(k)}=\frac{3}{8}\times \frac{1+{\xi }_{\mathrm{L}}{\eta }_{\mathrm{L}}{V}_{\mathrm{G}}^{(k)}}{1-{\eta }_{\mathrm{L}}{V}_{\mathrm{G}}^{(k)}}\times E_{\mathrm{M}}+\frac{5}{8}\times \frac{1+{\xi }_{\mathrm{W}}{\eta }_{\mathrm{W}}{V}_{\mathrm{G}}^{(k)}}{1-{\eta }_{\mathrm{W}}{V}_{\mathrm{G}}^{(k)}}\times E_{\mathrm{M}}$$

(1)

where

$${\eta }_{\mathrm{L}}=\frac{\left(E_{\mathrm{G}}/E_{\mathrm{M}}\right)-1}{\left(E_{\mathrm{G}}/E_{\mathrm{M}}\right)+{\xi }_{\mathrm{L}}},{\eta }_{\mathrm{W}}=\frac{\left(E_{\mathrm{G}}/E_{\mathrm{M}}\right)-1}{\left(E_{\mathrm{G}}/E_{\mathrm{M}}\right)+{\xi }_{\mathrm{W}}}$$

(2)

k = 1, 2, …, $N_L$, and $E_{\mathrm{G}}$, $E_{\mathrm{M}}$ represent Young’s modulus of GPLs and polymer matrix, respectively. ${\xi }_{\mathrm{L}}$ and ${\xi }_{\mathrm{W}}$ are the parameters associated with the geometry and size of the GPLs, defined as

$${\xi }_{\mathrm{W}}=2\left(\frac{w_{\mathrm{G}}}{h_{\mathrm{G}}}\right),{\xi }_{\mathrm{L}}=2\left(\frac{l_{\mathrm{G}}}{h_{\mathrm{G}}}\right)$$

(3)

with $l_{\mathrm{G}}$, $w_{\mathrm{G}}$, and $h_{\mathrm{G}}$ being the average length, width, and thickness of the GPLs. In addition, the volume fraction $V_{\mathrm{G}}^{\left(k\right)}$ of GPLs in the kth layer is defined as follows:

$$V_{\mathrm{G}}^{(k)}=\frac{f_{\mathrm{G}}^{(k)}}{f_{\mathrm{G}}^{(k)}+\left({\rho }_{\mathrm{G}}/{\rho }_{\mathrm{M}}\right)\left(1-f_{\mathrm{G}}^{(k)}\right)}$$

(4)

where $f_{\mathrm{G}}^{\left(k\right)}$ is GPLs weight fraction of the kth layer, ${\rho }_{\mathrm{M}}$ and ${\rho }_{\mathrm{G}}$ denote the mass densities of the polymer matrix and GPLs, respectively.

According to the rule of mixture, the effective density and Poisson’s ratio of the kth layer of the FG-GRC plates are defined as follows:

$${\rho }_c^{(k)}={\rho }_{\mathrm{G}}{V}_{\mathrm{G}}^{(k)}+{\rho }_{\mathrm{M}}\left(1-V_{\mathrm{G}}^{(k)}\right),{\nu }_c^{(k)}={\nu }_{\mathrm{G}}{V}_{\mathrm{G}}^{(k)}+{\nu }_{\mathrm{M}}\left(1-V_{\mathrm{G}}^{(k)}\right)$$

(5)

${\nu }_{\mathrm{G}}$ and ${\nu }_{\mathrm{M}}$ are Poisson’s ratios of GPLs and polymer matrix, respectively. Assuming $N_L$ as an even number, the GPLs weight fraction of the kth layer for the four distribution modes is as follows [43]:

$$f_{\mathrm{G}}^{(k)}=\{\begin{array}{ll}{f}_{\mathrm{G}}\hfill & \mathrm{UD}\hfill \\ 4f_{\mathrm{G}}\left(\frac{N_L+1}{2}-\left|k-\frac{N_L+1}{2}\right|\right)/\left(2+N_L\right)\hfill & \mathrm{FG}-\mathrm{O}\hfill \\ 4f_{\mathrm{G}}\left(\frac{1}{2}+\left|k-\frac{N_L+1}{2}\right|\right)/\left(2+N_L\right)\hfill & \mathrm{FG}-\mathrm{X}\hfill \\ 2k{f}_{\mathrm{G}}/\left(N_L+1\right)\hfill & \mathrm{FG}-\mathrm{A}\hfill \end{array}$$

(6)

in which $f_{\mathrm{G}}$ is the total weight fraction of GPLs in the FG-GRC plates.

3. The Plate Theories (RPT and IRPT)

Consider the rectangular plate shown in Figure 1, using the assumptions shown below:

(1)

In-plane displacement and transverse displacement consists of bending and shear components.

(2)

The bending component of the in-plane displacement is assumed to be analogous to the displacement that is given by classical plate theory.

(3)

The shear part of the in-plane displacement causes a parabolic change in the shear strain, thereby generating shear stress through the thickness of the plate so that the shear stress on the top and bottom surfaces disappears.

Based on the above three assumptions, the displacement field of RPT is defined as follows [32,41]:

$$\begin{array}{l}{u}_1(x,y,z,t)=u(x,y,t)-z\frac{\partial w_b}{\partial x}-f\left(z\right)\frac{\partial w_s}{\partial x}\\ u_2(x,y,z,t)=v(x,y,t)-z\frac{\partial w_b}{\partial y}-f\left(z\right)\frac{\partial w_s}{\partial y}\\ u_3\left(x,y,z,t\right)=w_b\left(x,y,t\right)+w_s\left(x,y,t\right)\end{array}$$

(7)

where $f\left(z\right)=-z/4+5z^3/3h^2$; $\left(u_1,u_2,u_3\right)$ are the displacements along the $\left(x,y,z\right)$ axis direction, respectively. $u$ and $v$ represent the displacements of a point on the middle surface of the plate along the x-axis and y-axis; $w_b$ is the bending part of the transverse displacement; and $w_s$ is the shear part of the transverse displacement. h indicates the thickness of the plate.

Considering the effect of thickness stretching based on RPT, we can obtain the displacement field of IRPT as follows [41]:

$$\begin{array}{l}{u}_1(x,y,z,t)=u(x,y,t)-z\frac{\partial w_b}{\partial x}-f\left(z\right)\frac{\partial w_s}{\partial x}\\ u_2(x,y,z,t)=v(x,y,t)-z\frac{\partial w_b}{\partial y}-f\left(z\right)\frac{\partial w_s}{\partial y}\\ u_3\left(x,y,z,t\right)=w_b\left(x,y,t\right)+w_s\left(x,y,t\right)+g\left(z\right)w_z\left(x,y,t\right)\end{array}$$

(8)

where $g\left(z\right)=1-f^{\prime }\left(z\right)$. $w_z$ is the unknown displacement function with respect to the effect of thickness stretching.

The linear relationship between strain and displacement is as follows:

$${\epsilon }_x=\frac{\partial u}{\partial x}-z\frac{{\partial }^2{w}_b}{\partial x^2}-f\left(z\right)\frac{{\partial }^2{w}_s}{\partial x^2}$$

(9)

$${\epsilon }_y=\frac{\partial v}{\partial y}-z\frac{{\partial }^2{w}_b}{\partial y^2}-f\left(z\right)\frac{{\partial }^2{w}_s}{\partial y^2}$$

(10)

$${\epsilon }_z=g^{\prime }\left(z\right)w_z$$

(11)

$${\gamma }_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}-2z\frac{{\partial }^2{w}_b}{\partial x\partial y}-2f\left(z\right)\frac{{\partial }^2{w}_s}{\partial x\partial y}$$

(12)

$${\gamma }_{xz}=g\left(z\right)\left(\frac{\partial w_s}{\partial x}+\frac{\partial w_z}{\partial x}\right)$$

(13)

$${\gamma }_{yz}=g\left(z\right)\left(\frac{\partial w_s}{\partial y}+\frac{\partial w_z}{\partial y}\right)$$

(14)

If the effect of thickness stretching is neglected, Equation (11) becomes ${\epsilon }_z=0$.

The relationship between strain and stress is as follows:

$${\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\tau }_{xy}\\ {\tau }_{xz}\\ {\tau }_{yz}\end{array}\right\}}^{(k)}=\left[\begin{array}{cccccc}{C}_{11}^{\left(k\right)}& C_{12}^{\left(k\right)}& C_{13}^{\left(k\right)}& 0& 0& 0\\ C_{21}^{\left(k\right)}& C_{22}^{\left(k\right)}& C_{23}^{\left(k\right)}& 0& 0& 0\\ C_{31}^{\left(k\right)}& C_{32}^{\left(k\right)}& C_{33}^{\left(k\right)}& 0& 0& 0\\ 0& 0& 0& C_{66}^{\left(k\right)}& 0& 0\\ 0& 0& 0& 0& C_{55}^{\left(k\right)}& 0\\ 0& 0& 0& 0& 0& C_{44}^{\left(k\right)}\end{array}\right]{\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ {\gamma }_{xy}\\ {\gamma }_{xz}\\ {\gamma }_{yz}\end{array}\right\}}^{\left(k\right)}$$

(15)

where $C_{ij}^{\left(k\right)}$ is the three-dimensional elastic constant as follows:

$$C_{11}^{\left(k\right)}=C_{22}^{\left(k\right)}=C_{33}^{\left(k\right)}=\frac{\left(1-{\nu }_c^{\left(k\right)}\right)E_c^{\left(k\right)}}{\left(1-2{\nu }_c^{\left(k\right)}\right)\left(1+{\nu }_c^{\left(k\right)}\right)}$$

(16)

$$C_{12}^{\left(k\right)}=C_{13}^{\left(k\right)}=C_{23}^{\left(k\right)}=\frac{{\nu }_c^{\left(k\right)}{E}_c^{\left(k\right)}}{\left(1-2{\nu }_c^{\left(k\right)}\right)\left(1+{\nu }_c^{\left(k\right)}\right)}$$

(17)

$$C_{44}^{\left(k\right)}=C_{55}^{\left(k\right)}=C_{66}^{\left(k\right)}=\frac{E_c^{\left(k\right)}}{2\left(1+{\nu }_c^{\left(k\right)}\right)}$$

(18)

If the thickness stretching effect is neglected, the stress–strain relationship is shown below.

$${\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\tau }_{xy}\\ {\tau }_{xz}\\ {\tau }_{yz}\end{array}\right\}}^{(k)}=\frac{E_c^{\left(k\right)}}{1-{\left[{\nu }_c^{\left(k\right)}\right]}^2}\left[\begin{array}{ccccc}1& {\nu }_c^{\left(k\right)}& 0& 0& 0\\ {\nu }_c^{\left(k\right)}& 1& 0& 0& 0\\ 0& 0& \frac{1-{\nu }_c^{\left(k\right)}}{2}& 0& 0\\ 0& 0& 0& \frac{1-{\nu }_c^{\left(k\right)}}{2}& 0\\ 0& 0& 0& 0& \frac{1-{\nu }_c^{\left(k\right)}}{2}\end{array}\right]{\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{xy}\\ {\gamma }_{xz}\\ {\gamma }_{yz}\end{array}\right\}}^{\left(k\right)}$$

(19)

The virtual strain energy of the system is as follows:

$$\delta U={\displaystyle {\int }_{\Omega }{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left({\sigma }_{xx}^{(k)}\delta {\epsilon }_{xx}^{(k)}+{\sigma }_{yy}^{(k)}\delta {\epsilon }_{yy}^{(k)}+{\sigma }_{zz}^{(k)}\delta {\epsilon }_{zz}^{(k)}+{\tau }_{xy}^{(k)}\delta {\gamma }_{xy}^{(k)}+{\tau }_{xz}^{(k)}\delta {\gamma }_{xz}^{(k)}+{\tau }_{yz}^{(k)}\delta {\gamma }_{yz}^{(k)}\right)}}}\mathrm{d}z\mathrm{d}x\mathrm{d}y$$

(20)

The variation in work carried out by the external forces q is as follows:

$$\delta V=-{\displaystyle {\int }_{\Omega }{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}q\left(\delta w_b+\delta w_s+g\left(z\right)\delta w_z\right)}\mathrm{d}z\mathrm{d}x\mathrm{d}y}}$$

(21)

where q represents the external load in the transverse direction.

The virtual kinetic energy of the system is as follows:

$$\delta K={\displaystyle {\int }_{\Omega }{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}{\rho }_c^{(k)}}}(\dot{u}\delta \dot{u}+\dot{v}\delta \dot{v}+\dot{w}\delta \dot{w})}\mathrm{d}z\mathrm{d}x\mathrm{d}y$$

(22)

In this paper, the dot-superscript convention denotes the differential with respect to the time variable.

If the effect of thickness stretching is ignored, then:

$$\delta U={\displaystyle {\int }_{\Omega }{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left({\sigma }_{xx}^{(k)}\delta {\epsilon }_{xx}^{(k)}+{\sigma }_{yy}^{(k)}\delta {\epsilon }_{yy}^{(k)}+{\tau }_{xy}^{(k)}\delta {\gamma }_{xy}^{(k)}+{\tau }_{xz}^{(k)}\delta {\gamma }_{xz}^{(k)}+{\tau }_{yz}^{(k)}\delta {\gamma }_{yz}^{(k)}\right)}}}\mathrm{d}z\mathrm{d}x\mathrm{d}y$$

(23)

$$\delta V=-{\displaystyle {\int }_{\Omega }{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}q\left(\delta w_b+\delta w_s\right)}\mathrm{d}z\mathrm{d}x\mathrm{d}y}}$$

(24)

$$\delta K={\displaystyle {\int }_{\Omega }{\displaystyle \sum _{k=1}^{N_L}{\displaystyle {\int }_{z_k}^{z_{k+1}}{\rho }_c^{(k)}}}(\dot{u}\delta \dot{u}+\dot{v}\delta \dot{v}+\dot{w}\delta \dot{w})}\mathrm{d}z\mathrm{d}x\mathrm{d}y$$

(25)

Hamilton’s principle is as follows:

$${\displaystyle {\int }_0^T\delta \left(U+V-K\right)\mathrm{d}t=0}$$

(26)

Substituting the expressions from Equations (20)–(22) into Equation (26), we obtain the governing equations under IRPT as follows:

$$\begin{array}{l}{A}_{11}\frac{{\partial }^{2}u}{\partial x^2}+A_{66}\frac{{\partial }^{2}u}{\partial y^2}+\left(A_{12}+A_{66}\right)\frac{{\partial }^{2}v}{\partial x\partial y}-B_{11}\frac{{\partial }^3{w}_b}{\partial x^3}-\left(B_{12}+2B_{66}\right)\frac{{\partial }^3{w}_b}{\partial x\partial y^2}-\\ B_{11}^s\frac{{\partial }^3{w}_s}{\partial x^3}-\left(B_{12}^s+2B_{66}^s\right)\frac{{\partial }^3{w}_s}{\partial x\partial y^2}+X_{13}\frac{\partial w_z}{\partial x}=I_0\ddot{u}-I_1\frac{\partial {\ddot{w}}_b}{\partial x}-J_1\frac{\partial {\ddot{w}}_s}{\partial x}\end{array}$$

(27)

$$\begin{array}{l}{A}_{22}\frac{{\partial }^{2}v}{\partial y^2}+A_{66}\frac{{\partial }^{2}v}{\partial x^2}+\left(A_{12}+A_{66}\right)\frac{{\partial }^{2}u}{\partial x\partial y}-B_{22}\frac{{\partial }^3{w}_b}{\partial y^3}-\left(B_{12}+2B_{66}\right)\frac{{\partial }^3{w}_b}{\partial x^2\partial y}-\\ B_{22}^s\frac{{\partial }^3{w}_s}{\partial y^3}-\left(B_{12}^s+2B_{66}^s\right)\frac{{\partial }^3{w}_s}{\partial x^2\partial y}+X_{23}\frac{\partial w_z}{\partial y}=I_0\ddot{v}-I_1\frac{\partial {\ddot{w}}_b}{\partial y}-J_1\frac{\partial {\ddot{w}}_s}{\partial y}\end{array}$$

(28)

$$\begin{array}{l}{B}_{11}\frac{{\partial }^{3}u}{\partial x^3}+\left(B_{12}+2B_{66}\right)\left(\frac{{\partial }^{3}u}{\partial x\partial y^2}+\frac{{\partial }^{3}v}{\partial x^2\partial y}\right)+B_{22}\frac{{\partial }^{3}v}{\partial y^3}-D_{11}\frac{{\partial }^4{w}_b}{\partial x^4}-D_{22}\frac{{\partial }^4{w}_b}{\partial y^4}-\\ 2\left(D_{12}+2D_{66}\right)\frac{{\partial }^4{w}_b}{\partial x^2\partial y^2}-D_{11}^s\frac{{\partial }^4{w}_s}{\partial x^4}-D_{22}^s\frac{{\partial }^4{w}_s}{\partial y^4}-2\left(D_{12}^s+2D_{66}^s\right)\frac{{\partial }^4{w}_s}{\partial x^2\partial y^2}+\\ Y_{13}\frac{{\partial }^2{w}_z}{\partial x^2}+Y_{23}\frac{{\partial }^2{w}_z}{\partial y^2}+q=I_0\left({\ddot{w}}_b+{\ddot{w}}_s\right)+J_0{\ddot{w}}_z+I_1\left(\frac{\partial \ddot{u}}{\partial x}+\frac{\partial \ddot{v}}{\partial y}\right)-I_2{\nabla }^2{\ddot{w}}_b-J_2{\nabla }^2{\ddot{w}}_s\end{array}$$

(29)

$$\begin{array}{l}{B}_{11}^s\frac{{\partial }^{3}u}{\partial x^3}+\left(B_{12}^s+2B_{66}^s\right)\left(\frac{{\partial }^{3}u}{\partial x\partial y^2}+\frac{{\partial }^{3}v}{\partial x^2\partial y}\right)+B_{22}^s\frac{{\partial }^{3}v}{\partial y^3}-D_{11}^s\frac{{\partial }^4{w}_b}{\partial x^4}-\\ D_{22}^s\frac{{\partial }^4{w}_b}{\partial y^4}-2\left(D_{12}^s+2D_{66}^s\right)\frac{{\partial }^4{w}_b}{\partial x^2\partial y^2}-H_{11}^s\frac{{\partial }^4{w}_s}{\partial x^4}-H_{22}^s\frac{{\partial }^4{w}_s}{\partial y^4}-\\ 2\left(H_{12}^s+2H_{66}^s\right)\frac{{\partial }^4{w}_s}{\partial x^2\partial y^2}+A_{55}^s\frac{{\partial }^2{w}_s}{\partial x^2}+A_{44}^s\frac{{\partial }^2{w}_s}{\partial y^2}+\left(Y_{13}^s+A_{55}^s\right)\frac{{\partial }^2{w}_z}{\partial x^2}+\\ \left(Y_{23}^s+A_{44}^s\right)\frac{{\partial }^2{w}_z}{\partial y^2}+q+N_{xx}^0\left(\frac{{\partial }^2{w}_b}{\partial x^2}+\frac{{\partial }^2{w}_s}{\partial x^2}\right)+N_{yy}^0\left(\frac{{\partial }^2{w}_b}{\partial y^2}+\frac{{\partial }^2{w}_s}{\partial y^2}\right)=\\ I_0\left({\ddot{w}}_b+{\ddot{w}}_s\right)+J_0{\ddot{w}}_z+J_1\left(\frac{\partial \ddot{u}}{\partial x}+\frac{\partial \ddot{v}}{\partial y}\right)-J_2{\nabla }^2{\ddot{w}}_b-K_2{\nabla }^2{\ddot{w}}_s\end{array}$$

(30)

$$\begin{array}{l}-X_{13}\frac{\partial u}{\partial x}-X_{23}\frac{\partial v}{\partial y}+Y_{13}\frac{{\partial }^2{w}_b}{\partial x^2}+Y_{23}\frac{{\partial }^2{w}_b}{\partial y^2}+\left(Y_{13}^s+A_{55}^s\right)\frac{{\partial }^2{w}_s}{\partial x^2}+\left(Y_{23}^s+A_{44}^s\right)\frac{{\partial }^2{w}_s}{\partial y^2}+\\ A_{55}^s\frac{{\partial }^2{w}_z}{\partial x^2}+A_{44}^s\frac{{\partial }^2{w}_z}{\partial y^2}-Z_{33}{w}_z+gq+g{N}_{xx}^0\left(\frac{{\partial }^2{w}_b}{\partial x^2}+\frac{{\partial }^2{w}_s}{\partial x^2}\right)+g{N}_{yy}^0\left(\frac{{\partial }^2{w}_b}{\partial y^2}+\frac{{\partial }^2{w}_s}{\partial y^2}\right)\\ =J_0\left({\ddot{w}}_b+{\ddot{w}}_s\right)+K_0{\ddot{w}}_z\end{array}$$

(31)

where

$$\left(A_{ij},A_{ij}^s,B_{ij},B_{ij}^s,D_{ij},D_{ij}^s,H_{ij}^s\right)={\displaystyle \sum _{k=1}^{NL}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left(1,g^2,z,f,z^2,fz,f^2\right)}{C}_{ij}{}^{\left(k\right)}dz}$$

(32)

$$\left(X_{ij},Y_{ij},Y_{ij}^s,Z_{ij}\right)={\displaystyle \sum _{k=1}^{NL}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left(g^{\prime },g^{\prime }z,g^{\prime }f,{g^{\prime }}^2\right)}{C}_{ij}{}^{\left(k\right)}dz}$$

(33)

$$\left(I_0,I_1,I_2\right)={\displaystyle \sum _{k=1}^{NL}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left(1,z,z^2\right)}{\rho }^{\left(k\right)}dz}$$

(34)

$$\left(J_0,J_1,J_2\right)={\displaystyle \sum _{k=1}^{NL}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left(g,f,zf\right)}{\rho }^{\left(k\right)}dz}$$

(35)

$$\left(K_0,K_2\right)={\displaystyle \sum _{k=1}^{NL}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left(g^2,f^2\right)}{\rho }^{\left(k\right)}dz}$$

(36)

4. Closed Solutions

Consider an FG-GRC plate with four sides simply supported, as shown in Figure 1. The transversal load q acts on the surface of the FG-GRC plate. The Navier solution is used to solve the closed solution of Equations (27)–(31). According to the Navier solution, the solution form is assumed to be as follows:

$$\begin{array}{l}u(x,y,t)={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{s=1}^S\left(U_{rs}{e}^{i\omega t}\cos \alpha x\sin \beta y\right)}}\hfill \\ v(x,y,t)={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{s=1}^S\left(V_{rs}{e}^{i\omega t}\sin \alpha x\cos \beta y\right)}}\hfill \\ w_b(x,y,t)={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{s=1}^S\left(W_{brs}{e}^{i\omega t}\sin \alpha x\sin \beta y\right)}}\hfill \\ \begin{array}{l}{w}_s(x,y,t)={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{s=1}^S\left(W_{srs}{e}^{i\omega t}\sin \alpha x\sin \beta y\right)}}\\ w_z(x,y,t)={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{s=1}^S\left(W_{zrs}{e}^{i\omega t}\sin \alpha x\sin \beta y\right)}}\end{array}\hfill \end{array}$$

(37)

where $U_{rs}$, $V_{rs}$, $W_{brs}$, $W_{srs}$, and $W_{zrs}$ are unknown coefficients, while r and s are the half-wave numbers in the x and y directions, respectively. ω is the natural frequency. In addition,

$$i=\sqrt{-1},\alpha =\frac{r\pi }{a},\beta =\frac{s\pi }{b}$$

(38)

The form of transversal load q is assumed to be as follows:

$$q(x,y)={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{s=1}^S{Q}_{rs}}}\sin \alpha x\sin \beta y$$

(39)

For the bending analysis in the following, it is necessary to consider a uniformly distributed transverse load $q\left(x,y\right)=F_0$, in which case [43],

$$Q_{rs}=4{\lambda }_{rs}{F}_0/(E_{\mathrm{M}}rs{\pi }^2)$$

(40)

where

$${\lambda }_{rs}=1-{\left(-1\right)}^r-{\left(-1\right)}^s+{\left(-1\right)}^{r+s}$$

(41)

Substituting Equations (37) and (39) into Equations (27)–(31), the closed solutions can be derived from the following equation.

$$\left(\left[\begin{array}{ccccc}{j}_{11}& j_{12}& j_{13}& j_{14}& j_{15}\\ j_{12}& j_{22}& j_{23}& j_{24}& j_{25}\\ j_{13}& j_{23}& j_{33}& j_{34}& j_{35}\\ j_{14}& j_{24}& j_{34}& j_{44}& j_{45}\\ j_{15}& j_{25}& j_{35}& j_{45}& j_{55}\end{array}\right]-{\omega }^2\left[\begin{array}{ccccc}{m}_{11}& 0& m_{13}& m_{14}& 0\\ 0& m_{22}& m_{23}& m_{24}& 0\\ m_{31}& m_{32}& m_{33}& m_{34}& m_{35}\\ m_{41}& m_{42}& m_{43}& m_{44}& m_{45}\\ 0& 0& m_{53}& m_{54}& m_{55}\end{array}\right]\right)\left\{\begin{array}{c}{U}_{rs}\\ V_{rs}\\ W_{brs}\\ W_{srs}\\ W_{zrs}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ Q_{rs}\\ Q_{rs}\\ 0\end{array}\right\}$$

(42)

where

$$\begin{array}{l}{j}_{11}=A_{11}{\alpha }^2+A_{66}{\beta }^2,j_{12}=\left(A_{12}+A_{66}\right)\alpha \beta ,j_{13}=-B_{11}{\alpha }^3-\left(B_{12}+2B_{66}\right)\alpha {\beta }^2\\ j_{14}=-B_{11}^s{\alpha }^3-\left(B_{12}^s+2B_{66}^s\right)\alpha {\beta }^2,j_{15}=-X_{13}\alpha ,j_{22}=A_{66}{\alpha }^2+A_{22}{\beta }^2\\ j_{23}=-B_2{\beta }^3-\left(B_{12}+2B_{66}\right){\alpha }^2\beta ,j_{24}=-B_{22}^s{\beta }^3-\left(B_{12}^s+2B_{66}^s\right){\alpha }^2\beta ,j_{25}=-X_{23}\beta \\ j_{33}=D_{11}{\alpha }^4+2\left(D_{12}+2D_{66}\right){\alpha }^2{\beta }^2+D_{22}{\beta }^4,j_{45}=\left(Y_{13}^s+A_{55}^s\right){\alpha }^2+\left(Y_{23}^s+A_{44}^s\right){\beta }^2\\ j_{34}=D_{11}^s{\alpha }^4+2\left(D_{12}^s+2D_{66}^s\right){\alpha }^2{\beta }^2+D_{22}^s{\beta }^4,j_{55}=A_{55}^s{\alpha }^2+A_{44}^s{\beta }^2+Z_{33}\\ j_{35}=Y_{13}{\alpha }^2+Y_{23}{\beta }^2,j_{44}=H_{11}^s{\alpha }^4+2\left(H_{12}^s+2H_{66}^s\right){\alpha }^2{\beta }^2+H_{22}^s{\beta }^4+A_{55}^s{\alpha }^2+A_{44}^s{\beta }^2\\ m_{11}=I_0,m_{13}=-\alpha I_1,m_{14}=-\alpha J_1,m_{22}=I_0,m_{23}=-\beta I_1,m_{24}=-\beta J_1\\ m_{33}=I_0+I_2\left({\alpha }^2+{\beta }^2\right),m_{34}=I_0+J_1\left({\alpha }^2+{\beta }^2\right),m_{35}=J_0,m_{44}=I_0+K_2\left({\alpha }^2+{\beta }^2\right)\\ m_{44}=I_0+K_2\left({\alpha }^2+{\beta }^2\right),m_{45}=J_0,m_{55}=K_0\end{array}$$

To solve for the bending deflection, just take ω = 0, and we can find $U_{rs}$, $V_{rs}$, $W_{brs}$, $W_{srs}$, and $W_{zrs}$, and then the deflection value using Equation (42). The free vibration frequency can be found by simply substituting Equation (42) for $Q_{rs}=0$.

If the effect of thickness stretching is neglected and the same procedure as IRPT is used, the closed solutions of RPT can be obtained from the following equation.

$$\left(\left[\begin{array}{cccc}{j}_{11}& j_{12}& j_{13}& j_{14}\\ j_{12}& j_{22}& j_{23}& j_{24}\\ j_{13}& j_{23}& j_{33}& j_{34}\\ j_{14}& j_{24}& j_{34}& j_{44}\end{array}\right]-{\omega }^2\left[\begin{array}{cccc}{m}_{11}& 0& m_{13}& m_{14}\\ 0& m_{22}& m_{23}& m_{24}\\ m_{13}& m_{23}& m_{33}& m_{34}\\ m_{14}& m_{24}& m_{34}& m_{44}\end{array}\right]\right)\left\{\begin{array}{c}{U}_{rs}\\ V_{rs}\\ W_{brs}\\ W_{srs}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ Q_{rs}\\ Q_{rs}\end{array}\right\}$$

(43)

where

$$\begin{array}{l}{j}_{11}=A\left({\alpha }^2+\frac{1-v_c^{\left(k\right)}}{2}{\beta }^2\right),j_{12}=\frac{1+v_c^{\left(k\right)}}{2}A\alpha \beta ,j_{22}=A\left(\frac{1-v_c^{\left(k\right)}}{2}{\alpha }^2+{\beta }^2\right)\\ j_{13}=-B\alpha \left({\alpha }^2+{\beta }^2\right),j_{14}=-B^s\alpha \left({\alpha }^2+{\beta }^2\right),j_{23}=-B\beta \left({\alpha }^2+{\beta }^2\right)\\ j_{24}=-B^s\beta \left({\alpha }^2+{\beta }^2\right),j_{33}=D{\left({\alpha }^2+{\beta }^2\right)}^2,j_{34}=D^s{\left({\alpha }^2+{\beta }^2\right)}^2\\ j_{44}=H{\left({\alpha }^2+{\beta }^2\right)}^2+A^s\left({\alpha }^2+{\beta }^2\right),m_{11}=m_{22}=I_0,m_{13}=-\alpha I_1,m_{14}=-\alpha J_1,m_{23}=-\beta I_1\\ m_{24}=-\beta J_1,m_{33}=I_0+I_2\left({\alpha }^2+{\beta }^2\right),m_{34}=I_0+J_2\left({\alpha }^2+{\beta }^2\right)\\ m_{44}=I_0+K_2\left({\alpha }^2+{\beta }^2\right)\\ \left(A,B,B^s,D,D^s,H\right)={\displaystyle \sum _{k=1}^{NL}{\displaystyle {\int }_{z_k}^{z_{k+1}}\left(1,z,f,z^2,zf,f^2\right)\frac{E_c^{\left(k\right)}}{1-v_c^{\left(k\right)}{}^2}}dz},A^s={\displaystyle \sum _{k=1}^{NL}{\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{g^2{E}_c^{\left(k\right)}}{2(1+v_c^{\left(k\right)})}}dz}\end{array}$$

The process of solving RPT is similar to solving IRPT.

5. Results and Discussion

In this section, the bending and free vibration behaviors of the FG-GRC plates are investigated in detail, with special attention on the effects of the distribution pattern of GPLs, weight fraction of GPLs, length-to-thickness ratio of FG-GRC plates, and a total number of layers $N_L$ on the mechanical behaviors of FG-GRC plates when the effect of thickness stretching are considered. It should be noted that since the reinforcing materials can be considered as uniformly distributed in the horizontal direction, several composite plates mentioned below (FG-GRC plates, Al/(ZrO₂)-1 plates and Al/Al₂O₃ plates) can be considered as quasi-isotropic material plates, and their material properties only produce gradient changes in the transverse thickness direction.

5.1. Static Bending

A functionally graded square plate composed of $\mathrm{Al}/\left({\mathrm{ZrO}}_2\right)\text{-}1$ is used as an example to verify the effectiveness of the current method. The material properties are shown below

$$\begin{array}{l}\mathrm{Al}:E_1=70\mathrm{GPa},{\nu }_1=0.3\\ \left({\mathrm{ZrO}}_2\right)\text{-}1:E_2=200\mathrm{GPa},{\nu }_2=0.3\end{array}$$

The plate has Al on the upper surface and $\left({\mathrm{ZrO}}_2\right)\text{-}1$ on the lower surface with an intermediate gradient distribution. The effective material properties are shown below [46]:

$$P(z)=P_1+\left(P_2-P_1\right){\left(\frac{1}{2}+\frac{z}{h}\right)}^n\left(n\ge 0\right)$$

(44)

where $P_1$ and $P_2$ represent the material properties of Al and $\left({\mathrm{ZrO}}_2\right)\text{-}1$, such as Young’s modulus, density, etc., respectively. n is the volume fraction index. The functionally graded plate is divided into $N_L$ layers of isotropic plates along the thickness direction. The total material properties of the functionally graded plates are superimposed on the material properties of each layer in the thickness direction, and the equivalent material property $P_{\mathrm{eq}}^{\left(k\right)}$ for the kth layer is defined as follows [47]:

$$P_{\mathrm{eq}}^{(k)}={\displaystyle {\int }_{z_k}^{z_{k+1}}\frac{P(z)}{z_{k+1}-z_k}} \mathrm{d}z,\left(k=1,2,\cdots ,N_L\right)$$

(45)

The dimensionless deflection is defined as follows and the uniformly distributed transverse load $F_0=500 \mathrm{kPa}$.

$$w_c=100w{E}_1{h}^3/\left[12\left(1-{\nu }_1^2\right)F_0{a}^4\right]$$

(46)

The effects of different truncation levels R and S and $N_L$ on the dimensionless deflection are given in

Table 1. Dimensionless central deflection values of $\mathrm{Al}/\left({\mathrm{ZrO}}_2\right)\text{-}1$ functionally graded square plates at different truncation levels. (n = 0.5).

		R = S = 1	R = S = 10	R = S = 20
NL = 10	RPT	0.2402	0.2313	0.2312
	IRPT	0.2550	0.2462	0.2462
NL = 20	RPT	0.2407	0.2318	0.2317
	IRPT	0.2556	0.2467	0.2467
NL = 30	RPT	0.2408	0.2319	0.2318
	IRPT	0.2557	0.2468	0.2468

. n is the volume fraction index. The results show that the dimensionless deflection is close to convergence when the truncation levels R and S of Equation (37) are taken as 10 for the different total numbers of layers. Therefore, the use of R = S = 10 for all subsequent bending deflection-related calculations balances accuracy and calculation costs.

The dimensionless central deflection values for $\mathrm{Al}/\left({\mathrm{ZrO}}_2\right)\text{-}1$ functionally graded square plates at $a/h=5$, $N_L=10$, and $R=S=10$ are given in

Table 2. Dimensionless central deflection values of $\mathrm{Al}/\left({\mathrm{ZrO}}_2\right)\text{-}1$ functionally graded square plates under a transverse uniform load. ($R=S=10,N_L=10,a/h=5$).

	n
	0	0.5	1	2
Ref. [48]	0.1703	0.2232	0.2522	0.2827
Ref. [49]	0.1671	0.2505	0.2905	0.3280
RPT	0.1716	0.2313	0.2715	0.3140
IRPT	0.1830	0.2462	0.2883	0.3327

, where n is the volume fraction index. By comparing with the finite element method [48] and the meshless local Petrov–Galerkin method [49], it can be seen that RPT and IRPT are in good agreement with the available literature in terms of results.

Next, consider an FG-GRC plate. The epoxy is selected for the matrix material, and its material properties are shown below.

$$\begin{array}{l}\mathrm{epoxy}:E_{\mathrm{M}}=3.0\mathrm{GPa},{\nu }_{\mathrm{M}}=0.34,{\rho }_{\mathrm{M}}=1.2{\mathrm{g}/\mathrm{cm}}^3\\ \mathrm{GPLs}:E_{\mathrm{G}}=1.01\mathrm{TPa},{\nu }_{\mathrm{G}}=0.186,{\rho }_{\mathrm{G}}=1.06{\mathrm{g}/\mathrm{cm}}^3\end{array}$$

The geometric dimensions of GPLs are: $l_{\mathrm{G}}=2.5\mathsf{\mu }\mathrm{m},w_{\mathrm{G}}=1.5\mathsf{\mu }\mathrm{m},h_{\mathrm{G}}=1.5\mathrm{nm}$. If not specified, the default is taken as $f_{\mathrm{G}}=1.0\%$ and $N_L=10$.

Figure 3 and Figure 4 investigated the effects of the total number of layers $N_L$ and the GPLs weight fraction $f_{\mathrm{G}}$ on the percentage deflection ratio $w_c/w_{\mathrm{M}}$ of the FG-GRC square plates under four different distribution patterns of GPLs using the plate theories IRPT and RPT, respectively. $w_c$ is the center deflection of FG-GRC plates and $w_{\mathrm{M}}$ is the center deflection of plates without GPLs reinforcement (matrix material plates).

It can be clearly seen that when the GPLs are asymmetrically distributed (FG-A), the difference between the results obtained from the two plate theories (RPT and IRPT) is very obvious. When the GPLs are symmetrically distributed within the FG-GRC plates (UD, FG-O, and FG-X), the difference between the results calculated by the two plate theories is very small. This is due to the asymmetric distribution of GPLs in the transverse thickness direction, resulting in a more significant transverse anisotropy in FG-A-type FG-GRC plates. That is to say, the more significant transverse anisotropy of the FG-A-type FG-GRC plates compared to the symmetrically distributed modes of GPLs leads to a significant difference in its deflection ratio in the calculated results of RPT and IRPT, and this difference is evident in both thick plates with $a/h=4$ and moderately thick plates with $a/h=10$.

It can also be seen from Figure 3 and Figure 4 that the difference in deflection ratio $w_c/w_{\mathrm{M}}$ decreases with decreasing a/h values for different distribution modes, and the effect of thickness stretching tends to increase as the total number of layers $N_L$ and the weight fraction $f_{\mathrm{G}}$ of GPLs increase. However, as the total number of layers $N_L$ and the weight fraction $f_{\mathrm{G}}$ of GPLs increase, the difference due to the effect of thickness stretching tends to stabilize.

5.2. Free Vibration

Firstly, metal-ceramic functionally graded materials are calculated to verify the validity of current methods. A simply supported ${\mathrm{Al}/\mathrm{Al}}_2{\mathrm{O}}_3$ functionally graded square plate with a total number of layers $N_L=10$ is considered. The material properties of $\mathrm{Al}$ and ${\mathrm{Al}}_2{\mathrm{O}}_3$ are shown in the following equations.

$$\begin{array}{l}\mathrm{Al}:E_1=70.0\times {10}^9{\mathrm{N}/\mathrm{m}}^2,{\nu }_1=0.3,{\rho }_1=2702{\mathrm{kg}/\mathrm{m}}^3\\ {\mathrm{Al}}_2{\mathrm{O}}_3:E_2=380\times {10}^9{\mathrm{N}/\mathrm{m}}^2,{\nu }_2=0.3,{\rho }_2=3800{\mathrm{kg}/\mathrm{m}}^3\end{array}$$

The effective material properties of the functionally graded plates are determined by Equations (45) and (46), which follow the same procedure as when bending is considered. n is the volume fraction index. The dimensionless natural frequency of the FG-GRC plates is defined by the following equation.

$${\omega }_c=\omega h\sqrt{{\rho }_2/E_2}$$

(47)

As shown in

Table 3. Dimensionless fundamental frequencies of ${\mathrm{Al}/\mathrm{Al}}_2{\mathrm{O}}_3$ functionally graded t square plates.

$\mathit{a}/\mathit{h}$		n
		0.0	0.5	1.0	4.0	10.0
2	Ref. [50]	0.9400	0.8232	0.7476	0.5994	0.5460
	RPT	0.9297	0.8126	0.7366	0.5919	0.5393
	IRPT	0.9420	0.8267	0.7529	0.6053	0.5479
5	Ref. [50]	0.2121	0.1819	0.1640	0.1383	0.1306
	RPT	0.2113	0.1811	0.1634	0.1375	0.1293
	IRPT	0.2122	0.1829	0.1662	0.1406	0.1310
10	Ref. [50]	0.05777	0.04917	0.04426	0.03811	0.03642
	RPT	0.05769	0.04914	0.04427	0.03795	0.03611
	IRPT	0.05777	0.04950	0.04495	0.03878	0.03656

, the results of RPT and IRPT are in good agreement with those of higher-order shear deformation plate theory [50].

Next, we will analyze the FG-GRC plates. The material parameters of the epoxy matrix materials and GPLs are the same as those in the bending analysis. Where the dimensionless natural frequencies are defined as shown below.

$${\omega }_c=\omega h\sqrt{{\rho }_{\mathrm{M}}/E_{\mathrm{M}}}$$

(48)

in which ${\omega }_c$ and ${\omega }_{\mathrm{M}}$ are the natural frequencies of FG-GRC square plates and pure epoxy square plates, respectively.

The dimensionless natural frequencies and relative frequency increase $\left({\omega }_c-{\omega }_{\mathrm{M}}\right)/{\omega }_{\mathrm{M}}$ of FG-GRC plates are given in

Table 4. Dimensionless natural frequencies of pure epoxy square plates and FG-GRC square plates with various mode shapes. ($N_L=10,f_{\mathrm{G}}=1.0\%$).

$\mathit{a}/\mathit{h}$	r, s	Theories	Pure Epoxy	UD	FG-O	FG-X	FG-A
4	1, 1	RPT	0.3146	0.6548 (108.1%)	0.5728 (82.1%)	0.6967 (121.5%)	0.6093 (93.7%)
		IRPT	0.3169	0.6595 (108.1%)	0.5784 (82.5%)	0.7006 (121.1%)	0.6250 (97.2%)
	2, 1	RPT	0.6584	1.3706 (108.2%)	1.2433 (88.8%)	1.3941 (111.7%)	1.2904 (96.0%)
		IRPT	0.6668	1.3879 (108.1%)	1.2663 (90.0%)	1.4069 (111.0%)	1.3272 (99.0%)
	2, 2	RPT	0.9294	1.9350 (108.2%)	1.7918 (92.8%)	1.9227 (106.9%)	1.8352 (97.5%)
		IRPT	0.9441	1.9654 (108.2%)	1.8352 (94.4%)	1.9444 (106.0%)	1.8903 (100.2%)
	3, 1	RPT	1.0862	2.2615 (108.2%)	2.1145 (94.7%)	2.2241 (104.8%)	2.1530 (98.2%)
		IRPT	1.1050	2.3003 (108.2%)	2.1719 (96.6%)	2.2514 (103.7%)	2.2189 (100.8%)
	3, 2	RPT	1.2970	2.7006 (108.2%)	2.5528 (96.8%)	2.6267 (102.5%)	2.5826 (99.1%)
		IRPT	1.3216	2.7514 (108.2%)	2.6309 (99.1%)	2.6621 (101.4%)	2.6633 (101.5%)
	3, 3	RPT	1.6034	3.3389 (108.2%)	3.1957 (99.3%)	3.2092 (100.1%)	3.2108 (100.2%)
		IRPT	1.6368	3.4079 (108.2%)	3.3068 (102.0%)	3.2577 (99.0%)	3.3130 (102.4%)
10	1, 1	RPT	0.0584	0.1216 (108.2%)	0.1023 (75.2%)	0.1366 (133.9%)	0.1118 (91.4%)
		IRPT	0.0585	0.1218 (108.2%)	0.1025 (75.2%)	0.1368 (133.8%)	0.1144 (95.6%)
	2, 1	RPT	0.1391	0.2895 (108.1%)	0.2470 (77.6%)	0.3188 (129.2%)	0.2674 (92.2%)
		IRPT	0.1396	0.2906 (108.2%)	0.2481 (77.7%)	0.3198 (129.1%)	0.2737 (96.1%)
	2, 2	RPT	0.2132	0.4437 (108.1%)	0.3828 (79.5%)	0.4809 (125.6%)	0.4111 (92.8%)
		IRPT	0.2143	0.4460 (108.1%)	0.3854 (79.8%)	0.4829 (125.3%)	0.4212 (96.5%)
	3, 1	RPT	0.2596	0.5402 (108.1%)	0.4691 (80.7%)	0.5803 (123.5%)	0.5015 (93.2%)
		IRPT	0.2612	0.5435 (108.1%)	0.4729 (81.0%)	0.5831 (123.2%)	0.5141 (96.8%)
	3, 2	RPT	0.3253	0.6770 (108.1%)	0.5930 (82.3%)	0.7190 (121.0%)	0.6302 (93.7%)
		IRPT	0.3277	0.6820 (108.1%)	0.5990 (82.8%)	0.7231 (120.7%)	0.6465 (97.3%)
	3, 3	RPT	0.4265	0.8877 (108.1%)	0.7870 (84.5%)	0.9284 (117.7%)	0.8294 (94.5%)
		IRPT	0.4304	0.8958 (108.1%)	0.7971 (85.2%)	0.9348 (117.2%)	0.8517 (97.9%)

Note: The values in parentheses indicate the relative frequency increase $\left({\omega }_c-{\omega }_{\mathrm{M}}\right)/{\omega }_{\mathrm{M}}$, in which ${\omega }_c$ and ${\omega }_{\mathrm{M}}$ are the natural frequencies of FG-GRC square plates and pure epoxy square plates, respectively.

. Each distribution mode of the FG-GRC square plates was calculated. The results of the calculations for the pure epoxy plates are also included in Table 4. The dimensionless natural frequencies and relative frequency increases in RPT and IRPT were compared. The role played by the effect of thickness stretching in the free vibration process of FG-GRC plates can be seen by comparing the results of RPT and IRPT calculations.

As seen in Table 4, the RPT underestimates the dimensionless natural frequency of the FG-GRC plate compared to the IRPT because the effect of thickness stretching is neglected. The calculated results for FG-A-type FG-GRC plates with asymmetric distribution of GPLs, including the dimensionless natural frequencies and relative frequency increases, are significantly affected by the effect of thickness stretching. This is because the transverse anisotropy of the FG-GRC plates with FG-A-type distribution is more significant than the other three distributions. It can also be seen from Table 4 that the FG-GRC plate with FG-X-type distribution has the best ability to improve the free vibration response of the plates among the four distribution modes. The UD-type distribution of FG-GRC plates is almost unaffected by the effect of thickness stretching, which is consistent with our expected results due to the isotropic nature of the UD-type distribution in the transverse displacement direction. We can also conclude that when the value of a/h decreases, that is, when the FG-GRC plates become thicker, the free vibration frequency of the FG-GRC plates increases significantly.

The effects of the total number of layers $N_L$ and the weight fraction $f_{\mathrm{G}}$ of GPLs on the fundamental frequency change in FG-GRC square plates at $a/h=4$ and $a/h=10$ are given in Figure 5 and Figure 6, respectively. The relative frequency change is equal to $\left({\omega }_c-{\omega }_{\mathrm{M}}\right)/{\omega }_{\mathrm{M}}$. Figure 5 is calculated by taking $f_{\mathrm{G}}=1.0\%$ and Figure 6 is calculated by taking $N_L=10$. All four distribution modes of the FG-GRC plate were considered and the results of the RPT and IRPT calculations are compared in Figure 5 and Figure 6.

Among the four distribution modes, GPLs symmetrically distributed UD-type, FG-O-type, and FG-X-type FG-GRC plates showed little difference in the calculated results under the two plate theories (RPT and IRPT). The calculated results for the FG-A-type FG-GRC plate with asymmetric distribution of GPLs are significantly different under RPT and IRPT. Compared with the IRPT results, the RPT results without considering the effect of thickness stretching are significantly lower, which is caused by the more significant anisotropy in the transverse displacement direction of the FG-A-type distribution than the other three GPLs symmetrically distributed FG-GRC plates. From Figure 5 and Figure 6, we can also see that the effect of thickness stretching becomes more significant as the total number of layers ($N_L$) and the weight fraction ($f_{\mathrm{G}}$) of GPLs increase. From Figure 3, Figure 4, Figure 5 and Figure 6, it can be seen that for the calculated results of bending deflection and free vibration fundamental frequencies of FG-GRC plates, the difference in results between the four distribution modes of GPLs decreases significantly as the value of a/h decreases.

6. Conclusions

Based on two types of plate theories considering (IRPT) and not considering (RPT) the thickness stretching effect, the effect of thickness stretching on the bending and free vibration behaviors of FG-GRC plates was investigated. The governing equations of the FG-GRC plates are derived by Hamilton’s principle. The Navier solution method is used to obtain the closed solutions for the percentage deflection ratio, dimensionless free frequency, and the relative frequency change in the FG-GRC rectangular simply supported plates. The effects of $a/h$ value, the total number of layers $N_L$, and weight fraction $f_{\mathrm{G}}$ of GPLs are investigated by comparing the two plate theories. It is found that:

(1) The effect of thickness stretching depends on the transverse anisotropy of FG-GRC plates. FG-GRC plates with asymmetric distribution of GPLs are most significantly affected by the effect of thickness stretching. Even for the moderately thick FG-GRC plates ($a/h=10$) with asymmetric distribution of GPLs, the effect of thickness stretching cannot be neglected when considering their bending and free vibration behaviors.

(2) When the value of $a/h$ decreases, the dimensionless natural frequencies of the FG-GRC plates increase significantly.

(3) The difference in bending and free vibration of FG-GRC plates due to different distribution patterns of GPLs decreases significantly as the value of a/h decreases.

(4) The effect of thickness stretching becomes more obvious as the total number of layers of FG-GRC plates ($N_L$) and the weight fraction of GPLs ($f_{\mathrm{G}}$) increase.