2.1. Distribution functions and macroscopic quantities

Since the characteristic temperatures of the rotational modes of typical diatomic molecules are less than 100 K (Bird Reference Bird1994), the translational and rotational DoF can be regarded as fully excited when the gas temperature is higher than 100 K. Therefore, it is a common choice to use constant values of DoF for these modes. Since the gap between two subsequent discrete levels is much lower for rotational energy than for vibrational energy, the rotational energy of gas molecules takes continuous values approximately. However, the vibrational DoF is temperature dependent, due to large energy gaps between discrete vibrational energy levels (Bird Reference Bird1994; Anderson Reference Anderson2019). Suppose there are $N$ energy levels allowed in the vibrational mode, then $N$ distribution functions $f_i(t,\boldsymbol {x},\boldsymbol {v},I_r)$ are needed to identify the states of molecular gas, where $t$ is the time, $\boldsymbol {x}$ is the spatial coordinates, $\boldsymbol {v}$ is the translational molecular velocity, $I_r$ is the rotational energy and $i=0,1,\ldots,N-1$ indicates the vibrational level with sensible energy $\varepsilon _i=ih\nu$ (the energy measured above zero-point energy, $h$ is the Planck number and $\nu$ is the vibrational frequency). Macroscopic variables, such as the molecular number density $n$, flow velocity $\boldsymbol {u}$, pressure tensor $\boldsymbol {P}$, temperatures $T_t,T_r,T_v$ and heat fluxes $\boldsymbol {q}_t,\boldsymbol {q}_r,\boldsymbol {q}_v$, are obtained by taking the moments of the distribution function and summing over all the $N$ vibrational states,

(2.1)\begin{equation} \left.\begin{array}{c} \displaystyle (n, n\boldsymbol{u},\boldsymbol{P})= \sum_{i}^{N}\int_{0}^{\infty}\int_{-\infty}^{\infty} (1,\boldsymbol{v},m\boldsymbol{c}\boldsymbol{c}) {f_i}\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}I_r, \\ \displaystyle \left(\dfrac{3}{2}k_BT_t,\dfrac{d_r}{2}k_BT_r, \dfrac{{d_v(T_v)}}{2}k_BT_v\right)= \dfrac{1}{n}\sum_{i}^{N}\int_{0}^{\infty} \int_{-\infty}^{\infty}{\left(\dfrac{1}{2}mc^2,I_r, \varepsilon_i\right)f_i}\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}I_r, \\ \displaystyle (\boldsymbol{q}_t,\boldsymbol{q}_r,\boldsymbol{q}_v)= \sum_{i}^{N}\int_{0}^{\infty}\int_{-\infty}^{\infty} \boldsymbol{c}\left(\dfrac{1}{2}mc^2,I_r,\varepsilon_i\right) {f_i}\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}I_r, \end{array}\right\} \end{equation}

where the parentheses are used to collect the set of variables defined in an analogical way and, thus, make the equations more compact; the subscripts $t, r, v$ denote the translational, rotational and vibrational components, respectively; $\boldsymbol {c}=\boldsymbol {v}-\boldsymbol {u}$ is the peculiar (thermal) velocity, $m$ is the molecular mass and $k_B$ is the Boltzmann constant; $d_r$ and $d_v(T_v)$ are the rotational and vibrational number of DoF, respectively. Note that the quantum mechanics must be applied for the vibration energy; thus, for a simple harmonic oscillator, the average energy $E_v$ in the case of the equilibrium Boltzmann distribution over the vibrational energy levels can be calculated as (Bird Reference Bird1994; Nagnibeda & Kustova Reference Nagnibeda and Kustova2009; Anderson Reference Anderson2019)

(2.2)\begin{equation} E_v = \frac{h\nu}{\exp{[h\nu(k_BT_v)^{{-}1}]}-1}. \end{equation}

The vibrational DoF $d_v$ is defined in the sense of an effective value such that $d_vk_BT_v/2=E_v$. Then, $d_v$ is obtained as a temperature-dependent value,

(2.3)\begin{equation} d_v(T_v)=\frac{2T_{{ref}}/T_v}{\exp({T_{{ref}}/T_v})-1}, \end{equation}

where $T_{{ref}}$ is the characteristic temperature of the active vibrational mode defined in terms of the oscillator frequency $T_{{ref}}=h\nu /k_B$. It shows that $d_v$ approaches $0$ when $T_v\ll T_{{ref}}$, while $d_v$ is around $1.54$ when $T_v=2T_{{ref}}$.

We also define the temperature $T_{tr}$ to be the equilibrium temperature between the translational and rotational modes, $T_{tv}$ to be the equilibrium temperature between the translational and vibrational modes, and $T$ to be the equilibrium temperature over all DoF,

(2.4a–c)\begin{equation} T_{tr}=\frac{3T_t+d_rT_r}{3+d_r}, \quad T_{tv}=\frac{3T_t+{{d_v(T_v)}}T_v}{3+{{d_v(T_{tv})}}}, \quad T=\frac{3T_t+d_rT_r+{{d_v(T_v)}}T_v}{3+d_r+{{d_v(T)}}}, \end{equation}

and the corresponding pressures are $(p_t, p_r, p_v, p, p_{tr}, p_{tv}) = nk_B(T_t, T_r, T_v, T, T_{tr}, T_{tv})$.

2.2. Generalized kinetic model with vibrational radiation transition

The distribution function $f_i(t,\boldsymbol {x},\boldsymbol {v},I_r)$ is changed due to the gas–gas interactions $J_{{gas}}$, the gas–photon interactions $J_{{photon}}$ that exchanges energy between vibrational mode and radiation field, and the streaming $\mathcal {D}$. Ignoring the momentum exchange between gas and photon (which is an acceptable approximation in most of the practical situations, since the momentum exchange rate is proportional to the ratio of radiative heat flux to the speed of light), the evolution of the distribution function under external body acceleration $\boldsymbol {a}$ is governed by

(2.5)\begin{equation} \underbrace{\frac{\partial{f_{i}}}{\partial{t}}+\boldsymbol{v} \boldsymbol{\cdot}\frac{\partial{f_{i}}}{\partial{\boldsymbol{x}}}+ \frac{\partial{(\boldsymbol{a}f_{i})}}{\partial {\boldsymbol{v}}}}_{\mathcal{D}f_i}=J_{{gas},i} + J_{{photon},i}. \end{equation}

The kinetic model of binary gas–gas collisions has been well established by the WCU equation. When the rotational mode is treated by classical mechanics, it reads

(2.6)\begin{equation} J_{{gas},i} = \sum_{i'j'}\sum_{j}\int_{-\infty}^{\infty} \int_{4{\rm \pi}}{\left(\frac{g_{i}g_{j}}{g_{i'}g_{j'}}f_{i'}f_{j'} -f_{i}f_{j}\right)|\boldsymbol{v}-\boldsymbol{v}_*| \sigma_{ij}^{i'j'}\,\mathrm{d}\varOmega\,\mathrm{d}\boldsymbol{v}_*}, \end{equation}

where $\boldsymbol {v}$ and $\boldsymbol {v}_*$ are the pre-collision velocities of the two molecules with vibrational states $i$ and $j$, respectively, the superscript $'$ indicates the state after collision; $g_i$ is the degeneracy of the vibrational states $i$, $\sigma _{ij}^{i'j'}$ is the scattering cross-section and $\varOmega$ is the solid angle. Since the total energy is conserved during the collision, such a transition occurs only when

(2.7)\begin{equation} |\boldsymbol{v}'-\boldsymbol{v}'_*|^2=|\boldsymbol{v}- \boldsymbol{v}_*|^2+\frac{4}{m}(e_{i}+e_{j}-e_{i'}-e_{j'})>0, \end{equation}

where e is the total internal energy. When the collision is admissible, the molecular velocity after the collision is

(2.8a,b)\begin{equation} \boldsymbol{v}'=\frac{\boldsymbol{v}+\boldsymbol{v}_\ast}{2} +\frac{|\boldsymbol{v}'-\boldsymbol{v}'_*|}{2}{\varOmega}, \quad \boldsymbol{v}'_\ast{=}\frac{\boldsymbol{v}+\boldsymbol{v}_\ast}{2} -\frac{|\boldsymbol{v}'-\boldsymbol{v}'_*|}{2}{\varOmega}. \end{equation}

Since each of the energy states $i,j,i',j'$ has $N$ possible choices, the computational cost will be $N^4$ times higher than that of the BCO, posing urgent needs to develop simplified gas kinetic models. Also, it should be noted that the collision is called elastic when $|\boldsymbol {v}'-\boldsymbol {v}'_*|=|\boldsymbol {v}-\boldsymbol {v}_*|$ since the kinetic energy is conserved. Otherwise, it is inelastic.

Molecular bound-bound radiation occurs as a result of radiative transitions between the quantized energy levels of gas molecules. These radiative events are determined by the molecular structure itself, but are not associated directly with the intermolecular collisions. For gas–photon interactions, the following three processes of vibrational energy transitions that are related to the radiation are considered in the present work:

(2.9)\begin{equation} \left. \begin{array}{c@{}} \text{spontaneous emission:}\ M_j \rightarrow M_i+h\nu_{ij}, \\ \text{stimulated emission:}\ M_j+h\nu_{ij} \rightarrow M_i+2h\nu_{ij}, \\ \text{absorption:}\ M_i+h\nu_{ij} \rightarrow M_j. \end{array} \right\} \end{equation}

Here $\nu _{ij}=(\varepsilon _j-\varepsilon _i)/h$ is the frequency of the photon absorbed or emitted during these transitions; $\varepsilon _j>\varepsilon _i$ are the corresponding energy of the vibrational energy levels $i$ and $j$ of the molecule $M$. Since the photon frequencies have $N(N-1)/2$ discretized values, the radiation field is described by the same number of radiation intensity functions $I_{\nu _{ij}}^{R}(t,\boldsymbol {x},\boldsymbol {\varOmega })$, which measure the energy fluxes per unit solid angle of photons propagating along the direction $\boldsymbol {\varOmega }$ with the frequency $\nu _{ij}$.

The transition changes the molecular population of the energy levels $i$ and $j$. The rates of change of distribution function $f_i$ during processes (2.9) can be calculated by the Einstein coefficients, $A_{ji},B_{ji},B_{ij}$ for the spontaneous emission, stimulated emission and absorption, which are denoted by the superscripts ${sp},{st}$ and ${ab}$, respectively,

(2.10)\begin{equation} \left. \begin{array}{c@{}} \displaystyle \left(\dfrac{\mathrm{d}{f_i}}{\mathrm{d}{t}}\right)^{{sp}}_j = \int_{4{\rm \pi}}{A_{ji}\,\mathrm{d}\varOmega}f_j, \\ \displaystyle \left(\dfrac{\mathrm{d}{f_i}}{\mathrm{d}{t}}\right)^{{st}}_j = \int_{4{\rm \pi}}{B_{ji}I_{\nu_{ij}}^{R}\,\mathrm{d}\varOmega}f_j, \\ \displaystyle \left(\dfrac{\mathrm{d}{f_i}}{\mathrm{d}{t}}\right)^{{ab}}_j ={-}\int_{4{\rm \pi}}{B_{ij}I_{\nu_{ij}}^{R}\,\mathrm{d}\varOmega}f_i. \end{array}\right\} \end{equation}

In the relativistic gas flow, the Doppler and aberration effects describe the shift of photon frequency and the change of propagation direction, respectively, and have a magnitude proportional to the ratio of flow velocity to the speed of light. When the non-relativistic gas flow is considered, these effects can be ignored. Therefore, based on (2.10) the rates of change of the distribution function $f_i$ due to all possible vibrational energy transitions, which give the gas–photon interaction term $J_{{photon},i}$, yield (Groppi & Spiga Reference Groppi and Spiga1999)

(2.11)\begin{align} J_{{photon},i} &=\sum_{j>i}\int_{4{\rm \pi}}{[A_{ji}f_j+ (B_{ji}f_j-B_{ij}f_i)I_{\nu_{ij}}^{R}]\,\mathrm{d}\varOmega}\nonumber\\ & \quad -\sum_{j< i}\int_{4{\rm \pi}}{[A_{ij}f_i+ (B_{ij}f_i-B_{ji}f_j)I_{\nu_{ji}}^{R}]\,\mathrm{d}\varOmega}. \end{align}

Next, we consider the evolution of radiation intensity due to the interaction with gas molecules $J_{\nu _{ij}}^{R}$. Ignoring the photon scattering (which is an accepted assumption when no particles or droplets are considered in the gas flow, and the radiation occupies a spectral range such that wavelengths are much larger compared with molecular size), the evolution of intensity $I_{\nu _{ij}}^{R}$ is

(2.12)\begin{equation} \frac{1}{c_l}\frac{\partial{I_{\nu_{ij}}^{R}}}{\partial{t}} +\boldsymbol{n} \boldsymbol{\cdot} \frac{\partial{I_{\nu_{ij}}^{R}}}{\partial {\boldsymbol{x}}}=J_{\nu_{ij}}^{R}, \end{equation}

where $c_l$ is the speed of light and $\boldsymbol {n}$ is the unit vector of photon propagation direction. When the radiation passes through the matter over distance $c_l\,\mathrm {d}{t}$, it is attenuated by a constant fraction $k_{\nu _{ij}}$, and it has a gain part $j_{\nu _{ij}}$ that does not depend on the photon intensity (Casto Reference Casto2004). Therefore, the rate of change of intensity $I_{\nu _{ij}}^{R}$ due to the gas–photon interaction can be written as

(2.13)\begin{equation} \left. \begin{array}{c@{}} \displaystyle \dfrac{1}{c_l}\left(\dfrac{\mathrm{d}{I_{\nu_{ij}}^{R}}}{\mathrm{d}{t}} \right)^{{loss}} ={-}k_{\nu_{ij}}I_{\nu_{ij}}^{R}, \\ \displaystyle \dfrac{1}{c_l}\left(\dfrac{\mathrm{d}{I_{\nu_{ij}}^{R}}}{\mathrm{d}{t}} \right)^{{gain}} = j_{\nu_{ij}}. \end{array} \right\} \end{equation}

According to processes (2.9), the gain part comes from the spontaneous emission alone, and the loss part is the difference between the absorption and stimulated emission, which indicates the decay of radiation intensity caused by the gas–photon interactions during propagation. Based on (2.10), the frequency-dependent absorptivity $k_{\nu _{ij}}$ and emissivity $j_{\nu _{ij}}$ are determined by the Einstein coefficients as

(2.14a,b)\begin{equation} k_{\nu_{ij}}=h\nu_{ij}\left(B_{ij}\int_{-\infty}^{\infty}{f_i\, \mathrm{d}\boldsymbol{v}}-B_{ji}\int_{-\infty}^{\infty}{f_j\, \mathrm{d}\boldsymbol{v}}\right), \quad j_{\nu_{ij}}= h\nu_{ij}A_{ji}\int_{-\infty}^{\infty}{f_j\,\mathrm{d} \boldsymbol{v}}. \end{equation}

Therefore, the rate of change of the photon intensity $I_{\nu _{ij}}^{R}$ due to gas–photon interactions is

(2.15)\begin{equation} J_{\nu_{ij}}^{R} = j_{\nu_{ij}}-k_{\nu_{ij}}I_{\nu_{ij}}^{R}. \end{equation}

According to Kirchhoff's law, the emissivity is related to the absorptivity as $j_{\nu _{ij}}=B^R_{\nu _{ij}}(T_v)k_{\nu _{ij}}$, where the Planck function $B^R_{\nu _{ij}}(T)$ is defined as

(2.16)\begin{equation} B^R_{\nu_{ij}}(T)=\frac{2h\nu_{ij}^3}{c_l^2} \frac{1}{\exp(h\nu_{ij}/k_BT)-1}. \end{equation}

Then, the evolution of intensity due to the interaction with gas molecules is

(2.17)\begin{equation} \frac{1}{c_l}\frac{\partial{I_{\nu_{ij}}^{R}}}{\partial{t}} +\boldsymbol{n} \boldsymbol{\cdot} \frac{\partial{I_{\nu_{ij}}^{R}}}{\partial {\boldsymbol{x}}}=k_{\nu_{ij}}(B^R_{\nu_{ij}}(T_v) -I_{\nu_{ij}}^{R}). \end{equation}

The macroscopic radiative temperature $T_R$ and heat flux $\boldsymbol {q}_R$ are defined as

(2.18)\begin{equation} (4\sigma_RT_R^4,\boldsymbol{q}_R)=\sum_{\nu_{ij}} \int_{4{\rm \pi}}{I_{\nu_{ij}}^{R}(1,\boldsymbol{n})}\,\mathrm{d}\varOmega, \end{equation}

where $\sigma _R=2{\rm \pi} ^5k_B^4/(15h^3c_l^2)$ is the Stefan–Boltzmann constant. Since the energy exchange between gas and photon is due to the radiation transition of vibrational energy, the radiative temperature $T_R$ approaches vibrational temperature $T_v$ at equilibrium state; hence, the equilibrium temperature in Planck function in the photon kinetic equation (2.17) is $T_v$.

Therefore, the governing equations of the rarefied molecular gas flow with radiation consist of the $N$ equations given by (2.5) of gas distribution functions with collision terms (2.6) and (2.11), and $N(N-1)/2$ equations given by (2.17) of radiative intensities. Then, simplifications of collision terms will be introduced to develop tractable model equations.

2.3. Kinetic model with RTA

2.3.1. Modified Rykov model for gas–gas collisions

Here we build a kinetic model to simplify the WCU collision operator (2.6) for gas–gas interactions based on the Rykov model, not only due to its greater freedom to reflect the relaxation process of heat fluxes, but also due to its much reduced computational complexity. In this model the elastic and inelastic collisions are considered separately with different relaxation times, which can be adjusted to give a correct bulk viscosity. Furthermore, the reference distribution functions to which the distribution function relaxes contain the heat fluxes, so that the thermal conductivity can be recovered. Although the Rykov model was initially proposed for diatomic gas without vibrational modes, it has been extended to polyatomic gas (Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015b) and gases with vibrational modes (Titarev & Frolova Reference Titarev and Frolova2018).

However, to recover the correct thermal relaxation rates other than thermal conductivities, the heat fluxes in the reference distribution functions in the original Rykov model have to be adjusted. Thus, the modified Rykov model for gas–gas collisions with discretized vibrational states becomes

(2.19)\begin{equation} J_{{gas},i} = \frac{g_{t,i}-f_i}{\tau} + \frac{g_{r,i}-g_{t,i}}{Z_r\tau} + \frac{g_{v,i}-g_{t,i}}{Z_{v}\tau} , \end{equation}

where $Z_r$ and $Z_v$ are the rotational and vibrational collision numbers, respectively. The reference distribution functions $g_{t,i}, g_{r,i}, g_{v,i}$ are expanded about the equilibrium distributions $E_t(T) E_r(T) E_{v,i}(T)$ in a series of orthogonal polynomials of peculiar velocity $\boldsymbol {c}$, rotational energy $I_r$, vibrational energy $\varepsilon _i$ and corresponding moments $\boldsymbol {q_t}, \boldsymbol {q_r}, \boldsymbol {q_v}$,

(2.20) \begin{gather} \left. \begin{array}{c@{}} \begin{aligned} g_{t,i} & = nE_t(T_t) E_r(T_r) E_{v,i}(T_v) \left[{1 + \dfrac{2m\boldsymbol{q}_t\boldsymbol{\cdot}{\boldsymbol{c}}}{15 {k_B}{T_t}{p_t}} \left(\dfrac{mc^2}{2k_BT_t}-\dfrac{5}{2}\right)} \right. \\ & \quad \left. +\,\dfrac{2m\boldsymbol{q}_r\boldsymbol{\cdot}{\boldsymbol{c}}}{d_rk_B T_tp_r}\left(\dfrac{I_r}{k_BT_r}-\dfrac{d_r}{2}\right) + \dfrac{2m\boldsymbol{q}_v\boldsymbol{\cdot}{\boldsymbol{c}}}{{{d_v(T_v)}}k_B T_tp_v}\left(\dfrac{\varepsilon_i}{k_BT_v}-\dfrac{{d_v (T_v)}}{2}\right)\right], \end{aligned}\\ \begin{aligned} g_{r,i} & = nE_t(T_{tr}) E_r(T_{tr}) E_{v,i}(T_v) \left[{1 + \dfrac{2m\boldsymbol{q}_0\boldsymbol{\cdot}{\boldsymbol{c}}}{15 {k_B}{T_{tr}}{p_{tr}}} \left(\dfrac{mc^2}{2k_BT_{tr}} -\dfrac{5}{2}\right)} \right. \\ & \quad \left. +\,\dfrac{2m\boldsymbol{q}_1\boldsymbol{\cdot}{\boldsymbol{c}}}{d_r k_BT_{tr}p_{tr}}\left(\dfrac{I_r}{k_BT_{tr}}- \dfrac{d_r}{2}\right) + \dfrac{2m\boldsymbol{q}_2 \boldsymbol{\cdot}{\boldsymbol{c}}}{{{d_v(T_v)}}k_BT_{tr}p_v} \left(\dfrac{\varepsilon_i}{k_BT_v}-\dfrac{{d_v(T_v)}}{2} \right)\right], \end{aligned}\\ \begin{aligned} g_{v,i} & =nE_t(T_{tv}) E_r(T_{r}) E_{v,i}(T_{tv}) \left[{1 + \dfrac{2m\boldsymbol{q}_0\boldsymbol{\cdot}{\boldsymbol{c}}}{15 {k_B}{T_{tv}}{p_{tv}}} \left(\dfrac{mc^2}{2k_BT_{tv}} -\dfrac{5}{2}\right)} \right. \\ & \quad \left. +\,\dfrac{2m\boldsymbol{q}_1\boldsymbol{\cdot}{\boldsymbol{c}}}{d_r k_BT_{tv}p_{r}}\left(\dfrac{I_r}{k_BT_{r}}-\dfrac{d_r}{2}\right) +\dfrac{2m\boldsymbol{q}_2\boldsymbol{\cdot}{\boldsymbol{c}}}{{{d_v(T_{tv})}}k_B T_{tv}p_{tv}}\left(\dfrac{\varepsilon_i}{k_BT_{tv}}-\dfrac{{d_v (T_{tv})}}{2}\right)\right], \end{aligned} \end{array}\right\} \end{gather}

with the equilibrium distribution functions,

(2.21) \begin{equation} \left. \begin{array}{c} \displaystyle E_t(T)={\left(\dfrac{m}{2{\rm \pi} k_BT}\right)}^{3/2} \exp{\left(-\dfrac{mc^2}{2k_BT}\right)}, \\ \displaystyle E_r(T)=\dfrac{I^{d_r/2-1}_{r}}{\varGamma(d_r/2)(k_BT)^{d_r/2}} \exp{\left(-\dfrac{I_r}{k_BT}\right)}, \\ \displaystyle E_{v,i}(T)=\dfrac{g_i}{\displaystyle \sum_j{g_j\exp {\left(-\dfrac{\varepsilon_j}{k_BT}\right)}}}\exp {\left(-\dfrac{\varepsilon_i}{k_BT}\right)}, \end{array}\right\} \end{equation}

where $\varGamma$ is the gamma function, ${\boldsymbol {q}_{0}}$, ${\boldsymbol {q}_{1}}$ and ${\boldsymbol {q}_{2}}$ in (2.20) are linear combinations of translational, rotational and vibrational heat fluxes, which will be designed to recover the exact thermal relaxation rates.

2.3.2. Reduced distribution functions and radiative intensity

In practical numerical simulations it is almost impossible to solve $N(N+1)/2$ equations even with RTA, especially when the number of vibrational energy levels $N$ considered is large. However, the complexity arising from the discrete vibrational energy can be eliminated with the reduced distribution technique (Chu Reference Chu1965; Mathiaud & Mieussens Reference Mathiaud and Mieussens2020). Therefore, we eliminate the rotational energy variable ${I_r}$ and vibrational energy level index $i$ by introducing the following reduced velocity distribution functions ${f_0, f_1, f_2}$:

(2.22)\begin{equation} (f_0,f_1,f_2)=\sum_{i}^{N}\int_{0}^{\infty} (1,I_r,\varepsilon_i)f_i (t,\boldsymbol{x},\boldsymbol{v},I_r) \,\mathrm{d}{I_r}. \end{equation}

Similarly, the total radiative intensity $I^{R}$ will be solved instead of the frequency-dependent ones, and the Planck function $B^{R}(T)$ summed over all frequencies are used,

(2.23) \begin{equation} \left. \begin{array}{c@{}} \displaystyle I^{R}=\sum_{\nu_{ij}}I^{R}_{\nu_{ij}}, \\ \displaystyle B^R(T)=\sum_{\nu_{ij}}B^R_{\nu_{ij}}(T)=\dfrac{1}{\rm \pi}\sigma_RT^4. \end{array} \right\} \end{equation}

Given the fact that the radiation transitions are not correlated with the translational motion of the gas molecules, the governing equations (2.5) and (2.17) with modified Rykov model (2.19) can be transferred to four coupled equations,

(2.24) \begin{equation} \left. \begin{array}{c@{}} \displaystyle \mathcal{D}f_0 = \dfrac{g_{0t}-f_0}{\tau} + \dfrac{g_{0r}-g_{0t}}{Z_r\tau} + \dfrac{g_{0v}-g_{0t}}{Z_v\tau}, \\ \displaystyle \mathcal{D}f_1 = \dfrac{g_{1t}-f_1}{\tau} + \dfrac{g_{1r}-g_{1t}}{Z_r\tau} + \dfrac{g_{1v}-g_{1t}}{Z_v\tau}, \\ \displaystyle \mathcal{D}f_2 = \dfrac{g_{2t}-f_2}{\tau} + \dfrac{g_{2r}-g_{2t}}{Z_r\tau} + \dfrac{g_{2v}-g_{2t}}{Z_v\tau} -\dfrac{f_0}{n}\int_{4{\rm \pi}}{(k^{e}B^R(T_v)-k^{ne}I^{R})\, \mathrm{d}\varOmega}, \\ \displaystyle \dfrac{1}{c_l}\dfrac{\partial{I^{R}}}{\partial{t}}+\boldsymbol{n} \boldsymbol{\cdot} \dfrac{\partial{I^{R}}}{\partial{\boldsymbol{x}}} = k^{e}B^R(T_v)-k^{ne}I^{R}, \end{array}\right\} \end{equation}

with the reduced reference velocity distribution functions

(2.25) \begin{equation} \left. \begin{array}{c@{}} \displaystyle g_{0t}=nE_t(T_t)\left[1+\dfrac{2m\boldsymbol{q}_t\boldsymbol{\cdot} {\boldsymbol{c}}}{15{k_B}{T_t}{p_t}} \left(\dfrac{mc^2}{2k_BT_t}- \dfrac{5}{2}\right)\right], \\ \displaystyle g_{0r}=nE_t(T_{tr})\left[1+\dfrac{2m\boldsymbol{q}_0 \boldsymbol{\cdot}{\boldsymbol{c}}}{15{k_B}{T_{tr}}{p_{tr}}} \left(\dfrac{mc^2}{2k_BT_{tr}}-\dfrac{5}{2}\right)\right], \\ \displaystyle g_{0v}=nE_t(T_{tv})\left[1+\dfrac{2m\boldsymbol{q}_0 \boldsymbol{\cdot}{\boldsymbol{c}}}{15{k_B}{T_{tv}}{p_{tv}}} \left(\dfrac{mc^2}{2k_BT_{tv}}-\dfrac{5}{2}\right)\right], \\ \displaystyle g_{1t}=\dfrac{d_r}{2}k_BT_rg_{0t}+\dfrac{m\boldsymbol{q}_r \boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_t}E_t(T_t), \\ \displaystyle g_{1r}=\dfrac{d_r}{2}k_BT_{tr}g_{0r}+\dfrac{m\boldsymbol{q}_1 \boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_{tr}}E_t(T_{tr}), \\ g_{1v}=\dfrac{d_r}{2}k_BT_{r}g_{0v}+\dfrac{m\boldsymbol{q}_1 \boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_{tv}}E_t(T_{tv}), \\ \displaystyle g_{2t}=\dfrac{{d_v(T_v)}}{2}k_BT_vg_{0t}+\dfrac{m\boldsymbol{q}_v \boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_t}E_t(T_t), \\ \displaystyle g_{2r}=\dfrac{{d_v(T_v)}}{2}k_BT_{v}g_{0r}+\dfrac{m\boldsymbol{q}_2 \boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_{tr}}E_t(T_{tr}), \\ \displaystyle g_{2v}=\dfrac{{d_v(T_{tv})}}{2}k_BT_{tv}g_{0v}+\dfrac{m \boldsymbol{q}_2\boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_{tv}}E_t(T_{tv}), \end{array}\right\} \end{equation}

and the effective absorptivities

(2.26) \begin{equation} \left. \begin{array}{c@{}} \displaystyle k^e=\dfrac{\displaystyle \sum_{\nu_{ij}}k_{\nu_{ij}} B^R_{\nu_{ij}}(T_v)}{\displaystyle \dfrac{1}{\rm \pi}\sigma_RT_v^4}, \\ \displaystyle k^{ne}=\dfrac{\displaystyle \sum_{\nu_{ij}}k_{\nu_{ij}}B^R_{\nu_{ij}} (T_R)}{\displaystyle \dfrac{1}{\rm \pi}\sigma_RT_R^4}. \end{array}\right\} \end{equation}

From (2.14a,b), it can be seen that the absorptivity $k_{\nu _{ij}}$ is proportional to the population of molecules at vibrational energy levels $i$ and $j$, which is determined by the number density as well as the vibrational DoF of the gas. When the gray model of photons is adopted, the absorptivity $k_{\nu _{ij}}=k_{{gray}}$ is frequency independent, thus, both $k^e$ and $k^{ne}$ reduce to $k_{{gray}}$.

The macroscopic quantities of gas flow in (2.1) and the radiative field in (2.18) can be calculated based on the reduced velocity distribution functions and intensity,

(2.27) \begin{equation} \left. \begin{array}{c@{}} \displaystyle (n, n\boldsymbol{u},\boldsymbol{P})=\int_{-\infty}^{\infty} (1,\boldsymbol{v},m\boldsymbol{c}\boldsymbol{c}){f_0} \,\mathrm{d}\boldsymbol{v}, \\ \displaystyle \left(\dfrac{3}{2}k_BT_t,\dfrac{d_r}{2}k_BT_r, \dfrac{{d_v(T_v)}}{2}k_BT_v\right)=\dfrac{1}{n} \int_{-\infty}^{\infty}{\left(\dfrac{1}{2}mc^2f_0,f_1,f_2 \right)}\,\mathrm{d}\boldsymbol{v}, \\ \displaystyle (\boldsymbol{q}_t,\boldsymbol{q}_r,\boldsymbol{q}_v) =\int_{-\infty}^{\infty}{\boldsymbol{c} \left(\dfrac{1}{2}mc^2f_0,f_1,f_2\right)}\,\mathrm{d}\boldsymbol{v}, \\ \displaystyle (4\sigma_RT_R^4,\boldsymbol{q}_R)= \int_{4{\rm \pi}}{I^{R}(1,\boldsymbol{n})}\,\mathrm{d}\varOmega. \end{array}\right\} \end{equation}

2.4. Kinetic model with BCO

In the RTA operator (2.19) all molecules relax with the same speed, which is not very physical, since in general molecules with larger peculiar velocity have larger collision probability and, hence, smaller relaxation time; in fact, when (2.19) is used, the temperature of the normal shock wave will be overpredicted in the upstream (Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015b). To circumvent this problem, by observing that the elastic collision term in the first equation of (2.24) is just the Shakhov-type approximation of the BCO for monatomic gas (Shakhov Reference Shakhov1968a,Reference Shakhovb), we replace the elastic collision term $(g_{0t}-f_0)/\tau$ back with the BCO in monatomic gas,

(2.28)\begin{equation} Q_B(f_0) = \int_{-\infty}^{\infty}\int_{4{\rm \pi}}{B(\cos{\theta},{| \boldsymbol{v}-\boldsymbol{v}_* |})[f_0(\boldsymbol{v}'_*)f_0(\boldsymbol{v}') -f_0(\boldsymbol{v}_*)f_0(\boldsymbol{v})] \,\mathrm{d}{\varOmega}}\,\mathrm{d}{\boldsymbol{v}_*}, \end{equation}

so that the relaxation time depends on the molecular velocity through the collision kernel $B(\cos {\theta },{|\boldsymbol {v}-\boldsymbol {v}_* |})$ that is determined by the intermolecular potential. Note that in (2.28), $\theta$ is the deflection angle of collision, $\boldsymbol {v}$ and $\boldsymbol {v}_*$ are the pre-collision velocities of the two molecules, while $\boldsymbol {v}'$ and $\boldsymbol {v}'_*$ are the corresponding post-collision velocities. When the inverse power-law potential is considered, the collision kernel is modelled as (Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2013; Wu, Reese & Zhang Reference Wu, Reese and Zhang2014)

(2.29)\begin{equation} B=\frac{5\sqrt{{\rm \pi}{m}k_BT_0}(4k_BT_0/m)^{(2\omega -1)/2}}{64{\rm \pi}\mu(T_0)\varGamma^2(9/4-\omega/2)}\sin^{(1 -2\omega)/2}\left(\frac{\theta}{2}\right) \cos^{(1-2\omega)/2}\left(\frac{\theta}{2}\right)| \boldsymbol{v}-\boldsymbol{v}_* |^{2(1-\omega)}, \end{equation}

where $\omega$ is the viscosity index in

(2.30)\begin{equation} \mu(T)=\mu(T_0)\left(\frac{T}{T_0}\right)^\omega. \end{equation}

Meanwhile, $g_{1t}$ and $g_{2t}$ are modified correspondingly (Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2015b), resulting in the following kinetic model for molecular gas:

(2.31) \begin{equation} \left. \begin{array}{c@{}} \displaystyle \mathcal{D}f_0 = Q_B(f_0) + \dfrac{g_{0r}-g_{0t}}{Z_r\tau} + \dfrac{g_{0v}-g_{0t}}{Z_v\tau}, \\ \displaystyle \mathcal{D}f_1 = \dfrac{g_{1t}'-f_1}{\tau} + \dfrac{g_{1r}-g_{1t}}{Z_r\tau} + \dfrac{g_{1v}-g_{1t}}{Z_v\tau}, \\ \displaystyle \mathcal{D}f_2 = \dfrac{g_{2t}'-f_2}{\tau} + \dfrac{g_{2r}-g_{2t}}{Z_r\tau} + \dfrac{g_{2v}-g_{2t}}{Z_v\tau} -\dfrac{f_0}{n}\int_{4{\rm \pi}}{(k^{e}B^R(T_v)-k^{ne}I^{R})\, \mathrm{d}\varOmega}, \\ \displaystyle \dfrac{1}{c_l}\dfrac{\partial{I^{R}}}{\partial{t}} + \boldsymbol{n} \boldsymbol{\cdot} \dfrac{\partial{I^{R}}}{\partial{\boldsymbol{x}}} = k^{e}B^R(T_v)-k^{ne}I^{R}. \end{array}\right\} \end{equation}

Here

(2.32) \begin{equation} \left. \begin{array}{c@{}} \displaystyle g_{1t}'=\dfrac{d_r}{2}k_BT_r[\tau Q(f_0)+f_0]+ \dfrac{m\boldsymbol{q}_r\boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_t}E_t(T_t), \\ \displaystyle g_{2t}'=\dfrac{{d_v(T_v)}}{2}k_BT_v[\tau Q(f_0)+f_0]+ \dfrac{m\boldsymbol{q}_v\boldsymbol{\cdot}{\boldsymbol{c}}}{k_BT_t}E_t(T_t). \end{array}\right\} \end{equation}

The kinetic model with the BCO is able to distinguish the role of intermolecular potentials (Sharipov & Bertoldo Reference Sharipov and Bertoldo2009; Takata & Funagane Reference Takata and Funagane2011; Wu et al. Reference Wu, Reese and Zhang2014,  Reference Wu, Liu, Zhang and Reese2015a), while the models based on the RTA do not have this capability, but keep the advantage of computational efficiency. To be specific, the collision terms of RTA can be solved by the discrete velocity method with a computational cost of $O(N_v^3)$, and the BCO can be solved by the fast spectral method (Wu et al. Reference Wu, White, Scanlon, Reese and Zhang2013) with a cost of $O(M^2N_v^3\ln {N_v})$, where $N_v$ and $M^2$ are the number of grid points in one velocity direction and the solid angle, respectively; both the computational costs are much smaller than solving the collision operator (2.6) in the WCU equation.

2.5. Determination of parameters from relaxation properties

So far, the kinetic model equations with the RTA (2.24) and the BCO (2.31) have been established, which contain the free parameters $\boldsymbol {q}_0$, $\boldsymbol {q}_1$, $\boldsymbol {q}_2$, $Z_r$, $Z_v$ and $\tau$. In addition to the shear viscosity and translational heat conductivity in a monatomic gas, the molecular gas possesses bulk viscosity and internal thermal conductivities. The essences of these new transport coefficients are the relaxation of internal temperature and heat fluxes. Therefore, the free parameters will be determined by the recovery of relaxation rates of shear stress, energy, and heat fluxes, which correspond to the recovery of shear viscosity, bulk viscosity and thermal conductivities of all modes, respectively.

We perform the Chapman–Enskog expansion and adjust the relevant relaxation coefficients to recover the transport coefficients in the continuum flow regime, so that the kinetic model and macroscopic fluid dynamics are consistent in this regime, which is the basic requirement of kinetic modelling. In highly rarefied gas flows, the definitions of transport coefficients lose validity, while the relaxation coefficients can still be determined by involving the effect of non-equilibrium.

Since the Shakhov model and the Boltzmann equation for monatomic gas have the same shear viscosity and translational thermal conductivity, it can be shown that the second model (2.31) has the same transport coefficients as the first model (2.24). Therefore, for the sake of simplicity, we use the first model with RTA to discuss the determination of free parameters.

2.5.1. Relaxation of shear stress

Consider a spatial-homogeneous system without external acceleration. Multiplying the first equation in (2.24) by $mc^2$ and $mc_ic_j$, and integrating them with respect to $\boldsymbol {v}$, yields

(2.33) \begin{equation} \left. \begin{array}{c@{}} \displaystyle \dfrac{\partial p_t}{\partial t} = \dfrac{p_{tr}-p_t}{Z_r\tau} + \dfrac{p_{tv}-p_t}{Z_v\tau}, \\ \displaystyle \dfrac{\partial p_{ij}}{\partial t} = \dfrac{p_{t}\delta_{ij}-p_{ij}}{\tau} + \dfrac{p_{tr}\delta_{ij}-p_t\delta_{ij}}{Z_r\tau} +\dfrac{p_{tv}\delta_{ij}-p_t\delta_{ij}}{Z_v\tau}. \end{array}\right\} \end{equation}

These two equations lead to the relaxation equation for the components of the non-equilibrium stress tensor,

(2.34)\begin{equation} \frac{\partial p_{ij}}{\partial t} ={-}\frac{1}{\tau}{p_{ij}}, \end{equation}

which indicates the relaxation time of shear stress is the mean collision time due to molecular translational motion. Therefore, the shear viscosity can be recovered by adjusting the relaxation time $\tau$ as (see Appendix A)

(2.35)\begin{equation} \mu=p_t\tau. \end{equation}

The dependence of shear viscosity on mean molecular collision time in molecular gas is the same as that in monatomic gas.

2.5.2. Relaxation of energy

During the contraction or expansion of gas, the work done by pressure is converted to translational energy immediately. However, in molecular gas the molecules exhibit internal relaxation that exchanges the translational and internal energies in a finite time, which gives rise to the resistance that opposes the volume change. This is the origin of bulk viscosity in dilute gases.

The energy relaxation can be obtained by multiplying the first equation in (2.24) by $\frac {1}{2}mc^2$ and integrating the first three equations with respect to $\boldsymbol {v}$,

(2.36) \begin{equation} \left. \begin{array}{c@{}} \displaystyle \dfrac{\partial T_t}{\partial t} =\dfrac{T_{tr}-T_{t}}{Z_r\tau} + \dfrac{T_{tv}-T_{t}}{Z_v\tau}, \\ \displaystyle \dfrac{\partial T_r}{\partial t} =\dfrac{T_{tr}-T_{r}}{Z_r\tau}, \\ \displaystyle \dfrac{\partial ({{d_v(T_v)}}T_v)}{\partial t} = \dfrac{{{d_v(T_{tv})}}T_{tv}-{{d_v(T_v)}}T_{v}}{Z_v\tau}- \dfrac{2}{nk_B}\int_{4{\rm \pi}}{(k^{e}B^R(T_v)-k^{ne}I^{R})\,\mathrm{d}\varOmega}. \end{array}\right\} \end{equation}

When the radiation is absent, (2.36) reduce to the Jeans–Landau–Teller equations for the rotational and vibrational relaxation at the macroscopic level,

(2.37a,b)\begin{equation} \frac{\mathrm{d}{T_r}}{\mathrm{d}{t}} = \frac{T_t-T_r}{\tau_r}, \quad \frac{\mathrm{d}{T_v}}{\mathrm{d}t} = \frac{T_t-T_v}{\tau_v}. \end{equation}

Given the rotational and vibrational collision times $\tau _r$ and $\tau _v$, we define the rotational and vibrational collision numbers as

(2.38a,b)\begin{equation} Z_r = \frac{3\tau_r}{(3+d_r)\tau}, \quad Z_v = \frac{3\tau_v}{(3+d_v)\tau}, \end{equation}

and find that the bulk viscosity can be derived based on the Chapman–Enskog expansion (see Appendix A),

(2.39)\begin{equation} \mu_b(T_t)=2p_t\tau\frac{(3+d_r)d_rZ_r+(3+d_v)d_v Z_v}{3(3+d_r+d_v)^2}. \end{equation}

It is shown that the ratio $\mu _b/\mu$ depends only on the numbers of internal DoF and the corresponding collision numbers. Larger $Z_r$ or $Z_v$ makes the energy exchange between translational and internal motions more difficult, thus leading to higher bulk viscosity.

During the radiation transitions of the vibrational mode, the energy exchange between gas and photon alters the vibrational energy immediately, while the corresponding changes in translational and rotational modes are delayed due to the finite time of internal energy relaxation. Therefore, a new type of bulk viscosity arises, denoted as $\mu _b^{R}$, which gives the resistance of the translational–internal energy changes due to the vibrational radiation transitions. Consider the energy conservation in the radiative gas,

(2.40)\begin{equation} \frac{nk_B}{2}\left(3\frac{\partial T_t}{\partial t} +d_r\frac{\partial T_r}{\partial t}+\frac{\partial ({{d_v(T_v)}}T_v)}{\partial t}\right) ={-}\int_{4{\rm \pi}}{(k^{e}B^R(T_v)-k^{ne}I^{R})\,\mathrm{d}\varOmega}. \end{equation}

When the deviations between equilibrium temperature $T$ and $T_t, T_r, T_v$ are small, the higher-order terms of $T-T_t$ can be ignored; thus, we have

(2.41)\begin{equation} p_t-p=2p_v\frac{3\tau_v+d_r(\tau_v-\tau_r)}{(3+d_r +d_v)^2}\frac{1}{T_v}\int_{4{\rm \pi}}{(k^{e}B^R(T_v) -k^{ne}I^{R})\,\mathrm{d}\varOmega}. \end{equation}

Then the bulk viscosity resisting the vibrational radiation transitions is obtained,

(2.42)\begin{equation} \mu_b^{R}(T_v)=2p_v\tau\frac{(3+d_r)(3+d_v)Z_v-(3+d_r)d_rZ_r}{3 (3+d_r+d_v)^2}. \end{equation}

Clearly, both $\mu _b$ and $\mu _b^{R}$ are determined by the internal collision numbers $Z_r$ and $Z_v$, which describe how rapidly the equipartition of kinetic energy among different modes can be reached. Unlike shear viscosity, the determination of bulk viscosity has always been a challenging task. Nevertheless, indirect techniques in the experiments, such as the absorption and dispersion of sound waves (Graves & Argrow Reference Graves and Argrow1999), and the Rayleigh–Brillouin scattering (Pan, Shneider & Miles Reference Pan, Shneider and Miles2005; Gu & Ubachs Reference Gu and Ubachs2013), provide ways to estimate the bulk viscosity. On the other hand, MD simulations are used to extract the internal collision numbers and bulk viscosity. To be specific, the collision numbers can be determined from the adiabatic/isothermal internal relaxation processes in homogeneous MD systems (Valentini, Zhang & Schwartzentruber Reference Valentini, Zhang and Schwartzentruber2012), and the bulk viscosity can be calculated through non-equilibrium MD simulations with the expansion or compression processes (Sharma & Kumar Reference Sharma and Kumar2019).

2.5.3. Relaxation of heat flux

The rotational and vibrational modes in molecular gas carry the thermal energy and contribute also to the heat flux, while the conductance can be quite different from that of the translational one. In the continuum flow limit the total thermal conductivity determines the gas dynamics in addition to the viscosity and diffusivity. However, the thermal conductivity of a single type mode may be important and even dominant when the gas is in a non-equilibrium state. For example, the mass flow rate in thermal transpiration is found to depend on the translational thermal conductivity of gas rather than the total thermal conductivity (Mason Reference Mason1963).

In the original Rykov model, the relaxation of translational heat flux is independent of the rotational one, and vice versa. However, due to the energy exchange between different modes, it is necessary to consider the fact that the relaxations of heat fluxes are coupled within all the DoF. In general, according to the WCU equation, the relaxation of translational, rotational and vibrational heat fluxes satisfies the following relation in a spatially homogeneous system (Mason & Monchick Reference Mason and Monchick1962):

(2.43) \begin{equation} \left[ \begin{array}{@{}c@{}} \partial{\boldsymbol{q}_{t}}/{\partial{t}} \\ \partial{\boldsymbol{q}_{r}}/{\partial{t}} \\ \partial{\boldsymbol{q}_{v}}/{\partial{t}} \end{array} \right] ={-}\frac{p_t}{\mu} \left[ \begin{array}{@{}ccc@{}} A_{tt} & A_{tr} & A_{tv} \\ A_{rt} & A_{rr} & A_{rv} \\ A_{vt} & A_{vr} & A_{vv} \end{array} \right] \left[ \begin{array}{@{}c@{}} \boldsymbol{q}_{t} \\ \boldsymbol{q}_{r} \\ \boldsymbol{q}_{v} \end{array} \right]. \end{equation}

Here the dimensionless relaxation rate $\boldsymbol {A}$ is a $3\times 3$ matrix encapsulating the dimensionless thermal relaxation rates. Accordingly, $\boldsymbol {q}_0$, $\boldsymbol {q}_1$, $\boldsymbol {q}_2$ in reference distributions (2.20) can be determined in terms of $\boldsymbol {q}_t$, $\boldsymbol {q}_r$, $\boldsymbol {q}_v$ and the thermal relaxation rates $\boldsymbol {A}$. To be specific, the first three equations in (2.24) are multiplied by $\frac {1}{2}mc^2\boldsymbol {c}$, $\boldsymbol {c}$ and $\boldsymbol {c}$, respectively, and then are integrated with respect to the molecular velocity $\boldsymbol {v}$, yielding

(2.44)\begin{equation} \begin{bmatrix} \boldsymbol{q}_{0} \\ \boldsymbol{q}_{1} \\ \boldsymbol{q}_{2} \end{bmatrix} = \begin{bmatrix} (2-3A_{tt})Z_{int}+1 & -3A_{tr}Z_{int} & -3A_{tv}Z_{int} \\ -A_{rt}Z_{int} & -A_{rr}Z_{int}+1 & -A_{rv}Z_{int} \\ -A_{vt}Z_{int} & -A_{vr}Z_{int} & -A_{vv}Z_{int}+1 \end{bmatrix} \begin{bmatrix} \boldsymbol{q}_{t} \\ \boldsymbol{q}_{r} \\ \boldsymbol{q}_{v} \end{bmatrix}, \end{equation}

where $Z_{int}=({1}/{Z_r}+{1}/{Z_v})^{-1}$.

When the gas flow is stationary, the thermal relaxation rates are related to the translational, rotational and vibrational thermal conductivities, $\kappa _t$, $\kappa _r$ and $\kappa _v$, respectively (see Appendix A),

(2.45) \begin{equation} \left[ \begin{array}{@{}c@{}} \kappa_t \\ \kappa_r \\ \kappa_v \end{array} \right] = \frac{k_B\mu}{2m} \left[ \begin{array}{@{}ccc@{}} A_{tt} & A_{tr} & A_{tv} \\ A_{rt} & A_{rr} & A_{rv} \\ A_{vt} & A_{vr} & A_{vv} \end{array} \right]^{{-}1} \left[ \begin{array}{@{}c@{}} 5 \\ d_r \\ d_v(T_v) \end{array} \right]. \end{equation}

It is shown that the presence of radiation transition does not affect thermal conductivities.

It will be convenient to use the following dimensionless (Eucken Reference Eucken1913) factor $f_{eu}$:

(2.46)\begin{equation} c_vf_{eu}\equiv\frac{\kappa}{\mu}=\frac{\kappa_t+\kappa_r+\kappa_v}{\mu}. \end{equation}

Here $\kappa$ is the total thermal conductivity and $c_v$ is the specific heat capacity at constant volume. Similarly, $f_t$, $f_r$ and $f_v$ represent the Eucken factors of the translational, rotational and vibrational modes, respectively,

(2.47a–c)\begin{equation} f_t = \frac{2}{3}\frac{m\kappa_t}{k_B\mu}, \quad f_r = \frac{2}{d_r}\frac{m\kappa_r}{k_B\mu}, \quad f_v = \frac{2}{d_v}\frac{m\kappa_v}{k_B\mu}. \end{equation}

Therefore, the Eucken factors are determined by the thermal relaxation rates as

(2.48) \begin{equation} \left[ \begin{array}{@{}c@{}} f_t \\ f_r \\ f_v \end{array} \right] = \left[ \begin{array}{@{}ccc@{}} 3A_{tt} & d_rA_{tr} & d_v(T_v)A_{tv} \\ 3A_{rt} & d_rA_{rr} & d_v(T_v)A_{rv} \\ 3A_{vt} & d_rA_{vr} & d_v(T_v)A_{vv} \end{array} \right]^{{-}1} \left[ \begin{array}{@{}c@{}} 5 \\ d_r \\ d_v(T_v) \end{array} \right]. \end{equation}

Clearly, the elements in matrix $\boldsymbol {A}$ cannot be fully determined even though all the Eucken factors $f_t, f_r, f_v$ (thermal conductivities $\kappa _t, \kappa _r, \kappa _v$ equivalently) are fixed. In other words, in molecular gas, having all the transport coefficients is not enough to exactly describe the relaxation of heat flux. Therefore, it is necessary to recover the thermal relaxation rates in the kinetic model correctly. Nevertheless, since the rate of translation–vibration energy exchange is usually much slower than that of translation–rotational energy exchange, the values of $A_{tv}, A_{vt}, A_{rv}, A_{vr}$ in matrix $\boldsymbol {A}$ are much smaller than the others and, thus, can be approximated to be zero practically (Mason & Monchick Reference Mason and Monchick1962). While the off-diagonal elements $A_{tr}, A_{rt}$ still have to be correctly recovered, which have been found to make uncertainty in predicting macroscopic gas dynamics (Li et al. Reference Li, Zeng, Su and Wu2021). We propose the method for extracting the thermal relaxation rates from MD simulations by monitoring the evolution of heat fluxes, which are initially generated with prescribed values by manipulating the velocities of all atoms in the systems (see the details in Appendix B).

2.6. The H theorem

With all the parameters determined by the relaxation properties, the entropy of both gas molecules and photons can be evaluated. We define the functional $\mathcal {H}$ according to Groppi & Spiga (Reference Groppi and Spiga1999),

(2.49)\begin{equation} \mathcal{H}(f_i,I_{\nu_{ij}}^{R}) = \sum_{i}^{N}\int_{0}^{\infty}\int_{-\infty}^{\infty} \ln(f_i){f_i}\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}I_r +\sum_{i,j}^{N}\int_{4{\rm \pi}}\ln(I_{\nu_{ij}}^{R'}){I_{\nu_{i j}}^{R'}}\,\mathrm{d}\varOmega, \end{equation}

where $I_{\nu _{ij}}^{R'}=B_{ij}I_{\nu _{ij}}^{R}/(A_{ji}+B_{ji}I_{\nu _{ij}}^{R})$. Then, with the kinetic model equations, the change of functional $\mathcal {H}$ can be obtained as

(2.50)\begin{align} \frac{\mathrm{d}\mathcal{H}}{\mathrm{d}t} &= \underbrace{\sum_{i}^{N}\iint\ln(f_i)J_{{gas},i}\, \mathrm{d}\boldsymbol{v}\,\mathrm{d}I_r}_{\left(\dfrac{\mathrm{d} \mathcal{H}}{\mathrm{d}t}\right)_{1}} \nonumber\\ &\quad +\underbrace{\sum_{i}^{N}\iint\ln(f_i)J_{{photon},i}\, \mathrm{d}\boldsymbol{v}\,\mathrm{d}I_r+\sum_{i,j}^{N} \int_{4{\rm \pi}}\ln\left(\frac{B_{ij}I_{\nu_{ij}}^{R}}{A_{ji} +B_{ji}I_{\nu_{ij}}^{R}}\right){J_{\nu_{ij}}^{R}}\, \mathrm{d}\varOmega}_{\left(\dfrac{\mathrm{d} \mathcal{H}}{\mathrm{d}t}\right)_{2}}. \end{align}

The term $({\mathrm {d}\mathcal {H}}/{\mathrm {d}t})_{1}$ indicates the change of $\mathcal {H}$ due to the intermolecular interactions. It should be noted that as long as the simplifications of the gas–gas collision operators are adopted, the evolution of entropy governed by the model equations is no longer the same as that of the WCU equation. As the extensions of the Shakhov model for monatomic gas and the Rykov model for diatomic gas, our kinetic models share similar properties in the evolution of entropy. In the literature, the Boltzmann H theorem for the Shakhov model has been proved in the linearized system (Shakhov Reference Shakhov1968b), while it has not been discussed for the Rykov model. We explicitly evaluate $({\mathrm {d}\mathcal {H}}/{\mathrm {d}t})_{1}$ for the RTA collision operator in a linearized molecular gas system that slightly deviates from an equilibrium state (see Appendix C),

(2.51) \begin{align} \left(\frac{\mathrm{d}{\mathcal{H}}}{\mathrm{d}{t}}\right)_{1} &={-}\frac{4}{{\tau}p_0^2v_m^2} \left[ \begin{array}{@{}c@{}} {\rm \Delta}\boldsymbol{q}_t \\ {\rm \Delta}\boldsymbol{q}_r \\ {\rm \Delta}\boldsymbol{q}_v \end{array} \right]^{{T}} \left[ \begin{array}{@{}ccc@{}} A_{tt}/5 & A_{tr}/5 & A_{tv}/5 \\ A_{rt}/d_r & A_{rr}/d_r & A_{rv}/d_r \\ A_{vt}/d_v & A_{vr}/d_v & A_{vv}/d_v \end{array} \right] \left[ \begin{array}{@{}c@{}} {\rm \Delta}\boldsymbol{q}_t \\ {\rm \Delta}\boldsymbol{q}_r \\ {\rm \Delta}\boldsymbol{q}_v \end{array} \right] \nonumber\\ & \quad-\frac{1}{{\tau}}\left(\frac{1}{p_0^2}{\rm \Delta} \boldsymbol{P}:{\rm \Delta}\boldsymbol{P}+\sum({HoT})^2\right), \end{align}

where ${\rm \Delta} \boldsymbol {q}_t, {\rm \Delta} \boldsymbol {q}_r, {\rm \Delta} \boldsymbol {q}_v$ are the perturbed heat fluxes, ${\rm \Delta} \boldsymbol {P}$ is the perturbed stress tensor and ${HoT}$ is the high-order term. The thermal relaxation rates $\boldsymbol {A}$ are positive definite, since the diagonal elements are positive and much larger than the magnitudes of off-diagonal elements. Therefore, it is demonstrated that $\mathcal {H}$ always decreases with time due to the gas–gas collisions in the linear system. Moreover, when the elastic collision operator is changed to BCO, $({\mathrm {d}\mathcal {H}}/{\mathrm {d}t})_{1}\leq 0$ still keeps valid.

In a nonlinear system, because the reference distribution functions constructed in our kinetic models only involve orthogonal polynomials in terms of the heat fluxes, the changes of functional $\mathcal {H}$ due to the shear stress and all the high-order terms are exactly the same as those governed by the standard BGK equation, which have been proven to be non-positive. Additionally, by recovering the correct thermal relaxation rates and thermal conductivities in the model equations, the change of functional $\mathcal {H}$ due to the heat flux should be the same as that of the Navier–Stokes–Fourier equations in the continuum limit. Also, the change of functional $\mathcal {H}$ is zero in the free molecular limit. Therefore, it is very likely that $({\mathrm {d}\mathcal {H}}/{\mathrm {d}t})_{1}\leq 0$ holds in all the flow regimes, though it is not explicitly calculated or rigorously proved. On the other hand, although the H theorem has been proven in the standard BGK and the ES-BGK models, the Shakhov and Rykov models perform better in many problems (Fei et al. Reference Fei, Liu, Liu and Zhang2020; Zeng et al. Reference Zeng, Li and Wu2022). Therefore, to achieve higher accuracy, we still choose to construct the gas–gas interaction terms based on the Rykov model.

The term $({\mathrm {d}\mathcal {H}}/{\mathrm {d}t})_{2}$ indicates the change of $\mathcal {H}$ due to the gas–photon interactions. By substituting the expressions of $J_{{photon},i}$ (2.11) and $J_{\nu _{ij}}^{R}$ (2.14a,b) and (2.15) into $({\mathrm {d}\mathcal {H}}/{\mathrm {d}t})_{2}$, we have

(2.52)\begin{align} \left(\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}t}\right)_{2} &= \sum_{i,j}^{N}\iiint{[A_{ji}f_i+(B_{ji}f_j-B_{ij}f_i)I_{\nu_{ij}}^{R}] \left[\ln(f_i)-\ln(f_j)+\ln\left(\frac{B_{ij}I_{\nu_{ij}}^{R}}{A_{ji} +B_{ji}I_{\nu_{ij}}^{R}}\right)\right]}\,\mathrm{d}\boldsymbol{v} \,\mathrm{d}I_r\mathrm{d}\varOmega \nonumber\\ &=\sum_{i,j}^{N}\iiint{(A_{ji}f_j+B_{ji}f_jI_{\nu_{ij}}^{R}) \left[1-\frac{B_{ij}f_iI_{\nu_{ij}}^{R}}{A_{ji}f_j +B_{ji}f_jI_{\nu_{ij}}^{R}}\right]\ln\left(\frac{B_{ij}f_i I_{\nu_{ij}}^{R}}{A_{ji}f_j+B_{ji}f_jI_{\nu_{ij}}^{R}} \right)}\,\mathrm{d}\boldsymbol{v}\,\mathrm{d}I_r\,\mathrm{d}\varOmega. \end{align}

Because of the positivity of Einstein coefficients and the convexity of the $\ln$ function, the gas–photon interactions always lead to the decrease of $\mathcal {H}$ with time $({\mathrm {d}\mathcal {H}}/{\mathrm {d}t})_{2}\leq 0$.

2.7. Dimensionless expressions

Let $L_0, T_0, n_0$ be the reference length, temperature and number density, respectively, then the most probable speed is $v_m=\sqrt {2k_BT_0/m}$ and the reference pressure is $p_0=n_0k_BT_0$. The dimensionless variables are introduced as

(2.53) \begin{equation} \left. \begin{array}{c@{}} \tilde{\boldsymbol{x}}=\boldsymbol{x}/L_0, \quad \tilde{T}=T/T_0, \quad \tilde{n}=n/n_0, \quad \tilde{t}=v_mt/L_0, \\ \tilde{\boldsymbol{v}}=\boldsymbol{v}/v_m, \quad \tilde{\boldsymbol{c}}=\boldsymbol{c}/v_m, \quad \tilde{p}=p/p_0, \quad \tilde{\boldsymbol{q}}=\boldsymbol{q}/(p_0v_m), \\ \tilde{f}_0=v_m^{3}f_0/n_0, \quad \tilde{f}_1=v_m^{3}f_1/p_0, \quad \tilde{f}_2=v_m^{3}f_2/p_0, \quad \tilde{I}^R=I^R/(p_0v_m). \end{array}\right\} \end{equation}

When the gray model is used for photon transport, the Knudsen numbers for gas flow and photon transport are defined, respectively, as

(2.54a,b)\begin{equation} {Kn}_{{gas}}=\frac{\mu(T_0)}{n_0L_0}\sqrt{\frac{\rm \pi}{2mk_BT_0}}, \quad {Kn}_{{photon}}=\frac{1}{k_{{gray}}L_0}. \end{equation}

The relative strength of the radiative to the convective heat transfer is given by the dimensionless parameter

(2.55)\begin{equation} \tilde{\sigma}_R=\frac{\sigma_RT_0^3}{n_0k_Bv_m}, \end{equation}

which is equivalent to the reciprocal of the Boltzmann number commonly used in radiation hydrodynamics (Casto Reference Casto2004). Therefore, the model equations are non-dimensionalized as

(2.56) \begin{equation} \left. \begin{array}{c@{}} \displaystyle \tilde{\mathcal{D}}\widetilde{f_0} = \tilde{Q}(\,\widetilde{f_0}) + \dfrac{\sqrt{\rm \pi}\tilde{n} \tilde{T}_t^{1-\omega}}{2{Kn}_{{gas}}}\left[\dfrac{\tilde{g}_{0r}- \tilde{g}_{0t}}{Z_r} + \dfrac{\tilde{g}_{0v}-\tilde{g}_{0t}}{Z_v}\right], \\ \displaystyle \tilde{\mathcal{D}}\widetilde{f_1} = \dfrac{\sqrt{\rm \pi}\tilde{n} \tilde{T}_t^{1-\omega}}{2{Kn}_{{gas}}}\left[(\tilde{g}_{1t}-\tilde{f}_1) +\dfrac{\tilde{g}_{1r}-\tilde{g}_{1t}}{Z_r} + \dfrac{\tilde{g}_{1v}- \tilde{g}_{1t}}{Z_v}\right], \\ \displaystyle \tilde{\mathcal{D}}\widetilde{f_2} = \dfrac{\sqrt{\rm \pi}\tilde{n} \tilde{T}_t^{1-\omega}}{2{Kn}_{{gas}}}\left[(\tilde{g}_{2t}-\tilde{f}_2) +\dfrac{\tilde{g}_{2r}-\tilde{g}_{2t}}{Z_r} + \dfrac{\tilde{g}_{2v} -\tilde{g}_{2t}}{Z_v}\right]\\ \hspace{-2pc}-\,\dfrac{\tilde{f}_0}{\tilde{n} {Kn}_{{photon}}}{\left(4\tilde{\sigma}_R\tilde{T}_v^4- \displaystyle\int_{4\pi}{\tilde{I}^{R}\,\mathrm{d}\varOmega}\right)}, \\ \displaystyle \dfrac{1}{\tilde{c}_l}\dfrac{\partial{\tilde{I}^{R}}}{\partial{\tilde{t}}} +\boldsymbol{n} \boldsymbol{\cdot} \dfrac{\partial{\tilde{I}^{R}}}{\partial {\tilde{\boldsymbol{x}}}} = \dfrac{1}{{Kn}_{{photon}}} \left(\dfrac{1}{\rm \pi}\tilde{\sigma}_R\tilde{T}_v^4-\tilde{I}^{R}\right), \end{array}\right\} \end{equation}

where $\tilde {\mathcal {D}}\widetilde {f_i}= {\partial {\widetilde {f_i}}}/{\partial {\tilde {t}}}+\boldsymbol {\tilde {v}} \boldsymbol {\cdot } ({\partial {\widetilde {f_i}}}/{\partial {\tilde {\boldsymbol {x}}}})+ {\partial {(\tilde {\boldsymbol {a}}\widetilde {f_i})}}/{\partial {\tilde {\boldsymbol {v}}}}$, $\tilde {c}_l=c_l/v_m$ does not affect the results of steady-state problems, and $\tilde {Q}(\,\tilde{f}_0)$ represents the dimensionless elastic collision operators

(2.57) \begin{equation} \tilde{Q}(\,\tilde{f}_0) = \left\{ \begin{array}{@{}l} \dfrac{\sqrt{\rm \pi}\tilde{n}\tilde{T}_t^{1-\omega}}{2{Kn}_{{gas}}} (\tilde{g}_{0r}-\tilde{g}_{0t}), \quad (\text{RTA}) \\ \begin{aligned} & \dfrac{1}{{Kn}_{{gas}}}\dfrac{5}{2^{7-\omega}\varGamma^2(9/4- \omega/2)}\iint\sin^{(1-2\omega)/2}\left(\dfrac{\theta}{2}\right) \cos^{(1-2\omega)/2}\left(\dfrac{\theta}{2}\right) \\ & \quad \times| \tilde{\boldsymbol{v}}-\tilde{\boldsymbol{v}}_* |^{2(1-\omega)} \boldsymbol{\cdot}[\,\tilde{f}_0(\tilde{\boldsymbol{v}}'_*)\tilde{f}_0(\tilde{\boldsymbol{v}}') -\tilde{f}_0(\tilde{\boldsymbol{v}}_*)\tilde{f}_0 (\tilde{\boldsymbol{v}})]\,\mathrm{d}{\varOmega}\, \mathrm{d}{\tilde{\boldsymbol{v}}_*}. \quad (\text{BCO}) \end{aligned} \end{array}\right. \end{equation}

The dimensionless reduced reference velocity distribution functions are

(2.58) \begin{equation} \left.\begin{array}{c} \displaystyle \tilde{g}_{0t}=\tilde{n}\left(\dfrac{1}{{\rm \pi}\tilde{T}_t}\right)^{3/2} \exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_t}\right)} \left[1+\dfrac{4\tilde{\boldsymbol{q}}_t\boldsymbol{\cdot}{\tilde{\boldsymbol{c}}}}{15 {\tilde{T}_t}{\tilde{p}_t}} \left(\dfrac{\tilde{c}^2}{\tilde{T}_t} -\dfrac{5}{2}\right)\right], \\ \displaystyle \tilde{g}_{0r}=\tilde{n}\left(\dfrac{1}{{\rm \pi}\tilde{T}_{tr}} \right)^{3/2}\exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_{tr}} \right)}\left[1+\dfrac{4\tilde{\boldsymbol{q}}_0\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{15{\tilde{T}_{tr}}{\tilde{p}_{tr}}} \left(\dfrac{\tilde{c}^2}{\tilde{T}_{tr}} -\dfrac{5}{2}\right)\right], \\ \displaystyle \tilde{g}_{0v}=\tilde{n}\left(\dfrac{1}{{\rm \pi}\tilde{T}_{tv}} \right)^{3/2}\exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_{tv}} \right)}\left[1+\dfrac{4\tilde{\boldsymbol{q}}_0\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{15{\tilde{T}_{tv}}{\tilde{p}_{tv}}} \left(\dfrac{\tilde{c}^2}{\tilde{T}_{tv}}-\dfrac{5}{2} \right)\right], \\ \tilde{g}_{1t}= \left\{ \begin{aligned} & \displaystyle \dfrac{d_r}{2}\tilde{T}_r\tilde{g}_{0t}+ \left(\dfrac{1}{{\rm \pi}\tilde{T}_t}\right)^{3/2}\exp{\left(- \dfrac{\tilde{c}^2}{\tilde{T}_t}\right)}\dfrac{2 \tilde{\boldsymbol{q}}_r\boldsymbol{\cdot}{\tilde{\boldsymbol{c}}}}{\tilde{T}_t}, \quad (\text{RTA})\\ & \displaystyle \dfrac{d_r}{2}\tilde{T}_r\left[\dfrac{2{Kn}_{{gas}}}{\sqrt{\rm \pi} \tilde{n}\tilde{T}_t^{1-\omega}} \tilde{Q}_B(\,\tilde{f}_0) +\tilde{f}_0\right]+\left(\dfrac{1}{{\rm \pi}\tilde{T}_t} \right)^{3/2}\exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_t} \right)}\dfrac{2\tilde{\boldsymbol{q}}_r\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{\tilde{T}_t}, \quad (\text{BCO}) \end{aligned}\right. \\ \displaystyle \tilde{g}_{1r}=\dfrac{d_r}{2}\tilde{T}_{tr}\tilde{g}_{0r}+ \left(\dfrac{1}{{\rm \pi}\tilde{T}_{tr}}\right)^{3/2} \exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_{tr}}\right)} \dfrac{2\tilde{\boldsymbol{q}}_1\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{\tilde{T}_{tr}}, \\ \displaystyle \tilde{g}_{1v}=\dfrac{d_r}{2}\tilde{T}_{r}\tilde{g}_{0v}+ \left(\dfrac{1}{{\rm \pi}\tilde{T}_{tv}}\right)^{3/2} \exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_{tv}}\right)} \dfrac{2\tilde{\boldsymbol{q}}_1\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{\tilde{T}_{tv}}, \\ \tilde{g}_{2t}= \left\{ \begin{aligned} & \dfrac{d_v(\tilde{T}_v)}{2}\tilde{T}_v\tilde{g}_{0t}+ \left(\dfrac{1}{{\rm \pi}\tilde{T}_t}\right)^{3/2}\exp {\left(-\dfrac{\tilde{c}^2}{\tilde{T}_t}\right)} \dfrac{2\tilde{\boldsymbol{q}}_v\boldsymbol{\cdot}{\tilde{\boldsymbol{c}}}}{\tilde{T}_t}, \quad (\text{RTA})\\ & \dfrac{d_v(\tilde{T}_v)}{2}\tilde{T}_v\left[\dfrac{2{Kn}_{{gas}}}{\sqrt{\rm \pi} \tilde{n}\tilde{T}_t^{1-\omega}} \tilde{Q}_B(\,\tilde{f}_0) +\tilde{f}_0\right]+\left(\dfrac{1}{{\rm \pi}\tilde{T}_t} \right)^{3/2}\exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_t} \right)}\dfrac{2\tilde{\boldsymbol{q}}_v\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{\tilde{T}_t}, \quad (\text{BCO}) \end{aligned}\right. \\ \displaystyle \tilde{g}_{2r}=\dfrac{{d_v(\tilde{T}_v)}}{2}\tilde{T}_{v} \tilde{g}_{0r}+\left(\dfrac{1}{{\rm \pi}\tilde{T}_{tr}}\right)^{3/2} \exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_{tr}}\right)} \dfrac{2\tilde{\boldsymbol{q}}_2\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{\tilde{T}_{tr}}, \\ \displaystyle \tilde{g}_{2v}=\dfrac{{d_v(\tilde{T}_{tv})}}{2}\tilde{T}_{tv} \tilde{g}_{0v}+\left(\dfrac{1}{{\rm \pi}\tilde{T}_{tv}}\right)^{3/2} \exp{\left(-\dfrac{\tilde{c}^2}{\tilde{T}_{tv}}\right)} \dfrac{2\tilde{\boldsymbol{q}}_2\boldsymbol{\cdot} {\tilde{\boldsymbol{c}}}}{\tilde{T}_{tv}}. \end{array}\right\} \end{equation}

When the gas species and relaxation rates are fixed, the solutions of (2.56) are determined by four dimensionless parameters: ${Kn}_{{gas}}$, ${Kn}_{{photon}}$, $\tilde {\sigma }_R$ and $T_0/T_{{ref}}$, where $\tilde {\sigma }_R$ gives the relative importance of radiative heat transfer in the system, and $T_0/T_{{ref}}$ indicates the degree of vibrational excitation. From the term of energy exchange between vibrational modes and the radiative field in the governing equations, it can be seen that when $\tilde {\sigma }_R/{Kn}_{{photon}}\rightarrow 0$, the interaction between gas and photon becomes negligible.

3.1. Relaxation rates extracted from DSMC

In the validation of the proposed kinetic models, we need to first compare the results of our models with the reference solutions using the same physical properties and system parameters to demonstrate that the model equations can give an accurate evolution of the system when all the necessary parameters have been known. In DSMC simulation for molecular gas flow, generally, not all the transport coefficients of a real gas can be recovered simultaneously in this method. However, the results of DSMC still can be regarded as reference solutions when we consider a virtual gas with the exact same relaxation rates realized in DSMC simulations. Therefore, we extract the thermal relaxation rates from the DSMC and apply them to our kinetic model to make a fair comparison. With the fixed shear viscosity and self-diffusion coefficient, the collision number $Z_r$ and $Z_v$ are the only parameters that affect the thermal relaxation rates in DSMC. As an example, we take the temperature-independent collision number $Z_r=2.667$, which is determined in DSMC by matching the experimentally measured thermal conductivity of nitrogen at room temperature (Li et al. Reference Li, Zeng, Su and Wu2021), and choose $Z_v=10Z_r$ as a typical value of vibrational collision number, as it is usually much larger than the rotational one. For the cases that $Z_v$ has an even higher order of magnitude, the vibrational modes can be treated as frozen and, thus, no longer need to be coupled with the other modes into the kinetic modelling.

Similar to the procedure of extracting thermal relaxation rates for the translational and rotational DoF from DSMC (Li et al. Reference Li, Zeng, Su and Wu2021), here a homogeneous system of nitrogen is simulated, which consists of $10^6$ simulation particles in a cubic cell of the volume $(10\ \textrm {nm})^3$. The periodic condition is applied at all boundaries. Binary collisions are described by the variable-soft-sphere model, and the system parameters and properties of nitrogen used in the simulations are $d_r=d_v=2$, $n_0=2.69\times 10^{25}\ \textrm {m}^{-3}$, $T_0=5000\ \textrm {K}$, $m=4.65\times 10^{-26}\ \textrm {kg}$, the molecular diameter is $d=4.11\times 10^{-10}\ \textrm {m}$, the viscosity index is $\omega =0.74$, the angular scattering parameter is $\alpha =1.36$ and the Schmidt number is $Sc=1/1.34$ (Bird Reference Bird1994). Initially, simulation particles with positive velocity in the $x_1$ direction follow the equilibrium distribution at $4500$ K, while those moving in the opposite direction follow the equilibrium distribution at $5500$ K, see figures 1(a) and 1(b), so that initial heat fluxes in all DoF are generated. Then the evolution of heat flux is monitored until the entire system reaches thermal equilibrium; see figure 1(c). An ensemble average is taken from 3000 independent runs to get the time derivative of heat flux in figure 1(d). Finally, the following relaxation rates are extracted by solving the linear regression problem (2.43) with the least squares method,

(3.1)\begin{equation} \left[ \begin{array}{ccc} A_{tt} & A_{tr} & A_{tv} \\ A_{rt} & A_{rr} & A_{rv} \\ A_{vt} & A_{vr} & A_{vv} \end{array} \right] = \left[ \begin{array}{ccc} 0.786 & -0.208 & 0.003 \\ -0.047 & 0.883 & -0.049 \\ -0.004 & -0.038 & 0.772 \end{array} \right]. \end{equation}

Hence, according to (2.48), we have $f_t=2.3635$, $f_r=1.3979$, $f_v=1.3825$ and $f_{eu}=1.807$. With these parameters, our kinetic model is uniquely determined.

Figure 1. Extraction of the thermal relaxation rates $\boldsymbol {A}$ in (2.43) from the DSMC simulation. Special distributions of (a) the molecular velocity and (b) rotational/vibrational energy (overlap with each other) are designed to generate an initial heat flux. (c) The evolution of heat fluxes and (d) their time derivatives are monitored until the system reaches thermal equilibrium.

3.2. Fourier flow

The heat transfer in the nitrogen gas between two parallel plates located at $x_2=0$ and $L_0$ are considered, where the temperature of the lower and upper plates are $T_w=T_l$ and $T_w=T_u$, respectively. The averaged number density of nitrogen is set to be $n_0$, and the characteristic length $L_0$ is chosen to be the distance between two plates. The Maxwell boundary conditions for molecular gas with a single accommodation coefficient $\alpha$ are adopted (Larina & Rykov Reference Larina and Rykov1986), so that the reflected distributions are

(3.2) \begin{align} \left. \begin{array}{c@{}} \displaystyle x_2=0, v_2\ge0: \quad f_0^{+} = \dfrac{n^{-}(x_2=0)}{n_0}E_t(T^{+}_t), \ f_1^{+}=\dfrac{d_r}{2}k_BT^{+}_rf_0^{+}, \ f_2^{+}=\dfrac{d_v}{2}k_BT^{+}_vf_0^{+}, \\ \displaystyle x_2=L_0, v_2\le0: \quad f_0^{+} = \dfrac{n^{-}(x_2=L_0)}{n_0}E_t(T^{+}_t), \ f_1^{+}=\dfrac{d_r}{2}k_BT^{+}_rf_0^{+}, \ f_2^{+}=\dfrac{d_v}{2}k_BT^{+}_vf_0^{+}, \end{array}\right\} \end{align}

where the superscripts $-$ and $+$ represent the incident and reflected molecules, respectively; $n^{-}$ is determined by the flux of incident number density of gas at the plates

(3.3)\begin{equation} n^{-}= \left|\left(\frac{2m{\rm \pi}}{k_BT_w}\right)^{1/2} \int_{v_n<0}v_2\,f^{-}_0\,\mathrm{d}\boldsymbol{v}\right|, \end{equation}

where $v_n$ is the molecular velocity component normal to the wall ($v_n=v_2$ at $x_2=0$ and $v_n=-v_2$ at $x_2=L_0$); $T^{-}_{t},T^{-}_{r},T^{-}_{v}$ are the temperatures of the incident molecules,

(3.4)\begin{equation} T^{-}_{t}=\frac{1}{2}\frac{\displaystyle{\int_{v_n<0} \frac{1}{2}mc^2v_n\,f_0^{-}\,\mathrm{d} \boldsymbol{v}}}{\displaystyle{\int_{v_n<0}v_n\,f_0^{-}\, \mathrm{d}\boldsymbol{v}}}, T^{-}_{r}=\frac{2}{d_r}\frac{\displaystyle{\int_{v_n<0}v_n\,f_ 1^{-}\,\mathrm{d}\boldsymbol{v}}}{\displaystyle{\int_{v_n<0} v_n\,f_0^{-}\,\mathrm{d}\boldsymbol{v}}}, T^{-}_{v}=\frac{2}{d_v}\frac{\displaystyle{\int_{v_n<0} v_n\,f_2^{-}\,\mathrm{d}\boldsymbol{v}}}{\displaystyle{\int_{v_n<0} v_n\,f_0^{-}\,\mathrm{d}\boldsymbol{v}}}. \end{equation}

Then the temperatures of reflected molecules $T^{+}_{t},T^{+}_{r},T^{+}_{v}$ are the linear combination of $T_w$ and $T^{-}_{t},T^{-}_{r},T^{-}_{v}$ as $T^{+}_{t}=\alpha T_w+(1-\alpha )T^{-}_{t}$, $T^{+}_{r}=\alpha T_w+(1-\alpha )T^{-}_{r}$ and $T^{+}_{v}=\alpha T_w+(1-\alpha )T^{-}_{v}$ (Larina & Rykov Reference Larina and Rykov1986).

3.2.1. Comparison with DSMC

In the comparison with DSMC simulations, the Knudsen numbers considered are ${Kn}_{{gas}}=0.1$ and 1. The temperature of the lower and upper plates are $T_l=0.8T_0$ and $T_u=1.2T_0$, respectively, and the accommodation coefficient $\alpha =1$. Numerical solutions of the kinetic models I and II, as well as the DSMC results, are shown in figure 2. For both ${Kn}_{{gas}}=0.1$ and ${Kn}_{{gas}}=1$, excellent agreement in the density, temperature and heat flux are observed. Meanwhile, profiles of translational, rotational and vibrational temperatures nearly overlap, although the relaxation times for different modes are different. Additionally, the rotational and vibrational heat fluxes are almost the same (figure 2f), due to the close values of the rotational and vibrational thermal conductivities. Thus, it is clearly seen that the values of collision numbers $Z_r$ and $Z_v$ do not have an influence on the distribution of macroscopic quantities for the steady-state planar Fourier flow.

Figure 2. Comparisons of the (a) density, (b) translational temperature, (c) rotational temperature,(d) vibrational temperature, (e) translational heat flux and (f) rotational/vibrational heat flux of nitrogen between kinetic model I (green lines), kinetic model II (red lines) and DSMC (blue circles) for the Fourier flows.

3.2.2. Comparison with experiments

Teagan & Springer (Reference Teagan and Springer1968) performed the experimental measurements of heat transfer and density distributions of argon/nitrogen between two parallel plates in the transition regime. The heat flux was determined by measuring the total electrical power input into the hot plate and subtracting the energy losses. The density profiles were measured by observing the luminescence produced by a high-energy electron beam traversed between the plates. In the experiments, the temperature of the cold plate is maintained at 288 K, while the temperature of the hot one is 296 K in the heat transfer measurements, and it is 368 K in the density distribution measurements to increase the accuracy. The accommodation coefficient was determined from the heat transfer measurement at low pressure (free molecular condition), where Knudsen's formula is applicable, and it gave the value $\alpha =0.76$ for nitrogen.

Both kinetic models I and II are used to predict the heat flux and density profiles of nitrogen at the experimental conditions. Since the temperature difference of the two plates is relatively small, the constant values of $\alpha =0.76$, $Z_r=2.667$ are used ($\alpha$ is obtained in experiments, and $Z_r$ is determined by matching the thermal conductivity of nitrogen at room temperature). The vibrational DoF is negligibly small in these conditions.

Figure 3 compares the results obtained by kinetic models and experimental data in a wide range of ${Kn}$. The total heat flux $q_{tot}$ in figure 3(a) is normalized by the value $q_{fm}$ at the free molecular limit, which is calculated at ${Kn}=1000$. The number density in figures 3(b) and 3(c) is normalized by the values at the midpoint between the plates. Good agreements between the numerical results and the experimental data are achieved. It is found that the maximum discrepancy in density profiles is less than 6 %, which is expected to be further reduced by using a more accurate gas–surface interaction model if accommodation coefficients of translational and rotational modes can be measured separately (Larina & Rykov Reference Larina and Rykov1986).

Figure 3. Comparisons of the normalized (a) total heat flux, (b,c) density distribution of nitrogen between kinetic model I (green squares/lines), kinetic model II (red diamonds/lines) and experimental data (Teagan & Springer Reference Teagan and Springer1968) (blue circles) for the planar heat transfer flows.

3.3. Couette flow

The configuration of the Couette flow is the same as the Fourier flow, while the temperatures of both plates are kept the same at $T_0$, and the speed of the lower and upper plates are $u_1=-v_m$ and $u_2=v_m$, respectively. Due to the symmetry, only half of the domain $(L_0/2\le x_2\le L_0)$ is simulated. The diffuse boundary condition at $x_2=L_0$ yields

(3.5)\begin{equation} v_2\le 0: \quad {f_0^{+}}=\frac{{n^{-}}(x_2=L_0)}{n_0}{E_t(T_0,\boldsymbol{u}_2)}, \quad {f_1^{+}}=\frac{d_r}{2}k_BT_0{f_0^{+}}, \quad {f_2^{+}}=\frac{d_v}{2}k_BT_0{f_0^{+}}, \end{equation}

where ${n^{-}}(x_2=L_0)$ is determined as the same way as (3.3), and $E_t(T_0,\boldsymbol {u}_2)$ is the equilibrium distribution function defined in (2.21) with temperature $T_0$ and macroscopic velocity $\boldsymbol {u}_2$. The symmetrical condition at $x_2=L_0/2$ reads

(3.6)\begin{equation} v_2\ge 0: \quad {f_0^{+}}={f_0^{-}}({-}v_1,-v_2,v_3), \quad {f_1^{+}}=\frac{d_r}{2}k_BT{f_0^{+}}, \quad {f_2^{+}}=\frac{d_v}{2}k_BT{f_0^{+}}. \end{equation}

The results from our kinetic models I and II as well as the DSMC simulation at ${Kn}_{{gas}}=0.5$ are shown in figure 4. The kinetic model II with BCO shows better accuracy, while the relative error given by kinetic model I in translational heat flux is around $7\,\%$. It can be seen that the vibrational temperature is much lower than the rotational one in the Couette flow, since the energy increase in internal DoF only comes from the exchange with translational DoF in this problem. Thus, a larger collision number leads to less increase in internal temperature at the same distance from the wall (due to the infrequent relaxation with the translational mode), and also contributes less to the heat flux.

Figure 4. Comparisons of the (a) density, (b) flow velocity, (c) temperature and (d) heat flux $q_1$ in the flow direction and $q_2$ perpendicular to the flow direction of nitrogen, between kinetic model I, kinetic model II and DSMC for the one-dimensional Couette flow at ${Kn}_{{gas}}=0.5$.

3.4. Creep flow driven by the Maxwell demon

The creep flow driven by the Maxwell demon is a thought test (Li et al. Reference Li, Zeng, Su and Wu2021), where each gas molecule is subjected to an external acceleration based on its kinetic energy,

(3.7)\begin{equation} a_1 = a_0\left(\frac{v_1^2}{v_m^2}-\frac{3}{2}\right). \end{equation}

That is, fast molecules are forced toward the positive direction, while slow molecules move in the opposite direction.

Consider the nitrogen flow driven by the Maxwell demon confined between two parallel plates with distance $L_0$ apart. To solve the force-driven flow, we choose small values of $a_0$ so that the gas flow deviates only slightly from the global equilibrium; the acceleration acting on the molecules is linearized, which results in

(3.8) \begin{equation} \left. \begin{array}{c} \displaystyle \dfrac{\partial{(\boldsymbol{a}f_0)}}{\partial{\boldsymbol{v}}} =\dfrac{2a_0L_0}{v_m^2}v_1E_t(T_0)\left(\dfrac{v_1^2}{v_m^2}- \dfrac{5}{2}\right),\\ \displaystyle \dfrac{\partial{(\boldsymbol{a}f_1)}}{\partial{\boldsymbol{v}}}= \dfrac{d_ra_0L_0}{v_m^2}v_1k_BT_0E_t(T_0)\left(\dfrac{v_1^2}{v_m^2} -\dfrac{5}{2}\right), \\ \displaystyle \dfrac{\partial{(\boldsymbol{a}f_2)}}{\partial{\boldsymbol{v}}}= \dfrac{d_va_0L_0}{v_m^2}v_1k_BT_0E_t(T_0)\left(\dfrac{v_1^2}{v_m^2} -\dfrac{5}{2}\right). \end{array}\right\} \end{equation}

The plates at rest are fully diffuse, then the boundary conditions are simply given by (3.2) and (3.3), but with the wall temperature replaced by $T_0$.

Figure 5 shows the good agreement between the solutions of kinetic model II and DSMC at ${Kn}_{{gas}}=1$. The rotational/vibrational heat flux is one/two orders of magnitude smaller than the translational heat flux and, hence, makes a negligible contribution to the total heat transfer in this problem. However, significant discrepancies are observed in both velocity and heat flux profiles predicted by kinetic model I. It demonstrates the importance of using correct velocity-dependent collision rates based on the intermolecular potential, as it is realized in kinetic model II.

Figure 5. Comparisons of the (a) velocity and (b) heat flux in the flow direction of nitrogen between kinetic model I (green lines), kinetic model II (red lines) and DSMC (blue circles) for one-dimensional creep flow driven by the Maxwell demon at ${Kn}_{{gas}}=1$. Both the flow velocity and the heat flux have been further normalized by ${2a_0L_0}/{v_m^2}$.

To assess the influence of the thermal relaxation rates on the creep flow, two more cases are conducted using kinetic model II by varying the values of the matrix $A$ but keeping the Eucken factors fixed. More specifically, the off-diagonal elements in $A$ in the two cases are set to be zero and double of those given by DSMC, respectively. The values of diagonal elements are calculated based on (2.48) using the fixed Eucken factors. Figure 6 shows that these relaxation rates affect the flow velocity and heat fluxes, despite the fact that the thermal conductivities are fixed. In particular, when the off-diagonal elements in $A$ are zero, the heat fluxes of different types of DoF are decoupled, so that the internal heat fluxes are exactly zero. These situations occur in many traditional kinetic models, such as the Rykov model and the ES-BGK model. This example demonstrates the importance of recovering the fundamental thermal relaxation process rather than the apparent thermal conductivities in rarefied gas flow simulations.

Figure 6. Same as figure 5, except that the off-diagonal elements in $A$ are set to be zero (blue), the values from DSMC (red) and double of those from DSMC (green), respectively.

3.5. Normal shock wave

In the simulations of the normal shock wave of nitrogen, the upstream number density $n_u=n_0$ and temperature $T_u=T_0$ are chosen to be the reference values, which also determine the characteristic length to be $L_0={16\mu (T_0)}/({5n_0\sqrt {2{\rm \pi} mk_BT_0}})$ and, hence, ${Kn}_{{gas}}={5{\rm \pi} }/{16}$ in this problem. The total length of the simulation domain is $90L_0$, so that the boundary conditions at both ends can be approximated by equilibrium states (the wave front is initially located at $x=0$),

(3.9) \begin{equation} \left. \begin{array}{c@{}} \displaystyle x={-}30L_0, v\ge 0: \quad f_0=\dfrac{n_{u}}{n_0}E_t(T_u), \quad f_1=\dfrac{d_r}{2}k_BT_uf_0, \quad f_2=\dfrac{d_v}{2}k_BT_uf_0, \\ \displaystyle x=60L_0, v\le 0: \quad f_0=\dfrac{n_{d}}{n_0}E_t(T_d), \quad f_1=\dfrac{d_r}{2}k_BT_df_0, \quad f_2=\dfrac{d_v}{2}k_BT_df_0, \end{array}\right\} \end{equation}

where the subscripts $u$ and $d$ represent the upstream and downstream ends, respectively. Given the Mach number, the macroscopic quantities at the downstream end are determined by the Rankine–Hugoniot relation.

3.5.1. Comparison with DSMC

Numerical results of both kinetic models I, II and DSMC are compared in figure 7, when the Mach number is ${Ma}=5$. As expected, kinetic model II reproduces the structure of normal shock waves with high accuracy. While kinetic model I significantly overestimates the temperature, heat flux and deviated pressure before the wave front, because the velocity dependence of the collision frequency is not involved in the RTA. The rotational and vibrational collision numbers, $Z_r$ and $Z_v$, which affect the energy exchange rate between internal and translational modes, play roles in the difference between rotational and vibrational temperatures. That is, the distance for the vibrational temperature to reach equilibrium is much longer than that for rotational modes. This is consistent with the fact that we set $Z_v=10Z_r$.

Figure 7. Comparisons of the (a) density and velocity, (b) temperature, (c) deviated pressure and (d) heat flux of nitrogen between kinetic model I (green lines), kinetic model II (red lines) and DSMC (blue circles) for a normal shock wave at ${Ma}=5$.

3.5.2. Comparison with experiments

Then we compare the results of our kinetic models with experimental data for normal shock waves with high enthalpy incoming nitrogen flow reported by Alsmeyer (Reference Alsmeyer1976). The experiments were performed in a conventional shock tube for argon and nitrogen with the Mach number range from 1.55 to 9 and from 1.5 to 10, respectively, and the density distribution of the shock wave was measured by the attenuation of an electron beam. The differentiated density profiles and reciprocal shock thickness may have an error of 4 %, as it is mentioned in Alsmeyer (Reference Alsmeyer1976).

According to the experimental conditions, the temperature of the upstream flow is 300 K. Thus, the temperature at the downstream end exceeds 2000 K when ${Ma}>6$, where the vibrational DoF of nitrogen cannot be ignored. Also, the temperature dependence of the vibrational DoF (2.3) and internal collision numbers in the modelling have to be considered to compare with experimental data. For diatomic molecules, the most commonly used empirical model for the temperature-dependent vibrational collision number $Z_v$ is given by (Millikan & White Reference Millikan and White1963)

(3.10)\begin{equation} Z_v=\frac{3{\rm \pi}}{4(3+d_v)}\left(\frac{C_1}{T_v^{\omega}}\right) \exp{(C_2T_v^{{-}1/3})}, \end{equation}

with parameters $C_1=9.1, C_2=220$. As for the rotational collision number, discrepancies for $Z_r$ between different experimental data are large (Carnevale & Larson Reference Carnevale and Larson1967; Annis & Malinauskas Reference Annis and Malinauskas1971). Parker (Reference Parker1959) presents an analytical derivation of $Z_r$, which shows good predictions of selected experimental data with appropriately chosen parameters. However, it gives a rather strong increase of $Z_r$ with temperature (Nyeland Reference Nyeland1967), and does not distinguish the difference between the translational and rotational temperatures (Lordi & Mates Reference Lordi and Mates1970). Nevertheless, MD simulations provide a reliable way to extract collision numbers for specific gas systems. Valentini et al. (Reference Valentini, Zhang and Schwartzentruber2012) studied $Z_r$ of nitrogen with the Lennard-Jones intermolecular potential using MD simulations, which shows weaker dependence on the gas temperature but a strong influence of difference between the translational and rotational temperatures. Also, a demonstrative model of $Z_r$ is provided by Valentini et al. (Reference Valentini, Zhang and Schwartzentruber2012), which is fitted by their MD results at the temperature below 2000 K,

(3.11)\begin{equation} \left. \begin{array}{c@{}} \displaystyle Z_r=\dfrac{\rm \pi}{4}Z_r^0\left[1-b_1 \left(1-\dfrac{T_r}{T_t}\right)\right], \\ Z_r^0=a_1T_t^{1/4}+a_2T_t^{{-}1/4}+a_3(T_t-1000\ \text{K}), \end{array}\right\} \end{equation}

where the parameters are $a_1=1.33868, a_2=-6.19992, a_3=-0.00107942,b_1=1$. When $T_t$ goes to 2000 K, $Z_r^0$ approaches a constant value, which will be used for higher $T_t$.

Because of the variations of the internal collision numbers from the upstream to the downstream end, the resultant changes in relaxation rates of heat flux have to be considered. We adopt the thermal relaxation rates $\boldsymbol {A}$ derived from the asymptotic expansion of the WCU equation (Mason & Monchick Reference Mason and Monchick1962; Gorji & Jenny Reference Gorji and Jenny2013),

(3.12)\begin{align} \boldsymbol{A}=\left[ \begin{array}{ccc} \displaystyle{\dfrac{2}{3}+\dfrac{5}{6}\left(\dfrac{d_r}{(3+d_r) Z_r}+\dfrac{d_v}{(3+d_v)Z_v}\right)} & \displaystyle{-\dfrac{5}{2(3+d_r)Z_r}} & \displaystyle{-\dfrac{5}{2(3+d_v)Z_v}} \\ \displaystyle {-\dfrac{d_r}{2(3+d_r)Z_r}} & \displaystyle{{Sc}+\dfrac{3}{2(3+d_r)Z_r}} & 0 \\ \displaystyle{-\dfrac{d_v}{2(3+d_v)Z_v}} & 0 & \displaystyle{{Sc}+\dfrac{3}{2(3+d_v)Z_v}} \end{array} \right]. \end{align}

The dimensionless reciprocal shock thickness, denoted as $\delta ^{-1}$, is defined as the maximum gradient of the density (Alsmeyer Reference Alsmeyer1976),

(3.13)\begin{equation} \delta^{{-}1} = \max\left(\frac{\partial n^*}{\partial x}\right)L_0, \end{equation}

where the density is normalized by $n^*=(n-n_d)/(n_u-n_d)$.

Figure 8 compares the reciprocal shock thickness and density profiles obtained by kinetic models and experiments, where excellent agreements can be observed. Although kinetic model I significantly overestimates the temperature before the wave front, as discussed previously, the density distribution and, hence, the shock thickness are not affected that much.

Figure 8. Comparisons of (a) the reciprocal shock thickness $\delta ^{-1}$ and (b,c) density profiles between kinetic model I (green squares/lines), kinetic model II (red diamonds/line) and experimental measurements reported by Alsmeyer (Reference Alsmeyer1976) (blues circles).

4.1. Fourier flow

Consider the heat transfer between two plates maintained at different temperatures $T_1=0.5T_0$ and $T_2=1.5T_0$, which emit photons with Planck distribution and reflect gas molecules diffusely. The Knudsen numbers of gas is ${Kn}_{{gas}}=0.1$, ${Kn}_{{photon}}$ varies from 0.1 to 100, $\tilde {\sigma }_R$ changes from 0.001 to 0.1 and $T_0/T_{{ref}}=2$.

It can be seen from figure 9 that the radiation significantly changes the profiles of vibrational temperature (figure 9b), especially in the low-temperature region, while its influence on translational temperature is relatively small (figure 9a). Also, the presence of radiation leads to a remarkable variation of vibrational heat flux (figure 9c), which indicates a strong coupling between molecular vibrational modes and photons. The radiative heat flux and its proportion to the total one are calculated and shown in figure 10(a) and 10(b), respectively. Firstly, it is observed that the radiative heat flux $q_R$ and its proportion to the total one increase with $\tilde {\sigma }_R$. This is because $\tilde {\sigma }_R$ indicates the relative strength of the radiative intensity in terms of the characteristic heat flux of gas flow. Secondly, when ${Kn}_{{photon}}$ approaches 100, the radiative heat flux is found to be linearly proportional to $\tilde {\sigma }_R$. The reason is that no interaction between gas molecules and photons occurs when ${Kn}_{{photon}}$ is large enough and, thus, the gas flow is transparent to radiation. Thirdly, when ${Kn}_{{photon}}=0.1$, the energy transported by the radiation field becomes negligible. This is due to the extremely short transport distance of radiative energy at small ${Kn}_{{photon}}$, which implies that the photons emitted by vibrational transitions are absorbed immediately. In summary, in the Fourier flow the radiative heat transport is important and even dominated when both $\tilde {\sigma }_R$ and ${Kn}_{{photon}}$ are large. The relative difference between the radiative heat flux from kinetic models I and II are less than $0.1\,\%$ in these cases.

Figure 9. Comparisons of the (a) translational temperature, (b) vibrational and radiative temperature, (c) vibrational and radiative heat flux of nitrogen in the Fourier flow between different radiation strengths, when ${Kn}_{{gas}}=0.1$ and $T_0/T_{{ref}}=2$. The results are obtained from kinetic model II.

Figure 10. The radiative heat flux change with ${Kn}_{{photon}}$ and $\tilde {\sigma }_R$ in the Fourier flow, when ${Kn}_{{gas}}=0.1$ and $T_0/T_{{ref}}=2$, (a) radiative heat flux, (b) ratio of radiative heat flux to the total one.

4.2. Couette flow

For the Couette flow between two plates with the same temperature $T_0$, the velocities of the lower and upper plates are $u_1=-v_m$ and $u_2=v_m$, respectively. Two vibrational collision numbers are considered $Z_v=2Z_r$ and $Z_v=10Z_r$ here. The Knudsen number of gas is ${Kn}_{{gas}}=0.1$, ${Kn}_{{photon}}$ varies from 0.1 to 100, $\tilde {\sigma }_R$ changes from 0.01 to 0.1 and $T_0/T_{{ref}}=2$.

Due to the symmetry, only half of the domain $(L_0/2\le x_2\le L_0)$ is shown in figure 11. A large vibrational collision number ($Z_v=10Z_r$) suppresses the rise of the vibrational temperature, and, hence, limits the radiation effect. When $Z_v=2Z_r$, the radiative effect lowers the vibrational temperature and heat flux significantly compared with those of the non-radiative gas flow. However, the radiation itself may contribute a non-negligible part to the total heat flux.

Figure 11. Comparisons of the (a,c) temperature, (b,d) heat flux of nitrogen Couette flow, when ${Kn}_{{photon}}=1$, $\tilde {\sigma }_R=0.1$, ${Kn}_{{gas}}=0.1$ and $T_0/T_{{ref}}=2$; (a,b) $Z_v=10Z_r$, (c,d) $Z_v=2Z_r$.

Figure 12 shows the variations of the ratio of radiative heat flux to the total one with respect to ${Kn}_{{photon}}$ and $\tilde {\sigma }_R$. Contrary to the results in the Fourier flow problems (figure 10), the proportion of radiative heat flux in the Couette flow is not a monotonical function of ${Kn}_{{photon}}$. With the increase of ${Kn}_{{photon}}$, the radiative heat flux firstly becomes significant due to the longer distance the photons can propagate before being absorbed by gas molecules. However, it decreases gradually when ${Kn}_{{photon}}$ further increases, since the rise of radiative temperature is generated by the gas–photon interaction in the Couette flow, and the lower emissivity of photons (larger ${Kn}_{{photon}}$) reduces the radiative temperature difference between the gas and the wall. As the result of the two opposite mechanisms, the radiative heat flux reaches the maximum value at an intermediate value of ${Kn}_{{photon}}$, which is around 1 in these cases.

Figure 12. The ratio of radiative heat flux to the total heat flux change with ${Kn}_{{photon}}$ and $\tilde {\sigma }_R$ in the Couette flow, when ${Kn}_{{gas}}=0.1$ and $T_0/T_{{ref}}=2$; (a) $Z_v=2Z_r$, (b) $Z_v=10Z_r$.

4.3. Normal shock wave

When the temperature-dependent vibrational DoF is considered, the specific heat ratio changes with the temperature across the shock wave structure, thus, the Rankine–Hugoniot relation for the upstream and downstream macroscopic quantities is no longer applicable. In this situation, consider the conservation of mass, momentum and energy of the normal shock wave,

(4.1)\begin{equation} \left. \begin{array}{c@{}} n_1u_1=n_2u_2, \\ mn_1u_1^2+p_1=mn_2u_2^2+p_2,\\ \displaystyle c_{p,1}T_1+\dfrac{1}{2}u_1^2+2{\rm \pi} \int{I^R_1\cos{\theta}\,\mathrm{d}\theta}=c_{p,2}T_2+ \dfrac{1}{2}u_2^2+2{\rm \pi}\int{I^R_2\cos{\theta}\,\mathrm{d}\theta}, \end{array}\right\} \end{equation}

where $c_p$ is the specific heat capacity with constant pressure, which relates to the specific heat ratio $\gamma (T)=(5+d_r+d_v(T))/(3+d_r+d_v(T))$ as $c_p=[\gamma /(\gamma -1)]k_B/m$. The heat flux due to radiation vanishes far away from the shock wave layer, then the (4.1) are solved,

(4.2a–c)\begin{equation} \frac{n_2}{n_1}=1-\frac{1}{\gamma(T_1){Ma}^2}\left(\frac{p_2}{p_1}-1\right), \quad \frac{u_2}{u_1}=\frac{n_1}{n_2}, \quad \frac{T_2}{T_1}=\frac{p_2}{p_1}\frac{n_1}{n_2}, \end{equation}

with

(4.3)\begin{align} \frac{p_2}{p_1}&=\frac{((\gamma(T_1)^2{Ma}^4+ \gamma(T_2)^2)(\gamma(T_1)-1)-2{Ma}^2\gamma(T_1)(\gamma(T_2)^2- \gamma(T_1)))^{1/2}}{(\gamma(T_1)-1)^{1/2}(\gamma(T_2)+1)} \nonumber\\ &\quad +\frac{\gamma(T_1){Ma}^2+1}{\gamma(T_2)+1}, \end{align}

where the Mach number is defined by the velocity of upstream flow ${Ma}=u_1/\sqrt {\gamma (T_1)k_BT_1/m}$.

Figure 13 shows the density, velocity, temperature and heat flux across the shock wave structure with different ${Kn}_{{photon}}$. The smaller ${Kn}_{{photon}}$ means a more frequent interaction between vibrational mode and radiation field and, thus, makes the radiative temperature very close to the vibrational one. When ${Kn}_{{photon}}=1$ (figure 13a–c), only the vibrational temperature and heat flux deviate from the non-radiative one slightly, while the other gas flow properties are almost unaffected. When ${Kn}_{{photon}}$ increases to 10 (figure 13d–f), the radiative temperature is more smoothly distributed between upstream and downstream ends (figure 13e), hence, significantly changing the gas flow properties and leading to a thicker normal shock wave.

Figure 13. Comparisons of the (a,d) density and velocity, (b,e) temperature and (c,f) heat flux distribution across the shock wave structure, when ${Ma}=2$, $\tilde {\sigma }_R=0.1$ and $T_0/T_{{ref}}=1$. Results are shown for (a–c) ${Kn}_{{photon}}=1$, (d–f) ${Kn}_{{photon}}=10$.

5.1. Effect of radiation

When ${Kn}_{{photon,ref}}=100$ and $\tilde {\sigma }_R=10$, the temperature and heat flux of each mode of the surrounding gas are shown in figure 14 for both non-radiative (dashed lines) and radiative (solid lines) gas flow. Due to the direct energy exchange between the vibrational mode of gas molecules and photons, the vibrational temperature is significantly reduced by radiation, while the influence on both the translational and rotational modes are relatively small. Therefore, the difference between the translational and vibrational temperatures enlarges in the radiative hypersonic gas flow, since both $\mu _b$ and $\mu _b^R$ contribute to the thermal non-equilibrium. To be specific, the bulk viscosity $\mu _b$ resisting the volume change makes the maximum vibrational temperature about 1/4 of the translational one. Meanwhile, the bulk viscosity $\mu _b^R$ resisting the radiation transition further lowers the maximum vibrational temperature to be 1/12 of the translational one. Therefore, the heat flux from the vibrational mode becomes negligible in this case.

Figure 14. A ${Ma}=15$ shock wave passing through the cylinder. The distribution of the (a) translational, (c) rotational, (e) vibrational temperatures (normalized by $T_0$), and the (b) translational, (d) rotational (f), vibrational heat fluxes (normalized by $n_0k_BT_0v_m$) around the cylinder solved by kinetic model equations when the incoming flow has ${Kn}_{{gas}}=0.05$ and $T_0=T_{{ref}}/2$. The dashed lines represent the results of the case without radiation, while the solid lines are the results of the case with ${Kn}_{{photon,ref}}=100$ and $\tilde {\sigma }_R=10$. Note that the colourbars for vibrational temperature and heat flux are not in the linear scale.

The photon absorptivity is proportional to the population of the vibrational energy state ((2.14a,b)), thus, it increases with the gas density and effective vibrational DoF. Figure 15(a) shows that, the order of magnitudes of the photon absorptivity increases significantly in the stagnation region and, hence, the local photon Knudsen number there becomes much lower. The vibrational/radiative temperatures and heat fluxes along the stagnation line are shown in figures 15(b) and 15(c). It can be seen that the variation of radiative temperature along the stagnation line is relatively small when compared with that of the vibrational temperature, since the photon Knudsen number is much larger than that of the gas over the entire domain. It is also found that the vibrational temperature quickly approaches the radiative temperatures on the surface of the cylinder. It is attributed to the high photon absorptivity (strong gas–photon interaction equivalently) near the stagnation point, which is caused by the high molecular density there. However, owing to the relatively high radiative strength $\tilde {\sigma }_R$, the radiative heat flux is one order of magnitude larger than that of the vibrational mode, and contributes around 32 % to the total heat flux over the cylinder surface.

Figure 15. The distribution of (a) the local photon Knudsen number and absorptivity (normalized by $1/L_0$), vibrational and radiative (b) temperature (normalized by $T_0$), and (c) heat flux (normalized by $n_0k_BT_0v_m$) along the stagnation line before the cylinder, when the incoming flow has ${Kn}_{{gas}}=0.05$, ${Ma}=15$, $T_0=T_{{ref}}/2$, ${Kn}_{{photon,ref}}=100$ and $\tilde {\sigma }_R=10$. The origin $x/L_0=0$ is located at the centre of the cylinder.

5.2. Influence of ${Kn}_{{photon}}$ and $\tilde {\sigma }_R$

By varying the reference photon Knudsen numbers ${Kn}_{{photon,ref}}$ from $10^{-1}$ to $10^3$, and the dimensionless radiative strength $\tilde {\sigma }_R$ from $10^{-1.5}$ to $10^1$, their influence on the radiative heat flux on the surface of the cylinder are systematically revealed. Figure 16 shows the magnitude of normal heat flux from both convection and radiation to the cylinder. Note that ${Kn}_{{photon,ref}}$ indicates the strength of gas–photon interaction at the reference state, which is determined by the molecular structure and number density.

Figure 16. The heat flux (normalized by $n_0k_BT_0v_m$) from convection (blue lines) and radiation (red lines) along the cylinder surface: (a) ${Kn}_{{photon,ref}}=10$, $\tilde {\sigma }_R=$1 (solid lines), 10 (dashed lines). (b) $\tilde {\sigma }_R=1$, ${Kn}_{{photon,ref}}=$0.1 (solid lines), 10 (dashed lines), 100 (dashed dot lines). Note that $\theta$ is the clockwise angle measured from the stagnation streamline. Also note that the two lines of convective heat flux when ${Kn}_{{photon,ref}}=0.1$ and 100 are overlapped in (b).

When ${Kn}_{{photon,ref}}$ is fixed (figure 16a), a higher value of $\tilde {\sigma }_R$ not only increases the relative radiative intensity and heat flux, but also leads to stronger coupling between the gas flow and the radiation field, and hence, decreases the convection heat flux. Therefore, the relative importance of radiative heat transfer becomes significant.

When $\tilde {\sigma }_R$ is fixed (figure 16b), the influence of ${Kn}_{{photon,ref}}$ in the windward region is different at different regimes of photon Knudsen number. Two mechanisms are summarized below.

1. (i) The gas–photon interaction becomes weaker with the increase of ${Kn}_{{photon,ref}}$ and, thus, makes the radiative temperature much lower than the gas temperature. Consequently, the radiative heat flux decreases while the convective heat flux increases. This mechanism dominates the influence of ${Kn}_{{photon,ref}}$ when the photon Knudsen number is relatively large, e.g. when ${Kn}_{{photon,ref}}$ changes from 10 to 100 shown in figure 16(b).

2. (ii) The mean free path of photons (reciprocal of absorptivity) increases with ${Kn}_{{photon,ref}}$, and then the transportation distance of energy carried by the photon before it is absorbed by the gas becomes longer. Therefore, the radiative heat flux increases. This mechanism surpasses the first one when the photon Knudsen number is relatively small, where the radiative temperature is almost the same as the vibrational temperature, e.g. when ${Kn}_{{photon,ref}}$ changes from 0.1 to 10 as shown in figure 16(b).

On the other hand, the radiative heat flux in the leeward side of the cylinder decreases monotonically when ${Kn}_{{photon,ref}}$ changes from 0.1 to 100 (figure 16b). The reason is that the local photon Knudsen number is much higher in the leeward region (due to the low gas density there) than that in the windward region, thus, the first mechanism always dominates. Therefore, due to the significant difference in the local photon Knudsen number on the two sides of the cylinder, the radiative heat flux becomes much more important than the convective one in the leeward region when ${Kn}_{{photon,ref}}$ is small. For example, the radiative heat flux is almost uniformly distributed along the surface of the cylinder when ${Kn}_{{photon,ref}}=0.1$.

To quantitatively measure the importance of the radiation in hypersonic gas flow as ${Kn}_{{photon,ref}}$ and $\tilde {\sigma }_R$ change, the ratio of radiative heat flux to the total one on the surface of the cylinder is calculated and shown in figure 17. The contour lines show two regimes of the influence from the two parameters: when ${Kn}_{{photon,ref}}\gg 1$, the proportion of the radiative heat flux is constant along $\tilde {\sigma }_R/{Kn}_{{photon,ref}}$; when ${Kn}_{{photon,ref}}\ll 1$, it is constant along $\tilde {\sigma }_R{Kn}_{{photon,ref}}$. Therefore, it gives clear criteria to determine the importance of radiative heat transfer in these problems.

Figure 17. The ratio of radiative heat flux to the total heat flux on the surface of the cylinder in hypersonic gas flow with the variation of ${Kn}_{{photon,ref}}$ and $\tilde {\sigma }_R$.