A quasi-dynamic one-equation model with joint constraints of kinetic energy and helicity fluxes for large eddy simulation of rotating turbulence

For the incompressible turbulent flows, the filtered Navier-Stokes equations are taken as follows: {subeqnarray} ∂¯ui∂xi &= 0 ,
∂¯ui∂t + ¯u_j∂¯ui∂xj = -1ρ∂¯p∂xi + ν∂2¯ui∂xj2 + ¯f_i -∂τij∂xj. where $\left(\bar{\cdot}\right)$ denotes spatial filtering at scale $\Delta$, $\bar{f}_{i}$ is the filtered forcing, and $\tau_{ij}=\overline{u_{i}u_{j}}-\bar{u}_{i}\bar{u}_{j}$ is the SGS stress tensor that needs to be modelled. The isotropic part of the SGS stress is often absorbed into pressure and eddy viscosity based SGS models based on the Boussinesq hypothesis take the form

$$\tau_{ij}-\frac{1}{3}\tau_{kk}\delta_{ij}=-2\nu_{sgs}\bar{S}_{ij},\quad\bar{S}_{ij}=\frac{1}{2}\left(\frac{\partial\bar{u}_{i}}{\partial x_{j}}+\frac{\partial\bar{u}_{j}}{\partial x_{i}}\right),$$

(1)

where $\delta_{ij}$ is the Kronecker-delta function.

In the absence of system rotation, the filtered forcing term in momentum equation (2.1) contains the Coriolis force term $-2\boldsymbol{\varOmega}\times\boldsymbol{\bar{u}}$, where $\bar{\boldsymbol{u}}$ is the filtered velocity and $\boldsymbol{\varOmega}$ is the rotating vector of the system. The pressure $\bar{p}$ is the total pressure containing both the static pressure and the centrifugal force.

According to Yan et al. (2020), the first and second channel of helicity flux $\varPi_{\Delta}^{H1},\varPi_{\Delta}^{H2}$ are expressed as

$$\varPi_{\Delta}^{H1}=-\tau_{ij}\bar{R}_{ij},\quad\varPi_{\Delta}^{H2}=-\gamma_{ij}\bar{\varOmega}_{ij}.$$

(2)

Here, $\bar{\varOmega}_{ij}=\frac{1}{2}\left(\partial\bar{u}_{i}/\partial x_{j}-\partial\bar{u}_{j}/\partial x_{i}\right)$, $\bar{R}_{ij}=\frac{1}{2}\left(\partial\bar{\omega}_{i}/\partial x_{j}+\partial\bar{\omega}_{j}/\partial x_{i}\right)$ is the large-scale rotation rate and the large-scale symmetric vorticity gradient, respectively. $\gamma_{ij}=\left(\overline{\omega_{i}u_{j}}-\bar{\omega}_{i}\bar{u}_{j}\right)-\left(\overline{\omega_{j}u_{i}}-\bar{\omega}_{j}\bar{u}_{i}\right)$ is the subgrid-scale vortex stretching stress.

The kinetic energy flux which represents energy transferred to small scales, is defined as

$$\varPi_{\Delta}^{E}=\tau_{ij}\bar{S}_{ij},$$

(3)

where $\bar{S}_{ij}=\frac{1}{2}\left(\partial\bar{u}_{i}/\partial x_{j}+\partial\bar{u}_{j}/\partial x_{i}\right)$ is the large-scale strain rate tensor.

Following the method of quasi-dynamic procedure, we derive SGS kinetic energy equation for both incompressible and compressible flows. For incompressible flow, the evolution equation of filtered sub-grid kinetic energy $k_{sgs}=\overline{u_{k}^{\prime}u_{k}^{\prime}}/2$ can be expressed as:

$$\begin{split}\frac{\partial k_{sgs}}{\partial t}=&-\frac{\partial}{\partial x_{j}}(k_{sgs}\bar{u}_{j})-\frac{1}{2}\frac{\partial}{\partial x_{j}}(\overline{u_{i}u_{i}u_{j}}-\overline{u_{i}u_{i}}\bar{u}_{j})-\frac{\partial}{\partial x_{j}}(\overline{pu_{j}}-\bar{p}\bar{u}_{j})\\ &+\frac{\partial}{\partial x_{j}}(\nu\frac{\partial k_{sgs}}{\partial x_{j}})+\frac{\partial}{\partial x_{j}}(\tau_{ij}\bar{u}_{i})\\ &-\nu(\overline{\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial u_{i}}{\partial x_{j}}}-\frac{\partial\bar{u}_{i}}{\partial x_{j}}\frac{\partial\bar{u}_{i}}{\partial x_{j}})-\tau_{ij}\frac{\partial\bar{u}_{i}}{\partial x_{j}}.\end{split}$$

(4)

To further simplify this equation, we use the analogy proposed by Knight et al. (1998):

$$\frac{1}{2}\left(\overline{u_{i}u_{i}u_{j}}-\overline{u_{i}u_{i}}\bar{u}_{j}\right)=\frac{1}{2}\left(\overline{u_{i}u_{i}u_{k}}-\bar{u}_{i}\bar{u}_{i}\bar{u}_{k}-2k_{sgs}\bar{u}_{k}\right)\approx\tau_{ik}\bar{u}_{i},$$

(5)

and we noticed that with the continuity equation, the term representing diffusion by pressure effect can be taken as

$$\frac{\partial}{\partial x_{j}}(\overline{pu_{j}}-\bar{p}\bar{u}_{j})=\overline{u_{j}\frac{\partial p}{\partial x_{j}}}-\bar{u}_{j}\frac{\partial\bar{p}}{\partial x_{j}}.$$

(6)

Finally, the simplified SGS kinetic energy equation can be written as:

$$\begin{split}\frac{\partial k_{sgs}}{\partial t}=&-\frac{\partial}{\partial x_{j}}(k_{sgs}\bar{u}_{j})-\left(\overline{u_{j}\frac{\partial p}{\partial x_{j}}}-\bar{u}_{j}\frac{\partial\bar{p}}{\partial x_{j}}\right)\\ &+\frac{\partial}{\partial x_{j}}(\nu\frac{\partial k_{sgs}}{\partial x_{j}})-\nu(\overline{\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial u_{i}}{\partial x_{j}}}-\frac{\partial\bar{u}_{i}}{\partial x_{j}}\frac{\partial\bar{u}_{i}}{\partial x_{j}})-\tau_{ij}\frac{\partial\bar{u}_{i}}{\partial x_{j}}.\end{split}$$

(7)

In the simulation of incompressible channel flow, which is our test case, the non-dimensional equation can be obtained by using half channel width $h$ and wall shear velocity $u_{\tau}$:

$$\begin{split}\frac{\partial k_{sgs}}{\partial t}=&-\frac{\partial}{\partial x_{j}}(k_{sgs}\bar{u}_{j})-\left(\overline{u_{j}\frac{\partial p}{\partial x_{j}}}-\bar{u}_{j}\frac{\partial\bar{p}}{\partial x_{j}}\right)\\ &+\frac{1}{Re_{\tau}}\frac{\partial^{2}k_{sgs}}{\partial x_{j}\partial x_{j}}-\frac{1}{Re_{\tau}}(\overline{\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial u_{i}}{\partial x_{j}}}-\frac{\partial\bar{u}_{i}}{\partial x_{j}}\frac{\partial\bar{u}_{i}}{\partial x_{j}})-\tau_{ij}\frac{\partial\bar{u}_{i}}{\partial x_{j}},\end{split}$$

(8)

where $Re_{\tau}=\frac{\rho u_{\tau}h}{\nu}$ is the wall shear Reynolds number and $\bar{u}=\bar{u}^{*}/u_{\tau}$ is actually the dimensionless velocity $u^{+}$.

In equation (8), there are three quantities that need to be modelled, they are the SGS stress tensor $\tau_{ij}$, the pressure dilatation $\varPi_{p}$ and the viscous dissipation $\epsilon$, respectively:

$$\tau_{ij},\quad\varPi_{p}=\left(\overline{u_{j}\frac{\partial p}{\partial x_{j}}}-\bar{u}_{j}\frac{\partial\bar{p}}{\partial x_{j}}\right),\quad\varepsilon=(\overline{\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial u_{i}}{\partial x_{j}}}-\frac{\partial\bar{u}_{i}}{\partial x_{j}}\frac{\partial\bar{u}_{i}}{\partial x_{j}}).$$

$(9a,b,c)$

Precise and accurate modelling of these unclosed quantities is crucial to the success of the new model. In the upcoming section, we will present modelling approaches for these unclosed quantities.

In this section, we will start by introducing the quasi-dynamic process and helicity flux dual channel constraint for incompressible flows and then extend the result to compressible flows.

A series expansion of filtered quantity is given by Bedford & Yeo (1993) as:

	$\displaystyle\overline{fg}-\bar{f}\bar{g}=$	$\displaystyle 2\alpha\frac{\partial\bar{f}}{\partial x_{k}}\frac{\partial\bar{g}}{\partial x_{k}}+\frac{1}{2!}(2\alpha)^{2}\frac{\partial^{2}\bar{f}}{\partial x_{k}\partial x_{l}}\frac{\partial^{2}\bar{g}}{\partial x_{k}\partial x_{l}}$		(10)
		$\displaystyle+\frac{1}{3!}(2\alpha)^{2}\frac{\partial^{3}\bar{f}}{\partial x_{k}\partial x_{l}\partial x_{m}}\frac{\partial^{2}\bar{g}}{\partial x_{k}\partial x_{l}\partial x_{m}}+\cdots,$

where

$$\alpha(y)=\int_{-\infty}^{\infty}x^{2}G(x,y)\mathrm{d}x.$$

(11)

For a box filter, an exact value $\alpha=\Delta^{2}/24$ can be obtained. In practice, this is set to be $2\alpha=C_{0}\Delta_{k}^{2}$, where $\Delta_{k}$ is the grid length scale in $k$ direction. Thus, the unclosed terms can be modelled, viz.

$$\tau_{ij}=C_{0}\Delta_{k}^{2}\frac{\partial\bar{u}_{i}}{\partial x_{k}}\frac{\partial\bar{u}_{j}}{\partial x_{k}}+\frac{1}{2!}(C_{0}^{2}\Delta_{k}^{2}\Delta_{l}^{2})\frac{\partial^{2}\bar{u}_{i}}{\partial x_{k}\partial x_{l}}\frac{\partial^{2}\bar{u}_{j}}{\partial x_{k}\partial x_{l}}+\cdots,$$

(12)

As it has been proved by previous successful applications of QKM model (Qi et al., 2022), it is enough for us to reserve only the first term for modelling $\tau_{ij}:$

$$\tau_{ij}\approx C_{0}\Delta_{k}^{2}\frac{\partial\bar{u}_{i}}{\partial x_{k}}\frac{\partial\bar{u}_{j}}{\partial x_{k}}.$$

(13)

As the SGS kinetic energy $k_{sgs}=\frac{1}{2}\tau_{kk}$, one can obtain

$$C_{0}=\frac{2k_{sgs}}{\Delta_{l}^{2}\dfrac{\partial\bar{u}_{k}}{\partial x_{l}}\dfrac{\partial\bar{u}_{k}}{\partial x_{l}}}.$$

(14)

For the unclosed quantities $\varPi_{p}$ and $\varepsilon$, we also apply this method to get

	$\displaystyle\varPi_{p}$	$\displaystyle\approx$	$\displaystyle C_{0}\Delta_{l}^{2}\frac{\partial\bar{p}}{\partial x_{l}}\frac{\partial\bar{u}_{k}}{\partial x_{l}\partial x_{k}},$		(15)
	$\displaystyle\varepsilon$	$\displaystyle\approx$	$\displaystyle C_{0}\Delta_{k}^{2}\frac{\partial^{2}\bar{u}_{i}}{\partial x_{j}\partial x_{k}}\frac{\partial^{2}\bar{u}_{i}}{\partial x_{j}\partial x_{k}}.$		(16)

Now eqaution (8) is closed and the value of $C_{0}$ can be solved dynamically. Through this so-called quasi-dynamic process, a more precise distribution of $k_{sgs}$ can be obtained. As mentioned, the deconvolution form of equation (13) shows a high correlation with the real value but is often numerically unstable. To maintain a strong correlation with the real values while ensuring numerical stability, we use the equation for SGS kinetic energy to derive precise SGS kinetic energy and determine the appropriate coefficients for the deconvolution model. We employ a joint-constraint approach to restrict the Smagorinsky model, which exhibits excellent numerical stability.

To tighten the constraints of the Smagorinsky model, the methodology leverages the kinetic energy flux and helicity flux. The primary objective is to ensure that the modified model demonstrates a strong correlation with both fluxes while maintaining numerical stability. Thus, the proposed approach employs a least square constraint method to determine the appropriate coefficient for the Smagorinsky model.

We denote the quantity related to the modified Smagorinksy model with a dynamic constant $C_{sm}$ with superscript $SM$ and quantities related to the deconvolution model with a dynamic coefficient $C_{0}$ that is solved through SGS kinetic energy equation with superscript $EQ$. The kinetic energy flux and dual-channel helicity flux can be expressed as

	$\displaystyle\varPi_{E}^{SM}=C_{sm}\Delta^{2}|\bar{S}|\bar{S}_{ij}\frac{\partial\bar{u}_{i}}{\partial x_{j}},\quad\varPi_{E}^{EQ}$	$\displaystyle=$	$\displaystyle C_{0}\Delta_{k}^{2}\frac{\partial\bar{u}_{i}}{\partial x_{k}}\frac{\partial\bar{u}_{j}}{\partial x_{k}}\frac{\partial\bar{u}_{i}}{\partial x_{j}},$		(17)
	$\displaystyle\varPi_{H1}^{SM}=-C_{sm}\Delta^{2}|\bar{S}|\bar{S}_{ij}\frac{\partial\bar{\omega}_{i}}{\partial x_{j}},\quad\varPi_{H1}^{EQ}$	$\displaystyle=$	$\displaystyle-C_{0}\Delta_{k}^{2}\frac{\partial\bar{u}_{i}}{\partial x_{k}}\frac{\partial\bar{u}_{j}}{\partial x_{k}}\frac{\partial\bar{\omega}_{i}}{\partial x_{j}},$		(18)
	$\displaystyle\varPi_{H2}^{SM}=-\nabla\times(C_{sm}\Delta^{2}|\bar{S}|\bar{S}_{ij})\frac{\partial\bar{u}_{i}}{\partial x_{j}},\quad\varPi_{H2}^{EQ}$	$\displaystyle=$	$\displaystyle-C_{0}\Delta_{k}^{2}\left(\frac{\partial\bar{u}_{j}}{\partial x_{k}}\frac{\partial\bar{\omega}_{i}}{\partial x_{k}}-\frac{\partial\bar{u}_{i}}{\partial x_{k}}\frac{\partial\bar{\omega}_{j}}{\partial x_{k}}\right)\frac{\partial\bar{u}_{i}}{\partial x_{j}}.$		(19)

To get the constraint for $C_{sm}$, the variation of the model coefficient is frozen. We can make the assumption that $\nabla C_{sm}$ is negligible so that:

$$\nabla\times(C_{sm}\Delta^{2}|\bar{S}|\bar{S}_{ij})\frac{\partial\bar{u}_{i}}{\partial x_{j}}=C_{sm}\nabla\times(\Delta^{2}|\bar{S}|\bar{S}_{ij})\frac{\partial\bar{u}_{i}}{\partial x_{j}}.$$

(20)

Now the deviation between the Smagorinsky model and the dynamic deconvolution model can be expressed as

	$\displaystyle\delta_{E}(C_{sm})$	$\displaystyle=\left(\varPi_{E}^{EQ}-\varPi_{E}^{SM}\right)^{2},$		(21)
	$\displaystyle\delta_{H1}(C_{sm})$	$\displaystyle=\left(\varPi_{H1}^{EQ}-\varPi_{H1}^{SM}\right)^{2},$
	$\displaystyle\delta_{H2}(C_{sm})$	$\displaystyle=\left(\varPi_{H2}^{EQ}-\varPi_{H2}^{SM}\right)^{2}.$

Since the total derivation $\delta=\delta_{E}+\delta_{H1}+\delta_{H2}$ can achieve minimum in a reasonable giving range of $C_{sm}$, take $\partial\delta/{\partial C_{sm}}=0$, we got the optimized $C_{sm}$:

	$\displaystyle C_{sm}$	$\displaystyle=$	$\displaystyle\frac{\varPi_{E}^{EQ}B_{1}+\varPi_{H1}^{EQ}B_{2}+\varPi_{H2}^{EQ}B_{3}}{B_{1}^{2}+B_{2}^{2}+B_{3}^{2}},$		(22)
	$\displaystyle B_{1}$	$\displaystyle=$	$\displaystyle-\Delta^{2}|\bar{S}|\bar{S}_{ij}\frac{\partial\bar{u}_{i}}{\partial x_{j}},$		(23)
	$\displaystyle B_{2}$	$\displaystyle=$	$\displaystyle-\Delta^{2}|\bar{S}|\bar{S}_{ij}\frac{\partial\bar{\omega}_{i}}{\partial x_{j}},$		(24)
	$\displaystyle B_{3}$	$\displaystyle=$	$\displaystyle-\nabla\times(\Delta^{2}|\bar{S}|\bar{S}_{ij})\frac{\partial\bar{u}_{i}}{\partial x_{j}}.$		(25)

With equation (22), the eddy viscosity can be obtained, i.e.

$$\nu_{sgs}=C_{sm}\Delta^{2}|\tilde{S}|.$$

(26)

For the range of clipping for $C_{sm}$, the relation recommend by Garnier et al. (2009) is adopted, i.e.

	$\displaystyle\nu_{sgs}+\nu\geqslant 0,$		(27)
	$\displaystyle C_{sm}\leqslant C_{max},$		(28)

where $C_{max}$ can be taken as the square of Smagorinksy constant $0.18^{2}$.

As the first test case of the new model, incompressible streamwise-rotating channel flow is widely studied in literature(Yang & Wang, 2018; Oberlack et al., 2006; Alkishriwi et al., 2008). Figure 1 shows the schematic diagram for the streamwise-rotating channel flow.

The rotating vector of the system is defined as $\boldsymbol{\varOmega}=\left(\varOmega_{x},0,0\right)$, where $\varOmega_{x}$ is set to be a constant so that the rotation is homogeneous. The flow is driven by constant streamwise pressure gradient $\boldsymbol{\varPi}=\left(\varPi_{x},0,0\right)$ and sufficient computation time (approximately $100h/u_{\tau}$) is taken to ensure the time duration is long enough for achieving a statistically stationary state. Periodic boundary conditions are applied to the streamwise and spanwise of the computational domain and the no-slip and impermeable boundary conditions are applied to the top and bottom walls. For LES, the dynamic Smagorinsky model (DSM) and the wall-adapting local eddy-viscosity (WALE) model are compared with the new model. The test filter width of the DSM model is chosen to be $2\Delta$, where $\Delta=(\Delta_{x}\Delta_{y}\Delta_{z})^{1/3}$ is the length scale of the mesh. For the WALE model, the model constant is chosen to be $C_{w}=\sqrt{10.6}C_{s}$ following Ducros et al. (1998), where $C_{s}$ is the Smagorinksy model constant.

As it is vitally important to control the numerical dissipation in LES(Chapelier & Lodato, 2017), a high-resolution pseudo-spectral method is used to solve the control equation. The grid near the wall is refined by using Cheyshev-Gauss-Lobatto points and flow variables are expanded into Chebyshev polynomials in the wall norm direction. Fourier series are used to expand flow variables in the streamwise and spanwise direction on uniform grids. More computational details can be found in the literature. To examine the effectiveness of the new model, a priori tests were carried out for channel flow with Reynolds number ($Re_{\tau}=u_{\tau}h/\nu$) fixed at 180 and the rotation number ($Ro_{\tau}=2\varOmega h/u_{\tau}$) fixed at 7.5. Then a posteriori analysis of the performance of the new model in higher Reynolds number ($Re_{\tau}=395$) is also presented. Table 1 shows the grid settings for direct numerical simulation (DNS) and LES in two different Reynolds numbers. The computational domain is set to $32\pi\times 2\times 8\pi$ in the streamwise, wall-normal, and spanwise direction, this computational domain is large enough to capture large-scale intermittency and most energetic eddy motions including Taylor-Götler vortices.

For the compressible flow simulation, applying a filter on the Navier-Stokes equation gives

	$\displaystyle\frac{\partial\bar{\rho}}{\partial t}+\frac{\partial\bar{\rho}\tilde{u}_{i}}{\partial x_{j}}$	$\displaystyle=$	$\displaystyle 0,$		(33)
	$\displaystyle\frac{\partial\bar{\rho}\tilde{u}_{i}}{\partial t}+\frac{\partial\bar{\rho}\tilde{u}_{i}\tilde{u}_{j}}{\partial x_{j}}$	$\displaystyle=$	$\displaystyle-\frac{\partial\bar{p}}{\partial x_{i}}+\frac{\partial\tilde{\sigma}_{ij}}{\partial x_{j}}-\frac{\partial\tau_{ij}}{\partial x_{j}},$		(34)
	$\displaystyle\frac{\partial\bar{\rho}\tilde{E}}{\partial t}+\frac{\partial(\bar{\rho}\tilde{E}+\bar{p})\tilde{u}_{j}}{\partial x_{j}}$	$\displaystyle=$	$\displaystyle-\frac{\partial\tilde{q}_{j}}{\partial x_{j}}+\frac{\partial\tilde{\sigma}_{ij}\tilde{u}_{i}}{\partial x_{j}}-\frac{\partial C_{p}Q_{j}}{\partial x_{j}}-\frac{\partial J_{j}}{\partial x_{j}},$		(35)

where $\left(\tilde{\cdot}\right)$ represents density-weighted (Favre) filtering ($\tilde{\phi}=\bar{\rho}\phi/\bar{\rho}$). Compare to the filtered equation of incompressible flow, there are several differences: The energy equation is needed here and the filtered total energy is written as

$$\bar{\rho}\tilde{E}=\bar{\rho}C_{v}\tilde{T}+\frac{1}{2}\bar{\rho}\tilde{u}_{i}\tilde{u}_{i}+\bar{\rho}k_{sgs},\quad\bar{\rho}k_{sgs}=\frac{1}{2}\bar{\rho}(\widetilde{u_{i}u_{i}}-\tilde{u}_{i}\tilde{u}_{i})$$

$(36a,b).$

The resolved heat flux $\tilde{q}_{j}$, the SGS heat flux $Q_{j}$ and the SGS turbulent diffusion term $J_{j}$ are expressed as

	$\displaystyle\tilde{q}_{j}$	$\displaystyle=$	$\displaystyle\frac{C_{p}\mu(\tilde{T})}{Pr}\frac{\partial\tilde{T}}{\partial x_{j}},$		(37)
	$\displaystyle Q_{j}$	$\displaystyle=$	$\displaystyle\bar{\rho}(\widetilde{u_{j}T}-\tilde{u}_{j}\tilde{T}),$		(38)
	$\displaystyle J_{j}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\bar{\rho}(\widetilde{u_{i}u_{i}u_{j}}-\widetilde{u_{i}u_{i}}\tilde{u}_{j}),$		(39)

where the molecular viscosity $\mu$ takes from Sutherland’s law.

The dual channels of helicity in compressible flows can be derived in a similar way as in incompressible flows. We can obtain the governing equation of large-scale helicity $\tilde{H}=\tilde{\boldsymbol{u}}\cdot\tilde{\boldsymbol{\omega}}$ as follows:

$$\frac{\partial\tilde{H}}{\partial t}+\frac{\partial\tilde{J}_{j}}{\partial x_{j}}=-\varPi_{\Delta}^{H1}-\varPi_{\Delta}^{H2}+\varPhi_{l}^{1}+\varPhi_{l}^{2}+D_{l}^{1}+D_{l}^{2},$$

(40)

where $\tilde{J}_{j}$ is the spatial transport term and $\varPhi_{l}$ is the pressure term, $D_{l}$ is the viscosity term. Here, superscript 1 is named the first channel originating from the Farve filtered velocity governing equation, and superscript 2 is named the second channel originating from the Farve filtered vorticity equation. Here, we mainly focus on the dual channels of helicity flux. The first and second channels of helicity flux are written respectively as

	$\displaystyle\varPi_{\Delta}^{H1}$	$\displaystyle=$	$\displaystyle-\bar{\rho}\tilde{\tau}_{ij}\frac{\partial}{\partial x_{j}}\left(\frac{\tilde{\omega}_{i}}{\bar{\rho}}\right),$		(41)
	$\displaystyle\varPi_{\Delta}^{H2}$	$\displaystyle=$	$\displaystyle-\epsilon_{ikl}\frac{\partial\bar{\rho}\tilde{\tau}_{lj}}{\partial x_{k}}\frac{\partial}{\partial x_{j}}\left(\frac{\tilde{u}_{i}}{\bar{\rho}}\right)-\bar{\rho}\tilde{\tau}_{ij}\frac{\partial}{\partial x_{j}}\left(\epsilon_{ikl}\tilde{u}_{k}\frac{\partial}{\partial x_{l}}\left(\frac{1}{\bar{\rho}}\right)\right).$		(42)

If we take the filtered density as constant, the equations above would recover to the form of incompressible flows. For the kinetic energy flux, it has the form same as in incompressible flows, i.e.

$$\varPi_{\Delta}^{E}=\tau_{ij}\tilde{S}_{ij}.$$

(43)

Similar to the SGS kinetic energy equation of incompressible flow, one can derive the SGS kinetic equation of compressible flow, written as

	$\displaystyle\frac{\partial\bar{\rho}k_{sgs}}{\partial t}=$	$\displaystyle-\frac{\partial\bar{\rho}k_{sgs}\tilde{u}_{j}}{\partial x_{j}}-\tau_{ij}\tilde{S}_{ij}-\varepsilon_{s}-\varepsilon_{d}+\varPi_{p}$		(44)
		$\displaystyle-\frac{\partial J_{j}}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}\left(\bar{\mu}\frac{\partial k_{sgs}}{\partial x_{j}}\right)+\frac{\partial}{\partial x_{j}}\left[\tau_{ij}\tilde{u}_{i}+\bar{\mu}\frac{\partial}{\partial x_{i}}\left(\frac{\tau_{ij}}{\bar{\rho}}\right)+R\tilde{q}_{j}\right],$

with the following notation:

	$\displaystyle\varepsilon_{s}$	$\displaystyle=$	$\displaystyle 2\bar{\mu}\left[\widetilde{\mathbb{S}_{ij}\mathbb{D}_{ij}}-\widetilde{\mathbb{S}}_{ij}\widetilde{\mathbb{D}}_{ij}\right],$		(45)
	$\displaystyle\varepsilon_{d}$	$\displaystyle=$	$\displaystyle\frac{\partial}{\partial x_{j}}\left[\frac{5}{3}\left(\bar{\mu}\widetilde{u_{j}\frac{\partial u_{k}}{\partial x_{k}}}-\bar{\mu}\tilde{u}_{j}\frac{\partial\tilde{u}_{k}}{\partial\tilde{x}_{k}}\right)\right],$		(46)
	$\displaystyle\varPi_{p}$	$\displaystyle=$	$\displaystyle\overline{p\frac{\partial u_{k}}{\partial x_{k}}}-\bar{p}\frac{\partial\tilde{u}_{k}}{\partial x_{k}}.$		(47)

Where $R=C_{p}-C_{v}$ is the specific gas constant and $\tilde{\mathbb{S}}_{ij},\tilde{\mathbb{D}}_{ij}$ are the traceless strain rate/velocity gradient tensor. Again, like the incompressible SGS kinetic energy equation, the equation of compressible SGS kinetic energy is not fully closed, several terms need to be modelled: the SGS stress tensor $\tau_{ij}$, the solenoidal dissipation $\varepsilon_{s}$, the dilatational dissipation $\varepsilon_{d}$ and the pressure dilatation $\varPi_{p}$. These terms will be modelled using the quasi-dynamic process method introduced in the next section.

As is mentioned in the derivation of the quasi-dynamic method of incompressible flow, by taking series expansion of the unclosed quantities, we can get the approximation of them, expressed as

	$\displaystyle\tau_{ij}$	$\displaystyle\approx$	$\displaystyle C_{0}\Delta_{k}^{2}\bar{\rho}\frac{\partial\tilde{u}_{i}}{\partial x_{k}}\frac{\partial\tilde{u}_{j}}{\partial x_{k}},$		(48)
	$\displaystyle\varepsilon_{s}$	$\displaystyle\approx$	$\displaystyle 2C_{0}\Delta_{k}^{2}\mu(\tilde{T})\frac{\partial\tilde{\mathbb{S}}_{ij}}{\partial x_{k}}\frac{\partial\tilde{\mathbb{D}}_{ij}}{\partial x_{k}},$		(49)
	$\displaystyle\varepsilon_{d}$	$\displaystyle\approx$	$\displaystyle\frac{5}{3}\frac{\partial}{\partial x_{j}}\left[C_{0}\Delta_{l}^{2}\mu(\tilde{T})\frac{\partial\tilde{u}_{j}}{\partial x_{l}}\frac{\partial^{2}\tilde{u}_{k}}{\partial x_{k}\partial x_{l}}\right],$		(50)
	$\displaystyle\varPi_{p}$	$\displaystyle\approx$	$\displaystyle C_{0}\Delta_{l}^{2}\frac{\partial\bar{p}}{\partial x_{l}}\frac{\partial^{2}\tilde{u}_{k}}{\partial x_{l}\partial x_{k}}.$		(51)

Here we note that all four unclosed quantities use the same $C_{0}$, this has been proved reasonable by both a prior and a posteriori tests (Qi et al., 2022).

From equation (48), we can find that

$$\tau_{kk}=C_{0}\Delta_{l}^{2}\bar{\rho}\frac{\partial\tilde{u}_{k}}{\partial x_{l}}\frac{\partial\tilde{u}_{k}}{\partial x_{l}},$$

(52)

and we have the relation $\tau_{kk}=2\bar{\rho}k_{sgs}$, so we can confirm the real value of $C_{0}$ as

$$C_{0}=\frac{2k_{sgs}}{\Delta^{2}_{l}\frac{\partial\tilde{u}_{k}}{\partial x_{l}}\frac{\partial\tilde{u}_{k}}{\partial x_{l}}}.$$

(53)

The kinetic energy flux $\varPi_{\Delta}^{E}=\tau_{ij}\tilde{S}_{ij}$, is another quantity used to constrain the model coefficient. Using eqaution (48), we can derive the deconvolution model representation of these quantities as follows:

	$\displaystyle\varPi_{\Delta}^{E}$	$\displaystyle=$	$\displaystyle C_{0}\Delta_{k}^{2}\bar{\rho}\frac{\partial\tilde{u}_{i}}{\partial x_{k}}\frac{\partial\tilde{u}_{j}}{\partial x_{k}}\tilde{S}_{ij},$		(54)
	$\displaystyle\varPi_{\Delta}^{H1}$	$\displaystyle=$	$\displaystyle-\bar{\rho}C_{0}\Delta_{k}^{2}\bar{\rho}\frac{\partial\tilde{u}_{i}}{\partial x_{k}}\frac{\partial\tilde{u}_{j}}{\partial x_{k}}\frac{\partial}{\partial x_{j}}\left(\frac{\tilde{\omega}_{i}}{\bar{\rho}}\right),$		(55)

	$\displaystyle\varPi_{\Delta}^{H2}=$	$\displaystyle-\epsilon_{ikl}\frac{\partial\bar{\rho}}{\partial x_{k}}\left(C_{0}\Delta_{m}^{2}\bar{\rho}\frac{\partial\tilde{u}_{l}}{\partial x_{m}}\frac{\partial\tilde{u}_{j}}{\partial x_{m}}\right)\frac{\partial}{\partial x_{j}}\left(\frac{\tilde{u}_{i}}{\bar{\rho}}\right)$		(56)
		$\displaystyle-\bar{\rho}C_{0}\Delta_{k}^{2}\bar{\rho}\frac{\partial\tilde{u}_{i}}{\partial x_{k}}\frac{\partial\tilde{u}_{j}}{\partial x_{k}}\frac{\partial}{\partial x_{j}}\left(\epsilon_{ikl}\tilde{u}_{k}\frac{\partial}{\partial x_{l}}\left(\frac{1}{\bar{\rho}}\right)\right).$

Again, adopting the least-square method, the constraint for the Smagorinksy model coefficient can be obtained:

	$\displaystyle C_{sm}$	$\displaystyle=$	$\displaystyle\frac{\varPi_{E}^{EQ}B_{1}+\varPi_{H1}^{EQ}B_{2}+\varPi_{H2}^{EQ}B_{3}}{B_{1}^{2}+B_{2}^{2}+B_{3}^{2}},$		(57)
	$\displaystyle B_{1}$	$\displaystyle=$	$\displaystyle-2\bar{\rho}\Delta^{2}|\tilde{S}|\tilde{\mathbb{S}}_{ij}\tilde{S}_{ij},$		(58)
	$\displaystyle B_{2}$	$\displaystyle=$	$\displaystyle-\bar{\rho}^{2}\Delta^{2}|\tilde{S}|\tilde{\mathbb{S}}_{ij}\frac{\partial\tilde{\omega}_{i}/\bar{\rho}}{\partial x_{j}},$		(59)
	$\displaystyle B_{3}$	$\displaystyle=$	$\displaystyle-\nabla\times(\bar{\rho}\Delta^{2}|\tilde{S}|\tilde{\mathbb{S}}_{ij})\frac{\partial\tilde{u}_{i}/\bar{\rho}}{\partial x_{j}}+\bar{\rho}^{2}\Delta^{2}|\tilde{S}|\tilde{\mathbb{S}}_{ij}\frac{\partial}{\partial x_{j}}\left(\epsilon_{ikl}\tilde{u}_{k}\frac{\partial}{\partial x_{l}}\left(\frac{1}{\bar{\rho}}\right)\right).$		(60)

With the new coefficient provided by equation (57), the anisotropic part of the Smagorinksy model can be determined, i.e.

$$\tau_{ij}^{mod}-\frac{1}{3}\delta_{ij}\tau_{kk}^{mod}=-2\bar{\rho}C_{sm}\Delta^{2}|\tilde{S}|\tilde{\mathbb{S}_{ij}}.$$

(61)

The isotropic part of the Smagorinksy model can be obtained directly from the relation $\tau_{kk}^{mod}=2\bar{\rho}k_{sgs}$.