*2.1. Basic Equations and Smagorinsky Subgrid Scale (SGS) Model*

The Favre-averaged filter and spatial filter are expressed as ¯ () and ˆ () respectively, and are used for the continuity and Navier–Stokes equations, which include the continuity equation, momentum equations, and equation of state. In Equation (2), the eddy viscosity assumption is used for SGS turbulence stress. A non-dimensionalization is applied to all variables by means of the streamwise velocity *U*<sup>0</sup> and chord length *C*. Since the boundary-fitted-grid is employed in our computation, general curvilinear coordinates have to be applied and is represented as *ξ*, *η*, and *ς*. Thus,

$$\frac{\partial \hat{\rho}}{\partial t} + \frac{1}{\int \hat{\partial} \tilde{\xi}^k} \left( J \hat{\rho} \hat{\mathcal{U}}^k \right) = 0,\tag{1}$$

$$\frac{\partial}{\partial t} \frac{(\not p \vec{u}\_i)}{\partial t} + \frac{1}{I} \frac{\partial}{\partial \not\_{\eth}^{\chi k}} \left( I \not p \vec{u}\_i \Omega^k \right) = -\frac{1}{I} \frac{\partial}{\partial \not\_{\eth}^{\chi k}} \left[ I \frac{\partial \not p}{\partial x\_i} \left( \not p + \frac{2}{3} \not p k\_{\text{sys}} \right) \right] + \frac{1}{I} \frac{\partial}{\partial \not\_{\eth}^{\chi k}} \left( I \frac{\partial \not p}{\partial x\_i} \mathcal{C}\_{ij} \right), \tag{2}$$

$$
\hat{p} = \hat{p}RT,\tag{3}
$$

where

$$
\sigma\_{\rm ij} = 2 \left( \frac{1}{R\varepsilon} + \hat{\rho} \nu\_{\rm sys} \right) \left( S\_{\rm ij} - \frac{1}{3} \delta\_{\rm ij} S\_{\rm kk} \right) \,\tag{4}
$$

where *J* represents the Jacobian of the coordinate transformation, *ρ* denotes the density, *U*ˆ *<sup>k</sup>* means the contravariant velocity, *ksgs* is the SGS kinetic energy, *δij* is the Kronecker symbol, T is the absolute temperature, R is the ideal gas constant, *Re* means the Reynolds number, as in *ρ*0*U*0/*ν*, *S*¯ *ij* denotes the strain rate tensor

$$\bar{S}\_{ij} = \frac{1}{2} \left( \frac{\partial \xi^k}{\partial x\_i} \frac{\partial \eta\_j}{\partial \xi^k} + \frac{\partial \xi^k}{\partial x\_j} \frac{\partial \eta\_i}{\partial \xi^k} \right) \tag{5}$$

and *νsgs* means the eddy viscosity since the effect of SGS turbulence on the large-scale motion is expected to be estimated by an SGS model.

Using the local equilibrium assumption, i.e., that the dissipation rate of SGS energy is in balance with the production rate, the dynamic Smagorinsky model can be obtained for LES:

$$\nu\_{\text{\ $}\$ } = \mathbb{C}\hat{\Lambda}^2 |\overline{\mathbb{S}}| \,\tag{6}$$

in which <sup>|</sup>*S*<sup>|</sup> denotes the norm of the strain rate tensor and is defined as 2*SijSij*; *C* is a variable to be dynamically determined from the grid-scale velocity field *u*¯.

Applying the test filter on the grid-filtered N-S equations, the Germano identity can be defined as

$$L\_{i\bar{j}} = T\_{i\bar{j}} - \overleftarrow{\tau}\_{i\bar{j}} = \overleftarrow{\Pi\_i}\overleftarrow{\tilde{\Pi}\_{\bar{j}}} - \overleftarrow{\tilde{\Pi}\_i}\overleftarrow{\tilde{\Pi}\_{\bar{j}}},\tag{7}$$

where *Lij* can be calculated based on the resolved scales, and *Tij* = *u iuj* − *uiu<sup>j</sup>* represents the residual turbulent stress at a test-filter scale Δ and can be given as

$$T\_{\rm ij} - \frac{1}{3} \delta\_{\rm ij} T\_{\rm kk} = -2C\tilde{\Delta}^2 |\tilde{\mathcal{S}}| \tilde{\mathcal{S}}\_{\rm ij}. \tag{8}$$

Upon substituting Equations (4) and (8) into Equation (7) and assuming Δˆ and C are constant inside the test filter, an equation for determining C was obtained:

$$L\_{ij} - \frac{1}{3} \delta\_{ij} L\_{kk} = -2\mathbb{C}\hat{\Lambda}^2 M\_{ij\prime} \tag{9}$$

where

$$\mathcal{M}\_{ij} = \frac{\widetilde{\Delta}^2}{\widetilde{\Delta}^2} |\widetilde{\mathfrak{S}}| \widetilde{\mathfrak{S}}\_{ij} - |\widetilde{\mathfrak{S}}| \dot{\mathfrak{S}}\_{ij}. \tag{10}$$

Minimization of the error of Equation (9) over all independent tensor components [29], and over some averaging region of statistical homogeneity, leads to

$$\mathcal{C} = -\frac{1}{2\dot{\Delta}^2} \frac{\langle L\_{ij}M\_{ij}\rangle}{\langle M\_{ij}M\_{ij}\rangle}.\tag{11}$$

Finally, substituting the dynamically determined variable *C* of Equation (11) into Equation (6), the eddy viscosity *νsgs* can be computed.

Through using the CIP-combined unified numerical procedure method [30], the fractional step method is applied for time-marching of Equation (2):

$$\left(\left(\hat{\rho}\mathfrak{u}\right)^{\mathrm{F}} = \frac{\Delta t}{2} \left(\frac{2}{\Delta t}(\hat{\rho}\mathfrak{u})^{\mathrm{n}} + 3\nabla \cdot \left[-\left(\hat{\rho}\mathfrak{u}\mathfrak{u}\right) + \tau\right]^{\mathrm{n}} - \nabla \cdot \left[-\left(\hat{\rho}\mathfrak{u}\mathfrak{u}\right) + \tau\right]^{\mathrm{n}-1}\right),\tag{12}$$

$$(\left(\not p\vec{u}\right)^{n+1} = \left(\not p\vec{u}\right)^{F} - \Delta t \nabla \not p^{n+1},\tag{13}$$

in which *τ* represents the viscous stress, Δ*t* the time increment, and *n* the time step count. The advancement method of time for Equation (1) can be written as:

$$
\boldsymbol{\phi}^{n+1} = \boldsymbol{\phi}^{n} - \Delta t \nabla \cdot (\boldsymbol{\phi} \boldsymbol{u})^{n+1}.\tag{14}
$$

Applying the divergence for Equation (13) and then substitution of this divergence into Equation (14) leads to the pressure elliptic equation for *p*ˆ*n*<sup>+</sup>1.

$$-\not{\rho}^{n+1} + \not{\rho}^n - \Delta t \nabla \cdot (\not{\rho} \vec{u})^F = -\Delta t^2 \nabla^2 \not{\rho}^{n+1}.\tag{15}$$

In order to take the compressibility into account for the pressure equation, Equation (3) can be written as the following for weakly compressible flows:

$$
\mathfrak{p}^{n+1} - \mathfrak{p}^n = (\mathfrak{p}^{n+1} - \mathfrak{p}^n)RT. \tag{16}
$$

Substituting Equation (16) into Equation (15) results in a pressure equation that accounts for the compressibility; thus,

$$\frac{1}{2}(\Delta t)^2 RT\nabla^2 \mathfrak{p}^{n+1} - \mathfrak{p}^{n+1} = (\Delta t)RT\nabla \cdot (\mathfrak{p}\mathfrak{u})^F - \mathfrak{p}^n. \tag{17}$$

Therefore, the second term in two sides of Equation (17) accounts for the compressibility effect. In a word, using Equation (17), *p*ˆ*n*+<sup>1</sup> is solved; and then, substituting *p*ˆ*n*+<sup>1</sup> into Equation (12), (*ρ*ˆ*u*¯)*n*+<sup>1</sup> is calculated; and finally, applying (1), *ρ*ˆ*n*+<sup>1</sup> is obtained. Note that he applicability of this procedure is limited with weakly compressible flows due to the fact that Δ*p*ˆ and Δ*ρ*ˆ change very slightly in one time step for the low Mach number turbulent flow.
