Theoretical and Numerical Solution for the Bending and Frequency Response of Graphene Reinforced Nanocomposite Rectangular Plates

Abstract

In this work, we study the vibration and bending response of functionally graded graphene platelets reinforced composite (FG-GPLRC) rectangular plates embedded on different substrates and thermal conditions. The governing equations of the problem along with boundary conditions are determined by employing the minimum total potential energy and Hamilton’s principle, within a higher-order shear deformation theoretical setting. The problem is solved both theoretically and numerically by means of a Navier-type exact solution and a generalized differential quadrature (GDQ) method, respectively, whose results are successfully validated against the finite element predictions performed in the commercial COMSOL code, and similar outcomes available in the literature. A large parametric study is developed to check for the sensitivity of the response to different foundation properties, graphene platelets (GPL) distribution patterns, volume fractions of the reinforcing phase, as well as the surrounding environment and boundary conditions, with very interesting insights from a scientific and design standpoint.

1. Introduction

Due to their outstanding thermal and mechanical properties, carbon-based nano-filler reinforced composites are widely applied in many engineering fields, such as civil, biomedical and automotive engineering [1,2,3,4,5,6]. In more detail, graphene platelets (GPLs) are increasingly introduced as carbon nano-fillers because of their relevant potentials in terms of high surface area, elasticity modulus, thermal conductivity, etc. GPLs, as one of novel nanosize reinforcements, have special properties, and their two-dimensional geometry enables them to be scattered in the matrix with less agglomeration, unlike the one-dimensional anisotropic ones. Due to their excellent mechanical, chemical, and physical properties, graphene-based composites demonstrate a wide range of applications in an engineering field, such as sensors, fuel cells, supercapacitors, and batteries. The addition of graphene as reinforcing agent in a polymer matrix, indeed, improves the overall performances and properties of composite materials, as largely demonstrated in the literature from researchers working in this area [7,8]. The primary interest of using graphene materials stems from their excellent mechanical, thermal, electrical and physicochemical properties with prosing results in all fields of technologies. For example, graphene represents one of the stiffest and most grounded materials, with an elastic modulus of $\sim 1\mathrm{TPa}$ and quality of $\sim 100\mathrm{GPa}$ [9,10,11]. By introducing 1 volume percent of graphene in a polymer matrix, the nanocomposite material reaches a conductivity of about $0.1\hspace{0.17em}{\mathrm{Sm}}^{-1}$ with adequate consequences for electrical applications, along with significant changes in quality and strength [12]. In such a context, several theories and computational models have been developed in the last decades in the field of GPL-reinforced media. Anamagh and Bediz [13], for example, studied the buckling and vibration response of GPL-reinforced rectangular plates with different boundary conditions, based on a spectral-Tchebychev model. Reddy et al. [14] surveyed the vibrational frequencies of the composite plates reinforced by GPLs, and investigated the effect of various parameters, primarily, boundary conditions, distribution patterns, geometry and weight fractions of GPLs, on the natural frequencies of the system. In addition, Qaderi and Ebrahimi [15] focused on the frequency response of GPL reinforced rectangular plates embedded on viscoelastic substrates, and their sensitivity to different damping coefficients. In line with the previous works, Song et al. [16] studied the free and forced vibration behavior of FG-GPL-reinforced (FG-GPLR) plates by applying a first-order shear deformation theory (FSDT), with a clear enhancement of the vibration performances even with the addition of small quantities of GPLs. Based on the Chebyshev—Ritz procedure, Yang et al. [17] investigated the natural frequencies and critical buckling loads of FG-GPLR nanocomposite plates in presence of different porosities levels. Among the recent literature, different continuum-based nonlocal models have been considered as effective methods to treat plate-like nanostructures and to avoid possible difficulties encountered during experimental characterizations or time-consuming computational atomistic simulations of nanotubes. In this context, some theoretical studies of the free vibration response of graphene sheets can be found in the recent works [18,19,20,21,22,23,24], based on different nonlocal theoretical assumptions, accounting for different small-scale parameters, geometrical properties, boundary and environmental conditions. It is also well-known that different substrates can surround a structural member, thus affecting its mechanical behavior and stability. Numerous engineering problems (e.g., heavy machines, pavement of roads, etc.), indeed, are modeled as structural members resting on an elastic medium [25]. The elastic substrates are commonly modeled as Winkler or Pasternak foundations by means of one or two parameters [26,27]. The effect of visco-Pasternak substrate on the nonlinear dynamic response of the FG-GPLRC rectangular plates can be found in the seminal works by Fan et al. [28], and by Liu et al. [29] along with a sensitivity study of the mechanical behavior to different foundation parameters and porosity distributions. Among further works, Gao et al. [30] analyzed the nonlinear vibrational frequencies of FG-GPLR porous plates embedded on a two-parameter-type elastic medium, where an increased porosity coefficient was found to reduce the overall stiffness of structures. The vibrational properties of FG rectangular plates resting on a two-parameter elastic substrate were also surveyed by Thai and Choi [31]. They demonstrated that an increased quantity of metal components can significantly increase the deformability in a structural system. Similarly, Zhou et al. [32] studied the frequency response of thick plates on elastic media, while checking for the effect of different parameters, namely, the foundation coefficients, boundary conditions and aspect ratios, on the structural stiffness. A FSDT was also proposed in [33] to assess the nonlinear vibrational frequency and dynamic behavior of FG-GPLR plates resting on a viscoelastic-Pasternak foundation, with a clear reduction of the structural capacity for increased compressive loads.

Starting with the available literature on the topic, the present work aims at determining a general thermo-elasticity solution to treat both the static and frequency problems of GPLRC rectangular plates under different boundary conditions and embedding foundations, as typically applied in many lightweight mechanical and biomedical components, as well as in membranes and flexible wearable sensors and actuators. Despite the available literature on plate-like nanostructures, usually based on nonclassical approaches, the proposed work explores the capability of a higher-order shear deformation plate formulation combined with a modified Halpin and Tsai model to handle the problem, and checks for the potentials of the generalized differential quadrature (GDQ) approach as high-performance numerical tool to solve the equations even with a reduced computational effort, in lieu of the most common continuum finite element methods from the literature. The governing equations are here derived by means of the Hamilton’s principle, accounting for a modified Halpin–Tsai model for the definition of the material properties and the effect of the dispersion in nanocomposites. The GDQ-based solution is here compared to the analytical once based on a Navier-type expansion, and numerically. An extensive study is performed systematically to analyze the impact of different parameters such as the distribution patterns and weight fractions of the reinforcement phase, complex environments, Winkler–Pasternak foundation coefficients, and Kerr substrate constants on the overall response of FG-GPLRC rectangular plates. Results of the present study would be useful for the design of advanced lightweight composite members in civil and mechanical engineering, due to the importance of nanofillers dispersion and the application of foundation structures. The proposed GDQ method represents an innovative computational tool for design purposes, due to its great capability to solve challenging problems, with high simplicity and accuracy. A further extension of the formulation accounts for the thermal buckling of nanocomposite members within a unified setting, as useful for coupled problems for which theoretical predictions are usually cumbersome to obtain.

2. Theoretical Formulation

Here, we consider a FG-GPLRC rectangular plate resting on an elastic Winkler–Pasternak and Kerr medium, whose geometry and dimensions are depicted in Figure 1. The GPLs reinforcement is assumed to be distributed either uniformly (GPL-UD) or in a functionally graded way throughout the thickness, with two symmetric patterns, GPL-X, and GPL-O, respectively.

2.1. Effective Material Properties

The material properties are here defined according to a modified Halpin–Tsai model, such that the effective Young’s modulus of the GPL/polymer composite $\overline{E}$ reads as follows [34]:

$$\overline{E}=\frac{3}{8}\frac{\left(1+{\xi }_L{\eta }_L{V}_{GPL}\right)}{\left(1-{\eta }_L{V}_{GPL}\right)}{E}_M+\frac{5}{8}\frac{\left(1+{\xi }_W{\eta }_W{V}_{GPL}\right)}{\left(1-{\eta }_W{V}_{GPL}\right)}{E}_M$$

(1)

where

$${\xi }_L=2\frac{L_{GPL}}{t_{GPL}},{\xi }_W=2\frac{W_{GPL}}{t_{GPL}},{\eta }_W=-\frac{1-\left(\frac{E_{GPL}}{E_M}\right)}{{\xi }_W+\left(\frac{E_{GPL}}{E_M}\right)},{\eta }_L=\frac{\left(\frac{E_{GPL}}{E_M}\right)-1}{{\xi }_L+\left(\frac{E_{GPL}}{E_M}\right)}$$

(2)

with $E_M$ and $E_{GPL}$ are the Young’s moduli of the polymer matrix and GPLs, respectively; $V_{GPL}$ is the GPL volume fraction, ${\xi }_L$ and ${\xi }_W$ are the parameters characterizing both the geometry and size of GPL nanofillers; $L_{GPL},W_{GPL}$ and $t_{GPL}$ are the average length, width, and thickness of GPLs, respectively.

In line with findings by Rafiee et al. [35], the effective Young’s modulus of GPL/polymer nanocomposites is well-approximated by the modified Halpin–Tsai model. The result determined by Equation (1), indeed, is only 2.7% higher than the experimental predictions. Based on the same rule of mixtures, the effective Poisson’s ratio and mass density read as follows:

$$\overline{\rho }={\rho }_{GPL}{V}_{GPL}+{\rho }_M\left(1-V_{GPL}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\overline{\nu }={\nu }_{GPL}{V}_{GPL}+{\nu }_M\left(1-V_{GPL}\right)$$

(3)

while, the effective shear modulus is defined as:

$$\overline{G}=\frac{\overline{E}}{2\left(1+\overline{\nu }\right)}$$

(4)

As also depicted in Figure 2, we select three different distribution patterns of GPLs along the thickness direction of the structure, whose analytical expressions take the following form [36]:

$$V_{GPL}\left(z_j\right)=\left\{\begin{array}{l}\left\{\begin{array}{l}2\hspace{0.17em}{V}_{GPL}^{*}\left(1-2\frac{Z}{h}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le Z\le \frac{h}{2}\hfill \\ 2\hspace{0.17em}{V}_{GPL}^{*}\left(1-2\frac{\left(h-Z\right)}{h}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{h}{2}\le Z\le h\hfill \end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}GPL-X\hfill \\ \left\{\begin{array}{l}2\hspace{0.17em}{V}_{GPL}^{*}\left(1-2\frac{\left(\frac{h}{2}-Z\right)}{h}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le Z\le \frac{h}{2}\hfill \\ 2\hspace{0.17em}{V}_{GPL}^{*}\left(1+2\frac{\left(\frac{h}{2}-Z\right)}{h}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{h}{2}\le Z\le h\hfill \end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}GPL-O\hfill \\ \left\{\begin{array}{l}{V}_{GPL}^{*}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le Z\le \frac{h}{2}\hfill \\ V_{GPL}^{*}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{h}{2}\le Z\le h\hfill \end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}UD\right.\hfill \end{array}\right.$$

(5)

being $V_{GPL}^{*}=\frac{{\Lambda }_{GPL}}{\left(\frac{{\rho }_{GPL}}{{\rho }_M}\right)\left(1-{\Lambda }_{GPL}\right)+{\Lambda }_{GPL}}$ and $z_j=\left(\frac{1}{2}+\frac{1}{2n}-\frac{j}{N_L}\right)h,\hspace{1em}j=1,2,3,\dots ,N_L$.

2.2. Displacement Field

As already mentioned in the introduction, we follow a higher order shear deformation theory (HSDT) to define the kinematic field of the structure, i.e., [37].

$$\begin{array}{l}u(x,y,z,t)=u_0(x,y,t)+z\hspace{0.17em}{u}_1(x,y,t)+z^2{u}_2(x,y,t)+z^3{u}_3(x,y,t)\\ v(x,y,z,t)=v_0(x,y,t)+z\hspace{0.17em}{v}_1(x,y,t)+z^2{v}_2(x,y,t)+z^3{v}_3(x,y,t)\\ w(x,y,z,t)=w_0(x,y,t)\end{array}$$

(6)

where $(u,v,w)$ refer to the axial displacement components of an arbitrary point $(x,y,z)$ within the domain; $(u_0,v_0,w_0)$ stand for the related components at the reference mid-plane; $(u_1,v_1,w_1)$ are the rotations of the normal about the y-, x-, and z-axis respectively; $u_2$,$v_2$,$u_3$,$v_3$ define the higher-order terms in the Taylor’s series expansion. Also, the non-null strain components are defined in Appendix A.

The constitutive relations for the elastic problem are expressed as:

$${\left\{\begin{array}{l}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{yz}\\ {\tau }_{xz}\\ {\tau }_{xy}\end{array}\right\}}^{(K)}={\left[\begin{array}{ccccc}{Q}_{11}& Q_{12}& 0& 0& 0\\ Q_{21}& Q_{22}& 0& 0& 0\\ 0& 0& Q_{44}& 0& 0\\ 0& 0& 0& Q_{55}& 0\\ 0& 0& 0& 0& Q_{66}\end{array}\right]}^{(K)}{\left\{\begin{array}{l}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{yz}\\ {\gamma }_{xz}\\ {\gamma }_{xy}\end{array}\right\}}^{(K)}$$

(7)

where the elastic constants are defined in Appendix A.

2.3. Hamilton’s Principle and Governing Equations

The fundamental equations of the problem are determined by applying the Hamilton’s principle, in the following variational energy form [38]:

$${\displaystyle {\int }_{t_1}^{t_2}\left(\delta {\Phi }_k-\delta {\Phi }_e-\left(\delta {\Phi }_{w1}+\delta {\Phi }_{w2}+\delta {\Phi }_{w3}+\delta {\Phi }_{w4}\right)\right)\mathrm{dt}=0}$$

(8)

where ${\Phi }_k$ and ${\Phi }_e$ stand for the kinetic and elastic energy, respectively, and the external work ${\Phi }_w$ is split as ${\Phi }_{w1}$, ${\Phi }_{w2}$, ${\Phi }_{w3}$_, and ${\Phi }_{w4}$ whose definition depends on the elastic Winkler–Pasternak and Kerr substrates, as well as on the mechanical loading, respectively. The above-mentioned quantities are defined in a variational form as:

$$\delta {\Phi }_k={\displaystyle \underset{V}{\int }\rho \left(\frac{\partial U}{\partial t}\frac{\partial \delta U}{\partial t}+\frac{\partial V}{\partial t}\frac{\partial \delta V}{\partial t}+\frac{\partial W}{\partial t}\frac{\partial \delta W}{\partial t}\right)dV}$$

(9)

$$\delta {\Phi }_e={\displaystyle \underset{V}{\int }\left({\sigma }_{xx}\delta {\epsilon }_{xx}+{\sigma }_{yy}\delta {\epsilon }_{yy}+{\sigma }_{zz}\delta {\epsilon }_{zz}+{\tau }_{xy}\delta {\gamma }_{xy}+{\tau }_{yz}\delta {\gamma }_{yz}+{\tau }_{xz}\delta {\gamma }_{xz}\right)dV}$$

(10)

In addition

$$\delta {\Phi }_{w1}=\begin{array}{c}{\displaystyle \underset{A}{\int }\left(-k_w{w}_o+k_p\left(\frac{{\partial }^2{w}_o}{\partial x^2}+\frac{{\partial }^2{w}_o}{\partial y^2}\right)\right)\delta w_{o}dA}\end{array}$$

(11)

$k_w$ and $k_p$ being the Winkler and Pasternak constants.

While

$$\delta {\Phi }_{w2}=\begin{array}{c}{\displaystyle \underset{A}{\int }\left(-\frac{k_l{k}_u}{k_l+k_u}{w}_o+\frac{k_s{k}_u}{k_l+k_u}\left(\frac{{\partial }^2{w}_o}{\partial x^2}+\frac{{\partial }^2{w}_o}{\partial y^2}\right)\right)\delta w_{o}dA}\end{array}$$

(12)

where $k_s$, $k_u$, $k_l$, refer to the shear layer, upper, and lower spring layers, respectively [39]. The last energy contribution related to the external load P acting on the top surface of the plate reads as follows [40]:

$$\delta {\Phi }_{w3}=-{\displaystyle \underset{A}{\int }P\delta w_0}dA$$

(13)

In addition, the conductive layer reinforced with GPLs satisfies the following Fourier heat conduction relation:

$${\nabla }^{2}T+R=\rho c\frac{\partial T}{\partial t}$$

(14)

In absence of a thermal generation, in steady-state conditions, it is:

$${\nabla }^{2}T=0$$

(15)

For a conductive layer reinforced with a UD or FG distribution of GPLs, we get the following relations:

$$\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}UD$$

(16a)

$$\frac{\partial }{\partial x}\left(K_c\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_c\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_c\frac{\partial T}{\partial z}\right)=0\hspace{0.17em}\mathrm{for}\mathrm{FG}$$

(16b)

whose thermal boundary conditions read as follows:

$$T(0,y,z)=0,\hspace{0.17em}T(a,y,z)=0,\hspace{0.17em}\hspace{0.17em}T(x,0,z)=0,\hspace{0.17em}\hspace{0.17em}T(x,b,z)=0,\hspace{0.17em}\hspace{0.17em}T(x,y,0)=T_1,\hspace{0.17em}\hspace{0.17em}T(x,y,h)=T_2$$

(17)

The last energy contribution related to the thermal load can be obtained as:

$$\delta \hspace{0.17em}{\Phi }_{w4}={\displaystyle \underset{A}{\int }\left((N_1^T)\frac{\partial w_0}{\partial x}\frac{\partial \delta w_0}{\partial x}+(N_2^T)\frac{\partial w_0}{\partial y}\frac{\partial \delta w_0}{\partial y}\right)}dA$$

(18)

with

$$N_1^{\mathrm{T}}={\displaystyle \underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}(Q_{11}}+Q_{12}){\alpha }_C(T-T_0)dz,N_2^{\mathrm{T}}={\displaystyle \underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}(Q_{21}}+Q_{22}){\alpha }_C(T-T_0)dz$$

(19)

and $T_0$ being the ambient temperature.

By substitution of Equations (9)–(13), and (18) into Equation (8), after a mathematical manipulation we get the following equations as presented in Appendix A (Equations (A4)–(A12)).

3. Thermal Field

To satisfy the thermal boundary conditions in Equation (17), we introduce a Fourier-type solution as follows:

$$T={\displaystyle \sum _{m=1}^{\infty }{\displaystyle \sum _{n=1}^{\infty }{T}_{mn}(z)\sin (P_{m}x)\hspace{0.17em}\sin (P_{n}y)\hspace{0.17em}\hspace{0.17em}.\hspace{0.17em}}}$$

(20)

where $P_m=m\pi /a$ and $P_n=n\pi /b$. Moreover, the thermal conductivity coefficients related to the GPLs distribution pattern are determined as:

$$UD:\hspace{0.17em}\frac{K_c}{K_m}=1+D$$

(21a)

$$GPL-X\hspace{0.17em}:\hspace{0.17em}\left\{\begin{array}{c}\frac{K_c}{K_m}=1+2D(1-2\frac{Z}{h})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le Z\le \frac{h}{2}\\ \frac{K_c}{K_m}=1+2D(-1+2\frac{Z}{h})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{h}{2}\le Z\le h\end{array}\right.$$

(21b)

$$GPL-O\hspace{0.17em}:\hspace{0.17em}\left\{\begin{array}{c}\frac{K_c}{K_m}=1+4D(\frac{Z}{h})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le Z\le \frac{h}{2}\\ \frac{K_c}{K_m}=1+4D(1-\frac{Z}{h})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{h}{2}\le Z\le h\end{array}\right.$$

(21c)

Inserting Equations (20) and (21a) into Equation (16a) and solving the equation analytically in its final form, we obtain the following expression of temperature gradient for GPLRC rectangular plates with a uniform distribution of GPLs.

$$T_{\mathrm{mn}}={\mathrm{C}}_{11}{\mathrm{e}}^{\sqrt{{\mathrm{A}}_{11}}Z}+C_{22}{\mathrm{e}}^{-\sqrt{{\mathrm{A}}_{11}}Z}$$

(22)

where ${\mathrm{C}}_{11}{,\mathrm{C}}_{22},$ and ${\mathrm{A}}_{11}$ are arbitrary constants determined with appropriate enforcement of the thermal surface boundary conditions (see more details in Appendix B).

At the same time, by combining Equations (20), (21b) and (21c), the heat conduction differential Equation (16b) reduces to the following hypergeometric equation:

$$\left({\mathrm{A}}_1\mathrm{Z}+{\mathrm{A}}_2\right)\frac{{\partial }^2{T}_{\mathrm{m}\mathrm{n}}\left(\mathrm{z}\right)}{\partial {\mathrm{Z}}^2}+\left({\mathrm{A}}_3\mathrm{Z}+{\mathrm{A}}_4\right)\frac{\partial T_{\mathrm{m}\mathrm{n}}\left(\mathrm{z}\right)}{\partial \mathrm{Z}}+\left({\mathrm{A}}_5\mathrm{Z}+{\mathrm{A}}_6\right)T_{\mathrm{m}\mathrm{n}}\left(\mathrm{z}\right)=0$$

(23)

with ${\mathrm{A}}_1{,\mathrm{A}}_2,\dots {,\mathrm{A}}_6$ being constant coefficients depending on the pattern of GPLs distribution (see Appendix B). The analytical solution of Equation (23) takes the following form:

$$T_{\mathrm{m}\mathrm{n}}={\mathrm{C}}_1{e}^{{\alpha }_1}\mathrm{K}\mathrm{u}\mathrm{m}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{M}({\beta }_1,{\beta }_2,{\beta }_3){({\mathrm{A}}_2+{\mathrm{A}}_1\mathrm{Z})}^{{\alpha }_2}+{\mathrm{C}}_2{\mathrm{e}}^{{\alpha }_1}\mathrm{K}\mathrm{u}\mathrm{m}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{U}({\beta }_1,{\beta }_2,{\beta }_3){({\mathrm{A}}_2+{\mathrm{A}}_1\mathrm{Z})}^{{\alpha }_2}$$

(24)

where the Kummer’s function, also known as the confluent hypergeometric function of the first kind, is a solution to a Kummer’s differential equation. In addition, KummerU and KummerM functions represent two special types of Kummer function.

As the thermal behavior of a structure depends on its thermo-mechanical properties, it is worth noticing in Figure 3 that the temperature distribution in thermoelastic solutions is completely different from a uniform or harmonic distribution.

4. Solution Procedure

4.1. Analytical Solution

Following a Navier-type procedure, we now introduce the analytical solution to the above-mentioned governing equations for simply supported FG-GNPRC plates, namely [37]:

$$\left(u_0,\hspace{0.17em}{u}_1,\hspace{0.17em}{u}_2,\hspace{0.17em}{u}_3\right)={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{\infty }\left(U_0,\hspace{0.17em}{U}_1,\hspace{0.17em}{U}_2,\hspace{0.17em}{U}_3\right)}}\cos \left(p_{m}x\right)\sin \left(p_{n}y\right)\exp \left(i\omega t\right)$$

(25a)

$$\left(v_0,\hspace{0.17em}{v}_1,\hspace{0.17em}{v}_2,\hspace{0.17em}{v}_3\right)={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{\infty }\left(V_0,\hspace{0.17em}{V}_1,\hspace{0.17em}{V}_2,\hspace{0.17em}{V}_3\right)}}\sin \left(p_{m}x\right)\cos \left(p_{n}y\right)\exp \left(i\omega t\right)$$

(25b)

$$w_0={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{\infty }{W}_0}}\sin \left(p_{m}x\right)\sin \left(p_{n}y\right)\exp \left(i\omega t\right)$$

(25c)

Substituting Equations (25a)–(25c) into Equations (A4)–(A12) (see Appendix A), it is possible to derive the following relations in matrix form, under the assumption $p=0$

$$\left(\left[K\right]-\left[M\right]{\omega }^2\right)\left\{\delta \right\}=\left\{0\right\}$$

(26)

with $\left[K\right]$ and $\left[M\right]$ being the stiffness and mass matrix, respectively, and ${\left\{\delta \right\}}^T=\left\{U_0,V_0,W_0,U_1,V_1,U_2,V_2,U_3,V_3\right\}$. Thus, the natural frequencies are determined by means of the following eigenvalue relation:

$$\left|\left[K\right]-\left[M\right]{\omega }^2\right|=0$$

(27)

4.2. Numerical Solution

The same problem is also solved numerically by means of the GDQ method, due to its capability to yield accurate solutions even with a reduced computational effort [41,42,43,44,45,46] while maintaining a certain flexibility when involving any kind of boundary condition along the structural edges. The proposed method allows solving of the problem in a strong form, by discretizing the derivatives of a function in the following form [41]:

$$\hspace{0.17em}{\left.\frac{\partial f}{\partial x}\right|}_{x=x_i,\hspace{0.17em}y=y_j}={\displaystyle \sum _{m=1}^{N_x}{\displaystyle \sum _{n=1}^{N_y}{A}_{im}^x}}{I}_{jn}^y{f}_{mn}$$

(28a)

$$\hspace{0.17em}{\left.\frac{\partial f}{\partial y}\right|}_{x=x_i,\hspace{0.17em}y=y_j}={\displaystyle \sum _{m=1}^{N_x}{\displaystyle \sum _{n=1}^{N_y}{I}_{im}^x}}{A}_{jn}^y{f}_{mn}$$

(28b)

$$\frac{\partial }{\partial r}\hspace{0.17em}\left({\left.\frac{\partial f}{\partial \theta }\right|}_{x=x_i,\hspace{0.17em}y=y_j}\right)={\displaystyle \sum _{m=1}^{N_x}{\displaystyle \sum _{n=1}^{N_y}{A}_{im}^x}}{A}_{jn}^y{f}_{mn}$$

(28c)

$$\hspace{0.17em}{\left.\frac{{\partial }^{2}f}{\partial x^2}\right|}_{x=x_i,\hspace{0.17em}y=y_j}={\displaystyle \sum _{m=1}^{N_x}{\displaystyle \sum _{n=1}^{N_y}{B}_{im}^x}}{I}_{jn}^y{f}_{mn}$$

(28d)

$$\hspace{0.17em}{\left.\frac{{\partial }^{2}f}{\partial y^2}\right|}_{x=x_i,\hspace{0.17em}y=y_j}={\displaystyle \sum _{m=1}^{N_x}{\displaystyle \sum _{n=1}^{N_y}{I}_{im}^x}}{B}_{jn}^y{f}_{mn}$$

(28e)

where $I_{im}^x$ and $I_{jn}^y$ are equal to one when $i=m$ and $j=n$, or equal to zero, otherwise. In addition, $A_{im}^x$, $A_{jn}^y$, $B_{im}^x$ and $B_{jn}^y$ are the weighting coefficients of the first and second-order derivatives along the $x$ and $y$-directions, respectively, defined as:

$$A_{im}^{\left(1\right)}=\left\{\begin{array}{l}\frac{\xi \left(x_i\right)}{\left(x_i-x_m\right)\xi \left(x_m\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}when\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i\ne m\hfill \\ -{\displaystyle \sum _{k=1,k\ne i}^{N_x}{A}_{ik}^{\left(1\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}when\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=m}\hfill \end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,m=1,2,\dots ,N_x$$

(29a)

$$A_{jn}^{\left(1\right)}\hspace{0.17em}=\left\{\begin{array}{l}\frac{\xi \left(y_j\right)}{\left(y_j-y_n\right)\xi \left(y_n\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}when\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j\ne n\hfill \\ -{\displaystyle \sum _{k=1,k\ne j}^{N_y}{A}_{jk}^{\left(1\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}when\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=n}\hfill \end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j,n=1,2,\dots ,N_y$$

(29b)

with

$$\xi \left(x_i\right)={\displaystyle \prod _{k=1,k\ne i}^{N_x}\left(x_i-x_k\right)}$$

(30a)

$$\xi \left(y_j\right)={\displaystyle \prod _{k=1,k\ne j}^{N_y}\left(y_j-y_k\right)}$$

(30b)

and

$$B_{im}^{\left(2\right)}=2\left(A_{ii}^{\left(1\right)}{A}_{im}^{\left(1\right)}-\frac{A_{im}^{\left(1\right)}}{\left(x_i-x_m\right)}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,m=1,2,\dots ,N_x\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}i\ne m$$

(31a)

$$B_{jn}^{\left(2\right)}=2\left(A_{jj}^{\left(1\right)}{A}_{jn}^{\left(1\right)}-\frac{A_{jn}^{\left(1\right)}}{\left(y_j-y_n\right)}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j,n=1,2,\dots ,N_y\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}j\ne n$$

(31b)

$$B_{ii}^{\left(2\right)}=-{\displaystyle \sum _{k=1,k\ne i}^{N_x}{B}_{ik}^{\left(2\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}i=1,2,\dots ,N_x},\hspace{0.17em}\hspace{0.17em}i=m$$

(31c)

$$B_{jj}^{\left(2\right)}=-{\displaystyle \sum _{k=1,k\ne j}^{N_y}{B}_{jk}^{\left(2\right)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}j=1,2,\dots ,N_y},\hspace{0.17em}\hspace{0.17em}j=n$$

(31d)

In what follows, we select a Chebyshev distribution of grid points within the domain defined as:

$$x_i=\frac{a}{2}\left(1-\cos \left(\frac{\left(i-1\right)}{\left(N_x-1\right)}\pi \right)\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2,3,\hspace{0.17em}\dots \hspace{0.17em},N_x$$

(32a)

$$y_j=\frac{b}{2}\left(1-\cos \left(\frac{\left(j-1\right)}{\left(N_y-1\right)}\pi \right)\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,3,\hspace{0.17em}\dots \hspace{0.17em},N_y$$

(32b)

Thus, the algebraic eigenvalue problem can be redefined in matrix form as:

$$\left\{\left[\begin{array}{cc}\left[M_{dd}\right]& \left[M_{db}\right]\\ \left[M_{bd}\right]& \left[M_{bb}\right]\end{array}\right]{\omega }_{mn}^2+\left[\begin{array}{cc}\left[K_{dd}\right]& \left[K_{db}\right]\\ \left[K_{bd}\right]& \left[K_{bb}\right]\end{array}\right]\right\}\left\{\begin{array}{l}{\delta }_d\\ {\delta }_b\end{array}\right\}=0$$

(33)

where we distinguish among inner and boundary grid-points, by means of subscripts $d$ and $b,$ respectively, whereas $\delta$ stands for the displacement vector. The natural frequencies of the problem are derived as solutions of Equation (33).

5. Minimum Total Potential Energy Principle

Now we apply the minimization procedure of the total energy associated to the structural system to study its static response [47], namely:

$$\delta \left({\Phi }_e+{\Phi }_{w1}+{\Phi }_{w2}+{\Phi }_{w3}+{\Phi }_{w4}\right)=0$$

(34)

which is combined with the energy quantities in Equations (10)–(13) and (18) to yield the following governing equations of GPL reinforced composite rectangular plates (see Equations (A34)–(A42) in Appendix B).

Bending Analysis

The analytical solution for a static problem stems from a Fourier-type series discretization of the mechanical force, as follows [48]:

$$P_{mn}={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{\infty }{q}_{mn}}}\sin \left(p_{m}x\right)\cos \left(p_{n}y\right)$$

(35)

in which $q_{mn}=\frac{4p_0}{mn{\pi }^2}\left(1-{\left(-1\right)}^n\right)\left(1-{\left(-1\right)}^m\right),$ and $p_0=0.1\hspace{0.17em}\mathrm{MPa}.$

By substitution of Equations (35) and (25a)–(25c) into Equations (A34)–(A42), with the assumption $\omega =0$, we get the following relation:

$$\left\{\left[F\right]+\left[K\right]\right\}\left\{\delta \right\}=0$$

(36)

which is solved in terms of the kinematic unknowns.

At the same time, based on a GDQ definition of the problem, the substitution of Equations (28a)–(28e) into Equations (A34)–(A42) leads to the following relation in matrix form:

$$\left\{\left[\begin{array}{cc}\left[F_{dd}\right]& \left[F_{db}\right]\\ \left[F_{bd}\right]& \left[F_{bb}\right]\end{array}\right]+\left[\begin{array}{cc}\left[K_{dd}\right]& \left[K_{db}\right]\\ \left[K_{bd}\right]& \left[K_{bb}\right]\end{array}\right]\right\}\left\{\begin{array}{l}{\delta }_d\\ {\delta }_b\end{array}\right\}=0$$

(37)

depending on the kinematic unknowns ${\delta }_d$ and ${\delta }_b$.

6. Results and Discussion

6.1. Validation

We now present the results from a large numerical investigation aimed at studying the static and vibrational response of GPLRC multilayer rectangular plates, with material properties as summarized in

Table 1. Material properties of the system (see Ref. [3]).

Polymer Epoxy (Matrix)	Graphene Platelets
${\nu }_m=0.34$	${\nu }_{GPL}=0.186$
${\rho }_m\hspace{0.17em}[\mathrm{k}\mathrm{g}/{\mathrm{m}}^3]=1.2\times {10}^3$	${\rho }_{GPL}\hspace{0.17em}\hspace{0.17em}[\mathrm{k}\mathrm{g}/{\mathrm{m}}^3]=1.06\times {10}^3$
$E_m\hspace{0.17em}[\mathrm{G}\mathrm{P}\mathrm{a}]=2.85\hspace{0.17em}$	$E_{GPL}\hspace{0.17em}[\mathrm{T}\mathrm{P}\mathrm{a}]=1.01$
	$l_{GPL}\hspace{0.17em}[\mathsf{\mu }\mathrm{m}]=2.5$
	$w_{GPL}\hspace{0.17em}[\mathsf{\mu }\mathrm{m}]=1.5$
	$t_{GPL}\hspace{0.17em}[\mathrm{n}\mathrm{m}]=1.5$

. After a preliminary convergence study, we test the performances of our proposed formulation with a comparative evaluation against the open literature or further numerical methods.

The GDQ numerical study starts considering the effect of an increased grid point distribution on the structural response in terms of dimensionless fundamental frequency and bending deflection, for two different boundary conditions (completely clamped and simply-supported), as visible in Figure 4 and Figure 5, respectively. Based on the plots in these figures it is worth noticing the very fast stabilization of results even with a reduced number of sampling points, whose rate of convergence maintains almost constant independently of the selected boundary conditions. As a further step we perform a parametric evaluation of the frequency response for GPL reinforced thick rectangular plates in terms of natural frequencies for different reinforcement distributions, accounting for different longitudinal and transverse modes, as summarized in

Table 2. Dimensionless natural frequency of a rectangular plate made by pure epoxy and different types of graphene distribution with various longitudinal and transverse mode shapes for ($a\times b\times h=$ 0.45 m × 0.45 m × 0.045 m),${\overline{\omega }}_{mn}={\omega }_{mn}h\sqrt{{\rho }_m/E_m},$ and ${\Lambda }_{GPL}=1\hspace{0.17em}\%$.

GPL-X	GPL-O	GPL-UD	Pure Epoxy	${\overline{\omega }}_{\mathit{m}\mathit{n}}$	(m, n)
0.1378	0.102	0.1216	0.0584	Ref. [16]	(1, 1)
0.1378	0.102	0.1216	0.0584	FSDT
0%	0%	0%	0%	Error%
0.3249	0.2456	0.2895	0.1391	Ref. [16]	(2, 1)
0.3249	0.2456	0.2895	0.1391	FSDT
0%	0%	0%	0%	Error%
0.4939	0.3796	0.4436	0.2132	Ref. [16]	(2, 2)
0.4939	0.3796	0.4436	0.2132	FSDT
0%	0%	0%	0%	Error%
0.5984	0.4645	0.54	0.2595	Ref. [16]	(3, 1)
0.5984	0.4645	0.54	0.2595	FSDT
0%	0%	0%	0%	Error%
0.7454	0.586	0.6767	0.3251	Ref. [16]	(3, 2)
0.7454	0.586	0.6767	0.3251	FSDT
0%	0%	0%	0%	Error%
0.969	0.7755	0.8869	0.4261	Ref. [16]	(3, 3)
0.969	0.7755	0.8869	0.4261	FSDT
0%	0%	0%	0%	Error%

. A FSDT is adopted in this case for comparative purposes with predictions by Song et al. [16], with a perfect matching for all the selected graphene distributions and mode shapes. Among the different GPL distributions, it seems that a GPL-X distribution predicts the highest vibrational frequencies of the system, whereas a pure epoxy material yields the lowest vibrational values.

We focus, now, on the statics of FG square plates, with the upper and lower surfaces made by a pure ceramic and metal, respectively. The mechanical properties of the system are assumed to vary along the thickness direction based on the following relation [49]:

$$\begin{gathered} E(z)=\left(E_c-E_m\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^p+E_m \\ \rho (z)=\left({\rho }_c-{\rho }_m\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^p+{\rho }_m \\ v(z)=\left(v_c-v_m\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^p+v_m \end{gathered}$$

(38)

where subscripts $c$ and $m$ refer to the ceramic and metal phases, respectively, and $p$ is the power-law index of the FG material. The mechanical properties of pure ceramic and metal phases are summarized in

Table 3. Material properties of the system.

Property	Value
Ec	380 GPa
Em	70 GPa
vc	0.3
vm	0.3

, where the ceramic phase features a meaningful higher stiffness than a pure metal. This means that an increased quantity of metal components in the structure would increase its flexibility, as specified in the deflection response of

Table 4. Effect of the volume fraction exponent of the FG material on the dimensionless deflections of SSSS FG square plates.

Present	Ref. [49]	Ref. [48]	p
0.4621	0.4666	0.4665	p = 0 (Ceramic)
0.9416	-	0.9287	1
1.2002	1.1908	1.194	2
1.3204	-	1.32	3
1.3869	1.3769	1.389	4
1.4342	-	1.4356	5
1.4741	1.4554	1.4727	6
1.5107	-	1.5049	7
1.5455	1.5157	1.5343	8
1.579	-	1.5617	9
1.6112	1.5695	1.5876	10
2.5085	-	2.5327	p = ∞ (Metal)

, for different exponents $p$, in line with findings by References [48,49]. As further validation, we compare our vibration results with finite element predictions as performed in COMSOL, for a square plate with constant thickness $h=10 \mathrm{mm}$, length 250 mm, and different GPL distribution profiles. As the GPLRC plate features an inhomogeneous behavior, the material properties are assumed to be graded in the thickness direction. In

Table 5. Comparison with FEM (GPL-UD).

S. No	Mode	$\omega \left(\mathrm{Hz}\right)$	$\mathrm{FEM} \left(\mathrm{Hz}\right)$	Relative Error%
S1	1	494.9	495.9	0.19%
S2	2	996.8	997.3	0.05%
S3	3	1453.0	1451.9	0.07%
S4	4	1766.1	1763.9	0.12%
S5	5	2180.2	2175.1	0.23%
S6	6	2752.5	2743.3	0.33%

,

Table 6. Comparison between FEM-based predictions and results from our formulation (GPL-X).

S. No	Mode	$\omega \left(\mathrm{Hz}\right)$	$\mathrm{FEM} \left(\mathrm{Hz}\right)$	Relative Error
S1	1	520.2	520.51	0.06%
S2	2	1046.4	1044.0	0.23%
S3	3	1523.8	1516.3	0.49%
S4	4	1840.6	1840.1	0.027%
S5	5	2283.2	2264.3	0.83%
S6	6	2879.9	2850.2	1.04%

and

Table 7. Comparison between FEM-based predictions and results from our formulation (GPL-O).

S. No	Mode	$\omega \left(\mathrm{Hz}\right)$	$\mathrm{FEM} \left(\mathrm{Hz}\right)$	Relative Error
S1	1	468.1	470.1	0.42%
S2	2	944.0	947.3	0.34%
S3	3	1377.3	1381.4	0.29%
S4	4	1665.5	1670.2	0.28%
S5	5	2069.4	2074.0	0.22%
S6	6	2615.0	2620.4	0.20%

, we compare the vibration results for the first six modes, and for three different reinforcement distributions, namely, a GPL–UD, –X, and –O profile, respectively, It is worth noticing the very good matching of results compared to finite elements, with a limited percentage error (lower than 1.04% in each case), despite the reduced computational effort. For a GPL-O distribution we also represent the related mode shapes in Figure 6, while keeping ${\Lambda }_{GPL}=0.1\hspace{0.17em}$(wt%), with a reasonable kinematic response with the selected simply-supported boundary condition.

6.2. Parametric Study

We perform, now, a parametric study focusing on the sensitivity of the frequency response $\Delta \omega \%=\left({\omega }_c-{\omega }_{epoxy}\right)/{\omega }_{epoxy}\times 100$ against the weight fraction of GPLs for each selected distribution pattern (Figure 7), with a clear beneficial effect on the structural stiffness and stability as ${\Lambda }_{GPL}$ increases within the material. Moreover, due to the presence of high normal stresses in the upper and lower sides of multilayered plates, a GPL-X distribution with an increased quantity of GPLs causes a hardening effect on the system, together with higher natural frequencies. On the other hand, a higher concentration of the reinforcing phase in the middle surface of the plate increases the structural deformability monotonically, whose variation rate is plotted in Figure 7, with a perfect agreement with predictions by Song et al. [16]. A further parameter affecting the frequency response is the number of layers $N_L$ in the structure, as plotted in Figure 8 for three different reinforcement distributions. Except for the uniform reinforcement case for which the response is unaffected by $N_L$, a small increase of this parameter can cause some hardening or softening effects on the structural stiffness, with a sharp increase or decrease of the frequency, for a GPL-X or GPL-O distribution, respectively. Even for these two distributions, the solutions stabilize for a number of layers equal or higher than 10, as also predicted by Song et al. [16]. The sensitivity of the response to the number of layers in the thickness direction is plotted in Figure 9 in terms of dimensionless vibrational frequency, for a square plate with $a/h=25$ and ${\Lambda }_{GPL}=0.3\%$ under a completely-clamped (CCCC) and simply-supported (SSSS) boundary condition. Based on a comparative evaluation of the plots in these two figures, the best stability response seems to be reached for a CCCC multilayered structure with $N_L=10$ and a GPL-X reinforcement distribution, in terms of vibration frequency (see Figure 9). At the same time, an increased number of layers more than 10 becomes deleterious for the overall stability of plates with a GPL-O type distribution, both for CCCC and SSSS boundary conditions. Figure 10 depicts the variation of the first vibration frequency against the thermal gradient $\Delta T$ of moderately thick square plates for different GPL weight fractions. Due to the variation of the thermoelastic properties of the system, together with the presence of initial internal stresses and strains in the structure by thermal attacks, this complex environment causes a combined hardening-softening impact on the system. More specifically, an increased thermal variation decreases the fundamental frequency of the system up to a certain value of ΔT, for which the fundamental frequency becomes zero, and the structure undergoes a static instability phenomenon. This critical temperature moves towards higher values for increased weight fractions of GPLs (from 0.2% up to 0.8%). It is worth also noticing that in the pre-divergence zone, by approaching the static instability phenomenon, an enhanced temperature causes a very fast reduction of the vibrational frequency of the system. In order to survey the effect of the Kerr foundation on the vibrational behavior and static instability of GPL-reinforced plates, we plot the variation of dimensionless first fundamental frequency versus the thermal variation for different substrate coefficients while selecting an SSSS and CCCC boundary condition, respectively, in Figure 11a,b. As the presence of an elastic foundation can vary the bending stiffness of a structure, we note that an enhanced value of the foundation constants improves the vibrational behavior of the system. Meanwhile, the divergence instability can be delayed by increasing the substrate coefficients values. In other words, the presence of an elastic medium gets higher critical values at which the static instability phenomenon takes place. In addition, clamped boundary constraints reduce the positive effect of foundations, with a less pronounced variation in the critical temperature corresponding to the static instability and natural vibrational frequencies of the system.

A further goal of the systematic analysis is also the evaluation of the maximum deflection of the structure, hereafter reported in dimensionless form. In Figure 12 we show the variation of this kinematic quantity versus the number of layers within a SSSS (Figure 12a) and CCCC (Figure 12b) laminated structure, accounting for the three different GPLs patterns. Unlike the UD of GPLs, the kinematic response seems to be clearly sensitive to the number of layers $N_L$ within a multilayered structure, with a monotone increase or decrease depending on whether a GPL–O or –X distribution is selected as reinforcing phase, with a plateau obtained in correspondence of $N_L=10$. With regard to the sensitivity of the deflection response to the reinforcement weight fraction, we plot the results in Figure 13a,b for a SSSS and CCCC boundary condition, respectively. Based on the plots in this figure, an increased amount of GPL nanofillers decreases monotonically the overall deformability of the structure. For example, the introduction of a small percentage of GPLs (equal to 0.5%) is able to reduce the deformability of the system up to a percentage of 360%, for a GPL-X symmetric distribution. This last reinforcement dispersion provides the highest stiffness in multilayered structures, for both a SSSS and CCCC boundary conditions, whereas the highest deformability is obtained for GPL-O distributions of the reinforcing phase under the same weight fraction assumptions. From a design standpoint, a GPL-X symmetric dispersion is desirable as the best reinforcing distribution among others, due to its capability to limit the structural deformability. At the same time, a further reduction in deformability can be obtained, accounting for the elastic properties of the surrounding medium, as plotted in Figure 14 and Figure 15 for a Winkler–Pasternak or Kerr elastic substrate, respectively. In both cases, indeed, elastic foundations with increased stiffness properties get lower deflections, while keeping fixed the GPL weight fraction within the structure. This reduction is even more pronounced for more relaxing boundary conditions as simply-supports, while assuming a Kerr medium in lieu of a Winkler–Pasternak-type foundation (compare the plots of Figure 14 and Figure 15).

7. Conclusions

In this paper we focused on the vibrational and static response of FG-GPLRC multilayer rectangular plates by combining a higher order formulation of shells together with a modified Halpin and Tsai model to account the effect of the dispersion in nanocomposites. The problem is here solved both theoretically with a Navier-type solution, and computationally, by means of the GDQ approach, as high-performance numerical tool. The proposed model is successfully validated in its accuracy against predictions from the literature and results from finite element formulations in the first part. Based on a parametric study, it seems that A GPL-X pattern in a multilayered member provides the highest fundamental frequencies and stiffness of the structure. These mechanical properties increase for an increased GPL weight fraction within the material. In addition, the vibration and kinematic results based on a uniform distribution of GPLs are always unaffected by the reinforcement weight fraction and number of layers within the structure. At the same time, the elastic foundations with increased stiffness properties reduce the overall deformability of multilayered GPL-reinforced structures, which confirms the importance of considering the correct mechanical performances of different substrates around a structural member for design purposes. It is also observed that the presence of a thermal environment reduces the structural efficiency and stiffness due to the introduction of an additional stress and strain field in the system. Meanwhile, elastic foundations with increased stiffness properties raise the critical temperature of multilayered structures while reducing their deformability. This study would provide useful scientific insights and an enhanced tool to engineers and designers for the development of novel and efficient composite structures and components, such as electronic circuits, sensors, or flexible electrodes for displays and solar cells.