## 2 Methodology
### Simulation Method  
The Reynolds-Averaged Navier-Stokes (RANS) method is used to analyze the fluid flow over the submarine in this study. Following this method, the instantaneous velocity is decomposed into its time-averaged velocity (\(\overline{u}\)) and fluctuating velocity, in which the fluctuating velocity is modeled as a function of the turbulence eddy viscosity and turbulence kinetic energy. The governing equations of the RANS method for incompressible flow have been introduced by [PERSON] (2011) and [PERSON] (2021), as follows:  
\[\frac{\partial\overline{u}_{i}}{\partial x_{i}}=0 \tag{14}\]  
\[\frac{\partial\overline{u}_{i}}{\partial t}+\overline{u}_{j}\frac{\partial \overline{u}_{i}}{\partial x_{j}}=-\frac{\partial\overline{p}}{\partial x_{i }}+\upsilon\frac{\partial^{2}\overline{u}_{i}}{\partial x_{j}\partial x_{j}} -\frac{\partial\tau_{ij}}{\partial x_{j}} \tag{15}\]  
The Reynolds stress, \(\tau_{ij}=\overline{u_{i}^{\prime}u_{j}^{\prime}}\), is modeled as the function of the eddy viscosity and kinetic energy based on the Boussinesq hypothesis. The turbulent model, \(k\)-\(w\) SST (Shear Stress Transport) (2011), for solving turbulent flow around the submarine is presented as follows:  
\[\frac{\partial k}{\partial t}+\overline{u}_{i}\frac{\partial k}{\partial x_{i }}=P_{k}-\beta^{*}\omega k+\frac{\partial}{\partial x_{j}}\bigg{(}(\upsilon+ \sigma_{k}\upsilon_{t})\frac{\partial k}{\partial x_{i}}\bigg{)} \tag{16}\]  
\[\frac{\partial\omega}{\partial t}+\overline{u}_{i}\frac{\partial\omega}{ \partial x_{i}}=\frac{\gamma}{\mu_{t}}P_{k}-\beta\omega^{2}+\frac{\partial}{ \partial x_{j}}\bigg{[}(\upsilon+\sigma_{\omega}\upsilon_{t})\frac{\partial \omega}{\partial x_{i}}\bigg{]}+2(1-F_{1})\frac{\sigma_{\omega^{2}}}{\omega} \frac{\partial k}{\partial x_{i}}\frac{\partial\omega}{\partial x_{i}} \tag{17}\]  
\[\mu_{t}=\frac{\rho a_{1}k}{max(a_{1}\omega,\Omega F_{2})} \tag{18}\]  
where \(P_{k}\) is the production rate of turbulent kinetic energy by the mean flow. The variables are defined as follows.  
\[P_{k}=\tau_{ij}\frac{\partial u_{i}}{\partial x_{j}}\]  
\[\tau_{ij}=\mu_{t}\bigg{(}2S_{ij}-\frac{2}{3}\frac{\partial u_{k}}{\partial x_ {k}}\delta_{ij}\bigg{)}-\frac{2}{3}\rho k\delta_{ij}\]  
\[S_{ij}=\frac{1}{2}\bigg{(}\frac{\partial u_{i}}{\partial x_{j}}+\frac{ \partial u_{j}}{\partial x_{i}}\bigg{)}\]  
\[\beta_{i}=F_{1}\beta_{i,1}+(1-F_{1})\beta_{i,2}\]  
\[\sigma_{k}=\frac{1}{F_{1}/\sigma_{k,1}+(1-F_{1})\sigma_{k,2}}\]  
\[\sigma_{\omega}=\frac{1}{F_{1}/\sigma_{\omega,1}+(1-F_{1})\sigma_{\omega,2}}\]  
\[F_{1}=\tanh\Big{(}\arg_{1}^{4}\Big{)}\]  
\[F_{2}=\tanh\Big{(}\arg_{2}^{2}\Big{)}\]  
\[arg_{1}=min\bigg{[}max\bigg{(}\frac{\sqrt{k}}{\beta^{*}\omega y};\frac{500 \upsilon}{y^{2}\omega}\bigg{)};\frac{4\rho\sigma_{\omega z}k}{CD_{k\omega}y ^{2}}\bigg{]}\]\[arg_{2}=max\left(2\frac{\sqrt{k}}{0.09\omega y};\frac{500\
u}{y^{2}\omega}\right)\]  
\[CD_{k\omega}=max\left(2\rho\sigma_{\omega,2}\frac{1}{\omega}\frac{\partial k}{ \partial x_{i}}\frac{\partial\omega}{\partial x_{i}},10^{-20}\right)\]  
with the model constants \(\sigma_{k,1}\) = 1.176, \(\sigma_{\omega,1}\) = 2.0, \(\sigma_{k,2}\) = 1.0, \(\sigma_{\omega,2}\) = 1.168, \(a_{1}\) = 0.31, \(\beta_{i,1}\) = 0.075, and \(\beta_{i,2}\) = 0.0828 [23].