# Modeling and simulation in supersonic three-temperature carbon dioxide
turbulent channel flow

Guiyu Cao<EMAIL_ADDRESS>Yipeng Shi<EMAIL_ADDRESS>Kun Xu
<EMAIL_ADDRESS>Shiyi Chen<EMAIL_ADDRESS>Academy for Advanced
Interdisciplinary Studies, Southern University of Science and Technology,
Shenzhen, Guangdong 518055, People’s Republic of China  Guangdong-Hong Kong-
Macao Joint Laboratory for Data-Driven Fluid Mechanics and Engineering
Applications, Southern University of Science and Technology, Shenzhen,
Guangdong 518055, People’s Republic of China  Department of Aeronautics &
Astronautics Engineering, Peking University, Beijing 100871, People’s Republic
of China  Department of Mathematics, Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, Hong Kong

###### Abstract

This paper pioneers the direct numerical simulation (DNS) and physical
analysis in supersonic three-temperature carbon dioxide ($CO_{2}$) turbulent
channel flow. $CO_{2}$ is a linear and symmetric triatomic molecular, with the
thermal non-equilibrium three-temperature effects arising from the
interactions among translational, rotational and vibrational modes under room
temperature. Thus, the rotational and vibrational modes of $CO_{2}$ are
addressed. Thermal non-equilibrium effect of $CO_{2}$ has been modeled in an
extended three-temperature BGK-type model, with the calibrated translational,
rotational and vibrational relaxation time. To solve the extended BGK-type
equation accurately and robustly, non-equilibrium high-accuracy gas-kinetic
scheme is proposed within the well-established two-stage fourth-order
framework. Compared with the one-temperature supersonic turbulent channel
flow, supersonic three-temperature $CO_{2}$ turbulence enlarges the ensemble
heat transfer of the wall by approximate $20\%$, and slightly decreases the
ensemble frictional force. The ensemble density and temperature fields are
greatly affected, and there is little change in Van Driest transformation of
streamwise velocity. The thermal non-equilibrium three-temperature effects of
$CO_{2}$ also suppress the peak of normalized root-mean-square of density and
temperature, normalized turbulent intensities and Reynolds stress. The
vibrational modes of $CO_{2}$ behave quite differently with rotational and
translational modes. Compared with the vibrational temperature fields, the
rotational temperature fields have the higher similarity with translational
temperature fields, especially in temperature amplitude. Current thermal non-
equilibrium models, high-accuracy DNS and physical analysis in supersonic
$CO_{2}$ turbulent flow can act as the benchmark for the long-term
applicability of compressible $CO_{2}$ turbulence.

###### keywords:

Carbon dioxide flow, Vibrational modes, Three-temperature effects, Supersonic
turbulent channel flows

##  1 Introduction

Mars exploration programs are currently experiencing a revival, such as
amazing Mars robotic helicopter ”Ingenuity” operating on Mars [1]. Mars’s
atmosphere consists of $95.32\%$ carbon dioxide ($CO_{2}$). For accurate
predictions of surface drag and heat flux on Martian vehicles, it is necessary
to take the peculiarities of $CO_{2}$ into account. Different with the
dominant diatomic gases nitrogen ($N_{2}$) and oxygen ($O_{2}$) on earth,
$CO_{2}$ is a linear and symmetric triatomic molecular, which has three
vibrational modes [2]. Triatomic molecular $CO_{2}$ is equipped with the
inherent thermal non-equilibrium multi-temperature effects arising from the
interactions among the translational, rotational and vibrational modes [3, 4].

Carbon dioxide is widely studied in the applications of physical chemistry and
fluid dynamics, i.e., $CO_{2}$-$N_{2}$ gas laser system [5], environmental
green-house problems [6], and Mars entry vehicles [7, 8, 9]. Complex molecular
structure and multiple internal energy relaxation mechanisms in $CO2$
significantly affect its physical and chemical properties [10]. In fluids
community, of special interest is the evaluation of bulk viscosity in $CO2$.
Stokes’ viscosity relation does not hold for $CO_{2}$ [11], as the bulk
viscosity can be thousands of times larger than the shear viscosity. Many
research works are engaged in experimental and theoretical studies in the bulk
viscosity of $CO_{2}$ [12, 13]. Recent experiments [14, 15] show that the
vibrational modes contribute dominantly to the bulk viscosity of $CO_{2}$, and
the bulk viscosity from rotational modes is only one half of its shear
viscosity approximately [15]. Thus, the vibrational modes should be modelled
carefully when simulating $CO_{2}$ flows. In view of the importance of
multiple internal energy modes, the multiple internal energy relaxation
mechanisms of $CO_{2}$ [3, 4] have been modeled and analyzed, which confirm
that equilibrium one-temperature gas flow description is not valid for
$CO_{2}$ flows even under room temperature (i.e., $300K$). To the author’s
knowledge, the thermal non-equilibrium physical models considering the multi-
temperature effects of $CO_{2}$ and its applications in turbulent flows are
seldom reported. For accurate predictions of $CO_{2}$ turbulence, it is
necessary to take the three-temperature effects of $CO_{2}$ into account.

In the past few decades, the gas-kinetic scheme (GKS) based on the Bhatnagar-
Gross-Krook (BGK) model [16, 17] has been developed systematically for the
computations from low speed flows to hypersonic ones [18, 19]. Based on the
time-dependent flux solver, including generalized Riemann problem solver and
GKS [20, 21], a reliable two-stage fourth-order framework was provided for
developing the high-order GKS (HGKS) into fourth-order accuracy. With the
advantage of finite volume GKS and HGKS, they have been naturally implemented
as a direct numerical simulation (DNS) tool in simulating turbulent flows [22,
23, 24], especially for compressible turbulence [25, 26]. Aiming to conduct
the large-scale DNS, a parallel in-house computational platform of HGKS has
been developed in uniform grids and curvilinear grids [27, 28], with high
efficiency, fourth-order accuracy and super robustness. In addition, with the
discretization of particle velocity space, a unified gas-kinetic scheme (UGKS)
[29, 30] and unified gas-kinetic wave particle method (UGKWP) [31, 32, 33]
have been developed for multi-scale physical transport problems. The well-
developed HGKS and multi-scale UGKS/UGKWP provides the solid foundation for
thermal non-equilibrium multi-temperature modeling and simulation in $CO_{2}$
flows. The multi-scale modeling and numerical framework can be applied in
multi-scale $CO_{2}$ flows, i.e., the Mar’s re-entry vehicles from rarefied to
continuum regimes. As a starter, current study focuses on the supersonic
$CO_{2}$ turbulence in the continuum regime.

In this paper, the vibrational modes of $CO_{2}$ are addressed, and the
translational, rotational and vibrational relaxation time of $CO_{2}$ are
calibrated. The three-temperature effects of $CO_{2}$ are modeled in an
extended three-temperature BGK-type model within the well-established kinetic
framework [34, 35, 36]. To achieve high-order accuracy in space and time for
simulating the supersonic $CO_{2}$ turbulence, the non-equilibrium high-
accuracy GKS has been constructed with the second-order kinetic flux, fifth-
order WENO-Z reconstruction [37], and two-stage fourth-order time
discretization [21]. One-temperature supersonic turbulent channel flow [38,
39] is simulated firstly to validate the numerical set-up with bulk Mach
number $Ma=3$ and bulk Reynolds number $Re=4880$. Considering the
translational, rotational, and vibrational specific heats at constant volume,
one-temperature supersonic turbulent channel flow of thermally perfect gas has
been studied [40]. With implementing the non-equilibrium high-accuracy GKS in
the large-scale parallel in-house platform [27, 28], for the first time, the
DNS in supersonic three-temperature $CO_{2}$ turbulent channel flow is
conducted. Compared with the one-temperature supersonic turbulent channel
flow, the three-temperature effects of $CO_{2}$ are analyzed. Numerical
simulation confirms the thermal non-equilibrium three-temperature performance
of $CO_{2}$. Both the maximum ensemble temperature and normalized r.m.s.
temperature sort from high to low is translational temperature, rotational
temperature, and vibrational temperature. Compared with the vibrational
temperature fields, the rotational temperature fields has the higher
similarity with translational temperature fields both in temperature amplitude
and its structure.

For physical modeling and numerical simulation in supersonic three-temperature
$CO_{2}$ turbulent channel flow, this paper is organized as follows. Section 2
addresses the internal energy modes of $CO_{2}$. Extended thermal non-
equilibrium three-temperature BGK-type model and corresponding non-equilibrium
high-accuracy GKS for $CO_{2}$ are included in Section 3. Numerical examples
and discussions are presented in Section 4. The last section is the conclusion
and remarks.

##  2 Internal energy modes of carbon dioxide

Thermal non-equilibrium three-temperature effects of $CO_{2}$ mainly arise
from the interactions among internal energy modes [4]. This section addresses
the vibrational modes of $CO_{2}$, and focuses on the translational,
rotational and vibrational relaxation time for the extended three-temperature
BGK-type model.

###  2.1 Rotational and vibrational modes

Carbon dioxide is a linear and symmetric triatomic molecular, with rotational
and vibrational internal energy modes. The characteristic rotational
temperature is defined as $\theta_{r}=h_{P}^{2}/(8\pi^{2}k_{B}I)$, where
$h_{P}$ is the Planck constant, $k_{B}$ the Boltzmann constant, $I$ the
molecular moment of inertia. For carbon dioxide, $\theta_{r}=0.56K$ can be
obtained [2], while the $\theta_{r}$ for $N_{2}$ and $O_{2}$ is $2.88K$ and
$2.08K$, respectively. Under room temperature, it is well known that the
rotational degrees of freedom (d.o.f.) are assumed to be excited completely
for $N_{2}$ and $O_{2}$. Since $CO_{2}$ is with the smaller characteristic
rotational temperature, the rotational d.o.f. of $CO_{2}$ are regarded as
complete excitation in current study (i.e., $CO_{2}$ gas temperature above
$300K$).

Carbon dioxide is equipped with three vibrational modes as one symmetric
stretching mode $\nu_{1}$, one double degenerated bending mode $\nu_{2}$, and
one asymmetric stretching mode $\nu_{3}$. The characteristic vibrational
temperature reads $\theta_{v}=h_{P}\tilde{\nu}c_{L}/k_{B}$, where
$\tilde{\nu}$ is the characteristic wavenumber, and $c_{L}$ the speed of light
in the vacuum. In experimental studies of $CO_{2}$, infrared spectrum gives
the wavelength for corresponding vibrational modes $\nu_{2}$ and $\nu_{3}$,
and Raman spectroscopy provides the wavelength for $\nu_{1}$ [2]. Table 1
presents the characteristic wavenumber, wavelength ($\lambda=1/\tilde{\nu}$),
characteristic vibrational temperature and corresponding degeneracy of
$CO_{2}$. The double degenerated bending modes $\nu_{2}$ are most likely to be
activated with characteristic vibrational temperature $\theta_{v}=959.66K$,
which is much lower than that of $N_{2}$ with $\theta_{v}=3521K$ and $O_{2}$
with $\theta_{v}=2256K$. Thus, the vibrational modes of $N_{2}$ and $O_{2}$ is
usually considered in high-temperature applications, i.e., re-entry vehicles
experiencing the temperature above $800K$ [41]. However, the excitation of
vibrational modes of $CO_{2}$ requires to be modeled and simulated even under
the room temperature [4, 15].

Vibrational mode | Wavenumber($\tilde{\nu}$)/$cm^{-1}$ | Wavelength($\lambda$)/$\mu m$ |  $\theta_{v}$/$K$ | Degeneracy
---|---|---|---|---
$\nu_{1}$ | $1388$ | $7.20$ | $1997.02$ | $1$
$\nu_{2}$ | $667$ | $14.99$ | $959.66$ | $2$
$\nu_{3}$ | $2349$ | $4.26$ | $3379.69$ | $1$
Table 1:  Parameters for three vibrational modes of $CO_{2}$.

With the assumption that there is a unique vibrational temperature at each
point in the flow fields, the translational internal energy $E_{t}$,
rotational internal energy $E_{r}$, and vibrational internal energy $E_{v}$
per unit mass of $CO_{2}$ read

$\displaystyle E_{t}$ $\displaystyle=\frac{N_{t}}{2}RT_{t},$ (1)
$\displaystyle E_{r}$ $\displaystyle=\frac{N_{r}}{2}RT_{r},$ (2)
$\displaystyle E_{v}$
$\displaystyle=R\sum_{i=1}^{3}g_{i}\frac{\theta_{v,i}}{e^{\theta_{v,i}/T_{v}}-1},$
(3)

with translational d.o.f. as $N_{t}$, rotational d.o.f. as $N_{r}$, and
vibrational d.o.f. as $N_{v}$, where $T_{t}$, $T_{r}$ and $T_{v}$ represent
translational temperature, rotational temperature and vibrational temperature,
respectively. In Eq.(3), $\theta_{v,i}$ is the characteristic vibrational
temperature of vibrational mode $\nu_{i}$, and $g_{i}$ is the corresponding
degeneracy of mode $\nu_{i}$ as shown in Table 1. The total specific heat at
constant volume $C_{V}$ can be obtained by the sum of three components as

$\displaystyle C_{V}=C_{V,t}+C_{V,r}+C_{V,v},$ (4)

with

$\displaystyle C_{V,t}$ $\displaystyle=\frac{N_{t}}{2}R,$ (5) $\displaystyle
C_{V,r}$ $\displaystyle=\frac{N_{r}}{2}R,$ (6) $\displaystyle C_{V,v}$
$\displaystyle=\sum_{i=1}^{3}g_{i}\frac{(\theta_{v,i}/T_{v})^{2}}{e^{\theta_{v,i}/T_{v}}+e^{-\theta_{v,i}/T_{v}}-2}R,$
(7)

where $C_{V,t}$, $C_{V,r}$ and $C_{V,v}$ denotes the componential specific
heat at constant volume for translational internal energy, rotational internal
energy and vibrational internal energy, respectively. In current study, the
supersonic $CO_{2}$ turbulence is considered above $300K$, thus, translational
d.o.f. as $N_{t}=3$ and rotational d.o.f. as $N_{r}=2$ are adopted.
Vibrational d.o.f. as $N_{v}$ can be obtained by the definition as
$C_{V,v}=N_{v}R/2$, and the specific heat ratio $\gamma$ is the function of
vibrational d.o.f. as

$\displaystyle N_{v}$
$\displaystyle=2\sum_{i=1}^{3}g_{i}\frac{(\theta_{v,i}/T_{v})^{2}}{e^{\theta_{v,i}/T_{v}}+e^{-\theta_{v,i}/T_{v}}-2},$
(8) $\displaystyle\gamma$
$\displaystyle=\frac{C_{P}}{C_{V}}=1+\frac{2}{5+N_{v}}.$ (9)

In Eq.(9), $C_{P}$ is the total specific heat at constant pressure. In the
high-temperature limit, note that the vibrational d.o.f. as $N_{v}$ approaches
to classical definition with
$2\sum_{i=1}^{3}g_{i}\frac{\theta_{v,i}/T_{v}}{e^{\theta_{vib,i}/T_{v}}-1}$.

For compressible wall-bounded turbulence, the Prandtl number $Pr$ plays a key
role in determining the statistical turbulent quantities [27], especially for
density and temperature fields. For $CO_{2}$, the $Pr$ can be calculated as
$Pr=\mu C_{P}/\kappa$, where shear viscosity $\mu$ and thermal conductivity
$\kappa$ depends on the translational temperature $T_{t}$ as subsequent
Eq.(10) and Eq.(12), and $C_{P}$ relies on the vibrational temperature as
Eq.(9). Figure 1 presents the comparisons on $\gamma$ and $Pr$ between
$CO_{2}$ and air. For air, the specific heat ratio $\gamma$ almost keeps as
$1.4$ up to $800K$, and $Pr$ is fixed at approximate $0.7$ between $300K$ and
$800K$. In terms of $CO_{2}$, experimental measurements on $\gamma$ and $Pr$
show the strong temperature-dependent behavior [42]. We find the numerical
profile of $\gamma$ in $CO_{2}$ with Eq.(9) agrees well with the experimental
one. The numerical profile of $Pr$ is calculated with the assumption of
$T_{t}=T_{v}$, and this assumption dose not hold for experimental measurement.
Thus, it is reasonable to find the discrepancy in $Pr$ between the numerical
profile and experimental result. In following simulation, without the
assumption of $T_{t}=T_{v}$ as shown in Figure 1, the $Pr$ of $CO_{2}$ depends
on the practical translational and vibrational temperatures. Notice that the
$\gamma$ and $Pr$ must be computed locally in each time step for each
computational grid [34].

Figure 1:  Comparisons on the specific heat ration $\gamma$ (left) and Prandtl
number $Pr$ (right) for $CO_{2}$ and air. The numerical profile of $Pr$ is
calculated with the assumption $T_{t}=T_{v}$. Experimental results are adapted
from the Ref [42].

###  2.2 Translational, rotational and vibrational relaxation time

For $CO_{2}$ gas flows, non-equilibrium translational, rotational, and
vibrational internal modes are equipped with different relaxation time. Aiming
to the following construction of thermal non-equilibrium three-temperature
BGK-type model, the translational, rotational and vibrational relaxation time
requires to be determined firstly. Translational relaxation time $\tau_{t}$ is
related with shear viscosity $\mu$, and the power law [42] gives the
approximation of shear viscosity as

$\displaystyle\mu(T_{t})=\mu_{0}(\frac{T_{t}}{T_{0}})^{n_{t}},$ (10)

where $T_{t}$ is the translational temperature, $n_{t}=0.79$ and
$\mu_{0}=1.370\times 10^{-5}kg/(m\cdot s)$ at $T_{0}=273K$ for $CO_{2}$ is
approximately valid between $209K$ and $1700K$. While, $n_{t}=0.666$ and
$\mu_{0}=1.716\times 10^{-5}kg/m\cdot s$ for air is approximately valid
between $210K$ and $1900K$. From Chapman-Enskog expansion [17, 19], the
translational relaxation time $\tau_{t}$ of $CO_{2}$ can be obtained by
$\tau_{t}=\mu/p$, where $p$ is the pressure. Thus, the corresponding power law
of translational relaxation time $\tau_{t}$ reads

$\displaystyle\tau_{t}=\tau_{t0}(\frac{T_{t}}{T_{0}})^{n_{t}},$ (11)

where $\tau_{t0}=1.35\times 10^{-10}s$ at $T_{0}=273K$ with the atmosphere
pressure $101.325\times 10^{3}Pa$. The power law also gives the approximation
of thermal conductivity $\kappa$ [42] as

$\kappa(T_{t})=\kappa_{0}(\frac{T_{t}}{T_{0}})^{n_{\kappa}},$ (12)

where $n_{\kappa}=1.30$ and $\kappa_{0}=1.46\times 10^{-2}W/(m\cdot K)$ at
$T_{0}=273K$ for $CO_{2}$ is approximately valid between $180K$ and $700K$.
While, $n_{\kappa}=0.81$ and $\kappa_{0}=2.41\times 10^{-2}W/(m\cdot K)$ for
air is valid between $210K$ and $2000K$. The thermal conductivity in Eq.(12)
gets involved with the calculation of $Pr$ as shown in Figure 1 and following
numerical simulation in $CO_{2}$ turbulence. For comparison, the Sutherland
law for shear viscosity and thermal conductivity of $CO_{2}$ and air, as well
as the curve-fit expression in the shear viscosity of $CO_{2}$ [7] over a much
wider range of translational temperature has been provided in Appendix A.

The rotational relaxation time $\tau_{r}$ and vibrational relaxation time
$\tau_{v}$ are related with the bulk viscosity $\eta_{b}$ [13], which is given
by

$\eta_{b}=\eta_{b,c}+\eta_{b,int},$ (13)

where $\eta_{b,c}$ represents the effect of elastic collisions, and
$\eta_{b,int}$ accounts for the inelastic collisional contribution of the
internal d.o.f. [43]. For dilute gas, $\eta_{b,c}$ can be neglected, resulting
in

$\displaystyle\eta_{b}$
$\displaystyle\approx\eta_{b,int}=(\gamma-1)^{2}\sum_{i=1}^{N}\frac{C_{V,l}}{R}p\tau_{l},$
(14)

where $\gamma$ is the specific heat ratio as Eq.(9), $R$ the gas constant,
$C_{V,l}$ the specific heat at constant volume as Eq.(5) - Eq.(7), and
$\tau_{l}$ the relaxation time for the $l$th internal energy mode. Assuming
that the rotational and vibrational modes relax independently with a single
rotational relaxation time $\tau_{r}$ and vibrational relaxation time
$\tau_{v}$, the bulk viscosity as Eq.(14) can be rewritten as

$\displaystyle\eta_{b}$ $\displaystyle=\eta_{b}^{r}+\eta_{b}^{v},$ (15)

with

$\displaystyle\eta_{b}^{r}$
$\displaystyle=(\gamma-1)^{2}\frac{N_{r}}{2}p\tau_{r},$ (16)
$\displaystyle\eta_{b}^{v}$
$\displaystyle=(\gamma-1)^{2}\frac{C_{V,v}}{R}p\tau_{v}.$ (17)

$\eta_{b}^{r}$ denotes the componential bulk viscosity arises from the
rotational modes, and $\eta_{b}^{v}$ the componential bulk viscosity results
from the vibrational modes.

The power law of $CO_{2}$ for rational relaxation time is calibrated by
Eq.(16), with the combination of experimental data [15] and Parker equation
[44, 45]. As reported in experimental study [15], the vibrational modes remain
frozen when measuring the bulk viscosity of $CO_{2}$ with Rayleigh-Brillouin
light scattering spectroscopy at $532nm$. In such circumstance, the rotational
modes contribute solely to the bulk viscosity, and Eq.(15) reduces to
$\eta_{b}=\eta_{b}^{r}$. Thence, the seminal estimated experimental data of
$\eta_{b}^{r}$ using Hammond-Wiggings hydrodynamic model [15] is utilized to
calibrate the $\tau_{r}$. Power law of rotational relaxation time $\tau_{r}$
is given by

$\displaystyle\tau_{r}=\tau_{r0}(\frac{T_{t}}{T_{0}})^{n_{r}},$ (18)

where $n_{r}=1.59$ and $\tau_{r0}=2.99\times 10^{-10}s$ at $T_{0}=273K$ is
approximately valid between $250K$ and $2000K$ with the atmosphere pressure.
The details in calibrating the rotational relaxation time $\tau_{r}$ can be
find in Appendix B. With the well-known calibration of $p\tau_{v}$ [13], the
vibrational relaxation time can be calibrated as

$\displaystyle\tau_{v}=\tau_{v0}(\frac{T_{t}}{T_{0}})^{n_{v}},$ (19)

where $n_{v}=-1.353$ and $\tau_{v0}=1.89\times 10^{-6}s$ at $T_{0}=800K$ is
approximately valid between $300K$ and $1700K$ with the atmosphere pressure.

Figure 2:  Comparisons on translational relaxation time $\tau_{t}$, rotational
relaxation time $\tau_{r}$ and vibrational relaxation time $\tau_{v}$ for
$CO_{2}$ and air.

As power laws of Eq.(11), Eq.(18) and Eq.(19) present, all relaxation times
only depend on the translational temperature. Figure 2 presents the
comparisons on translational, rotational and vibrational relaxation time for
$CO_{2}$ and air. In Figure 2, the translational relaxation time of air is
slightly larger than that of $CO_{2}$. We clearly observe that the vibrational
relaxation time $\tau_{v}$ can be approximately thousands of times longer than
that of $\tau_{t}$ and $\tau_{r}$. Current quantitative calibration on
$\tau_{t}$, $\tau_{r}$ and $\tau_{v}$ deviates from the previous theoretical
calculations using kinetic theory [10], while the qualitative behaviors are
similar. It should be noted that the bulk viscosity of $CO_{2}$ is mainly
dominated by the $\eta_{b}^{v}$ arising from the vibrational modes [15]. Thus,
present $\tau_{v}$ indeed gives the large bulk viscosity with Eq.(17), which
is consistent with the findings of inherent large bulk viscosity of $CO_{2}$
[11, 43].

##  3 Three-temperature kinetic model and numerical scheme for carbon dioxide

In this section, an extended thermal non-equilibrium three-temperature BGK-
type model for $CO_{2}$ is introduced with well-calibrated relaxation time.
The three-temperature kinetic equation is going to be solved by the proposed
finite-volume non-equilibrium high-accuracy GKS.

###  3.1 Three-temperature BGK-type model

In terms of one-temperature equilibrium gas flow, the BGK model [16] has been
well proposed. For thermal non-equilibrium three-temperature diatomic gas
flow, the extended BGK-type model [34, 46] reads

$\displaystyle\frac{\partial f}{\partial t}+u_{i}\frac{\partial f}{\partial
x_{i}}=\frac{f^{v}-f}{\tau_{t}}+\frac{f^{r}-f^{v}}{\tau_{r}}+\frac{g-f^{r}}{\tau_{v}}\equiv\frac{f^{v}-f}{\tau_{t}}+Q_{v},$
(20)

where $f(\bm{x}_{i+1/2,j_{m},k_{n}},t,\bm{u},\xi_{r},\xi_{v})$ is the number
density of molecules at position $(x_{1},x_{2},x_{3})^{T}$ and time $t$, with
particle velocity $\bm{u}=(u_{1},u_{2},u_{3})^{T}$ and internal energy
$(\xi_{r},\xi_{v})^{T}$. The left hand side of Eq.(20) represents the free
streaming of molecules in space, and the right side denotes the collision
term. In Eq.(20), two intermediate equilibrium states $f^{v}$ and $f^{r}$, and
Maxwellian distribution $g$ are introduced with three temperatures $T_{t}$,
$T_{r}$ and $T_{v}$ as

$\displaystyle f^{v}$
$\displaystyle=\rho(\frac{\lambda_{t}}{\pi})^{\frac{N_{t}}{2}}e^{-\lambda_{t}(u_{i}-U_{i})^{2}}(\frac{\lambda_{r}}{\pi})^{\frac{N_{r}}{2}}e^{-\lambda_{r}\xi_{r}^{2}}(\frac{\lambda_{v}}{\pi})^{\frac{N_{v}}{2}}e^{-\lambda_{v}\xi_{v}^{2}},$
(21) $\displaystyle f^{r}$
$\displaystyle=\rho(\frac{\lambda_{t}}{\pi})^{\frac{N_{t}}{2}}e^{-\lambda_{t}(u_{i}-U_{i})^{2}}(\frac{\lambda_{r}}{\pi})^{\frac{N_{r}}{2}}e^{-\lambda_{r}\xi_{r}^{2}},$
(22) $\displaystyle g$
$\displaystyle=\rho(\frac{\lambda_{t}}{\pi})^{\frac{N_{t}}{2}}e^{-\lambda_{t}(u_{i}-U_{i})^{2}},$
(23)

where $\rho$ is the density, $\bm{U}$ denotes three-dimension velocities
$(U_{1},U_{2},U_{3})^{T}$, $\lambda_{t}=m_{0}/(2k_{B}T_{t})$ is related to the
translational temperature $T_{t}$, $\lambda_{r}=m_{0}/(2k_{B}T_{r})$ and
$\lambda_{v}=m_{0}/(2k_{B}T_{v})$ account for the rotational temperature
$T_{r}$ and vibrational temperature $T_{v}$, respectively. Above extended BGK-
type kinetic model has the similarity with the two relaxation time BGK models
for gases with internal d.o.f. [47, 48]. For triatomic molecule $CO_{2}$, we
reasonably adopt above three-temperature extend BGK-type model, with the well-
calibrated relaxation time as Eq.(11), Eq.(18) and Eq.(19). The right-hand-
side collision operator as Eq.(20) contains three terms, corresponding to the
elastic collision ($f\to f^{v}$) and inelastic collisions ($f^{v}\to f^{r}$
and $f^{r}\to g$). As shown in Figure 2, for multiple relaxation precess of
$CO_{2}$, notice that the inelastic collision takes longer time than that of
elastic collision. The additional term $Q_{v}$ resulting from the inelastic
collisions accounts for the internal energy exchange among the translational,
rotational and vibrational internal energy.

The relation between density $\rho$, momentum $\rho\bm{U}$, total energy $\rho
E$, rotational internal energy $\rho E_{r}$, and internal vibrational energy
$\rho E_{v}$ are determined by taking moments of the intermediate equilibrium
distribution function $f^{v}$ as

$\displaystyle\bm{Q}\equiv\int\bm{\psi_{v}}f^{v}\text{d}\Xi_{v}=\begin{pmatrix}\rho,\
\rho\bm{U},\ \rho E,\ \rho E_{r},\ \rho E_{v}\end{pmatrix}^{T},$ (24)

with the vector of extended collision invariants
$\bm{\psi_{v}}=(1,u_{1},u_{2},u_{3},\frac{1}{2}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2}+\xi_{r}^{2}+\xi_{v}^{2}),\frac{1}{2}\xi_{r}^{2},\frac{1}{2}\xi_{v}^{2})^{T}$
and
$\text{d}\Xi_{v}=\text{d}u_{1}\text{d}u_{2}\text{d}u_{3}\text{d}\xi_{r}\text{d}\xi_{v}$.
In Eq.(24), total energy is $E=\bm{U}^{2}/2+E_{t}+E_{r}+E_{v}$. Eq.(20)
introduces the new rotational and vibrational temperatures $T_{r}$ and
$T_{v}$, thus, the constraints of rotational and vibrational internal energy
relaxation have to be imposed on the extended kinetic model to self-
consistently determine all unknowns. Since only mass, momentum and total
energy are conserved during molecule collisions, the compatibility condition
for the collision term turns into

$\displaystyle\bm{S}\equiv\int(\frac{f^{v}-f}{\tau_{t}}+Q_{v})\bm{\psi_{v}}\text{d}\Xi_{v}=(0,0,0,0,0,S_{r},S_{v})^{T}.$
(25)

Source terms $S_{r}$ and $S_{v}$ are from the internal energy exchange among
translational, rotational and vibrational modes during inelastic collision.
These source terms cannot be derived from the BGK model itself. These two
source terms for the rotational internal energy and vibrational internal
energy can be modeled through the Landau-Teller-Jeans-type relaxation model
[34, 36], which read

$\displaystyle S_{r}$ $\displaystyle=\frac{(\rho E_{r})^{eq}-\rho
E_{r}}{\tau_{r}},$ (26) $\displaystyle S_{v}$ $\displaystyle=\frac{(\rho
E_{v})^{eq}-\rho E_{v}}{\tau_{v}}.$ (27)

For $CO_{2}$, the rotational and vibrational relaxation time $\tau_{r}$ and
$\tau_{v}$ have been determined by Eq.(18) and Eq.(19). The left unknown
equilibrium rotational internal energy $(\rho E_{r})^{eq}$ and vibrational one
$(\rho E_{v})^{eq}$ are determined by the assumption
$T_{v}=T_{r}=T_{t}=T^{eq}$ [34, 36], i.e.,

$\displaystyle T_{eq}$
$\displaystyle\equiv\frac{N_{t}T_{t}+N_{r}T_{r}+N_{v}T_{v}}{N_{t}+N_{r}+N_{v}},$
(28) $\displaystyle(\rho E_{r})^{eq}$ $\displaystyle=\frac{N_{r}}{2}\rho
RT_{eq},$ (29) $\displaystyle(\rho E_{v})^{eq}$
$\displaystyle=\frac{N_{v}}{2}\rho RT_{eq},$ (30)

with the translational d.o.f. as $N_{t}$, rotational d.o.f. as $N_{r}$ and
vibrational d.o.f. as $N_{v}$. Up to this point, the thermal non-equilibrium
three-temperature BGK-type model for $CO_{2}$ is completed with the well-
determined relaxation time $\tau_{t}$, $\tau_{r}$, and $\tau_{v}$, as well as
the modeling source terms $S_{r}$ and $S_{v}$.

Using the intermediate equilibrium state $f^{v}$, with the frozen of
rotational and vibrational internal energy exchange, the $1$st-order Chapman-
Enskog expansion [17] gives

$\displaystyle f=f^{v}-\tau_{t}(\frac{\partial{f^{v}}}{\partial
t}+u_{i}\frac{\partial f^{v}}{\partial x_{i}}),$ (31)

from which the corresponding non-equilibrium three-temperature macroscopic
governing equations in three-dimensions can be derived [36]. The technical
details in the standard derivation of three-temperature macroscopic equations
are similar as that of deriving two-temperature macroscopic equations [49],
which are omitted in this paper. The key finding is that two additional
equations get involved with above source terms Eq.(26) and Eq.(27), which
govern the evolutionary dynamics of rotational internal energy $\rho E_{r}$
and vibrational internal energy $\rho E_{v}$. Thus current thermal non-
equilibrium three-temperature macroscopic system goes beyond the one-
temperature supersonic turbulent channel flow of thermally perfect gas [40].
The non-equilibrium three-temperature macroscopic equations give the fixed
Prandtl number $Pr=1$. In the numerical simulation of supersonic $CO_{2}$
turbulent channel flow, the heat flux through the cell interface will be
corrected to obtain the targeted $Pr$ to any realistic value [42]. The three-
temperature BGK-type equation as Eq.(20) is going to be solved by following
non-equilibrium high-accuracy GKS, as the numerical fluxes at cell interfaces
are evaluated based on the time-dependent gas distribution solution.

###  3.2 Non-equilibrium high-accuracy gas-kinetic scheme

For finite-volume non-equilibrium GKS, the key procedure is updating the
macroscopic flow variables inside each control volume through the numerical
fluxes. In this section, the spatial and temporal high-accuracy non-
equilibrium GKS is proposed within the two-stage fourth-order framework [21].

Taking moments of the extended three-temperature BGK-type model as Eq.(20) and
integrating with respect to control volume, 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},$
(32)

where $\bm{Q}_{ijk}$ is the cell averaged macroscopic variables as Eq.(24),
$\bm{S}_{ijk}$ is the cell averaged source term as Eq.(25). 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 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$
(33)
$\displaystyle=\sum_{m,n=1}^{2}\omega_{mn}\int\bm{\psi_{v}}u_{1}f(\bm{x}_{i+1/2,j_{m},k_{n}},t,\bm{u},\xi_{r},\xi_{v})\text{d}\Xi_{v}\Delta
x_{2}\Delta x_{3},$

where $\bm{n}$ is the outer normal direction. The Gaussian quadrature is used
over the cell interface for Eq.(33), 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 numerical fluxes for Eq.(20), the secondary relaxation term
$Q_{v}$ is splitly taken into account as the source term $\bm{S}_{ijk}$ in
Eq.(32).

Without considering $Q_{v}$, the gas distribution function
$f(\bm{x}_{i+1/2,j_{m},k_{n}},t,\bm{u},\xi_{r},\xi_{v})$ in the local
coordinate can be obtained by the integral solution of Eq.(20) as

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

where $\bm{x}^{\prime}=\bm{x}_{i+1/2,j_{m},k_{n}}-\bm{u}(t-t^{\prime})$ is the
trajectory of molecule on grids, $f_{0}$ the initial gas distribution
function, and $f^{v}$ the corresponding intermediate equilibrium state in the
form of Eq.(21). Along the line of GKS [18], for the multi-dimensional kinetic
solver, $f^{v}$ and $f_{0}$ can be constructed as

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

and

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

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

$\displaystyle\int\bm{\psi_{v}}f^{v}_{0}\text{d}\Xi=\int_{u_{1}>0}\bm{\psi_{v}}f^{v}_{l}\text{d}\Xi_{v}+\int_{u_{1}<0}\bm{\psi_{v}}f^{v}_{r}\text{d}\Xi_{v}.$
(37)

Substituting $f^{v}$ and $f_{0}$ into Eq.(34), 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_{r},\xi_{v})$
$\displaystyle=(1-e^{-t/\tau_{t}})f^{v}_{0}+((t+\tau_{t})e^{-t\tau_{t}}-\tau_{t})(\overline{a}_{1}u_{1}+\overline{a}_{2}u_{2}+\overline{a}_{3}u_{3})f^{v}_{0}$
(38) $\displaystyle+(t-\tau_{t}+\tau_{t}e^{-t\tau})\overline{A}f^{v}_{0}$
$\displaystyle+e^{-t/\tau_{t}}f^{v}_{l}[1-(\tau_{t}+t)(a_{1}^{l}u_{1}+a_{2}^{l}u_{2}+a_{3}^{l}u_{3})-\tau_{t}A_{l}]H(u_{1})$
$\displaystyle+e^{-t/\tau_{t}}f^{v}_{r}[1-(\tau_{t}+t)(a_{1}^{r}u_{1}+a_{2}^{r}u_{2}+a_{3}^{r}u_{3})-\tau_{t}A_{r}](1-H(u_{1})).$

With the relation of macroscopic variables and intermediate equilibrium
distribution function $f^{v}$, 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 Eq.(38)
can be determined and details are presented in Appendix C. It is noticed that
Eq.(38) provides a gas evolution process from kinetic scale to hydrodynamic
scale on 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
[18, 35], and has advantages in comparison with traditional hydrodynamic
solver in which the Riemann solver [50] and central difference are used for
the inviscid and viscous terms splitly. For Prandtl number fix, similar to Ref
[18], the total energy flux $\bm{F}{(\rho E)}$ in Eq.(33) is modified as
$\bm{F}^{new}{(\rho E)}=\bm{F}{(\rho E)}+(1/Pr-1)q$, where the time-dependent
heat flux through the cell interface can be evaluated precisely by
$q=\int(u-U)\\{[(u_{i}-U_{i})^{2}+\xi_{r}^{2}+\xi_{v}^{2}]/2\\}fd\Xi_{v}$. As
presented in Section 2.1, $Pr$ depends on both the translational and
vibrational temperatures, which is computed locally in each time step for each
computational grid.

With the time-dependent kinetic flux as Eq.(38), the second-order accuracy can
be achieved by one step integration. To achieve high-order accuracy in space
and time, the fifth-order WENO-Z spatial reconstruction [37] and two-stage
fourth-order time discretization [20, 21] are implemented. For source term in
Eq.(32), the one-step forward Euler method is applied in two-stage updating
process to guarantee the robustness. Notice that the $T_{v}$ and $N_{v}$ as
Eq.(3) and Eq.(8) should be calculated in the center of control volume for
source term, as well as at the cell interface when calculating non-equilibrium
fluxes in Appendix C. Newton–Raphson method is utilized to compute the $T_{v}$
at $n+1$ step. With the initial guess $T_{v}$ at $n$ step, several iterations
are enough to obtain the convergent solution, i.e., with convergence error
$\epsilon_{T_{v}}\equiv||T_{v}^{*}-T_{v}^{**}||/T_{v}^{n}\leq 10^{-8}$ where
$T_{v}^{*}$ and $T_{v}^{**}$ are the successive iterated value during the
iteration process. $||\phi||$ denotes the absolute value of $\phi$. The
present non-equilibrium high-accuracy GKS has been constructed with the
second-order kinetic flux as Eq.(38), fifth-order WENO-Z reconstruction, two-
stage fourth-order time discretization. Current non-equilibrium high-accuracy
GKS is well implemented in the in-house parallel computational platform for
turbulence simulation [26, 28], and the DNS in supersonic $CO_{2}$ turbulent
channel flow is presented subsequently.

##  4 Numerical simulation and discussion

Numerical simulations in supersonic turbulent channel flows with bulk Mach
number $Ma=3.0$ and bulk Reynolds number $Re=4880$ are implemented in this
section. The benchmark as equilibrium one-temperature supersonic turbulent
channel flow is validated with HGKS [27, 28] firstly. The supersonic thermal
non-equilibrium three-temperature $CO_{2}$ turbulent channel flow is simulated
with detailed physical analysis.

###  4.1 One-temperature supersonic turbulent channel flow

The computational studies of supersonic turbulent channel flow [38, 40, 51,
39] have been extensively carried out to study the compressible turbulent
boundary layer. One-temperature supersonic turbulent channel flow with bulk
Mach number $Ma=3.0$ and bulk Reynolds number $Re=4880$ [38, 39] is firstly
used to validate the high-accuracy of in-house HGKS solver with non-uniform
grids [27]. In the computation, the physical domain is $(x,y,z)\in[0,4\pi
H]\times[-H,H]\times[0,4\pi H/3]$ and the computational domain takes
$(\xi,\eta,\zeta)\in[0,4\pi H]\times[0,3\pi H]\times[0,4\pi H/3]$. The
coordinate transformation is given as previous studies [27, 28]. The periodic
boundary conditions are used in streamwise $X$-direction and spanwise
$Z$-directions, and the non-slip and isothermal boundary conditions are
utilized in wall-normal $Y$-direction. In what follows, note that $X$-, $Y$\-
and $Z$\- directions are equivalent as $x_{1}$-,$x_{2}$\- and $x_{3}$\-
directions in Section 3.2. $(U_{1},U_{2},U_{3})^{T}$ is re-expressed in
$(U,V,W)^{T}$ for convenient comparisons with refereed studies.

The turbulent channel flow is initiated with uniform density $\rho=1$, and the
initial streamwise velocity $U(y)$ profile is given by the perturbed
Poiseuille flow profile $U(y)=1.5(1-y^{2})+\text{white noise}$, where the
white noise is added with $10\%$ amplitude of local streamwise velocity. The
spanwise and wall-normal velocity is initiated with white noise. The initial
uniform pressure is computed through the corresponding intial bulk Mach number
and Reynolds number. The targeted non-dimensional parameters bulk Mach number
$Ma$ and bulk Reynolds number $Re$ are defined as $Ma=U_{b}/c_{w}$,
$Re=\rho_{b}U_{b}H/\mu_{w}$. The bulk velocity $U_{b}$ and bulk density
$\rho_{b}$ are given by $U_{b}=\int_{-H}^{H}U(y)\text{d}y$ and
$\rho_{b}=\int_{-H}^{H}\rho(y)\text{d}y$. $H=1$ is the half height of the
channel, $\mu_{w}$ the wall molecule viscosity. $c_{w}=\sqrt{\gamma RT_{w}}$
is the wall sound speed, $T_{w}$ the wall temperature and $R$ the gas
constant. $T_{w}$ is set to $1$ in current simulation. In current on-
temperature validation case, the shear viscosity $\mu$ adopts the exact same
power law $\mu(T)\propto T^{0.7}$ [38], where $T$ is the equilibrium
temperature. The plus unit $Y^{+}$ and plus velocity $U^{+}$ are defined as
$Y^{+}=\rho u_{\tau}y/\mu$, $U^{+}=U/u_{\tau}$ with the friction velocity
$u_{\tau}=\sqrt{\tau_{w}/\rho_{w}}$, the wall shear stress
$\tau_{w}=\mu_{w}\partial U/\partial y\big{|}_{w}$, and the wall density
$\rho_{w}$. The friction Mach number $Ma_{\tau}$ and the friction Reynolds
number $Re_{\tau}$ are given by $Ma_{\tau}=u_{\tau}/{c_{w}}$ and
$Re_{\tau}=H/\delta_{v}$ with $\delta_{v}=\mu_{w}/(\rho_{w}u_{\tau})$. The
heat flux $q_{w}$ and the non-dimensional heat flux $B_{q}$ of the wall are
defined as $q_{w}=-\kappa\partial T/\partial y\big{|}_{w}$,
$B_{q}=q_{w}/(\rho_{w}C_{P}u_{\tau}T_{w})$. In the one-temperature validation,
the fixed Prandtl number $Pr=0.70$ is used as refereed simulation [38], which
is close to the $Pr$ of air as shown in Figure 1. These statistical quantities
are used to quantitatively validate the performance of HGKS and non-
equilibrium high-accuracy GKS subsequently.

Case | Physical domain | $N_{x}\times N_{y}\times N_{z}$ |  $\Delta Y^{+}_{min}$/$Y^{+}_{N10}$ | $\Delta X^{+}$ | $\Delta Z^{+}$
---|---|---|---|---|---
Ref1 | $4\pi H\times 2H\times 4\pi H/3$ | $144\times 90\times 60$ | 0.20/17 | 39 | 24
Ref2 | $4\pi H\times 2H\times 4\pi H/3$ | $400\times 210\times 320$ | 0.65/- | 14.32 | 5.96
$G_{1}$ | $4\pi H\times 2H\times 4\pi H/3$ | $128\times 128\times 128$ | 0.52/12.94 | 43.36 | 14.45
Table 2:  Supersonic one-temperature turbulent channel flow: numerical
parameters of current validation case $G_{1}$ and the reference simulations
[38, 39]. ”-” means that the data can not be find in the refereed paper.

Figure 3:  Supersonic one-temperature turbulent channel flow: the profiles of
normalized r.m.s. of density $\rho_{rms}^{{}^{\prime}}/\langle\rho\rangle$,
r.m.s. of Mach number $Ma_{rms}^{{}^{\prime}}$ and normalized r.m.s. of
temperature $T_{rms}^{{}^{\prime}}/\langle T\rangle$ (upper), and the profiles
of normalized turbulence intensities $U_{rms}^{{}^{\prime}}/\langle
u_{\tau}\rangle$, $V_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$,
$W_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$ and Reynolds stress
$-\langle\rho U^{{}^{\prime}}V^{{}^{\prime}}\rangle/\langle\tau_{w}\rangle$
(lower).

In this one-temperature supersonic turbulent channel flow simulation, the
details of numerical parameters are given in Table 2. The numerical cases of
DNS in refereed paper[38] and [39] are denoted as Ref1 and Ref2, and case
$G_{1}$ is implemented by one-temperature HGKS [28]. The spectral method and
high-order difference scheme is used in Ref1 and Ref2, respectively. $\Delta
Y^{+}_{min}$ is the first grid center space off the wall in the wall-normal
direction, and $Y^{+}_{N10}$ is the plus unit for the first ten points (grid
center) off the wall. $\Delta X^{+}$ and $\Delta Z^{+}$ are the equivalent
plus unit for uniform streamwise and spanwise grids, respectively. As shown in
Table 2, the grid resolution of case $G_{1}$ meet the requirement for DNS
[38]. The constant moment flux is used to determine the external force [28] in
transition and the fully developed turbulence periods. The supersonic
turbulent channel flow takes longer time to transit than that of near
incompressible turbulent channel flow with $Ma=0.1$. During the fully
developed turbulence period, $680$ characteristic periodic time $H/U_{b}$ is
used to obtain the statistically stationary turbulence. In what follows, the
ensemble average of $\phi$ over time and the $X$\- and $Z$-directions is
represented by $\langle\phi\rangle$. The fluctuation of $\phi$ is denoted by
$\phi^{{}^{\prime}}=\phi-\langle\phi\rangle$, and the root-mean square
(r.m.s.) of $\phi$ is defined as
$\phi_{rms}^{{}^{\prime}}=\sqrt{\langle(\phi-\langle\phi\rangle)^{2}\rangle}$,
where $\phi$ represents the density, temperature and velocity, etc. To further
quantify the performance of HGKS in one-temperature supersonic turbulent
channel flow, the profiles of normalized r.m.s. of density
$\rho_{rms}^{{}^{\prime}}/\langle\rho\rangle$, r.m.s. of Mach number
$Ma_{rms}^{{}^{\prime}}$, and normalized r.m.s. temperature
$T_{rms}^{{}^{\prime}}/\langle T\rangle$, and the profiles of normalized
turbulence intensities (r.m.s. of velocities as $U_{rms}^{{}^{\prime}}/\langle
u_{\tau}\rangle$, $V_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$,
$W_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$) and Reynolds stress
$-\langle\rho U^{{}^{\prime}}V^{{}^{\prime}}\rangle/\langle\tau_{w}\rangle$
are presented in Figure 3. In order to account for the mean property of
variations caused by compressibility, the Van Driest (VD) transformation [52]
is proposed for the mean velocity, i.e., density-weighted velocity
${\left\langle U\right\rangle}_{VD}^{+}=\int_{0}^{{\left\langle
U\right\rangle}^{+}}(\left\langle\rho\right\rangle/\left\langle\rho_{w}\right\rangle)^{1/2}\text{d}{\left\langle
U\right\rangle}^{+}$. Compared with the refereed DNS solution [38], the well-
agreed performance of case $G_{1}$ confirms the high-accuracy flow-fields has
been obtained by current HGKS for one-temperature supersonic turbulent channel
flow. The small deviations between the case $G_{1}$ and Ref1 may result from
the numerical solutions in different governing equations. For VD
transformation of streamwise velocity, the solution from case $G_{1}$ matches
well with the DNS in very fine grids [39] with fixed Prandtl number $Pr=0.72$.
Subsequently, non-equilibrium high-accuracy GKS for three-temperature model
proposed in Section 3.2 is used to simulate supersonic thermal non-equilibrium
three-temperature $CO_{2}$ turbulent channel flow.

###  4.2 Three-temperature supersonic $CO_{2}$ turbulent channel flow

For the first time, the DNS in supersonic thermal non-equilibrium three-
temperature $CO_{2}$ turbulent channel flow is implemented. The numerical
setup is same as one-temperature validation case $G_{1}$ in Section 4.1. The
bulk Mach number and bulk Reynolds number take $Ma=3.0$ and $Re=4880$. The
three-temperature initial flow fields restarts from the one-temperature fully
developed turbulence, with initializing three temperatures as
$T_{t}=T_{r}=T_{v}$. With absolute wall temperature $T_{w}=1$, the equivalent
wall temperature $T_{we}=300K$ is adopted to compute the practical physical
temperatures when determining $\tau_{r}$ and $\tau_{v}$. The isothermal
boundary condition is utilized for translational temperature $T_{t}$. To the
author’s knowledge, there is no report on wall boundary condition for
$CO_{2}$. In current study, assuming rotational and vibrational modes of $CO2$
do not exchange internal energy with the wall, thus, adiabatic boundary
condition is used for rotational and vibrational temperatures in wall-normal
$Y$-direction. The realistic wall boundary conditions for rotational and
vibrational internal energy of $CO_{2}$ deserve to be explored by seminal
experimental measurements and theoretical studies.

Figure 4:  Supersonic thermal non-equilibrium three-temperature $CO_{2}$
turbulent channel flow: the time history of mean wall density $[\rho_{w}]$
(left), and the time history of mean total energy $[\rho E]$, mean rotational
internal energy $[\rho E_{r}]$ and mean vibrational internal energy $[\rho
E_{v}]$ (right).

Figure 5:  Supersonic thermal non-equilibrium three-temperature $CO_{2}$
turbulent channel flow: contours of translational, rotational, and vibrational
temperatures.

To balance the wall shear stress, the constant moment flux is used to
determine the external force of supersonic thermal non-equilibrium three-
temperature $CO_{2}$ turbulent channel flow. In Figure 4, note that the mean
average of $\phi$ over whole computational domain (or the two wall planes) is
represented by $[\phi]$. After the long running time, the statistical mean
variables as $[\rho_{w}]$, $[\rho E]$, $[\rho E_{r}]$ and $[\rho E_{v}]$ in
Figure 4 oscillate in a narrow range, indicating that the three-temperature
$CO_{2}$ turbulent channel flow reaches the fully developed state.
Specifically, the mean total energy $[\rho E]$ of one-temperature case $G_{1}$
is approximate $94\%$ of current three-temperature $CO_{2}$. The larger mean
total energy $[\rho E]$ of three-temperature $CO_{2}$ is reasonable, since the
$CO_{2}$ turbulence is equipped with the additional excited vibrational
internal energy. Figure 5 shows three-dimension contours of translational,
rotational, and vibrational temperatures of supersonic $CO_{2}$ turbulent
flows in the fully developed turbulence period. We clearly observe that the
contour of vibrational temperature distinguishes from the translational and
rotational temperature contours, confirming the thermal non-equilibrium
performance.

Case | ${\left\langle\rho_{w}\right\rangle}$ | ${\left\langle u_{\tau}\right\rangle}/U_{b}$ | ${\left\langle Ma_{\tau}\right\rangle}$ | ${\left\langle Re_{\tau}\right\rangle}$ | ${\left\langle B_{q}\right\rangle}$ | ${\left\langle\rho_{c}\right\rangle}$ | ${\left\langle T_{c}^{eq}\right\rangle}$
---|---|---|---|---|---|---|---
Ref1 | 2.388 | 0.0387 | 0.116 | 451 | 0.137 | 0.952 | 2.490
G1 | 2.407 | 0.0376 | 0.113 | 442 | 0.137 | 0.948 | 2.521
C1 | 2.080 | 0.0400 | 0.120 | 406 | 0.127 | 0.966 | 2.084
Table 3:  Supersonic thermal non-equilibrium three-temperature $CO_{2}$
turbulent channel flow: ensemble quantities at the wall and the center plane
of channel.

Figure 6:  Supersonic thermal non-equilibrium three-temperature $CO_{2}$
turbulent channel flow: the profiles of ensemble density $\langle\rho\rangle$,
ensemble Mach number $\langle Ma\rangle$, and ensemble equilibrium temperature
$\langle T^{eq}\rangle$ (left), and the profiles of ensemble streamwise
velocity profiles ${\left\langle U\right\rangle}^{+}$ and VD transformation of
streamwise velocity ${\left\langle U\right\rangle}_{VD}^{+}$ (right).

Figure 7:  Supersonic thermal non-equilibrium three-temperature $CO_{2}$
turbulent channel flow: the profiles of normalized r.m.s. of density
$\rho_{rms}^{{}^{\prime}}/\langle\rho\rangle$, r.m.s. of Mach number
$Ma_{rms}^{{}^{\prime}}$ and normalized r.m.s. of equilibrium temperature
${(T^{eq})}^{{}^{\prime}}_{rms}/\langle T^{eq}\rangle$ (left), and the
profiles of normalized turbulence intensities $U_{rms}^{{}^{\prime}}/\langle
u_{\tau}\rangle$, $V_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$,
$W_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$ and normalized Reynolds stress
$-\langle\rho U^{{}^{\prime}}V^{{}^{\prime}}\rangle/\langle\tau_{w}\rangle$
(right).

To quantitatively analyze the three-temperature performance of supersonic
$CO_{2}$ turbulence, $900$ characteristic periodic time $H/U_{b}$ is used to
obtain the statistically stationary three-temperature turbulence. The key
ensemble quantities at the wall and the center plane of channel are presented
in Table 3. Running case $C_{1}$ denotes the supersonic thermal non-
equilibrium three-temperature $CO_{2}$ case. Again, ensemble quantities from
case $G_{1}$ agree well with these of refereed solution in Ref1 [38]. Compared
with the solutions of one-temperature case $G_{1}$, the ensemble friction
velocity $\langle u_{\tau}\rangle$, the ensemble friction Mach number
${\left\langle Ma_{\tau}\right\rangle}$ and the ensemble center density
${\left\langle\rho_{c}\right\rangle}$ (density at the center plane of channel,
namely at $Y=0$ plane) of case $C_{1}$ are slightly larger. The non-
dimensional heat flux $B_{q}$ of the wall of case $C_{1}$ is slightly smaller
than that of case $G_{1}$. While the ensemble wall density
${\left\langle\rho_{w}\right\rangle}$, the ensemble friction Reynolds number
$Re_{\tau}$, and the ensemble central equilibrium temperature ${\left\langle
T_{c}^{eq}\right\rangle}$ of case $C_{1}$ reduce dramatically. Frictional
force and heat transfer of the wall is of special interest in the long-term
applicability of compressible $CO_{2}$ turbulence. The ensemble frictional
force $\langle\tau_{w}\rangle$ of case $G_{1}$ and case $C_{1}$ is $3.41\times
10^{-3}$ and $3.32\times 10^{-3}$, respectively. While, the ensemble heat flux
$\langle q_{w}\rangle$ of $G_{1}$ and $C_{1}$ is $-0.0309$ and $-0.0396$,
respectively. It is concluded that the thermal non-equilibrium three-
temperature effects of $CO_{2}$ enlarge the ensemble heat transfer by $20\%$,
and slightly decrease the ensemble frictional force.

For the first-order ensemble statistical quantities in supersonic three-
temperature $CO_{2}$ turbulent channel flow, the profiles of key ensemble
quantities are presented in Figure 6. The equilibrium temperature $T^{eq}$ is
defined as Eq.(28). Compared with the one-temperature case $G_{1}$, we observe
the large discrepancies in the near-wall density profiles (approximate
$||Y/H||\geq 0.8$) and the equilibrium temperature profiles in the off-wall
region (approximate $||Y/H||\leq 0.95$). Both the ensemble wall density and
central equilibrium temperature decrease by approximate $15\%$. It can be
concluded that the thermal non-equilibrium three-temperature effects of
supersonic $CO_{2}$ affect its ensemble thermal quantities greatly. In Figure
6, the ensemble streamwise velocity in the log-law region of case $C_{1}$ is
lower than that of case $G_{1}$, indicating that the ${\left\langle
U\right\rangle}^{+}$ in log-law region is suppressed by the non-equilibrium
three-temperature effects. There is little change in the ensemble Mach number
${\left\langle Ma\right\rangle}$ and ensemble VD transformation of streamwise
velocity ${\left\langle U\right\rangle}_{VD}^{+}$, showing that the VD
transformation [52] still works well for supersonic $CO_{2}$ turbulent channel
flow. Figure 7 shows that the peak normalized r.m.s. of density
$\rho_{rms}^{{}^{\prime}}/\langle\rho\rangle$ and equilibrium temperature
${(T^{eq})}^{{}^{\prime}}_{rms}/\langle T^{eq}\rangle$ of supersonic three-
temperature $CO_{2}$ are suppressed near the wall region, and the
corresponding peak locations are much closer to the wall than these of case
$G_{1}$. While, the peak of r.m.s. of Mach number $Ma_{rms}^{{}^{\prime}}$ in
case $C_{1}$ is larger than that of case $G_{1}$. Compared with one-
temperature supersonic turbulent channel flow, Figure 7 shows that the
normalized turbulent intensities of three-temperature $CO_{2}$ are suppressed
above the $Y^{+}\approx 20$, $Y^{+}\approx 100$, and $Y^{+}\approx 40$ regions
for $U_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$,
$V_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$, and
$W_{rms}^{{}^{\prime}}/\langle u_{\tau}\rangle$, respectively.
Correspondingly, the normalized Reynolds stress is suppressed above the
$Y^{+}\approx 120$ region. In supersonic turbulent channel flow, it can be
concluded that the thermal non-equilibrium three-temperature effects of
$CO_{2}$ suppress the peak of normalized r.m.s. of density and temperature,
normalized turbulent intensities and Reynolds stress.

Figure 8:  Supersonic thermal non-equilibrium three-temperature $CO_{2}$
turbulent channel flow: the ensemble translational, rotational, and
vibrational temperature $\langle T_{t}\rangle$, $\langle T_{r}\rangle$,
$\langle T_{v}\rangle$ (left), and the normalized r.m.s. of translational,
rotational, and vibrational temperature $(T_{t})^{{}^{\prime}}_{rms}/\langle
T_{t}\rangle$, $(T_{r})^{{}^{\prime}}_{rms}/\langle T_{r}\rangle$,
$(T_{v})^{{}^{\prime}}_{rms}/\langle T_{v}\rangle$ (right).

Figure 9:  Supersonic thermal non-equilibrium three-temperature $CO_{2}$
turbulent channel flow: slices of translational, rotational, and vibrational
temperatures at $Y^{+}=20.3$ (left column) and $Y^{+}=40.6$ (right column).

Be of special interest in the non-equilibrium performance of three internal
energy modes of supersonic three-temperature $CO_{2}$ turbulent channel flow.
Figure 8 shows the ensemble average of translational, rotational and
vibrational temperatures. At the wall, the three ensemble temperatures sort
from the high to low as $\langle T_{v}\rangle$, $\langle T_{r}\rangle$,
$\langle T_{t}\rangle$, and the opposite order of three temperatures is
observed in the off-wall region (approximate $||Y/H||\leq 0.95$). Figure 8 as
well shows that the peak of normalized r.m.s. of translational temperature
$(T_{t})^{{}^{\prime}}_{rms}/\langle T_{t}\rangle$ and rotational temperature
$(T_{r})^{{}^{\prime}}_{rms}/\langle T_{r}\rangle$ is almost co-located at
$Y^{+}\approx 15$. While the peak location of normalized r.m.s. of vibrational
temperature $(T_{v})^{{}^{\prime}}_{rms}/\langle T_{v}\rangle$ is far from the
wall at $Y^{+}\approx 30$. The peak value of normalized r.m.s. of temperatures
sorting from the high to low is translational temperature, rotational
temperature, and vibrational temperature. Additionally, it is observed that
$(T_{v})^{{}^{\prime}}_{rms}/\langle T_{v}\rangle$ is larger than
$(T_{t})^{{}^{\prime}}_{rms}/\langle T_{t}\rangle$ and
$(T_{r})^{{}^{\prime}}_{rms}/\langle T_{r}\rangle$ above the $Y^{+}\approx 30$
region. Figure 9 clearly shows the streamwise low-temperature and high-
temperature ribbon-like regions for all three temperatures. At $Y^{+}=20.3$,
the low-temperature ribbon-like regions are dominated in vibrational
temperature fields, while the high-temperature ones dominate the translational
temperature fields at $Y^{+}=40.6$. Compared with the vibrational temperature
fields, Figure 9 obviously shows the rotational temperature fields have the
higher similarity with translational temperature fields in temperature
amplitude. As shown in Figure 2, for $CO_{2}$, the rotational relaxation time
$\tau_{r}$ is much closer with the translational relaxation time $\tau_{t}$.
However, the vibrational relaxation time $\tau_{v}$ is thousands of times
larger than that of $\tau_{t}$. Thus, the $T_{t}$ and $T_{r}$ are more likely
to be equilibrium through the interaction between translational and rotational
modes with closer relaxation time. Meanwhile, the much longer relaxation
process of vibrational modes may account for its small correlation with
translational modes and rotational modes.

##  5 Conclusion and remarks

The present paper focuses on thermal non-equilibrium three-temperature effects
of carbon dioxide ($CO_{2}$) in supersonic turbulent channel flow. Essential
ingredient has been addressed for compressible $CO_{2}$ turbulent flows,
namely, the thermal non-equilibrium interactions among translational,
rotational, and vibrational modes. The three vibrational modes of $CO_{2}$ are
addressed, and the double degenerated bending modes $\nu_{2}$ is equipped with
the characteristic vibrational temperature $\theta_{v}=959.66K$, which is much
lower than that of $N_{2}$ and $O_{2}$. Thus, the excitation of vibrational
modes of $CO_{2}$ requires to be modeled and simulated carefully even under
the room temperature. The translational, rotational and vibrational relaxation
time of $CO_{2}$ are calibrated. Then, $CO_{2}$ is modeled in an extended
three-temperature BGK-type model within the well-established kinetic
framework. To achieve high-order accuracy in space and time for simulating
$CO_{2}$ turbulence, the non-equilibrium high-accuracy GKS has been
constructed with the second-order kinetic flux, fifth-order WENO-Z
reconstruction, and two-stage fourth-order time discretization.

With implementing the non-equilibrium high-accuracy GKS in the large-scale
parallel in-house solver, the DNS of supersonic $CO_{2}$ turbulent channel
flow is conducted. Compared with the one-temperature supersonic turbulent
channel flow, the three-temperature effects of $CO_{2}$ are analyzed. The
ensemble frictional force and ensemble heat flux of the wall, as well as the
typical ensemble and fluctuating turbulent quantities of supersonic $CO_{2}$
turbulent channel flow are investigated. Thermal non-equilibrium three-
temperature effects of $CO_{2}$ enlarges the ensemble heat transfer of the
wall by approximate $20\%$, and slightly decreases the ensemble frictional
force. The ensemble density and temperature fields are greatly affected, and
both the ensemble wall density and central equilibrium temperature decrease by
approximate $15\%$. The ensemble streamwise velocity in log-law region is
suppressed. There is little change in VD transformation of streamwise
velocity, which shows that the VD transformation still works well for
supersonic $CO_{2}$ turbulent channel flow. We observe that the peak of
normalized r.m.s. of density and temperature, normalized turbulent intensities
and Reynolds stress are suppressed in supersonic three-temperature $CO_{2}$
turbulent flow.

Numerical simulation confirms the thermal non-equilibrium three-temperature
performance of $CO_{2}$. The streamwise low-temperature and high-temperature
ribbon-like regions are clearly observed for all three temperatures near the
wall region. The vibrational modes of $CO_{2}$ behave quite differently with
rotational and translational modes. The peak positions of normalized r.m.s. of
translational temperature and rotational temperature are almost co-located.
Both the maximum ensemble temperature and normalized r.m.s. temperature sort
from high to low is translational temperature, rotational temperature, and
vibrational temperature. Compared with the vibrational temperature fields, the
rotational temperature fields have the higher similarity with translational
temperature fields both in temperature amplitude and its structure. The much
longer relaxation process of vibrational modes of $CO_{2}$ may account for its
small correlation with translational modes and rotational modes.

In the future, we expect to explore the interdisciplinary studies on $CO_{2}$
transition in supersonic/hypersonic flat plate, i.e., the overheating
phenomenon at the late transitional period [53]. The multi-scale numerical
framework UGKS and UGKWP [29, 54] also provide the solid foundation in further
multi-scale $CO_{2}$ flows, i.e., the Mars re-entry vehicles from rarefied to
continuum regimes. In addition, the realistic wall boundary conditions for
rotational and vibrational internal energy of $CO_{2}$ deserve to be explored
for more delicate simulations in wall-bounded $CO_{2}$ turbulence.

## Ackonwledgement

Many thanks to Dr. Zhao Wang for helpful discussions on the wall boundary
conditions of $CO_{2}$ gas flows. This research is supported by the National
Numerical Windtunnel project, the Department of Science and Technology of
Guangdong Province (2020B1212030001). The authors would like to acknowledge
TaiYi supercomputers in the SUSTech for providing high performance
computational resources.

## Appendix A. Sutherland law and curve-fit shear viscosity of $CO_{2}$

For dilute $CO_{2}$ gas flows, Sutherland law [42] gives the shear viscosity
as

$\mu(T_{t})=\mu_{0}(\frac{T_{t}}{T_{0}})^{\frac{3}{2}}\frac{T_{0}+S}{T_{t}+S},$
(39)

where the reference shear viscosity $\mu_{0}=1.370\times 10^{-5}kg/(m\cdot
s)$, $T_{0}=273K$ and $S=222K$ is approximately valid between $190K$ and
$1700K$. In terms of air, it is noted that $\mu_{0}=1.716\times
10^{-5}kg/m\cdot s$, $T_{0}=273K$ and $S=111K$ for air is approximately valid
between $210K$ and $1900K$. Sutherland law gives the thermal conductivity [42]
as

$\kappa(T_{t})=\kappa_{0}(\frac{T_{t}}{T_{0}})^{\frac{3}{2}}\frac{T_{0}+S}{T_{t}+S},$
(40)

where the reference thermal conductivity $\kappa_{0}=1.46\times
10^{-2}W/(m\cdot K)$, $T_{0}=273K$ and $S=1800K$ is approximately valid
between $180K$ and $700K$. In terms of air, it is noted that
$\kappa_{0}=2.41\times 10^{-2}W/(m\cdot K)$, $T_{0}=273K$ and $S=194$ is
approximately valid between $160K$ and $2000K$.

It is reported that the shear viscosity of $CO_{2}$ [7] is well approximated
over a wide range of temperature up to $20000K$. The curve-fit expression of
shear viscosity reads

$\mu(T_{t})=0.1exp\\{(AlnT_{t}+B)lnT_{t}+C\\},$ (41)

where $A=-0.01952739$, $B=1.047818$, and $C=-14.32212$. In Eq.(41), shear
viscosity $\mu$ is in the unit $kg/(m\cdot s)$.

## Appendix B. Calibrated rotational relaxation time of $CO_{2}$

This appendix provides the power law for rotational relaxation time $\tau_{r}$
of $CO_{2}$. For bulk viscosity $\eta_{b}^{r}$ arising from the rotational
modes, the calibration data of $CO_{2}$ [15] is utilized as Table 4. Based on
the least-square method, with $T_{0}=273K$, calibration data in Table 4
provides the $\tau_{r0}=2.99\times 10^{-10}s$ and $n_{r}=4.47$. It should be
noted that the data of bulk viscosity in high temperature is not adequate.

$T_{t}/K$ | $258.05$ | $274.36$ | $293.24$ | $312.80$ | $330.73$ | $353.15$
---|---|---|---|---|---|---
$\eta_{b}^{r}/\mu$ | $0.226$ | $0.180$ | $0.188$ | $0.198$ | $0.200$ | $0.191$
Table 4:  Calibration data [15] for determining the rotational relaxation time
$\tau_{r}$ of $CO_{2}$.

To absorb the high-temperature information of $\tau_{r}$, the power $n_{r}$
should be corrected based on the refereed calculation of $CO_{2}$ [45]. The
rotational relaxation time is used as $\tau_{r}=Z_{r}\tau_{t}$, where the
rotational collision number $Z_{r}$ is given by Parker equation [44] as

$\displaystyle
Z_{r}=\frac{Z_{r}^{\infty}}{1+(\pi^{3/2}/2)\sqrt{T^{\ast}/T_{t}}+(\pi+\pi^{2}/4)(T^{\ast}/T_{t})}.$
(42)

$T^{\ast}$ is the characteristic temperature of inter-molecular potential, and
$Z_{r}^{\infty}$ is the limiting value. The values $Z_{r}^{\infty}=20.39$ and
$T^{\ast}=91.5K$ are used for $CO_{2}$ [45]. Thence, with the fixed
$T_{0}=273K$ and $\tau_{r0}=2.99\times 10^{-10}s$, the
$Zr\tau_{t}=\tau_{r0}(T/T_{0})^{n_{r}}$ gives the corrected power $n_{r}=1.59$
at $T=2000K$.

## Appendix C. Connection between macroscopic variables and mesoscopic
coefficients

The connection between the spatial derivatives of macroscopic flow variables
and the expansion of intermediate equilibrium distribution function $f^{v}$ as
Eq.(21) reads

$\frac{\partial\bm{Q}}{\partial x_{i}}=\int\frac{\partial f^{v}}{\partial
x_{i}}\bm{\psi_{v}}\text{d}\Xi_{v}\equiv\int
af^{v}\bm{\psi_{v}}\text{d}\Xi_{v},$ (43)

where $a$ denotes the spatial mesoscopic coefficients in Eq.(38) as

$\displaystyle
a=\bm{a}^{T}\bm{\tilde{\psi}_{v}}=a_{1}+a_{2}u_{1}+a_{3}u_{2}+a_{4}u_{3}+a_{5}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2})+a_{6}\xi_{r}^{2}+a_{7}\xi_{v}^{2}.$
(44)

Eq.(43) can be rewritten into following linear system

$\frac{1}{\rho}\frac{\partial\bm{Q}}{\partial
x_{i}}=\Big{(}\frac{1}{\rho}\int\bm{\psi_{v}}\otimes\bm{\tilde{\psi}_{v}}^{T}f^{v}\text{d}\Xi_{v}\Big{)}\bm{a}\equiv\bm{M}\bm{a},$
(45)

Each component of $(a_{1},…,a_{6},a_{7})^{T}$ in Eq.(45) can be determined
uniquely

$\left\\{\begin{aligned} a_{7}&=\frac{4\lambda_{v}^{2}}{N_{v}}B_{6},\\\
a_{6}&=\frac{4\lambda_{r}^{2}}{N_{r}}B_{5},\\\
a_{5}&=\frac{4\lambda_{t}^{2}}{3}(B_{4}-U_{1}B_{1}-U_{2}B_{2}-U_{3}B_{3}-B_{5}-B_{6}),\\\
a_{4}&=2\lambda_{t}B_{3}-2U_{3}a_{5},\\\
a_{3}&=2\lambda_{t}B_{2}-2U_{2}a_{5},\\\
a_{2}&=2\lambda_{t}B_{1}-2U_{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}-(U_{1}^{2}+U_{2}^{2}+U_{3}^{2}+\frac{3}{2\lambda_{t}})a_{5}-\frac{N_{r}}{2\lambda_{r}}a_{6}-\frac{N_{v}}{2\lambda_{v}}a_{7},\end{aligned}\right.$
(46)

with

$\left\\{\begin{aligned} B_{1}&=\frac{1}{\rho}[\frac{\partial(\rho
U_{1})}{\partial x_{i}}-U_{1}\frac{\partial\rho}{\partial x_{i}}],\\\
B_{2}&=\frac{1}{\rho}[\frac{\partial(\rho U_{2})}{\partial
x_{i}}-U_{2}\frac{\partial\rho}{\partial x_{i}}],\\\
B_{3}&=\frac{1}{\rho}[\frac{\partial(\rho U_{3})}{\partial
x_{i}}-U_{3}\frac{\partial\rho}{\partial x_{i}}],\\\
B_{4}&=\frac{1}{\rho}[\frac{\partial(\rho E)}{\partial
x_{i}}-E\frac{\partial\rho}{\partial x_{i}}],\\\
B_{5}&=\frac{1}{\rho}[\frac{\partial(\rho E_{r})}{\partial
x_{i}}-E_{r}\frac{\partial\rho}{\partial x_{i}}],\\\
B_{6}&=\frac{1}{\rho}[\frac{\partial(\rho E_{v})}{\partial
x_{i}}-E_{v}\frac{\partial\rho}{\partial x_{i}}].\end{aligned}\right.$ (47)

For the temporal mesoscopic coefficient in Eq.(38), the relation between
temporal derivatives of macroscopic variables and $f^{v}$ can be written as

$\frac{\partial\bm{Q}}{\partial t}=\int\frac{\partial f^{v}}{\partial
t}\bm{\psi_{v}}\text{d}\Xi_{v}\equiv\int Af^{v}\bm{\psi_{v}}\text{d}\Xi_{v},$
(48)

with

$\displaystyle
A=\bm{A}^{T}\bm{\tilde{\psi}_{v}}=A_{1}+A_{2}u_{1}+A_{3}u_{2}+A_{4}u_{3}+A_{5}(u_{1}^{2}+u_{2}^{2}+u_{3}^{2})+A_{6}\xi_{r}^{2}+A_{7}\xi_{v}^{2}.$
(49)

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

$\displaystyle\int(\frac{\partial{f^{v}}}{\partial t}+u_{i}\frac{\partial
f^{v}}{\partial x_{i}})\bm{\psi_{v}}\text{d}\Xi_{v}=\bm{0}.$ (50)

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

## References

  * [1] Mars helicopter - nasa mars.  https://mars.nasa.gov/technology/helicopter/, 2022.
  * [2] Peter Atkins and Julio de Paula.  Physical chemistry. eight ed, 2006.
  * [3] EV Kustova and EA Nagnibeda.  On a correct description of a multi-temperature dissociating co2 flow.  Chemical physics, 321(3):293–310, 2006.
  * [4] E Kustova, M Mekhonoshina, and A Kosareva.  Relaxation processes in carbon dioxide.  Physics of Fluids, 31(4):046104, 2019.
  * [5] Raymond L Taylor and Steven Bitterman.  Survey of vibrational relaxation data for processes important in the c o 2-n 2 laser system.  Reviews of Modern Physics, 41(1):26, 1969.
  * [6] RJ Stouffer, Syukuro Manabe, and K Bryan.  Interhemispheric asymmetry in climate response to a gradual increase of atmospheric co2.  Nature, 342(6250):660–662, 1989.
  * [7] GRAHAM CANDLER.  Computation of thermo-chemical nonequilibrium martian atmospheric entry flows.  In 5th Joint Thermophysics and Heat Transfer Conference, page 1695, 1990.
  * [8] Y-K Chen, WD Henline, and ME Tauber.  Mars pathfinder trajectory based heating and ablation calculations.  Journal of Spacecraft and Rockets, 32(2):225–230, 1995.
  * [9] Robert D Braun and Robert M Manning.  Mars exploration entry, descent and landing challenges.  In 2006 IEEE Aerospace Conference, pages 18–pp. IEEE, 2006.
  * [10] Elena Kustova and Mariia Mekhonoshina.  Models for bulk viscosity in carbon dioxide.  In AIP Conference Proceedings, volume 2132, page 150006. AIP Publishing LLC, 2019.
  * [11] L Tisza.  Supersonic absorption and stokes’ viscosity relation.  Physical Review, 61(7-8):531, 1942.
  * [12] Xingguo Pan, Mikhail N Shneider, and Richard B Miles.  Power spectrum of coherent rayleigh-brillouin scattering in carbon dioxide.  Physical Review A, 71(4):045801, 2005.
  * [13] Mark S Cramer.  Numerical estimates for the bulk viscosity of ideal gases.  Physics of fluids, 24(6):066102, 2012.
  * [14] Seyed Hossein Jamali, Mariette De Groen, Othonas A Moultos, Remco Hartkamp, Thijs JH Vlugt, Wim Ubachs, and Willem Van De Water.  Rayleigh-brillouin light scattering spectra of co2 from molecular dynamics.  The Journal of Chemical Physics, 151(6):064201, 2019.
  * [15] Yuanqing Wang, Wim Ubachs, and Willem Van De Water.  Bulk viscosity of co2 from rayleigh-brillouin light scattering spectroscopy at 532 nm.  The Journal of Chemical Physics, 150(15):154502, 2019.
  * [16] Prabhu Lal Bhatnagar, Eugene P Gross, and Max Krook.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  Physical review, 94(3):511, 1954.
  * [17] Sydney Chapman and TG Cowling.  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, 1970.
  * [18] Kun Xu.  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, 2001.
  * [19] Kun Xu.  Direct modeling for computational fluid dynamics: construction and application of unified gas-kinetic schemes.  World Scientific, 2015.
  * [20] Jiequan Li and Zhifang Du.  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, 2016.
  * [21] Liang Pan, Kun Xu, Qibing Li, and Jiequan Li.  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, 2016.
  * [22] Song Fu and Qibing Li.  Numerical simulation of compressible mixing layers.  International journal of heat and fluid flow, 27(5):895–901, 2006\.
  * [23] Wei Liao, Yan Peng, and Li-Shi Luo.  Gas-kinetic schemes for direct numerical simulations of compressible homogeneous turbulence.  Physical Review E, 80(4):046702, 2009.
  * [24] G Kumar, Sharath S Girimaji, and J Kerimo.  Weno-enhanced gas-kinetic scheme for direct simulations of compressible transition and turbulence.  Journal of Computational Physics, 234:499–523, 2013.
  * [25] Guiyu Cao, Liang Pan, and Kun Xu.  Three dimensional high-order gas-kinetic scheme for supersonic isotropic turbulence i: criterion for direct numerical simulation.  Computers & Fluids, 192:104273, 2019.
  * [26] Guiyu Cao, Liang Pan, and Kun Xu.  Three dimensional high-order gas-kinetic scheme for supersonic isotropic turbulence ii: Coarse-graining analysis of compressible ksgs budget.  Journal of Computational Physics, 439:110402, 2021.
  * [27] Guiyu Cao, Kun Xu, Liang Pan, and Shiyi Chen.  High-order gas-kinetic scheme in general curvilinear coordinate for iles of compressible wall-bounded turbulent flows.  arXiv preprint arXiv:2107.08609, 2021.
  * [28] Guiyu Cao, Liang Pan, and Kun Xu.  High-order gas-kinetic scheme with parallel computation for direct numerical simulation of turbulent flows.  Journal of Computational Physics, 448:110739, 2022.
  * [29] Kun Xu and Juan-Chen Huang.  A unified gas-kinetic scheme for continuum and rarefied flows.  Journal of Computational Physics, 229(20):7747–7764, 2010.
  * [30] Juan-Chen Huang, Kun Xu, and Pubing Yu.  A unified gas-kinetic scheme for continuum and rarefied flows iii: Microflow simulations.  Communications in Computational Physics, 14(5):1147–1173, 2013\.
  * [31] Yajun Zhu, Chang Liu, Chengwen Zhong, and Kun Xu.  Unified gas-kinetic wave-particle methods. ii. multiscale simulation on unstructured mesh.  Physics of Fluids, 31(6):067105, 2019.
  * [32] Chang Liu, Yajun Zhu, and Kun Xu.  Unified gas-kinetic wave-particle methods i: Continuum and rarefied gas flow.  Journal of Computational Physics, 401:108977, 2020.
  * [33] Xiaojian Yang, Wei Shyy, and Kun Xu.  Unified gas-kinetic wave–particle method for gas–particle two-phase flow from dilute to dense solid particle limit.  Physics of Fluids, 34(2):023312, 2022.
  * [34] Chunpei Cai, Danny D Liu, and Kun Xu.  One-dimensional multiple-temperature gas-kinetic bhatnagar-gross-krook scheme for shock wave computation.  AIAA journal, 46(5):1054–1062, 2008.
  * [35] Guiyu Cao, Hualin Liu, and Kun Xu.  Physical modeling and numerical studies of three-dimensional non-equilibrium multi-temperature flows.  Physics of Fluids, 30(12):126104, 2018.
  * [36] Hualin Liu, Guiyu Cao, and Weifang Chen.  Multiple-temperature gas-kinetic scheme for type iv shock/shock interaction.  Communication in Computational Physics, page 853, 2021.
  * [37] Marcos Castro, Bruno Costa, and Wai Sun Don.  High order weighted essentially non-oscillatory weno-z schemes for hyperbolic conservation laws.  Journal of Computational Physics, 230(5):1766–1792, 2011.
  * [38] Gary N Coleman, John Kim, and Robert D Moser.  A numerical study of turbulent supersonic isothermal-wall channel flow.  Journal of Fluid Mechanics, 305:159–183, 1995.
  * [39] Peng Zhang and Zhenhua Xia.  Contribution of viscous stress work to wall heat flux in compressible turbulent channel flows.  Physical Review E, 102(4):043107, 2020.
  * [40] XiaoPing Chen, XinLiang Li, and ZuChao Zhu.  Effects of dimensional wall temperature on velocity-temperature correlations in supersonic turbulent channel flow of thermally perfect gas.  SCIENCE CHINA Physics, Mechanics & Astronomy, 62(6):1–11, 2019\.
  * [41] John David Anderson.  Hypersonic and high temperature gas dynamics.  Aiaa, 2000.
  * [42] Frank M White and Joseph Majdalani.  Viscous fluid flow, volume 3.  McGraw-Hill New York, 2006.
  * [43] George Emanuel.  Bulk viscosity of a dilute polyatomic gas.  Physics of Fluids A: Fluid Dynamics, 2(12):2252–2254, 1990.
  * [44] JG Parker.  Rotational and vibrational relaxation in diatomic gases.  The Physics of Fluids, 2(4):449–462, 1959.
  * [45] E Kustova, I Alekseev, and L Tan.  Investigation of shock wave structure in co2 based on the continuum and dsmc approaches.  In Journal of Physics: Conference Series, volume 1959, page 012032\. IOP Publishing, 2021.
  * [46] Zhao Wang, Hong Yan, Qibing Li, and Kun Xu.  Unified gas-kinetic scheme for diatomic molecular flow with translational, rotational, and vibrational modes.  Journal of Computational Physics, 350:237–259, 2017.
  * [47] TF Morse.  Kinetic model for gases with internal degrees of freedom.  The physics of fluids, 7(2):159–169, 1964.
  * [48] Henning Struchtrup.  The bgk model for an ideal gas with an internal degree of freedom.  Transport theory and statistical physics, 28(4):369–385, 1999.
  * [49] Guiyu Cao.  High-order gas-kinetic schemes for turbulence modeling and simulation.  Hong Kong University of Science and Technology (Hong Kong), 2021.
  * [50] Eleuterio F Toro.  Riemann solvers and numerical methods for fluid dynamics: a practical introduction.  Springer Science & Business Media, 2013.
  * [51] Ming Yu, Chun-Xiao Xu, and Sergio Pirozzoli.  Genuine compressibility effects in wall-bounded turbulence.  Physical Review Fluids, 4(12):123402, 2019.
  * [52] PG Huang and Gary N Coleman.  Van driest transformation and compressible wall-bounded flows.  AIAA journal, 32(10):2110–2113, 1994.
  * [53] Kenneth J Franko and Sanjiva K Lele.  Breakdown mechanisms and heat transfer overshoot in hypersonic zero pressure gradient boundary layers.  Journal of Fluid Mechanics, 730:491–532, 2013.
  * [54] Yipei Chen, Yajun Zhu, and Kun Xu.  A three-dimensional unified gas-kinetic wave-particle solver for flow computation in all regimes.  Physics of Fluids, 32(9):096108, 2020.