## 3 Natural Element (NE) Bending and Free Vibration Approximation  
In the NE approximation, a uniform NEM grid is constructed on the mid-surface \(\varpi\) of the FG porous plate. Referring to a previous paper [27], the NEM grid consists of a finite number of grid points called nodes and Delaunay triangles. For the NE approximation  
Figure 2: Uniform, center-, lower-, upper-biased porosity distributions [12]. of the (3,3,2) hierarchical model (9), its actual and virtual displacements \(\mathbf{u}^{h}\) and \(\mathbf{v}^{h}\) are expanded as (\(\alpha,\beta=x,y,z\))  
\[u^{h}_{\alpha}(\mathbf{x})=\sum_{m=0}^{q_{x}}\left(\sum_{J=1}^{N}U^{m}_{a,J}\phi_{J} (x,y)\right)\cdot\left(\frac{2z}{h}\right)^{m} \tag{10}\]  
\[v^{h}_{\beta}(\mathbf{x})=\sum_{\ell=0}^{q_{\beta}}\left(\sum_{I=1}^{N}V^{\ell}_{ \beta,I}\phi_{I}(x,y)\right)\cdot\left(\frac{2z}{h}\right)^{\ell} \tag{11}\]  
with \(q_{x}=q_{y}=3\) and \(q_{z}=2\). Here, \(\phi_{I}(x,y)\) are Laplace interpolation (\(L/I\)) functions [32; 33] defined on the NEM grid and \(U^{m}_{a,J}\) indicates the nodal values to be determined. Letting \(q_{x}=q_{y}=q_{z}=q\) for the concise expression of mathematical formula, the virtual linear elastic strain vector \(\mathbf{\varepsilon}\left(\mathbf{v}^{h}\right)\) and the actual Cauchy stress vector \(\sigma\left(\mathbf{u}^{h}\right)\) are approximated as  
\[\mathbf{\varepsilon}\left(\mathbf{v}^{h}\right)=\sum_{\ell=0}^{q}\sum_{I=1}^{N}\mathbf{L} _{\infty}\phi_{I}V^{\ell}_{I}\cdot\left(\frac{2z}{h}\right)^{\ell}=\sum_{\ell =0}^{q}\sum_{I=1}^{N}\mathbf{B}^{\ell}_{I}\mathbf{V}^{\ell}_{I}\cdot\left(\frac{2z}{h }\right)^{\ell} \tag{12}\]  
\[\mathbf{\sigma}\left(\mathbf{u}^{h}\right)=\mathbf{E}\mathbf{\varepsilon}\left(\mathbf{u}^{h} \right)=\mathbf{E}\sum_{m=0}^{q}\sum_{J=1}^{N}\mathbf{B}^{m}_{I}\mathbf{U}^{m}_{J}\cdot \left(\frac{2z}{h}\right)^{m} \tag{13}\]  
in which the gradient-like matrix \(\mathbf{L}_{\infty}\) and the partial differential matrix \(\mathbf{B}^{m}_{I}\) are respectively defined by  
\[\mathbf{L}_{\infty}= \begin{bmatrix}\partial_{x}&0&0&\partial_{y}&0&m/z\\ 0&\partial_{y}&0&\partial_{x}&m/z&0\\ 0&0&m/z&0&\partial_{y}&\partial_{x}\end{bmatrix} ^{T} \tag{14}\]  
\[\mathbf{B}^{m}_{I}=\mathbf{L}_{\infty}\phi_{I}\] \[=\left[ \begin{array}{cccccc}\phi_{J,x}&0&0&\phi_{J,y}&0&m\phi_{ \tilde{\gamma}I}/z\\ 0&\phi_{J,y}&0&\phi_{J,x}&m\phi_{\tilde{\gamma}I}/z&0\\ 0&0&m\phi_{\tilde{\gamma}I}/z&0&\phi_{\tilde{\gamma}I},&\phi_{\tilde{\gamma}I },x\end{array} \right]^{T} \tag{15}\]  
with \(\partial_{x}=\partial/\partial x\) and \(\phi_{J,x}=\partial\phi_{\tilde{\gamma}I}/\partial x\). \(E\) denotes the \((6\times 6)\) matrix of the elastic modulus \(E\) and the Poisson's ratio \(\
u\). Introducing Equations (12) and (13) into the weak form of static equilibrium [27] of FG porous plates leads to  
\[[\mathbf{K}]^{\ell m}_{I\tilde{\gamma}}\{\mathbf{U}\}^{\ell}_{I}=\{\mathbf{F}\}^{m}_{I} \tag{16}\]  
with the stiffness matrix \([\mathbf{K}]\) and the load vector \(\{\mathbf{F}\}\) defined by  
\[[\mathbf{K}]^{\ell m}_{II}= \tag{17}\]  
\[\{\mathbf{F}\}^{m}_{I}=\int_{\partial\Omega_{N}}\widetilde{t}_{a}\phi_{I}ds\cdot(2 z^{*}/h)^{m}dA \tag{18}\]  
Here, \(\widetilde{t}\) indicates the external force acting on the natural boundary \(\partial\Omega_{N}\) located at the vertical position \(z^{*}\). The subscripts \(FI\) and \(RI\) in Equation (17) denote the full numerical integration using 7 Gauss points and the reduced integration using only 1 Gauss point. The reduced integration is used to overcome shear locking [34; 35] for the bending-prevailed thin elastic structures, for which the material constant matrix \(\mathbf{E}\) is divided into  
\[[\mathbf{E}]= \begin{bmatrix}\mathbf{E}_{1}&0\\ 0&\mathbf{E}_{2}\end{bmatrix} ,\,[\mathbf{E}_{1}]= \begin{bmatrix}C_{1}&C_{2}&C_{2}\\ C_{2}&C_{1}&C_{2}\\ C_{2}&C_{2}&C_{1}\end{bmatrix}