Kinetic modelling of rarefied gas flows with radiation

Abstract

Two kinetic models are proposed for high-temperature rarefied (or non-equilibrium) gas flows with internal degrees of freedom and radiation. One of the models uses the Boltzmann collision operator to model the translational motion of gas molecules, which has the ability to capture the influence of intermolecular potentials, while the other adopts the relaxation time approximations, which has higher computational efficiency. In our kinetic model equations, not only the transport coefficients such as the shear/bulk viscosity and thermal conductivity but also their underlying relaxation processes are recovered. The non-equilibrium dynamics of gas flow and radiation are tightly coupled, where the transport properties of gas molecules and photons are correlatively dependent. The proposed kinetic models are validated by the direct simulation Monte Carlo method in several non-radiative rarefied gas flows (e.g. the normal shock wave, Fourier flow, Couette flow and the creep flow driven by the Maxwell demon), and the experimental data of planar heat transfer and normal shock waves in nitrogen. Then, the rarefied gas flows with strong radiation are studied based on the kinetic models, not only in the above one-dimensional gas flows, but also in the two-dimensional radiative hypersonic flow passing a cylinder. The characteristics of heat transfer in the tightly coupled fields of gas and radiation are systematically investigated, particularly the influence of the non-equilibrium photon transport and their interactions with gas molecules are revealed. It is found that the radiation makes a profound contribution to the total heat transfer in radiative hypersonic flow at an intermediate photon Knudsen number.

1. Introduction

The non-equilibrium dynamics of molecular (diatomic/polyatomic) gas is commonly encountered in aerospace engineering. For example, at a high Mach number, the air surrounding an aircraft decelerates and heats up rapidly after compression by shock waves, which causes a strong conversion from the translational energy into the internal energy (i.e. rotational, vibrational and electronic energies). The temperature may reach thousands of degrees Kelvin and leads to significant changes in the physical and chemical properties of the gas (Ivano & Gimelshein 1998; Anderson 2019). Meanwhile, the thermal radiation, which is induced by transitions between excited states of gas molecules, makes a significant contribution to the overall heat load on hypersonic aircraft. For instance, in the re-entry into Mars’ atmosphere, the shock layer radiation can constitute a larger portion than convective heating in the stagnation region of spacecraft (Edquist et al. 2014), and the order of magnitude analysis suggests that the absolute radiative power at the stagnation point scales as the eighth power of entry velocity (da Silva & Beck 2011). Therefore, accurate predictions of molecular gas flow and radiative energy transfer are essential for developing thermal protection systems for aircraft entering planetary atmospheres. Additionally, radiative gas dynamics also play an important role in the interpretation of spectrometer measurements (Horvath et al. 2010) and technologies of gas dynamics lasers.

Under the assumption of thermodynamic equilibrium, the traditional Navier–Stokes–Fourier equations with proper boundary conditions, involving velocity slip and temperature jump, are used to predict the thermal environment and aerodynamic characteristics of the aircraft in the near-continuum regime. The influence of internal degrees of freedom (DoF) is taken into account by the variations of heat capacity and transport properties of molecular gas (Malik & Anderson 1991). On the other hand, when thermodynamic non-equilibrium occurs, gases with different temperatures associated with various relaxation processes need to be considered. Several sets of Navier–Stokes-type equations have been developed with multi-temperatures of different types of kinetic modes (Taylor & Bitterman 1969; Lee 1984; Park 1985; Colonna et al. 2006; Bruno & Giovangigli 2011; Aoki et al. 2020), and applied in the simulations of thermal non-equilibrium gas dynamics (Kustova et al. 2011; Armenise, Reynier & Kustova 2016). In addition, studies of radiative gas flow are largely oriented towards hypersonic re-entry vehicles (Vincent & Traugott 1971). Since the transport properties of gas molecules and photons depend on each other, the fields of gas flow and radiation have to be determined self-consistently. Conventionally, the coupling of gas dynamics and radiative heat transfer is achieved by integrating the radiative heat transfer as a heat source into the total energy conservation equation.

Because the macroscopic models are obtained at small Knudsen numbers from the gas kinetic theory, they are only applicable in the near-continuum flow regime, where the mean free path of gas molecules is much smaller than the characteristic flow length. However, the gas could be in highly thermal non-equilibrium in many realistic situations, such as the re-entry of aircraft into the atmosphere, where the gas flow changes from the continuum to the free molecular regime. Therefore, the treatment based on gas kinetic theory from a mesoscopic perspective is inevitable, as the molecular dynamics (MD) simulation at the microscopic level is limited to small spatial and temporal domains. The fundamental equation in gas kinetic theory is the Boltzmann equation, but it is only rigorously established for monatomic gas. For the molecular gas, its internal DoF pose difficulties in the kinetic modelling. The heuristic way to describe the molecular gas dynamics in all flow regimes is using the Wang-Chang–Uhlenbeck (WCU, Wang-Chang & Uhlenbeck 1951) equation, which treats the internal DoF quantum mechanically and assigns each internal energy level an individual velocity distribution function. Besides, a general framework for the kinetic modelling of molecular gases is proposed recently based on a set of allowed internal states endowed with a suitable measure (Borsoni, Bisi & Groppi 2022). However, the complexity and excessive computational burden prevent their engineering applications.

The direct simulation Monte Carlo (DSMC) method (Bird 1994) is prevailing in simulating the rarefied gas dynamics. Although it is proven that DSMC is equivalent to the Boltzmann equation for monatomic gas (Wagner 1992), there are some drawbacks when it is applied to radiative molecular gas flows. First, the bulk viscosity and the thermal conductivities cannot be recovered simultaneously. The reason lies in its phenomenological collision model of Borgnakke & Larsen (1975), which realizes the correct exchange rate between the translational and internal energies to exactly recover the bulk viscosity (Haas et al. 1994; Gimelshein, Gimelshein & Lavin 2002); however, it cannot guarantee that the thermal conductivity, or its translational and internal components, is recovered at the same time (Wu et al. 2020; Li et al. 2021). Second, the DSMC method has been coupled with the photon Monte Carlo method to simulate radiative gas flow, by calculating the rates of radiative heating/cooling at each time step (Sohn et al. 2012). Because of the high sensitivity of radiation rates on temperature, however, the fluctuation of temperature sampled in simulation cells may lead to significant instabilities (Prem et al. 2019). Third, DSMC is computationally costly in the simulation of low-speed or low-Knudsen-number flows due to its intrinsic stochastic nature. For instance, it has been found that the computational cost increases as ${Ma}^{-2}$ (${Ma}$ is the Mach number) when the flow speed is approaching zero (Hadjiconstantinou et al. 2003). Therefore, the multiscale feature in non-equilibrium hypersonic flows passing spacecraft makes the DSMC method time consuming and even intractable in some cases.

Alternatively, kinetic models are proposed to imitate as closely as possible the behaviour of the WCU equation, and multiscale deterministic methods are developed to solve those kinetic models. The Bhatnagar–Gross–Krook (BGK) type kinetic models (Bhatnagar, Gross & Krook 1954; Holway 1966; Shakhov 1968a), which replace the Boltzmann collision operator (BCO) with a single relaxation approximation, has achieved notable success in the modelling of monatomic rarefied gas flows. These kinetic models have been extended to model molecular gas by introducing additional internal energy variables in the distribution function (Morse 1964; Rykov 1975; Andries et al. 2000; Rahimi & Struchtrup 2016; Tantos et al. 2016; Wang et al. 2017; Bernard, Iollo & Puppo 2019), as well as the gas mixture of polyatomic molecules (Pirner 2018). Besides, the Fokker–Planck models have been proposed (Gorji & Jenny 2013; Mathiaud & Mieussens 2020), which take advantage of the continuous distribution functions in terms of stochastic velocity processes to speed up the stochastic particle methods.

These models however do not reduce to the Boltzmann equation for monatomic gases when the translational–internal energy exchange is absent. Therefore, these models cannot distinguish the influence of different intermolecular potentials. For example, the uncertainties caused by different intermolecular potentials have been demonstrated in the calculation of the thermal creep slip on diffuse walls (Loyalka 1990), the thermal creep and Poiseuille flow (Sharipov & Bertoldo 2009; Takata & Funagane 2011; Wu et al. 2015a), the viscous slip of the Couette flow (Su et al. 2019), and the slip coefficient and mass flow rate of thermal transpiration (Wang, Su & Wu 2020). On the other hand, all these kinetic model equations concern only the transport coefficients, such as the thermal conductivity and bulk viscosity, while their fundamental relaxation processes are not captured, which are found to be important in rarefied molecular gas dynamics. For example, the relaxation rates of heat flux can significantly affect the creep flow driven by a molecular velocity-dependent external force (Li et al. 2021). Thus, it is necessary to tackle the two difficulties when building a gas kinetic model for molecules with rational and vibrational DoF.

Additionally, the kinetic model of molecular gas flow with radiative heat transfer is still far from being well developed. Groppi & Spiga (1999) proposed a kinetic model with radiation transitions between energy levels by the absorption and emission of photons, as an extension of the WCU equation with a radiative field. However, it is unlikely to be solved practically due to its even much higher complexity than the WCU equation. Therefore, the tractable kinetic model equation incorporating radiative heat transfer is highly desirable and urgently needed.

Hence, the present work is dedicated to developing general kinetic models of molecular gas in radiative rarefied flow, which include translational, rotational and discrete vibrational modes with correct rates of relaxation processes. The rest of the paper is organized as follows. In § 2 the kinetic models with both relaxation time approximation (RTA) and BCO are proposed, and the transport coefficients and their intrinsic relation to relaxation rates are discussed. In § 3 the kinetic models without a radiation field are validated by DSMC in typical rarefied gas flows and experimental data of planar heat transfer and normal shock waves. Then, in § 4 the influence of radiative heat transfer is examined by solving the kinetic models. Furthermore, the kinetic models are applied to solve hypersonic flow passing a cylinder where the radiative heat transfer becomes essential in § 5. Finally, conclusions are presented in § 6.

2. Kinetic model

In gas kinetic theory the distribution function is used to describe the status of the dilute gaseous system at the mesoscopic level. The evolution of the distribution function of molecular gas is governed by the WCU equation, which is too complicated to be applied to realistic problems. Therefore, kinetic models are urgently needed to simplify the collision operator. A fundamental requirement in constructing a kinetic model is that all the transport coefficients are consistent with those obtained from the Boltzmann equation for monatomic gas or the WCU equation for molecular gas. In dilute gases, due to the excitation of rotational and vibrational DoF in molecular gas, additional relaxation processes occur between different types of energies, which lead to exclusively transport coefficients in molecular gas such as the bulk viscosity and internal thermal conductivity. The recovery of these transport coefficients in the kinetic model is crucial to accurately describe rarefied gas dynamics. For instance, the modelling of the shock wave requires correct bulk viscosity due to its high compressibility, while the modelling of thermal transpiration requires the recovery of translational thermal conductivity, rather than the total thermal conductivity (Mason 1963; Porodnov, Kulev & Tuchvetov 1978; Loyalka & Storvick 1979; Wu et al. 2020; Li et al. 2021).

Well-known kinetic models are the stochastic (Borgnakke & Larsen 1975) model in DSMC and the deterministic (Rykov 1975) and ellipsoidal-statistical (ES) BGK models (Holway 1966; Andries et al. 2000), with the emphasis to recover the transport coefficients, rather than the rates of the essential relaxation processes. To be specific, although the total thermal conductivity can be recovered in the ES-BGK model, this model cannot give correct translational and internal thermal conductivities, respectively; the Rykov model can recover each component of the thermal conductivity, and therefore, has flexibility in the simulation of thermal transpiration, but in the rarefied flow driven by the Maxwell demon, the flow velocity and heat flux are incorrect (Li et al. 2021; Zeng, Li & Wu 2022).

When radiation is concerned, transitions between different energy states of a gas molecule lead to the absorption and emission of photons, and therefore, the interactions between gas and photons have to be taken into account. The kinetic equations of the gas and photon are coupled with each other through gas–photon collision terms (Groppi & Spiga 1999), which are related to the radiation intensity and photon absorptivity that depends on the gas distribution function. Thus, simplifications of gas and photon kinetic equations are required to reduce the complexity of collision terms, both analytically and computationally.

We start from a generalized kinetic model based on each discrete vibrational energy level, which includes not only the intermolecular collisions but also the interactions between vibrational transition and radiation. By introducing the reduced velocity distribution function and total radiative intensity, the evolution of the gas flow and radiative fields is governed by four coupled equations, three for the distribution functions of gas and one for photon intensity. Meanwhile, we modify the Rykov model to recover the correct thermal relaxation rates and then give a kinetic model based on the RTA. Furthermore, by adopting the BCO for elastic intermolecular collisions, the second kinetic model, which can accurately distinguish the intermolecular potentials, is proposed.

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 1994), 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 1994; Anderson 2019). 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 1994; Nagnibeda & Kustova 2009; Anderson 2019)

(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 1999)

(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 2004). 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. 2015b) and gases with vibrational modes (Titarev & Frolova 2018).

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 1965; Mathiaud & Mieussens 2020). 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. 2015b). 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 1968a,b), 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. 2013; Wu, Reese & Zhang 2014)

(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. 2015b), 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 2009; Takata & Funagane 2011; Wu et al. 2014,  2015a), 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. 2013) 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 1999), and the Rayleigh–Brillouin scattering (Pan, Shneider & Miles 2005; Gu & Ubachs 2013), 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 2012), and the bulk viscosity can be calculated through non-equilibrium MD simulations with the expansion or compression processes (Sharma & Kumar 2019).

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 1963).

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 1962):

(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 1913) 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 1962). 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. 2021). 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 (1999),

(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 1968b), 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. 2020; Zeng et al. 2022). 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 2004). 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. Validation of the kinetic model

When the radiation is absent, the accuracy of the proposed kinetic models (2.24) and (2.31) are evaluated by comparing the numerical solutions of one-dimensional Fourier flow, Couette flow, thermal creep flow and the normal shock wave in nitrogen with constant vibrational DoF to DSMC solutions, and also compared with the experimental data of Fourier flow and normal shock waves. We use kinetic models I and II to represent the model equations with RTA and BCO, respectively. The kinetic model equations are solved by the discretized velocity method with the fast spectral method for the BCO (Wu et al. 2013,  2014), while DSMC simulations are conducted using the open source code SPARTA (Plimpton et al. 2019).

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. 2021), 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. 2021), 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 1994). 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 1986), 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 1986).

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 (1968) 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 1986).

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 1968) (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. 2021), 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 (1976). 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 (1976).

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 1963)

(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 1967; Annis & Malinauskas 1971). Parker (1959) 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 1967), and does not distinguish the difference between the translational and rotational temperatures (Lordi & Mates 1970). Nevertheless, MD simulations provide a reliable way to extract collision numbers for specific gas systems. Valentini et al. (2012) 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. (2012), 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 1962; Gorji & Jenny 2013),

(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 1976),

(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 (1976) (blues circles).

4. Rarefied molecular gas flow with radiation

The influence of thermal radiation on the total heat transfer is investigated by solving the kinetic model equation with BCO. Three typical problems are considered: one-dimensional Fourier flow, Couette flow and the normal shock wave in nitrogen flow. The vibrational DoF change with temperature as (2.3), and the effective absorptivity in the photon gray model $k_{{gray}}$ is assumed to be constant for simplicity. The transport coefficients and relaxation rates of gas are fixed as in § 3, while the parameters ${Kn}_{{photon}}$ and $\tilde {\sigma }_R$ are varied to investigate the importance of radiation in rarefied molecular gas flows.

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. Applications to two-dimensional hypersonic radiative flows

One of the most important situations in which the radiative energy interacts strongly with the rarefied gas flow is radiative shock wave passing obstacles. The gas is heated sufficiently for the radiative flux to be comparable to the flux of kinetic energy. With the ability to perform coupled simulations of rarefied gas dynamics and radiative fields, the effect of radiation on the flow structure of hypersonic non-equilibrium can be understood. Since our objective here is to examine the significance of radiative transfer, for the sake of computational simplicity, kinetic model I with RTA is adopted. The photon absorptivity is still approximated by the gray model, but its value is proportional to the product of gas number density and vibrational DoF, as indicated in (2.14a,b). When the spectrum absorptivities are available, the extension to the non-gray model of the radiation field is described by (2.26).

We consider the hypersonic gas flow with density $n_0$ at ${Ma}=15$ passing a cylinder with diameter $L_0$, the temperatures of both the incoming flow and surfaces of the cylinder are maintained at $T_0=T_{{ref}}/2$. The Knudsen number of the incoming gas flow is ${Kn}_{{gas}}=0.05$, the Knudsen number of photons at reference state ${Kn}_{{photon,ref}}$, when the gas density is $n_0$ and the vibrational DoF is 1, varies from 0.1 to 1000, and the dimensionless radiative strength $\tilde {\sigma }_R$ changes from 0.1 to 10.

The discretized velocity method is applied to solve the kinetic model equations. The simulation domain has a radius ten times larger than that of the cylinder ($r=L_0/2$). Only the upper half-domain is used in the simulation due to symmetry, which is divided into $40\times 60$ structured quadrilateral meshes with refinement near the cylinder surface. The reduced two-dimensional molecular velocity space is truncated by $[-25\sqrt {2}v_m/2,25\sqrt {2}v_m/2]\times [-7\sqrt {2}v_m,7\sqrt {2}v_m]$, and discretized uniformly by $120\times 50$ velocity points. The properties of photons over the solid angle are obtained by the Gauss–Legendre integral with $48\times 48$ discretized points. In the numerical implementation a finite volume method with a second-order reconstruction scheme is adopted. The convention fluxes are evaluated implicitly by the point relaxation technique, while the collision terms are calculated with the Venkata limiter.

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$.

6. Conclusions

In summary, we have proposed two tractable kinetic models to describe the high-temperature rarefied gas flows with radiation, where the gas flow and radiation are coupled self-consistently from the mesoscopic perspective for the first time in literature. Both the two kinetic models can recover not only the transport coefficients, but also the underlying relaxation processes. While they differ in dealing with the elastic intermolecular collisions: one uses the BCO that has the ability to capture the influence of intermolecular potentials; the other adopts RTA that has higher computational efficiency.

The accuracy of the proposed models has been assessed by comparing with DSMC simulations for one-dimensional non-radiative Fourier flow, Couette flow, creep flow driven by the Maxwell demon and normal shock wave. The kinetic model with the BCO demonstrates its excellent accuracy in all these benchmarks, while the one with RTA shows a slight discrepancy in the normal shock wave and creep flows, due to its ignoring of the velocity-dependence of intermolecular collision rates. Moreover, excellent agreements between the solutions of kinetic models and the experimental data of planar heat transfer and normal shock waves in nitrogen are obtained.

Then, the study on the effect of radiation has been conducted by solving the kinetic models in both one-dimensional typical rarefied gas flows and two-dimensional radiative hypersonic flow passing a cylinder. It is found that the presence of radiation can double the heat load on the obstacle surface compared with that with convective heat transfer only. Additionally, the influences of the photon Knudsen number and relative radiation strength are systemically investigated, and the parameter regions determining the importance of radiative heat transfer are displayed at different regimes of the photon Knudsen number.

Although only the radiation transitions of the vibrational mode are involved in the coupling of gas and radiation field in this work, a full consideration of all possible radiation transitions can be achieved straightforwardly in the framework of our general kinetic models. Meanwhile, an accurate prediction of the radiative environment requires a detailed understanding and reliable data of the radiation transitions at the molecular level, which can be acquired from the mathematical models (Zalogin et al. 2001; Babou et al. 2009), as well as experimental measurements in relevant conditions (Reynier 2021). In future work these realistic parameters will be incorporated into our general kinetic models to study the high-temperature rarefied (non-equilibrium) gas flow with strong radiation.