*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{cases} u(x,y,z,t) = u\_0(x,y,t) + z \, 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 \, 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{cases} \tag{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 midplane; (*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*<sup>3</sup> 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:

$$\begin{Bmatrix} \sigma\_{\mathbf{x}} \\ \sigma\_{y} \\ \tau\_{yz} \\ \tau\_{xz} \\ \tau\_{xy} \end{Bmatrix}^{(K)} = \begin{bmatrix} 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{bmatrix}^{(K)} \begin{Bmatrix} \varepsilon\_{\mathbf{x}} \\ \varepsilon\_{y} \\ \gamma\_{yz} \\ \gamma\_{xz} \\ \gamma\_{xy} \end{Bmatrix}^{(K)} \tag{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]:

$$\int\_{t\_1}^{t\_2} \left( \delta \Phi\_k - \delta \Phi\_{\varepsilon} - \left( \delta \Phi\_{w1} + \delta \Phi\_{w2} + \delta \Phi\_{w3} + \delta \Phi\_{w4} \right) \right) \text{d}t = 0\tag{8}$$

where Φ*<sup>k</sup>* and Φ*<sup>e</sup>* stand for the kinetic and elastic energy, respectively, and the external work Φ*<sup>w</sup>* is split as Φ*w*1**,** Φ*w*2, Φ*w*3*,* and Φ*w*<sup>4</sup> 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 = \int\_V \rho \left( \frac{\partial \mathcal{U}}{\partial t} \frac{\partial \delta \mathcal{U}}{\partial t} + \frac{\partial \mathcal{V}}{\partial t} \frac{\partial \delta \mathcal{V}}{\partial t} + \frac{\partial \mathcal{W}}{\partial t} \frac{\partial \delta \mathcal{W}}{\partial t} \right) dV \tag{9}$$

$$\delta\Phi\_{\varepsilon} = \int\_{V} \left( \sigma\_{xx}\delta\varepsilon\_{xx} + \sigma\_{yy}\delta\varepsilon\_{yy} + \sigma\_{zz}\delta\varepsilon\_{zz} + \tau\_{xy}\delta\gamma\_{xy} + \tau\_{yz}\delta\gamma\_{yz} + \tau\_{xz}\delta\gamma\_{xz} \right) dV \tag{10}$$

In addition

$$
\delta\Phi\_{w1} = \int\_{A} \left( -k\_w w\_o + k\_p \left( \frac{\partial^2 w\_a}{\partial x^2} + \frac{\partial^2 w\_a}{\partial y^2} \right) \right) \delta w\_o dA \tag{11}
$$

*kw* and *kp* being the Winkler and Pasternak constants.

While

$$
\delta\Phi\_{\rm w2} = \int\_{A} \left( -\frac{k\_l k\_u}{k\_l + k\_u} w\_0 + \frac{k\_s k\_\mu}{k\_l + k\_u} \left( \frac{\partial^2 w\_0}{\partial x^2} + \frac{\partial^2 w\_0}{\partial y^2} \right) \right) \delta w\_0 dA \tag{12}
$$

where *ks*, *ku*, *kl*, 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\_{\text{uv3}} = -\int\_A P \delta w\_0 dA \tag{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} \tag{14}
$$

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

$$
\nabla^2 T = 0 \tag{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 \qquad\qquad\text{for}\quad \text{lUD}\tag{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 \quad \text{for} \quad \text{FG} \tag{16b}$$

whose thermal boundary conditions read as follows:

$$T(0, y, z) = 0, \ T(a, y, z) = 0, \ T(\mathbf{x}, 0, z) = 0, \ T(\mathbf{x}, b, z) = 0, \ T(\mathbf{x}, y, 0) = T\_1, \ T(\mathbf{x}, y, h) = T\_2 \tag{17}$$

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

$$\delta \Phi\_{\text{w4}} = \int\_{A} \left( (N\_1^T) \frac{\partial w\_0}{\partial \mathbf{x}} \frac{\partial \delta w\_0}{\partial \mathbf{x}} + (N\_2^T) \frac{\partial w\_0}{\partial y} \frac{\partial \delta w\_0}{\partial y} \right) dA \tag{18}$$

with

$$\mathbf{N}\_1^T = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} (Q\_{11} + Q\_{12}) \mathbf{a}\_\mathbb{C} (T - T\_0) dz,\\ \mathbf{N}\_2^T = \int\_{-\frac{\hbar}{2}}^{\frac{\hbar}{2}} (Q\_{21} + Q\_{22}) \mathbf{a}\_\mathbb{C} (T - T\_0) \, dz \tag{19}$$

and *T*<sup>0</sup> 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 Fouriertype solution as follows:

$$T = \sum\_{m=1}^{\infty} \sum\_{n=1}^{\infty} T\_{mn}(z) \sin(P\_m x) \sin(P\_n y) \ . \tag{20}$$

where *Pm* = *mπ*/*a* and *Pn* = *nπ*/*b*. Moreover, the thermal conductivity coefficients related to the GPLs distribution pattern are determined as:

$$
\hbar lD : \frac{\mathbb{K}\_{\mathfrak{c}}}{\mathbb{K}\_{m}} = 1 + D \tag{21a}
$$

$$GPL - X : \begin{cases} \frac{K\_c}{K\_m} = 1 + 2D(1 - 2\frac{Z}{h}) & 0 \le Z \le \frac{h}{2} \\\frac{K\_c}{K\_m} = 1 + 2D(-1 + 2\frac{Z}{h}) & \frac{h}{2} \le Z \le h \end{cases} \tag{21b}$$

$$GPL-O: \begin{cases} \frac{K\_c}{K\_m} = 1 + 4D(\frac{Z}{h}) & 0 \le Z \le \frac{h}{2} \\\frac{K\_c}{K\_m} = 1 + 4D(1 - \frac{Z}{h}) & \frac{h}{2} \le Z \le h \end{cases} \tag{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\_{\rm mn} = \mathbb{C}\_{11}\mathbf{e}^{\sqrt{\mathcal{N}\_{11}}\mathbf{Z}} + \mathbb{C}\_{22}\mathbf{e}^{-\sqrt{\mathcal{N}\_{11}}\mathbf{Z}} \tag{22}$$

where C11 , C22, and A11 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:

$$(\mathbf{A}\_1 \mathbf{Z} + \mathbf{A}\_2) \frac{\partial^2 T\_{\rm mn}(\mathbf{z})}{\partial \mathbf{Z}^2} + (\mathbf{A}\_3 \mathbf{Z} + \mathbf{A}\_4) \frac{\partial T\_{\rm mn}(\mathbf{z})}{\partial \mathbf{Z}} + (\mathbf{A}\_5 \mathbf{Z} + \mathbf{A}\_6) T\_{\rm mn}(\mathbf{z}) = 0 \tag{23}$$

with A1, A2, ...,A6 being constant coefficients depending on the pattern of GPLs distribution (see Appendix B). The analytical solution of Equation (23) takes the following form:

$$T\_{\rm mm} = \mathbb{C}\_1 \mathbf{e}^{\mathbb{A}\_1} \text{KummerM}(\beta\_1, \beta\_2, \beta\_3)(\mathbf{A}\_2 + \mathbf{A}\_1 \mathbf{Z})^{\mathbb{A}\_2} + \mathbb{C}\_2 \mathbf{e}^{\mathbb{A}\_1} \text{KummerU}(\beta\_1, \beta\_2, \beta\_3)(\mathbf{A}\_2 + \mathbf{A}\_1 \mathbf{Z})^{\mathbb{A}\_2} \tag{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.

**Figure 3.** Different temperature distributions, when b = 10 h, a = b, Λ*GPL* = 0.3 (wt%), GPL-X, T1 = 300, T2 = 400. (**a**) Representation of the effect of different distribution, (**b**) temperature variation for different temperature boundary T2.

#### **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(\left(\mu\_{0},\mu\_{1},\mu\_{2},\mu\_{3}\right)\right) = \sum\_{n=1}^{\infty} \sum\_{m=1}^{\infty} \left(\mathcal{U}\_{0},\mathcal{U}\_{1},\mathcal{U}\_{2},\mathcal{U}\_{3}\right) \cos(p\_{m}\mathbf{x})\sin(p\_{n}y)\exp(i\omega t)\tag{25a}$$

$$\left(\left(v\_{0\prime}, v\_{1\prime}, v\_{2\prime}, v\_{3}\right) = \sum\_{n=1}^{\infty} \sum\_{m=1}^{\infty} \left(V\_{0\prime}, V\_{1\prime}, V\_{2\prime}, V\_{3}\right) \sin(p\_m \mathbf{x}) \cos(p\_n y) \exp(i\omega t) \tag{25b}$$

$$w\_0 = \sum\_{n=1}^{\infty} \sum\_{m=1}^{\infty} \mathcal{W}\_0 \sin(p\_m x) \sin(p\_n y) \exp(i\omega t) \tag{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) \{ \delta \} = \{ 0 \} \tag{26}$$

with [*K*] and [*M*] being the stiffness and mass matrix, respectively, and {*δ*}*<sup>T</sup>* <sup>=</sup> {*U*0, *<sup>V</sup>*0, *<sup>W</sup>*0, *U*1, *V*1, *U*2, *V*2, *U*3, *V*3}. Thus, the natural frequencies are determined by means of the following eigenvalue relation:

$$\left| [K] - [M]\omega^2 \right| = 0 \tag{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–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]:

$$\left. \frac{\partial f}{\partial x} \right|\_{x=x\_i, y=y\_j} = \sum\_{m=1}^{N\_x} \sum\_{n=1}^{N\_y} A\_{im}^x I\_{ju}^y f\_{mn} \tag{28a}$$

$$\left. \frac{\partial f}{\partial y} \right|\_{x=x\_i, y=y\_j} = \sum\_{m=1}^{N\_x} \sum\_{n=1}^{N\_y} I\_{im}^x A\_{ju}^y f\_{mn} \tag{28b}$$

$$\frac{\partial}{\partial r} \left( \left. \frac{\partial f}{\partial \theta} \right|\_{x=x\_i, y=y\_j} \right) = \sum\_{m=1}^{N\_x} \sum\_{n=1}^{N\_y} A\_{im}^x A\_{jn}^y f\_{mn} \tag{28c}$$

$$\left. \frac{\partial^2 f}{\partial x^2} \right|\_{x=x\_i, y=y\_j} = \sum\_{m=1}^{N\_x} \sum\_{n=1}^{N\_y} B\_{im}^x I\_{jn}^y f\_{mn} \tag{28d}$$

$$\left. \frac{\partial^2 f}{\partial y^2} \right|\_{x=x\_i, y=y\_j} = \sum\_{m=1}^{N\_x} \sum\_{n=1}^{N\_y} I\_{im}^x B\_{jn}^y f\_{mn} \tag{28e}$$

where *I<sup>x</sup> im* and *I y jn* are equal to one when *i* = *m* and *j* = *n*, or equal to zero, otherwise. In addition, *A<sup>x</sup> im*, *<sup>A</sup><sup>y</sup> jn*, *<sup>B</sup><sup>x</sup> im* and *<sup>B</sup><sup>y</sup> jn* are the weighting coefficients of the first and second-order derivatives along the *x* and *y*-directions, respectively, defined as:

$$A\_{im}^{(1)} = \begin{cases} \frac{\tilde{\xi}(x\_i)}{(x\_i - x\_m)\tilde{\xi}(x\_m)} & \text{when} \quad i \neq m\\ -\sum\_{k=1, k \neq i}^{N\_\Gamma} A\_{ik}^{(1)} & \text{when} \quad i = m \end{cases} \quad i, m = 1, 2, \dots, N\_\Gamma \tag{29a}$$

$$A\_{jn}^{(1)} = \begin{cases} \frac{\xi(y\_j)}{\left(y\_j - y\_n\right)\overline{\xi}(y\_n)} & \text{when} \ j \neq n\\ -\sum\_{\substack{N\_y\\k=1, k\neq j}}^{N\_y} A\_{jk}^{(1)} & \text{when} \ j = n \end{cases} \quad j, n = 1, 2, \dots, N\_y \tag{29b}$$

with

$$\xi(\mathbf{x}\_i) = \prod\_{k=1, k \neq i}^{N\_x} (\mathbf{x}\_i - \mathbf{x}\_k) \tag{30a}$$

$$\xi(y\_j) = \prod\_{k=1, k \neq j}^{N\_y} (y\_j - y\_k) \tag{30b}$$

and

$$B\_{im}^{(2)} = 2\left(A\_{ii}^{(1)}A\_{im}^{(1)} - \frac{A\_{im}^{(1)}}{(\mathbf{x}\_i - \mathbf{x}\_m)}\right) \quad i, m = 1, 2, \dots, N\_\mathbf{x} \; , \; i \neq m \tag{31a}$$

$$B\_{ji}^{(2)} = 2\left(A\_{jj}^{(1)}A\_{ji}^{(1)} - \frac{A\_{ji}^{(1)}}{(y\_j - y\_n)}\right) \quad j, n = 1, 2, \dots, N\_{\mathcal{Y}} \; , \; j \neq n \tag{31b}$$

$$B\_{ii}^{(2)} = -\sum\_{k=1, k \neq i}^{N\_x} B\_{ik}^{(2)} \quad , i = 1, 2, \dots, N\_{\mathbf{x}\_i} \text{ i } = m \tag{31c}$$

$$B\_{jj}^{(2)} = -\sum\_{\substack{k=1, k \neq j}}^{N\_y} B\_{jk}^{(2)} \quad j = 1, 2, \dots, N\_{y\_{\prime}} \ j = n \tag{31d}$$

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

$$\mathbf{x}\_{i} = \frac{a}{2} \left( 1 - \cos \left( \frac{(i-1)}{(N\_{\mathbf{x}} - 1)} \pi \right) \right) \qquad i = 1, 2, 3, \dots, N\_{\mathbf{x}} \tag{32a}$$

$$\mathbf{y}\_{j} = \frac{b}{2} \left( 1 - \cos\left(\frac{(j-1)}{\left(N\_{y} - 1\right)}\pi\right) \right) \quad j = 1, 2, 3, \dots, N\_{y} \tag{32b}$$

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

$$\left\{ \left[ \begin{array}{c} \left[ M\_{dd} \right] \\ \left[ M\_{bd} \right] \end{array} \begin{array}{c} \left[ M\_{db} \right] \\ \left[ M\_{bb} \right] \end{array} \right] \omega\_{mn}^{2} + \left[ \begin{array}{c} \left[ K\_{dd} \right] \\ \left[ K\_{bd} \right] \end{array} \begin{array}{c} \left[ K\_{db} \right] \\ \left[ K\_{bb} \right] \end{array} \right] \right\} \left\{ \begin{array}{c} \delta\_{d} \\ \delta\_{b} \end{array} \right\} = 0 \tag{33}$$

where we distinguish among inner and boundary grid-points, by means of subscripts *d* and *b*, respectively, whereas *δ* 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(\Phi\_{\mathfrak{e}} + \Phi\_{\mathfrak{w}1} + \Phi\_{\mathfrak{w}2} + \Phi\_{\mathfrak{w}3} + \Phi\_{\mathfrak{w}4}) = 0 \tag{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\_{\rm III} = \sum\_{n=1}^{\infty} \sum\_{m=1}^{\infty} q\_{\rm III} \sin(p\_m \mathbf{x}) \cos(p\_n \mathbf{y}) \tag{35}$$

in which *qmn* = <sup>4</sup>*p*<sup>0</sup> *mnπ*<sup>2</sup> 1 − (−1) *<sup>n</sup>*<sup>1</sup> <sup>−</sup> (−1) *m* , and *p*<sup>0</sup> = 0.1 MPa.

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

$$\{\{F\} + [K]\} \{\delta\} = 0 \tag{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}{c} \left[ F\_{dd} \right] \\ \left[ F\_{bd} \right] \end{array} \right] \begin{array}{c} \left[ F\_{db} \right] \\ \left[ F\_{bb} \right] \end{array} \right\} + \left[ \begin{array}{c} \left[ K\_{dd} \right] \\ \left[ K\_{bd} \right] \end{array} \begin{array}{c} \left[ K\_{db} \right] \\ \left[ K\_{bb} \right] \end{array} \right] \right\} \left\{ \begin{array}{c} \delta\_{d} \\ \delta\_{b} \end{array} \right\} = 0 \tag{37}$$

depending on the kinematic unknowns *δ<sup>d</sup>* and *δ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. 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.

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


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 Figures 4 and 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. 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.

**Figure 4.** Convergence study of the first fundamental frequency, when assuming a GPL-UD, *b*/*a* = 5, *a*/*h* = 10, and Λ*GPL* = 0.5%: (**a**) CCCC boundary conditions, (**b**) SSSS boundary conditions.

**Figure 5.** Convergence study of the bending deflection, when assuming a GPL-UD, *p*<sup>0</sup> = 105, *b*/*a* = 5, *a*/*h* = 10, and Λ*GPL* = 0.5%: (**a**) CCCC boundary conditions, (**b**) SSSS boundary conditions.

**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* × *b* × *h* = 0.45 m × 0.45 m × 0.045 m), *ωmn* = *ωmnh* -*ρm*/*Em*, and Λ*GPL* = 1 %.


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{aligned} E(z) &= (E\_{\mathfrak{c}} - E\_m) \left( \frac{z}{\mathfrak{h}} + \frac{1}{2} \right)^p + E\_m \\ \rho(z) &= (\rho\_{\mathfrak{c}} - \rho\_m) \left( \frac{z}{\mathfrak{h}} + \frac{1}{2} \right)^p + \rho\_m \\ v(z) &= (v\_{\mathfrak{c}} - v\_m) \left( \frac{z}{\mathfrak{h}} + \frac{1}{2} \right)^p + v\_m \end{aligned} \tag{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, 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, 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 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 Tables 5–7, 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 Λ*GPL* = 0.1 (wt%), with a reasonable kinematic response with the selected simply-supported boundary condition.

**Table 3.** Material properties of the system.


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


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


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



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

**Figure 6.** The first six mode shapes of a GPLRC rectangular plate with a GPL-O distribution pattern.

#### *6.2. Parametric Study*

We perform, now, a parametric study focusing on the sensitivity of the frequency response Δ*ω*% = *ω<sup>c</sup>* − *ωepoxy* /*ωepoxy* × 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 Λ*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 *NL* 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 *NL*, 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 Λ*GPL* = 0.3% under a completely-clamped (CCCC) and simplysupported (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 *NL* = 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 Δ*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 *NL* 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 *NL* = 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 Figures 14 and 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 simplysupports, while assuming a Kerr medium in lieu of a Winkler–Pasternak-type foundation (compare the plots of Figures 14 and 15).

**Figure 7.** Relative frequency variation vs. the GPL weight fraction, for different GPL distribution patterns (*b*/*a* = 1, *a*/*h* = 10, and SSSS boundary condition).

**Figure 8.** Relative frequency variation vs. the number of layers, for different GPL distribution patterns (*b*/*a* = 1, *a*/*h* = 10, and SSSS boundary condition).

**Figure 9.** Dimensionless fundamental frequency vs. the number of layers, for different GPL distribution patterns (*b*/*a* = 1, *a*/*h* = 25, Δ*GPL* = 0.3%): (**a**) CCCC boundary conditions, (**b**) SSSS boundary conditions.

**Figure 10.** Dimensionless fundamental frequency vs. thermal variation, for different GPL weight fractions (*a* = *b*, *a*/*h* = 25, and GPL-X).

**Figure 11.** Dimensionless fundamental frequency vs. thermal variation, for different Kerr substrate coefficients (*a* = *b*, *a*/*h* = 25, ΛGPL = 0.3%, Kl = 105, and GPL-X): (**a**) SSSS boundary conditions, (**b**) CCCC boundary conditions.

**Figure 12.** Dimensionless deflection vs. number of layers, for different distribution patterns of reinforcement (b/a = 1, *a*/*h* = 25, Δ*GPL* = 0.3%, *p*<sup>0</sup> = 105): (**a**) SSSS boundary conditions, (**b**) CCCC boundary conditions.

**Figure 13.** Dimensionless deflection changes vs. GPL weight fraction, for different distribution patterns of reinforcement (*b*/*a* = 1, *a*/*h* = 10, *p*<sup>0</sup> = 105): (**a**) SSSS boundary conditions, (**b**) CCCC boundary conditions.

**Figure 14.** Dimensionless deflection vs. GPL weight fraction, for different elastic foundation coefficients (*b*/*a* = 1, *<sup>a</sup>*/*<sup>h</sup>* <sup>=</sup> 10, *GPL* <sup>−</sup> *<sup>X</sup>*, *<sup>p</sup>*<sup>0</sup> <sup>=</sup> <sup>10</sup>5): (**a**) SSSS boundary conditions, (**b**) CCCC boundary conditions.

**Figure 15.** Dimensionless deflection vs. GPL weight fraction, for different Kerr foundation coefficients (*b*/*a* = 1, *<sup>a</sup>*/*<sup>h</sup>* <sup>=</sup> 10, *GPL* <sup>−</sup> *<sup>X</sup>*, *<sup>p</sup>*<sup>0</sup> <sup>=</sup> <sup>10</sup><sup>5</sup> and KI <sup>=</sup> <sup>10</sup>5 : (**a**) SSSS boundary conditions, (**b**) CCCC boundary conditions.

#### **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.

**Author Contributions:** Conceptualization, M.S., A.F., R.D. and F.T.; methodology, M.S., A.F., R.D. and F.T.; software, M.S. and A.F.; validation, M.S., A.F., R.D. and F.T.; formal analysis, M.S. and A.F.; investigation, M.S., A.F., R.D. and F.T.; writing—original draft preparation, M.S. and A.F.; writing—review and editing, R.D. and F.T.; visualization, M.S., A.F., R.D. and F.T.; supervision, R.D. and F.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The non-null strain components are defined as

where

$$\begin{Bmatrix} \varepsilon\_{xx} \\ \varepsilon\_{yy} \\ \gamma\_{yz} \\ \gamma\_{zx} \\ \gamma\_{xy} \end{Bmatrix} = \begin{Bmatrix} \varepsilon\_{xx}^{0} \\ \varepsilon\_{yy}^{0} \\ \gamma\_{yz}^{0} \\ \gamma\_{xy}^{0} \end{Bmatrix} + z \begin{Bmatrix} k\_{xx}^{0} \\ k\_{yy}^{0} \\ k\_{yz}^{0} \\ k\_{xy}^{0} \end{Bmatrix} + z^{2} \begin{Bmatrix} \varepsilon\_{xx}^{\*} \\ \varepsilon\_{yy}^{\*} \\ \gamma\_{yz}^{\*} \\ k\_{xy}^{\*} \end{Bmatrix} + z^{3} \begin{Bmatrix} k\_{xx}^{\*} \\ k\_{yy}^{\*} \\ k\_{yz}^{\*} \\ k\_{xy}^{\*} \end{Bmatrix} \tag{A1}$$

$$\begin{Bmatrix}\begin{Bmatrix}\varepsilon\_{xx}^{0}\\ \varepsilon\_{yy}^{0}\\ \gamma\_{xz}^{0}\\ \gamma\_{xz}^{0}\\ \gamma\_{xy}^{0}\end{Bmatrix}\\ \begin{array}{c}\varepsilon\_{yx}^{0}\\ \gamma\_{yz}^{0}\\ \gamma\_{zx}^{0}\end{array}\end{Bmatrix} = \left\{ \begin{array}{c}\frac{\mathrm{\dot{\boldsymbol{\alpha}}}{\mathrm{\dot{\alpha}}}{\mathrm{\dot{\alpha}}}\\ \frac{\mathrm{\dot{\alpha}}{\mathrm{\dot{\alpha}}}\\ \mathrm{\dot{\alpha}}\\ \frac{\mathrm{\dot{\alpha}}{\mathrm{\dot{\alpha}}}\\ \mathrm{\dot{\alpha}}\end{array}\end{Bmatrix} \left\{ \begin{array}{c}k\_{xx}^{0}\\ k\_{yy}^{0}\\ k\_{yz}^{0}\\ k\_{zx}^{0}\\ \mathrm{\dot{\alpha}}\\ \end{array}\right\} + \left\{\begin{array}{c}\frac{\mathrm{\dot{\alpha}}{\mathrm{\dot{\alpha}}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \end{array}\right\} + \left\{\begin{array}{c}\frac{\mathrm{\dot{\alpha}}{\mathrm{\dot{\alpha}}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \end{array}\right\} = \left\{ \begin{array}{c}\frac{\mathrm{\dot{\alpha}}{\mathrm{\dot{\alpha}}}\\ \frac{\mathrm{\dot{\alpha}}{\mathrm{\dot{\alpha}}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \end{array}\right\} + \left\{ \begin{array}{c}\frac{\mathrm{\dot{\alpha}}{\mathrm{\dot{\alpha}}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\dot{\alpha}}\\ \mathrm{\$$

The elastic constants used in Equation (7) are introduced as

$$Q\_{11} = \frac{E\_{\mathbb{C}}}{1 - \nu\_{\mathbb{C}}^2},\ Q\_{12} = \frac{\nu\_{\mathbb{C}} E\_{\mathbb{C}}}{1 - \nu\_{\mathbb{C}}^2},\ Q\_{13} = Q\_{23} = Q\_{12},\ Q\_{33} = Q\_{22} = Q\_{11},\ Q\_{44} = Q\_{55} = Q\_{66} = G\_{\mathbb{C}}\tag{A3}$$

The set of governing associated to the Hamilton's principle takes the following form

$$
\delta\mathfrak{u}\_0: \frac{\partial N\_{xx}}{\partial \mathbf{x}} + \frac{\partial N\_{xy}}{\partial y} = I\_0 \ddot{\mathfrak{u}}\_0 + I\_1 \ddot{\mathfrak{u}}\_1 + I\_2 \ddot{\mathfrak{u}}\_2 + I\_3 \ddot{\mathfrak{u}}\_3 \tag{A4}
$$

$$
\delta v\_0 : \frac{\partial N\_{yy}}{\partial y} + \frac{\partial N\_{xy}}{\partial x} = I\_0 \ddot{v}\_0 + I\_1 \ddot{v}\_1 + I\_2 \ddot{v}\_2 + I\_3 \ddot{v}\_3 \tag{A5}
$$

$$\begin{array}{c} \delta w\_{0}: \quad \frac{\partial Q\_{x}}{\partial x} + \frac{\partial Q\_{y}}{\partial y} - N\_{1}^{T} \frac{\partial^{2} w\_{0}}{\partial x^{2}} - N\_{2}^{T} \frac{\partial^{2} w\_{0}}{\partial y^{2}} + k\_{p} \frac{\partial^{2} w\_{0}}{\partial x^{2}} + k\_{p} \frac{\partial^{2} w\_{0}}{\partial y^{2}} - k\_{\partial} w\_{0} \\ \quad - \frac{k\_{l} k\_{u}}{k\_{l} + k\_{u}} w\_{0} + \frac{k\_{s} k\_{u}}{k\_{l} + k\_{u}} \frac{\partial^{2} w\_{0}}{\partial x^{2}} + \frac{k\_{s} k\_{u}}{k\_{l} + k\_{u}} \frac{\partial^{2} w\_{0}}{\partial y^{2}} - P w\_{0} = I\_{0} \ddot{w}\_{0} \end{array} \tag{A6}$$

$$
\delta u\_1 : \frac{\partial M\_x}{\partial \mathbf{x}} + \frac{\partial M\_{xy}}{\partial y} - Q\_x = I\_1 \ddot{u}\_0 + I\_2 \ddot{u}\_1 + I\_3 \ddot{u}\_2 + I\_4 \ddot{u}\_3 \tag{A7}
$$

$$
\delta v\_1 : \frac{\partial M\_y}{\partial y} + \frac{\partial M\_{xy}}{\partial x} - Q\_y = I\_1 \ddot{v}\_0 + I\_2 \ddot{v}\_1 + I\_3 \ddot{v}\_2 + I\_4 \ddot{v}\_3 \tag{A8}
$$

$$\delta\mu\_2: \frac{\partial N\_x^\*}{\partial \mathbf{x}} + \frac{\partial N\_{xy}^\*}{\partial y} - 2S\_x = I\_2\ddot{u}\_0 + I\_3\ddot{u}\_1 + I\_4\ddot{u}\_2 + I\_5\ddot{u}\_3\tag{A9}$$

$$\delta v\_2 : \frac{\partial N\_y^\*}{\partial y} + \frac{\partial N\_{xy}^\*}{\partial x} - 2S\_y = I\_2 \ddot{v}\_0 + I\_3 \ddot{v}\_1 + I\_4 \ddot{v}\_2 + I\_5 \ddot{v}\_3 \tag{A10}$$

$$\delta\mu\_3: \frac{\partial M\_x^\*}{\partial x} + \frac{\partial M\_{xy}^\*}{\partial y} - 3Q\_x^\* = I\_3\ddot{u}\_0 + I\_4\ddot{u}\_1 + I\_5\ddot{u}\_2 + I\_6\ddot{u}\_3\tag{A11}$$

$$\delta v\_{\Im} \colon \frac{\partial M\_y^\*}{\partial y} + \frac{\partial M\_{xy}^\*}{\partial x} - 3Q\_y^\* = I\_3\ddot{v}\_0 + I\_4\ddot{v}\_1 + I\_5\ddot{v}\_2 + I\_6\ddot{v}\_3 \tag{A12}$$

where .. ( ◦ ) refers to the acceleration field, and *Ii* are the mass inertias. The corresponding boundary conditions are defined as

$$\begin{array}{ll} \delta u\_{0} = 0 & \text{or} & N\_{\text{x}} \ n\_{x} + N\_{xy} n\_{y} = 0 \\ \delta v\_{0} = 0 & \text{or} & N\_{\text{y}} \ n\_{y} + N\_{xy} n\_{x} = 0 \\ \delta w\_{0} = 0 & \text{or} & N\_{\text{x}} \ n\_{x} + Q\_{y} u\_{y} + Q\_{x} n\_{x} + N\_{1}^{T} \frac{\partial w\_{0}}{\partial x} n\_{x} + N\_{2}^{T} \frac{\partial w\_{0}}{\partial y} n\_{y} \\ & & -k\_{p} \frac{\partial w\_{0}}{\partial x} n\_{x} - k\_{p} \frac{\partial w\_{0}}{\partial y} n\_{y} - \frac{k\_{x} k\_{y}}{k\_{l} + k\_{x}} \frac{\partial w\_{0}}{\partial x} n\_{x} - \frac{k\_{x} k\_{y}}{k\_{l} + k\_{x}} \frac{\partial w\_{0}}{\partial y} n\_{y} = 0 \\\\ \delta u\_{1} = 0 & \text{or} & M\_{\text{x}} \ n\_{x} + M\_{xy} n\_{y} = 0 \\ \delta v\_{1} = 0 & \text{or} & M\_{\text{y}} \ n\_{y} n\_{x} + M\_{xy} n\_{x} = 0 \end{array} \tag{A14}$$

$$\begin{array}{ll} \delta u\_2 = 0 & \text{or} \quad \begin{array}{l} N\_x^\* \ n\_x + N\_{xy}^\* n\_y = 0 \\ N\_y^\* \ n\_y + N\_{xy}^\* n\_x = 0 \end{array} . \end{array} \tag{A15}$$

$$\begin{array}{ll}\delta\mu\_3 = 0 & \text{or} & M\_{\text{x}}^{\*}\ n\_{\text{x}} + M\_{\text{xy}}^{\*}\ n\_{\text{y}} = 0\\ \delta\upsilon\_3 = 0 & \text{or} & M\_{\text{y}}^{\*}\ n\_{\text{y}}\upsilon\_{\text{x}} + M\_{\text{xy}}^{\*}\upsilon\_{\text{x}} = 0 \end{array} \tag{A16}$$

where

⎡ ⎣ Nx N∗ x Ny N∗ y Nxy N∗ xy ⎤ <sup>⎦</sup> <sup>=</sup> *h* 2 − *h* 2 ⎧ ⎨ ⎩ *σ*x *σ*y *σ*xy ⎫ ⎬ ⎭ # 1 z<sup>2</sup> \$ dz , ⎧ ⎨ ⎩ Mx My Mxy ⎫ ⎬ ⎭ = *h* 2 − *h* 2 ⎧ ⎨ ⎩ *σ*x *σ*y *σ*xy ⎫ ⎬ ⎭ zdz , ⎧ ⎨ ⎩ M∗ x M∗ y M∗ xy ⎫ ⎬ ⎭ = *h* 2 − *h* 2 ⎧ ⎨ ⎩ *σ*x *σ*y *σ*xy ⎫ ⎬ ⎭z3dz , ! Qx Qy " = *h* 2 − *h* 2 \* *σ*xz *<sup>σ</sup>*yz , dz ,! Sx Sy " = *h* 2 − *h* 2 \* *σ*xz *<sup>σ</sup>*yz , z dz , (I0, I1, I2, I3, I4, I5, I6) = z *ρ* 1, z1, z2, z3, z4, z5, z6 dz. (A17)

and

*Nxx* = *A*<sup>11</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>B</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>C</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>A</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>B</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>C</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>D</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>3</sup> *∂y Nxy* <sup>=</sup> *<sup>A</sup>*44-*∂u*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>0</sup> *∂x* <sup>+</sup> *<sup>B</sup>*44-*∂u*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>1</sup> *∂x* <sup>+</sup> *<sup>C</sup>*44-*∂u*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>2</sup> *∂x* <sup>+</sup> *<sup>D</sup>*44-*∂u*<sup>3</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>3</sup> *∂x Nyy* = *A*<sup>21</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>B</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>C</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>A</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>B</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>C</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>D</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>3</sup> *∂y Qy* <sup>=</sup> *<sup>A</sup>*55-*∂w*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* + *v*<sup>1</sup> + *B*55(2*v*2) + *C*55(3*v*3) *Qx* <sup>=</sup> *<sup>A</sup>*66-*∂w*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* + *u*<sup>1</sup> + *B*66(2*u*2) + *C*66(3*u*3) *Mx* = *B*<sup>11</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>C</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>E</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>B</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>C</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>D</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>E</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>3</sup> *∂y Mxy* <sup>=</sup> *<sup>B</sup>*44-*∂u*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>0</sup> *∂x* <sup>+</sup> *<sup>C</sup>*44-*∂u*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>1</sup> *∂x* <sup>+</sup> *<sup>D</sup>*44-*∂u*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>2</sup> *∂x* <sup>+</sup> *<sup>E</sup>*44-*∂u*<sup>3</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>3</sup> *∂x My* = *B*<sup>21</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>C</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>E</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>B</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>C</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>D</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>E</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>3</sup> *∂y N*∗ *<sup>x</sup>* = *<sup>C</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>E</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>F</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>C</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>D</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>E</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>F</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>3</sup> *∂y N*∗ *xy* <sup>=</sup> *<sup>C</sup>*44-*∂u*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>0</sup> *∂x* <sup>+</sup> *<sup>D</sup>*44-*∂u*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>1</sup> *∂x* <sup>+</sup> *<sup>E</sup>*44-*∂u*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>2</sup> *∂x* <sup>+</sup> *<sup>F</sup>*44-*∂u*<sup>3</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>3</sup> *∂x N*∗ *<sup>y</sup>* = *<sup>C</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>E</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>F</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>C</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>D</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>E</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>F</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>3</sup> *∂y M*∗ *<sup>x</sup>* = *<sup>D</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>E</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>F</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>Y</sup>*<sup>11</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>E</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>F</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>Y</sup>*<sup>12</sup> *<sup>∂</sup>v*<sup>3</sup> *∂y M*∗ *xy* <sup>=</sup> *<sup>D</sup>*44-*∂u*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>0</sup> *∂x* <sup>+</sup> *<sup>E</sup>*44-*∂u*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>1</sup> *∂x* <sup>+</sup> *<sup>F</sup>*44-*∂u*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>2</sup> *∂x* <sup>+</sup> *<sup>Y</sup>*44-*∂u*<sup>3</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>∂</sup>v*<sup>3</sup> *∂x M*∗ *<sup>y</sup>* = *<sup>D</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>E</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>F</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>2</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>Y</sup>*<sup>21</sup> *<sup>∂</sup>u*<sup>3</sup> *<sup>∂</sup><sup>x</sup>* <sup>+</sup> *<sup>D</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>0</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>E</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>F</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>2</sup> *<sup>∂</sup><sup>y</sup>* <sup>+</sup> *<sup>Y</sup>*<sup>22</sup> *<sup>∂</sup>v*<sup>3</sup> (A18)

Based on relations (A13)–(A16), different boundary conditions can be set as follows

$$\begin{array}{c} \mathbf{C} \mathbf{A} \mathbf{m} \mathbf{p} \mathbf{e} \mathbf{ (C)} \mathbf{e} \mathbf{d} \mathbf{g} \mathbf{e} \mathbf{s} \\\quad \left\{ \begin{array}{ll} \mathbf{x} = \mathbf{0}, \mathbf{a} \\\ y = \mathbf{0}, \mathbf{b} \end{array} \rightarrow \begin{cases} u\_0 = 0 & u\_1 = 0 \\ v\_0 = 0 & v\_1 = 0 \\\ w\_1 = 0 & v\_2 = 0 \end{cases} \nu\_3 = 0 \\\quad \text{(A19)} \\\hline \\ w\_0 = 0 \end{array} \tag{A19}$$

*∂y*

$$\begin{array}{rcl} \textbf{Simply (S) edges:} \begin{cases} \textbf{x} = \textbf{0, a} & \rightarrow \begin{cases} N\_{\textbf{xx}} = M\_{\textbf{xx}} = N\_{\textbf{xx}}^{\*} = M\_{\textbf{xx}}^{\*} = \textbf{0} \\\ \textbf{y} = \textbf{0, b} & \rightarrow \begin{cases} N\_{\textbf{yy}} = M\_{\textbf{yy}} = N\_{\textbf{yy}}^{\*} = M\_{\textbf{yy}}^{\*} = \textbf{0} \end{cases} \end{array} . \end{array} . \end{array} \tag{A20}$$

#### **Appendix B**

The parameters in Equations (21a)–(21c) are defined as follows

$$D = \frac{1}{3} V\_{GPL}^\* \left\{ \frac{2}{H + \frac{1}{\frac{K\_x}{K\_W} - 1}} + \frac{1}{\frac{1}{2}(1 - H) + \frac{1}{\frac{K\_x}{K\_W} - 1}} \right\},\tag{A21}$$

$$H = \frac{\text{Ln}\left[\rho\left(\rho + \sqrt{\rho^2 - 1}\right)\right]}{\sqrt{\left(\rho^2 - 1\right)^3}} - \frac{1}{\rho^2 - 1},\tag{A22}$$

$$K\_{\chi} = \frac{K\_{\chi}}{\frac{2R\_{k}K\_{\chi}}{L} + 1},\tag{A23}$$

$$K\_z = \frac{K\_\%}{\frac{2R\_kK\_\chi}{t} + 1},\tag{A24}$$

$$
\rho = \frac{l\_{GPL}}{t\_{GPL}}.\tag{A25}
$$

The constants in Equation (23) are determined as

$$\mathbf{C}\_{11} = \frac{\begin{vmatrix} T\_1 & T\_2 \\ 1 & e^{-\sqrt{A\_1}h} \end{vmatrix}}{\begin{vmatrix} 1 & 1 \\ e^{\sqrt{A\_1}h} & e^{-\sqrt{A\_1}h} \end{vmatrix}}, \quad \mathbf{C}\_{22} = \frac{\begin{vmatrix} T\_2 & T\_1 \\ e^{\sqrt{A\_1}h} & 1 \end{vmatrix}}{\begin{vmatrix} 1 & 1 \\ e^{\sqrt{A\_1}h} & e^{-\sqrt{A\_1}h} \end{vmatrix}}, \quad A\_{11} = P\_m^2 + P\_{n\prime}^2 \tag{A26}$$

where

*FG* − *X*: *Z* < *<sup>h</sup>* <sup>2</sup> <sup>⇒</sup> *<sup>A</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>4*<sup>D</sup> <sup>h</sup>* ; *<sup>A</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>+</sup> <sup>2</sup>*<sup>D</sup>* ; *<sup>A</sup>*<sup>3</sup> <sup>=</sup> 0 ; *<sup>A</sup>*<sup>4</sup> <sup>=</sup> <sup>−</sup>4*<sup>D</sup> <sup>h</sup>* ; *<sup>A</sup>*<sup>5</sup> <sup>=</sup> <sup>−</sup>*A*<sup>11</sup> <sup>4</sup>*<sup>D</sup> <sup>h</sup>* ; *A*<sup>6</sup> = (1 + 2*D* )*A*<sup>11</sup> *Z* > *<sup>h</sup>* <sup>2</sup> <sup>⇒</sup> *<sup>A</sup>*<sup>1</sup> <sup>=</sup> <sup>4</sup>*<sup>D</sup> <sup>h</sup>* ; *<sup>A</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> <sup>2</sup>*<sup>D</sup>* ; *<sup>A</sup>*<sup>3</sup> <sup>=</sup> 0 ; *<sup>A</sup>*<sup>4</sup> <sup>=</sup> <sup>4</sup>*<sup>D</sup> <sup>h</sup>* ; *<sup>A</sup>*<sup>5</sup> <sup>=</sup> *<sup>A</sup>*<sup>11</sup> <sup>4</sup>*<sup>D</sup> <sup>h</sup>* ; *A*<sup>6</sup> = (1 − 2*D* )*A*<sup>11</sup> (A27)

*FG* − *O* :

 $Z < \frac{h}{2} \implies A\_1 = \frac{4D}{h};\ A\_2 = 1;\ A\_3 = 0$  ;\ A\_4 = \frac{4D}{h};\ A\_5 = A\_{11} $\frac{4D}{h};\ A\_6 = A\_{11} \tag{A28}$  $Z > \frac{h}{2} \implies A\_1 = -\frac{4D}{h};\ A\_2 = 1 + 4D$  ;\ A\_3 = 0 ;\ A\_4 = -\frac{4D}{h};\ A\_5 = -\frac{4D}{h} $A\_{11};\ A\_6 = (1 + 4D)A\_{11}$ 

More details about coefficients in Equations (25a)–(25c) are defined in the following

$$\alpha\_1 = -\frac{\left(A\_3 + \sqrt{-4A\_5A\_1 + A\_3^2}\right)}{2A\_1}Z\tag{A29}$$

$$\mu\_2 = \frac{A\_1^{\cdot^2} - A\_1 A\_4 + A\_2 A\_3}{A\_1^{\cdot^2}} \tag{A30}$$

$$\begin{array}{ll} \beta\_{1} &= \frac{-2A\_{6}A\_{1}^{2} - 2A\_{1}^{2}\sqrt{-4A\_{5}A\_{1} + A\_{3}^{2}} - 2A\_{2}A\_{5}A\_{1}}{2A\_{1}^{2}\sqrt{-4A\_{5}A\_{1} + A\_{3}^{2}}} \\ &+ \frac{-A\_{1}A\_{3}A\_{4} + A\_{1}A\_{4}\sqrt{-4A\_{5}A\_{1} + A\_{3}^{2}} + A\_{2}A\_{3}^{2} - A\_{2}A\_{3}\sqrt{-4A\_{5}A\_{1} + A\_{3}^{2}}}{2A\_{1}^{2}\sqrt{-4A\_{5}A\_{1} + A\_{3}^{2}}} \\ \end{array} \tag{A31}$$

$$\beta\_2 = \frac{2A\_1^{\cdot^2} - A\_1 A\_4 + A\_2 A\_3}{A\_1^{\cdot^2}} \tag{A32}$$

$$\beta\_3 = \frac{\sqrt{-4A\_5A\_1 + A\_3^2}(A\_2 + A\_1Z)}{A\_1^2} \tag{A33}$$

The governing equations of GPLRC plates determined by means of minimum total potential energy principle in Section 5 are defined as

$$
\delta u\_0 : \frac{\partial N\_{xx}}{\partial x} + \frac{\partial N\_{xy}}{\partial y} = 0 \,, \tag{A34}
$$

$$
\delta v\_0 : \frac{\partial N\_{yy}}{\partial y} + \frac{\partial N\_{xy}}{\partial x} = 0,\tag{A35}
$$

$$\delta w\_{0} : \begin{array}{l} \frac{\partial Q\_{\text{t}}}{\partial x} + \frac{\partial Q\_{\text{y}}}{\partial y} - N\_{1}^{T} \frac{\partial^{2} w\_{0}}{\partial x^{2}} - N\_{2}^{T} \frac{\partial^{2} w\_{0}}{\partial y^{2}} + k\_{p} \frac{\partial^{2} w\_{0}}{\partial x^{2}} + k\_{p} \frac{\partial^{2} w\_{0}}{\partial y^{2}} - k\_{\text{uv}} w\_{0} \\ - \frac{k\_{l} k\_{u}}{k\_{l} + k\_{u}} w\_{0} + \frac{k\_{v} k\_{u}}{k\_{l} + k\_{u}} \frac{\partial^{2} w\_{0}}{\partial x^{2}} + \frac{k\_{v} k\_{u}}{k\_{l} + k\_{u}} \frac{\partial^{2} w\_{0}}{\partial y^{2}} - P w\_{0} = 0, \end{array} \tag{A36}$$

$$
\delta u\_1 : \frac{\partial M\_x}{\partial x} + \frac{\partial M\_{xy}}{\partial y} - Q\_x = 0,\tag{A37}
$$

$$
\delta v\_1 : \frac{\partial M\_y}{\partial y} + \frac{\partial M\_{xy}}{\partial x} - Q\_y = 0,\tag{A38}
$$

$$
\delta\mu\_2: \frac{\partial N\_x^\*}{\partial x} + \frac{\partial N\_{xy}^\*}{\partial y} - 2S\_x = 0 \,, \tag{A39}
$$

$$
\delta v\_2 : \frac{\partial N\_y^\*}{\partial y} + \frac{\partial N\_{xy}^\*}{\partial x} - 2S\_y = 0 \,, \tag{A40}
$$

$$
\delta\mathfrak{u}\_3 : \frac{\partial M\_x^\*}{\partial \mathfrak{x}} + \frac{\partial M\_{xy}^\*}{\partial y} - \mathfrak{H}\_x^\* = 0 \,, \tag{A41}
$$

$$
\delta v\_3 : \frac{\partial M\_y^\*}{\partial y} + \frac{\partial M\_{xy}^\*}{\partial x} - \mathfrak{H}\_y^\* = 0 \,, \tag{A42}
$$

The boundary conditions and static quantities associated with the problem are the same as defined in Equations (A13)–(A16).

#### **References**

