## 2 Theoretical Formulation
### 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_{e}-(\delta\[PERSON]_{w1}+ \delta\Phi_{w2}+\delta\Phi_{w3}+\delta\Phi_{w4})\right)\text{dt}=0 \tag{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}=\int\limits_{V}\rho\bigg{(}\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}\bigg{)}dV \tag{9}\]  
\[\delta\Phi_{e}=\int\limits_{V}\big{(}\sigma_{xx}\delta\varepsilon_{xx}+\sigma _{yy}\delta\varepsilon_{yy}+\sigma_{zz}\delta\varepsilon_{zz}+\tau_{xy} \delta\tau_{xy}+\tau_{yz}\delta\gamma_{yz}+\tau_{xz}\delta\gamma_{xz}\big{)}dV \tag{10}\]  
In addition  
\[\delta\Phi_{w1}=\int\limits_{A}\Big{(}-k_{w}w_{o}+k_{p}\Big{(}\frac{\partial^ {2}w_{o}}{\partial x^{2}}+\frac{\partial^{2}w_{o}}{\partial y^{2}}\Big{)} \Big{)}\delta w_{o}dA \tag{11}\]  
\(k_{w}\) and \(k_{p}\) being the Winkler and Pasternak constants.  
While  
\[\delta\Phi_{w2}=\int\limits_{A}\Big{(}-\frac{k_{s}k_{u}}{k_{l}+k_{w}}w_{o}+ \frac{k_{k}k_{y}}{k_{l}+k_{w}}\Big{(}\frac{\partial^{2}w_{o}}{\partial x^{2}} +\frac{\partial^{2}w_{o}}{\partial y^{2}}\Big{)}\Big{)}\delta w_{o}dA \tag{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}=-\int\limits_{A}P\delta w_{0}dA \tag{13}\]  
In addition, the conductive layer reinforced with GPLs satisfies the following Fourier heat conduction relation:  
\[\
abla^{2}T+R=\rho c\frac{\partial T}{\partial t} \tag{14}\]  
In absence of a thermal generation, in steady-state conditions, it is:  
\[\
abla^{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 UD \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(x,0,z)=0,\ T(x,b,z)=0,\ T(x,y,0)=T_{1},\ T(x,y,h)=T _{2} \tag{17}\]  
The last energy contribution related to the thermal load can be obtained as:  
\[\delta\,\Phi_{\
u 4}=\int\limits_{A}\left((N_{1}^{\rm T})\frac{\partial w_{0}}{ \partial x}\frac{\partial\delta w_{0}}{\partial x}+(N_{2}^{\rm T})\frac{ \partial w_{0}}{\partial y}\frac{\partial\delta w_{0}}{\partial y}\right)dA \tag{18}\]  
with  
\[N_{1}^{\rm T}=\int\limits_{-\frac{h}{2}}^{\frac{h}{2}}(Q_{11}+Q_{12})\alpha_{ C}(T-T_{0})dz,\ N_{2}^{\rm T}=\int\limits_{-\frac{h}{2}}^{\frac{h}{2}}(Q_{21}+Q_{22} )\alpha_{C}(T-T_{0})\ dz \tag{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)).