## 2 Theoretical Definition of the Problem
### Governing Equations  
The governing equations of the problem are derived here from the principle of virtual work, as follows:  
\[\delta U-\delta V_{g}=0 \tag{5}\]  
where \(U\) is the total strain potential energy of the plate defined on the domain \(\Omega\) as:  
\[U=\frac{1}{2}\int_{\Omega}\sigma_{ij}\varepsilon_{ij}d\Omega \tag{6}\]  
and \(V_{g}\) is the potential energy related to geometry, which takes the following form:  
\[V_{g}=\frac{1}{2}\int_{\Omega}\left\{\begin{array}{l}P_{x}\Bigg{[}\left( \frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial w}{\partial x} \right)^{2}+\left(\frac{\partial w}{\partial x}\right)^{2}\Bigg{]}+P_{y} \Bigg{[}\left(\frac{\partial u}{\partial y}\right)^{2}+\left(\frac{\partial w }{\partial y}\right)^{2}+\left(\frac{\partial v}{\partial y}\right)^{2} \Bigg{]}+\\ P_{xy}\Bigg{(}\frac{\partial w}{\partial x}\frac{\partial w}{\partial y}+ \frac{\partial v}{\partial x}\frac{\partial v}{\partial y}-\frac{\partial u }{\partial x}\frac{\partial u}{\partial y}\Bigg{)}\end{array}\right\}d\Omega \tag{7}\]  
By substitution of Equations (6) and (7) into Equation (5), the following relation is obtained:  
\[\int_{0}^{h}\int_{0}^{a}\int_{0}^{b}\left[\begin{array}{l}\left(\sigma_{xx} \delta\varepsilon_{xx}+\sigma_{yy}\delta\varepsilon_{yy}+\sigma_{zz}\delta \varepsilon_{zz}+\sigma_{yz}\delta\gamma_{yz}+\sigma_{xz}\delta\gamma_{xz}+ \sigma_{xy}\delta\gamma_{xy}\right)\\ -P_{x}\bigg{(}\frac{\partial u}{\partial x}\frac{\partial\delta u}{\partial x }+\frac{\partial v}{\partial x}\frac{\partial\delta v}{\partial x}+\frac{ \partial w}{\partial x}\frac{\partial\delta w}{\partial x}\bigg{)}-P_{y} \bigg{(}\frac{\partial u}{\partial y}\frac{\partial\delta u}{\partial y}+ \frac{\partial v}{\partial y}\frac{\partial\delta v}{\partial y}+\frac{ \partial w}{\partial y}\frac{\partial\delta w}{\partial y}\bigg{)}\\ -P_{xy}\bigg{(}\frac{\partial\delta w}{\partial x}\frac{\partial w}{\partial y }+\frac{\partial\delta w}{\partial y}\frac{\partial w}{\partial x}+\frac{ \partial\delta w}{\partial x}\frac{\partial w}{\partial y}+\frac{\partial \delta w}{\partial y}\frac{\partial w}{\partial x}-\frac{\partial\delta u}{ \partial x}\frac{\partial u}{\partial y}-\frac{\partial\delta u}{\partial y} \frac{\partial u}{\partial x}\bigg{)}\end{array}\right]dv=0 \tag{8}\]  
where the constitutive relations for FG-saturated porous plates in 3D poroelasticity can be defined according to [PERSON]'s constitutive law, as \([\sigma_{ij}]=[C][\varepsilon_{ij}]\). For FG-saturated porous rectangular plates, the elasticity matrix reads as follows:  
\[C=\frac{E(1-\
u)}{(1+\
u)(1-2\
u)}\left(\begin{array}{cccccc}1&\frac{\
u}{ 1-\
u}&\frac{\
u}{1\,\frac{\
u}{\
u}}&0&0&0\\ \frac{\
u}{1-\
u}&1&\frac{1}{1-\
u}&0&0&0\\ \frac{\
u}{1-\
u}&\frac{\
u}{1-\
u}&1&0&0&0\\ 0&0&0&\frac{1-2\
u}{2(1-\
u)}&0&0\\ 0&0&0&0&\frac{1-2\
u}{2(1-\
u)}&0\\ 0&0&0&0&0&\frac{1-2\
u}{2(1-\
u)}\end{array}\right)+\left(\begin{array}{ ccccccc}\overline{M}a^{2}&0&0&0&0&0\\ 0&\overline{M}a^{2}&0&0&0&0\\ 0&0&\overline{M}a^{2}&0&0&0\\ 0&0&\overline{M}a^{2}&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\end{array}\right)= \tag{9}\]  
The elasticity modulus, \(E\), is assumed to vary along the \(z\)-direction, whereas the Poisson's ratio, \(\
u\), remains constant.