Kinetic modelling of rarefied gas mixtures with disparate mass in strong non-equilibrium flows

Abstract

The simulation of rarefied gas flow based on the Boltzmann equation is challenging, especially when the gas mixtures have disparate molecular masses. In this paper, a computationally tractable kinetic model is proposed for monatomic gas mixtures, to mimic the Boltzmann collision operator as closely as possible. The intra- and inter-collisions are modelled separately using relaxation approximations, to correctly recover the relaxation time scales that could span several orders of magnitude. The proposed kinetic model preserves the accuracy of the Boltzmann equation in the continuum regime by recovering four critical transport properties of a gas mixture: the shear viscosity, the thermal conductivity, the coefficients of diffusion and the thermal diffusion. While in the rarefied flow regimes, the kinetic model is found to be accurate when comparing its solutions with those from the direct simulation Monte Carlo method in several representative cases (e.g. one-dimensional normal shock wave, Fourier flow and Couette flow, two-dimensional supersonic flow passing a cylinder and nozzle flow into a vacuum), for binary mixtures with a wide range of mass ratios, species concentrations and different intermolecular potentials. Pronounced separations in species properties have been observed, and the flow characteristics of gas mixtures in shock waves are found to change as the molecular mass ratio increases from 10 to 1000.

1. Introduction

The dynamics of rarefied gas mixtures has long been an important issue, and one of the particular interests lies in the non-equilibrium phenomena of disparate-mass mixtures widely encountered in plasma physics, aerospace engineering and chip industry. For example, during the re-entry of a vehicle into the planetary atmosphere at a significantly high Mach number, plasma comprising ions, electrons, as well as neutral species, is formed between the shock and the vehicle surface, where the mass ratio between the mixture components can be as large as $10^3$ to $10^5$ (Brun 2012). In the low-pressure environment of an extreme ultraviolet (EUV) lithography system, hydrogen is commonly employed as a clean gas to effectively inhibit the diffusion of pollutant gas molecules (e.g. hydrocarbon) generated by photoresists, whose molecular mass is tens to hundreds times greater than that of clean gases (Teng et al. 2023). Similar situations are also encountered in the design of particle exhaust systems in nuclear fusion devices (Tantos et al. 2024).

In addition to the well-known rarefaction effects occurring when the mean free path of gas molecules is comparable to the characteristic flow length, multiscale non-equilibrium also exists temporally in rarefied gas mixtures with disparate mass. From the mesoscopic perspective, the light molecules have a larger thermal velocity than the heavy ones, thus leading to dispersed relaxation time scales of molecular collisions. According to Grad (1960), the light molecules will reach equilibrium among themselves first through intra-species collisions, then the heavy molecules reach their own equilibrium and, finally, all the species approach the common state through inter-species collisions. The large difference in relaxation time generates difficulties in the simulation of gas mixture flows. The Navier–Stokes-type equations involving a common flow temperature and approximated diffusion velocities of each species are adequate only when all these relaxation times are much smaller than the characteristic time of gas flow (Nagnibeda & Kustova 2009). Otherwise, the gas kinetic theory needs to be adopted to capture the non-equilibrium behaviours (Agrawal, Singh & Ansumali 2020; Sawant, Dorschner & Karlin 2020). Although the Boltzmann equation is rigorously established for gas mixtures consisting of monatomic molecules at a mesoscopic level, it is practically difficult to apply in realistic applications, due to the numerical complexity of the high-dimensional nonlinear integral collision operator, especially for mixtures with mass disparity. Even using the fast spectral method, the computational cost for each binary collision operator will be increased by the square root of the mass ratio, and hence, the numerical simulation is only conducted at a molecular mass ratio less than 36 (Wu et al. 2015b).

A well-acknowledged method for the simulations of rarefied gas flows is the direct simulation Monte Carlo (DSMC) method (Bird 1994), which models the kinetic processes of a collection of simulated particles. It has been proven by Wagner (1992) that the DSMC method is equivalent to the Boltzmann equation for a monatomic gas, as the number of simulation particles approaches infinity. Although applicable in all flow regimes, the DSMC method is computationally costly in the simulation of flows with low Knudsen numbers (${Kn}$), since the cell size and time step should be smaller than the mean free path and mean collision times of gas molecules, respectively. In a gas mixture with disparate mass, the mean collision time of different types of collisions spans multiple scales, thereby restricting the time step in the DSMC method to be smaller than the fastest relaxation time and significantly decelerating the numerical evolution of the system. In some simulation cases of collisional plasma, the electron mass is increased by three orders of magnitude to achieve an acceptable simulation time but sacrifice the accuracy (Farbar & Boyd 2010). In addition, the concentration of components in the gas mixture could vary significantly in many of the realistic applications. For instance, the ultraviolet radiation from nitric oxide is of particular interest in the flow fields surrounding hypersonic re-entry vehicles, while the mole fraction of nitric oxide is typically less than $10^{-5}$ (Erdman et al. 1993). Thus, the conventional DSMC method with equally weighted particles has enormous difficulties in terms of either huge computational costs or significant statistical noise. The differentially weighted schemes, which although solve the disparate mole fraction problem, face conservation issues during each collision (Boyd 1996), and involve a process of creating and destroying fractions of particles (Alves et al. 2018). Therefore, the multiscale feature in gas mixtures with disparate mass and concentrations makes the DSMC method time consuming and even intractable.

Therefore, it is important to develop kinetic model equations with much-simplified collision operators to imitate as closely as possible the behaviour of the Boltzmann equation, and multiscale numerical methods to solve those kinetic models. For a single-species monatomic gas, the Bhatnagar–Gross–Krook (BGK) model equation replaces the Boltzmann collision operator with a relaxation approximation (Bhatnagar, Gross & Krook 1954) to achieve high computational efficiency and lay the foundation for more sophisticated models. However, the main drawback of the BGK model is the incorrect Prandtl number produced by its single-relaxation rates for both stress and heat flux. This issue has been addressed by the modified kinetic models, e.g. the ellipsoidal statistical BGK (ES-BGK) model (Holway 1966; Andries et al. 2000; Mathiaud, Mieussens & Pfeiffer 2022) and the Shakhov model (Shakhov 1968a,b), both of which can reproduce correct shear viscosity and thermal conductivity simultaneously in the continuum flow regime. Together with the multiscale numerical methods, these kinetic models have found applications in many engineering problems (Liu, Zhu & Xu 2020; Su et al. 2020a; Pfeiffer, Garmirian & Gorji 2022; Zeng, Su & Wu 2023a; Liu et al. 2024).

However, the extension of the single-species kinetic models to gas mixtures is a non-trivial task. First, collisions between different species of molecules lead to exchanges of momentum and energy between mixture components, exhibiting notable disparities in relaxation times due to variations in species properties. Second, apart from the component-specific shear viscosities and thermal conductivities, a multi-species gas also possesses effective mixture viscosity and thermal conductivity, as well as diffusion and thermal-diffusion coefficients that correspond to the Fick and Soret effects, respectively (Chapman & Cowling 1970), which have to be recovered by the kinetic models in the continuum limit. Previously proposed kinetic models using BGK approaches can be classified into two types (Pirner 2021). One uses a single-relaxation term involving both inter- and intra-species collisions, and the other one has a sum of collision terms modelling each type of collision individually. For the single-relaxation models, the model of Andries, Aoki & Perthame (2002) is the most widely applied one, which reduces to the single-species BGK model for mechanically identical components that cannot recover shear viscosity and thermal conductivity simultaneously. Besides, the diffusion coefficient is not correctly captured. Later, more adjustable parameters are introduced into the single-relaxation model by adopting ES-BGK and Shakhov-type operators, thus, thermal conductivity (Brull 2015) and diffusion coefficients (Groppi, Monica & Spiga 2011; Todorova & Steijl 2019) can be recovered. It should be noted that, however, the proposed models are incapable of modelling the thermally induced flow of non-Maxwell molecules, since the thermal-diffusion effect is not reproduced, which is an important transport phenomenon in gas mixtures (Sharipov 2024). Although the single-relaxation models are found to be accurate for mixtures with a small mass ratio (Pfeiffer, Mirza & Nizenkov 2021) and also easy to be extended to polyatomic gases with internal energy (Bisi & Travaglini 2020; Todorova, White & Steijl 2020) and chemical reactions (Groppi & Spiga 2004; Bisi, Monaco & Soares 2018), they are not able to distinguish the multiple scales of relaxation times and different types of interactions. On the other hand, several multi-relaxation models have been proposed (Morse 1964; Hamel 1965; Greene 1973; Haack, Hauck & Murillo 2017; Klingenberg, Pirner & Puppo 2017; Bobylev et al. 2018; Bisi et al. 2022), which mainly differ in the construction of the auxiliary properties in inter-species collision terms. Because of the complexity of the collision operators, these models were basically derived mathematically but rarely applied to realistic problems (Tantos, Varoutis & Day 2021; Bisi et al. 2022). More importantly, the recovery of transport properties from the multi-relaxation models in the continuum limit is usually overlooked, and the determination of the adjustable parameters (e.g. the multiple relaxation times) is still an open question. It is also noted that, the models proposed by McCormack (1973) and Kosuge (2009) replace the Boltzmann collision operator by the use of polynomial expansions in the molecular velocity, where the model coefficients are determined by matching the moments of the model collision operator to the Boltzmann collision operator. The polynomial models correctly recover all the transport coefficients of a gas mixture and show good performance in slightly to moderately non-equilibrium flows (Ho et al. 2016; Tantos & Valougeorgis 2018), while accurate intermolecular potentials are required when applying to realistic applications. However, as noted by Kosuge (2009) himself, the polynomial model is not well suited for simulating strong non-equilibrium mixture flows, where the reference velocity distributions in the collision terms may exhibit remarkably negative values, leading to an unphysical prediction of the macroscopic quantities.

Despite the great effort made in the past decades, establishing accurate and computationally tractable kinetic models for gas mixtures with disparate mass is still of significant challenge. The present work is dedicated to developing a kinetic model based on the idea of multi-relaxation models for rarefied monatomic gas mixtures with disparate molecular mass, which not only recovers the transport properties of a gas mixture including shear viscosity, thermal conductivity, diffusion coefficient and thermal-diffusion coefficients, but also correctly captures the multiscale relaxation rates of different collision processes. More importantly, with the deterministic numerical methods and multiscale schemes, the proposed kinetic model can be used to solve gas mixture flows with disparate mass and mole fractions efficiently and accurately in all flow regimes.

The rest of the paper is organized as follows. In § 2 the kinetic model is proposed with all the adjustable parameters determined by the transport properties of the gas, and the transport coefficients and hydrodynamic equations in the continuum limit are given. In § 3, to make a consistent comparison to the DSMC method, the kinetic model parameters are obtained based on the DSMC collision model. The accuracy of the proposed kinetic model is assessed by the DSMC method in several one-dimensional and two-dimensional problems in §§ 4 and 5, respectively, and the flow characteristics of gas mixtures with disparate mass are discussed. Finally, conclusions are presented in § 6.

2. Kinetic model equation

The gas kinetic theory describes the status of a gaseous system in the phase space using the velocity distribution functions, whose evolution is governed by the Boltzmann equations when only binary collisions are considered. Because of the unaffordable computational cost of solving the $N^2$ number of Boltzmann operators for an $N$-component gas mixture, kinetic model equations with simplified collision operators are highly demanded. Theoretically, several requirements need to be followed in the kinetic modelling of monatomic gas mixtures: (i) the collision terms satisfy the conservations of mass, momentum and energy, and restore the equilibrium state for an isolated system; (ii) the relaxation rates for different types of intermolecular interaction can be correctly captured; (iii) all the transport properties given by the model equation are consistent with those obtained from the Boltzmann equation in the hydrodynamic limit; (iv) the momentum and energy exchange during inter-species collisions are close to those obtained from the Boltzmann equation; and (v) the kinetic model complies with the indifferentiability principle and H-theorem.

However, it is impractical to build such an ideal kinetic model to satisfy all the requirements with affordable computational cost. Even for a single-species gas, although the ES-BGK model has been proven to keep the non-negativity of the velocity distribution functions and satisfy the H-theorem (Andries et al. 2000), the Shakhov model is found to be more accurate for many strong non-equilibrium problems due to its better approximations of high-order moments of molecular velocity distributions (Chen, Xu & Cai 2015; Fei et al. 2020; Yuan & Wu 2022; Park et al. 2024). Moreover, despite the fact that the model of Andries et al. (2002) complies with the indifferentiability principle unconditionally, it recovers only one transport coefficient and ignores the multiple relaxation time scales.

Therefore, to achieve a balance between the accuracy and computational burden for a gas mixture, we build our kinetic model using multi-relaxation operators to distinguish different types of binary interactions. The tunable parameters in the kinetic model are determined by recovering the transport properties in the continuum flow regime, including viscosity, thermal conductivity, diffusion coefficients and thermal-diffusion coefficients. Additionally, the indifferentiability principle can be satisfied in the near-equilibrium condition, where the proposed kinetic model reduces to the Shakhov model for a mixture with identical components. The H-theorem is not proven for the present model in this paper.

Note that our kinetic model employs expansions of molecular thermal velocity in the reference distribution functions and, thus, differs from the previous multi-relaxation ES-BGK models. Also, our model is different to those proposed by McCormack (1973) and Kosuge (2009), where polynomial expansions around the equilibrium state of the entire mixture are adopted. Thus, they may not accurately capture the velocity distributions of individual gas components in a mixture with substantial differences in velocity and temperature among the species, and can result in notably negative values for the reference velocity distributions. In contrast, our model differs in its construction of the reference velocity distributions, which are expansions based on the equilibrium states of species temperatures and auxiliary velocities of each respective binary collision. Thus, the proposed model is expected to perform better in strong non-equilibrium flows; see Appendix C for an example.

2.1. Kinetic description of monatomic gas mixture

We consider an $N$-components mixture of monatomic gases with the velocity distribution functions $f_s(\boldsymbol {x},\boldsymbol {v},t)$ describing their mesoscopic states, where $s$ indicates the species, $t$ is the time, $\boldsymbol {x}\in \mathbb {R}^3$ is the spatial coordinates and $\boldsymbol {v}\in \mathbb {R}^3$ is the molecular velocity. Since all the collisions, either intra-species or inter-species, are binary, the distribution function for species $s$ under external force $\boldsymbol {F}_s$ evolve according to the Boltzmann equation,

(2.1)\begin{equation} \underbrace{\frac{\partial{f_{s}}}{\partial{t}}+ \boldsymbol{v}\boldsymbol{\cdot}\frac{\partial{f_{s}}}{\partial{\boldsymbol{x}}}+ \frac{\boldsymbol{F}_s}{m_s}\boldsymbol{\cdot} \frac{\partial{f_{s}}}{\partial{\boldsymbol{v}}}}_{\mathcal{D}f_{s}} = \sum_{r=1}^N Q_{sr}(\,f_s,f_r), \quad s=1,2,\ldots,N, \end{equation}

where $m_s$ is the molecular mass; the left-hand side is known as the streaming term $\mathcal {D}f_{s}$ and the right-hand side is a sum over all binary Boltzmann collision operators $Q_{sr}$ between molecules of species $s$ and $r$,

(2.2)\begin{equation} Q_{sr}(f_s,f_r) = \int_{\mathbb{R}^3}\int_{4{\rm \pi}} |\boldsymbol{v}_{*}|\sigma_{sr}\left(|\boldsymbol{v}_{*}|, \boldsymbol{\varOmega}\boldsymbol{\cdot}\frac{\boldsymbol{v}_{*}}{|\boldsymbol{v}_{*}|} \right)[f_{s}(\boldsymbol{v}')f_{r}(\boldsymbol{w}')- f_{s}(\boldsymbol{v})f_{r}(\boldsymbol{w})]\, \mathrm{d}\boldsymbol{\varOmega}\,\mathrm{d}\boldsymbol{w}. \end{equation}

Here, $\sigma _{sr}$ is the differential scattering cross-section depending on the intermolecular potential between the two species and $\boldsymbol {\varOmega }$ is the unit vector of the solid angle; $\boldsymbol {v}$ and $\boldsymbol {w}$ are the pre-collision velocities of the two molecules of species $s$ and $r$, respectively. Hence, $\boldsymbol {v}_{*}=\boldsymbol {v}-\boldsymbol {w}$ is the relative velocity, which determines the post-collision velocities $\boldsymbol {v}'$ and $\boldsymbol {w}'$ of the collision pair as

(2.3a,b)\begin{equation} \boldsymbol{v}' = \boldsymbol{v} - \frac{2m_r}{m_s+m_r}(\boldsymbol{v}_{*} \boldsymbol{\cdot}\boldsymbol{\varOmega})\boldsymbol{\varOmega}, \quad \boldsymbol{w}' = \boldsymbol{w} + \frac{2m_s}{m_s+m_r}(\boldsymbol{v}_{*} \boldsymbol{\cdot}\boldsymbol{\varOmega})\boldsymbol{\varOmega}. \end{equation}

The macroscopic variables of each species $s$, namely, the number density $n_s$, mass density $\rho _s$, flow velocity $\boldsymbol {u_s}$, temperatures $T_s$, pressure tensor $\boldsymbol{\mathsf{P}}_s$ and heat flux $\boldsymbol {q}_s$, are obtained by taking the moments of the respective velocity distribution function $f_s$,

(2.4)\begin{equation} \left.\begin{array}{@{}c@{}} n_s = \langle 1,f_s\rangle , \quad \rho_s = \langle m_s,f_s\rangle , \quad \rho_s\boldsymbol{u}_s = \langle m_s\boldsymbol{v},f_s\rangle , \quad \tfrac{3}{2}n_sk_BT_s = \langle \tfrac{1}{2}m_s |\boldsymbol{v}-\boldsymbol{u}_s|^2,f_s\rangle , \\ \boldsymbol{\mathsf{P}}_s = \langle m_s(\boldsymbol{v}-\boldsymbol{u}_s)(\boldsymbol{v}-\boldsymbol{u}_s),f_s\rangle , \quad \boldsymbol{q}_s = \langle \tfrac{1}{2}m_s|\boldsymbol{v}-\boldsymbol{u}_s|^2(\boldsymbol{v}- \boldsymbol{u}_s),f_s\rangle, \end{array}\right\} \end{equation}

where $k_B$ is the Boltzmann constant; and the operator $\langle h,\psi \rangle$ is defined as an integral of $h\psi$ over the velocity space,

(2.5)\begin{equation} \langle h,\psi\rangle \equiv \int_{\mathbb{R}^3}{h\psi}\,\mathrm{d}\boldsymbol{v}. \end{equation}

Then, the corresponding macroscopic quantities for the mixture, the number density $n$, mass $\rho$, flow velocity $\boldsymbol {u}$, temperatures $T$, pressure tensor $\boldsymbol{\mathsf{P}}$ and heat flux $\boldsymbol {q}$, are given by

(2.6) \begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle n = \sum_s\langle 1,f_s\rangle = \sum_{s}n_s, \\ \displaystyle \rho = \sum_s\langle m_s,f_s\rangle = \sum_{s}\rho_s, \\ \displaystyle \rho\boldsymbol{u} = \sum_s\langle m_s\boldsymbol{v},f_s\rangle = \sum_{s}\rho_s\boldsymbol{u}_s, \\ \displaystyle \dfrac{3}{2}nk_BT = \sum_{s}\left\langle \dfrac{1}{2}m_sc^2,f_s \right\rangle =\sum_{s}\dfrac{3}{2}n_sk_BT_s + \dfrac{1}{2}\sum_{s}\rho_s |\boldsymbol{u}_s-\boldsymbol{u}|^2, \\ \displaystyle \boldsymbol{\mathsf{P}} = \sum_{s}\langle m_s \boldsymbol{c}\boldsymbol{c},f_s\rangle = \sum_{s}\boldsymbol{\mathsf{P}}_s + \sum_{s}\rho_s (\boldsymbol{u}_s-\boldsymbol{u})(\boldsymbol{u}_s-\boldsymbol{u}), \\ \begin{aligned} \boldsymbol{q} & = \sum_{s}\left\langle \dfrac{1}{2}m_sc^2\boldsymbol{c},f_s\right\rangle = \sum_{s}\boldsymbol{q}_s + \sum_{s}\dfrac{3}{2}n_sk_BT_s(\boldsymbol{u}_s-\boldsymbol{u})\\ & \quad +\dfrac{1}{2}\sum_{s}\rho_s|\boldsymbol{u}_s-\boldsymbol{u}|^2 (\boldsymbol{u}_s-\boldsymbol{u}) + \sum_{s}\boldsymbol{\mathsf{P}}_s \boldsymbol{\cdot}(\boldsymbol{u}_s-\boldsymbol{u}), \end{aligned} \end{array}\right\} \end{equation}

where $\boldsymbol {c}=\boldsymbol {v}-\boldsymbol {u}$ is the peculiar velocity with respect to the mixture velocity $\boldsymbol {u}$, and therefore, the diffusion velocities $\boldsymbol {u}_s-\boldsymbol {u}$ contribute to the mixture temperature, stress and heat flux. Note that the scalar $c$ denotes the magnitude of the vector $\boldsymbol {c}$, with $c=|\boldsymbol {c}|$, and this notation is consistently used for both the molecular velocity $\boldsymbol {v}$ and the macroscopic velocity $\boldsymbol {u}$ throughout this paper.

2.2. Kinetic model with multi-relaxation collision operators

The proposed kinetic model for gas mixtures adopts relaxation time approximations for each pair of gas components individually to simplify the Boltzmann collision operator (2.2), and thus, the evolution of distribution functions $f_{s}$ can be written as

(2.7)\begin{equation} {\mathcal{D}f_{s}} = {\sum_{r=1}^N\frac{1}{\tau_{sr}}(g_{sr}-f_s)}, \quad s=1,2,\ldots,N, \end{equation}

where $s=r$ indicates a intra-species collision operator, $s\neq r$ are inter-species collisions; $\tau _{sr}$ is the corresponding relaxation time and $g_{sr}$ is the reference distribution function constructed in the form

(2.8) \begin{align} g_{sr}&=\hat{n}_{sr}\left(\frac{m_s}{2{\rm \pi} k_B{T}_{s}}\right)^{3/2} \exp\left(-\frac{m_s(\boldsymbol{v}-\hat{\boldsymbol{u}}_{sr})^2}{2k_B {T}_{s}}\right) \times \left[1+\frac{\hat{T}_{sr}-T_s}{T_s} \left(\frac{m_s(\boldsymbol{v}-\hat{\boldsymbol{u}}_{sr})^2}{2k_B {T}_{s}}-\frac{3}{2}\right)\right.\nonumber\\ &\quad +\left.\frac{2m_s\hat{\boldsymbol{q}}_{sr} \boldsymbol{\cdot}(\boldsymbol{v}-\hat{\boldsymbol{u}}_{sr})}{5\hat{n}_{sr} k_B^2\hat{T}_{s}^2}\left(\frac{m_s(\boldsymbol{v}- \hat{\boldsymbol{u}}_{sr})^2}{2k_B\hat{T}_{s}}- \frac{5}{2}\right)\right], \end{align}

with $\hat {n}_{sr}, \hat {T}_{sr}, \hat {\boldsymbol {u}}_{sr}, \hat {\boldsymbol {q}}_{sr}$ being the auxiliary parameters.

Construction of the auxiliary parameters is the most crucial task in building the kinetic model. As the essential constraints, the conservations of mass, momentum and energy have to be guaranteed during any binary collisions. For the intra-species collision operators, the conservations are the same as those for a single-species gas, i.e.

(2.9a–c)\begin{align} \left\langle 1,\frac{1}{\tau_{ss}}(g_{ss}-f_{s})\right\rangle = 0, \quad \left\langle m_s\boldsymbol{v},\frac{1}{\tau_{ss}}(g_{ss}-f_{s})\right\rangle = 0, \quad \left\langle \frac{1}{2}m_sv^2,\frac{1}{\tau_{ss}}(g_{ss}-f_{s})\right\rangle = 0. \end{align}

Therefore, the auxiliary number density, flow velocity and temperature in intra-species collision operators can be simply determined by the macroscopic properties of each species as $\hat {n}_{ss}=n_s$, $\hat {\boldsymbol {u}}_{ss}=\boldsymbol {u}_s$, $\hat {T}_{ss}=T_s$. Thus, the reference distribution function $g_{ss}$ in the intra-species collision term reduces to that in the Shakhov model.

The inter-species collision operators conserve the number density of each inert species, and the total momentum and energy of the collision pairs of species $s$ and $r$:

(2.10) \begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \left\langle 1,\dfrac{1}{\tau_{sr}}(g_{sr}-f_{s})\right\rangle = 0, \quad \left\langle 1,\dfrac{1}{\tau_{rs}}(g_{rs}-f_{r})\right\rangle = 0, \\ \displaystyle \left\langle m_s\boldsymbol{v},\dfrac{1}{\tau_{sr}}(g_{sr}-f_{s})\right\rangle + \left\langle m_r\boldsymbol{v},\dfrac{1}{\tau_{rs}}(g_{rs}-f_{r})\right\rangle = 0, \\ \displaystyle \left\langle \dfrac{1}{2}m_sv^2,\dfrac{1}{\tau_{sr}}(g_{sr}-f_{s})\right\rangle + \left\langle \dfrac{1}{2}m_rv^2,\dfrac{1}{\tau_{rs}}(g_{rs}-f_{r})\right\rangle = 0. \end{array}\right\} \end{equation}

Then, the auxiliary number density is obtained as $\hat {n}_{sr} = n_s$ and $\hat {n}_{rs} = n_r$, while the auxiliary velocities and temperatures $\hat {\boldsymbol {u}}_{sr}, \hat {\boldsymbol {u}}_{rs}, \hat {T}_{sr}, \hat {T}_{rs}$ cannot be uniquely determined, but yield the constraints

(2.11)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \dfrac{\rho_s}{\tau_{sr}}\hat{\boldsymbol{u}}_{sr} +\dfrac{\rho_r}{\tau_{rs}}\hat{\boldsymbol{u}}_{rs} = \dfrac{\rho_s}{\tau_{sr}}{\boldsymbol{u}}_{s} +\dfrac{\rho_r}{\tau_{rs}}{\boldsymbol{u}}_{r}, \\ \displaystyle\dfrac{1}{\tau_{sr}}\left[\dfrac{3}{2}n_sk_B(\hat{T}_{sr}- T_s)+\dfrac{1}{2}\rho_s(\hat{u}^2_{sr}-u^2_s)\right]\\ \quad + \dfrac{1}{\tau_{rs}}\left[\dfrac{3}{2}n_rk_B(\hat{T}_{rs}-T_r)+ \dfrac{1}{2}\rho_r(\hat{u}^2_{rs}-u^2_r)\right]=0. \end{array}\right\} \end{equation}

Therefore, we further impose the following assumptions to determine the auxiliary velocities and temperatures in inter-species collision operators:

(2.12)\begin{equation} \left.\begin{array}{@{}c@{}} \hat{\boldsymbol{u}}_{sr}-\hat{\boldsymbol{u}}_{rs} = (1-a_{sr})({\boldsymbol{u}}_{s}-{\boldsymbol{u}}_{r}) -b_{sr}(\boldsymbol{\nabla}{\ln T_s}+\boldsymbol{\nabla}{\ln T_r}), \\ \hat{T}_{sr}-\hat{T}_{rs} = (1-c_{sr})(T_{s}-T_{r}) -d_{sr}|{\boldsymbol{u}}_{s}-{\boldsymbol{u}}_{r}|^2. \end{array}\right\} \end{equation}

Here $a_{sr}=a_{rs}, b_{sr}=-b_{rs}, c_{sr}=c_{rs}, d_{sr}=-d_{rs}$ are the adjustable parameters, which describe how rapidly the equilibrium among different gas components can be achieved through inter-species collisions. To be specific, the assumptions imposed to determine the auxiliary velocities and temperatures are designed to recover the correct momentum and energy relaxation, respectively. Since the momentum relaxation is closely related to the diffusion phenomenon, the term with parameters $a_{sr}$ is presented to account for the ordinary diffusion arising from concentration gradients, that is, the Fick effect. The term with parameters $b_{sr}$ is introduced to encapsulate the thermal diffusion induced by temperature gradients, known as the Soret effect. Note that, compared with those in the literature (such as Haack et al. 2017; Bobylev et al. 2018), (2.12) gives a more general expression for the relations between the auxiliary and macroscopic properties, by adding the temperature gradient to phenomenologically model the thermally induced flow. Then, combined with (2.11), the auxiliary velocities and temperatures are given by

(2.13)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \hat{\boldsymbol{u}}_{sr} ={\boldsymbol{u}}_{s}- \dfrac{\rho_r\tau_{sr}}{\rho_s\tau_{rs}+\rho_r\tau_{sr}}\boldsymbol{{X}}_{sr},\\ \displaystyle \hat{T}_{sr} = T_s-\dfrac{n_r\tau_{sr}}{n_s\tau_{rs}+n_r \tau_{sr}}{Y}_{sr}-\dfrac{\rho_s\rho_r\tau_{sr}\tau_{rs}\boldsymbol{{X}}_{sr} \boldsymbol{\cdot}[\boldsymbol{{X}}_{sr}-2({\boldsymbol{u}}_{s}- {\boldsymbol{u}}_{r})]}{3k_B(n_s\tau_{rs}+ n_r\tau_{sr})(\rho_s\tau_{rs}+\rho_r\tau_{sr})}, \end{array}\right\} \end{equation}

with

(2.14)\begin{equation} \left.\begin{array}{@{}c@{}} \boldsymbol{{X}}_{sr} = a_{sr}({\boldsymbol{u}}_{s}-{\boldsymbol{u}}_{r}) +b_{sr}(\boldsymbol{\nabla}{\ln T_s}+\boldsymbol{\nabla}{\ln T_r}), \\ {Y}_{sr} = c_{sr}(T_{s}-T_{r})+d_{sr} |{\boldsymbol{u}}_{s}-{\boldsymbol{u}}_{r}|^2. \end{array}\right\} \end{equation}

In addition, the auxiliary properties $\hat {\boldsymbol {q}}_{sr}$ are constructed to adjust the relaxation rates of heat fluxes as

(2.15)\begin{equation} \hat{\boldsymbol{q}}_{sr} = (1-{Pr}_{sr}){\boldsymbol{q}}_{s}+\gamma_{sr} ({\boldsymbol{q}}_{sr}-{\boldsymbol{q}}_{s}), \end{equation}

where ${Pr}_{sr}$ is an effective Prandtl number giving the thermal relaxation of species $s$ due to collisions with species $r$, which reduces to ${Pr}_{ss}=2/3$ for the intra-species collisions of a monatomic gas; ${\boldsymbol {q}}_{sr}$ is defined as the heat flux of species $s$ measured relative to auxiliary velocity $\hat {\boldsymbol {u}}_{sr}$,

(2.16)\begin{equation} {\boldsymbol{q}}_{sr} = \langle \tfrac{1}{2}m_s |\boldsymbol{v}-\hat{\boldsymbol{u}}_{sr}|^2(\boldsymbol{v}- \hat{\boldsymbol{u}}_{sr}),f_s\rangle , \end{equation}

and thus, (2.15) yields $\hat {\boldsymbol {q}}_{ss} = (1-{Pr}_{ss}){\boldsymbol {q}}_{s}$ for the intra-species collision term; $\gamma _{sr}=-\gamma _{rs}$ is a dimensionless coefficient taking into account the Dufour effects caused by the diffusive thermal conductivity. Based on the asymptotic analysis of the Boltzmann equation (Chapman & Cowling 1970), it is an effect inverse to thermal diffusion and, hence, $\gamma _{sr}$ is not an independent parameter but can be determined by the other adjustable parameters, as shown in the following section.

2.3. Continuum flow limit

The multi-relaxation kinetic model equations have been constructed for a monatomic gas mixture, which contains several adjustable parameters: relaxation rates $\tau _{sr}$, effective Prandtl number ${Pr}_{sr}$ and the set of coefficients $a_{sr}$, $b_{sr}$, $c_{sr}$, $d_{sr}$, $\gamma _{sr}$ used for calculating auxiliary properties. Here, we perform the Chapman–Enskog analysis to get the asymptotic limit of the proposed model equation and determine the relevant parameters to recover the transport coefficients, so that the kinetic model and macroscopic fluid dynamics are consistent in the continuum flow regime, which is the basic requirement to achieve accurate kinetic modelling.

Without loss of generality, we consider a binary mixture of monatomic gases. The extension to a multi-species mixture is expected to be relatively easy, due to the fact that our kinetic model is built on each binary collision pair. Therefore, although the expressions of the transport coefficients of the whole mixture become more complex for mixtures with more than two components, only the properties of binary mixtures are required to determine all the model parameters in each possible binary collision term.

In the continuum limit, when all the relaxation times are considerably smaller than the characteristic time of the gas flow, only the species number density $n_s$, total momentum $\rho \boldsymbol {u}$ and energy $\frac {3}{2}nk_BT+\frac {1}{2}\rho u^2$ are the collisional invariants. Then, the set of equations for the conserved macroscopic variables $n_s, \boldsymbol {u}, T$ can be obtained by taking momentums of the kinetic equations (2.7) and summing over the species ($s=1,2$) for momentum and energy equations,

(2.17) \begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \dfrac{\partial{n_s}}{\partial{t}} + \boldsymbol{\nabla}\boldsymbol{\cdot} (n_s\boldsymbol{u}) + \boldsymbol{\nabla}\boldsymbol{\cdot} (n_s(\boldsymbol{u}_s-\boldsymbol{u})) = 0, \\ \displaystyle \dfrac{\partial}{\partial{t}}(\rho\boldsymbol{u}) + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u}\boldsymbol{u}) + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{P}} = n_1\boldsymbol{F}_1+n_2\boldsymbol{F}_2, \\ \begin{aligned} & \dfrac{\partial}{\partial{t}} \left(\dfrac{3}{2}nk_BT+\dfrac{1}{2}\rho u^2\right) +\boldsymbol{\nabla}\boldsymbol{\cdot}\left[\left(\dfrac{3}{2}nk_BT+\dfrac{1}{2}\rho u^2\right)\boldsymbol{u}\right] \\ & \quad +\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\mathsf{P}}\boldsymbol{\cdot}\boldsymbol{u}) + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q} = n_1\boldsymbol{F}_1\boldsymbol{\cdot}\boldsymbol{u}_1+n_2\boldsymbol{F}_2\boldsymbol{\cdot}\boldsymbol{u}_2, \end{aligned} \end{array}\right\} \end{equation}

where $\boldsymbol {u}_s$ is the species macroscopic velocity defined in (2.4), $\boldsymbol{\mathsf{P}}$ and $\boldsymbol {q}$ are the mixture pressure tensor and heat flux, respectively, given in (2.6). To close the above set of equations, it is necessary to find the approximations to distribution functions $f_s$, and hence, the properties $\boldsymbol {u}_s, \boldsymbol{\mathsf{P}}, \boldsymbol {q}$ can be expressed as functions of the macroscopic variables $n_s, \boldsymbol {u}, T$, giving the constitutive relations and transport properties of the gas mixture (see the details in Appendix A).

The Navier–Stokes-type equation for a gas mixture is then obtained by the second approximation of $f_s$, and the species velocities $\boldsymbol {u}_s^{{NS}}$, stress tensor $\boldsymbol{\mathsf{P}}^{{NS}}$ and heat flux $\boldsymbol {q}^{{NS}}$ are given by

(2.18) \begin{align} \left.\begin{array}{@{}c@{}} \displaystyle \boldsymbol{u}_1^{{NS}} =\boldsymbol{u} -\dfrac{\rho_1\tau_{21}+\rho_2\tau_{12}}{a_{12}} \dfrac{p}{\rho\rho_1}\boldsymbol{d}_{12}- \dfrac{2b_{12}\rho_2}{a_{12}\rho}\boldsymbol{\nabla}{\ln T}, \\ \displaystyle \boldsymbol{u}_2^{{NS}} = \boldsymbol{u}+ \dfrac{\rho_1\tau_{21}+\rho_2\tau_{12}}{a_{12}} \dfrac{p}{\rho\rho_2}\boldsymbol{d}_{12}+\dfrac{2b_{12}\rho_1}{a_{12} \rho}\boldsymbol{\nabla}{\ln T}, \\ \displaystyle \boldsymbol{\mathsf{P}}^{{NS}} = nk_BT\boldsymbol{\mathsf{I}} - k_BT\left(\dfrac{\tau_{11}\tau_{12}}{\tau_{11}+\tau_{12}}n_1+ \dfrac{\tau_{22}\tau_{21}}{\tau_{22}+\tau_{21}}n_2\right) \left(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{\mathrm{T}}- \dfrac{2}{3}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}\boldsymbol{\mathsf{I}}\right) , \\ \begin{aligned} \boldsymbol{q}^{{NS}} & = \dfrac{5}{2}k_BT(n_1\boldsymbol{u}_1^{{NS}} +n_2\boldsymbol{u}_2^{{NS}}-n\boldsymbol{u}) + \dfrac{2b_{12}\rho_1\rho_2}{\rho_1\tau_{21}+\rho_2\tau_{12}} ({\boldsymbol{u}}_{1}^{{NS}}-{\boldsymbol{u}}_{2}^{{NS}}) \\ & \quad -\left[\left(\dfrac{n_1}{m_1}\dfrac{\tau_{11}\tau_{12}}{{Pr}_{12} \tau_{11}+{Pr}_{11}\tau_{12}}+\dfrac{n_2}{m_2}\dfrac{\tau_{22} \tau_{21}}{{Pr}_{21}\tau_{22}+{Pr}_{22}\tau_{21}}\right) \dfrac{5k_B^2T}{2} \right.\\ & \quad - \left.\dfrac{4b_{12}^2\rho_1\rho_2}{a_{12} (\rho_1\tau_{21}+\rho_2\tau_{12})T}\right]\boldsymbol{\nabla} T, \end{aligned} \end{array}\right\} \end{align}

where $p=nk_BT$ is the pressure, $\boldsymbol{\mathsf{I}}$ is a $3\times 3$ identity matrix and $\boldsymbol {d}_{12}$ is given as

(2.19)\begin{equation} \boldsymbol{d}_{12} = \boldsymbol{\nabla} \left(\frac{n_1}{n}\right) +\frac{n_1n_2(m_2-m_1)}{n\rho}\boldsymbol{\nabla}{\ln p} -\frac{\rho_1\rho_2}{\rho p}\left(\frac{\boldsymbol{F}_1}{m_1}-\frac{\boldsymbol{F}_2}{m_2}\right), \end{equation}

which yields $\boldsymbol {d}_{12}=-\boldsymbol {d}_{21}$, since $\boldsymbol {\nabla } ({n_1}/{n}) =-\boldsymbol {\nabla } ({n_2}/{n})$.

2.3.1. Viscosity

The shear viscosity $\mu$ of a binary mixture can be obtained from the off-diagonal components of the non-equilibrium stress tensor $\boldsymbol{\mathsf{P}}$ in (2.18), which leads to

(2.20)\begin{equation} \mu = k_BT\left(\frac{\tau_{11}\tau_{12}}{\tau_{11}+ \tau_{12}}n_1+\frac{\tau_{22}\tau_{21}}{\tau_{22}+\tau_{21}}n_2\right). \end{equation}

Clearly, the mixture shear stress depends on all the relaxation times. It is known that the dependence of species shear viscosity $\mu _s$ on mean collision time $\tau _{ss}$ of intra-species collisions is $\mu _s=n_sk_BT\tau _{ss}$. However, the relaxation time $\tau _{sr} (s\neq r)$ in the kinetic model is no longer a mean molecular collision time of inter-species interactions, but measures the time scale approaching the reference states $g_{sr}$ for the component $s$ due to the collisions with component $r$. In other words, $n_r\tau _{sr}$ should be smaller than $n_s\tau _{rs}$ when $m_s< m_r$, because the collisions between molecules with disparate mass have a more significant influence on the light one than on the heavy one. For example, by considering a mixture consisting of ions and electrons with mass $m_i$ and $m_e$, respectively, the approximated relation between the inter-species relaxation times given in the literature is $n_i\tau _{ei}/n_e\tau _{ie}=m_e/m_i\ll 1$ (Bellan 2006).

Therefore, here we define

(2.21)\begin{equation} \phi_{sr}=\frac{n_s\tau_{ss}}{n_r\tau_{sr}}, \end{equation}

to quantify the ratio between the relaxation times of each intra- and inter-species collision for species $s$, where $\phi _{sr}$ is only determined by the intermolecular interactions and temperature, but independent of the concentrations of the components. Then the mixture viscosity (2.20) can be rewritten in the form

(2.22)\begin{equation} \mu = \frac{\mu_1}{1+\dfrac{n_2}{n_1}\phi_{12}} + \frac{\mu_2}{1+\dfrac{n_1}{n_2}\phi_{21}}, \end{equation}

which is a linear combination of the species viscosities of the mixture components and shares the exact same form as that given by Wilke's mixture rule (Wilke 1950). Nevertheless, the more accurate values of $\phi _{sr}$ can be obtained by fitting the mixture viscosity measured experimentally.

2.3.2. Diffusion

In a gas mixture the ordinary diffusion coefficient $D_{sr}$ is usually measured when the gas mixture is uniform in temperature and pressure and without external forces acting on the gas molecules; while the thermal diffusion with coefficient $D_{T,sr}$ also contributes to the diffusion velocity when there is a temperature gradient present. In general, the mass flux $J_s$ of the species $s$ in a binary mixture caused by gradients of concentration, pressure and temperature is given by (Hirschfelder, Curtiss & Bird 1954)

(2.23)\begin{equation} J_s =\rho_s(\boldsymbol{u}_s-\boldsymbol{u}) ={-}\frac{n^2}{\rho}m_sm_r[D_{sr}\boldsymbol{d}_{sr}-D_{T,sr}\boldsymbol{\nabla}{\ln T}]. \end{equation}

Meanwhile, the diffusion velocity in the proposed kinetic model is obtained as (2.18) from the Chapman–Enskog method,

(2.24)\begin{equation} \boldsymbol{u}_s-\boldsymbol{u} ={-}\frac{\rho_1\tau_{21}+ \rho_2\tau_{12}}{a_{sr}}\frac{p}{\rho_1\rho_2}\boldsymbol{d}_{sr}- \frac{2b_{sr}}{a_{sr}}\boldsymbol{\nabla}{\ln T}. \end{equation}

Therefore, the binary diffusion coefficient $D_{12}$ is

(2.25)\begin{equation} D_{12} = \frac{k_BT}{m_1m_2n}\frac{\rho_1\tau_{21}+\rho_2\tau_{12}}{a_{12}}, \end{equation}

which determines the parameter $a_{12}=a_{21}$ when the relaxation times are known from viscosity. Also, the thermal-diffusion coefficient $D_{T,12}$, as well as the thermal-diffusion ratio $k_{T,12}=D_{T,12}/D_{12}$ are obtained:

(2.26a,b)\begin{equation} D_{T,12} = \frac{2b_{12}n_1n_2}{a_{12}n^2}, \quad k_{T,12} =\frac{2b_{12}\rho_1\rho_2}{p(\rho_1 \tau_{21}+\rho_2\tau_{12})}. \end{equation}

Thus, the parameter $b_{12}=-b_{21}$ is adjusted to match the thermal-diffusion ratio and recover the Soret effect in the continuum limit. The Soret effect describes mass separation due to a temperature gradient and it may be positive or negative depending on the mass difference and intermolecular potentials (Chapman & Cowling 1970). In general, thermal diffusion becomes stronger in molecules with a larger mass difference and, thus, is crucial in the modelling of gas mixtures with significant mass disparity. Meanwhile, although the absolute value of the thermal-diffusion ratio $k_{T,12}$ usually has an order of magnitude less than $10^{-1}$ for a neutral gas mixture, it can be greatly increased in ionized gases (Chapman 1958).

2.3.3. Thermal conductivity

In addition to the direct transfer of kinetic energy during collisions, the diffusional migration of molecules also carries thermal energy and contributes to the heat transport in a gas mixture, which is known as the Dufour effect. Based on the Chapman–Enskog method, the mixture heat flux $\boldsymbol {q}$ is obtained as

(2.27)\begin{align} \boldsymbol{q} &= \frac{5}{2}k_BT[n_1(\boldsymbol{u}_1- \boldsymbol{u})+n_2(\boldsymbol{u}_2-\boldsymbol{u})] + \frac{5}{2}k_BTA({\boldsymbol{u}}_{1}-{\boldsymbol{u}}_{2})\nonumber\\ &\quad -\frac{5}{2}k_BT\left(\frac{n_1}{m_1} \frac{\tau_{11}\tau_{12}}{{Pr}_{12}\tau_{11}+{Pr}_{11}\tau_{12}}+ \frac{n_2}{m_2}\frac{\tau_{22}\tau_{21}}{{Pr}_{21}\tau_{22}+{Pr}_{22}\tau_{21}} - A\frac{2b_{12}}{a_{12}k_BT}\right)k_B\boldsymbol{\nabla} T. \end{align}

with

(2.28)\begin{equation} A = \gamma_{12} \left(\frac{1}{m_1}\frac{\tau_{11} \tau_{12}}{{Pr}_{12}\tau_{11}+{Pr}_{11}\tau_{12}}+ \frac{1}{m_2}\frac{\tau_{22}\tau_{21}}{{Pr}_{21}\tau_{22} +{Pr}_{22}\tau_{21}}\right)\frac{a_{12}\rho_1\rho_2}{\rho_1 \tau_{21}+\rho_2\tau_{12}}. \end{equation}

It can be seen that there are three factors contributing to the total heat flux. (i) The first term occurs since the heat flux is measured relative to the mixture flow velocity, instead of the species flow velocity, and thus, represents energy carried by the molecular flux in the presence of diffusion $\boldsymbol {u}_s-\boldsymbol {u}$. (ii) The second term with binary diffusion velocity $\boldsymbol {u}_1-\boldsymbol {u}_2$ arises as an inverse process to thermal diffusion, and has a coefficient $\frac {5}{2}k_BTA=nk_BTk_{T,12}$ based on the asymptotic analysis of the original Boltzmann equation (Chapman & Cowling 1970), leading to

(2.29)\begin{equation} \gamma_{12} =\left(\frac{1}{m_1}\frac{\tau_{11}\tau_{12}}{{Pr}_{12} \tau_{11}+{Pr}_{11}\tau_{12}}+\frac{1}{m_2}\frac{\tau_{22} \tau_{21}}{{Pr}_{21}\tau_{22}+{Pr}_{22}\tau_{21}}\right)^{{-}1} \frac{4b_{12}}{5a_{12}k_BT}. \end{equation}

(iii) The third term $-\kappa \boldsymbol {\nabla } T$ is generated by a temperature gradient, where $\kappa$ is the ordinary thermal conductivity of the mixture that is usually measured experimentally in the absence of any diffusion velocity:

(2.30)\begin{align} \kappa = \left(\frac{n_1}{m_1}\frac{\tau_{11}\tau_{12}}{{Pr}_{12} \tau_{11}+{Pr}_{11}\tau_{12}}+\frac{n_2}{m_2}\frac{\tau_{22} \tau_{21}}{{Pr}_{21}\tau_{22}+{Pr}_{22}\tau_{21}}\right) \frac{5k_B^2T}{2} - \frac{4b_{12}^2\rho_1\rho_2}{a_{12} (\rho_1\tau_{21}+\rho_2\tau_{12})T}. \end{align}

In analogy with the parameter $\phi _{sr}$ used to measure the ratio between the relaxation times of intra- and inter-species collisions, we also define

(2.31)\begin{equation} \varphi_{sr}=\frac{{Pr}_{sr}}{{Pr}_{ss}}, \end{equation}

to represent the ratio of the thermal relaxation rates. Given the species thermal conductivity $\kappa _s ={5n_sk_B^2T\tau _{ss}}/{2m_s{Pr}_{ss}}$, the mixture thermal conductivity (2.30) can be rewritten in the form

(2.32)\begin{equation} \kappa =\frac{\kappa_1}{1+\dfrac{n_2}{n_1}\phi_{12}\varphi_{12}} +\frac{\kappa_2}{1+\dfrac{n_1}{n_2}\phi_{21}\varphi_{21}} - \frac{D_{12}k_T^2n^3k_B}{n_1n_2}, \end{equation}

where the last term indicates the effect of thermal diffusion on the thermal conductance of a gas mixture. Similar to the viscosity, the values of $\varphi _{sr}$ can be determined by matching the mixture thermal conductivity measured at different proportions of gas components.

2.4. Inter-species energy relaxation

We have shown that in the continuum flow limit the temperatures of different components stay the same up to the Navier–Stokes approximation of the proposed model, when all the relaxation times are considerably smaller than the characteristic time of gas flow. Thus, the transport coefficients in the continuum flow limit are not affected by the parameters $c_{sr}$ and $d_{sr}$. However, these parameters determine the auxiliary temperatures and, hence, the energy relaxation rates between different species, which may have a significant impact in strong non-equilibrium situations. In the kinetic modelling of a single-species gas, we have found that having all the transport coefficients is not enough to exactly describe the underlying relaxation processes (Li et al. 2021, 2023; Zeng, Li & Wu 2022). The same is true for the energy relaxation during inter-species collisions of a gas mixture.

Therefore, we determine the parameters $c_{sr}$ and $d_{sr}$ for calculating auxiliary temperature by imposing that the energy exchange rates of the inter-species collision operator coincide with that of the Boltzmann collision operator,

(2.33)\begin{equation} \left\langle\frac{1}{2}m_sv^2,\frac{1}{\tau_{sr}}(g_{sr}-f_{s})\right\rangle = \left\langle \frac{1}{2}m_sv^2,Q_{sr}\right\rangle . \end{equation}

The exchange rates of the kinetic model (left-hand side) can be calculated straightforwardly, while that of the Boltzmann collision operator may only be explicitly evaluated for the Maxwellian intermolecular potential,

(2.34)\begin{align} &\frac{n_s}{\tau_{sr}}\left[\frac{3}{2}k_B(\hat{T}_{sr}-T_s)+ \frac{m_s}{2}(\hat{\boldsymbol{u}}_{sr}-\boldsymbol{u}_s)^2\right]\nonumber\\ &\quad =\lambda_{sr}\frac{m_sm_r}{(m_s+m_r)^2}n_sn_r [3k_B(T_r-T_s)+m_r(\boldsymbol{u}_r-\boldsymbol{u}_s)^2], \end{align}

where $\lambda _{sr}=\lambda _{rs}$ are constants related to the integral of collision cross-sections. By substituting the auxiliary velocity and temperature (2.13) into (2.34), we immediately get

(2.35)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle c_{sr} =\dfrac{2\lambda_{sr}m_sm_r}{(m_s+m_r)^2}(n_s\tau_{rs}+n_r\tau_{sr}),\\ \displaystyle d_{sr} =\dfrac{\lambda_{sr}m_sm_r}{3k_B(m_s+m_r)^2}[{\lambda_{sr} (n_r\rho_r\tau_{sr}^2-n_s\rho_s\tau_{rs}^2)} -{2(\rho_r\tau_{sr}-\rho_s\tau_{rs})}]. \end{array}\right\} \end{equation}

Note that the same form of (2.34) can be obtained for non-Maxwellian intermolecular potentials with non-constant $\lambda _{sr}$, when it is subject to the restriction that the distribution functions are Maxwellian at different temperatures but with small diffusion velocities (Morse 1964). Therefore, we calculate the parameters $c_{sr}$ and $d_{sr}$ using (2.35) for any type of intermolecular potential.

Furthermore, variable $\lambda _{sr}$ can be approximated by matching the momentum exchange rates of the inter-species collision operator with that of the Boltzmann collision operator for the Maxwellian intermolecular potential,

(2.36)\begin{equation} \left\langle m_sv,\frac{1}{\tau_{sr}}(g_{sr}-f_{s})\right\rangle = \langle m_sv,Q_{sr}\rangle , \end{equation}

which leads to

(2.37)\begin{equation} \frac{n_sm_s}{\tau_{sr}}(\hat{\boldsymbol{u}}_{sr}-{\boldsymbol{u}}_{s}) =\lambda_{sr}\frac{m_sm_r}{m_s+m_r}n_sn_r({\boldsymbol{u}}_{r}- {\boldsymbol{u}}_{s}). \end{equation}

Therefore, $\lambda _{sr}$ can be evaluated by parameter $a_{sr}$ as

(2.38)\begin{equation} \lambda_{sr} = \frac{a_{sr}(m_s+m_r)}{\rho_s\tau_{rs}+\rho_r\tau_{sr}}, \end{equation}

which is thus determined by the binary diffusion coefficient as given by (2.25).

2.5. Indifferentiability principle

The indifferentiability principle states that, when all the molecules are mechanically identical (e.g. they have the same mass and scattering cross-section), the model equation reduces to a single one by adding the distribution functions (Garzó, Santos & Brey 1989), although in reality one does not have this situation unless one uses two Boltzmann equations to simulate single-species flow. This property holds for the Boltzmann equation because of the bilinearity of its operators. It is however non-trivial for a model equation to inherit since the operators constructed are usually highly nonlinear. Historically, several kinetic models using a single-relaxation collision operator have been proven to fulfil this principle. It should be noted that some of them require a condition that the diffusion velocities vanish for the indifferentiable molecules, when the models contain parameters to recover the Fick law (Brull, Pavan & Schneider 2012; Todorova & Steijl 2019).

Therefore, we also adopt the assumption that $\boldsymbol {u}_s=\boldsymbol {u}$ for all the indifferentiable species $s$, to demonstrate that our multi-relaxation model with the linearised collision operator complies with the indifferentiability principle. Consider a gas mixture system that slightly deviates from an equilibrium state with flow velocity $\boldsymbol {u}$, temperature $T$ and number density $n_s$ of each component, the reference distribution in the linearised collision operator is given as

(2.39) \begin{align} g_{sr}^{linear} &= f_{s}^{eq}\left[1 + \frac{m_s (\hat{\boldsymbol{u}}_{sr}-\boldsymbol{u}) \boldsymbol{\cdot}\boldsymbol{c}}{k_BT} + \frac{\hat T_{sr}-T}{T} \left(\frac{m_sc^2}{2k_B{T}}-\frac{3}{2}\right)\right.\nonumber\\ &\quad +\left. \frac{2m_s\hat{\boldsymbol{q}}_{sr}\boldsymbol{\cdot}\boldsymbol{c}}{5{n}_{s} k_B^2{T}^2}\left(\frac{m_sc^2}{2k_B{T}}-\frac{5}{2}\right)\right], \end{align}

with

(2.40)\begin{equation} f_s^{eq} = n_s\left(\frac{m_s}{2{\rm \pi} k_BT}\right)^{3/2} \exp{\left(-\frac{m_s(\boldsymbol{v}-{\boldsymbol{u}})^2}{2k_BT}\right)}. \end{equation}

For the indifferentiable molecules, we have (i) $m_s=m$ and $n_r\tau _{sr}$ is constant for any $s$ and $r$ due to the identical mass and scattering cross-sections, respectively; (ii) $\hat {\boldsymbol {u}}_{sr}=\boldsymbol {u}$ based on (2.13) with the thermal-diffusion coefficient vanishing; (iii) $\hat T_{sr}+\hat T_{rs}=T_s+T_r$ and $nT=\sum n_sT_s$ from total energy conservation (2.11) and calculation of the mixture temperature (2.6), respectively; (iv) $\hat {\boldsymbol {q}}_{sr}=(1-{Pr})\boldsymbol {q}_{s}$ from (2.15). Then, the sum of the kinetic equations over all species yields

(2.41)\begin{align} {\mathcal{D}\left(\sum_{s=1}^Nf_{s}\right)} &= \sum_{s=1}^N{\sum_{r=1}^N\frac{1}{\tau_{sr}} (g_{sr}^{linear}-f_s)} \nonumber\\ &= \frac{1}{\tau}\left\{\sum_sf_s^{eq}\left[1+\frac{2m(1 -{Pr}){\left(\sum_s\boldsymbol{q}_{s}\right)}\boldsymbol{\cdot} \boldsymbol{c}}{5{n}k_B^2{T}^2}\left(\frac{mc^2}{2k_B{T}}- \frac{5}{2}\right)\right]-\sum_{s=1}^Nf_{s}\right\}, \end{align}

where $\tau =n_r\tau _{sr}/n$ is the overall relaxation time. Clearly, the kinetic model equation reduces to the Shakhov model of single-species monatomic gas with the distribution function $f=\sum f_{s}$.

2.6. Dimensionless forms

Let $L_0, T_0, n_0, m_0$ be the reference length, temperature, number density and mass, respectively, then the most probable speed is

(2.42)\begin{equation} v_m=\sqrt{\frac{2k_BT_0}{m_0}}, \end{equation}

and the reference pressure is $p_0=n_0k_BT_0$. The dimensionless variables are introduced as

(2.43)\begin{equation} \left.\begin{array}{@{}c@{}} \tilde{\boldsymbol{x}}=\boldsymbol{x}/L_0, \quad \tilde{n}=n/n_0, \quad \tilde{m}=m/m_0, \quad\tilde{T}=T/T_0, \\ \tilde{\boldsymbol{v}}=\boldsymbol{v}/v_m, \quad \tilde{\boldsymbol{c}}=\boldsymbol{c}/v_m, \quad \tilde{t}=v_mt/L_0, \quad \tilde{\tau}=v_m\tau/L_0, \\ \tilde{p}=p/p_0, \quad \tilde{\boldsymbol{q}}=\boldsymbol{q}/(p_0v_m), \quad \tilde{f}_s=v_m^{3}f_s/n_0. \end{array}\right\} \end{equation}

The Knudsen numbers ${Kn}_{s}$ of each species $s$ is defined as

(2.44)\begin{equation} {Kn}_{s}=\frac{\mu_s(T_0)}{n_0L_0}\sqrt{\frac{\rm \pi}{2m_sk_BT_0}}. \end{equation}

It is noted that the species-specific Knudsen numbers are correlated as ${Kn}_{r}={Kn}_{s}\beta ^{\mu }_{rs}/\sqrt {\beta ^m_{rs}}$, with $\beta ^{\mu }_{rs}=\mu _r(T_0)/\mu _s(T_0)$ being the viscosity ratio at the reference temperature and $\beta ^m_{rs}=m_r/m_s$ the mass ratio. Therefore, the dimensionless relaxation times can be written in terms of the Knudsen numbers,

(2.45a,b)\begin{equation} \tilde{\tau}_{ss} = 2{Kn}_{s}\sqrt{\frac{\tilde{m}_s}{\rm \pi}} \frac{\tilde{T}_s^{\omega_s-1}}{\tilde{n}_s}, \quad \tilde{\tau}_{sr} = \tilde{\tau}_{ss}\phi_{sr}^{{-}1} \frac{\tilde{n}_s}{\tilde{n}_r}, \end{equation}

where $\omega _s$ is the viscosity index of species $s$ in

(2.46)\begin{equation} \mu_s(T)=\mu_s(T_0)\left(\frac{T}{T_0}\right)^\omega_s. \end{equation}

Then, the kinetic model equations are non-dimensionalised as

(2.47)\begin{equation} \frac{\partial{\tilde{f}_{s}}}{\partial{\tilde{t}}}+\tilde{\boldsymbol{v}} \boldsymbol{\cdot} \frac{\partial{\tilde{f}_{s}}}{\partial{\tilde{\boldsymbol{x}}}}+ \frac{\tilde{\boldsymbol{F}}_s}{\tilde{m}_s}\boldsymbol{\cdot} \frac{\partial{\tilde{f}_{s}}}{\partial{\tilde{\boldsymbol{v}}}} =\sqrt{\frac{\rm \pi}{\tilde{m}_s}}\frac{\tilde{T}_s^{1-\omega_s}}{2{Kn}_{s}} \left[\tilde{n}_s(\tilde{g}_{ss}-\tilde{f}_s)+\sum_{r \neq s}\tilde{n}_r\phi_{sr}(\tilde{g}_{sr}-\tilde{f}_s)\right], \end{equation}

with the dimensionless reference velocity distribution function

(2.48) \begin{align} \tilde{g}_{sr}&=\widetilde{{n}}_{sr}\left(\frac{\tilde{m}_s}{{\rm \pi} \tilde{T}_{s}}\right)^{3/2}\exp\left(-\frac{\tilde{m}_s(\tilde{\boldsymbol{v}}- \tilde{\hat{\boldsymbol{u}}}_{sr})^2}{\tilde{T}_{s}}\right) \times \left[1+\frac{\tilde{\hat{T}}_{sr}-\widetilde{T_s}}{\widetilde{T_s}} \left(\frac{\tilde{m}_s(\tilde{\boldsymbol{v}}- \tilde{\hat{\boldsymbol{u}}}_{sr})^2}{\tilde{T}_{s}}- \frac{3}{2}\right)\right.\nonumber\\ &\quad +\left.\vphantom{1+\frac{\tilde{\hat{T}}_{sr}-\widetilde{T_s}}{\widetilde{T_s}} \left(\frac{\tilde{m}_s(\tilde{\boldsymbol{v}}- \tilde{\hat{\boldsymbol{u}}}_{sr})^2}{\tilde{T}_{s}}- \frac{3}{2}\right)}\frac{4\tilde{m}_s\tilde{\hat{\boldsymbol{q}}}_{sr} \boldsymbol{\cdot}(\tilde{\boldsymbol{v}}-\tilde{\hat{\boldsymbol{u}}}_{sr})}{5 \tilde{\hat{n}}_{sr}\tilde{T}_{s}^2}\left(\frac{\tilde{m}_s (\tilde{\boldsymbol{v}}-\tilde{\hat{\boldsymbol{u}}}_{sr})^2}{\tilde{T}_{s}}- \frac{5}{2}\right)\right]. \end{align}

It clearly shows that the strengths of intra- and inter-species collisions are indicated by the magnitudes of $\tilde {n}_s$ and $\tilde {n}_r\phi _{sr}$, respectively. Therefore, in a gas mixture with a large disparity in concentration or mass, the intra-species collisions become dominant for the component $s$ with a major proportion of number density ($\tilde {n}_s \gg \tilde {n}_r$) or significantly heavier mass ($\phi _{sr} \ll 1$ when $m_r\ll m_s$).

3. Determination of parameters

Given the relationship between the model parameters and transport coefficients of the mixtures as (2.22), (2.25), (2.26a,b) and (2.32), the adjustable parameters can be uniquely determined by the experimentally measured properties of gas mixtures directly, without the knowledge of any intermolecular potentials that is also constructed to approximate the real gas properties.

Nevertheless, to validate our kinetic model in various rarefied flow problems of gas mixtures with a wide range of mass ratios and different types of molecular interactions, we compare the solutions of our model with DSMC results for virtual gases with well-defined intermolecular potentials. The transport properties of a gas in DSMC simulations are the result of the transfer of mass, momentum and energy through particle movement and collision dynamics of corresponding collision models. Given the information of any intermolecular potentials, the DSMC method with the variable-soft-sphere (VSS) collision model (Koura & Matsumoto 1991) captures the viscosity and diffusion cross-sections simultaneously by adjusting two parameters, the viscosity index $\omega$ and the angular scattering parameter $\alpha$, and thus, provides reference solutions consistent with those from the Boltzmann equation for monatomic gases. For example, the hard-sphere molecules have $\omega =0.5$ and $\alpha =1$, and the Maxwell molecules have $\omega =1.0$ and $\alpha =2.14$. Therefore, to make a consistent comparison between the results from the proposed kinetic model and the DSMC method, we determine the parameters in the model equation by matching the transport properties of a gas mixture from the VSS model applied in the DSMC simulations.

The transport coefficients for both simple gases and gas mixtures can be approximated by the Chapman–Enskog solutions of the Boltzmann equation. It should be noted that, the first approximations of the transport coefficients for non-Maxwell gases may have differences compared with the exact values. Based on the Chapman–Enskog solutions to higher order in Sonine polynomials (Tipton, Tompson & Loyalka 2009a,b), for a gas mixture consisting of hard-sphere molecules, the discrepancies in mixture viscosity and diffusion coefficient are usually limited, namely, less than 3 % for a wide range of molecular mass, sizes and mole fractions. On the other hand, the discrepancies can be non-negligible in thermal conductivity, that is, around 10 % in some mixtures, and even significant for thermal-diffusion coefficients (higher than 20 %). In a binary gas mixture of monatomic molecules, the first approximations of the transport coefficients, denoted by $[{\cdot }]_1$, are given by (Chapman & Cowling 1970)

(3.1)\begin{align} \left.\begin{array}{@{}c@{}} \displaystyle {[\mu_s]}_1 = \dfrac{5k_BT}{8\varOmega_s^{(2)}(2)}, \\ \displaystyle {[\kappa_s]}_1 = \dfrac{75k_B^2T}{32m_s\varOmega_s^{(2)}(2)}, \\ \displaystyle {[D_{12}]}_1 = \dfrac{3E}{2n(m_1+m_2)}, \\ \displaystyle {[k_T]}_1 = 5C\dfrac{n_1n_2}{n}\dfrac{n_1S_1-n_2S_2}{n_1^2Q_1+n_2^2Q_2+n_1n_2Q_{12}}, \\ \displaystyle {[\mu]}_1 = \dfrac{n_1^2R_1+n_2^2R_2+n_1n_2R'_{12}}{n_1^2R_1/[\mu_1]_1+n_2^2R_2/ {[\mu_2]}_1+n_1n_2R_{12}}, \\ \displaystyle {[\kappa]}_1 = \dfrac{n_1^2Q_1[\kappa_1]_1+n_2^2Q_2[\kappa_2]_1 +n_1n_2Q'_{12}}{n_1^2Q_1+n_2^2Q_2+n_1n_2Q_{12}}, \end{array}\right\} \end{align}

where $\varOmega _{s}$ is the integral of intra-species collisions, the variables $S, Q, R$ can be expressed in terms of species viscosities $[\mu _s]_1$ and mass fraction $M_s=m_s/(m_1+m_2)$,

(3.2)\begin{gather} \left.\begin{array}{@{}c@{}} \displaystyle S_1 = \dfrac{M_1^2E}{[\mu_1]_1}-M_2(3(M_2-M_1)+4M_1A), \\ \displaystyle S_2 = \dfrac{M_2^2E}{[\mu_2]_1}-M_1(3(M_1-M_2)+4M_2A), \\ \displaystyle Q_1 = \dfrac{M_1E}{[\mu_1]_1}(6M_2^2+5M_1^2-4M_1^2B+8M_1M_2A), \\ \displaystyle Q_2 = \dfrac{M_2E}{[\mu_2]_1}(6M_1^2+5M_2^2-4M_2^2B+8M_2M_1A), \\ \displaystyle Q_{12} = 3(M_1-M_2)^2(5-4B)+4M_1M_2A(11-4B)+\dfrac{2M_1M_2E^2}{[\mu_1]_1[\mu_2]_1}, \\ \displaystyle Q'_{12} = \dfrac{15k_BE}{2(m_1+m_2)}\left(\dfrac{M_1E}{[\mu_1]_1} +\dfrac{M_2E}{[\mu_2]_1}+(11-4B-8A)M_1M_2\right), \end{array}\right\} \end{gather}

(3.3)\begin{gather} \left.\begin{array}{@{}c@{}} \displaystyle R_1 = \dfrac{2}{3}+\dfrac{M_1A}{M_2}, \quad R_2 = \dfrac{2}{3}+\dfrac{M_2A}{M_1}, \\ \displaystyle R_{12} = \dfrac{E}{2[\mu_1]_1[\mu_2]_1}+\dfrac{4A}{3EM_1M_2}, \\ \displaystyle R'_{12} = \dfrac{E}{2[\mu_1]_1}+\dfrac{E}{2[\mu_2]_1} +2\left(\dfrac{2}{3}-A\right). \end{array}\right\} \end{gather}

Here, $A, B, C, E$ are functions of integrals $\varOmega _{12}$ of inter-species collisions,

(3.4)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle A =\dfrac{\varOmega_{12}^{(2)}(2)}{5\varOmega_{12}^{(1)}(1)}, \quad B = \dfrac{5\varOmega_{12}^{(1)}(2)-\varOmega_{12}^{(1)}(3)}{5 \varOmega_{12}^{(1)}(1)}, \\ \displaystyle C = \dfrac{2\varOmega_{12}^{(1)}(2)}{5\varOmega_{12}^{(1)}(1)}-1, \quad E = \dfrac{k_BT(m_1+m_2)^2}{8m_1m_2\varOmega_{12}^{(1)}(1)}. \end{array}\right\} \end{equation}

Based on the VSS collision model, the $\varOmega$ integrals can be calculated as given by Stephani, Goldstein & Varghese (2012),

(3.5)\begin{equation} \left.\begin{array}{@{}c@{}} \displaystyle \varOmega_{s}^{(2)}(2) = \dfrac{4\alpha_s}{(\alpha_s+1)(\alpha_s+2)}\dfrac{\rm \pi}{2} \left(\dfrac{k_BT}{{\rm \pi} m_{s}}\right)^{1/2}\left(\dfrac{5}{2}- \omega_{s}\right)\left(\dfrac{7}{2}-\omega_{s}\right) \left(\dfrac{T_{0}}{T}\right)^{\omega_{s}-1/2}d_{s,ref}^2, \\ \displaystyle \varOmega_{12}^{(1)}(1) = \dfrac{2}{\alpha_{12}+1} \dfrac{\rm \pi}{2}\left(\dfrac{k_BT}{2{\rm \pi} m_{12}}\right)^{1/2} \left(\dfrac{5}{2}-\omega_{12}\right)\left(\dfrac{T_{0}}{T} \right)^{\omega_{12}-1/2}d_{12,ref}^2, \\ \displaystyle \varOmega_{12}^{(1)}(2) = \varOmega_{12}^{(1)}(1)\left(\dfrac{7}{2}-\omega_{12}\right), \\ \displaystyle \varOmega_{12}^{(1)}(3) = \varOmega_{12}^{(1)}(1) \left(\dfrac{7}{2}-\omega_{12}\right)\left(\dfrac{9}{2}-\omega_{12}\right), \\ \displaystyle \varOmega_{12}^{(2)}(2) = \dfrac{4\alpha_{12}}{(\alpha_{12}+1)(\alpha_{12}+2)} \dfrac{\rm \pi}{2}\left(\dfrac{k_BT}{2{\rm \pi} m_{12}}\right)^{1/2} \left(\dfrac{5}{2}-\omega_{12}\right)\left(\dfrac{7}{2}-\omega_{12}\right) \left(\dfrac{T_{0}}{T}\right)^{\omega_{12}-1/2}d_{12,ref}^2, \end{array}\right\} \end{equation}

where $m_{12}=m_1m_2/(m_1+m_2)$ is the reduced mass; $d_{ref}$ is the reference collision diameter, $\omega$ is the viscosity index, $\alpha$ is the angular scattering parameter in the VSS model, with subscripts $s$ and $12$ indicating the values for intra- and inter-species collisions, respectively.

Therefore, given the collision parameters in the VSS model, the first approximations of the transport properties of a gas mixture can be explicitly evaluated. Then the parameters $\phi _{sr}$ measuring the ratio between relaxation times of intra- and inter-species collisions can be obtained by following two straightforward steps.

1. (i) First, the mixture viscosities are calculated as a function of mole concentrations based on (3.1)–(3.5).

2. (ii) Second, we find the optimized values of $\phi _{12}$ and $\phi _{21}$ in expression (2.22) by the least square method to fit the mixture viscosities given by the first step.

Similarly, the ratio of thermal relaxation rates $\varphi _{sr}$ can be fitted by matching the mixture thermal conductivity using (2.32), once $\phi _{sr}$ is determined. Also, the parameters $a_{sr}$ and $b_{sr}$ are calculated based on the diffusion coefficients (2.25) and thermal-diffusion ratio (2.26a,b), respectively.

The value of $\phi _{sr}=n_s\tau _{ss}/n_r\tau _{sr}$ is not only essential to recover the shear viscosity of the gas mixture, but also signifies the relaxation time scales of different types of binary collisions. Then, to explore the dependence of the relaxation time ratio on mass disparity, we obtain $\phi _{sr}$ across a wide range of mass ratios $m_2/m_1=1 \sim 10^4$ for Maxwell gas mixtures with a fixed reference diameter ratio $d_2/d_1=1$ and hard-sphere gas mixtures with $d_2/d_1=2$, as shown in figure 1. It is found that for lighter species (typically smaller in diameter), $\phi _{12}=n_1\tau _{11}/n_2\tau _{12}$ is greater than 1 and also increases gradually with mass ratio (figure 1a), which implies an important role of inter-species collisions with heavy gas molecules. On the other hand, $\phi _{21}=n_2\tau _{22}/n_1\tau _{21}$ rapidly decreases as the mass disparity grows, and roughly scales as $({m_2}/{m_1})^{-0.59}$ when the mass ratio exceeds 20 (figure 1b). Consequently, in a gas mixture with disparate mass, the inter-species collision term in our kinetic model will have a negligible impact on the dynamics of the heavier species unless its mole fraction is significantly small. Note that, knowing the relaxation time scale of each collision term is also crucial for developing a multi-temperature hydrodynamic equation for mixtures with disparate mass, since the corresponding constitutive relations depend on the orders of magnitude of the inter-species relaxation rates (Bisi, Groppi & Martalò 2021).

Figure 1. The relaxation time ratios (a) $\phi _{12}$ and (b) $\phi _{21}$ fitted by the first approximations of the mixture viscosities with the mass ratio $m_2/m_1$ varies from 1 to $10^4$, for Maxwell gas mixtures with a fixed reference diameter ratio $d_2/d_1=1$ and hard-sphere gas mixtures with $d_2/d_1=2$.

In the present work, three types of binary gas mixtures are considered in the following simulations: mixtures 1 and 2 consist of Maxwell gas molecules possessing identical reference diameters and mass ratios of 10 and 1000, respectively; mixture 3 composites hard-sphere molecules characterized by a mass ratio of 100 and a diameter ratio of 2. All the parameters associated with inter-species collisions in the VSS model (namely, $\omega _{12}$, $\alpha _{12}$ and $d_{12}$) are determined simply through the arithmetic averaging of the corresponding parameters of the individual components in each mixture. The fitting parameters $\phi _{sr}$ and $\varphi _{sr}$ for relaxation rate ratios are given in table 1. Due to the discrepancy by the first approximation of mixture thermal conductivity for hard-sphere molecules as mentioned above, we determine the values of $\varphi _{sr}$ for mixture 3 based on the thermal conductivities calculated from the DSMC method directly, to make a more fair comparison for model validation. It is noteworthy that the fitted values of $\phi _{sr}$ and $\varphi _{sr}$ for the considered mixtures are independent of temperature, although this is generally not the case for arbitrary mixtures containing gas molecules with different interaction potentials.

Table 1. The constituents of the three binary mixtures considered in the present work, and the corresponding parameters $\phi _{sr}$ and $\varphi _{sr}$ in the kinetic model fitted by matching the mixture viscosity and thermal conductivity from the VSS collision model, respectively.

4. Numerical results of one-dimensional problems

In this section the accuracy of our kinetic model is assessed by the DSMC method in a one-dimensional normal shock wave, Fourier flow and Couette flow of the binary gas mixtures listed in table 1. We compare not only the average properties of the mixture but also those of the individual components, which is crucial for accurately describing mixture flows, as the different species in the mixture can vary significantly in concentration, velocity and temperature in non-equilibrium flows.

The DSMC simulations are conducted using the open-source code SPARTA (Plimpton et al. 2019). We use uniform spatial cells in the one-dimensional simulations, with each cell size $L_{{cell}}^{{DSMC}}$ set to approximately 1/10 of the mean free path of the species with a larger collision diameter. Note that for the Fourier and Couette flows with a Knudsen number exceeding 0.1, 100 cells are used. The time step applied in the simulations is $L_{{cell}}^{{DSMC}}/5v_m^{{light}}$, where $v_m^{{light}}$ is the most probable speed of the lighter molecules. On average, we use 100 simulation particles per cell on the upstream side in normal shock wave problems, and 200 simulation particles per cell in Fourier and Couette flow simulations.

On the other hand, to reduce the computational cost of solving kinetic model equations, the velocity distribution functions are dimensionally reduced to be quasi-one-dimensional in velocity space for the normal shock wave and Fourier flow, by introducing functions $f_{s,x1}$ and $f_{s,x2}$:

(4.1a,b)\begin{equation} f_{s,x1} = \int_{\mathbb{R}^2}{f_s}\,\mathrm{d}{v_y}\,\mathrm{d}{v_z}, \quad f_{s,x2} =\int_{\mathbb{R}^2}{(v_y^2+v_z^2)f_s}\,\mathrm{d}{v_y} \,\mathrm{d}{v_z}. \end{equation}

Similarly, $f_{s}$ can be reduced to be quasi-two-dimensional in velocity space for the Couette flow, by introducing functions $f_{s,xy1}$ and $f_{s,xy2}$:

(4.2a,b)\begin{equation} f_{s,xy1} =\int_{\mathbb{R}^1}{f_s}\,\mathrm{d}{v_z}, \quad f_{s,xy2} =\int_{\mathbb{R}^1}{v_z^2f_s}\,\mathrm{d}{v_z}. \end{equation}

The macroscopic variables are then calculated by taking moments of the reduced distribution functions. The kinetic model equations are solved by the discretized velocity method (Zeng et al. 2023b; Liu et al. 2024; Zeng, Li & Wu 2024). A finite difference method with a second-order upwind scheme is adopted in the numerical implementation, where we use the forward Euler scheme for the time derivative and apply implicit treatment to the convection term and distribution functions in the collision term. Non-uniform spatial grids are employed, with mesh refinement applied around the centre of the shock waves and in the vicinity of the solid walls in Fourier and Couette flows. Since the thermal velocity of the lighter gas can be much higher than that of the heavier one by the square root of the mass ratio, the individual velocity space of each component is used with different truncations and discretization. Note that the model equation for each species is concerned only with its own velocity distribution function, thus eliminating the need for interpolation within the velocity space in the numerical scheme implementation.

4.1. Normal shock wave

We investigate the structure of normal shock waves under various free-stream conditions. The conditions are defined by the species mole fraction $\chi _s$ in the upstream, as well as the Mach number (${Ma}$) calculated based on the speed of sound ${v}_{mix}=\sqrt {5k_BT_u/3m_{mix}}$ in the upstream, where $T_u$ is the upstream temperature and $m_{mix}=\sum {m_s\chi _s}$ is the averaged mass of the mixture. Given the upstream conditions, the macroscopic quantities at the downstream end are determined by the classical Rankine–Hugoniot relation. The mass of the lighter species (denoted as species 1), the mixture number density $n_u$ and the temperature $T_u$ of the upstream flow are taken as reference values, namely, $m_0=m_1$, $n_0=n_u$, $T_0=T_u$. The characteristic length $L_0$ is set equal to the mean free path of the lighter species in the upstream flow, thus leading to ${Kn}_1=1$.

The numerical investigations of normal shock wave structures in gas mixtures based on the Boltzmann equation and the DSMC method have been conducted in the literature, taking into account various intermolecular interactions, including hard-sphere and ab initio potentials. We validate the proposed kinetic model by comparing its results with existing published data (Kosuge, Aoki & Takata 2001; Raines 2002; Sharipov & Dias 2018), as given in Appendix B. This comparison not only demonstrates the accuracy of our model but also highlights its applicability in accommodating diverse intermolecular potentials. Furthermore, we apply both the proposed kinetic model and the model of Kosuge (2009) to simulate a weak and a strong shock wave with ${Ma}=2$ and 10, respectively. The comparative analysis with DSMC simulations indicates that both models can accurately capture the structure of the weak shock wave, while our model outperforms the Kosuge model in modelling strong non-equilibrium shock waves, as detailed in Appendix C.

We then simulate the binary mixtures with large mass ratios (up to 1000) to determine the limits of our kinetic model. The simulation domain $[-L_x,L_x]$ is selected to ensure that the upstream and downstream are in equilibrium states, and then the specific values of $L_x$ are chosen as $40L_0,500L_0,30L_0$ for mixtures 1, 2 and 3, respectively, due to the significant differences in the properties of these gas mixtures and the upstream conditions. Numerical results of both the kinetic model and the DSMC method are compared in figures 2–4 for mixtures 1, 2 and 3, respectively. We present the normalized values for number density, flow velocity, temperature and dimensionless heat flux for each species and also the mixture, under the specified conditions with different concentrations $\chi _1=0.1,0.5,0.9$ and ${Ma}=3,4,5$. For mixtures 1 and 2 consisting of Maxwell gases, excellent agreements between the results of our kinetic model and the DSMC method are achieved (figures 2 and 3), even when the mass ratio is as high as 1000. For mixture 3 consisting of hard-sphere molecules, the number density and velocity of the lighter species predicted by the kinetic model deviate from DSMC results, when the mole fraction of the lighter species is small ($\chi _1=0.1$; see figures 4a and 4d). However, the average properties of the mixture remain highly accurate, since the lighter species, due to its low concentration, has a minimal impact on the overall behaviour of the mixture. Generally, three reasons may contribute to the possible deviation predicted by the kinetic model for mixtures consisting of non-Maxwell gases: (i) the parameters $c_{sr}$ and $d_{sr}$ accounting for the energy relaxation are derived approximately for non-Maxwell molecules; (ii) the term with parameter $b_{sr}$ is designed to capture the correct thermal-diffusion phenomena in the continuum limit, and thus, may have a discrepancy in strong non-equilibrium cases; (iii) the velocity-dependent collision frequency for non-Maxwell molecules is not recovered in the BGK-type operators (Yuan & Wu 2022). Nevertheless, as shown in figure 4, the kinetic model gives good overall agreement with DSMC simulations for hard-sphere gas mixtures. Therefore, the agreement suggests that the kinetic model can predict accurate results for realistic molecular models, whose behaviour usually lies between that of hard-sphere and Maxwell molecules.

Figure 2. Comparisons of the normalized (a–c) number density, (d–f) flow velocity, (g–i) temperature and (j–l) dimensionless heat flux of the gas mixture between the kinetic model (lines) and the DSMC method (symbols) for the normal shock wave at ${Ma}=5$. The binary mixture consists of Maxwell molecules with a mass ratio $m_2/m_1=10$, a diameter ratio $d_2/d_1=1$ and the mole fraction of light species $\chi _1=0.1,0.5,0.9$.

Figure 3. Comparisons of the normalized (a–c) number density, (d–f) flow velocity, (g–i) temperature and (j–l) dimensionless heat flux of the gas mixture between the kinetic model (lines) and the DSMC method (symbols) for the normal shock wave at ${Ma}=3$. The binary mixture consists of Maxwell molecules with a mass ratio $m_2/m_1=1000$, a diameter ratio $d_2/d_1=1$ and the mole fraction of light species $\chi _1=0.1,0.5,0.9$.

Figure 4. Comparisons of the normalized (a–c) number density, (d–f) flow velocity, (g–i) temperature and (j–l) dimensionless heat flux of the gas mixture between the kinetic model (lines) and the DSMC method (symbols) for the normal shock wave at ${Ma}=4$. The binary mixture consists of hard-sphere molecules with a mass ratio $m_2/m_1=100$, a diameter ratio $d_2/d_1=2$ and the mole fraction of light species $\chi _1=0.1,0.5,0.9$.

Although the shock waves in mixtures with large mass ratios (e.g. 32.8 for Helium–Xenon) have been studied in the literature, and have shown unique characteristics absent in those composed of similar gas molecules. The shock structures can be significantly altered by substantial mass disparity. By comparing the shock structures in a wide range of mass ratios and species concentrations, the following features can be observed and correctly captured by our kinetic model.

1. (i) The shock wave thickness of a gas mixture is markedly thicker than that of a pure gas, especially when the mass ratio is large, as the mixture viscosity and diffusivity become stronger and the relaxation between components gets slower. For example, as shown in figure 3, the Maxwell gas mixture with a mass ratio of 1000 and diameter ratio of 1 form a significantly large transition zone from the upstream to the downstream, which spans several hundreds of the molecular mean free path.

2. (ii) In mixture 1 with a mass ratio of 10 (moderate mass difference), a pronounced temperature overshoot of the heavy species (higher than the downstream temperature) is observed when the heavier gas has only a small proportion (figure 2i), which has been shown in the literature from kinetic modelling (Bird 1968; Kosuge et al. 2001; Sharipov & Dias 2018) and hydrodynamic equations (Schmidt, Seiler & Wörner 1984). However, when the mass ratio increases to 1000, the temperature overshoot gradually vanishes when $\chi _1=0.9$. As illustrated in figure 3, a comparable or even larger proportion of the heavier gas is required for the temperature overshoot to occur.

3. (iii) Flow velocity undershoot of the lighter gas happens in a mixture with a large mass ratio (100 and 1000) and a small proportion of lighter species ($\chi _1=0.1$) (figures 3d and 4d), which confirms the phenomena predicted by the multi-temperature hydrodynamic equations on weak shock (Goldman & Sirovich 1969). In these situations, the lighter gas decelerates to a velocity lower than that of the downstream flow, and thus, is compressed to a density above the downstream one.

4. (iv) A two-stage shock structure with distinct gradients of the mixture properties can be observed. In mixture 1 with a mass ratio of 10 and a large concentration of light species ($\chi _1=0.9$), as shown in the third column of figure 2, the properties of the shock exhibit a very steep change on the upstream side, while followed by a sudden change in the form of a long tail downstream. This phenomenon happened in a mixture with a moderate mass difference in experiments of Helium–Xenon (Gmurczyk, Tarczynski & Walenta 1979) and found by hydrodynamic equations (Ruyev, Fedorov & Fomin 2002). However, when the mass ratio becomes significantly large in mixtures 2 and 3 (figures 3 and 4), this phenomenon disappears. On the contrary, the two-stage shock structure with an opposite trend occurs when there is only a small proportion of light species in the mixture ($\chi _1=0.1$). This structure consists of a smooth change in gas properties on the upstream side followed by a sudden and dramatic change on the downstream side.

4.2. Planar Fourier flow

A stationary rarefied mixture confined between two infinite parallel plates located at $x=0$ and $x=L_0$ is considered. The surfaces of the plates have different temperatures $T_w=T_0 \pm \Delta T$, and reflect the gas molecules in a fully diffuse way. The Knudsen number is determined in terms of the averaged number density of the mixture $n_0$, the temperature $T_0$ and the distance between the two plates $L_0$. A variety of cases were considered for the mixtures with different mole fractions $\chi _s$, temperature difference $\Delta T$ and Knudsen numbers, while a selection of the representative results is shown in figures 5 and 6. It can be seen that the solutions given by the kinetic model agree well with the DSMC results.

Figure 5. Comparisons of the normalized number density (first column), dimensionless temperature (second column) and heat flux (third column) of the gas mixture between the kinetic model (lines) and the DSMC method (symbols) for the planar Fourier flow with $\Delta T=0.2T_0$. The binary mixture consists of Maxwell molecules with a diameter ratio $d_2/d_1=1$ and the mole fraction of light species $\chi _1=0.5$. Each row belongs to a specific flow condition: (a–c) mixture 1 ($m_2/m_1=10$), ${Kn}_1=0.1$; (d–f) mixture 1, ${Kn}_1=1$; (g–i) mixture 2 ($m_2/m_1=1000$), ${Kn}_1=0.1$; (j–l) mixture 2, ${Kn}_1=1$.

Figure 6. Comparisons of the normalized number density (a,d), dimensionless temperature (b,e) and heat flux (c, f) of the gas mixture between the kinetic model (lines) and the DSMC method (symbols) for the Fourier flow, when ${Kn}_1=0.1$ and $\Delta T=0.8T_0$. The binary mixture consists of hard-sphere molecules with a mass ratio $m_2/m_1=100$, a diameter ratio $d_2/d_1=2$ and the mole fraction of light species $\chi _1=0.1$ (a–c) and 0.9 (d–f). The black lines are the kinetic model solutions without thermal-diffusion effects ($b_{12}=0$).

The temperatures of the components in mixture 1 ($m_2/m_1=10$) stay close when ${Kn}$ is up to 1 (figures 5a–5f), while pronounced temperature separation and concentration variation can be observed for mixture 2 ($m_2/m_1=1000$) when ${Kn}=0.1$ (figures 5g–5l). It is noteworthy that all these Maxwell molecules under consideration have the same size of mean free path, due to their identical diameter and interaction potential. However, the inter-species relaxation is much slower for the gases with disparate mass, although the spatial Knudsen number is the same. In other words, the mixtures with larger mass ratios may exhibit significant non-equilibrium phenomena even at small ${Kn}$, and hence, shrink the applicable range of the hydrodynamic description of the mixtures. The same observation has been found in a previous work solving the linearised Boltzmann equation and McCormack model (Ho et al. 2016).

Unlike Maxwell gas mixtures, where the thermal-diffusion effect is absent, mixtures of hard-sphere molecules have a significant thermal-diffusion effect. As shown in figure 6 (${Kn}_1=0.1$ and $\Delta T=0.8T_0$), although the temperature of the two species remains the same, the concentration ratio between components varies across the domain due to the temperature gradient, and the heavy gas molecules tend to concentrate in the cold region. Meanwhile, the highly nonlinear feature arising from the pronounced temperature difference between the two plates is accurately captured. We also solve the kinetic model without the thermal-diffusion effect by setting parameter $b_{12}=0$ (other transport properties and relaxation rates remain unchanged), corresponding results are shown in figure 6 by black lines. Clearly, the species separation phenomenon cannot be reproduced, thus leading to an incorrect prediction of concentration and, hence, the heat flux of each species.

Additionally, we also compare our kinetic model with the DSMC results reported by Strapasson & Sharipov (2014), where the heat transfer in a He–Ar mixture confined between two parallel plates are investigated across a wide range of Knudsen numbers. Figure 7 presents the dimensionless heat flux $q^*$ for a He–Ar mixture with hard-sphere molecules as functions of the rarefaction parameter $\delta$ ranging from 0.01 to 40, for various mole fractions of helium $\chi _{\textrm {He}}=0.25,0.5,0.75$, and temperature differences between two plates $\Delta T/T_0=0.1,0.75$. The maximum relative error, observed in the moderate rarefaction regime, is consistently below 2.8 %. Note that the rarefaction parameter $\delta$ is related to the Knudsen number through the equation

(4.3)\begin{equation} {Kn}_s = \frac{\sqrt{\rm \pi}\mu_s}{2\mu_{mix}}\sqrt{\frac{m_{mix}}{m_s}}\frac{1}{\delta}, \end{equation}

where $\mu _{mix}$ is the viscosity of the mixture. The dimensionless heat flux $q^*$ is defined as

(4.4)\begin{equation} q^* = \frac{q}{p_0}\sqrt{\frac{m_{mix}}{2k_BT_0}}\frac{T_0}{\Delta T}. \end{equation}

Figure 7. Comparison of the dimensionless heat flux between the kinetic model (lines) and the DSMC method (symbols) (Strapasson & Sharipov 2014) for He–Ar mixture with hard-sphere molecules at rarefaction parameter $\delta$ ranging from 0.01 to 40, with mole fractions $\chi _{\textrm {He}}=0.25,0.5,0.75$ and temperature differences between the two plates $\Delta T/T_0=0.1,0.75$. The corresponding parameters in the kinetic model fitted by matching the mixture viscosity and thermal conductivity are $\phi _{12}=2.149$, $\phi _{21}=0.3775$, $\varphi _{12}=1.373$, $\varphi _{21}=1.642$.

4.3. Planar Couette flow

The Couette flow shares the same configuration as that of the Fourier flow, but the temperatures of both plates are kept the same at $T_0$, and the lower and upper plates move along the $y$ direction with velocity $u_{1,y}=-v_w$ and $u_{2,y}=v_w$, respectively, where $v_w=\sqrt {2k_BT_0/m_{mix}}$. Due to the symmetry, only half of the domain $(L_0/2\le {}x_2\le {}L_0)$ is used in the simulations.

The simulation results of our kinetic model and the DSMC method are shown in figure 8 for the three mixtures with $\chi _1=0.5$, when ${Kn}_1=1$ for mixture 1 and ${Kn}_1=0.1$ for the others. Good agreements are found for all the macroscopic properties. We analyse the mixture shear stress predicted by two methods across various mole fractions $\chi _1$ and Knudsen numbers ${Kn}_1$ (table 2). Note that the actual plate's velocities vary significantly for mixtures with different mass ratios and mole fractions, and the boundary velocity is supersonic or even hypersonic for the heavier species, but subsonic for lighter gas in most cases (exclude the one for mixture 1 with $\chi _1=0.9$). Nevertheless, the maximum relative error of shear stress in all the considered cases remains below $3.2\,\%$. It should be noted that the moderate rarefaction level extends to a wider range of Knudsen numbers for a mixture with a larger mass difference, because of the presence of multiscale species ${Kn}$ and relaxation times of inter-species collisions. Importantly, the data demonstrate a consistent level of accuracy for our kinetic model, regardless of variations in the mass ratio of the gas species.

Figure 8. Comparisons of the (a–c) normalized number density and dimensionless (d–f) flow velocity, (g–i) temperature and (j–l) heat flux between the kinetic model (lines) and the DSMC method (symbols) for the Couette flow with the mole fraction of light species $\chi _1=0.5$, when ${Kn}_1=1$ (mixture 1 in the first column) and 0.1 (mixture 2 in the second column and mixture 3 in the third column).

Table 2. The dimensionless shear stress $P_{xy}/n_0k_BT_0$ of the mixture calculated by the kinetic model equation for the Couette flow, with the mole fraction $\chi _1=0.1,0.5,0.9$ and ${Kn}_1=0.1,1,10$. The values in parentheses are the relative errors between the results given by the kinetic model and the DSMC method.

5. Numerical results for two-dimensional problems

In this section the kinetic model is further applied to solve two-dimensional mixture flows, that is, a supersonic mixture passing a circular cylinder and a gas mixture flowing through a nozzle. Also, the results are compared with the DSMC simulations to evaluate the accuracy of our kinetic model.

In the two-dimensional DSMC simulations, we utilize orthogonal grids with refinement in areas of high gas density, ensuring that the cell size does not exceed one-third of the local mean free path of the species with the larger collision diameter. The time step is calculated as the minimum cell size divided by $5v_m^{{light}}$. Each cell contains 50 simulation particles on average.

The discretized velocity method (Zeng et al. 2023b; Liu et al. 2024; Zeng et al. 2024) is applied to solve the model equation that is reduced to be quasi-two-dimensional in velocity space using (4.2a,b). The spatial domains simulated in all the following cases are discretized by structured quadrilateral meshes with refinement near the surfaces, and the grid size is maintained smaller than the local mean free path of the gas molecules, to ensure a reliable solution of the kinetic model equation. The finite volume method is used in the numerical implementation, where the convention fluxes are evaluated implicitly by the lower-upper symmetric Gauss-Seidel technique, and the collision terms are calculated with the Venkata limiter. Further details on the numerical scheme implementation solving the kinetic equation of gas mixtures are available in Zeng et al. (2024).

5.1. Supersonic flow around a circular cylinder

We consider the supersonic gas flow with density $n_0$ at ${Ma}_{\infty }=3$ passing a cylinder with diameter $L_0$. The temperatures of both the free stream and isothermal surfaces of the cylinder are maintained at $T_0$. The Knudsen number of the lighter species in the free stream is ${Kn}_{1}=0.5$. Based on the diameter ratio of the mixtures’ components, the Knudsen number of the heavier species in mixtures 1 and 2 is ${Kn}_{2}=0.5$ as well, while that in mixture 3 of hard-sphere molecules is ${Kn}_{2}=0.125$. Besides, the mole fraction considered in this problem is $\chi _1=0.5$ for all mixtures. The simulations are conducted only in the upper half-domain $[-L,L]\times [0,L]$ due to symmetry, with $L/L_0=6,12,10$ for mixtures 1, 2 and 3, respectively.

The detailed flow fields about species number density and flow velocity of the surrounding gas are presented in figure 9 for the gas mixture with a mass ratio of 1000, where the result given by our kinetic model matches the DSMC data well. Figure 10 compares kinetic model results of the windward side number density, flow velocity, temperature and heat flux along the stagnation line with those solved by the DSMC method, and the overall agreement is very good.

Figure 9. Comparisons of the dimensionless (a) number density and (b) flow velocity between the results of the kinetic model (black lines) and the DSMC method (background contours) for a supersonic mixture flow passing a cylinder, when the mole fraction of the free stream is $\chi _1=0.5$, ${Ma}_{\infty }=3$ and ${Kn}_{1}=0.5$. The binary mixture consists of Maxwell molecules with a mass ratio $m_2/m_1=1000$ and a diameter ratio $d_2/d_1=1$.

Figure 10. Comparisons of the dimensionless (a–c) number density, (d–f) flow velocity in the $x$ direction, (g–i) temperature and (j–l) heat flux in the $x$ direction, along the windward side stagnation line between the results of the kinetic model and the DSMC method for a supersonic mixture flow passing a cylinder, when the mole fraction of the free stream is $\chi _1=0.5$, ${Ma}_{\infty }=3$ and ${Kn}_{1}=0.5$. The first, second and third columns correspond to mixture 1, 2 and 3, respectively.

Despite having the same free-stream flow velocity that exceeds the sound speed of the mixture, individual species in a gas mixture with disparate mass experience vastly different flow characteristics. This is due to their distinct species Mach numbers ${Ma}_{\infty,s}$, which are defined based on their own individual sound speeds $\sqrt {5k_BT_0/3m_{s}}$. For instance, ${Ma}_{\infty,1}=0.13$ and ${Ma}_{\infty,2}=4.24$ for the mixture shown in figure 9. Therefore, the lighter component forms a subsonic flow field and exhibits significantly less compression compared with the heavier species. As shown by the density profiles in figure 10 for the mixture with a mass ratio of 1000, the number density of the heavier gas gets nearly 10 times that of the lighter one near the stagnation point.

The overall properties of the shock are mainly determined by the heavier species, because of its higher number density and molecular mass, as long as the components’ mole fractions in the free stream are comparable. However, it is found that, as the mass ratio changes from 10 to 1000 for mixtures 1 and 2, the thickness of the shock wave and the peak values of the mixture temperature increase only slightly (10 %). Comparatively, mixture 3 has a thinner shock structure since its effective Knudsen number is lower (${Kn}_{2}={Kn}_{1}/4$). Note that the actual free-stream velocities of the mixtures vary significantly due to the distinct average mixture mass, though ${Ma}_{\infty }$ keeps constant.

Figure 11 shows the pressure and heat flux along the surface of the cylinder. The pressure predicted by the kinetic model matches the DSMC results very well, while the heat fluxes given by the two methods have a 6 % relative difference around the windward side stagnation region for mixture 3 consisting of hard-sphere molecules. Interestingly, the forces acting on the object are found to be very close in the three types of mixtures, despite the disparity in average mixture mass, and even the intermolecular potential. Particularly, for the flow of mixtures 1 and 2, the aerodynamic forces are nearly the same and, hence, roughly independent of the mass ratio, which is the only different dimensionless variable in the two flows. Also, the force is found to be insensitive to the intermolecular potential and molecular diameter ratio (mixture 3). On the other hand, the values of heat flux on the surface are not only inversely scaled by the square root of the mass ratio, but also significantly influenced by the intermolecular potential and molecular diameter ratio.

Figure 11. Comparisons of the dimensionless (a–c) pressure and (d–f) heat flux along the surface of the cylinder between the results of the kinetic model and the DSMC method for the mixtures, when the mole fraction of the free stream is $\chi _1=0.5$, ${Ma}_{\infty }=3$ and ${Kn}_{1}=0.5$. Here $p_x$ and $p_y$ denote the pressures exerted on the surface along the $x$ and $y$ directions, respectively. The angle is measured from the windward to the leeward side. The first, second and third columns correspond to mixture 1, 2 and 3, respectively.

5.2. Nozzle flow

Our kinetic model is applied to simulate a two-dimensional rarefied gas mixture flowing through a nozzle into the vacuum. The structure of the nozzle is shown in figure 12, which has a straight channel with width $L_0$, a converging section shrinking the width to $L_0/2$ at the throat and a diverging section. At the inlet of the nozzle ($x=0$), the flow was assumed to be considerably subsonic ${Ma}_{in}=0.05$ with Maxwellian velocity distributions at temperature $T_0$. The gas molecules are reflected on the cold walls of the nozzle ($T_w=T_0/2$) with complete thermal accommodation, and then go through the outlet ($x=3L_0$) into the vacuum. The Knudsen number, defined in terms of the gas properties at the inlet, is ${Kn}_1=0.1$ for lighter species for all the mixtures.

Figure 12. The dimensionless (a,c) number density and (b,d) local Mach number solved by the kinetic model (black lines) and the DSMC method (background contours) for a gas mixture flowing through a nozzle, when ${Kn}_{1}=0.1$ at the inlet. The inlet mole fraction is $\chi _1=0.05$ for mixture 2 (a,b) and $\chi _1=0.5$ for mixture 3 (c,d).

We consider a very small proportion of lighter molecules ($\chi _1=0.05$) mixed with heavier ones for mixture 2 and a half-half mix ($\chi _1=0.5$) of hard-sphere molecules for mixture 3. Figure 12 shows the species number density and local Mach number distribution in the nozzle solved by the kinetic model and the DSMC method. Also, figure 13 shows the density, temperature and heat flux along the centreline of the nozzle. The good accuracy of our model equation is demonstrated.

Figure 13. Comparisons of the dimensionless (a,d) flow velocity, (b,e) temperature and (c, f) heat flux along the centreline between the results of the kinetic model and the DSMC method for a gas mixture flowing through a nozzle, when ${Kn}_{1}=0.1$ at the inlet. The inlet mole fraction is $\chi _1=0.05$ for mixture 2 (a–c) and $\chi _1=0.5$ for mixture 3 (d–f).

It is found that, when a gas mixture with disparate mass flows through a nozzle, the density and temperatures of each component do not experience remarkable separation across the nozzle. However, the species velocities at the outlet become noticeably different. As illustrated in figures 13(a) and 13(d), mixture 2, having a mass ratio of 1000, exhibits an outlet velocity for the lighter species that is 20 times higher than that of the heavier one; similarly, mixture 3 shows a velocity ratio of 3.7 between its lighter and heavier components. Consequently, the significant diffusion velocities make a noticeable contribution to the total heat flux of the mixture, particularly in the diverging section of the nozzle, where the heat flux due to conductance becomes negligible. On the other hand, the two components achieve close values of species Mach number ${Ma}_{s,local}$ calculated using their respective local sound speeds, as plotted in figures 12(b) and 12(d), which also indicates similar degrees of compression of the disparate species flowing through the nozzle. Meanwhile, the relatively higher ${Ma}_{s,local}$ of the heavier components primarily arises from inter-species collisions with lighter molecules, which accelerate the heavier ones.

We also calculate the flow rate of individual components passing through the nozzle outlet, as plotted in figure 14. The molecular number flux is found to be comparable for the mixture components, while the heavier species will hence dominate the mass flow rate due to its much higher molecular mass. Compared with DSMC results, the mass flow rate predicted by the kinetic model yields less than 1.7 % relative difference for all the species in mixture 2 and 3.2 % relative difference for those in mixture 3.

Figure 14. Comparisons of the dimensionless molecular number flux of each species passing through the outlet of the nozzle between the results of the kinetic model (lines) and the DSMC method (symbols). The results of mixtures 2 and 3 are given in (a) and (b), respectively.

6. Conclusions

In summary, we have proposed a new kinetic model for monatomic gas mixtures, which can describe the dynamics of rarefied gas flow with disparate molecular mass, and reduce to the Shakhov model for a single-species gas when the components are mechanically identical and the diffuse velocity vanishes. The tunable parameters in the kinetic model are uniquely determined by the transport properties of the gas mixture, and thus, the shear viscosity, thermal conductivity, diffusion coefficient and thermal-diffusion coefficients can be recovered by the model equation in the continuum limit.

The accuracy of the proposed models has been assessed by comparing with DSMC simulations for various binary gas mixture flows in representative problems, including the one-dimensional Fourier flow, Couette flow and normal shock waves, as well as the two-dimensional supersonic flow passing a cylinder and nozzle flow into a vacuum. A wide range of mass ratios, species concentrations and different intermolecular potentials have been considered. The kinetic model demonstrates its high accuracy not only for predicting the average mixture properties but also for capturing the individual flow fields of components.

The proposed kinetic model benefits from the following features that enable its applicability in the modelling of a rarefied gas mixture with disparate mass.

1. (i) The model equation is constructed using a sum of relaxation operators imitating each type of binary collision individually; thus, it can correctly capture the multiscale relaxation rates inherent in different collision processes. We reveal the dependence of the inter-species relaxation time scales on mass difference, which differ for the lighter and heavier components by orders of magnitude in a mixture with disparate mass.

2. (ii) All the transport coefficients can be recovered simultaneously by the kinetic model in the continuum limit. Particularly, the model goes beyond previous BGK-type models by incorporating the thermal-diffusion effect for mixtures composed of non-Maxwell gases. The concentration variation is correctly captured by our model in simulating the Fourier flow of hard-sphere gas mixtures. Therefore, the kinetic model can be applied to investigate the gas separation processes driven by temperature gradients, such as those occurring in a Knudsen compressor (Takata, Sugimoto & Kosuge 2007).

3. (iii) More importantly, the model equation exhibits minimal loss of accuracy when the components’ mass ratio increases from 10 to 1000, when compared with the original Boltzmann equation. This remarkable consistency suggests the accuracy of the kinetic model for the gas mixtures with even larger mass disparities, and particularly, its potential extension to the kinetic modelling of plasmas.

Gas mixtures with disparate mass possess substantial velocity and temperature non-equilibrium due to significantly slow inter-species relaxations, thus forming the following unique flow characteristics, which are most evident in supersonic mixture flows.

1. (i) Previous observations of temperature overshoots and two-stage structures in normal shock waves are found at low concentrations of the heavier gas in mixtures with moderate mass ratios. However, mixtures with significant mass disparity exhibit these phenomena under different conditions (figures 3 and 4). Specifically, the temperature overshoot occurs when the heavier molecules are present in equal and even lower amounts compared with the lighter molecules. Additionally, the shock wave displays a reversed two-stage structure, characterized by a smooth and expansive upstream region followed by a steep change in properties on the downstream side.

2. (ii) The supersonic mixture around an object leads to the coexistence of a subsonic lighter gas flow and a super/hypersonic heavier gas flow, hence posing a dramatic temperature difference between the components. Interestingly, the aerodynamic force acting on the cylinder is found to be independent of the mass ratio and insensitive to the intermolecular potential, while the heat transfer to the cylinder can be significantly affected by these factors.

Last but not least, a gas mixture with disparate mass usually exhibits multiscale features both spatially and temporally, as well as significant concentration differences between its components. These characteristics bring unaffordable computational costs when solving the Boltzmann equation or conducting DSMC simulations in such mixture flows. On the other hand, the deterministic numerical methods with multiscale schemes that solve the kinetic equations have shown their excellent performance in multiscale problems. For example, the general synthetic iterative scheme developed in recent years can asymptotically preserve the Navier–Stokes equation in the continuum limit (thus removing the constraint on the spatial cell size), and find the steady-state solution of a kinetic equation within dozens of iterations at any Knudsen number (Su et al. 2020a; Su, Zhu & Wu 2020b; Liu et al. 2024). Therefore, with the computationally tractable kinetic model proposed here and the multiscale numerical methods, the gas mixture flows with disparate mass can be accurately and efficiently solved, leading to promising engineering applications, e.g. the particle exhaust system in nuclear fusion devices (Tantos et al. 2024) and the gas dynamics locker in EUV lithography (Teng et al. 2023).

It is crucial to recognize that the model is equally applicable to time-dependent scenarios, given that it accurately captures the relaxation rates of stress, energy and heat flux during intra- and inter-species collisions. Then the multiscale numerical scheme can also be applied to solve the unsteady gas flows (Zeng et al. 2023a). This method stands out as particularly beneficial when contrasted with the DSMC method, which can become excessively time consuming in handling time-dependent conditions of gas mixtures with disparate mass. On the other hand, the proposed kinetic model is solely applicable to dilute gas mixtures, where only binary collisions are considered. In the case of dense gas mixtures, the proposed model can be further extended to dense gas mixtures with additional collision terms describing non-local collision effects, based on our previous work in modelling dense gas flow for a single-species gas (Wu, Zhang & Reese 2015c). Thus, it is expected to have broad implications, potentially solving complex phenomena such as the kinetic processes within non-equilibrium plasmas in inertial confinement fusion. Additionally, although the model is currently designed for monatomic gases, with our experience in the kinetic modelling of single species with internal degrees of freedom (Wu et al. 2015a; Li et al. 2021, 2023), the kinetic models for multi-species gas mixtures with internal degrees of freedom are expected to be established in the near future.