# Non-equilibrium time-relaxation kinetic model for compressible turbulence
modeling

Guiyu Cao1,2  Liang Pan3  Kun Xu2,4  Minping Wan5,6<EMAIL_ADDRESS>Shiyi Chen1,5,6<EMAIL_ADDRESS>1Academy for Advanced Interdisciplinary
Studies, Southern University of Science and Technology, Shenzhen, Guangdong
518055, PR China 2Department of Mathematics, Hong Kong University of Science
and Technology, Clear Water Bay, Kowloon, Hong Kong 3Laboratory of Mathematics
and Complex Systems, School of Mathematical Sciences, Beijing Normal
University, Beijing 100875, PR China 4Shenzhen Research Institute, Hong Kong
University of Science and Technology, Shenzhen, Guangdong 518057, PR China
5Department of Mechanics and Aerospace Engineering, Southern University of
Science and Technology, Shenzhen, Guangdong 518055, PR China 6Guangdong-Hong
Kong-Macao Joint Laboratory for Data-Driven Fluid Mechanics and Engineering
Applications, Southern University of Science and Technology, Shenzhen,
Guangdong 518055, PR China

###### Abstract

For the first time, the non-equilibrium time-relaxation kinetic model (NTRKM)
is proposed for compressible turbulence modeling on unresolved grids. Within
the non-equilibrium time-relaxation framework, NTRKM is extended in the form
of modified Bhatnagar-Gross-Krook model. Based on the first-order Chapman-
Enskog expansion, NTRKM connects with the six-variable macroscopic governing
equations. The first five governing equations correspond to the conservative
laws in mass, momentum and total energy, while the sixth equation governs the
evolution of unresolved turbulence kinetic energy $K_{utke}$. The unknowns in
NTRKM, including turbulent relaxation time and source term, are determined by
essential gradient-type assumption and standard dynamic modeling approach.
Current generalized kinetic model on unresolved grids consequently offers a
profound mesoscopic understanding for one-equation subgrid-scale turbulence
kinetic energy $K_{sgs}$ model in compressible large eddy simulation. To solve
NTRKM accurately and robustly, a non-equilibrium gas-kinetic scheme is
developed, which succeeds the well-established gas-kinetic scheme for
simulating Navier-Stokes equations. Three-dimensional decaying compressible
isotropic turbulence and temporal compressible plane mixing layer on
unresolved grids are simulated to evaluate the generalized kinetic model and
non-equilibrium gas-kinetic scheme. The performance of key turbulent
quantities up to second-order statistics confirms that NTRKM is comparable
with the widely-used eddy-viscosity Smagorinsky model (SM) and dynamic
Smagorinsky model (DSM). Specifically, compared with the DNS solution in
temporal compressible plane mixing layer, the performance of NTRKM is much
closer with DSM and better than SM. This study provides a workable approach
for compressible turbulence modeling on unresolved grids, enriching the
understanding of turbulence modeling within the non-equilibrium time-
relaxation framework.

###### keywords:

time-relaxation kinetic model, compressible turbulence modeling, non-
equilibrium gas-kinetic scheme

##  1 Introduction

Turbulence modeling on unresolved grids is an extremely challenging issue in
turbulence community for decades (Pope, 2001). With the rapid increasing of
computational power, the large eddy simulation (LES) (Manabe et al., 1965;
Lilly, 1967) gradually becomes the tractable workhorse for high-fidelity
unsteady turbulence simulation. To simulate turbulent flows on unresolved
grids, LES solves the filtered Navier-Stokes (NS) equations with resolvable
large-scale turbulent structures explicitly, while the unresolved structures
are modeled through subgrid-scale (SGS) models (Sagaut, 2006; Garnier et al.,
2009). The widely-used eddy-viscosity LES models in physical space mainly
include zero-equation Smagorinsky-class models and one-equation SGS turbulence
kinetic energy (TKE) models (Schumann, 1975; Yoshizawa & Horiuti, 1985).

Smagorinsky model (SM) proposed by Manabe et al. (1965) belongs to zero-
equation eddy-viscosity model. SM models unresolved turbulent structures
through gradient-type assumption between the SGS stress and the resolved
velocity gradient. In practice, SM requires to adjust model coefficients
according to the flow types, and suffers the dissipative performance near the
wall, as well as SGS effect does not disappear in the laminar flow region
(Deardorff, 1970). To deal with those drawbacks of SM, dynamic Smagorinsky
model (DSM) (Germano et al., 1991; Moin et al., 1991; Lilly, 1992; Meneveau et
al., 1996) has been proposed for incompressible and compressible turbulence.
DSM allows the modeling coefficients to be computed locally on the basis of
dynamic approaches. In recent two decades, modified zero-equation eddy-
viscosity models have been constructed to be comparable with DSM while still
keep the simple algebraic form, such as wall-adapting local eddy-viscosity
model (Nicoud & Ducros, 1999), Vreman-type model (Vreman, 2004),
$\sigma$-model (Nicoud et al., 2011) and anisotropic minimum-dissipation model
(Rozema et al., 2015). By far, zero-equation eddy-viscosity models are the
most commonly used class of LES models (Moser et al., 2021)

Another branch to model the unresolved turbulent structures is deriving and
modeling SGS turbulence kinetic energy equation. Schumann (1975); Yoshizawa &
Horiuti (1985) have pioneered one-equation SGS turbulence kinetic energy
models to incorporate history and non-local effects through transport equation
of SGS turbulence kinetic energy $K_{sgs}$. These one-equation SGS turbulence
kinetic energy models can be analogous to the workable one-equation Reynolds
averaged Navier-Stokes (RANS) eddy-viscosity models (Wilcox et al., 1998). As
the grid filter width is taken as the characteristic modeling length scale,
only the SGS turbulence kinetic energy equation is required to determine the
eddy viscosity in LES models. One-equation SGS turbulence kinetic energy
models have been extensively applied in incompressible LES (Krajnović &
Davidson, 2002; De Stefano et al., 2008), which have shown better performance
in the prediction of turbulent flows. Compared with the well-established
research on compressibility correction for the unresolved TKE equation in the
RANS simulation (Sarkar et al., 1991; Wilcox et al., 1998), there only exists
limited work on one-equation SGS turbulence kinetic energy models for
compressible LES (Yoshizawa, 1986; Pomraning & Rutland, 2002; Chai & Mahesh,
2012; Cao et al., 2021). Considering the compressibility effects can hardly be
modeled in zero-equation eddy-viscosity model (Garnier et al., 2009), one-
equation SGS turbulence kinetic energy models indeed offer the great promise
in modeling compressible turbulent flows.

In the past decades, the gas-kinetic scheme (GKS) based on the Bhatnagar-
Gross-Krook (BGK) model (Bhatnagar et al., 1954; Chapman & Cowling, 1970) has
been developed systematically for the computations from low speed flows to
hypersonic ones (Xu, 2001, 2015). Based on the time-dependent flux solver,
including generalized Riemann problem solver and gas-kinetic scheme (Li & Du,
2016; Pan et al., 2016b), a reliable framework was provided for developing the
GKS into fourth-order accuracy. More importantly, the high-order gas-kinetic
scheme (HGKS) is as robust as the second-order one and works perfectly from
the subsonic to hypersonic viscous heat conducting flows (Cao et al., 2018).

With the advantage of finite volume GKS and HGKS, they have been naturally
implemented in simulating turbulent flows, especially in the compressible
regime. For practical turbulent flows, Hou et al. (1996); Chen et al. (2003)
pioneered the turbulent relaxation time $\tau_{t}$ for BGK-type models in
turbulence modeling within the equilibrium time-relaxation framework.
Following the concept of turbulent relaxation time, the second-order GKS/HGKS
coupled with $k-\omega$ SST model (Jiang & Qian, 2012; Righi, 2016; Cao et
al., 2019b), S-A model (Pan et al., 2016a), Vreman-type LES model, and the
hybrid RANS-LES method (Tan et al., 2018) have been performed in simulating
high-Reynolds number turbulence. These practical simulations have confirmed
the accuracy and robustness of second-order GKS/HGKS coupled with traditional
eddy-viscosity model. In terms of low-Reynolds number turbulent flows, the
second-order GKS/HGKS have been directly used as a DNS tool in simulating the
canonical benchmarks (Fu & Li, 2006; Liao et al., 2009; Kumar et al., 2013;
Cao et al., 2019a, 2021, 2022), such as compressible mixing layer,
compressible homogeneous turbulence, turbulent channel flows, etc. HGKS shows
special advantages in the supersonic turbulence due to its accuracy and super
robustness, i.e., the supersonic isotropic turbulence with initial turbulent
Mach number $Ma_{t0}=2.0$ has been simulated successfully (Cao et al., 2021).
Aiming to conduct the large-scale DNS, a parallel in-house code of HGKS has
been developed (Cao et al., 2022). Large-scale DNS up to $1024^{3}$ grids
shows that the computational cost of HGKS is comparable with the high-order
finite difference method (Bogey & Bailly, 2004).

To construct one-equation SGS turbulence kinetic energy model for compressible
LES within the time-relaxation framework, Cao et al. (2019a, 2021) have
systematically studied the high-fidelity DNS and delicate coarse-graining
analysis on decaying compressible isotropic turbulence. This paper aims to
complete the compressible one-equation $K_{sgs}$ model for LES. Firstly, non-
equilibrium time-relaxation kinetic model (NTRKM) is extended in the form of
modified BGK model. NTRKM can offer an mesoscopic understanding for transport
equation of the compressible SGS turbulence kinetic energy. To reasonably
maintain the accurate and robust numerical performance of HGKS, finite volume
non-equilibrium gas-kinetic scheme is designed when solving NTRKM. Comparable
with the widely-used eddy-viscosity SM and DSM, the decaying compressible
isotropic turbulence (DCIT) (Samtaney et al., 2001; Cao et al., 2019a) and
temporal compressible plane mixing layer (TCPML) (Sandham & Reynolds, 1991;
Vreman et al., 1997; Pantano & Sarkar, 2002) are simulated to assess current
generalized kinetic model and corresponding non-equilibrium gas-kinetic
scheme. These two cases are main engines to drive the development of
compressible turbulence models.

The organization of this paper is as follows. In §2, NTRKM for compressible
turbulence modeling is presented. §3 constructs finite volume non-equilibrium
gas-kinetic scheme for NTRKM. Posteriori tests on DCIT and TCPML are conducted
in §4. Conclusion and discussion are drawn in §5.

##  2 Non-equilibrium time-relaxation kinetic model for compressible
turbulence modeling

In this section, NTRKM on unresolved grids will be proposed. The first-order
Chapman-Enskog expansion provides the link between NTRKM and corresponding
macroscopic governing equations with six macroscopic variables.
Phenomenologically, the unknown turbulent relaxation time and source term can
be modeled through the gradient-type assumption and dynamic modeling approach.

###  2.1 Bhatnagar-Gross-Krook time-relaxation kinetic model

For molecular transport and collision, the simplification of Boltzmann
equation without external force is given by the BGK model (Bhatnagar et al.,
1954)

$\displaystyle\frac{\partial f}{\partial t}+u_{i}\frac{\partial f}{\partial
x_{i}}=\frac{g-f}{\tau},$ (1)

where $f(\bm{x},t,\bm{u},\xi)$ is the number density of molecular at position
$\bm{x}=(x_{1},x_{2},x_{3})^{T}$ and molecular velocity
$\bm{u}=(u_{1},u_{2},u_{3})^{T}$ at time $t$ with internal degrees of freedom
$\xi$. The relation between distribution function $f(\bm{x},t,\bm{u},\xi)$ and
macroscopic variables, such as mass, momentum and total energy can be obtained
by taking moments in molecular velocity of $f(\bm{x},t,\bm{u},\xi)$ (Xu, 2001,
2015). The left hand side of BGK model denotes the free transport process, and
the right hand side is the time-relaxation collision term. The collision term
in BGK model shows simple relaxation process from $f(\bm{x},t,\bm{u},\xi)$ to
a local equilibrium state $g$, with a molecular relaxation time $\tau$ which
is related to the molecular viscosity $\mu$ and heat conduction coefficient
$\kappa$ (see Appendix B (Xu, 2015)). The local equilibrium state $g$ is a
Maxwellian distribution

$\displaystyle
g=\rho(\frac{\lambda}{\pi})^{\frac{N+3}{2}}e^{-\lambda[(u_{i}-U_{i})^{2}+\xi^{2}]},$
(2)

where $\rho$ is the density, $\lambda=m_{o}/(2k_{B}T)$ as $m_{o}$ is the
molecular mass, $k_{B}$ the Boltzmann constant and $T$ the temperature,
$U_{i}$ denotes the macroscopic velocity in $x_{i}$ direction. For three-
dimensional equilibrium diatomic gas, the total number of degree of freedom in
$\xi$ is $N=2$, accounting for the two rotational modes
$\xi^{2}=\xi_{1}^{2}+\xi_{2}^{2}$. The specific heat ratio $\gamma$ is
determined by $\gamma=(N+5)/(N+3)$. Zeroth-order Chapman-Enskog expansion
(Chapman & Cowling, 1970) with $f=g$ offers the Euler equations. NS equations
can be derived with first-order truncation of Chapman-Enskog expansion

$\displaystyle f=g-\tau(\frac{\partial{g}}{\partial t}+u_{i}\frac{\partial
g}{\partial x_{i}}).$ (3)

For Euler and NS equations, the second-order and high-order GKS based on BGK
model (see equation (1)) has been systematically developed (Xu, 2001; Pan et
al., 2016b). The well-established second-order GKS/HGKS presents its accurate
and robust numerical performance from low speed flows to hypersonic ones (Xu,
2015).

\begin{overpic}[width=173.44534pt]{intro_numerical_turbulence_dns}
\put(-7.0,92.0){$(a)$} \end{overpic}

\begin{overpic}[width=174.3121pt]{intro_numerical_turbulence_les}
\put(-7.0,92.0){$(b)$} \end{overpic}

Figure 1:  Illustration of two-dimensional turbulence structure on $(a)$
resolved $8^{2}$ grids (DNS eliminates turbulence models entirely) and $(b)$
corresponding unresolved $1$ grid (turbulence modeling is required). Black
lines form grids and bule lines represent ”turbulent eddies”.

###  2.2 Non-equilibrium time-relaxation kinetic model

Numerically, the unresolved state or resolved state on a numerical cell
depends on the ratio of spatial-temporal resolution of numerical simulation to
the local characteristic scale of flow field. Spatial-temporal resolution is
mainly determined by the grid resolution and corresponding time step
(determined by CFL condition (Courant et al., 1928)), as well as the accuracy
of numerical scheme. In terms of characteristic scale, setting the boundary
layer as an example, thickness of boundary layer is reasonable characteristic
scale in laminar boundary layer (White & Majdalani, 2006), while smallest eddy
scale as Kolmogorov scale acts as intrinsic characteristic scale in turbulent
boundary layer Kim et al. (1987). When the numerical spatial-temporal
resolution is not adequate for resolving the local characteristic-scale
structures, the turbulence modeling is required. By contrast, DNS resolves
full scales above Kolmogorov scale, eliminating turbulence models entirely.

Figure 1 illustrates the comparison of two-dimensional turbulence structure in
”turbulent eddies” on resolved $8^{2}$ grids (see figure 1$(a)$) and
corresponding unresolved $1$ grid (see figure 1$(b)$). As demonstrated in
figure 1$(b)$, the unresolved grid means that the grid and corresponding time
step is not fine enough to resolve the local smallest turbulence structure
with fixed numerical scheme. Unresolved grids definitely leads to the lost of
turbulent information due to inevitable space and time averaging process when
updating the macroscopic variables (similar as the averaging process in finite
volume scheme (Xu, 2015)). The key point for turbulence modeling is to model
the unresolved turbulence structure through additional non-trivial quantities
on unresolved grids, i.e., non-trivial turbulent frequency governed by the
stochastic differential equation (Pope, 2001). The trivial quantities are
mass, momentum and total energy which are governed by the conservative laws,
without contributing non-trivial information to the unresolved turbulent
structures. Subsequently, on unresolved grids (see figure 1$(b)$), the non-
trivial unresolved turbulence kinetic energy $K_{utke}$ and its quantitative
dynamic evolution will be proposed for modeling the unresolved turbulent
process.

Figure 2:  $gk_{u}(k_{u})$ from coarsening DNS (Cao et al., 2021) at
$t/\tau_{to}=0.5$, $8^{3}$ resolved grids are coarsened to $1$ unresolved
grids with Box filter (Vreman et al., 1994). $\tau_{to}$ is the large-eddy-
turnover time.

To model turbulence on unresolved grids, the non-equilibrium double time-
relaxation kinetic model is proposed in the form of modified BGK model

$\displaystyle\frac{\partial f}{\partial t}+u_{i}\frac{\partial f}{\partial
x_{i}}=\frac{f^{eq}-f}{\tau+\tau_{t}}+\frac{g-f^{eq}}{\tau_{\ast}}\equiv\frac{f^{eq}-f}{\tau+\tau_{t}}+Q_{s},$
(4)

where $f(\bm{x},t,\bm{u},\xi,k_{u})$ is the generalized number density of
molecular on unresolved grids, $k_{u}$ the sample-space variable corresponding
to unresolved TKE $K_{utke}$, $\tau_{t}$ the turbulent relaxation time,
$\tau_{\ast}$ the newly-defined relaxation time from $f^{eq}$ to $g$, and
$Q_{s}$ the secondary relaxation term. Turbulence equilibrium state $f^{eq}$
is introduced with Maxwellian distribution $g$ for resolved flow variables and
Gamma distribution $gk_{u}$ for unresolved TKE as

$\displaystyle
f^{eq}=g\bm{\cdot}gk_{u}=\rho(\frac{\lambda}{\pi})^{\frac{N+3}{2}}e^{-\lambda[(u_{i}-U_{i})^{2}+\xi^{2}]}\bm{\cdot}\frac{1}{\Gamma(\alpha)}(\frac{\alpha}{K_{utke}})^{\alpha}k_{u}^{\alpha-1}e^{-\frac{\alpha
k_{u}}{K_{utke}}},$ (5)

where $\bm{\cdot}$ denotes the multiplication. We assume $gk_{u}(k_{u})$ is
the Gamma distribution, with non-negative shape parameter $\alpha$, mean
$Mean(k_{u})=K_{utke}$, variance $Var(k_{u})=K_{utke}^{2}/\alpha$, where
$K_{utke}$ is the total unresolved TKE on unresolved grids. In Jayesh-Pope
model (Pope, 2001), Gamma distribution is the stationary distribution of
turbulence frequency for statistically stationary isotropic turbulence. In
NTRKM, the distribution of $K_{utke}$ on unresolved grids is chosen as Gamma
distribution (see equation (5)), intuitively. Thence, the double-relaxation
process is named as non-equilibrium kinetic model, as the unresolved
turbulence information $K_{utke}$ participates in the non-equilibrium
relaxation process. Validity of Gamma distribution on $K_{utke}$ is conducted
through the following coarse-graining process. Based on the previous DNS study
on $384^{3}$ grids (Cao et al., 2021), the distribution of $K_{utke}$ through
coarsening DNS solution is presented in figure 2. $K_{utke}$ is normalized as
$1$ in figure 2. Gamma distribution with the parameter $\alpha=0.50$ is the
canonical distribution. Gamma distribution equipped with the parameter
$\alpha=0.39$ fits well with sample-space $k_{u}$ from coarsening DNS data. In
thick-tail region of figure 2, there exists apparent deviation between the
fitted Gamma distribution and coarsening DNS data. This deviation implies that
the intense events in the compressible isotropic turbulence (Wang et al.,
2013) are hard to be modeled by Gamma distribution. The optimal choice of
distribution for $K_{utke}$ on unresolved grids still requires to be
investigated. However, the form of $gk_{u}$ does not affect the evolution of
$K_{utke}$ with the subsequent finite volume non-equilibrium gas-kinetic
scheme. Since non-equilibrium gas-kinetic scheme acts as a hydrodynamic
solver, only the total $K_{utke}$ gets involved in the updating process
instead of $k_{u}$ in equation (5). By contrast, when the kinetic solver is
applied in updating the distribution function $f(\bm{x},t,\bm{u},\xi,k_{u})$
directly on unresolved grids, i.e., unified gas-kinetic scheme (UGKS) (Xu &
Huang, 2010), the form of $gk_{u}$ will contribute to the evolution of
$f(\bm{x},t,\bm{u},\xi,k_{u})$. If NTRKM is solved by kinetic solver, the
distribution of $K_{utke}$ requires to be modeled much carefully. As a
starter, turbulence equilibrium state $f^{eq}$ (see equation (5)) has been
proposed for constructing the non-trivial quantity on unresolved grid, namely,
depicting the $k_{u}$ for unresolved ”turbulent eddies” as illustrated in
figure 1$(b)$.

The relation between macroscopic variables as mass $\rho$, momentum $(\rho
U_{1},\rho U_{2},\rho U_{3})$, total energy $\rho E$, and unresolved
turbulence kinetic energy $\rho K_{utke}$ with the generalized distribution
function $f(\bm{x},t,\bm{u},\xi,k_{u})$ on unresolved grids is given by

$\displaystyle\bm{Q}=\int\bm{\psi}f\text{d}\Xi=\begin{pmatrix}\rho,\rho
U_{1},\rho U_{2},\rho U_{3},\rho E,\rho K_{utke}\end{pmatrix}^{T},$ (6)

where
$\bm{\psi}=(1,u_{1},u_{2},u_{3},\displaystyle\frac{1}{2}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+\xi^{2})+k_{u},k_{u})^{T}$
and
$\text{d}\Xi=\text{d}u_{1}\text{d}u_{2}\text{d}u_{3}\text{d}\xi\text{d}k_{u}$.
As the $K_{utke}$ is introduced to model unresolved turbulent process
quantitatively, one more constraint has to be imposed on current NTRKM to
self-consistently determine all unknowns. This additional constraint is the
$K_{utke}$ relaxation. Since only mass, momentum and total energy are
conserved during collisions, the compatibility condition for the collision
term becomes

$\displaystyle\bm{S}=\int\bm{\psi}(\frac{f^{eq}-f}{\tau+\tau_{t}}+Q_{s})\text{d}\Xi=(0,0,0,0,0,S_{t})^{T}.$
(7)

Unknown source term $S_{t}$ in equation (7) can be modeled through relaxation
model. This relaxation process is analogy to well-established non-equilibrium
kinetic model for multi-temperature flows (Xu et al., 2008) as

$\displaystyle S_{t}=\frac{\rho(K_{utke}^{eq}-K_{utke})}{\tau_{\ast}}.$ (8)

Conceptually, the equilibrium unresolved TKE $K_{utke}^{eq}$ and relaxation
time $\tau_{\ast}$ in equation (8) can be modeled on unresolved grids.
However, these two unknowns require profound priori knowledge and physical
understanding of turbulence within the non-equilibrium time-relaxation
framework. As far as the authors know, modeling $K_{utke}^{eq}$ and
$\tau_{\ast}$ in equation (8) directly is pretty challenging for current stage
of turbulence studies. To overcome this barrier, a comparison between the
$\rho K_{utke}$ equation derived from the first-order Chapman-Enskog expansion
on NTRKM and compressible $K_{sgs}$ equation from the compressible LES (Cao et
al., 2021) will be conducted. Consequently, source term $S_{t}$ can be modeled
in an alternative standard paradigm. After modeling $S_{t}$, the dynamic
evolution of non-trivial quantity $K_{utke}$ can be determined by equation (4)
and equation (7) quantitatively.

In contrast to the BGK model (1), the right-hand-side collision operator in
NTRKM equation (4) contains two terms corresponding to two-level collisions on
unresolved grids. The relaxation process has been extended as $f\to f^{eq}\to
g$, and the process from $f^{eq}\to g$ may take a much longer time
$\tau_{\ast}$ than that of process from $f\to f^{eq}$ by $(\tau+\tau_{t})$. On
unresolved grids, for the first collision process, the information of
unresolved $K_{utke}$ is memorized by the intermediate turbulence equilibrium
distribution $f^{eq}$; for the second collision process, $K_{utke}$ is
released into resolved kinetic energy and internal energy by $\tau_{\ast}$.
The total energy assignment in two-level collision process can be classified
as intermediate turbulence equilibrium state $\rho E=\rho\bm{U}^{2}/2+\rho
e+\rho K_{utke}$, and Maxwellian equilibrium state $\rho
E=\rho\bm{U}^{2}/2+\rho e$. $\bm{U}=(U_{1},U_{2},U_{3})^{T}$ is the resolved
macroscopic velocity vector with the definition of
$\bm{U^{2}}=U_{1}^{2}+U_{2}^{2}+U_{3}^{2}$. $e=(N+3)RT/2$ is the internal
energy in which $R$ is the gas constant. As figure 1$(b)$ approaches to figure
1$(a)$, the unresolved grids approach to resolved ones, meaning that the grid
and time step is fine enough to resolve the smallest spatial-temporal
characteristic structures. Under this approaching process, the unresolved
turbulent information approaches to disappear, so the non-trivial modeling
information on unresolved turbulence structures would be eliminated
automatically. This limit implies $K_{utke}\to 0$, turbulent relaxation time
$\tau_{t}\to 0$, turbulence equilibrium state $f^{eq}\to g$, and
$\tau_{\ast}\to\tau$. Therefore, this limit leads the NTRKM as equation (4) to
be consistent with the BGK model as equation (1) on resolved grid. This
asymptotic process also indicates turbulent relaxation time $\tau_{t}$ depends
on grid resolution and unresolved $K_{utke}$, which shed lights on the
modeling unknown turbulent relaxation time $\tau_{t}$.

###  2.3 Models for turbulent relaxation time $\tau_{t}$ and source term
$S_{t}$

To overcome the barrier of modeling the turbulent relaxation time $\tau_{t}$
and source term $S_{t}$ directly in kinetic model, the corresponding
macroscopic governing equations from the NTRKM will be derived by Chapman-
Enskog expansion. NTRKM to macroscopic equations can be regarded as a
projection, which essentially bridge the unknowns in NTRKM with the particular
terms in macroscopic governing equations.

Using the turbulence equilibrium state denoted as equation (5), with the
frozen of $K_{utke}$ exchange, the first-order Chapman-Enskog expansion gives

$\displaystyle f=f^{eq}-(\tau+\tau_{t})(\frac{\partial{f^{eq}}}{\partial
t}+u_{i}\frac{\partial f^{eq}}{\partial x_{i}}).$ (9)

From which the corresponding macroscopic governing equations in three
dimension have been derived for the first time, as shown in Appendix A, namely

$\displaystyle\rho_{,t}+(\rho U_{j})_{,j}=0,$ (10) $\displaystyle(\rho
U_{i})_{,t}+(\rho U_{i}U_{j}+p\delta_{ij})_{,j}=(\tau_{ij})_{,j},$ (11)
$\displaystyle(\rho E)_{,t}+((\rho
E+p)U_{j})_{,j}=(U_{i}\tau_{ij}+q_{j})_{,j},$ (12) $\displaystyle(\rho
K_{utke})_{,t}+(\rho K_{utke}U_{j})_{,j}=(U_{j}\tau_{st}+qk_{j})_{,j}+S_{t},$
(13)

where $p$ is pressure related to the resolved temperature $p=\rho
RT=\rho/(2\lambda)$, the total energy $\rho E=\rho(\bm{U}^{2}+3RT+NRT)/2+\rho
K_{utke}$, and $\delta_{ij}$ is the Kronecker symbol. The viscous stress term
in equation (11) is denoted by

$\displaystyle\tau_{ij}$
$\displaystyle=(\mu+\mu_{t})(U_{i,j}+U_{j,i}-\frac{2}{3}U_{k,k}\delta_{ij})+\eta_{t}U_{k,k}\delta_{ij}-\frac{2}{N+3}\tau_{st}\delta_{ij},$
(14)

with

$\displaystyle\tau_{st}$ $\displaystyle=(\tau+\tau_{t})S_{t},$ (15)

where molecular viscosity $\mu=\tau p$, turbulent eddy viscosity
$\mu_{t}=\tau_{t}p$, turbulent bulk viscosity
$\eta_{t}=2N(\mu+\mu_{t})/[3(N+3)]$, and the last term $\tau_{st}$ results
from the source term $S_{t}$. Typically, beyond the filtered compressible
momentum equation in NS equations (Chai & Mahesh, 2012), the generalized
viscous stress $\tau_{ij}$ on unresolved grids in equation (14) contains
additional terms. These terms are turbulent bulk viscosity term and term
related to the energy interaction from source term $S_{t}$. Unresolved
turbulence structure contributes to the generalized viscous stress so that
$\tau_{ij}$ on unresolved grids deviates from the linear constitutive relation
in equation (14). The heat conduction term in equation (12) reads

$\displaystyle q_{j}=(\kappa+\kappa_{t})T_{,j}+qk_{j},$ (16)

where molecular thermal conductivity is $\kappa=(N+5)\tau pR/2$, turbulent
thermal conductivity $\kappa_{t}=(N+5)\tau_{t}pR/2$. Appendix A shows that
Prandtl number $Pr=1$, and turbulent Prandtl number $Pr_{t}=1$. When
recovering the realistic laminar and turbulent Prandtl number, a similar
modification in the energy transport (Xu, 2001) should be implemented. The
modification will be presented briefly in §3. As presented in equation (13) ,
the term $qk_{j}$ is related to the governing equation of $K_{utke}$ as

$\displaystyle qk_{j}=(\mu+\mu_{t})(K_{utke})_{,j}.$ (17)

In summary, the first five governing equations in equations (10)-(12)
correspond to the conservative laws in mass, momentum and total energy with
generalized constitutive relationship of stress as equations (14)-(15) and
heat conduction term as equations (16)-(17). While the sixth equation (13)
governs the evolution of $K_{utke}$. So far, the unclosed terms in NTRKM are
the turbulent relaxation time $\tau_{t}$, and the source term $S_{t}$. In the
following part, $\tau_{t}$ and $S_{t}$ will be modeled through the comparison
amomg equations (10)-(13) and transport equations of the compressible SGS
turbulence kinetic energy $K_{sgs}$.

With the Favre filtering process (i.e., $\widetilde{U_{i}}$ denotes Favre
average of $U_{i}$), the compressible $\overline{\rho}K_{sgs}$ transport
equation (Cao et al., 2021) for compressible LES has been derived as

$\left\\{\begin{aligned}
&(\overline{\rho}K_{sgs})_{,t}+(\overline{\rho}K_{sgs}\widetilde{U}_{j})_{,j}=P_{sgs}-D_{sgs}+\Pi_{sgs}+T_{sgs},\\\
&P_{sgs}=-\tau_{ij}^{sgs}\widetilde{S}_{ij},\\\
&D_{sgs}=\overline{\sigma_{ij}U_{i,j}}-\overline{\sigma}_{ij}\widetilde{U}_{i,j},\\\
&\Pi_{sgs}=\overline{pU_{k,k}}-\overline{p}\widetilde{U}_{k,k},\\\
&T_{sgs}=[-\frac{1}{2}\overline{\rho}(\widetilde{U_{i}U_{i}U_{j}}-\widetilde{U_{i}U_{i}}\tilde{U}_{j})+\tau_{ij}^{sgs}\widetilde{U}_{i}\\\
&\qquad+(\overline{\sigma_{ij}U_{i}}-\overline{\sigma}_{ij}\widetilde{U}_{i})-\overline{\rho}R(\widetilde{TU_{j}}-\tilde{T}\widetilde{U}_{j})]_{,j},\end{aligned}\right.$
(18)

where $P_{sgs}$ is the production term, $D_{sgs}$ the total dissipation term,
$\Pi_{sgs}$ the pressure-dilation transfer, and last term $T_{sgs}$ the sum of
SGS diffusion terms. In equation (18), SGS turbulence kinetic energy
$\overline{\rho}K_{sgs}=\tau_{kk}^{sgs}/2=\overline{\rho}(\widetilde{U_{k}U_{k}}-\widetilde{U}_{k}\widetilde{U}_{k})/2$,
SGS stress
$\tau_{ij}^{sgs}=\overline{\rho}(\widetilde{U_{i}U_{j}}-\widetilde{U}_{i}\widetilde{U}_{j})$,
$\widetilde{S}_{ij}=(\widetilde{U}_{i,j}+\widetilde{U}_{j,i})/2$, and
$\sigma_{ij}=\mu\big{(}U_{i,j}+U_{j,i}-2U_{k,k}\delta_{ij}/3)$. The total SGS
dissipation rate $D_{sgs}$ can be decomposed into two parts, namely, SGS
solenoidal dissipation rate $\varepsilon_{s}^{sgs}$ and SGS dilational
dissipation rate $\varepsilon_{d}^{sgs}$ as

$\left\\{\begin{aligned} D_{sgs}&=\epsilon_{s}^{sgs}+\epsilon_{d}^{sgs},\\\
\varepsilon_{s}^{sgs}&=\overline{\mu}(\widetilde{\omega_{i}\omega_{i}}-\tilde{\omega}_{i}\tilde{\omega}_{i}),\\\
\varepsilon_{d}^{sgs}&=\frac{4}{3}\overline{\mu}(\widetilde{U_{k,k}^{2}}-\widetilde{U}_{k,k}^{2}),\end{aligned}\right.$
(19)

where $\omega_{i}=\epsilon_{ijk}U_{k,j}$ is the vorticity and
$\tilde{\omega}_{i}=\epsilon_{ijk}\widetilde{U}_{k,j}$ the filtered one with
the permutation symbol $\epsilon_{ijk}$. Comparing the governing equation of
$\rho K_{utke}$ in equation (13) with the exact $\overline{\rho}K_{sgs}$
equation in equation (18), it is seen that source term $S_{t}$ in equation
(13) is the net effect of SGS production term, SGS dissipation term and the
SGS pressure-dilation transfer as

$\displaystyle S_{t}\equiv P_{sgs}-D_{sgs}+\Pi_{sgs}.$ (20)

Consequently, current NTRKM as equation (4) provides an mesoscopic
understanding in transport equation of the compressible SGS turbulence kinetic
energy $K_{sgs}$ as equation (18). The non-trivial quantity $K_{utke}$ is
proposed for modeling the unresolved turbulence structures (see figure
1$(b)$), and the governing equation of $K_{utke}$ is responsible for the
evolution of unresolved turbulent process. This macroscopic description is
consistent with the projection of NTRKM, namely, the double-relaxation kinetic
model can be regarded as the mesoscopic understanding of one-equation SGS
$K_{sgs}$ model. We stress that the $K_{utke}$ on unresolved grids will be
defaultly regarded as SGS turbulence kinetic energy $K_{sgs}$ on filtered
grid. Especially, in the following modeling and simulations, grid filter width
adopts as the effective grid length of control volume (see §4.1 and §4.2) in
finite-volume numerical scheme, so $K_{sgs}$ can be treated as $K_{utke}$ by
default. Similarly, the SGS variables are treated equivalently as unresolved
variables without special statement.

As presented in equation (14), the connection between eddy viscosity $\mu_{t}$
and turbulent relaxation time $\tau_{t}$ is given by $\tau_{t}=\mu_{t}/p$. As
shown in figure 1$(b)$, $\tau_{t}$ can be explained as the relaxation time for
the turbulent eddies (Chen et al., 2003). The larger turbulent relaxation time
originates from the strong non-equilibrium turbulence process, i.e., eddies
transport and collision on unresolved grids. Following the seminal modeling
strategy (Yoshizawa, 1986; Chai & Mahesh, 2012), turbulent relaxation time
$\tau_{t}$ and SGS stress $\tau_{ij}^{sgs}$ can be modeled as

$\displaystyle\tau_{t}=\frac{C_{s}\overline{\Delta}\overline{\rho}K_{utke}^{\frac{1}{2}}}{p},$
(21)
$\displaystyle\tau_{ij}^{sgs}=-2C_{s}\overline{\Delta}\overline{\rho}K_{utke}^{\frac{1}{2}}\widetilde{S}_{ij}^{\ast}+\frac{2}{3}\overline{\rho}K_{utke}\delta_{ij},$
(22)

where $C_{s}$ is the model coefficient, $\overline{\Delta}$ the grid filter
width,
$\widetilde{S}_{ij}^{\ast}=\widetilde{S}_{ij}-\widetilde{S}_{kk}\delta_{ij}/3$
the traceless tensor of $\widetilde{S}_{ij}$. When correcting the total energy
transport to recover the realistic turbulent Prandtl number $Pr_{t}$, the
dynamic Prandtl number $Pr_{t}$ can be modeled as

$q_{j}=-C_{s}\overline{\Delta}\overline{\rho}K_{utke}^{\frac{1}{2}}\tilde{T}_{,j}/Pr_{t}.$
(23)

As shown in equation (21), with the aid of essential gradient-type assumption,
turbulent relaxation time $\tau_{t}$ has been closed in NTRKM. After modeling
the SGS stress $\tau_{ij}^{sgs}$, the SGS production term in source term in
equation (20) is modelled correspondingly. For the left unknowns in source
term $S_{t}$, the models of SGS dissipation rate and SGS pressure-dilation
transfer read (Chai & Mahesh, 2012)

$\displaystyle\epsilon_{s}^{sgs}=\frac{C_{\epsilon
s}\overline{\rho}K_{utke}^{\frac{3}{2}}}{\overline{\Delta}},$ (24)
$\displaystyle\epsilon_{d}^{sgs}=\frac{C_{\epsilon
d}\overline{\rho}Ma_{k}^{2}K_{utke}^{\frac{3}{2}}}{\overline{\Delta}},$ (25)
$\displaystyle\Pi_{sgs}=C_{\Pi}\overline{\Delta}^{2}\overline{p_{,j}}(\widetilde{U}_{k})_{j,k},$
(26)

where $Ma_{k}^{2}=2K_{utke}/(\gamma RT)$ is the unresolved TKE Mach number. In
terms of determining the unknown model coefficients, current paper follows the
standard dynamic approach (Germano et al., 1991; Moin et al., 1991; Chai &
Mahesh, 2012). In equations (21) - (26), model coefficients $C_{s}$, $C_{\Pi}$
and $Pr_{t}$ can be dynamically computed through Germano identity (Germano et
al., 1991; Moin et al., 1991). Additionally, $C_{\epsilon s}$ and $C_{\epsilon
d}$ can be obtained by the analogy between the grid-filter-level SGS stress
and the resolved stress across the test filter level (Menon & Kim, 1996; Chai
& Mahesh, 2012). The detailed derivation of all dynamic model coefficients and
necessary remarks are presented in Appendix B.

Turbulent relaxation time $\tau_{t}$ and source term $S_{t}$ have been modeled
on the basis of equations (20)-(26) with essential gradient-type assumption
and standard dynamic approaches. In the subsequent section, instead of solving
equations (10)-(13) with the traditional finite-volume hydrodynamic solver,
the NTRKM as equations (4) is solved directly with the flux function provided
by the time-dependent integral solution in the spirit of well-established gas-
kinetic scheme (Xu, 2001, 2015).

##  3 Non-equilibrium gas-kinetic scheme for generalized kinetic model

In this section, to maintain the accurate and robust numerical performance of
HGKS (Pan et al., 2016b; Cao et al., 2018), the finite volume non-equilibrium
gas-kinetic scheme is proposed to solve NTRKM.

For finite volume method, the key procedure is updating the macroscopic flow
variables inside each control volume through the numerical fluxes. Taking
moments of the NTRKM as equation (4) and integrating with respect to control
volume on unresolved grids, the finite volume scheme can be expressed as

$\displaystyle\frac{\text{d}(\bm{Q}_{ijk})}{\text{d}t}=-\frac{1}{|\Omega_{ijk}|}\sum_{s=1}^{6}\mathbb{F}_{s}(t)+\bm{S}_{ijk},$
(27)

where $\bm{Q}_{ijk}$ is the cell averaged macroscopic variables as equation
(6), $\bm{S}_{ijk}$ is cell averaged source term as equation (7) with $S_{t}$
modeled through equations (20)-(26). The control volume
$\Omega_{ijk}=[(x_{1})_{i}-\Delta x_{1}/2,(x_{1})_{i}+\Delta
x_{1}/2]\bm{\cdot}[(x_{2})_{j}-\Delta x_{2}/2,(x_{2})_{j}+\Delta
x_{2}/2]\bm{\cdot}[(x_{3})_{k}-\Delta x_{3}/2,(x_{3})_{k}+\Delta x_{3}/2]$,
$|\Omega_{ijk}|$ is the volume of $\Omega_{ijk}$ and $\mathbb{F}_{s}(t)$ is
the time-dependent numerical flux across the cell interface $\Sigma_{s}$. The
numerical flux $\mathbb{F}_{s}(t)|_{x_{1}}$ in $x_{1}$ direction (at cell
interface $(x_{1})_{i+1/2}$) is given as example

$\displaystyle\mathbb{F}_{s}(t)|_{x_{1}}$
$\displaystyle=\iint_{\Sigma_{s}|_{x_{1}}}\bm{F}(\bm{Q})\bm{\cdot}\bm{n}\text{d}\sigma$
(28)
$\displaystyle=\sum_{m,n=1}^{2}\omega_{mn}\int\bm{\psi}u_{1}f(\bm{x}_{i+1/2,j_{m},k_{n}},t,\bm{u},\xi,k_{u})\text{d}\Xi\Delta
x_{2}\Delta x_{3},$

where $\bm{n}$ is the outer normal direction. The Gaussian quadrature is used
over the cell interface for equation (28), where $\omega_{mn}$ is the
quadrature weight,
$\bm{x}_{i+1/2,j_{m},k_{n}}=[(x_{1})_{i+1/2},(x_{2})_{j_{m}},(x_{3})_{k_{n}}]^{T}$,
and $[(x_{2})_{j_{m}},(x_{3})_{k_{n}}]$ is the quadrature point of cell
interface $[(x_{2})_{j}-\Delta x_{2}/2,(x_{2})_{j}+\Delta
x_{2}/2]\bm{\cdot}[(x_{3})_{k}-\Delta x_{3}/2,(x_{3})_{k}+\Delta x_{3}/2]$.
When constructing the numerical fluxes, the secondary relaxation term $Q_{s}$
in equation (4) is not considered, and the effect of $Q_{s}$ is taken into
account as the source term in equation (27). The gas distribution function
$f(\bm{x}_{i+1/2,j_{m},k_{n}},t,\bm{u},\xi,k_{u})$ in the local coordinate can
be obtained by the integral solution of equation (4) as

$\displaystyle f(\bm{x}_{i+1/2,j_{m},k_{n}},t,\bm{u},\xi,k_{u})$
$\displaystyle=\frac{1}{(\tau+\tau_{t})}\int_{0}^{t}f^{eq}(\bm{x}^{\prime},t^{\prime},\bm{u},\xi,k_{u})e^{-(t-t^{\prime})/(\tau+\tau_{t})}\text{d}t^{\prime}$
(29) $\displaystyle+e^{-t/(\tau+\tau_{t})}f_{0}(-\bm{u}t,\xi,k_{u}),$

where $\bm{x}^{\prime}=\bm{x}_{i+1/2,j_{m},k_{n}}-\bm{u}(t-t^{\prime})$ is the
trajectory of molecular on unresolved grids, $f_{0}$ the initial gas
distribution function, and $f^{eq}$ the corresponding turbulence equilibrium
state in the form of equation (5). Along the line of GKS (Xu, 2001; Cao et
al., 2018), for the multi-dimensional kinetic solver, $f^{eq}$ and $f_{0}$ can
be constructed as

$\displaystyle
f^{eq}=f^{eq}_{0}(1+\overline{a}_{1}x_{1}+\overline{a}_{2}x_{2}+\overline{a}_{3}x_{3}+\overline{A}t),$
(30)

and

$\displaystyle
f_{0}=\begin{cases}f^{eq}_{l}[1+(a_{1}^{l}x_{1}+a_{2}^{l}x_{2}+a_{3}^{l}x_{3})-(\tau+\tau_{t})(a_{1}^{l}u_{1}+a_{2}^{l}u_{2}+a_{3}^{l}u_{3}+A_{l})],&x\leq
0,\\\
f^{eq}_{r}[1+(a_{1}^{r}x_{1}+a_{2}^{r}x_{2}+a_{3}^{r}x_{3})-(\tau+\tau_{t})(a_{1}^{r}u_{1}+a_{2}^{r}u_{2}+a_{3}^{r}u_{3}+A_{r})],&x>0,\end{cases}$
(31)

where $f^{eq}_{l}$ and $f^{eq}_{r}$ are the initial gas distribution functions
on both sides of a cell interface $\Sigma_{s}$. $f^{eq}_{0}$ is the initial
turbulence equilibrium state located at the cell interface, which can be
determined through the compatibility condition

$\displaystyle\int\bm{\psi}f^{eq}_{0}\text{d}\Xi=\int_{u_{1}>0}\bm{\psi}f^{eq}_{l}\text{d}\Xi+\int_{u_{1}<0}\bm{\psi}f^{eq}_{r}\text{d}\Xi.$
(32)

Substituting $f^{eq}$ (see equation (30)) and $f_{0}$ (see equation (31)) into
equation (29), the time-dependent gas distribution function at the Gaussian
point is evaluated as

$\displaystyle
f(\bm{x}_{i+1/2,j_{m},k_{n}},t,\bm{u},\xi,k_{u})=(1-e^{-t/(\tau+\tau_{t})})f^{eq}_{0}$
(33)
$\displaystyle+[(t+\tau+\tau_{t})e^{-t\tau}-\tau](\overline{a}_{1}u_{1}+\overline{a}_{2}u_{2}+\overline{a}_{3}u_{3})f^{eq}_{0}$
$\displaystyle+[t-(\tau+\tau_{t})+(\tau+\tau_{t})e^{-t(\tau+\tau_{t})}]\overline{A}f^{eq}_{0}$
$\displaystyle+e^{-t/(\tau+\tau_{t})}f^{eq}_{l}[1-(\tau+\tau_{t}+t)(a_{1}^{l}u_{1}+a_{2}^{l}u_{2}+a_{3}^{l}u_{3})-(\tau+\tau_{t})A_{l}]H(u)$
$\displaystyle+e^{-t/(\tau+\tau_{t})}f^{eq}_{r}[1-(\tau+\tau_{t}+t)(a_{1}^{r}u_{1}+a_{2}^{r}u_{2}+a_{3}^{r}u_{3})-(\tau+\tau_{t})A_{r}](1-H(u)).$

With the relation of macroscopic variables and turbulence equilibrium
distribution function $f^{eq}$, the spatial mesoscopic coefficients
$\overline{a}_{1}$, $a_{1}^{l}$, $\cdots$, $a_{3}^{l}$, $a_{3}^{r}$ and
temporal mesoscopic coefficients $\overline{A}$, $A_{l}$, $A_{r}$ in equation
(33) can be determined and details are presented in Appendix C. Equation (33)
provides a gas evolution process from kinetic scale to hydrodynamic scale on
unresolved grids, where both inviscid and viscous fluxes are recovered from a
time-dependent and multi-dimensional gas distribution function at a cell
interface. This flux function couples the inviscid and all dissipative terms
(Xu, 2001; Cao et al., 2018), and has advantages in comparison with
traditional hydrodynamic solver in which the Riemann solver and central
difference are used for the inviscid and viscous terms. For Prandtl number
fix, both the laminar Prandtl number $Pr$ and turbulent Prandtl number
$Pr_{t}$ should be taken into consideration. Total energy flux $\bm{F}{(\rho
E)}$ in equation (28) should be modified as $\bm{F}^{new}{(\rho
E)}=\bm{F}{(\rho E)}+\\{(\mu Pr_{t}+\mu_{t}Pr)/[PrPr_{t}(\mu+\mu_{t})]-1\\}q$,
where the time-dependent heat flux can be evaluated precisely by
$q=\int(u-U)\\{[(u_{i}-U_{i})^{2}+\xi^{2}]/2+k_{u}\\}fd\Xi$.

The second-order accuracy in time can be achieved by one step integration,
with the time-dependent kinetic flux as equation (33). To achieve high-order
accuracy in space and time, the fifth-order WENO-Z spatial reconstruction
(Castro et al., 2011) and two-stage fourth-order time discretization (Li & Du,
2016; Pan et al., 2016b) are implemented. The characteristic reconstruction is
applied to improve the robustness for compressible flows with strong
discontinuities (i.e., DCIT) (Pan & Xu, 2020). The characteristic variables
are defined as $\bm{Q}_{c}=\bm{R}^{-1}\bm{Q}$, where $\bm{R}$ is the right
eigenmatrix of Jacobian matrix at Gaussian quadrature point and details are
given in Appendix D. When dealing with compressible flows without strong
discontinuities, such as TCPML, the linear WENO spatial reconstruction based
on conservative variables is adopted. For source term in equation (27), the
one-step forward Euler method is applied in two-stage updating process to
guarantee the robustness. Therefore, the finite volume non-equilibrium gas-
kinetic scheme has been constructed with the second-order kinetic flux, fifth-
order WENO-Z reconstruction, two-stage fourth-order time discretization and
one-step forward Euler method for source term. The current non-equilibrium
gas-kinetic scheme has been well implemented in the in-house platform for
turbulence simulation (Cao et al., 2022), and the posteriori tests on
compressible turbulent flows will be presented in the following section.

##  4 Posteriori tests

In this section, the decaying compressible isotropic turbulence (Samtaney et
al., 2001; Cao et al., 2019a) and temporal compressible plane mixing layer
(Sandham & Reynolds, 1991; Vreman et al., 1997; Pantano & Sarkar, 2002) are
regarded as cornerstones to assess NTRKM and non-equilibrium gas-kinetic
scheme.

###  4.1 Decaying compressible isotropic turbulence

DCIT (Samtaney et al., 2001) is the building-block case to demonstrate the
performance of modeling on compressible turbulence. For the flow with
discontinuities, we have

$\displaystyle\tau+\tau_{t}=\frac{\mu+\mu_{t}}{p}+C_{num}|\frac{p_{l}-p_{r}}{p_{l}+p_{r}}|\text{d}t,$
(34)

where $p$ is the pressure at the cell interface, $p_{l}$ and $p_{r}$ the
pressure on the left and right sides of the cell interface. $\text{d}t$ is the
time step, and a fixed $C_{num}=2.5$. The reason for including artificial
dissipation through the additional term in the molecular relaxation time
$\tau$ and the turbulent relaxation time $\tau_{t}$ is to improves the
numerical stability. As the earlier remark states, the grid filter width
adopts as the grid length of control volume for DCIT, i.e.,
$\overline{\Delta}=\Delta x_{1}=\Delta x_{2}=\Delta x_{3}$ on equivalent
spaced grids. Additionally, test filter width $\widehat{\overline{\Delta}}$ is
set to twice the grid filter width $\overline{\Delta}$ when determining the
dynamic model coefficients, namely
$\widehat{\overline{\Delta}}=2\overline{\Delta}$. DNS for DCIT has been well
studied using HGKS systematically (Cao et al., 2019a, 2021). In this section,
following previous DNS set-up, LES on unresolved grids will be conducted
directly.

Case | grid size | $K_{0}$ | $I_{0}$ | $\kappa_{max}\eta_{0}$
---|---|---|---|---
DNS | $512^{3}$ | 0.5055 | 0 | 3.6
R1 | $128^{3}$ | 0.4931 | 0.0151 | 0.90
Table 1:  Numerical parameters for DCIT.

In the computation, the initial Taylor microscale Reynolds number is
$Re_{\lambda 0}=72$ and the initial turbulent Mach number is fixed at
$Ma_{t0}=0.6$. The detailed initial conditions are set as previous work (Cao
et al., 2019a), and the periodic boundary condition for six macroscopic
variables is used. A three-dimensional solenoidal random initial velocity
field can be generated by a specified spectrum as

$\displaystyle E(\kappa)=A_{0}\kappa^{4}\exp(-2\kappa^{2}/\kappa_{0}^{2}),$
(35)

with the fixed $A_{0}=0.00013$ and $\kappa_{0}=8$. After generating the
initial velocity field on $512^{3}$ resolved grids, the filtered velocity
fields can be obtained on unresolved grids, i.e., filtered flow fields on
$128^{3}$ grids. When filtering velocity field, the positive definite kernel
of Box filter is adopted to guarantee the positive unresolved TKE (Vreman et
al., 1994), thus the initial pointwise $K_{utke0}$ on unresolved grids can be
obtained as the initial condition for equation (13). Table 1 shows the
numerical parameters for DCIT of DNS and R1, where $\kappa_{max}$ is maximum
resolved wave number and $\eta_{0}$ is the Kolmogorov length scale.
$\left\langle K_{0}\right\rangle$ is the initial ensemble resolved TKE in
which $\left\langle\cdot\right\rangle$ denotes the spatial average on the
whole computational domain. Turbulence intensity $I_{0}$ denotes the ratio of
initial ensemble unresolved $\left\langle K_{utke0}\right\rangle$ to the
initial ensemble resolved TKE as $I_{0}=\left\langle
K_{utke0}\right\rangle/\left\langle K_{0}\right\rangle$. Table 1 shows that
the grid resolution meets the DNS criterion $\kappa_{max}\eta_{0}\geq 2.71$
for DCIT (Cao et al., 2019a). Obviously, the grid resolution of R1 is not
adequate for DNS, which is regarded as compressible LES on unresolved grids.
The representative key statistical quantities, including the resolved root-
mean-square density fluctuation $\rho_{rms}$ and resolved turbulence kinetic
energy $K$, are given by

$\displaystyle\rho_{rms}=\left\langle(\rho-\left\langle\rho\right\rangle)^{2}\right\rangle^{\frac{1}{2}},$
(36) $\displaystyle K=\frac{1}{2}\rho\bm{U}^{2}.$ (37)

The ensemble budget of resolved $K$ is computed, which can be described
approximately by (Sarkar et al., 1991)

$\displaystyle\frac{\text{d}\left\langle
K\right\rangle}{\text{d}t}=-\left\langle\varepsilon\right\rangle+\left\langle
pU_{k,k}\right\rangle,$ (38)
$\displaystyle\varepsilon=\varepsilon_{s}+\varepsilon_{d},$ (39)

where $\varepsilon_{s}=\mu\omega_{i}\omega_{i}$ is the resolved solenoidal
dissipation rate, $\displaystyle\varepsilon_{d}=4\mu U_{k,k}^{2}/3$ the
resolved dilational dissipation rate without considering bulk viscosity,
$pU_{k,k}$ the resolved pressure-dilation transfer.

\begin{overpic}[width=214.63977pt]{cit_ntrkm_cs} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_ntrkm_kutke}
\put(-2.0,75.0){$(b)$} \end{overpic}

Figure 3:  Time history of $(a)$ dynamic coefficient $C_{s}$ and $(b)$
normalized ensemble unresolved $\left\langle
K_{utke}\right\rangle/\left\langle K_{utke0}\right\rangle$ for case R1 with
NTRKM.

Figure 3 shows the time history of dynamic coefficient $C_{s}$ as equation
(21) for turbulent relaxation time $\tau_{t}$ and normalized ensemble
$\left\langle K_{utke}\right\rangle/\left\langle K_{utke0}\right\rangle$.
$\tau_{cit}=t/\tau_{to}$ is the normalized time and $\tau_{to}$ is the large-
eddy-turnover time (Cao et al., 2019a). Firstly, model coefficient $C_{s}$ is
presented to validate the implementation of dynamic modeling approach as
Appendix B.2. The empirical model coefficient $C_{s}$ is recommended as a
fixed value $0.05$ (Yoshizawa & Horiuti, 1985). Figure 3$(a)$ shows that the
dynamic coefficient $C_{s}$ in NTRKM fluctuates between $[0.06,0.12]$ for case
R1. Current dynamic approach shows that $C_{s}$ reasonably depends on the grid
resolution and the evolution of flow fields. In figure 3$(b)$, the normalized
ensemble unresolved $\left\langle K_{utke}\right\rangle/\left\langle
K_{utke0}\right\rangle$ increases approximately within $\tau_{cit}\leq 1.5$
and decrease consecutively, which behaves similarly as previous literature
(Chai & Mahesh, 2012). In equation (18), the SGS production term represents
the inter-scale transfer associated with the interaction of resolved and
unresolved scales. The SGS dissipation terms act as the sink of $K_{utke}$ in
source term $S_{t}$. The evolution of $\left\langle
K_{utke}\right\rangle/\left\langle K_{utke0}\right\rangle$ implies that the
ensemble forward resolved energy cascade dominates at the early stage, and
then the SGS dissipation terms dominate. Figure 3$(b)$ indicates that the
intrinsic equilibrium assumption on $K_{utke}$, such as $\left\langle
S_{t}\right\rangle\approx 0$ for zero-equation eddy-viscosity LES models
(Lilly, 1967; Germano et al., 1991; Moin et al., 1991) is not valid, which
confirms that the evolution of $K_{utke}$ on unresolved grids is crucial for
compressible LES modeling (Yoshizawa & Horiuti, 1985; Chai & Mahesh, 2012).

\begin{overpic}[width=214.63977pt]{cit_ntrkm_kutke_sketch}
\put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_ntrkm_st}
\put(-2.0,75.0){$(b)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_ntrkm_production}
\put(-2.0,75.0){$(c)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_ntrkm_d1}
\put(-2.0,75.0){$(d)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_ntrkm_d2}
\put(-2.0,75.0){$(e)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_ntrkm_ptheta}
\put(-2.0,75.0){$(f)$} \end{overpic}

Figure 4:  Three-dimension contours of $(a)$ unresolved $K_{utke}$, $(b)$
source term $S_{t}$, and the components of source term $(c)$ $P_{sgs}$, $(d)$
$\varepsilon_{s}^{sgs}$, $(e)$ $\varepsilon_{d}^{sgs}$, $(f)$ $\Pi_{sgs}$ at
$\tau_{cit}=0.5$ for case R1 with NTRKM.

Figure 4 shows the contours of unresolved $K_{utke}$, source term $S_{t}$, and
the components of source term $P_{sgs}$, $\varepsilon_{s}^{sgs}$,
$\varepsilon_{d}^{sgs}$, $\Pi_{sgs}$ at $\tau_{cit}=0.5$ for case R1. Figure
4$(a)$ illustrates the contour of unresolved $K_{utke}$, and figure 4$(b)$
confirms the ensemble positive source term $S_{t}$ which accounts for the
increase of unresolved $K_{utke}$ in figure 3$(b)$. More specifically, figure
4$(c)$ presents that the magnitude of ensemble unresolved production rate at
$\tau_{cit}=0.5$ is larger than that of ensemble unresolved dissipation rate
and ensemble unresolved pressure-dilation transfer. Qualitatively, the
negative unresolved dissipation rate and the high similarity between
unresolved solenoidal dissipation rate and unresolved dilational dissipation
rate are observed in figure 4$(d)(e)$. Figure 4$(f)$ shows that the magnitude
and portion of negative unresolved pressure-dilation transfer are larger than
the positive ones, which behave similarly as delicate priori coarse-graining
analysis of compressible unresolved TKE budget (Cao et al., 2021).

\begin{overpic}[width=214.63977pt]{cit_les_rhorms} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_les_k}
\put(-2.0,75.0){$(b)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_les_dtotal}
\put(-2.0,75.0){$(c)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_les_d1}
\put(-2.0,75.0){$(d)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_les_d2}
\put(-2.0,75.0){$(e)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cit_les_ptheta}
\put(-2.0,75.0){$(f)$} \end{overpic}

Figure 5:  Time history of $(a)$ normalized resolved root-mean-square density
fluctuation $\rho_{rms}/Ma_{t0}^{2}$, $(b)$ normalized ensemble resolved
turbulence kinetic energy $\left\langle K\right\rangle/\left\langle
K_{0}\right\rangle$, $(c)$ ensemble resolved dissipation rate
$\left\langle\varepsilon\right\rangle$, $(d)$ ensemble resolved solenoidal
dissipate rate $\left\langle\varepsilon_{s}\right\rangle$, $(e)$ ensemble
resolved dilational dissipation rate
$\left\langle\varepsilon_{d}\right\rangle$, and $(f)$ ensemble resolved
pressure-dilation transfer $\left\langle p\theta\right\rangle$ for case R1.

To evaluate the performance of NTRKM and non-equilibrium gas-kinetic scheme,
figure 5 shows the key resolved statistical quantities as equations (36)-(39)
for case R1. Figure 5 shows that the key statistical quantities of NTRKM are
comparable with those from the widely-used SM (Manabe et al., 1965) and DSM
(Moin et al., 1991). SM and DSM are dealt with the equilibrium time-relaxation
framework (Cao et al., 2019b), which are also implemented in the in-house DNS
code (Cao et al., 2022) with modifying the $\tau$ to $\tau+\tau_{t}$ on
unresolved grids. Among these three compressible LES models, the heat flux is
not modified, namely, $Pr=1$ and $Pr_{t}=1$ are treated fairly on all
simulations. In terms of density-weighted SM, the eddy viscosity takes
$\mu_{t}=(C_{sm}\overline{\Delta})^{2}\overline{\rho}|\widetilde{S}|$ with
magnitude
$|\widetilde{S}|=(2\widetilde{S}_{ij}\widetilde{S}_{ij})^{\frac{1}{2}}$.
Coefficient $C_{sm}$ is fixed as $0.17$ deduced from the scaling law of TKE
spectrum in high-Reynolds number isotropic turbulence (Lilly, 1967). For
density-weighted DSM, the eddy viscosity is determined as
$\mu_{t}=C_{dsm}\overline{\Delta}^{2}\overline{\rho}|\widetilde{S}|$, and the
dynamic coefficient $C_{dsm}$ is computed through the dynamic technique
presented in Appendix B.1. The dynamic coefficient $C_{dsm}$ in density-
weighted DSM evolves between $[0.016,0.025]$ for case R1. Thus, the dynamic
coefficient $C_{dsm}$ is slightly smaller than the empirical static
coefficient $C_{sm}^{2}\approx 0.029$ in current DCIT, which accounts for the
comparable performance between density-weighted SM and density-weighted DSM on
unresolved grids. Overall, the performance of key statistical turbulent
quantities shows that current NTRKM is comparable with the widely-used eddy-
viscosity SM and DSM. Numerical performance of NTRKM and corresponding non-
equilibrium gas-kinetic scheme (see §3) offers the confidence for simulating
practical turbulence on unresolved grids.

###  4.2 Temporal compressible plane mixing layer

For temporal compressible plane mixing layer (Sandham & Reynolds, 1991; Vreman
et al., 1997; Pantano & Sarkar, 2002), practical simulations on unresolved
grids are conducted to further assess the performance of NTRKM and the non-
equilibrium gas-kinetic scheme (see §3). For the flow without strong
discontinuities, the collision time is given by

$\displaystyle\tau+\tau_{t}=\frac{\mu+\mu_{t}}{p}.$ (40)

For TCPML, the grid filter width adopts as the effective grid length of
control volume, i.e., $\overline{\Delta}=(\Delta x_{1}\times\Delta
x_{2}\times\Delta x_{3})^{1/3}$. Test filter width
$\widehat{\overline{\Delta}}$ still keeps to twice the grid filter width
$\overline{\Delta}$ in determining the dynamic model coefficients, namely
$\widehat{\overline{\Delta}}=2\overline{\Delta}$. In this section, DNS in
TCPML will be validated firstly. Then, LES studies restart from the filtered
DNS solution on unresolved grids.

TCPML is initialized by a hyperbolic tangent profile for the streamwise
velocity (Arun et al., 2019) as

$\left\\{\begin{aligned} U_{1}&=\frac{1}{2}\Delta
U\tanh(\frac{-x_{2}}{2\delta_{\theta 0}}),\\\ U_{2}&=0,\\\
U_{3}&=0,\end{aligned}\right.$ (41)

where $\Delta U=U_{lo}-U_{up}$, and initial momentum thickness $\delta_{\theta
0}=1$ is adopted. As equation (41) presents, two equal and opposite streamwise
velocities are simulated as $-U_{up}=U_{lo}=1$. With the unity Prandtl number,
Crocco-Busemann relation (Sandham, 1990) gives the initial temperature profile
as

$\displaystyle\frac{T}{T_{\infty}}=1+Ma_{c}^{2}\frac{\gamma-1}{2}(1-U_{1}^{2}).$
(42)

The initial density is set to a uniform value $\rho_{\infty}$. The convective
Mach number $Ma_{c}=0.75$ and initial vorticity thickness-based Reynolds
number $Re_{\omega 0}=640$ are simulated as

$\displaystyle Ma_{c}=\frac{\Delta U}{2\sqrt{\gamma RT_{\infty}}},$ (43)
$\displaystyle Re_{\omega 0}=\frac{\rho_{\infty}\Delta U\delta_{\omega
0}}{\mu_{\infty}},$ (44)

where $\mu_{\infty}$ is the reference viscosity corresponding to reference
temperature $T_{\infty}$, and viscosity is determined by power law as
$\mu(T)=\mu_{\infty}(T/T_{\infty})^{0.67}$ (Sandham & Reynolds, 1991). With
the uniform initial density, initial vorticity thickness can be estimated as
$\delta_{\omega 0}=\Delta U/|\partial U_{1}/\partial x_{2}|_{max}$ in which
the maximum of denominator is reached in the centre plane. $|\cdot|$
represents the absolute value. The momentum thickness $\delta_{\theta}$ is
defined as

$\displaystyle\delta_{\theta}=\frac{1}{\rho_{\infty}(\Delta
U)^{2}}\int_{-\infty}^{\infty}[\left\langle\rho
U_{1}(x_{2})\right\rangle-\left\langle\rho
U_{lo}\right\rangle][\left\langle\rho U_{up}\right\rangle-\left\langle\rho
U_{1}(x_{2})\right\rangle]\text{d}x_{2},$ (45)

where $\left\langle\cdot\right\rangle$ represents the plane average along the
streamwise and spanwise directions, and $\delta_{\omega 0}\approx
4\delta_{\theta 0}$ since the finite transverse domain $[-L_{1}/2,L_{1}/2]$.
For TCPML, the turbulent stress tensor $R_{ij}$ and anisotropy stress tensor
$b_{ij}$ read

$\displaystyle b_{ij}=\frac{R_{ij}-\frac{2}{3}K_{R}\delta_{ij}}{2K_{R}},$ (46)
$\displaystyle R_{ij}=\frac{\left\langle\rho
U_{i}^{{}^{\prime}}U_{j}^{{}^{\prime}}\right\rangle}{\left\langle\rho\right\rangle},$
(47)

where $U_{i}^{{}^{\prime}}=U_{i}-\left\langle\rho
U_{i}\right\rangle/\left\langle\rho\right\rangle$, $K_{R}$ the so-called
resolved turbulence kinetic energy as $K_{R}=R_{ii}/2$. Anisotropy stress
tensor $b_{ij}$ is an important characteristic of turbulence, especially for
advanced turbulence closures (Pantano & Sarkar, 2002). In following
statistical process, $R_{ij}$ is integrated across mixing layer within
$[-\delta_{\omega}(\tau_{ml}),\delta_{\omega}(\tau_{ml})]$, while $b_{ij}$ is
integrated within $[-4\delta_{\theta}(\tau_{ml}),4\delta_{\theta}(\tau_{ml})]$
with normalized time $\tau_{ml}=\Delta Ut/\delta_{\theta 0}$. To accelerate
the transition process, the initial condition is specified by adding a random
number to density, temperature, transverse and spanwise velocities at each
mesh point (Sandham & Reynolds, 1991), i.e.,
$\rho_{p}=\rho_{\infty}+A_{s1}r_{d}e^{-x_{2}/(2\delta_{\theta 0})^{2}}$.
$r_{d}$ is a random number uniformly distributed between $[-0.5,0.5]$ and the
amplitude $A_{s1}=0.2$. In terms of streamwise velocity, besides the random
number, the artificial sinusoidal-type perturbation has been added as

$\displaystyle U_{1p}$
$\displaystyle=U_{1}+U_{1}[A_{s1}+A_{s2}sin(\gamma_{s1}x_{2})(B_{s1}+B_{s2})]r_{d}e^{-x_{2}/(2\delta_{\theta
0})^{2}},$ (48)

where $A_{s2}=0.6$, $\gamma_{s1}=0.25$, and
$B_{s1}=A_{s3}[cos(\gamma_{s2}x_{1})+cos(2\gamma_{s2}x_{1})+cos(4\gamma_{s2}x_{1})]$,
$B_{s2}=A_{s4}[cos(\gamma_{s2}x_{1})cos(\gamma_{s2}x_{3})+cos(2\gamma_{s2}x_{1})cos(2\gamma_{s2}x_{3})+cos(4\gamma_{s2}x_{1})cos(4\gamma_{s2}x_{3})$
with $A_{s3}=0.2$, $A_{s4}=0.4$ and $\gamma_{s2}=0.235$. The initial condition
for primitive variables $(\rho_{p},U_{1p},U_{2p},U_{3p},T_{p})^{T}$ can be
obtained for DNS. The computational domain is discretized uniformly in three
directions. Boundary conditions in the homogeneous streamwise and spanwise
directions are periodic. In the transverse direction, the non-reflective
boundary condition of conservative variables is given according to one-
dimensional Riemann invariants (Toro, 2013), whereas the outlet boundaries are
used for the unresolved $K_{utke}$ (see equation (13)).

Case | $L_{1}\times L_{2}\times L_{3}$ | $N_{1}\times N_{2}\times N_{3}$ | $Ma_{c}$ | $Re_{\omega 0}$
---|---|---|---|---
Ref1 | $172\delta_{\theta 0}\times 129\delta_{\theta 0}\times 86\delta_{\theta 0}$ | $256\times 192\times 128$ | 0.70 | 640
Ref2 | $314\delta_{\theta 0}\times 157\delta_{\theta 0}\times 78.5\delta_{\theta 0}$ | $512\times 256\times 128$ | 0.75 | 640
DNS | $314\delta_{\theta 0}\times 157\delta_{\theta 0}\times 78.5\delta_{\theta 0}$ | $576\times 384\times 192$ | 0.75 | 640
M1 | $314\delta_{\theta 0}\times 157\delta_{\theta 0}\times 78.5\delta_{\theta 0}$ | $144\times 96\times 48$ | 0.75 | 640
Table 2:  Numerical parameters for TCPML.

\begin{overpic}[width=214.63977pt]{cml_dns_vort400} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_dns_vort1400}
\put(-2.0,75.0){$(b)$} \end{overpic}

Figure 6:  Contours in magnitude of vorticity at $(a)$ $\tau_{ml}=400$ and $(b)$ $\tau_{ml}=1400$ for DNS. Case | $\dot{\delta}_{\theta}/\dot{\delta}_{inc}$ | $Re_{\omega}$ | $Ma_{t}$ | $(b_{11},b_{12},b_{22})$
---|---|---|---|---
Ref1 | 0.675 | 7790 | - | (0.15, 0.15, -0.10)
Ref2 | 0.589 | 8160 | 0.30 | (0.13, 0.13, -0.12)
DNS | 0.588 | 8052 | 0.28 | (0.14, 0.16, -0.07)
Table 3:  Key quantities for DNS in TCPML at $\tau_{ml}=1400$.

\begin{overpic}[width=214.63977pt]{cml_dns_theta} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_dns_bij}
\put(-2.0,75.0){$(b)$} \end{overpic}

Figure 7:  Time history of $(a)$ normalized momentum thickness
$\delta_{\theta}/\delta_{\theta 0}$ and $(b)$ evolution of anisotropy stress
tensor $b_{ij}$ for DNS.

\begin{overpic}[width=214.63977pt]{cml_dns_r11} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_dns_r22}
\put(-2.0,75.0){$(b)$} \end{overpic}

Figure 8:  Profiles of normalized turbulent stress $(a)$
$(R_{11})^{1/2}/\Delta U$ and $(a)$ $(R_{22})^{1/2}/\Delta U$ for experiment
(Elliott & Samimy, 1990), Ref1 (Pantano & Sarkar, 2002) and DNS.

Table 2 shows the numerical parameters for TCPML. Ref1 was simulated by high-
order finite difference method (Pantano & Sarkar, 2002) with a smaller
$Ma_{c}=0.70$ and Ref2 was simulated by WENO-enhanced GKS (Arun et al., 2019).
The size effect of computational domain and grid convergence studies have been
well studied in Ref2. Compared with the WENO-enhanced GKS (Kumar et al.,
2013), the current in-house code is equipped with the genuine spatial-temporal
HGKS (Pan et al., 2016b; Cao et al., 2022). Table 2 shows that the same
computational domain adopted and much finer grid are used by the current HGKS,
which definitely guarantee the resolution of DNS. Figure 6 shows the contours
in magnitude of vorticity $\|\omega_{mg}\|=\sqrt{2\omega_{i}\omega_{i}}$ at
$\tau_{ml}=400$ and $\tau_{ml}=1400$. Against the performance of
$\|\omega_{mg}\|$ at transitional stage ($\tau_{ml}=400$, see figure 6$(a)$),
the magnitude of vorticity not only enlarges thicker but also behaves more
intermittently at the self-similarity stage ($\tau_{ml}=1400$, see figure
6$(b)$).

It is well known that compressibility suppresses the mixing layer growth rate
$\dot{\delta}_{\theta}=\text{d}\delta_{\theta}/\text{d}\tau_{ml}$. In table 3,
the normalized growth rate $\dot{\delta}_{\theta}/\dot{\delta}_{inc}=0.589$
agrees well with that in Ref2 (Arun et al., 2019), and reasonably smaller than
that of Ref1 (smaller $Ma_{c}=0.70$ corresponding to larger normalized
momentum thickness) (Pantano & Sarkar, 2002). Incompressible growth rate
$\dot{\delta}_{inc}=0.016$ is chosen for the hyperbolic tangent profile as
equation (41). Growth rate in this paper is computed by the least-square
method within $\tau_{ml}\in[1000,1300]$ as shown in figure 7$(a)$. Table 3
also presents the vorticity thickness-based Reynolds number $Re_{\omega}$ and
the turbulent Mach number $Ma_{t}$ at $\tau_{ml}=1400$. Here, turbulent Mach
number is defined as $Ma_{t}^{2}=2K_{R}/(\gamma RT_{\infty})$. Table 3 shows
that $Re_{\omega}$ and $Ma_{t}$ at the center plane are in good agreement for
all cases. (In Ref1 and Ref2, the ending of the simulation may be
$\tau_{ml}=600$). Figure 7$(b)$ shows the time history of anisotropy stress
tensor $b_{ij}$ as equation (46). We observe the the well-matched quasi-
stationary profiles during the self-similarity stage. Table 3 shows the
components of anisotropy stress tensor $(b_{11},b_{12},b_{22})$ at
$\tau_{ml}=1400$. The large deviation in $b_{22}$ can be attributed to the
differences in setting up the initial perturbation field (Arun et al., 2019).
More specifically, figure 8 shows the profiles of normalized turbulent stress
$(R_{11})^{1/2}/\Delta U$ and $(R_{22})^{1/2}/\Delta U$ (see equation (47)).
Refereed solutions suggest the envelop for normalized turbulent stress in
TCPML correspond to $Ma_{c}=0.75$. The reasonable deviation originates from
the different convective Mach number, where the data of Ref1 and experiment
(Elliott & Samimy, 1990) correspond to $Ma_{c}=0.70$ and $Ma_{c}=0.64$,
respectively. Overall, the current DNS results agree well with refereed
numerical simulations. After obtaining the high-fidelity flow fields from DNS,
M1 with NTRKM, SM and DSM will be conducted on unresolved grids subsequently.

Table 2 shows that M1 with NTRKM, SM and DSM are conducted on unresolved
uniform grids $144\times 96\times 48$. Box filter is used to generate the
initial (restarted) six-variable flow field (see equations 10-13) from the DNS
solution at $\tau_{ml}=400$, i.e., $4^{3}$ resolved grids are coarsen to $1$
unresolved grid. The same computational domain and boundary conditions are
applied as the DNS. In terms of SM, to guarantee the numerical stability and
improve the dissipative behavior, Vreman et al. (1997) recommended
$C_{sm}=0.1$ for compressible mixing layer. When implementing NTRKM, the
minimum unresolved TKE is set as $\left\langle K_{utke0}\right\rangle/10000$,
where the initial ensemble unresolved $\left\langle
K_{utke0}\right\rangle=0.00058$. Figure 9$(a)$ shows the initial pointwise
unresolved TKE. The initial unresolved $K_{utke}$ is in a small magnitude and
restricted in a narrow region. From transitional stage to self-similarity
stage (at $\tau_{ml}=1400$), figure 9$(b)$ shows that the magnitude of
unresolved $K_{utke}$ increase obviously, as well as entrain to a much wider
region similar as the figure 6$(b)$. Again, the intrinsic equilibrium
assumption on $K_{utke}$, such as $\left\langle S_{t}\right\rangle\approx 0$
for zero-equation eddy-viscosity LES models (Lilly, 1967; Germano et al.,
1991; Moin et al., 1991) may not hold.

Figure 10$(a)$ shows the time history of ensemble unresolved $\left\langle
K_{utke}\right\rangle$ from NTRKM and ensemble resolved kinetic energy
$\left\langle K\right\rangle$. The ensemble unresolved TKE from NTRKM
increases, while the ensemble resolved kinetic energy $\left\langle
K\right\rangle$ decrease in the dissipative system. Figure 10$(b)$ also shows
the evolution of normalized momentum thickness $\delta_{\theta}/\delta_{\theta
0}$. Figure 10 shows that the performance of statistical quantities from NTRKM
is much closer with that from DSM. In table 4, the normalized growth rate
$\dot{\delta}_{\theta}/\dot{\delta}_{inc}=0.611$ from NTRKM agrees well with
that from DSM, slightly larger than that of DNS. However, SM overestimates the
normalized growth rate up to $20\%$. At the center plane, table 4 also shows
the vorticity thickness-based Reynolds number $Re_{\omega}$, and the turbulent
Mach

\begin{overpic}[width=214.63977pt]{cml_ntrkm_ksgs400} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_ntrkm_ksgs1400}
\put(-2.0,75.0){$(b)$} \end{overpic}

Figure 9:  Contours of unresolved $K_{utke}$ at $(a)$ $\tau_{ml}=400$ and $(b)$ $\tau_{ml}=1400$ for case M1 with NTRKM. Case M1 | $\dot{\delta}_{\theta}/\dot{\delta}_{inc}$ | $Re_{\omega}$ | $Ma_{t}$ | $(b_{11},b_{12},b_{22})$
---|---|---|---|---
NTRKM | 0.611 | 11912 | 0.27 | (0.13, 0.16, -0.08)
SM | 0.694 | 10106 | 0.29 | (0.12, 0.16, -0.04)
DSM | 0.619 | 11814 | 0.30 | (0.15, 0.15, -0.08)
Table 4:  Key quantities for LES in TCPML at $\tau_{ml}=1400$.

\begin{overpic}[width=214.63977pt]{cml_les_ek} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_les_theta}
\put(-2.0,75.0){$(b)$} \end{overpic}

Figure 10:  Time history of $(a)$ ensemble unresolved $\left\langle
K_{utke}\right\rangle$ (only from NTRKM) and ensemble resolved kinetic energy
$\left\langle K\right\rangle$, and $(b)$ normalized momentum thickness
$\delta_{\theta}/\delta_{\theta 0}$ for case M1.

number $Ma_{t}$ from LES at $\tau_{ml}=1400$. The vorticity thickness-based
Reynolds numbers $Re_{\omega}$ based on three LES models are larger than that
of DNS, where the turbulent Mach number $Ma_{t}$ from LES agrees well with
that of DNS. Components of anisotropy stress tensor at $\tau_{ml}=1400$ agree
well with each other in $b_{11}$ and $b_{12}$ (see equation (46)). The large
deviation appears in $b_{22}$, as NTRKM solution coincides with the DSM
solution. Compared with DNS results in table 3, NTRKM and DSM outperform the
SM on $b_{22}$.

Figure 11 presents the profiles of normalized turbulent stress
$(R_{11})^{1/2}/\Delta U$, $(R_{22})^{1/2}/\Delta U$, $(R_{33})^{1/2}/\Delta
U$, and $(R_{12})^{1/2}/\Delta U$ (see equation (47)). The normalized
turbulent stress from three models shows quite small deviations. Against the
DNS solution as shown in figure 8, the peak magnitude of turbulent stress as
$(R_{11})^{1/2}/\Delta U$ and $(R_{22})^{1/2}/\Delta U$ are smaller in LES
simulations, indicating that LES on unresolved grids underestimate the
turbulent fluctuations. In TCPML, the performance of key turbulent quantities
up to second-order statistics confirm that current NTRKM is comparable with
the widely-used eddy-viscosity models. The results of NTRKM and DSM are much
closer to DNS solution, and outweighing the SM guided with recommended
modeling coefficient.

\begin{overpic}[width=214.63977pt]{cml_les_r11} \put(-2.0,75.0){$(a)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_les_r22}
\put(-2.0,75.0){$(b)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_les_r33}
\put(-2.0,75.0){$(c)$}
\end{overpic}\begin{overpic}[width=214.63977pt]{cml_les_r12}
\put(-2.0,75.0){$(d)$} \end{overpic}

Figure 11:  Profiles of normalized turbulent stress $(a)$
$(R_{11})^{1/2}/\Delta U$, $(b)$ $(R_{22})^{1/2}/\Delta U$, $(c)$
$(R_{33})^{1/2}/\Delta U$ and $(d)$ $(R_{12})^{1/2}/\Delta U$ for case M1.

##  5 Concluding remarks

We propose the non-equilibrium time-relaxation kinetic model for compressible
turbulence modeling for the first time. The key idea is constructing the non-
trivial turbulent quantities and its corresponding dynamic evolution
quantitatively on unresolved grids. Within the non-equilibrium time-relaxation
framework, NTRKM is constructed in the form of modified BGK model on
unresolved grids. Based on the $1$st-order Chapman-Enskog expansion, NTRKM
connects with the six-variable macroscopic governing equations.
Phenomenologically, the unknown turbulent relaxation time and source term in
NTRKM are determined by gradient-type assumption and dynamic modeling
approach. Therefore, the non-equilibrium kinetic model provides an profound
mesoscopic understanding for transport equation of the compressible SGS
turbulence kinetic energy. To solve the NTRKM accurately and robustly, finite
volume non-equilibrium gas-kinetic scheme is developed in the spirt of well-
established gas-kinetic scheme. DCIT and TCPML are simulated as benchmarks to
evaluate current non-equilibrium kinetic model and non-equilibrium gas-kinetic
scheme. The performance of key statistical turbulent quantities up to second-
order statistics confirms that current NTRKM is comparable with the widely-
used eddy-viscosity models. As expected, the performance of NTRKM is much
closer with DSM and outperforming SM. Present work not only points an
alternative way for compressible turbulence modeling on unresolved grids, but
also opens the great possibilities to simulate multi-scale flow physics within
the non-equilibrium time-relaxation framework.

NTRKM will be further implemented for the more practical compressible
turbulent flows to validate its strength, such as compressible wall-bounded
turbulent flows and shock-boundary layer interaction (Chen et al., 2012,
2017). Compared with incompressible turbulence, compressible turbulent flows
are more complex due to the non-linear coupling of velocity, density and
pressure fields. Consequently, compressible LES models are much difficult to
construct than the incompressible ones. Be of scientific interest, the optimal
physical distribution of $K_{utke}$ and modeling the source term directly as
equation (8) on unresolved grids, still requires the continuous effort in
understanding the non-equilibrium properties of turbulence. To that extent, a
complete time-relaxation compressible turbulence modeling without any ad-hoc
technique from macroscopic turbulence models may be achieved. In addition, if
NTRKM can be solved by multi-scale kinetic solver, such as UGKS (Xu & Huang,
2010), the non-equilibrium multi-scale fluxes may improve its performance
correspondingly. These challenging topics deserve to be explored in the
subsequent studies.

Acknowledgements. Thank Z.Y. Wang and Dr. W. Xu at HKUST for helpful
discussions. The authors would like to thank TaiYi supercomputers in the
SUSTech for providing high performance computational resources.

Funding. This research is supported by the National Natural Science Foundation
of China (91852114, 11772281 and 11701038), the National Numerical Windtunnel
project, the Department of Science and Technology of Guangdong Province
(2020B1212030001), and the Fundamental Research Funds for the Central
Universities.

Declaration of interests. The authors report no conflict of interest.

##  Appendix A Connection between non-equilibrium time-relaxation kinetic
model and macro governing equations

This appendix provides the details for the derivation of corresponding
macroscopic governing equations on unresolved grids based on the NTRKM.
Derivation of the Euler equations and the NS equations from the BGK model can
be found in the Appendix B (Xu, 2015). Similar as molecular relaxation time
$\tau=\epsilon\hat{\tau}$ in refereed derivation, the turbulent relaxation
time is rewritten as $\tau_{t}=\epsilon\hat{\tau_{t}}$ in which $\epsilon$ is
a small dimensionless quantity. Suppose that $f^{eq}$ has a Taylor series
expansion about point ($\bm{x}$, $t$). Since $\tau$ and $\tau_{t}$ depends on
the local thermodynamic variables and unresolved $K_{utke}$, and these depend
on the moments of $f^{eq}$, we may assume that $\tau$ and $\tau_{t}$ are
consequently $\hat{\tau}$ and $\hat{\tau_{t}}$ can be expanded about the point
($\bm{x}$, $t$). Now consider the formal solution of the NTRKM for
$f(\bm{x},t,\bm{u},\xi,k_{u})$, supposing that $f^{eq}$ is known, it can be
shown that $f(\bm{x},t,\bm{u},\xi,k_{u})$ has an expansion in powers of
$\epsilon$.

We derive the terms in this expansion from the formal solution of
$f(\bm{x},t,\bm{u},\xi,k_{u})$, by putting the power expansion
$f=f_{0}+\epsilon f_{1}+\epsilon^{2}f_{2}+\cdots$, and
$(\tau+\tau_{t})=\epsilon(\hat{\tau}+\hat{\tau_{t}})$ into the NTRKM directly.
Let $D_{\bm{u}}=\partial/\partial t+u_{i}\partial/\partial x_{i}$, and write
the NTRKM as $\epsilon(\hat{\tau}+\hat{\tau_{t}})D_{\bm{u}}f+f-f^{eq}=0$, with
handling the source term $Q_{s}$ in equation (4) splitly. An expansion of this
equation in powers of $\epsilon$ yields

$f=g-\epsilon(\hat{\tau}+\hat{\tau_{t}})D_{\bm{u}}f^{eq}+\epsilon^{2}(\hat{\tau}+\hat{\tau_{t}})D_{\bm{u}}((\hat{\tau}+\hat{\tau_{t}})D_{\bm{u}}f^{eq})+\cdots.$
(49)

According to compatibility condition, after dividing by
$\epsilon(\hat{\tau}+\hat{\tau_{t}})$, gives

$\displaystyle\int\bm{\psi}D_{\bm{u}}f^{eq}\text{d}\Xi=\epsilon\int\bm{\psi}D_{\bm{u}}[(\hat{\tau}+\hat{\tau_{t}})D_{\bm{u}}f^{eq}]\text{d}\Xi+\mathcal{O}{(\epsilon^{2})},$
(50)

where
$\bm{\psi}=(1,u_{1},u_{2},u_{3},\displaystyle\frac{1}{2}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+\xi^{2})+k_{u},k_{u})^{T}$,
and
$\text{d}\Xi=\text{d}u_{1}\text{d}u_{2}\text{d}u_{3}\text{d}\xi\text{d}k_{u}$.
The integral on the left-hand-side $\mathcal{L}_{\alpha}$ and the right-hand-
side $\mathcal{R}_{\alpha}$ of equation (50) can be defined as

$\displaystyle\mathcal{L}_{\alpha}=\epsilon\mathcal{R}_{\alpha}+\mathcal{O}{(\epsilon^{2})},$
(51)

which shows that $\mathcal{L}_{\alpha}$ is at least $\mathcal{O}{(\epsilon)}$,
obviously. Therefore, in reducing the $\mathcal{R}_{\alpha}$ on the right side
of equation (51), which is already $\mathcal{O}{(\epsilon)}$, we can drop
$\mathcal{O}{(\epsilon)}$ quantities and their derives. To simplify the
notation, let

$\displaystyle<\psi_{\alpha}(\cdot)>\equiv\int\psi_{\alpha}f^{eq}\text{d}\Xi,\quad\alpha=1,2,3,4,5,6,$
(52)

where $<\psi_{\alpha}(\cdot)>$ denotes the moments of $f^{eq}$ on
$\psi_{\alpha}$, and $\psi_{\alpha}$ is the component of $\bm{\psi}$.
Therefore, $\mathcal{L}_{\alpha}$ and $\mathcal{R}_{\alpha}$ are rewritten as

$\displaystyle\mathcal{L}_{\alpha}$
$\displaystyle=<\psi_{\alpha}>_{,t}+<\psi_{\alpha}u_{l}>_{,l},$ (53)
$\displaystyle\mathcal{R}_{\alpha}$
$\displaystyle=\\{(\hat{\tau}+\hat{\tau_{t}})[<\psi_{\alpha}u_{k}>_{,t}+<\psi_{\alpha}u_{k}u_{l}>_{,l}]\\}_{,k}+\mathcal{O}{(\epsilon)},$

where $k$ and $l$ are subscript of molecular velocity, taking from $1$ to $3$
for three-dimension derivation. The macroscopic governing equations in Euler-
type and NS-type can be obtained by truncating the equation (51) up to the
order of $\mathcal{O}{(1)}$ and $\mathcal{O}{(\epsilon)}$, respectively. The
Euler-type macroscopic equations can be derived straightforwardly by the
Chapman-Enskog expansion up to zeroth order, which are used to simplify the
derivation of NS-type equations. While, this appendix focuses on the NS-type
macroscopic equations by the Chapman-Enskog expansion up to first order
truncation of $\tau+\tau_{t}$.

###  A.1 Continuity equation

Continuity equation is derived straightforwardly as

$\displaystyle\rho_{,t}+(\rho U_{k})_{,k}=0,$ (54)

which can be used to simplify the momentum equations, the total energy
equations, and the unresolved TKE equation.

###  A.2 Momentum equation

To simplify the time derivative of pressure in the momentum equations, we
introduce the following procedure firstly. For total energy equation, the left
side $\mathcal{L}_{5}$ of equation (51) can be grouped as

$\displaystyle\mathcal{L}_{5}$
$\displaystyle=\frac{1}{2}U_{n}^{2}[{\rho_{,t}+(\rho U_{k})_{,k}}]+\rho
U_{n}U_{n,t}+\rho U_{k}U_{n}U_{n,k}+U_{k}p_{,k}$ (55)
$\displaystyle+\frac{N+3}{2}[p_{,t}+U_{k}p_{,k}]+\frac{N+5}{2}pU_{k,k}+{(\rho
K_{utke})_{,t}+(\rho K_{utke}U_{k})_{,k}},$

where $N$ is the total number of degrees of freedom in $\xi$. The first term
is $U_{n}^{2}\mathcal{L}_{1}/2$ which is $\mathcal{O}(\epsilon^{2})$, and next
three are $U_{n}\mathcal{L}_{n}$, and are therefore $\mathcal{O}(\epsilon)$.
Then $\mathcal{L}_{5}$ can be rewritten as

$\displaystyle\mathcal{L}_{5}$
$\displaystyle=\frac{N+3}{2}[p_{,t}+U_{k}p_{,k}]+\frac{N+5}{2}pU_{k,k}+U_{n}\mathcal{L}_{n}$
(56) $\displaystyle+{(\rho K_{utke})_{,t}+(\rho
K_{utke}U_{k})_{,k}}+\mathcal{O}{(\epsilon^{2})}.$

Based on the Chapman-Enskog expansion up to $0$th order, unresolved TKE
equation $\mathcal{L}_{6}$ in Euler-type can be written as

$\displaystyle\mathcal{L}_{6}=(\rho K_{utke})_{,t}+(\rho
K_{utke}U_{k})_{,k}-S_{t}+\mathcal{O}{(\epsilon)},$ (57)

which can be used to simplify the time derivative of $(\rho K_{utke})_{,t}$ in
the following derivation. Combining equation (56) and equation (57), we get

$\displaystyle
p_{,t}+U_{k}p_{,k}=-\frac{N+5}{N+3}pU_{k,k}-\frac{2S_{t}}{N+3}+\mathcal{O}(\epsilon),$
(58)

which can be used to simplify the time derivative of pressure $p_{,t}$ in the
following derivation.

For the right hand sides of the momentum equations in equation (51),
considering $\mathcal{R}_{j}=[(\hat{\tau}+\hat{\tau_{t}})F_{jk}]_{,k}$, we get

$\displaystyle F_{jk}$
$\displaystyle=<u_{j}u_{k}>_{,t}+<u_{j}u_{k}u_{l}>_{,l}$ (59)
$\displaystyle=U_{j}[{{(\rho U_{k})_{,t}+[(\rho
U_{k}U_{l})+p\delta_{kl}]_{,l}}}]+\rho U_{k}U_{j,t}+(p\delta_{jk})_{,t}$
$\displaystyle+(\rho
U_{k}U_{l}+p\delta_{kl})U_{j,l}+(U_{l}p\delta_{jk}+U_{k}p\delta_{jl})_{,l},$

where the fact that all odd moments in $w_{k}$ vanish has been used. Here
$w_{k}=u_{k}-U_{k}$ is the peculiar velocity. The term in square brackets
multiplying $U_{j}$ is $\mathcal{L}_{k}$, i.e. it is $\mathcal{O}(\epsilon)$,
and can therefore be ignored. Then, after gathering terms with coefficients
$U_{k}$ and $p$, we have

$\displaystyle F_{jk}=U_{k}[{{\rho U_{j,t}+\rho
U_{l}U_{j,l}+p_{,j}}}]+p[U_{k,j}+U_{j,k}+U_{l,l}\delta_{jk}]+\delta_{jk}[{{p_{,t}+U_{l}p_{,l}}}]+\mathcal{O}(\epsilon).$
(60)

The coefficient of $U_{k}$ can be simplified by the following equation

$\displaystyle\mathcal{L}_{i}=\rho U_{i,t}+\rho
U_{k}U_{i,k}+p_{,i}+\mathcal{O}(\epsilon^{2}),$ (61)

where equation (61) is derived by multiplying the continuity equation by
$U_{i}$ and subtracting the result from $\mathcal{L}_{i}$ ($i=2,3,4$). To
eliminate $p_{,t}$ from the last term we use the equation (58) for
$\mathcal{L}_{5}$. Finally, decompose the tensor $U_{k,j}$ into its dilation
and shear parts in the usual way, which gives

$\displaystyle
F_{jk}=p[U_{k,j}+U_{j,k}-\frac{2}{3}U_{l,l}\delta_{jk}]+\frac{2}{3}\frac{N}{(N+3)}pU_{l,l}\delta_{jk}-\frac{2S_{t}}{N+3}\delta_{jk}+\mathcal{O}(\epsilon).$
(62)

The second term is due to bulk viscosity involves energy sharing between
transnational and internal degrees of freedom of the molecular (Xu, 2015;
Cramer, 2012), and the last term resulted from the energy interaction between
the resolved kinetic energy and the unresolved TKE.

###  A.3 Total energy equation

Analogy to derive the NS total energy equation, we write
$\mathcal{R}_{5}=\\{(\hat{\tau}+\hat{\tau_{t}})N_{k}\\}_{,k}$ with

$\displaystyle
N_{k}=<u_{k}(\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u})>_{,t}+<u_{k}u_{l}(\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u})>_{,l},$
(63)

where the $k_{u}$ is the sample-space variable corresponding to unresolved
$K_{utke}$. $N_{k}$ can be decomposed into $N_{k}=N_{k}^{(1)}+N_{k}^{(2)}$,
where

$\displaystyle
N_{k}^{(1)}=[U_{k}(\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u})]_{,t}+[U_{k}<u_{l}(\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u})>]_{,l},$
(64)

and

$\displaystyle
N_{k}^{(2)}=<w_{k}(\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u})>_{,t}+<w_{k}u_{l}(\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u})>_{,l}.$
(65)

For $N_{k}^{(1)}$, we have

$\displaystyle N_{k}^{(1)}$
$\displaystyle=U_{k}[{{<\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u}>_{,t}+<u_{l}(\frac{u_{n}^{2}+\xi^{2}}{2}+k_{u})>_{,l}}}]$
(66) $\displaystyle+[\frac{1}{2}\rho U_{n}^{2}+\frac{N+3}{2}p+\rho
K_{utke}]U_{k,t}+[U_{l}(\frac{1}{2}\rho U_{n}^{2}+\frac{N+5}{2}p+\rho
K_{utke})]U_{k,l}.$

The coefficient of $U_{k}$ in the equation above is $\mathcal{L}_{5}$, and
therefore can be dropped, and the remaining terms can be rewritten as

$\displaystyle N_{k}^{(1)}=[\frac{1}{2}\rho U_{n}^{2}+\frac{N+3}{2}p+\rho
K_{utke}][U_{k,t}+U_{l}U_{k,l}]+pU_{l}U_{k,l}+\mathcal{O}(\epsilon).$ (67)

According to equation (61) to replace $U_{k,t}$, we get

$\displaystyle N_{k}^{(1)}$
$\displaystyle=-[\frac{1}{2}U_{n}^{2}+\frac{N+3}{2}\frac{p}{\rho}+K_{utke}]p_{,k}+pU_{l}U_{k,l}+\mathcal{O}(\epsilon).$
(68)

For $N_{k}^{(2)}$, remembering that moments odd in $w_{k}$ vanish, we have

$\displaystyle N_{k}^{(2)}$
$\displaystyle=<U_{n}w_{n}w_{k}>_{,t}+<U_{l}U_{n}w_{n}w_{k}>_{,l}$ (69)
$\displaystyle+\frac{1}{2}<U_{n}^{2}w_{k}w_{l}>_{,l}+<w_{k}w_{l}(\frac{w_{n}^{2}+\xi^{2}}{2}+k_{u})>_{,l}$
$\displaystyle=(pU_{k})_{,t}+(pU_{k}U_{l})_{,l}+\frac{1}{2}(U_{n}^{2}p)_{,k}+\frac{N+5}{2}(\frac{p^{2}}{\rho})_{,k}+(pK_{utke})_{,k}.$

$N_{k}^{(2)}$ can be rewritten as

$\displaystyle N_{k}^{(2)}$
$\displaystyle=p[U_{k,t}+U_{l}U_{k,l}+U_{k}U_{l,l}+U_{l}U_{l,k}]$ (70)
$\displaystyle+U_{k}(p_{,t}+U_{l}p_{,l})+\frac{1}{2}U_{n}^{2}p_{,k}+\frac{N+5}{2}(\frac{p^{2}}{\rho})_{,k}+(pK_{utke})_{,k}.$

The $p_{,t}$ and $U_{k,t}$ can be replaced by equation (58) and equation (61),
respectively. Hence

$\displaystyle N_{k}^{(2)}$
$\displaystyle=p[-\frac{p_{,k}}{\rho}+U_{k}U_{l,l}+U_{l}U_{l,k}]+U_{k}[{{-\frac{N+5}{N+3}pU_{l,l}-\frac{2S_{t}}{N+3}}}]$
(71)
$\displaystyle+\frac{1}{2}U_{n}^{2}p_{,k}+\frac{N+5}{2}(\frac{p^{2}}{\rho})_{,k}+(pK_{utke})_{,k}+\mathcal{O}(\epsilon).$

Finally, $N_{k}$ can be obtained by summing $N_{k}^{(1)}$ and $N_{k}^{(2)}$ up

$\displaystyle N_{k}$
$\displaystyle=p[U_{l}(U_{k,l}+U_{l,k})-\frac{2}{N+3}U_{k}U_{l,l}]-U_{k}\frac{2S_{t}}{N+3}$
(72)
$\displaystyle+{\frac{N+5}{2}p(\frac{p}{\rho})_{,k}+(pK_{utke})_{,k}}+\mathcal{O}(\epsilon).$

###  A.4 Unresolved turbulence kinetic energy equation

For unresolved $K_{utke}$ equation, we write
$\mathcal{R}_{6}=\\{(\hat{\tau}+\hat{\tau_{t}})Z\\}_{,k}$, where

$\displaystyle Z$ $\displaystyle=<ku_{k}>_{,t}+<ku_{k}u_{l}>_{,l}$ (73)
$\displaystyle=K_{utke}[(\rho U_{k})_{,t}+(\rho
U_{k}U_{l}+p\delta_{kl})_{,l}]+\rho U_{k}K_{utke,t}+K_{utke,l}[\rho
U_{k}U_{l}+p\delta_{kl}].$

The term in square brackets is $\mathcal{L}_{k}$, i.e.
$\mathcal{O}(\epsilon)$, and can be dropped. Equation (57) subtracts the
multiplication of the continuity equation (54) by $K_{utke}$ gives

$\displaystyle\mathcal{L}_{6}=\rho K_{utke,t}+\rho
U_{l}K_{utke,l}-S_{t}+\mathcal{O}(\epsilon^{2}).$ (74)

Gathering terms with coefficients $U_{k}$ and $p$, and replacing $K_{utke,t}$
through equation (74), we have

$\displaystyle Z=U_{k}S_{t}+pK_{utke,k}+\mathcal{O}(\epsilon).$ (75)

All time derivatives have now been removed from $\mathcal{R}_{\alpha}$, and
the remaining steps in deriving corresponding macroscopic governing equations
for NTRKM may be summarized briefly as

  * •

Drop $\mathcal{O}(\epsilon^{2})$ in equation (51),

  * •

Combine $\epsilon$ and $\hat{\tau}+\hat{\tau_{t}}$ to recover
$\tau+\tau_{t}=\epsilon(\hat{\tau}+\hat{\tau_{t}})$,

  * •

Define the molecular dynamic viscosity as $\mu=\tau p$, and turbulent eddy
viscosity is recovered by $\mu_{t}=\tau_{t}p$,

  * •

Define the molecular thermal conductivity $\kappa=(N+5)\tau p/2$, and the
turbulent thermal conductivity $\kappa_{t}=(N+5)\tau_{t}p/2$.

Finally, corresponding macroscopic governing equations to NTRKM can be
rewritten as equations (10)-(13) in §2.3.

##  Appendix B Dynamic approach to determine modeling coefficients

Model coefficients $C_{s}$, $Pr_{t}$, and $C_{\Pi}$ can be dynamically
computed through Germano identity (Germano et al., 1991; Lilly, 1992), which
assumes the similarity of SGS quantities between the grid filter width
$\overline{\Delta}$ and the test filter width $\widehat{\overline{\Delta}}$.
In finite volume framework, explicit filter is not used, while the grid length
of control volume itself acts as the grid filter width, and the projection
process when updating the macroscopic variables can be regarded as the
filtering process. For any term
$\phi=\overline{\beta_{1}\beta_{2}}-\overline{\beta}_{1}\overline{\beta}_{2}$
on grid filter level, assuming that
$\Phi=\widehat{\overline{\beta_{1}\beta_{2}}}-\widehat{\overline{\beta}}_{1}\widehat{\overline{\beta}}_{2}$
still holds on test filter level. Then, the resolved tensor (or vector/scalar)
is defined as $L=\Phi-\widehat{\phi}$. Assume $\phi$ is modeled by the linear
constitutive relationship $\phi=Cm$, where $m$ is a function of the resolved
quantities. At the test filter level, $\Phi=CM$, $M$ takes similar form to $m$
but is a function of the test-filtered quantities. Plugging the linear model
for $\Phi$ and $\phi$, the Germano identity reads

$\displaystyle
L=\widehat{\overline{\beta}_{1}\overline{\beta}_{2}}-\widehat{\overline{\beta}}_{1}\widehat{\overline{\beta}}_{2}=C(M-\widehat{m}).$
(76)

To avoid computational instability, $C$ can be optimized by the least-square
method (Lilly, 1992).

###  B.1 Dynamic coefficient for density-weighted Smagorinsky model

For dynamic density-weighted SM (Moin et al., 1991), the model coefficient
$C_{dsm}$ is determined by

$\left\\{\begin{aligned} L_{ij}^{*}&=-2C_{dsm}M_{ij},\\\
M_{ij}&=\widehat{\overline{\Delta}}^{2}\widehat{\overline{\rho}}|\widehat{\widetilde{S}}|\widehat{\widetilde{S}}^{*}_{ij}-{\overline{\Delta}}^{2}\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{\overline{\rho}|\widetilde{S}|\widetilde{S}^{*}_{ij}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
\overline{\rho}|\widetilde{S}|\widetilde{S}^{*}_{ij}\\\
\rule{-4.30554pt}{0.0pt}\end{array},\\\
C_{dsm}&=-\frac{L_{ij}^{*}M_{ij}}{2{M_{ij}M_{ij}}},\end{aligned}\right.$ (77)

where
$L_{ij}=\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{{{\overline{\rho{U_{i}}}\,\overline{\rho{U_{j}}}}}/{\overline{\rho}}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
{{\overline{\rho{U_{i}}}\,\overline{\rho{U_{j}}}}}/{\overline{\rho}}\\\
\rule{-4.30554pt}{0.0pt}\end{array}-\widehat{\overline{\rho{U_{i}}}}\,\widehat{\overline{\rho{U_{j}}}}/{\widehat{\overline{\rho}}}$,
$L_{ij}^{*}=L_{ij}-L_{kk}\delta_{ij}/3$,
$\widetilde{S}_{ij}=(\widetilde{U}_{i,j}+\widetilde{U}_{j,i})/2$,
$\widetilde{S}_{ij}^{\ast}=\widetilde{S}_{ij}-\widetilde{S}_{kk}\delta_{ij}/3$,
$|\widetilde{S}|=(2\widetilde{S}_{ij}\widetilde{S}_{ij})^{\frac{1}{2}}$,
$|\widehat{\widetilde{S}}|=(2\widehat{\widetilde{S}}_{ij}\widehat{\widetilde{S}}_{ij})^{\frac{1}{2}}$.

###  B.2 Dynamic coefficient for non-equilibrium time-relaxation kinetic model

In NTRKM, $\tau_{ij}^{sgs}$ is modeled by the gradient-type eddy viscosity
model as equation (22), and the coefficient $C_{s}$ is determined by

$\left\\{\begin{aligned} L_{ij}^{*}&=-2C_{s}M_{ij},\\\
M_{ij}&=\widehat{\overline{\Delta}}\widehat{\overline{\rho}}\widehat{\overline{K}}_{utke}^{\frac{1}{2}}\widehat{\widetilde{S}}^{*}_{ij}-{\overline{\Delta}}\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{\overline{\rho}K_{utke}^{\frac{1}{2}}\widetilde{S}^{*}_{ij}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
\overline{\rho}K_{utke}^{\frac{1}{2}}\widetilde{S}^{*}_{ij}\\\
\rule{-4.30554pt}{0.0pt}\end{array},\\\
C_{s}&=-\frac{L_{ij}^{*}M_{ij}}{2{M_{ij}M_{ij}}},\end{aligned}\right.$ (78)

where
$L_{ij}=\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{{{\overline{\rho{U_{i}}}\,\overline{\rho{U_{j}}}}}/{\overline{\rho}}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
{{\overline{\rho{U_{i}}}\,\overline{\rho{U_{j}}}}}/{\overline{\rho}}\\\
\rule{-4.30554pt}{0.0pt}\end{array}-\widehat{\overline{\rho{U_{i}}}}\,\widehat{\overline{\rho{U_{j}}}}/{\widehat{\overline{\rho}}}$,
$L_{ij}^{*}=L_{ij}-L_{kk}\delta_{ij}/3$, and
$\widehat{\overline{K}}_{utke}=\widehat{K}_{utke}+(\widehat{\widetilde{U}_{k}\widetilde{U}_{k}}-\widehat{\widetilde{U}}_{k}\widehat{\widetilde{U}}_{k})/2$.
When modifying the intrinsic turbulent Prandtl number $Pr_{t}=1$ in NTRKM, the
realistic turbulent Prandtl number can be computed dynamically as

$\left\\{\begin{aligned} L_{j}&=-C_{s}M_{j}/Pr_{t},\\\
M_{j}&=\widehat{\overline{\Delta}}\widehat{\overline{\rho}}\widehat{\overline{K}}_{utke}^{\frac{1}{2}}\widehat{\widetilde{T}}_{,j}-\overline{\Delta}\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{\overline{\rho}K_{utke}^{\frac{1}{2}}\widetilde{T}_{,j}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
\overline{\rho}K_{utke}^{\frac{1}{2}}\widetilde{T}_{,j}\\\
\rule{-4.30554pt}{0.0pt}\end{array},\\\
Pr_{t}&=-\frac{C_{s}M_{j}M_{j}}{L_{j}M_{j}},\end{aligned}\right.$ (79)

where
$L_{j}=R(\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{{\overline{{\rho}T}\,\overline{{\rho}U_{j}}}/{\overline{\rho}}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
{\overline{{\rho}T}\,\overline{{\rho}U_{j}}}/{\overline{\rho}}\\\
\rule{-4.30554pt}{0.0pt}\end{array}-\widehat{\overline{{\rho}T}}\,\widehat{\overline{{\rho}U_{j}}}/{\widehat{\overline{\rho}}})$.
To model the pressure-dilation transfer, following the series expansion (Chai
& Mahesh, 2012), $C_{\Pi}$ reads

$\left\\{\begin{aligned} L&=C_{\Pi}M,\\\
M&=\widehat{\overline{\Delta}}^{2}\widehat{\overline{p}}_{,j}(\widehat{\overline{\rho
U_{k}}}/{\widehat{\overline{\rho}}})_{j,k}-\overline{\Delta}^{2}\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{\overline{p}_{,j}(\widetilde{U}_{k})_{j,k}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
\overline{p}_{,j}(\widetilde{U}_{k})_{j,k}\\\
\rule{-4.30554pt}{0.0pt}\end{array},\\\
C_{\Pi}&=\frac{L}{M},\end{aligned}\right.$ (80)

where
$L=\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{\overline{p}(\overline{\rho
U_{k}}/{\overline{\rho}})_{,k}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
\overline{p}(\overline{\rho U_{k}}/{\overline{\rho}})_{,k}\\\
\rule{-4.30554pt}{0.0pt}\end{array}-\widehat{\overline{p}}(\widehat{\overline{\rho
U_{k}}}/\widehat{\overline{\rho}})_{,k}$.

The coefficients $C_{\epsilon s}$ and $C_{\epsilon d}$ are obtained by the
analogy between the grid-filter-level SGS stress and the resolved stress
across the test filter level (Menon & Kim, 1996), where $L=CM$ instead of
$L=C(M-\widehat{m})$. To model the solenoidal dissipation rate, $C_{\epsilon
s}$ is obtained by $C_{\epsilon s}=L/M$, where

$\left\\{\begin{aligned}
L&=\widehat{\overline{\mu}\widetilde{w}_{i}\widetilde{w}_{i}}-\widehat{\overline{\mu}}\widehat{\widetilde{w}}_{i}\widehat{\widetilde{w}}_{i},\\\
M&=\widehat{\overline{\rho}}\widehat{\overline{K}}_{utke}^{\frac{3}{2}}/{\widehat{\overline{\Delta}}}.\end{aligned}\right.$
(81)

To model the dilational dissipation rate, $C_{\epsilon d}$ is obtained by
$C_{\epsilon d}=L/M$, where

$\left\\{\begin{aligned}
L&=4[\begin{array}[]{c}\stretchto{\scaleto{\scalerel*[\widthof{\overline{\mu}\widetilde{U}_{k,k}^{2}}]{\kern-0.5pt\bigwedge\kern-0.5pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}\\\
\overline{\mu}\widetilde{U}_{k,k}^{2}\\\
\rule{-4.30554pt}{0.0pt}\end{array}-\widehat{\overline{\mu}}\widehat{\widetilde{U}}_{k,k}^{2}]/3,\\\
M&=\widehat{\overline{\rho}}^{2}\widehat{\overline{K}}_{utke}^{\frac{5}{2}}/({\gamma\widehat{\overline{p}}\,\widehat{\overline{\Delta}}}).\end{aligned}\right.$
(82)

###### Remark 1

For DCIT, the global dynamic coefficient is obtained through the ensemble
average in the whole computational domain. In terms of TCPML, the local
dynamic coefficient is obtained through the plane averaging in the streamwise
and spanwise directions.

###### Remark 2

When modeling SGS pressure-dilation transfer, the dynamic denominator in
$C_{\Pi}$ may approach to a pretty tiny value, which cause a large value for
model coefficient $C_{\Pi}$. To smoothen the spurious behavior of $C_{\Pi}$,
the $C_{\Pi}$ is limited in a bound $[-0.1,0.1]$ to guarantee the numerical
robustness. When $C_{\Pi}$ goes beyond the bound, it would be modified as
$C_{\Pi}=C_{\Pi}^{lim}\bm{\cdot}sign(C_{\Pi})$, where $C_{\Pi}^{lim}=0.05$ and
$sign(\cdot)$ is the sign function. Both for DCIT and TCPML, numerical tests
indicate that the key statistical quantities are not sensitive to the
$C_{\Pi}^{lim}$, as long as the numerical simulations are stable.

###### Remark 3

For macroscopic varialbe $\beta$, the filtered macroscopic variable and Favre-
filtered variable are denoted as $\overline{\beta}$ and $\tilde{\beta}$ in
equations (18) - (26) and Appendix B, respectively. By default, macroscopic
variable $\beta$ in the rest part of this paper, such as the derivation of
macroscopic governing equations from NTRKM in Appendix A represents resolved
value without additional filtering symbols. Since resolved conservative
variables on unresolved grids can be regarded as the filtered conservative
variables, the corresponding Favre-filtered primitive variables also can be
obtained. In summary, the resolved variables in current paper can be treated
as the customary filtered variables equivalently.

##  Appendix C Connection between macroscopic variables and mesoscopic
coefficients

The connection between the spatial derivatives of macroscopic flow variables
and the expansion of turbulence equilibrium distribution function $f^{eq}$
reads

$\frac{\partial\bm{Q}}{\partial x_{i}}=\int\frac{\partial f^{eq}}{\partial
x_{i}}\bm{\psi}\text{d}\Xi\equiv\int af^{eq}\bm{\psi}\text{d}\Xi,$ (83)

where $a$ denotes the spatial mesoscopic coefficients in equation (33) as

$a=\bm{a}^{T}\bm{\psi}=a_{1}+a_{2}u_{1}+a_{3}u_{2}+a_{4}u_{3}+a_{5}[\frac{1}{2}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+\xi^{2})+k_{u}]+a_{6}k_{u}.$
(84)

Equation (83) can be rewritten into following linear system

$\frac{1}{\rho}\frac{\partial\bm{Q}}{\partial
x_{i}}=\Big{(}\frac{1}{\rho}\int\bm{\psi}\otimes\bm{\psi}^{T}f^{eq}\text{d}\Xi\Big{)}\bm{a}\triangleq\bm{M}\bm{a},$
(85)

Each component of $(a_{1},…,a_{6})^{T}$ in equation (85) can be determined
uniquely

$\left\\{\begin{aligned} a_{6}&=\frac{4}{\rho K_{utke}^{2}}B_{5}-a_{5},\\\
a_{5}&=\frac{8\lambda^{2}}{\rho(N+3)}(B_{4}-U_{1}B_{1}-U_{2}B_{2}-U_{3}B_{3}-B_{5}),\\\
a_{4}&=\frac{2\lambda}{\rho}B_{3}-U_{3}a_{5},\\\
a_{3}&=\frac{2\lambda}{\rho}B_{2}-U_{2}a_{5},\\\
a_{2}&=\frac{2\lambda}{\rho}B_{1}-U_{1}a_{5},\\\
a_{1}&=\frac{1}{\rho}\frac{\partial\rho}{\partial
x_{i}}-U_{1}a_{2}-U_{2}a_{3}-U_{3}a_{4}-Ea_{5}-K_{utke}a_{6},\end{aligned}\right.$
(86)

where

$\left\\{\begin{aligned} B_{1}&=\frac{\partial(\rho U_{1})}{\partial
x_{i}}-U_{1}\frac{\partial\rho}{\partial x_{i}},\\\ B_{2}&=\frac{\partial(\rho
U_{2})}{\partial x_{i}}-U_{2}\frac{\partial\rho}{\partial x_{i}},\\\
B_{3}&=\frac{\partial(\rho U_{3})}{\partial
x_{i}}-U_{3}\frac{\partial\rho}{\partial x_{i}},\\\ B_{4}&=\frac{\partial(\rho
E)}{\partial x_{i}}-E\frac{\partial\rho}{\partial x_{i}},\\\
B_{5}&=\frac{\partial(\rho K_{utke})}{\partial
x_{i}}-K_{utke}\frac{\partial\rho}{\partial x_{i}}.\end{aligned}\right.$ (87)

For the temporal mesoscopic coefficient in equation (33), the relation between
temporal derivatives of macroscopic variables and turbulence equilibrium
distribution can be written as

$\frac{\partial\bm{Q}}{\partial t}=\int\frac{\partial f^{eq}}{\partial
t}\bm{\psi}\text{d}\Xi\equiv\int Af^{eq}\bm{\psi}\text{d}\Xi,$ (88)

where

$A=\bm{A}^{T}\bm{\psi}=A_{1}+A_{2}u_{1}+A_{3}u_{2}+A_{4}u_{3}+A_{5}[\frac{1}{2}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+\xi^{2})+k_{u}]+A_{6}k_{u}.$
(89)

The temporal derivatives of macroscopic variables can be given according to
the compatibility condition as

$\displaystyle\int(\frac{\partial{f^{eq}}}{\partial t}+u_{i}\frac{\partial
f^{eq}}{\partial x_{i}})\bm{\psi}\text{d}\Xi=\bm{0}.$ (90)

In a similar way, the above components $(A_{1},…,A_{6})^{T}$ in equation (88)
can be determined uniquely.

##  Appendix D Eigenstructure for characteristic reconstruction

The Jacobian matrix $\bm{J}$ for the flux $\bm{F}(\bm{Q})$ of the hyperbolic
part in equations (10)-(13) is given by

$\displaystyle\bm{J}=\begin{pmatrix}0&1&0&0&0&0\\\
\frac{\hat{\gamma}}{2}\bm{U^{2}}-U_{1}^{2}&(3-\gamma)U_{1}&-\hat{\gamma}U_{2}&-\hat{\gamma}U_{3}&\hat{\gamma}&-\hat{\gamma}\\\
-U_{1}U_{2}&U_{2}&U_{1}&0&0&0\\\ -U_{1}U_{3}&U_{3}&0&U_{1}&0&0\\\
U_{1}(\frac{\gamma-2}{2}\bm{U^{2}}-\frac{c^{2}}{\gamma-1}-K_{utke})&H-\hat{\gamma}U_{1}^{2}+K_{utke}&-\hat{\gamma}U_{1}U_{2}&-\hat{\gamma}U_{1}U_{3}&\gamma
U_{1}&-\hat{\gamma}U_{1}\\\ -U_{1}K_{utke}&K_{utke}&0&0&0&U_{1}\\\
\end{pmatrix},$ (91)

where $H=\bm{U^{2}}/2+c^{2}/(\gamma-1)$, and $\hat{\gamma}=\gamma-1$. The
eigenvalues of the quasi one-dimensional system are

$\lambda_{1}=U_{1}-c,\quad\lambda_{2}=U_{1},\quad\lambda_{3}=U_{1},\quad\lambda_{4}=U_{1},\quad\lambda_{5}=U_{1},\quad\lambda_{6}=U_{1}+c,$
(92)

where $c$ is the local sound speed. The matrix of corresponding right
eigenmatrix is

$\displaystyle\bm{R}=\begin{pmatrix}1&1&0&0&0&1\\\
U_{1}-c&U_{1}&0&0&0&U_{1}+c\\\ U_{2}&U_{2}&1&0&0&U_{2}\\\
U_{3}&U_{3}&0&1&0&U_{3}\\\
H-U_{1}c+K_{utke}&\frac{1}{2}\bm{U^{2}}&U_{2}&U_{3}&1&H+U_{1}c+K_{utke}\\\
K_{utke}&0&0&0&1&K_{utke}\\\ \end{pmatrix},$ (93)

and the inverse matrix of $\bm{R}$ is given by

$\displaystyle\bm{R}^{-1}=\frac{\hat{\gamma}}{2c^{2}}\begin{pmatrix}\frac{1}{2}\bm{U^{2}}+\frac{U_{1}c}{\hat{\gamma}}&-(U_{1}+\frac{c}{\hat{\gamma}})&-U_{2}&-U_{3}&1&-1\\\
-\bm{U^{2}}+\frac{2c^{2}}{\hat{\gamma}}&2U_{1}&2U_{2}&2U_{3}&-2&2\\\
-\frac{2U_{2}c^{2}}{\hat{\gamma}}&0&\frac{2c^{2}}{\hat{\gamma}}&0&0&0\\\
-\frac{2U_{3}c^{2}}{\hat{\gamma}}&0&0&\frac{2c^{2}}{\hat{\gamma}}&0&0\\\
-\bm{U^{2}}K_{utke}&2U_{1}K_{utke}&2U_{2}K_{utke}&2U_{3}K_{utke}&-2K_{utke}&2K_{utke}+\frac{2c^{2}}{\hat{\gamma}}\\\
\frac{1}{2}\bm{U^{2}}-\frac{U_{1}c}{\hat{\gamma}}&-U_{1}+\frac{c}{\hat{\gamma}}&-U_{2}&-U_{3}&1&-1\end{pmatrix}.$
(94)

## References

  * Arun et al. (2019) Arun, S, Sameen, A, Srinivasan, B & Girimaji, SS 2019 Topology-based characterization of compressibility effects in mixing layers. Journal of Fluid Mechanics 874, 38–75.
  * Bhatnagar et al. (1954) Bhatnagar, Prabhu Lal, Gross, Eugene P & Krook, Max 1954 A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems. Physical review 94 (3), 511.
  * Bogey & Bailly (2004) Bogey, Christophe & Bailly, Christophe 2004 A family of low dispersive and low dissipative explicit schemes for flow and noise computations. Journal of Computational physics 194 (1), 194–214.
  * Cao et al. (2018) Cao, Guiyu, Liu, Hualin & Xu, Kun 2018 Physical modeling and numerical studies of three-dimensional non-equilibrium multi-temperature flows. Physics of Fluids 30 (12), 126104\.
  * Cao et al. (2019a) Cao, Guiyu, Pan, Liang & Xu, Kun 2019a Three dimensional high-order gas-kinetic scheme for supersonic isotropic turbulence i: criterion for direct numerical simulation. Computers & Fluids 192, 104273.
  * Cao et al. (2021) Cao, Guiyu, Pan, Liang & Xu, Kun 2021 Three dimensional high-order gas-kinetic scheme for supersonic isotropic turbulence ii: coarse-graining analysis of compressible ksgs budget. Journal of Computational Physics p. 110402.
  * Cao et al. (2022) Cao, Guiyu, Pan, Liang & Xu, Kun 2022 High-order gas-kinetic scheme with parallel computation for direct numerical simulation of turbulent flows. Journal of Computational Physics 448, 110739.
  * Cao et al. (2019b) Cao, Guiyu, Su, Hongmin, Xu, Jinxiu & Xu, Kun 2019b Implicit high-order gas kinetic scheme for turbulence simulation. Aerospace Science and Technology 92, 958–971.
  * Castro et al. (2011) Castro, Marcos, Costa, Bruno & Don, Wai Sun 2011 High order weighted essentially non-oscillatory weno-z schemes for hyperbolic conservation laws. Journal of Computational Physics 230 (5), 1766–1792.
  * Chai & Mahesh (2012) Chai, Xiaochuan & Mahesh, Krishnan 2012 Dynamic-equation model for large-eddy simulation of compressible flows. Journal of Fluid Mechanics 699, 385–413.
  * Chapman & Cowling (1970) Chapman, Sydney & Cowling, TG 1970 The mathematical theory of non-uniform gases: An account of the kinetic theory of viscosity, thermal conduction and diffusion in gases. cambridge mathematical library. Cambridge University Press 1, 27–52.
  * Chen et al. (2003) Chen, Hudong, Kandasamy, Satheesh, Orszag, Steven, Shock, Rick, Succi, Sauro & Yakhot, Victor 2003 Extended boltzmann kinetic equation for turbulent flows. Science 301 (5633), 633–636.
  * Chen et al. (2017) Chen, Liang, Xiao, Zuoli, Shi, Yipeng & Chen, Shiyi 2017 Constrained large-eddy simulation of supersonic turbulent boundary layer over a compression ramp. Journal of Turbulence 18 (8), 781–808.
  * Chen et al. (2012) Chen, Shiyi, Xia, Zhenhua, Pei, Suyang, Wang, Jianchun, Yang, Yantao, Xiao, Zuoli & Shi, Yipeng 2012 Reynolds-stress-constrained large-eddy simulation of wall-bounded turbulent flows. Journal of Fluid Mechanics 703, 1–28.
  * Courant et al. (1928) Courant, Richard, Friedrichs, Kurt & Lewy, Hans 1928 Über die partiellen differenzengleichungen der mathematischen physik. Mathematische annalen 100 (1), 32–74.
  * Cramer (2012) Cramer, Mark S 2012 Numerical estimates for the bulk viscosity of ideal gases. Physics of fluids 24 (6), 066102\.
  * De Stefano et al. (2008) De Stefano, Giuliano, Vasilyev, Oleg V & Goldstein, Daniel E 2008 Localized dynamic kinetic-energy-based models for stochastic coherent adaptive large eddy simulation. Physics of Fluids 20 (4), 045102.
  * Deardorff (1970) Deardorff, James W 1970 A numerical study of three-dimensional turbulent channel flow at large reynolds numbers. Journal of Fluid Mechanics 41 (2), 453–480.
  * Elliott & Samimy (1990) Elliott, Gregory S & Samimy, M 1990 Compressibility effects in free shear layers. Physics of Fluids A: Fluid Dynamics 2 (7), 1231–1240.
  * Fu & Li (2006) Fu, Song & Li, Qibing 2006 Numerical simulation of compressible mixing layers. International journal of heat and fluid flow 27 (5), 895–901.
  * Garnier et al. (2009) Garnier, Eric, Adams, Nikolaus & Sagaut, Pierre 2009 Large eddy simulation for compressible flows. Springer Science & Business Media.
  * Germano et al. (1991) Germano, Massimo, Piomelli, Ugo, Moin, Parviz & Cabot, William H 1991 A dynamic subgrid-scale eddy viscosity model. Physics of Fluids A: Fluid Dynamics 3 (7), 1760–1765.
  * Hou et al. (1996) Hou, Shuling, Sterling, J, Chen, Shiyi & Doolen, GD 1996 A lattice boltzmann subgrid model for high reynolds number flows, fields inst. In Comm. Citeseer.
  * Jiang & Qian (2012) Jiang, Jin & Qian, Yuehong 2012 Implicit gas-kinetic bgk scheme with multigrid for 3d stationary transonic high-reynolds number flows. Computers & Fluids 66, 21–28.
  * Kim et al. (1987) Kim, John, Moin, Parviz & Moser, Robert 1987 Turbulence statistics in fully developed channel flow at low reynolds number. Journal of fluid mechanics 177, 133–166.
  * Krajnović & Davidson (2002) Krajnović, Siniša & Davidson, Lars 2002 A mixed one-equation subgrid model for large-eddy simulation. International journal of heat and fluid flow 23 (4), 413–425.
  * Kumar et al. (2013) Kumar, G, Girimaji, Sharath S & Kerimo, J 2013 Weno-enhanced gas-kinetic scheme for direct simulations of compressible transition and turbulence. Journal of Computational Physics 234, 499–523.
  * Li & Du (2016) Li, Jiequan & Du, Zhifang 2016 A two-stage fourth order time-accurate discretization for lax–wendroff type flow solvers i. hyperbolic conservation laws. SIAM Journal on Scientific Computing 38 (5), A3046–A3069.
  * Liao et al. (2009) Liao, Wei, Peng, Yan & Luo, Li-Shi 2009 Gas-kinetic schemes for direct numerical simulations of compressible homogeneous turbulence. Physical Review E 80 (4), 046702\.
  * Lilly (1967) Lilly, Douglas K 1967 The representation of small-scale turbulence in numerical simulation experiments. IBM Form pp. 195–210.
  * Lilly (1992) Lilly, Douglas K 1992 A proposed modification of the germano subgrid-scale closure method. Physics of Fluids A: Fluid Dynamics 4 (3), 633–635.
  * Manabe et al. (1965) Manabe, Syukuro, Smagorinsky, Joseph & Strickler, Robert F 1965 Simulated climatology of a general circulation model with a hydrologic cycle. Mon. Wea. Rev 93 (12), 769–798.
  * Meneveau et al. (1996) Meneveau, Charles, Lund, Thomas S & Cabot, William H 1996 A lagrangian dynamic subgrid-scale model of turbulence. Journal of fluid mechanics 319, 353–385.
  * Menon & Kim (1996) Menon, S & Kim, W-W 1996 High reynolds number flow simulations using the localized dynamic subgrid-scale model. In 34th Aerospace Sciences Meeting and Exhibit, p. 425.
  * Moin et al. (1991) Moin, Parviz, Squires, Kyle, Cabot, W & Lee, Sangsan 1991 A dynamic subgrid-scale model for compressible turbulence and scalar transport. Physics of Fluids A: Fluid Dynamics 3 (11), 2746–2757.
  * Moser et al. (2021) Moser, Robert D, Haering, Sigfried W & Yalla, Gopal R 2021 Statistical properties of subgrid-scale turbulence models. Annual Review of Fluid Mechanics 53, 255–286.
  * Nicoud & Ducros (1999) Nicoud, Franck & Ducros, Frédéric 1999 Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, turbulence and Combustion 62 (3), 183–200.
  * Nicoud et al. (2011) Nicoud, Franck, Toda, Hubert Baya, Cabrit, Olivier, Bose, Sanjeeb & Lee, Jungil 2011 Using singular values to build a subgrid-scale model for large eddy simulations. Physics of fluids 23 (8), 085106.
  * Pan et al. (2016a) Pan, Dongxin, Zhong, Chengwen, Li, Ji & Zhuo, Congshan 2016a A gas-kinetic scheme for the simulation of turbulent flows on unstructured meshes. International Journal for Numerical Methods in Fluids 82 (11), 748–769.
  * Pan & Xu (2020) Pan, Liang & Xu, Kun 2020 High-order gas-kinetic scheme with three-dimensional weno reconstruction for the euler and navier-stokes solutions. Computers & Fluids 198, 104401\.
  * Pan et al. (2016b) Pan, Liang, Xu, Kun, Li, Qibing & Li, Jiequan 2016b An efficient and accurate two-stage fourth-order gas-kinetic scheme for the euler and navier–stokes equations. Journal of Computational Physics 326, 197–221.
  * Pantano & Sarkar (2002) Pantano, C & Sarkar, S 2002 A study of compressibility effects in the high-speed turbulent shear layer using direct simulation. Journal of Fluid Mechanics 451, 329.
  * Pomraning & Rutland (2002) Pomraning, Eric & Rutland, Christopher J 2002 Dynamic one-equation nonviscosity large-eddy simulation model. AIAA journal 40 (4), 689–701.
  * Pope (2001) Pope, Stephen B 2001 Turbulent flows.
  * Righi (2016) Righi, Marcello 2016 A gas-kinetic scheme for turbulent flow. Flow, Turbulence and Combustion 97 (1), 121–139.
  * Rozema et al. (2015) Rozema, Wybe, Bae, Hyun J, Moin, Parviz & Verstappen, Roel 2015 Minimum-dissipation models for large-eddy simulation. Physics of Fluids 27 (8), 085107.
  * Sagaut (2006) Sagaut, Pierre 2006 Large eddy simulation for incompressible flows: an introduction. Springer Science & Business Media.
  * Samtaney et al. (2001) Samtaney, Ravi, Pullin, Dale I & Kosović, Branko 2001 Direct numerical simulation of decaying compressible turbulence and shocklet statistics. Physics of Fluids 13 (5), 1415–1430.
  * Sandham & Reynolds (1991) Sandham, ND & Reynolds, WC 1991 Three-dimensional simulations of large eddies in the compressible mixing layer. Journal of Fluid Mechanics 224, 133–158.
  * Sandham (1990) Sandham, Neil David 1990 A numerical investigation of the compressible mixing layer. PhD thesis, Stanford University.
  * Sarkar et al. (1991) Sarkar, Sutanu, Erlebacher, Gordon, Hussaini, M Yousuff & Kreiss, Heinz Otto 1991 The analysis and modelling of dilatational terms in compressible turbulence. Journal of Fluid Mechanics 227, 473–493.
  * Schumann (1975) Schumann, Ulrich 1975 Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli. Journal of computational physics 18, 376–404.
  * Tan et al. (2018) Tan, Shuang, Li, Qibing, Xiao, Zhixiang & Fu, Song 2018 Gas kinetic scheme for turbulence simulation. Aerospace Science and Technology 78, 214–227.
  * Toro (2013) Toro, Eleuterio F 2013 Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media.
  * Vreman (2004) Vreman, AW 2004 An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications. Physics of fluids 16 (10), 3670–3681.
  * Vreman et al. (1994) Vreman, Bert, Geurts, Bernard & Kuerten, Hans 1994 Realizability conditions for the turbulent stress tensor in large-eddy simulation. Journal of Fluid Mechanics 278, 351–362.
  * Vreman et al. (1997) Vreman, Bert, Geurts, Bernard & Kuerten, Hans 1997 Large-eddy simulation of the turbulent mixing layer. Journal of fluid mechanics 339, 357–390.
  * Wang et al. (2013) Wang, Jianchun, Yang, Yantao, Shi, Yipeng, Xiao, Zuoli, He, XT & Chen, Shiyi 2013 Cascade of kinetic energy in three-dimensional compressible turbulence. Physical review letters 110 (21), 214505.
  * White & Majdalani (2006) White, Frank M & Majdalani, Joseph 2006 Viscous fluid flow, , vol. 3. McGraw-Hill New York.
  * Wilcox et al. (1998) Wilcox, David C & others 1998 Turbulence modeling for CFD, , vol. 2. DCW industries La Canada, CA.
  * Xu (2001) Xu, Kun 2001 A gas-kinetic bgk scheme for the navier–stokes equations and its connection with artificial dissipation and godunov method. Journal of Computational Physics 171 (1), 289–335.
  * Xu (2015) Xu, Kun 2015 Direct modeling for computational fluid dynamics: construction and application of unified gas-kinetic schemes. World Scientific.
  * Xu et al. (2008) Xu, Kun, He, Xin & Cai, Chunpei 2008 Multiple temperature kinetic model and gas-kinetic method for hypersonic non-equilibrium flow computations. Journal of computational physics 227 (14), 6779–6794.
  * Xu & Huang (2010) Xu, Kun & Huang, Juan-Chen 2010 A unified gas-kinetic scheme for continuum and rarefied flows. Journal of Computational Physics 229 (20), 7747–7764.
  * Yoshizawa (1986) Yoshizawa, Akira 1986 Statistical theory for compressible turbulent shear flows, with the application to subgrid modeling. The Physics of fluids 29 (7), 2152–2164.
  * Yoshizawa & Horiuti (1985) Yoshizawa, Akira & Horiuti, Kiyosi 1985 A statistically-derived subgrid-scale kinetic energy model for the large-eddy simulation of turbulent flows. Journal of the Physical Society of Japan 54 (8), 2834–2839.