## 4 Mathematical Modeling  
According to the First-order Shear Deformation Theory (FSDT), the three-dimensional displacements (\(u_{1},u_{2},u_{3}\)) are defined as follows [13; 38]:  
\[ \begin{array}{l}u_{1}=u+\zeta\beta_{1}\\ u_{2}=v+\zeta\beta_{2}\\ u_{3}=w\end{array}  \tag{15}\]  
where \(u\), \(v\), and \(w\) are the displacement components, and \(\beta_{1}\) and \(\beta_{2}\) are the rotations of the tangents to the middle surface about \(\alpha_{1}\) and \(\alpha_{2}\) axes, respectively. These five parameters are the degrees of freedom of the problem and can be collected in the corresponding vector \(U\). The equations of motion can be derived through the [PERSON]'s variational principle, applied within a generic time interval bounded by \(t_{0},t_{1}\):  
\[\int_{t_{0}}^{t_{1}}(\delta U+\delta W-\delta K)dt=0, \tag{16}\]  
where \(\delta U\), \(\delta K\), and \(\delta W\) are the variations of strain energy and kinetic energy, as well as the work done by the external forces, respectively. The strain energy of an elastic body occupying a volume \(V\) is defined as follows:  
\[U=\frac{1}{2}\int(\sigma_{ij}\epsilon_{ij})dV. \tag{17}\]  
The three-dimensional strain components are defined as  
\[\epsilon_{11}=\epsilon_{11}+\zeta k_{11} \tag{18}\]  
\[\epsilon_{22}=\epsilon_{22}+\zeta k_{22} \tag{19}\]  
\[\epsilon_{12}=\frac{1}{2}[\epsilon_{12}+\zeta k_{12}] \tag{20}\]  
\[\gamma_{13}=\beta_{1}+\frac{\partial w}{\partial\alpha_{1}}-\frac{u}{R_{1}} \tag{21}\]  
\[\gamma_{23}=\beta_{2}+\frac{\partial w}{\partial\alpha_{2}}-\frac{v}{R_{2}}, \tag{22}\]  
where  
\[e_{11}=\frac{1}{(1+\zeta/R_{1})}\bigg{[}\frac{\partial u}{\partial\alpha_{1}} +\frac{w}{R_{1}}\bigg{]} \tag{23}\]  
\[e_{22}=\frac{1}{(1+\zeta/R_{2})}\bigg{[}\frac{\partial v}{\partial\alpha_{2}} +\frac{w}{R_{2}}\bigg{]} \tag{24}\]  
\[e_{12}=\frac{1}{(1+\zeta/R_{2})}\bigg{[}\frac{\partial u}{\partial\alpha_{2}} \bigg{]}+\frac{1}{(1+\zeta/R_{1})}\bigg{[}\frac{\partial v}{\partial\alpha_{ 1}}\bigg{]} \tag{25}\]  
\[U=\frac{1}{2}\int(\sigma_{ij}\epsilon_{ij})dV. \tag{17}\]0  
\[U=\frac{1}{2}\int(\sigma_{ij}\epsilon_{ij})dV. \tag{17}\]1  
\[U=\frac{1}{2}\int(\sigma_{ij}\epsilon_{ij})dV. \tag{17}\]2  
under the assumption \(A_{1}=A_{2}=1\), related to the correspondence between curvilinear abscissae and principal coordinates. The stress-strain relationships at a given point within the FGM porous doubly-curved shell element take the following form:  
\[\left\{ \begin{array}{c}\sigma_{11}\\ \sigma_{22}\\ \sigma_{12}\\ \sigma_{23}\\ \sigma_{13}\end{array} \right\}=\left[ \begin{array}{cccccc}C_{11}&C_{12}&0&0& 0\\ C_{12}&C_{22}&0&0&0\\ 0&0&C_{66}&0&0\\ 0&0&0&C_{44}&0\\ 0&0&0&0&C_{55}\end{array} \right]\left\{ \begin{array}{c}\epsilon_{11}\\ \epsilon_{22}\\ \gamma_{12}\\ \gamma_{23}\\ \gamma_{13}\end{array} \right\}, \tag{29}\]  
where the three-dimensional stress components are indicated on the left. Here, \(C_{ij}\) are the elastic coefficients included in the corresponding elasticity matrix. They are defined below as  
\[C_{11}=\frac{E(\zeta)}{1-\
u(\zeta)^{2}},C_{12}=\
u(\zeta)C_{11},C_{22}=C_{11}, \ C_{44}=C_{55}=C_{66}=\frac{E(\zeta)}{2(1+\
u(\zeta))}. \tag{30}\]  
The kinetic energy \(K\) can be written as follows:  
\[K=\frac{1}{2}\int_{V}\rho(\zeta)\dot{U}^{T}\dot{U}dV, \tag{31}\]  
where \(U\) is the displacement vector, which collects the three-dimensional displacement components introduced above. The dot notation used in the expression above stands for the time derivative. In other words, \(\dot{U}\) collects the velocity components. On the contrary, the variations of the work done by external forces is given by  
\[\delta W=\int_{\alpha_{2}}\int_{\alpha_{1}}(q_{1}\delta u+q_{2}\delta v+q_{3} \delta w)d\alpha_{1}d\alpha_{2}, \tag{32}\]  
where \(q_{1}\), \(q_{2}\), and \(q_{3}\) are the external forces along the \(\alpha_{1}\), \(\alpha_{2}\), and \(\zeta\) directions, respectively.