1 Introduction↩︎

In their seminal 1954 paper [1], Bhatnagar, Gross and Krook (BGK) proposed a phenomenological model describing kinetic processes in a pure gas, and two years later, Gross and Krook extended this result to gas mixtures [2]. Even though neither of these models follows from the first principles, they are believed to provide a qualitatively correct approximation of the Boltzmann kinetic equation, and a lot of work has been done to generalize and extend the BGK approach. In application to mixtures, the effort has mostly gone into making the BGK model more adaptable, so that it would be able to describe a wide range of real fluids (e.g., Refs. [3]–[17]).

Note, however, that the multispecies BGK model has never been tested for compliance with the Onsager reciprocal relations, which impose certain constraints on the transport coefficients. Models derived from the first principles satisfy them automatically, whereas phenomenological models may or may not do so. An example of a non-compliant model can be viewed in the Enskog theory of dense fluids [18], and an example of a compliant one, in the so-called modified Enskog theory [19], [20]. The noncompliance with the Onsager relations casts doubt on the model’s physical relevance, and it is no coincidence that the modified Enskog theory has eventually been shown to follow from the first principles for a fluid of hard spheres [21], [22].

As demonstrated in the present paper, the most common version of the BGK model (which includes the original result of Gross and Krook [2] as a particular case) does not comply with the Onsager reciprocal relations. According to one of those, the coefficient of the temperature gradient in the mass flux should be inter-linked in a certain way with the coefficient of the density gradient in the heat flux. According to the BGK model, however, the former is zero, whereas the latter is proportional to the coefficient of the density gradient in the mass flux – hence, cannot be zero. Not only does this undermine the physical relevance of the model, this also makes the BGK model impossible to calibrate – i.e., choose the values of the parameters involved to ensure that the transport properties of the fluid under consideration are described correctly.

In Sec. 2 of this paper, one the most general BGK-type models will be formulated, and in Sec. 3, it will be shown to not comply with the Onsager relations. Other BGK-type models are briefly discussed in Sec. 4. For simplicity, only binary (two-species) mixtures will be considered, but the resulting conclusions apply to the general case as well.

2 Formulation: the BGK model for binary mixtures↩︎

Consider a mixture of two monatomic gases, described by the distribution functions \(f_{i}(t,\mathbf{r},\mathbf{v})\), where \(i\) is the species number, \(t\) is the time, \(\mathbf{r}\) is the position vector, and \(\mathbf{v}\), the molecular velocity. The macroscopic number density \(n_{i}\), velocity \(\mathbf{V}_{i}\), and temperature \(T_{i}\) of the \(i\)-th species are given by\[n_{i}={\displaystyle\int} f_{i}\mathrm{d}^{3}\mathbf{v}, \label{2461}\tag{1}\] \[n_{i}\mathbf{V}_{i}={\displaystyle\int} \mathbf{v}f_{i}\mathrm{d}^{3}\mathbf{v}, \label{2462}\tag{2}\] \[3n_{i}T_{i}={\displaystyle\int} m_{i}\left\vert \mathbf{v}-\mathbf{V}_{i}\right\vert ^{2}f_{i}\mathrm{d}^{3}\mathbf{v}, \label{2463}\tag{3}\] where \(m_{i}\) is the molecular mass, and \(T_{i}\) is measured in energy units (so that the Boltzmann constant equals unity).

The most general form of the multispecies BGK model (e.g., [15], [23]) consists in\[\begin{align} \frac{\partial f_{1}}{\partial t}+\mathbf{v}\cdot\mathbf{\nabla}f_{1} & =\nu_{11}\left( M_{1}-f_{1}\right) +\nu_{12}\left( M_{12}-f_{1}\right) ,\tag{4}\\ \frac{\partial f_{2}}{\partial t}+\mathbf{v}\cdot\mathbf{\nabla}f_{2} & =\nu_{22}\left( M_{2}-f_{2}\right) +\nu_{21}\left( M_{21}-f_{2}\right) , \tag{5} \end{align}\] where \(\nu_{ij}\) are the frequencies of collisions between the molecules of \(i\)-th and \(j\)-th species,\[M_{i}=n_{i}\left( \frac{m_{i}}{2\pi T_{i}}\right) ^{3/2}\exp\left( -\frac{m_{i}\left\vert \mathbf{v}-\mathbf{V}_{i}\right\vert ^{2}}{2T_{i}}\right) , \label{2466}\tag{6}\] \[\begin{align} M_{12} & =\left( \frac{m_{1}}{2\pi T_{12}}\right) ^{3/2}n_{1}\exp\left( -\frac{m_{1}\left\vert \mathbf{v}-\mathbf{V}_{12}\right\vert ^{2}}{2T_{12}}\right) ,\tag{7}\\ M_{21} & =\left( \frac{m_{2}}{2\pi T_{21}}\right) ^{3/2}n_{2}\exp\left( -\frac{m_{2}\left\vert \mathbf{v}-\mathbf{V}_{21}\right\vert ^{2}}{2T_{21}}\right) , \tag{8} \end{align}\] are various Maxwellian distributions, and\[\begin{align} \mathbf{V}_{12} & =\mathbf{V}_{1}+\beta_{1}\left( \mathbf{V}_{2}-\mathbf{V}_{1}\right) ,\tag{9}\\ \mathbf{V}_{21} & =\mathbf{V}_{2}+\beta_{2}\left( \mathbf{V}_{1}-\mathbf{V}_{2}\right) , \tag{10} \end{align}\] \[\begin{align} T_{12} & =T_{1}+\alpha_{1}\left( T_{2}-T_{1}\right) +\gamma_{1}\left\vert \mathbf{V}_{1}-\mathbf{V}_{2}\right\vert ^{2},\tag{11}\\ T_{21} & =T_{2}+\alpha_{2}\left( T_{1}-T_{2}\right) +\gamma_{2}\left\vert \mathbf{V}_{2}-\mathbf{V}_{1}\right\vert ^{2}. \tag{12} \end{align}\] Note that the parameters \(\nu_{ij}\), \(\alpha_{i}\), \(\beta_{i}\), and \(\gamma_{i}\) may depend on the macroscopic characteristics \(n_{1}\), \(n_{2}\), \(\mathbf{V}_{1}\), \(\mathbf{V}_{2}\), etc. – hence, may vary with \(t\) and \(\mathbf{r}\), but not with \(\mathbf{v}\). Various particular cases of model (1 )–(12 ) have been examined in Refs. [2], [3], [5]–[8], [11], [24].

Eqs. (1 )–(12 ) form a closed set for \(f_{1}(t,\mathbf{r},\mathbf{v})\) and \(f_{2}(t,\mathbf{r},\mathbf{v})\). One can readily show that they conserve mass – i.e., satisfy\[\frac{\partial n_{1}}{\partial t}+\mathbf{\nabla}\cdot\left( n_{1}\mathbf{V}_{1}\right) =0,\qquad\frac{\partial n_{2}}{\partial t}+\mathbf{\nabla}\cdot\left( n_{2}\mathbf{V}_{2}\right) =0. \label{24613}\tag{13}\] As for the momentum and energy, Eqs. (1 )–(12 ) do not conserve them automatically, but only subject to the following constraints:\[\alpha_{1}=\frac{\alpha}{\nu_{12}},\qquad\alpha_{2}=\frac{\alpha}{\nu_{21}}, \label{24614}\tag{14}\] \[\beta_{1}=\frac{\beta}{\nu_{12}m_{1}},\qquad\beta_{2}=\frac{\beta}{\nu _{21}m_{2}}, \label{24615}\tag{15}\] \[\begin{align} \gamma_{1} & =\frac{1}{3\nu_{12}}\left( \beta-\frac{\beta^{2}}{\nu _{12}m_{1}}+3\gamma\right) ,\tag{16}\\ \gamma_{2} & =\frac{1}{3\nu_{21}}\left( \beta-\frac{\beta^{2}}{\nu _{21}m_{2}}-3\gamma\right) , \tag{17} \end{align}\] where the coefficients \(\alpha\), \(\beta\), and \(\gamma\) may depend on \(\mathbf{r}\) and \(t\). The above constraints are equivalent to those derived in Refs. [15], [23], albeit presented in a different form.

Given constraints (14 )–(17 ), Eqs. (1 )–(12 ) imply that\[\begin{gather} \frac{\partial\left( m_{1}n_{1}\mathbf{V}_{1}+m_{2}n_{2}\mathbf{V}_{2}\right) }{\partial t}\\ +\mathbf{\nabla}\cdot{\displaystyle\int} \mathbf{v}\otimes\mathbf{v}\left( m_{1}f_{1}+m_{2}f_{2}\right) \mathrm{d}^{3}\mathbf{v}=0, \label{24618} \end{gather}\tag{18}\] \[\begin{gather} \frac{\partial}{\partial t}\left( \frac{3}{2}n_{1}T_{1}+\frac{m_{1}\left\vert \mathbf{V}_{1}\right\vert ^{2}}{2}+\frac{3}{2}n_{2}T_{2}+\frac{m_{1}\left\vert \mathbf{V}_{2}\right\vert ^{2}}{2}\right) \\ +\mathbf{\nabla}\cdot{\displaystyle\int} \frac{\left\vert \mathbf{v}\right\vert ^{2}}{2}\mathbf{v}\left( m_{1}f_{1}+m_{2}f_{2}\right) \mathrm{d}^{3}\mathbf{v}=0, \label{24619} \end{gather}\tag{19}\] which reflect the momentum and energy conservation, respectively.

3 Transport fluxes under the diffusion approximation↩︎

3.1 The standard hydrodynamics↩︎

Within the framework of the Enskog–Chapman approach (e.g., Ref. [25], chapter 6), the mass and heat fluxes are given by\[\mathbf{J}_{i}=-m_{i}n_{i}\left( \sum_{j}D_{ij}\mathbf{d}_{j}+B_{i}\frac{\mathbf{\nabla}T}{T}\right) ,\label{3461}\tag{20}\] \[\mathbf{Q}=-\kappa\mathbf{\nabla}T-\sum_{i}C_{i}\mathbf{d}_{i},\label{3462}\tag{21}\] where\[\begin{gather} \mathbf{d}_{j}=\mathbf{\nabla}\frac{n_{j}}{n_{1}+n_{2}}\\ +\left( \frac{n_{j}}{n_{1}+n_{2}}-\frac{\rho_{j}}{m_{1}n_{1}+m_{2}n_{2}}\right) \frac{\mathbf{\nabla}p}{p},\label{3463} \end{gather}\tag{22}\] is the “diffusion driving force”of the \(i\)-th species, and the pressure is\[p=\left( n_{1}+n_{2}\right) T.\label{3464}\tag{23}\] \(D_{ij}\), \(B_{i}\), \(C_{i}\), and \(\kappa\) are the transport coefficients: \(D_{ij}\) is the diffusivity, \(B_{i}\) is the thermodiffusivity (it describes the Soret effect, i.e., the mass flux due to a temperature gradient), \(C_{i}\) describes the Dufour effect (i.e., heat flux due to a concentration gradient), and \(\kappa\) is the thermal conductivity.

Most importantly, \(C_{i}\) is linked to \(B_{i}\) via one of the Onsager reciprocal relations,\[C_{i}=pB_{i}.\] and the other Onsager relation requires that the diffusivity matrix be symmetric,\[D_{ij}=D_{ji}.\] In addition to the above relations, the coefficients \(D_{ij}\) and \(B_{i}\) should also satisfy\[\sum_{i}D_{ij}=0,\qquad\sum_{i}B_{i}=0\] (see Ref. [25]). Thus, for a binary mixture, one can express \(D_{ij}\), \(B_{i}\), and \(C_{i}\) through only two coefficients – say, \(D\) and \(B\) – so that\[D_{11}=\frac{\rho_{2}}{\rho_{1}}D,\qquad D_{22}=\frac{\rho_{1}}{\rho_{2}}D,\label{3465}\tag{24}\] \[D_{21}=D_{12}=-D,\label{3466}\tag{25}\] \[B_{1}=\frac{B}{m_{1}n_{1}},\qquad B_{2}=-\frac{B}{m_{2}n_{2}},\label{3467}\tag{26}\] \[C_{1}=\frac{\left( n_{1}+n_{2}\right) T}{m_{2}n_{2}}B,\qquad C_{2}=-\frac{\left( n_{1}+n_{2}\right) T}{m_{1}n_{1}}B.\label{3468}\tag{27}\] In this paper, expressions (20 )–(27 ) will be used under the diffusion approximation – which includes the isobaricity assumption (more details given later) – i.e., \(p\approx\operatorname{const}\). Thus, expressions (20 )–(27 ) yield\[\begin{gather} \mathbf{J}_{1}\approx-D\frac{\left( m_{1}n_{1}+m_{2}n_{2}\right) \left( n_{2}\mathbf{\nabla}n_{1}-n_{1}\mathbf{\nabla}n_{2}\right) }{\left( n_{1}+n_{2}\right) ^{2}}\\ +B\frac{\mathbf{\nabla}n_{1}+\mathbf{\nabla}n_{2}}{n_{1}+n_{2}},\label{3469} \end{gather}\tag{28}\] \[\begin{gather} \mathbf{J}_{2}\approx-D\frac{\left( m_{1}n_{1}+m_{2}n_{2}\right) \left( n_{1}\mathbf{\nabla}n_{2}-n_{2}\mathbf{\nabla}n_{1}\right) }{\left( n_{1}+n_{2}\right) ^{2}}\\ -B\frac{\mathbf{\nabla}n_{1}+\mathbf{\nabla}n_{2}}{n_{1}+n_{2}},\label{34610} \end{gather}\tag{29}\] \[\mathbf{Q}\approx-B\frac{T\left( n_{1}\mathbf{\nabla}n_{2}-n_{2}\mathbf{\nabla}n_{1}\right) }{m_{1}n_{1}m_{2}n_{2}\left( n_{1}+n_{2}\right) ^{2}}-\kappa\frac{\mathbf{\nabla}n_{1}+\mathbf{\nabla}n_{2}}{\left( n_{1}+n_{2}\right) ^{2}}.\label{34611}\tag{30}\] In the next subsection, these expressions will be compared to their BGK counterparts. This is, generally, how the latter could be calibrated, so that its coefficients are related to the measured values of \(D\), \(B\), and \(\kappa\) of the gas mixture under consideration.

3.2 The multispecies BGK model↩︎

To derive the hydrodynamic approximation of a kinetic model, one should assume that the spatial scale of the solution exceeds the length \(l\) of the free path, and the solution’s temporal scale exceeds \(l/v\) where \(v\) is the mean velocity. Mathematically, these assumptions amount to ‘stretching’ the coordinates and time – i.e., replacing\[\frac{\partial}{\partial t}\rightarrow\varepsilon\frac{\partial}{\partial t},\qquad\mathbf{\nabla}\rightarrow\varepsilon\mathbf{\nabla}, \label{34612}\tag{31}\] where \(\varepsilon\) is a small parameter. One should then assume that the distribution function is nearly Maxwellian, with the velocity \(\mathbf{V}\) and temperature \(T\) being the same for all the species, i.e., \[\begin{gather} f_{i}=n_{i}\left( \frac{m_{i}}{2\pi T}\right) ^{3/2}\exp\left( -\frac{m_{i}\left\vert \mathbf{v}-\mathbf{V}\right\vert ^{2}}{2T}\right) \\ +\varepsilon f_{i}^{(1)}+\mathcal{O}(\varepsilon^{2}). \label{34613} \end{gather}\tag{32}\] In application to a BGK-type model, the hydrodynamic approximation was considered in Ref. [9], who also let\[\begin{align} \mathbf{V}_{i} & =\mathbf{V}+\varepsilon\mathbf{V}_{i}^{(1)}+\mathcal{O}(\varepsilon^{2}),\tag{33}\\ T_{i} & =T+\varepsilon T_{i}^{(1)}+\mathcal{O}(\varepsilon^{2}), \tag{34} \end{align}\] while leaving \(n_{i}\) nonexpanded. Substituting (32 )–(34 ) into the rescaled versions of Eqs. (4 )–(5 ), one can find \(f_{i}^{(1)}\) and then use Eqs. (2 )–(3 ) to find \(\mathbf{V}_{i}^{(1)}\) and \(T_{i}^{(1)}\), while Eq. (1 ) for \(n_{i}\) does not seem to be needed. Such a non-straightforward procedure was chosen in Ref. [9] because of the highly-nonlinear structure of the BGKmodel, making the straightforward calculation of higher-order corrections, such as the transport fluxes, cumbersome.

In the present paper, a slightly different approach is employed, where the transport coefficients are calculated under the diffusion approximation instead of the hydrodynamic one. The difference between the two approximations is two-fold. Firstly, the diffusion flow is slow, so that scaling (31 ) should be replaced with\[\frac{\partial}{\partial t}\rightarrow\varepsilon^{2}\frac{\partial}{\partial t},\qquad\mathbf{\nabla}\rightarrow\varepsilon\mathbf{\nabla.} \label{34616}\tag{35}\] Secondly, the diffusion flow is weak – so that expansions (32 )–(33 ) should be replaced with\[\begin{gather} f_{i}=n_{i}^{(0)}\left( \frac{m_{i}}{2\pi T}\right) ^{3/2}\exp\left( -\frac{m_{i}\left\vert \mathbf{v}\right\vert ^{2}}{2T}\right) \\ +\varepsilon f_{i}^{(1)}+\mathcal{O}(\varepsilon^{2}). \label{34617} \end{gather}\tag{36}\] \[\mathbf{V}_{i}=\varepsilon\mathbf{V}_{i}^{(1)}+\mathcal{O}(\varepsilon ^{2}),\qquad T_{i}=T+\varepsilon T_{i}^{(1)}+\mathcal{O}(\varepsilon^{2}), \label{34618}\tag{37}\] and the density should also be expanded, \[n_{i}=n_{i}^{(0)}+\varepsilon n_{i}^{(1)}+\mathcal{O}(\varepsilon^{2}). \label{34619}\tag{38}\] Under such an approximation, the transport fluxes emerge from the leading order of the expansion.

Having rescaled the BGK equations (4 )–(5 ) according to (35 ), one should substitute into them expansions (36 )–(38 ). \(\mathbf{V}_{ij}\) and \(T_{ij}\) should also be expanded [similarly to how \(\mathbf{V}_{i}\) and \(T_{i}\) are expanded in (37 )], as well as all the coefficients,

where\[M_{i}^{(0)}=n_{i}^{(0)}\left( \frac{m_{i}}{2\pi T}\right) ^{3/2}\exp\left( -\frac{m_{i}\left\vert \mathbf{v}\right\vert ^{2}}{2T}\right) .\] Substituting these expressions into the leading order of Eqs. (1 )–(3 ), one can verify that Eq. (1 ) is satisfied identically, and Eqs. (2 )–(3 ), (9 )–(12 ) yield\[\mathbf{\nabla}\left( n_{1}^{(0)}T^{(0)}+n_{2}^{(0)}T^{(0)}\right) =0, \label{34620}\tag{39}\] \[\begin{align} \mathbf{V}_{1}^{(1)} & =\mathbf{V}^{(1)}-\frac{\mathbf{\nabla}\left( n_{1}^{(0)}T^{(0)}\right) }{2\beta^{(0)}n_{1}^{(0)}n_{2}^{(0)}},\tag{40}\\ \mathbf{V}_{2}^{(1)} & =\mathbf{V}^{(1)}-\frac{\mathbf{\nabla}\left( n_{2}^{(0)}T^{(0)}\right) }{2\beta^{(0)}n_{1}^{(0)}n_{2}^{(0)}}, \tag{41} \end{align}\] \[T_{12}^{(1)}=T_{1}^{(1)}=T_{21}^{(1)}=T_{2}^{(1)}=T^{(1)},\] where \(\mathbf{V}^{(1)}(\mathbf{r},t)\) and \(T^{(1)}(\mathbf{r},t)\) are undetermined functions (neither will appear in the final expressions for the fluxes). Note that \(\beta_{1}\) and \(\beta_{2}\) which appear in the original set have been expressed through \(\beta\) using (15 ).

Physically, Eq. (39 ) reflects the isobaric nature of diffusion and heat conduction (the same result follows from the standard hydrodynamic equations or any other model). Introducing the leading-order pressure \(p^{(0)}\) (which may depend only on \(t\)), one can rewrite (39 ) in the form\[T^{(0)}=\frac{p^{(0)}}{n_{1}^{(0)}+n_{2}^{(0)}}.\label{34624}\tag{42}\] Next, introduce the mass-averaged velocity,\[\mathbf{\bar{V}}=\frac{m_{1}n_{1}\mathbf{V}_{1}+m_{2}n_{2}\mathbf{V}_{2}}{m_{1}n_{1}+m_{2}n_{2}},\] and the diffusion flux of the \(i\)-th species,\[\mathbf{J}_{i}=m_{i}n_{i}\left( \mathbf{V}_{i}-\mathbf{\bar{V}}\right) .\] Using (40 )–(41 ), one can calculate \(\mathbf{J}_{i}\) and then use (42 ) to eventually obtain

Comparing this expression to its ‘correct’ counterpart (30 ) with \(B=0\), one can see that the two results coincide only in the limit \(\beta\rightarrow\infty\) which makes the whole diffusive flux equal zero, not only its thermodiffusive part.

One way or another, no such value of the parameter \(\beta\) exists that makes the BGK fluxes satisfy the Onsager reciprocal relation – neither for the general case nor for a fluid with zero thermodiffusivity.

4 Concluding remarks↩︎

It should be emphasized that the results of the present paper apply to some, but not all, of the existing BGK-type models. Apart from the model examined above, they apply to that of Refs. [9], [10], which consists of Eqs. (1 )–(10 ), but with Eqs. (11 )–(12 ) replaced with\[\begin{align} T_{12} & =T_{1}+\alpha_{1}\left( T_{2}-T_{1}\right) +\gamma_{1}\left( \left\vert \mathbf{V}_{1}\right\vert ^{2}-\left\vert \mathbf{V}_{2}\right\vert ^{2}\right) ,\\ T_{21} & =T_{2}+\alpha_{2}\left( T_{1}-T_{2}\right) +\gamma_{2}\left( \left\vert \mathbf{V}_{2}\right\vert ^{2}-\left\vert \mathbf{V}_{1}\right\vert ^{2}\right) . \end{align}\] Even though these expressions differ from their counterparts examined here, the expressions for \(\mathbf{V}_{12}\) and \(\mathbf{V}_{21}\) are still the same, and this is enough for noncompliance with the Onsager relations. As for models where the collision frequencies \(\nu_{ij}\) depend on the molecular velocity (e.g., [14], [16]), those need to be tested separately. The present results do not cover them.

Note also that, even though the models examined in this paper do not formally satisfy the Onsager relations, they satisfy them asymptotically in the limit \[\frac{m_{1}}{m_{2}}\rightarrow0.\label{34626}\tag{44}\] To understand why, observe that the ratio of the first to second terms of heat flux (43 ) is proportional to \(m_{1}/m_{2}\) – hence, condition (44 ) allows one to neglect the first term. After that, expression (43 ) matches the standard heat flux expression (21 ) with \(C_{i}=0\), and so the corresponding Onsager relation holds. Note also that asymptotic limit (44 ) is important physically, as it describes ionized plasma (where the mass of electrons is indeed much smaller than that of ions). The asymptotic compliance with the Onsager relations occurs also in the limit \(m_{1}/m_{2}\rightarrow1\), in which case the first term in expression (43 ) vanishes.

One should not assume, however, that BGK-type models cannot satisfy the Onsager relations exactly. The model proposed in Ref. [24], for example, does have the correct transport properties – and also satisfies the so-called indifferentiability principle (i.e., if the molecules of the species have identical mechanical parameters, the distribution function of the mixture satisfies the single-species BGK equation). Unfortunately, this model does not seem to comply with the H theorem, as pointed out in Ref. [9].

Overall, a ‘perfect’ multispecies kinetic model should satisfy the following requirements:

1. conservation of mass, momentum, and energy;

2. H theorem;

3. indifferentiability principle;

4. positivity of the temperature and concentration;

5. ability to represent fluids with arbitrary values of the Prandtl and Schmidt numbers, and an arbitrary ratio of the bulk and shear viscosities;

6. Onsager reciprocal relations.

So far, none of the existing BGK-type models has been shown to comply with all of the above requirements (see, for example, the review sections of Refs. [13], [17]). This does not mean, however, that a fully compliant model does not exist in principle, and so one should hope that such will be developed in the future.

Finally, note that requirements 5–6 of the above list are particularly important for end users – i.e., researchers who need a practical tool to work with applications (like the present author, who looks for a tool to model evaporation of water into air). These requirements allow one to calibrate the model, so that its transport properties match those of the fluid under consideration.

References↩︎