*2.4. Mathematical Model*

The governing equations are the conservation of mass, momentum (Navier–Stokes), and energy equations, which are described as follows in that order:

$$\frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_i} (\rho u\_i) = 0 \tag{2}$$

$$\frac{\partial}{\partial t}(\rho u\_i) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho u\_i u\_j) = -\frac{\partial P}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} - \frac{2}{3} \delta\_{ij} \frac{\partial u\_k}{\partial \mathbf{x}\_k} \right) \right] + \frac{\partial}{\partial \mathbf{x}\_j} \left( -\rho \left| u\_i' u\_j' \right| \right) \tag{3}$$

$$\frac{\partial}{\partial t}(\rho E) + \frac{\partial}{\partial \mathbf{x}\_j} \{u\_j(\rho E + P)\} = \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mathbf{k}\_{eff}) \frac{\partial T}{\partial \mathbf{x}\_j} \right] + \frac{\partial}{\partial \mathbf{x}\_j} \left[ u\_i \mu\_{eff} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} - \frac{2}{3} \delta\_{ij} \frac{\partial u\_k}{\partial \mathbf{x}\_k} \right) \right] \tag{4}$$
 
$$i\_r \ \mathbf{j}\_r \ k = 1, \ \mathbf{2}, \ \mathbf{3}$$

Here, ρ is the fluid density, *u* is the fluid velocity, *P* is the fluid pressure, and μ is the fluid viscosity; *E* is the specific internal energy, *keff* is the e ffective thermal conductivity, and μ*eff* is the e ffective dynamic viscosity.

To solve the governing equations, a density-based solver built using the finite volume method was employed. The implicit Roe's flux di fference scheme was used as the spatial discretization scheme. The least-squares cell-based method was selected for gradient calculation. Furthermore, the flow was discretized using the second-order upwind scheme. The turbulent kinetic energy and specific dissipation rate were set according to the first-order upwind scheme. Although the first-order upwind scheme may not yield results of greater accuracy compared with the second-order upwind scheme, the first-order upwind scheme easily converges and minimizes the computational cost.

For the time integration, the first-order implicit scheme was used to process the unsteady simulation. The unsteady constant-time stepping with a time step size of 8 × 10−<sup>5</sup> was adopted, and the maximum number of iterations per time step was selected as 20 to reach a residual convergence of 10−4. Additionally, the Courant number was set at 2 in all the simulations.

A suitable turbulence model can enhance the precision and reliability of the CFD simulation. Direct numerical simulation (DNS) and large eddy simulation (LES) can e ffectively analyze the complex fluid-flow model, especially for the turbulence structure [16–20]. However, these turbulence models require a high computational cost. The Reynolds-averaged Navier–Stokes (RANS) model can balance the computational cost and performance e fficiency of the numerical analysis. A few RANS models such as *k* − ε, *k* − ω, SST γ, transition SST, et cetera. are widely used for compressible flow simulations. In this study, the pod operates in the transonic to supersonic speed regimes. Therefore, to specify the turbulent flow, the SST *k* − ω model was applied. This is a hybrid model combining the advantages of the Wilcox (standard) *k* − ω and the *k* − ε models. The two variables, the turbulence kinetic energy *k* - m2/s2 and the specific dissipation rate ω - s<sup>−</sup><sup>1</sup> , were respectively determined by the following two equations [15,21,22]:

$$\frac{\partial(\rho k)}{\partial t} + \frac{\partial(\rho u\_j k)}{\partial \mathbf{x}\_j} = P\mathbf{x} - \boldsymbol{\beta}^\* \rho \boldsymbol{\omega} k + \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mu + \frac{\mu\_t}{\sigma\_k}) \frac{\partial k}{\partial \mathbf{x}\_j} \right] \tag{5}$$

$$\frac{\partial(\rho\omega)}{\partial t} + \frac{\partial(\rho u\_{\circ}\omega)}{\partial \mathbf{x}\_{\circ}} = \frac{\gamma}{\nu\_{\text{I}}}P\_{\text{K}} - \beta\rho\omega^{2} + \frac{\partial}{\partial \mathbf{x}\_{\circ}} \left[ \left( \mu + \frac{\mu\_{\text{I}}}{\sigma\_{\text{\omega}}} \right) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\circ}} \right] + 2\rho(1 - F\_{\text{I}}) \frac{1}{\sigma\_{\text{\omega},2}} \frac{1}{\omega} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\circ}} \frac{\partial \omega}{\partial \mathbf{x}\_{\circ}} \tag{6}$$

where

$$P\_K = \pi\_{i\bar{j}} \frac{\partial u\_{\bar{i}}}{\partial x\_{\bar{j}}} \, \tag{7}$$

$$
\tau\_{ij} = 2\mu\_l S\_{ij} - \frac{2}{3}\rho k \delta\_{ij}.\tag{8}
$$

where <sup>τ</sup>*ij* denotes the Reynolds stresses (kgm−1s−2), *Sij* denotes the mean rate of deformation component (s−1), and δ*ij* denotes the Kronecker delta function.

$$F\_1 = \tanh\left\{ \left( \min \left[ \max \left( \frac{\sqrt{k}}{\beta^\* \omega y}, \frac{500\nu}{y^2 \omega} \right), \frac{4\rho \sigma\_{\omega,2} k}{\overline{\text{CD}}\_{k\omega} y^2} \right] \right)^4 \right\},\tag{9}$$

$$\text{CD}\_{k\omega} = \max\left(2\rho \frac{1}{\sigma\_{\omega,2}} \frac{1}{\omega} \frac{\partial k}{\partial \mathbf{x}\_i} \frac{\partial \omega}{\partial \mathbf{x}\_j}, 10^{-10}\right) \tag{10}$$

$$
\sigma\_k = \frac{1}{F\_1/\sigma\_{k,1} + (1 - F\_1)/\sigma\_{k,2}} \, ^\prime \tag{11}
$$

$$
\sigma\_{\omega} = \frac{1}{F\_1/\sigma\_{\omega,1} + (1 - F\_1)/\sigma\_{\omega,2}}.\tag{12}
$$

σ*k* and σω are the turbulent Prandtl numbers for *k* and ω, respectively. μ*t* (kg/ms) is the turbulence viscosity, which is calculated by the following equations:

$$\frac{1}{\rho}\mu\_t = \frac{a\_1 k}{\max(a\_1 \omega, SF\_2)},\tag{13}$$

$$F\_2 = \tanh\left(\max\left(2\frac{\sqrt{k}}{\beta^\*\omega y}, \frac{500\nu}{y^2\omega}\right)^2\right). \tag{14}$$

In Equation (13), the term *S* = -2*SijSij*1/2 is the invariant measure of the strain rate, and *a*1 is a constant equal to 0.31. The other constant values in the above equations are given as follows [15]:

$$
\beta^\* = 0.09,
$$

$$
\sigma\_{k,1} = 2.0, \,\sigma\_{k,2} = 1.0,
$$

$$
\sigma\_{\omega,1} = 2.0, \,\sigma\_{\omega,2} = 1.168.
$$

This turbulence model has been widely used to simulate the aerodynamic characteristics of high-speed trains [6,8,10,23–25]; it has improved the accuracy and reliability of free shear flows and transonic shockwaves predictions [15]. For more details of this model, readers are directed to [6,21,22].
