2. The distribution function and macroscopic variables

The kinetic theory of gases provides a statistical description of the motion of molecules in terms of the distribution function, $f(\boldsymbol {x}, \boldsymbol {c}, t)$, which is the probability density of finding a particle in an infinitesimal volume element $\textrm {d} \boldsymbol {x}$ containing the point $\boldsymbol {x}$ with a velocity in the range of $\textrm {d} \boldsymbol {c}$ about the point $\boldsymbol {c}$, where $\boldsymbol {x}$ and $\boldsymbol {c}$ are position and molecular velocity, respectively (Chapman & Cowling Reference Chapman and Cowling1970). Similarly, for an $\mathcal {N}$-component mixture the description is provided in terms of $f_i(\boldsymbol {x}, \boldsymbol {c}, t)\, \textrm {d} \boldsymbol {x}\, \textrm {d} \boldsymbol {c}$, which is the probability of finding a particle of the $i$th type in the neighbourhood of $(\boldsymbol {x}, \boldsymbol {c})$. The relevant macroscopic quantities can then be found by taking the appropriate moment of the distribution function. The component number density $n_i$ and the mixture number density $n$ are

(2.1a,b)\begin{equation} n_{i} = \langle 1, f_{i} \rangle, \quad n = \sum_{i} n_{i}, \end{equation}

where the summation is over all components and the $\langle \phi _1, \phi _2\rangle$ operator is

(2.2)\begin{equation} \langle \phi_1(\boldsymbol{c}_{i}), \phi_2(\boldsymbol{c}_{i}) \rangle = \int _{-\infty}^{\infty} \phi_1(\boldsymbol{c}_{i}) \phi_2(\boldsymbol{c}_{i})\, \textrm{d} \boldsymbol{c}_{i}. \end{equation}

A convention of explicit summation over all components is used for mixture quantities. As an example, the mixture density is defined as $\rho = \sum _{i} \rho _i$, where the component mass density is defined as $\rho _i = m_in_i$ with $m_i$ being the mass of each particle of the $i$th component. Similarly, in $D$ dimensions, the momentum $\rho \boldsymbol {u}$, the energy $E$ and the temperature $T$ of the mixture are defined as

(2.3a–c)\begin{equation} \rho\boldsymbol{u} = \sum_{i} \langle m_{i}\boldsymbol{c}_{i}, f_{i} \rangle,\quad E = \sum_{i} \left\langle \frac{m_{i}{c}_{i}^2}{2}, f_{i} \right\rangle, \quad \frac{D}{2}nk_BT = \sum_{i} \left\langle \frac{m_{i}(\boldsymbol{c}_{i} - \boldsymbol{u})^2}{2}, f_{i} \right\rangle, \end{equation}

where $k_B$ is the Boltzmann constant. Similar to the single-component case, the component velocity $\boldsymbol {u}_{i}$ and component temperature $T_{i}$ are

(2.4a,b)\begin{equation} \rho_{i}\boldsymbol{u}_{i} = \langle m_{i}\boldsymbol{c}_{i}, f_{i} \rangle,\quad \frac{D}{2}n_{i}k_BT_{i} = \left\langle \frac{m_{i}(\boldsymbol{c}_{i} - \boldsymbol{u}_{i})^2}{2}, f_{i} \right\rangle. \end{equation}

Similarly the pressure $p$, the stress $\sigma _{\alpha \beta }$ and the heat flux $q_{\alpha }$ are

(2.5a–c)\begin{equation} p = \sum_{i} n_ik_BT, \quad \sigma_{\alpha\beta} = \sum_{i} \left\langle m_i\overline{\xi_{i\alpha} \xi_{i\beta}}, f_i \right\rangle, \quad q_{\alpha} = \sum_{i} \left\langle m_i\frac{\xi_{i}^2}{2}\xi_{i\alpha}, f_i \right\rangle, \end{equation}

where $\xi _{i\alpha } = c_{i \alpha } - u_{\alpha }$ and with $\overline {{\boldsymbol{\mathsf{A}}}_{\alpha \beta }}$ indicating the traceless part of the tensor.

At equilibrium, the distribution function attains the Maxwell–Boltzmann distribution form

(2.6)\begin{equation} f^{MB}_{i} = n_{i}\left(\frac{m_{i}}{2 {\rm \pi} k_B T } \right)^{D/2} \exp{\left(-\frac{m_{i}}{2k_B T} (\boldsymbol{c_{i}} - \boldsymbol{u})^2 \right)}. \end{equation}

The component velocities and temperatures assume the values of their mixture counterparts, while the stress and heat flux become zero at equilibrium, that is

(2.7a–d)\begin{equation} \boldsymbol{u}_{i}[\,f^{MB}] = \boldsymbol{u}, \quad T_{i}[f^{MB}] = T, \quad \sigma_{\alpha\beta}[f^{MB}] = 0, \quad q_{\alpha}[f^{MB}] = 0. \end{equation}

The kinetic theory of gases also extends the idea of entropy present in statistical mechanics to non-equilibrium situations. This is achieved via the $H$-function defined as

(2.8)\begin{equation} H = \sum_{i} \int \textrm{d} \boldsymbol{c}_{i} \, (\,f_{i}\ln f_{i} - f_{i}). \end{equation}

It can be shown that the $H$-function steadily decreases as the system progresses in time and at equilibrium attains a form similar to the Sackur–Tetrode expression of entropy in thermodynamics (Chapman & Cowling Reference Chapman and Cowling1970), as $S_B = -k_B H[f^{MB}]$, wherein

(2.9)\begin{equation} H[f^{MB}] = \sum_{i} n_i \left[\frac{D}{2}\log\frac{2{\rm \pi} k_BT}{m_i} - \log n_i + \frac{D}{2} \right], \end{equation}

which shows that kinetic theory is consistent with features of statistical mechanics. This completes the description of the various relevant macroscopic quantities that are calculated from the distribution function and their behaviour at equilibrium.

3. The Boltzmann equation

For the case of dilute gases, the time evolution of the distribution function is described by the Boltzmann equation (Chapman & Cowling Reference Chapman and Cowling1970). For the single-component dilute gas, the Boltzmann equation has the form

(3.1)\begin{equation} \partial _t \,f (\boldsymbol{x}, \boldsymbol{c}, t ) + c_{\alpha}\partial_{x_{\alpha}} \,f(\boldsymbol{x}, \boldsymbol{c}, t ) = \Omega, \end{equation}

where $\Omega$ is a term that accounts for the change due to collisions between particles. The Boltzmann collisional operator, $\Omega ^{B}$, quantifies the change in the distribution function from all possible binary collisions, and is expressed as

(3.2)\begin{equation} \Omega^{B} = \iiint [w^{\prime}f (\boldsymbol{x}, \boldsymbol{ c}^{\prime}, t ) f(\boldsymbol{x}, \boldsymbol{c}^{\prime}_1, t ) - w\,f (\boldsymbol{x}, \boldsymbol{c}, t ) f(\boldsymbol{x}, \boldsymbol{c}_1, t )] \,\textrm{d}\boldsymbol{c}_{1}\,\textrm{d}\boldsymbol{c}^{\prime}\, \textrm{d}\boldsymbol{c}^{\prime}_{1}, \end{equation}

where $w \equiv w(\boldsymbol {c}^{\prime }, \boldsymbol {c}^{\prime }_{1}; \boldsymbol {c}, \boldsymbol {c}_{1})$ is the probability of the colliding pair to transition from the velocities $(\boldsymbol {c}, \boldsymbol {c}_{1})$ to $(\boldsymbol {c}^{\prime },\boldsymbol {c}^{\prime }_{1})$ and vice versa for $w^{\prime }$, and it can be shown that $w = w^{\prime }$ (Lifschitz & Pitaevskii Reference Lifschitz and Pitaevskii1981). The first term of the integrand represents the increase (gain) in the value of the distribution function $f(\boldsymbol {x}, \boldsymbol {c}, t)$ and similarly the second term represents the decrease (loss). In order to satisfy the conservation of momentum and energy, the velocity pairs must satisfy

(3.3a,b)\begin{equation} \boldsymbol{c} + \boldsymbol{c}_{1} = \boldsymbol{c}^{\prime} + \boldsymbol{c}^{\prime}_{1}, \quad c^2 + c^{2}_1 = c^{\prime 2} + c^{\prime 2}_{1}. \end{equation}

Therefore, by integrating over all possible $(\boldsymbol {c}_{1}, \boldsymbol {c}^{\prime }, \boldsymbol {c}^{\prime }_{1})$, the total change in $f(\boldsymbol {x}, \boldsymbol {c}, t)$ from collisions can be calculated. Further, by considering appropriate moments and integrating over the velocity space, the dynamics of various macroscopic quantities can be derived (Grad Reference Grad1949; Chapman & Cowling Reference Chapman and Cowling1970).

The Boltzmann equation is a highly complex integro-differential equation and hence does not lend itself to analysis even for simple boundary value problems. Therefore, approximations are made to the collisional term in order to find solutions to the equation. A simple, yet highly useful, model is the BGK approximation (Bhatnagar et al. Reference Bhatnagar, Gross and Krook1954), where the Boltzmann collisional operator is modelled as an approach to the equilibrium distribution function. The BGK collisional operator is

(3.4)\begin{equation} \Omega^{BGK} = \frac{1}{\tau_{BGK}} (f^{MB} - f ), \end{equation}

where $\tau _{BGK}$ is the relaxation time.

A variant of this approach is the ES-BGK model (Holway Reference Holway1966), which has the form

(3.5)\begin{equation} \Omega^{ES} = \frac{1}{\tau_{ES}} (f^{ES} - f ), \end{equation}

where $\tau _{ES}$ is the relaxation time associated with this model and $f^{ES}$ is the anisotropic Gaussian with the form

(3.6)\begin{equation} f^{ES} = \frac{n}{\sqrt{\textrm{det}[2{\rm \pi}\lambda_{\alpha\beta}]}} \exp\left(-\frac{1}{2}\xi_{\alpha}\lambda_{\alpha\beta}^{-1}\xi_{\beta}\right), \end{equation}

where $\textrm {det}[~]$ is the determinant and $\lambda _{\alpha \beta }$ is

(3.7)\begin{equation} \lambda_{\alpha\beta} = \frac{k_B T}{m}\delta_{\alpha\beta} + b\frac{\sigma_{\alpha\beta}}{\rho}, \end{equation}

where the $b$ parameter is used to tune the Prandtl number, as opposed to the BGK model where the Prandtl number is set to $1$.

Another method of approximation is the Fokker–Planck operator (Lebowitz et al. Reference Lebowitz, Frisch and Helfand1960), wherein the approach to equilibrium is modelled as drift and diffusion dynamics

(3.8)\begin{equation} \Omega^{FP} = \frac{1}{\tau_{{FP}}} \partial _{c_{\alpha}} \left(\xi _{\alpha}\,f + \frac{k_BT}{m}\frac{\partial f}{\partial c_{\alpha}} \right) = \frac{1}{\tau_{FP}}\frac{k_BT}{m}\partial_{c_{\alpha}} (f \partial_{c_{\alpha}} (\ln f - \ln f^{MB} ) ), \end{equation}

which is essentially the diffusion dynamics in velocity space, with $\xi _{\alpha }$ acting as the drift coefficient, $k_BT/m$ assumes the role of diffusion coefficient and $\tau ^{-1}_{{FP}}$ is the friction constant. It has been shown recently that this approximation is a useful alternative to methods such as discrete simulation Monte Carlo (DSMC) for simulating moderately high-$Kn$ flows.

The Boltzmann equation and its aforementioned approximations are shown to have the following properties

1. (a) Conservation laws: As binary collisions do not change the mass, momentum or energy of the system, we have

   (3.9)\begin{equation} \langle\Omega, \{m, m\boldsymbol{c}, mc^2/2\}\rangle = \{0, \boldsymbol{0}, 0\}. \end{equation}
   Using this result and calculating appropriate moments of the Boltzmann equation, the conservation laws are
   (3.10)\begin{equation} \left.\begin{array}{c@{}} \partial_{t} \rho + \partial_{\alpha}\rho u_{\alpha} = 0,\\ \partial_{t} \rho u_{\alpha} + \partial_{\beta}(\rho u_{\alpha}u_{\beta}+p \delta_{\alpha \beta})+\partial_{\beta} \sigma_{\alpha \beta}= 0, \\ \partial_{t}E+\partial_{\alpha} (( E+p) u_{\alpha} +\sigma_{\alpha \gamma} u_{\gamma} )+\partial_{\alpha} q_{\alpha} = 0,\end{array}\right\} \end{equation}
   which are in accordance with the macroscopic laws of conservation.

2. (b) Zero of the collision: When the collisions between particles do not affect the state of the system, it reaches a state of equilibrium:

   (3.11)\begin{equation} \Omega = 0 \quad\implies\quad f = f^{MB}. \end{equation}
   The converse is also true, i.e. when $f = f^{MB}$ then $\Omega = 0$.

3. (c) $H$-theorem: The Boltzmann equation extends the idea of entropy to non-equilibrium situations. This is highlighted from the evolution of the $H$-function,

   (3.12)\begin{equation} \partial _t H + \partial _{\alpha} J_{\alpha} = -\sigma^{S}, \end{equation}
   where the $H$-function is $\langle (\ln f - 1), f \rangle$, $J_{\alpha }$ is the entropy flux term and $\sigma ^{S} = \langle \Omega , \ln f \rangle$ is the entropy generation term. The Boltzmann collisional operator, the BGK and Fokker–Planck approximations ensure that
   (3.13)\begin{equation} \sigma^{S} \geq 0. \end{equation}
   Entropy production is greater than or equal to zero and hence the Boltzmann equation for rarefied gases is in accordance with the laws of thermodynamics. It is also noted that entropy production is zero at equilibrium, that is, when $f = f^{MB}$.

This concisely summarizes the features of the Boltzmann equation and its approximations for the single-component case. In the following section, we briefly explain the Boltzmann equation for binary mixtures, and some well-known approximations, and outline their important features.

4. Kinetic modelling of binary mixtures

The dynamics of binary mixtures is fundamentally different from the single-component case, as the two components exchange momentum and energy through collisions. Hence the Boltzmann equation for binary mixtures considers the different collisional possibilities, as schematically shown in figure 1. The changes in distribution of component A arise from A–A and A–B type collisions, and vice versa for component B. As a result, the mixture momentum and energy are conserved as opposed to their componentwise counterparts, i.e.

(4.1a,b)\begin{equation} m_{A}\boldsymbol{c}_{A} + m_{B}\boldsymbol{c}_{B} = m_{A}\boldsymbol{c}^{\prime}_{A} + m_{B}\boldsymbol{c}^{\prime}_{B}, \quad m_{A}c_{A}^2 + m_{B}c^{2}_{B} = m_{A} c^{\prime 2}_{A} + m_{B} c^{\prime 2}_{B}, \end{equation}

where $m_{A}$ and $m_{B}$ are the mass of A-type and B-type particles, respectively. Here $\boldsymbol {c}$ and $\boldsymbol {c^{\prime }}$ with a subscript denote the pre-collision and post-collision velocities of the respective particles. The Boltzmann equation for binary mixtures is (Chapman & Cowling Reference Chapman and Cowling1970)

(4.2)\begin{equation} \left.\begin{array}{c@{}} \dfrac{\partial f_{A}}{\partial t} + c_{A\alpha}\dfrac{\partial f_{A}}{\partial x_{\alpha}} = \Omega^{{B}}_{A} = \underbrace{\Omega^{{B}}(f_{A}, f_{A})}_{{self\textit{-}collision}} + \underbrace{\Omega^{{B}}(f_{A}, f_{B})}_{{cross\textit{-}collision}},\\ \dfrac{\partial f_{B}}{\partial t} + c_{B\alpha}\dfrac{\partial f_{B}}{\partial x_{\alpha}} = \Omega^{{B}}_{B} = \underbrace{\Omega^{{B}}(f_{B}, f_{B})}_{{self\textit{-}collision}} + \underbrace{\Omega^{{B}}(f_{B}, f_{A})}_{{cross\textit{-}collision}},\end{array}\right\} \end{equation}

where the right-hand side of each equation is the change in distribution of the respective component arising from self-collisions represented by $\Omega ^{{B}}(f_{i}, f_{i})$ and cross-collisions represented by $\Omega ^{{B}}(f_{i}, f_{j})$. Similar to the single-component case, the evolution equation of various macroscopic variables can be derived using this equation. The collisional operator has the following properties, which should ideally be satisfied by its approximations.

1. (a)

   Figure 1.

   The three types of collisional possibilities: A–A, A–B and B–B.

   

   Conservation laws: The mass of individual species as well as the total momentum and energy of the mixture are conserved, as binary collisions do not contribute any change to these quantities. This is represented as

   (4.3)\begin{equation} \left\langle \Omega^{{B}}_{i}, \left\{ m_{i}, \sum_{i=A,B} m_{i}\boldsymbol{c}_{i}, \sum_{i=A,B} m_{i}\frac{c^{2}_{i}}{2} \right\} \right\rangle = \{ 0, \boldsymbol{0}, 0 \} , \end{equation}
   using which the conservation laws can be calculated similarly to the single-component case. However, the componentwise momentum and energy are not conserved, as the two components exchange momentum and energy between themselves through cross-collisions (A–B type collisions). It is, in fact, these collisions that facilitate the relaxation of the component momentum and energy to the mixture momentum and energy (Hamel Reference Hamel1965). The mixture variables adhere to the conservation laws as mentioned in (3.10).

2. (b) Equilibrium: Similar to the single-component case, the system reaches a state of statistical equilibrium. The distribution of any component $i$ at equilibrium is

   (4.4)\begin{equation} f^{MB}_{i} = n_{i}\,\left(\frac{m_{i}}{2 {\rm \pi} k_B T } \right)^{3/2} \exp{\left(-\frac{m_{i}}{2 k_B T} (\boldsymbol{c}_{i}-\boldsymbol{u})^2 \right)}. \end{equation}
   The converse is also true, i.e. when $\Omega ^{B}_{i} = 0$ then the distribution function attains the form $f = f^{MB}_{i}$.

3. (c) H-theorem: The Boltzmann collision kernel for binary mixtures satisfies the $H$-theorem, that is, $\sigma ^{(s)} \geq 0$ similar to the single-component case.

4. (d) Indifferentiability principle: The equations should adhere to the indifferentiability principle, i.e. the equation should converge to the single-component case for $m_A = m_B$.

Similar to the single-component case, the corresponding BGK collision kernel for a binary mixture is

(4.5)\begin{equation} \Omega^{BGK} = \frac{1}{\tau}(f^{MB}_{i}(\rho_i, \boldsymbol{u}, T) - f_{i}). \end{equation}

The fundamental drawback with such a model is that there is only a single relaxation rate for all quantities, whereas, for the case of a binary mixture, there are two important time scales present in the system – the rate of mass diffusion and the rate of momentum diffusion. The dimensionless parameter that is used to characterize these time scales is known as the Schmidt number and is defined as (Bergman et al. Reference Bergman, Incropera, DeWitt and Lavine2011)

(4.6)\begin{equation} Sc = \frac{\text{viscous diffusion rate}}{\text{mass diffusion rate}} = \frac{\mu}{\rho D_{AB}}, \end{equation}

where $\mu$ is the viscosity, $\rho$ is the density and $D_{AB}$ is the mass diffusion coefficient. For such an approximation, $Sc = 1$ for all cases, and hence does not manage to accurately describe the system. Therefore, the collision kernel should be approximated in a manner capable of preserving these different time scales.

Various approaches to correct this defect exist. In order to deal with multiple time scales, the basic idea of fast–slow decomposition of motions near the quasi-equilibrium was introduced (Gorban & Karlin Reference Gorban and Karlin1994a,Reference Gorban and Karlinb; Levermore Reference Levermore1996). In accordance with this idea, the relaxation to equilibrium is modelled as a two-step process, where ‘fast’ relaxation happens from the initial to the quasi-equilibrium state and ‘slow’ relaxation happens from the quasi-equilibrium state to the final equilibrium state. In the context of multiple time scales, the quasi-equilibrium models are a simple alternative (Gorban & Karlin Reference Gorban and Karlin1994a,Reference Gorban and Karlinb; Levermore Reference Levermore1996) to the BGK approximation which can effectively incorporate multiple time scales of the system. The collision kernel for the quasi-equilibrium model is

(4.7)\begin{equation} \Omega^{{QE}}_{i} = \frac{1}{\tau_1}(f^{*}_{i}(M^{{quasi\text{-}slow}}, M^{{slow}}) - f_{i}) + \frac{1}{\tau _2}(f^{{MB}}_{i}(M^{{slow}}) - f^{*}_{i}(M^{{quasi\text{-}slow}}, M^{{slow}})), \end{equation}

where $f^{*}_{i}(M^{{quasi\text {-}slow}}, M^{{slow}})$ is the quasi-equilibrium distribution function and is a function of the quasi-slow moments $M^{{quasi\text {-}slow}}$ and the slow moments $M^{{slow}}$ (Gorban & Karlin Reference Gorban and Karlin1994b; Levermore Reference Levermore1996). The idea is that the system moves towards a state of quasi-equilibrium, where the quasi-slow moments relax first, and then proceeds towards equilibrium, where the slow moments react; a visual description of the idea is presented in figure 2. Similar to the approach adopted by Arcidiacono et al. (Reference Arcidiacono, Karlin, Mantzaras and Frouzakis2007), two possible forms for the quasi-equilibrium distribution are constructed: for low $Sc$ where mass diffusion occurs at a higher rate as compared to momentum diffusion, and vice versa for high $Sc$.

Figure 2.

Trajectory of the distribution function in a quasi-equilibrium model. The various concentric circles depict different levels of $H$.

The two distinct regimes of the Schmidt number can be explained at a kinetic level in the following manner.

1. (a) When the self-collisions (A–A and B–B types) dominate and have shorter time scales, the low-Schmidt-number regime is observed. In this case, the individual species first relax to the Maxwell–Boltzmann distribution corresponding to the component velocity and temperature, following which they relax to the Maxwell–Boltzmann distribution with mixture velocity and temperature at longer time scales.

2. (b) When the cross-collisions (A–B type) have a shorter characteristic time scale, the system has a higher Schmidt number, as the particles of different species collide at a higher rate and exchange momentum and energy, thereby initially relaxing to the mixture velocity and temperature and over a long time attain the Maxwell–Boltzmann distribution.

For the first case, the physically relevant quasi-slow variables are

(4.8)\begin{equation} M^{{quasi\text{-}slow}} = \{\rho_{i}, \rho_{i}\boldsymbol{u}_{i}, n_{i}k_B T_{i}\}, \end{equation}

which imposes the following condition on the quasi-equilibrium distribution function $f^{*}_{i}$:

(4.9)\begin{equation} \left\langle \left\{m_{i}, m_{i}\boldsymbol{c}_{i}, m_{i}\frac{(\boldsymbol{c}_{i} - \boldsymbol{u}_{i})^2}{2} \right\}, f^{*}_{i} \right\rangle = \left\{\rho _{i}, \rho _{i} \boldsymbol{u}_{i}, \frac{3}{2}n_{i}k_B T_{i} \right\}. \end{equation}

By minimizing the $H$-function as defined in (2.8) under these constraints, the form of quasi-equilibrium for the low-Schmidt-number limit, $f^{*(L)}_{i}$,is (Arcidiacono et al. Reference Arcidiacono, Mantzaras, Ansumali, Karlin, Frouzakis and Boulouchos2006)

(4.10)\begin{equation} f^{*(L)}_{i} = n_{i}\left(\frac{m_{i}}{2{\rm \pi} k_B T_{i}}\right)^{3/2}\exp{\left(-\frac{m_{i}(c_{i} - u_{i})^2}{2k_BT_{i}}\right)}. \end{equation}

Similarly, for the second case where the momentum diffuses faster, the set of constraints under which the $H$-function is to be minimized is

(4.11)\begin{equation} \left\langle \left\{m_{i}, m_{i}\boldsymbol{c}_{i}, \sum_{i=A,B} m_{i}\xi_{i\alpha}\xi_{i\beta}\right\}, f^{*}_{i} \right\rangle = \{\rho _{i}, \rho _{i} \boldsymbol{u}, n\theta_{\alpha\beta} \}, \end{equation}

where

(4.12)\begin{equation} \theta _{\alpha\beta} = \frac{1}{n}\sum_{i=A,B}\left\langle m_{i}\xi_{i\alpha}\xi_{i\beta}, f_{i} \right\rangle. \end{equation}

The quasi-equilibrium distribution function for the high-Schmidt-number limit, $f^{*(H)}_{i}$, is (Arcidiacono et al. Reference Arcidiacono, Mantzaras, Ansumali, Karlin, Frouzakis and Boulouchos2006)

(4.13)\begin{equation} f^{*(H)}_{i} = n_{i}\left(\frac{m_{i}}{2{\rm \pi} |\theta_{\alpha\beta}|} \right)^{{3}/{2}} \exp\left(\frac{-m_{i}\xi_{i\alpha}\theta^{-1}_{\alpha\beta}\xi_{i\beta}}{2}\right), \end{equation}

where $|\theta _{\alpha \beta }|$ is the determinant.

These two distinct forms of quasi-equilibrium can be used to build two different collision kernels based on the Fokker–Planck approximation, which can be solved for binary mixtures.

5. Quasi-equilibrium models for Fokker–Planck formulation

The Fokker–Planck approximation to the Boltzmann equation of the form (3.8) involves only a single time scale and therefore is not well suited for modelling systems with multiple time scales. To resolve this problem, Gorji & Jenny (Reference Gorji and Jenny2012) proposed a model that could account for the different collision frequencies and with a cubic drift term that enabled correct calculation of the transport coefficients. We propose to extend the Fokker–Planck approximation for binary mixtures by incorporating the concept of quasi-equilibrium models in a manner that would correctly represent the multiple time scales present in the system and its approach to equilibrium.

From (3.8), the Fokker–Planck approximation is

(5.1)\begin{equation} \Omega^{FP} = \frac{1}{\tau_{ FP}}\frac{k_BT}{m}\partial_{c_{\alpha}} (f\partial_{c_{\alpha}} (\ln f - \ln f^{MB} )). \end{equation}

This form of the Fokker–Planck model better illustrates the approach of the Maxwell–Boltzmann distribution, similar to the BGK model. In order to build a quasi-equilibrium like model with multiple time scales, we extend the Fokker–Planck model for it to have a similar form. Here, the approach to equilibrium is defined as a two-step process wherein the first term represents a logarithmic approach to quasi-equilibrium and the second to equilibrium:

(5.2)\begin{equation} \Omega^{FP}_{i} = \left(\frac{1}{\tau_{1}} - \frac{1}{\tau_{2}}\right) {\boldsymbol{\mathsf{A}}}_{\alpha\beta}\partial_{c_{i\beta}} (f_i \partial_{c_{i\alpha}} (\ln f_i - \ln f_{i}^{*} ) ) + \frac{1}{\tau_{2}}\frac{k_BT}{m_{i}}\partial_{c_{i\alpha}} (f_i \partial_{c_{i\alpha}} (\ln f_i - \ln f_{i}^{MB} ) ), \end{equation}

where $\tau _1$ and $\tau _2$ are the characteristic time scales associated with the approach to quasi-equilibrium and equilibrium, and ${{\boldsymbol{\mathsf{A}}}_{\alpha \beta }}$ is the diffusion coefficient that relaxes the system to the quasi-equilibrium state, for the low-Schmidt-number dynamics; ${{\boldsymbol{\mathsf{A}}}_{\alpha \beta }}$ can be chosen as ${\boldsymbol{\mathsf{A}}}_{\alpha \beta } = (k_{B} T_{i}/ m_{i}) \delta _{\alpha \beta }$. Using the form of $f^{*}_{i}$ presented in (4.10), the collision kernel for the low-Schmidt-number limit is

(5.3)\begin{align} \Omega_{i}^{FP(L)} = \frac{1}{\tau_{1}}\partial_ {c_{i\alpha}} \left( (c_{i\alpha} - u_{i\alpha}) f_{i} + \frac{k_{B} T_{i}}{m_{i}} \frac{\partial f_{i}}{\partial c_{i\alpha}}\right) + \frac{1}{\tau_{2}}\partial_ {c_{i\alpha}} \left( (u_{i\alpha} - u_{\alpha}) f_{i} + \frac{k_{B} {\rm \Delta} T}{m_{i}} \frac{\partial f_{i}}{\partial c_{i\alpha}}\right), \end{align}

where ${\rm \Delta} T = T - T_{i}$ is the difference in component and mixture temperatures.

Similarly, for the high-Schmidt-number dynamics, ${\boldsymbol{\mathsf{A}}}_{\alpha \beta }$ is taken as ${\boldsymbol{\mathsf{A}}}_{\alpha \beta } = \theta _{\alpha \beta }$, and the collision kernel for the high-Schmidt-number limit is

(5.4)\begin{equation} \Omega^{FP(H)}_{i} = \frac{1}{\tau_{1}}\partial_ {c_{i\alpha}} \left( (c_{i\alpha} - u_{\alpha}) f_{i} + \frac{ \theta_{\alpha \beta}}{m_{i}} \frac{\partial f_{i}}{\partial c_{i\beta}}\right) + \frac{1}{\tau_{2}} \partial_ {c_{i\alpha}} \left( \left( \frac{k_{B} T \delta_{\alpha \beta}}{m_{i}} - \frac{ \theta_{\alpha \beta}}{m_{i}} \right)\frac{\partial f_{i}}{\partial c_{i\beta}}\right). \end{equation}

For this model to be considered canonical, it must satisfy the properties of collision as mentioned in § 4. By integrating over the velocity space $\boldsymbol {c}_{i}$, it can be verified that the quasi-equilibrium Fokker–Planck model satisfies the constraints of (4.3). The evolution equations for component mass, mixture momentum and energy are the same as the conservation laws mentioned in (3.10). Furthermore, the component momentum and energy equations in relaxation form are

(5.5)\begin{equation} \left.\begin{array}{c@{}} \left\langle \Omega^{{FP}(L)}_{i}, m_{i}c_{i\alpha} \right\rangle = \dfrac{1}{\tau_{2}} ( \rho_{i} u_{\alpha} - \rho_{i} u_{i\alpha} ), \\ \left\langle \Omega^{{FP}(L)}_{i}, \dfrac{m_{i}c_{i}^{2}}{2} \right\rangle = \dfrac{1}{\tau_{2}} ( \rho_{i} u_{i\alpha} (u_{\alpha} - u_{i\alpha}) + D k_{B} n_{i} (T - T_{i}) ), \end{array}\right\} \end{equation}

for the low-Schmidt-number case. Similarly, for the high-Schmidt-number case the relaxation equations for component momentum and energy are

(5.6)\begin{equation} \left.\begin{array}{c@{}} \left\langle \Omega^{{FP}(H)}_{i}, m_{i}c_{i\alpha} \right\rangle = \dfrac{1}{\tau_{1}} ( \rho_{i} u_{\alpha} - \rho_{i} u_{i\alpha} ), \\ \left\langle \Omega^{{FP}(H)}_{i}, \dfrac{m_{i}c_{i}^{2}}{2} \right\rangle = \dfrac{1}{\tau_{1}} (\rho_{i} u_{i\alpha} (u_{\alpha} - u_{i\alpha}) + D k_{B} n_{i} (T - T_{i}) ). \end{array}\right\} \end{equation}

If $\tau _1 \leq \tau _2$, the component velocities equilibrate faster in the second case than in the first, which is as expected since the second model is applicable for the high-Schmidt-number regime wherein the viscous diffusion rate dominates the mass diffusion rate.

For the proposed model, the expression for entropy generation ($\sigma ^S$) is

(5.7)\begin{equation} \sigma^{S} = \frac{1}{\tau_{{eff}}} \sum_{i} - D n_{i}+ \int \frac{1}{f_{i}} \frac{\partial f_{i}}{\partial c_{i\alpha} }{\boldsymbol{\mathsf{A}}} _{\alpha\beta}\frac{\partial f_{i}}{\partial c_{i\beta}} + \frac{1}{\tau_{2}} \sum_{i} - D n_{i} + \int \frac{k_{B} T}{m_{i}}\frac{1}{f_{i}} \frac{\partial f_{i}}{\partial c_{i\alpha} } \frac{\partial f_{i}}{\partial c_{i\alpha}}, \end{equation}

where $\tau _{{eff}} = \tau _{2}\tau _{1}/(\tau _{2} - \tau _{1})$ and ${\boldsymbol{\mathsf{A}}}_{\alpha \beta } = (k_BT/m_i) \delta _{\alpha \beta }$ for the low-Schmidt-number case and ${\boldsymbol{\mathsf{A}}}_{\alpha \beta } = \theta _{\alpha \beta }/m_i$ for the high-Schmidt-number case. Following Singh & Ansumali (Reference Singh and Ansumali2015), (5.7) can be rewritten as

(5.8)\begin{equation} \sigma^{S} = \sum_{i} \underbrace{\int f_{i} \frac{\partial \ln ( f_{i}/ f^{*}_{i} )} {\partial c_{i\alpha}} {\boldsymbol{\mathsf{A}}}_{\alpha\beta} \frac{\partial \ln ( f_{i}/ f^{*}_{i} )}{\partial c_{i\beta}} \,\textrm{d}\boldsymbol{c}_{i}} _{positive} + \sum_{i} \underbrace{\int \frac{k_{B}T}{m_{i}}f_{i} \left(\frac{\partial \ln ( f_{i}/ f^{*}_{i} )}{\partial c_{i\alpha}}\right)^{2} \textrm{d}\boldsymbol{c}_{i}}_{positive}, \end{equation}

which suggests that

(5.9)\begin{equation} \sigma^S \geq 0, \quad \forall \ \tau_{1} \leq \tau_{2}. \end{equation}

Therefore, the proposed model satisfies the $H$-theorem for $\tau _1 \leq \tau _2$.

An important condition for $\Omega ^{FP}$ to be considered valid is that the zero of collision must imply that the distribution function has attained a Maxwell–Boltzmann form. For the present model, the zero of collision, i.e. $\Omega _{i}^{FP} = 0$, implies

(5.10a,b)\begin{equation} \left \langle \Omega^{{FP}(1)}_{i}, m_{i}c_{i\alpha} \right \rangle = 0, \quad \left \langle \Omega^{{FP}(1)}_{i}, \frac{m_{i}c_{i}^{2}}{2} \right \rangle = 0, \end{equation}

then as per (5.5) equilibrium $u_{i\alpha } = u_{\alpha }$ and $T_{i} = T$; therefore $\Omega _{i}^{FP(L)} = 0$ then reduces to

(5.11)\begin{equation} \partial_{c_{i\alpha}}\left((c_{i\alpha} - u_{\alpha})f_{i} + \frac{k_{B}T}{m_{i}}\frac{\partial f_{i}}{\partial c_{i\alpha}}\right) =0. \end{equation}

Integrating (5.11) with respect to the velocity space and using the fact that the distribution function and its derivatives tend to zero at infinity, we have

(5.12)\begin{equation} (c_{i\alpha} - u_{\alpha})f_{i} + \frac{k_{B}T}{m_{i}}\frac{\partial f_{i}}{\partial c_{i\alpha}} = 0. \end{equation}

Solving (5.12), we get the Maxwell–Boltzmann distribution as the solution. To find the equilibrium distribution function for the high-Schmidt-number case, we first note that

(5.13)\begin{equation} \sum_{i}\left\langle \Omega^{{FP}(H)}_{i}, m_{i}\xi_{i\alpha}\xi_{i\beta} \right\rangle = \frac{2}{\tau_{2}} \left(\frac{k_BT}{m_i} - \frac{\theta_{\alpha\beta}}{m_i}\right), \end{equation}

which suggests that $\Omega ^{{FP}(H)} = 0 \implies \theta _{\alpha \beta }/m_i = (k_BT/m_i)\delta _{\alpha \beta }$. Hence, $\Omega ^{{FP}(H)} = 0$ then reduces to (5.11), the solution for which is the Maxwell–Boltzmann distribution.

Furthermore, the model must be consistent with the indifferentiability principle, i.e. one must be able to recover the Fokker–Planck approximation for the single-component case. In the case where $\tau _1 = \tau _2 = \tau$ and $m_A = m_B = m$, the Fokker–Planck collision kernel for binary mixtures outlined in (5.2) reduces to the approximation for the single-component case, indicating that the proposed model abides by theindifferentiability principle.

As demonstrated above, the proposed model does indeed satisfy the conservation laws, the $H$-theorem, the zero of collision and the indifferentiability principle. Thus, this model is an acceptable approximation to the Boltzmann equation for binary mixtures.

6. Transport coefficients

In order to obtain the transport coefficients, we perform the Chapman–Enskog expansion, wherein the time derivative, distribution function and other relevant variables are represented as a series, with ${Kn}$ acting as the smallness parameter (Chapman & Cowling Reference Chapman and Cowling1970). The distribution function is expressed in series form as

(6.1)\begin{equation} f_{i} = f^{{MB}}_{i} + {Kn}\,f_{i}^{(1)} + {Kn}^2\,f_{i}^{(2)} + \cdots, \end{equation}

with the following constraints imposed on $f_{i}$:

(6.2)\begin{equation} \left\langle f_{i}^{(n)}, \left\{ m_{i}, \sum_{i=A,B} m_{i}c_{i\alpha}, \sum_{i=A,B} m_{i}c_{i}^2/2 \right\} \right\rangle = \{0, \boldsymbol{0}, 0\}, \quad \forall \ n \geq 1. \end{equation}

These constraints ensure that component density, mixture momentum and energy are slow moments. The higher-order moments in series form are

(6.3)\begin{equation} \left.\begin{array}{c@{}} \sigma_{\alpha\beta} = {Kn}\,\sigma _{\alpha\beta}^{(1)} + {Kn}^2\,\sigma_{\alpha\beta}^{(2)} + \cdots, \\ q_{\alpha} = {Kn}\,q_{\alpha}^{(1)} + {Kn}^2\,q_{\alpha}^{(2)} + \cdots,\end{array}\right\} \end{equation}

as the stress and heat flux are zero at equilibrium and are expected to be functions of slow moments otherwise. The time derivative is expressed in series form as

(6.4)\begin{equation} \partial _t = \partial _t^{(0)} + {Kn}\,\partial_t^{(1)} + {Kn}^2\,\partial_t^{(2)} + \cdots. \end{equation}

The time derivative of $f_i$ at zeroth order is computed using (Liboff Reference Liboff2003)

(6.5)\begin{equation} \partial_{t}^{(0)} f_i^{MB}(\rho_i, \boldsymbol{u}, T) = \frac{\partial f_i^{MB}}{\partial \rho_i} \partial_{t}^{(0)} \rho_i + \frac{\partial f_i^{MB}}{\partial \boldsymbol{u}} \partial_{t}^{(0)} \boldsymbol{u} + \frac{\partial f_i^{MB}}{\partial T} \partial_{t}^{(0)} T, \end{equation}

where the expression for the time derivatives of the conserved variables can be calculated from the conservation laws mentioned in (3.10)

In order to find an expression for the viscosity, we first calculate the stress evolution equation. For the first model it has the form

(6.6)\begin{align} &\partial_{t}\sigma_{\alpha \beta} + \partial_{\gamma} (\sigma_{\alpha \beta} u_{\gamma})+ 2 p \overline{\partial_{\alpha} u_{\beta}}+2\overline{\sigma_{\alpha\gamma}\partial_{\gamma} u_{\beta}}+\partial_{\gamma}Q_{\alpha \beta \gamma}+ \frac{4}{D+2} \overline{\partial_{\alpha}q_{\beta}} \nonumber\\ &\quad=- \frac{2}{\tau_{1}} \left(\sigma_{\alpha \beta} + \rho \overline{u_{\alpha} u_{\beta}} - \sum_{i=A,B}\rho_{i}\overline{ u_{i\alpha}u_{i\beta}}\right) - \frac{2}{\tau_{2}} \left( \sum_{i=A,B} \rho_{i} \overline{ u_{i\alpha} u_{i\beta}} - \rho \overline{ u_{\alpha} u_{\beta}}\right), \end{align}

where $Q_{\alpha \beta \gamma } = \sum _{i=A,B} \langle m_{i}\overline {\xi _{i\alpha }\xi _{i\beta }\xi _{i\gamma }} \rangle$. Retaining terms up to $\mathcal {O}({Kn})$, the stress evolution equation yields

(6.7)\begin{equation} 2 p \overline{\partial_{\alpha}u_{\beta}} = -\frac{2\sigma_{\alpha\beta}^{(1)}}{\tau_{1}} , \end{equation}

and comparing with the Navier–Stokes law for the stress tensor, we have

(6.8)\begin{equation} \mu = \frac{p\tau_1}{2}. \end{equation}

Similarly, for the second model the right-hand side of the stress evolution equation is

(6.9)\begin{equation} \sum_{i=A,B} \langle m_{i}\overline{\xi_{i\alpha}\xi_{i\beta}}, \Omega^{{FP}(2)}_{i} \rangle = - \frac{2}{\tau_{2}} \sigma_{\alpha \beta}. \end{equation}

Hence, the expression for viscosity for this model is

(6.10)\begin{equation} \mu = \frac{p \tau_{2}}{2}. \end{equation}

Similarly, the expression for the diffusion coefficient can be calculated by considering the relaxation of diffusion flux defined as

(6.11)\begin{equation} V_{\alpha} = m_{AB}(u_{A\alpha} - u_{B\alpha}), \end{equation}

where $m_{AB} = (\rho _A \rho _B)/\rho$. Diffusive flux essentially quantifies the difference between the momentum of a given component and the momentum of the mixture. The series expansion for this quantity is

(6.12)\begin{equation} V_{\alpha} = {Kn}\,V^{(1)}_{\alpha} + {Kn}^{2}\,V^{(2)}_{\alpha} + \cdots . \end{equation}

Similar to stress and heat flux, at equilibrium the diffusive flux attains zero values, as the momenta of both components relax to the mixture momentum. In order to calculate the expression for $V_{\alpha }$, we write the expression for individual component velocities. For the first model, we have

(6.13)\begin{equation} \left.\begin{array}{c@{}} \partial _{t} \rho_{A} u_{A\alpha} + \partial _{\alpha} P_{A\alpha\beta} = \dfrac{1}{\tau_2}(\rho u_{\alpha} - \rho_{A}u_{A\alpha}),\\ \partial _{t} \rho_{B} u_{B\alpha} + \partial _{\alpha} P_{B\alpha\beta} = \dfrac{1}{\tau_2}(\rho u_{\alpha} - \rho_{B}u_{B\alpha}),\end{array} \right\} \end{equation}

where $P_{i\alpha \beta } = \langle m_{i}c_{i\beta }c_{i\beta } \rangle$ and at equilibrium attains the value $P_{i\alpha \beta } = p_{i}\delta _{\alpha \beta } + \rho _{i}u_{\alpha }u_{\beta }$. After subtracting one equation from the other and collecting terms up to $\mathcal {O}(Kn)$, we have

(6.14)\begin{equation} \partial_t^{(0)} (\rho_A - \rho_B) u_{\alpha} + \partial_{\beta} [ (n_A - n_B)k_BT_0 \delta_{\alpha\beta} + (\rho_A - \rho_B)u_{\alpha}u_{\beta} ]= -\frac{2}{\tau_2}V_{\alpha\beta}^{(1)}.\end{equation}

The temporal derivatives are replaced using

(6.15a,b)\begin{equation} \partial_{t}^{(0)} \rho_{i} = -\partial_{\alpha}(\rho_{i}u_{\alpha}), \quad \partial_{t}^{(0)} \rho u_{\alpha} = -\partial_{\beta} (nk_BT_0 \delta_{\alpha\beta} + \rho u_{\alpha}u_{\beta}). \end{equation}

After some rearrangement (6.14) takes the form

(6.16)\begin{equation} V_{\alpha}^{(1)} = \tau_2 [ Y_{A}\partial_{\alpha}p - p\partial_{\alpha}X_{A} - X_{A}\partial_{\alpha}p ], \end{equation}

where $X_i = n_i/n$ is the component mole fraction and $Y_i = \rho _{i}/\rho$ is the component mass fraction. Rearranging (6.16) we have

(6.17)\begin{equation} \partial_{\alpha} X_{A} = - \frac{V_{\alpha}^{(1)}}{ \tau_{2}p} + (Y_{A} - X_{A}) \frac{\partial_{\alpha} p}{p}. \end{equation}

This has the same form as the Stefan–Maxwell equation (Bergman et al. Reference Bergman, Incropera, DeWitt and Lavine2011), which governs the diffusion in multi-component systems, and for binary mixtures is

(6.18)\begin{equation} \partial_{\alpha} X_{A} = \frac{X_{A} X_{B}}{D_{AB}} \frac{V_{\alpha}}{m_{AB}} + (Y_{A} - X_{A}) \frac{\partial_{\alpha} p}{p}. \end{equation}

Comparing (6.17) with the Stefan–Maxwell equation, we get the following expression for the diffusion coefficient:

(6.19)\begin{equation} D_{AB} = X_{A} X_{B} \frac{p}{m_{AB}} \tau_{2}. \end{equation}

The Schmidt number can now be computed as

(6.20)\begin{equation} Sc = \frac{\mu}{\rho D_{AB}} = \frac{\tau_{1}}{2\tau_{2}} \frac{m_{AB}}{X_{A} X_{B}}\frac{1}{\rho} = \frac{\tau_{1}}{ 2\tau_{2}} \frac{Y_{A} Y_{B}}{X_{A} X_{B}}. \end{equation}

The existence of the $H$-theorem for this model suggests that $\tau _1 \leq \tau _2$, hence

(6.21)\begin{equation} Sc \leq \frac{Y_{A} Y_{B}}{2 X_{A} X_{B}}. \end{equation}

This model has an upper limit on the Schmidt number and this is in accordance with the characteristics of the quasi-equilibrium distribution. Similarly, for the second model, the Schmidt number is calculated as

(6.22)\begin{equation} Sc = \frac{\mu}{\rho D_{AB}} = \frac{\tau_{2}}{2\tau_{1}} \frac{m_{AB}}{X_{A} X_{B}}\frac{1}{\rho} = \frac{\tau_{2}}{ 2\tau_{1}} \frac{Y_{A} Y_{B}}{X_{A} X_{B}}, \end{equation}

and, since the limitation $\tau _1 \leq \tau _2$ exists, as consistent with the hypothesis, there is a lower bound on the Schmidt number, which is

(6.23)\begin{equation} Sc \geq \frac{Y_{A} Y_{B}}{2 X_{A} X_{B}}. \end{equation}

Hence, both models in conjunction can cover a large range of Schmidt numbers.

7. Numerical scheme

A Fokker–Planck equation that describes the evolution of the probability density function of the random variable $\boldsymbol {\eta }$ is of the form

(7.1)\begin{align} \frac{\textrm{d}\mathcal{P}(\boldsymbol{\eta}, t)}{\textrm{d} t} &= -\Lambda^{(1)}_{\alpha}(\boldsymbol{\eta},t)\frac{\partial \mathcal{P}(\boldsymbol{\eta}, t)}{\partial \eta_{\alpha}} + \frac{\zeta^{(1)}_{\alpha\beta}(\boldsymbol{\eta},t)}{2}\frac{\partial^{2} \mathcal{P}(\boldsymbol{\eta}, t)}{\partial \eta_{\alpha}\partial \eta_{\beta}} \nonumber\\ &\quad- \Lambda^{(2)}_{\alpha}(\boldsymbol{\eta},t)\frac{\partial \mathcal{P}(\boldsymbol{\eta}, t)}{\partial \eta_{\alpha}} + \frac{\zeta^{(2)}_{\alpha\beta}}{2}(\boldsymbol{\eta},t)\frac{\partial^{2} \mathcal{P}(\boldsymbol{\eta}, t)}{\partial \eta_{\alpha}\partial \eta_{\beta}}, \end{align}

where $\boldsymbol {\Lambda }^{(i)}$ are the drift terms and $\boldsymbol {\zeta }^{(i)}$ are the diffusion coefficients. This form of the Fokker–Planck equation is equivalent to the Langevin equation (Risken Reference Risken1996)

(7.2)\begin{equation} \dot{\eta}_{\alpha} = h^{(1)}_{\alpha} (\boldsymbol{\eta}, t) + g^{(1)}_{\alpha \beta} (\boldsymbol{\eta}, t) \Gamma_{\beta}(t) + h^{(2)}_{\alpha}(\boldsymbol{\eta}, t) + g^{(2)}_{\alpha\beta} (\boldsymbol{\eta}, t ) \Gamma_{\beta}^{\prime}(t), \end{equation}

where $\boldsymbol {h}^{(i)}$ are the drift terms, $\boldsymbol {g}^{(i)}$ are the diffusion coefficients, and $\boldsymbol {\Gamma }, \boldsymbol {\Gamma }^{\prime }$ are Gaussian distributed random numbers that have the following properties:

(7.3a,b)\begin{equation} \langle \Gamma_{\alpha}(t) \rangle = 0, \quad \langle \Gamma_{\alpha}(t) \Gamma_{\beta}(t^{\prime}) \rangle = \delta(t-t^{\prime}) \delta_{\alpha\beta}. \end{equation}

Under these conditions, the following relations hold (Risken Reference Risken1996):

(7.4)\begin{equation} \left.\begin{array}{c@{}} \Lambda^{(1)}_{\alpha} = h^{(1)}_{\alpha}(\boldsymbol{\eta}, t), \quad \zeta^{(1)}_{\alpha\beta} = g^{(1)}_{\alpha\gamma}g^{(1)}_{\gamma\beta},\\ \Lambda^{(2)}_{\alpha} = h^{(2)}_{\alpha}(\boldsymbol{\eta}, t), \quad \zeta^{(2)}_{\alpha\beta} = g^{(2)}_{\alpha\gamma}g^{(2)}_{\gamma\beta}. \end{array} \right\} \end{equation}

The central idea is that the solution to the Fokker–Planck equation is approximated by considering an ensemble of trajectories generated by the Langevin dynamics. In this case, a large number of particles have their positions and velocities updated using (7.2). We now discuss the numerical scheme for the two cases.

7.1. Low-Schmidt-number limit

For the first model, the equivalent Langevin equations are

(7.5)\begin{equation} \left.\begin{array}{c@{}} \dfrac{{\textrm{d}x}_{\alpha}}{\textrm{d}t} = c_{i\alpha},\\ \dfrac{\textrm{d}c_{i\alpha}}{\textrm{d}t} = - \left(\dfrac{1}{\tau_{{eff}}} \right) ( c_{i\alpha} - u_{i\alpha}) - \dfrac{1}{\tau_{2}}( c_{i\alpha} - u_{\alpha}) + \sqrt{\dfrac{2 k_{B} T_{i}}{m_{i}}}\,\textrm{d} W_{\alpha} + \sqrt{\dfrac{2 k_{B} T}{m_{i}}} \,\textrm{d} W_{\alpha}^{\prime},\end{array}\right\} \end{equation}

where $\textrm {d}W_{\alpha }$ and $\textrm {d}W_{\alpha }^{\prime }$ denote random forces with the following statistics:

(7.6a,b)\begin{equation} \langle\textrm{d}W_{\alpha} \rangle = 0, \quad \langle\textrm{d}W_{\alpha}^{\prime}\rangle = 0, \quad \langle\textrm{d}W_{\alpha} \,\textrm{d} W_{\alpha}^{\prime}\rangle = 0. \end{equation}

More specifically, $\textrm {d}W = W(t+{\rm \Delta} t) - W(t)$ is the standard Weiner process, where $W(t)$ is a rapidly changing random force with mean and variance as (Gardiner Reference Gardiner1985)

(7.7a,b)\begin{equation} \langle \textrm{d} W_{\alpha}(t)\rangle = 0, \quad \langle \textrm{d} W_{\alpha} \,\textrm{d} W_{\beta}\rangle = \textrm{d}t \,\delta_{\alpha \beta}. \end{equation}

Thus, the detailed binary collision description is approximated by a random collision with a heat bath in the model.

These Langevin equations can be solved efficiently using the stochastic version of the Verlet algorithm. For the present model, the discretization scheme we have used is (Kloeden & Platen Reference Kloeden and Platen2013; Singh & Ansumali Reference Singh and Ansumali2015)

(7.8)\begin{align} \left.\begin{array}{c@{}} x^{(1)}_{\alpha} = x_{\alpha}(t) + \dfrac{1}{2} c_{i\alpha}(t) {\rm \Delta} t,\\ c_{i\alpha}(t + {\rm \Delta} t) = c_{i\alpha}(t) - \left(\dfrac{\vartheta_{1}}{1 + \vartheta_{1}/2}\right) (c_{i\alpha}(t) - u_{i\alpha}) - \left(\dfrac{\vartheta_{2}}{1 + \vartheta_{2}/2} \right) (c_{i\alpha}(t) - u_{\alpha})\\ + \dfrac{\sqrt{2\mathcal{D}^{(1)}_{i}\vartheta_{1}}}{1 + \vartheta_{1}/2}\phi_{\alpha} + \dfrac{\sqrt{2\mathcal{D}^{(2)}_{i}\vartheta_{2}}}{1 + \vartheta_{2}/2}\phi^{\prime}_{\alpha}, \\ x_{\alpha}(t + {\rm \Delta} t) = x_{\alpha}^{(1)} + \dfrac{1}{2} c_{i\alpha}(t + {\rm \Delta} t) {\rm \Delta} t, \end{array} \right\} \end{align}

where $\vartheta _{1} = {\rm \Delta} t/\tau _{{eff}}$, $\vartheta _{2} = {\rm \Delta} t/\tau _{2}$, $\phi _{\alpha }$ and $\phi _{\alpha }^{\prime }$ are Gaussian random numbers with mean zero and variance unity, and $\mathcal {D}^{(1)}_{i}$ and $\mathcal {D}^{(2)}_{i}$ are $k_BT_i/m_i$ and $k_BT/m_i$, respectively. The recently proposed Molecular Dice algorithm (Agrawal, Bhattacharya & Ansumali Reference Agrawal, Bhattacharya and Ansumali2018) was used to generate these Gaussian random numbers, which indicated considerable increase in efficiency without any loss of accuracy. This scheme works efficiently for small time steps such that $\max \{\vartheta _{1}, \vartheta _{2}\} \leq 0.001$.

7.2. High-Schmidt-number limit

The formulation for this model remains largely unchanged and the equivalent Langevin equations are

(7.9)\begin{equation} \left.\begin{array}{c@{}} \dfrac{{\textrm{d}x}_{\alpha}}{\textrm{d}t} = c_{i\alpha},\\ \dfrac{\textrm{d}c_{i\alpha}}{\textrm{d}t} = - \left(\dfrac{1}{\tau_{{eff}}} \right) ( c_{i\alpha} - u_{\alpha}) - \dfrac{1}{\tau_{2}}( c_{i\alpha} - u_{\alpha}) + \sqrt{2} \theta^{\prime}_{i\alpha\beta} \,\textrm{d}W_{\beta} + \sqrt{\dfrac{2 k_{B} T}{m_{i}}} \,\textrm{d}W_{\alpha}^{\prime}, \end{array}\right\} \end{equation}

where $\theta ^{\prime }_{i\alpha \gamma }\theta ^{\prime }_{i\gamma \beta } = \theta _{i\alpha \beta }/m_{i}$ and $\theta ^{\prime }_{\alpha \beta }$ can be calculated by using Cholesky decomposition of $\theta _{\alpha \beta }/m_{i}$. The discretization scheme for this model is

(7.10)\begin{equation} \left.\begin{array}{c@{}} x^{(1)}_{\alpha} = x_{\alpha}(t) + \dfrac{1}{2} c_{i\alpha}(t) {\rm \Delta} t,\\ c_{i\alpha}(t + {\rm \Delta} t) = c_{i\alpha}(t) - \left(\dfrac{\vartheta_{1}}{1 + \vartheta_{1}/2}\right) (c_{i\alpha}(t) - u_{\alpha}) - \left(\dfrac{\vartheta_{2}}{1 + \vartheta_{2}/2} \right) (c_{i\alpha}(t) - u_{\alpha})\\ + \dfrac{\sqrt{2\vartheta_{1}}\theta^{\prime}_{i\alpha\beta}}{1 + \vartheta_{1}/2}\phi_{\beta} + \dfrac{\sqrt{2\mathcal{D}^{(2)}_{i}\vartheta_{2}}}{1 + \vartheta_{2}/2}\phi^{\prime}_{\alpha},\\ x_{\alpha}(t + {\rm \Delta} t) = x_{\alpha}^{(1)} + \dfrac{1}{2} c_{i\alpha}(t + {\rm \Delta} t) {\rm \Delta} t. \end{array}\right\} \end{equation}

In order to validate the numerical scheme, we started with a mixture with $m_B/m_A = 2$ with $N = 10^5$ particles in a single periodic box. For model I, the velocities of the lighter particles were initialized uniformly in the range $[0,1)$ and the heavier particles in the range $[0,2)$. For model II, the velocities of the lighter particles were initialized with a Gaussian distribution with mean $4$ and variance $10$, and the heavier particles were Gaussian distributed with mean $1$ and variance $1$. The plots of energy of the two components and the mixture with time averaged over an ensemble of $15$ trajectories and the distribution of velocities in a particular direction, for both cases, are shown in figures 3 and 4.

Figure 3.

Plot of the distribution of velocities of the light and heavy components at equilibrium for (a) model I and (b) model II.

Figure 4.

Plot of the ratio of energy at time $t$ to the initial energy ($E(t)/E_0$) versus time for individual components and the mixture for (a) model I and (b) model II.

8. Simulation results

In this section, we present the results for three benchmark problems: Graham's law for effusion, Couette flow and binary diffusion.

8.1. Graham's law for effusion

Effusion is a process wherein gas molecules escape through a small hole. The length parameter of this hole is much smaller than the mean free path of the gas, i.e. $d \ll \lambda _{{mfp}}$. A sketch of the process is shown in figure 5. The number flux of the gas through this small hole is

(8.1)\begin{equation} \Phi_{i} = \langle c_{iz}, f(\boldsymbol{c}_{i}) \rangle, \end{equation}

where $\Phi _{i}$ is the number flux and $c_{iz}$ is the molecular velocity in the direction perpendicular to the plane of the hole. By integrating over velocity space, facilitated by a shift to the spherical coordinate system, the expression for $\Phi _{i}$ is

(8.2)\begin{equation} \Phi_{i} = \frac{P}{\sqrt{2{\rm \pi} m_{i}k_BT}}, \end{equation}

where $P$ is the pressure and $T$ is the temperature of the gas. Then, for a well-mixed binary mixture, the ratio of the fluxes is (Mason & Kronstadt Reference Mason and Kronstadt1967)

(8.3)\begin{equation} \frac{\Phi_{A}}{\Phi_{B}} = \sqrt{\frac{m_{B}}{m_{A}}}. \end{equation}

Figure 5.

The lighter particles (in this case filled circles) escape through the hole faster than the heavier particles, with a factor proportional to the square root of their mass ratios.

We simulated this system for three mass ratios, i.e. $m_{B}/m_{A} = 4$, $m_{B}/m_{A} = 16$ and $m_{B}/m_{A} = 100$. The boundary conditions in the transverse directions were taken to be periodic while maintaining constant pressure in the system. The results have been plotted in figure 6. As can be seen, the simulations are in excellent agreement with the analytical solution.

Figure 6.

Model I was used to simulate a set-up that could mimic Graham's law for effusion. The results observed are in very good agreement with the expected behaviour, for all three cases.

8.2. Couette flow

The set-up of the problem is simple: fluid between two plates is sheared in opposite directions with equal magnitudes. A sketch of the problem is shown in figure 7. In order to validate the model, we calculate the global stress tensor defined as (Sharipov, Cumin & Kalempa Reference Sharipov, Cumin and Kalempa2004)

(8.4)\begin{equation} \Pi = -\frac{v_0}{2UP_0}P_{xy}, \end{equation}

where $P_0 = nk_BT_0$ is the reference pressure. This quantity is calculated in the entire range of rarefaction parameter, $\delta$, which is the inverse of the Knudsen number and is defined as

(8.5)\begin{equation} \delta = \frac{HP_0}{\mu v_0}, \end{equation}

where $\mu$ is the mixture viscosity and $v_0$ is the characteristic molecular velocity of the mixture, defined as

(8.6)\begin{equation} v_0 = \sqrt{\frac{2k_BT_0}{m_0}}, \end{equation}

where $m_0 = C_0m_A + (1-C_0)m_B$, with $C_0$ being the concentration of the lighter component. We simulated the system for three mixtures, neon–argon (Ne–Ar), helium–argon (He–Ar) and helium–xenon (He–Xe), for rarefaction parametersin the range $[0.01, 40]$ for three different concentrations $(0.1,0.5,0.9)$. The value for $\Pi$ was computed by averaging over $10^5$ iterations for each parameter and the results are tabulated in table 1. The error bar (standard deviation) was of the same order for all parameters and in the range $[0.00478, 0.00544]$. The results were found to be in good agreement with the reported results of Sharipov et al. (Reference Sharipov, Cumin and Kalempa2004). This indicates that the proposed method is indeed capable of simulating flows in a wide range of Knudsen numbers.

Figure 7.

A representative sketch of the Couette flow set-up. Two walls with a separation $H$ are sheared in opposite directions with velocity $U/2$.

Table 1.

The $\Pi$ values for Ne–Ar, He–Ar and He–Xe mixtures for three different compositions.

8.3. Binary diffusion

The profile of the mixture in this set-up is determined by the step function

(8.7)\begin{equation} \left.\begin{array}{c@{}} X_A = 90\,\%, \quad X_B = 10\,\% \quad \text{if } x < 0, \\ X_B = 10\,\%, \quad X_B = 90\,\% \quad \text{if } x \geq 0, \end{array}\right\} \end{equation}

where the mass ratio of the components was chosen to be $m_B/m_A = 5$. The step function is used instead of a smooth profile as it is a more severe check for the numerical scheme. Under the assumption that, at infinity, the initial concentrations remains unchanged, this problems yields the analytical solution (Bergman et al. Reference Bergman, Incropera, DeWitt and Lavine2011)

(8.8)\begin{equation} X_i = \left[\frac{1}{2} + \frac{{\rm \Delta} X_i}{2} \, \mathrm{erf} \left(\frac{x}{\sqrt{4D_{AB}t}} \right)\right], \end{equation}

where $D_{AB}$ is the diffusion constant. The simulation was done for $20\,000$ time steps, and the plots for both components compared against their respective analytical solutions are plotted in figure 8. The simulation results were very close to the analytical solution. This exercise proves that the value of $D_{AB}$ set by the numerical scheme is accurate.

Figure 8.

Plot of the concentrations after 20 000 time steps of the component (a) A and (b) B, compared with the analytical solution given by (8.8).

9. Outlook

We have developed a new thermodynamically consistent Fokker–Planck approximation to the Boltzmann equation for binary gas mixtures, based on quasi-equilibrium models. These models were subjected to numerical experiments, like Graham's law, Couette flow and binary diffusion, and it was determined that the algorithm is capable of simulating flow for a wide range of Knudsen numbers and diffusion coefficients.

Finally, we wish to point out that this Fokker–Planck model of binary mixtures needs to be extended for multi-component mixtures and possibly to reactive systems due to its practical relevance related to reaction–diffusion systems. This extension, though it is conceptually straightforward, does require some new technical ingredients. We leave the development of such a model to future work, but we would like to highlight some of these technicalities before concluding. A straightforward model for an $\mathcal {N}$-component mixture model can be formulated by introducing separate relaxation times and quasi-equilibrium for each component in the proposed collision approximation (5.2). One would expect that most characteristics of the current model will be preserved by such an extension. The tunable relaxation time is obtained by comparing the hydrodynamic limit of the kinetic equation with its macroscopic counterpart, the Stefan–Maxwell diffusion equation for an isothermal system. However, as the diffusion coefficient is an $\mathcal {N}\times \mathcal {N}$ symmetric matrix while the number of tunable parameters is only $\mathcal {N}$, such a representation is not complete. This issue can be resolved by utilizing the mixture-averaged diffusion approximation (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot1960). Such formulations are correct only asymptotically when there is either only one bulk component and $\mathcal {N}-1$ trace species, or when all species but one have nearly the same diffusion coefficient. A more detailed model should introduce $\mathcal {N}(\mathcal {N}-1)/2$ independent relaxation times. We leave it to future work to explore these issues and possibly extend the model to also incorporate chemical reactions in the kinetic equation framework.