2.1. Composition and equation of state of ideal gas mixture

We begin with introducing some nomenclature and notation. Let us consider a mixture composed of $M$ ideal gases. The composition is described by the species densities $\rho _a$, $a=1,\ldots , M$, while the mixture density is

(2.1)\begin{equation} \rho=\sum_{a=1}^{M}\rho_a. \end{equation}

Equivalently, the mixture composition is defined by the mixture density $\rho$ and the $M-1$ independent mass fractions $Y_a$,

(2.2a,b)\begin{equation} Y_a=\frac{\rho_a}{\rho},\quad \sum_{a=1}^{M} Y_a = 1. \end{equation}

With the molar mass of the component $m_a$, the mean molar mass $m$ depends on the composition,

(2.3)\begin{equation} \frac{1}{{m}}=\sum_{a=1}^M\frac{Y_a}{m_a}.\end{equation}

The equation of state provides a relation between the pressure $P$, the temperature $T$ and the composition,

(2.4)\begin{equation} P=\rho R T.\end{equation}

Here, $R$ is the specific gas constant that contains the information about the composition of the gas by way of the mean molar mass $m$,

(2.5)\begin{equation} R=\frac{R_U}{m},\end{equation}

where $R_U\approx 8.314\ \textrm {kJ}\ \textrm {K}^{-1} \boldsymbol {\cdot } \textrm {kmol}^{-1}$ is the universal gas constant. Thus, for a mixture of ideal gases, the specific gas constant $R$ is a function of local composition and changes in space and time. The pressure of an individual component $P_a$ is related to the pressure of the mixture $P$ through Dalton's law of partial pressures as

(2.6)\begin{equation} P_a=X_a P,\end{equation}

where the mole fraction of a component $X_a$ is related to its mass fraction $Y_a$ as

(2.7a,b)\begin{equation} X_a=\left(\frac{m}{m_a}\right)Y_a, \quad \sum_{a=1}^{M} X_a = 1.\end{equation}

A consequence of Dalton's law of partial pressure is that $\sum _{a=1}^{M}P_a=P$. Combined with the equation of state (2.4), the partial pressure $P_a$ takes the form

(2.8)\begin{equation} P_a=\rho_a R_a T,\end{equation}

where $R_a$ is the specific gas constant of the component,

(2.9)\begin{equation} R_a=\frac{R_U}{m_a}. \end{equation}

With these thermodynamic relations in mind, we proceed to setting up the kinetic equations that recover the Stefan–Maxwell diffusion in the macroscopic limit.

2.2. Kinetic equation for the species

In this section we set-up kinetic equations which recover the Stefan–Maxwell diffusion model for an $M$-component ideal gas mixture. Each component is described by a set of populations $f_{ai}$, $a=1,\ldots ,M$, corresponding to the discrete velocities $\boldsymbol {c}_i$, $i=0,\ldots , Q-1$. The proposed kinetic equation for each of the species $a$ is written as

(2.10)\begin{equation} \partial_t \,f_{ai} + \boldsymbol{c}_{i} \boldsymbol{\cdot}\boldsymbol{\nabla} f_{ai} = \sum_{b\!{\ne} a}^M \frac{1}{\theta_{ab}} \left[ \left( \frac{f_{ai}^{eq}-f_{ai}}{\rho_a} \right) - \left( \frac{f_{bi}^{eq}-f^*_{bi}}{\rho_b} \right) \right].\end{equation}

Here, $\rho _a$ is the density of the component $a$, which is defined as the zeroth moment of populations $f_{ai}$,

(2.11)\begin{equation} \rho_a=\sum_{i=0}^{Q-1}f_{ai}.\end{equation}

Furthermore, a symmetric set of relaxation parameters $\theta _{ab}=\theta _{ba}$ shall be related to the binary diffusion coefficients in the following. The equilibrium $f_{ai}^{eq}$ and the quasi-equilibrium $f^*_{ai}$ populations will be fully defined in § 2.4. Here, we only need to specify the conditions for the low-order moments thereof. To that end, let us introduce the partial momenta $\rho _a\boldsymbol {u}_a$ as first moments of the species’ populations,

(2.12)\begin{equation} \rho_a\boldsymbol{u}_a=\sum_{i=0}^{Q-1} f_{ai}\boldsymbol{c}_i.\end{equation}

The quasi-equilibrium populations $f_{ai}^*$ must satisfy the following set of constraints,

(2.13)\begin{gather} \sum_{i=0}^{Q-1} f_{ai}^{*} = \rho_a, \end{gather}

(2.14)\begin{gather}\sum_{i=0}^{Q-1} f_{ai}^{*}\boldsymbol{c}_i = \rho_a \boldsymbol{u}_a. \end{gather}

The momenta of the components sum up to the mixture momentum,

(2.15)\begin{equation} \sum_{a=1}^{M} \rho_a \boldsymbol{u}_a = \rho \boldsymbol{u}.\end{equation}

The equilibrium populations $f_{ai}^{eq}$ have to verify the following set of constraints:

(2.16)\begin{gather} \sum_{i=0}^{Q-1} f_{ai}^{eq} = \rho_a, \end{gather}

(2.17)\begin{gather}\sum_{i=0}^{Q-1} f_{ai}^{eq}\boldsymbol{c}_i = \rho_a \boldsymbol{u}, \end{gather}

(2.18)\begin{gather}\sum_{i=0}^{Q-1} f_{ai}^{eq}\boldsymbol{c}_i\otimes \boldsymbol{c}_i = P_a\boldsymbol{I} +\rho_a\boldsymbol{u}\otimes \boldsymbol{u}. \end{gather}

In (2.18) the partial pressure $P_a$ (2.8) depends on the temperature $T$, which is obtained from the mixture kinetic equations of § 3. Finally, the quasi-equilibrium distribution must match the equilibrium if the species velocity equals the velocity of the mixture, $\boldsymbol {u}_a=\boldsymbol {u}$:

(2.19)\begin{equation} f^*_{ai}(\boldsymbol{u}_a)|_{\boldsymbol{u}_a=\boldsymbol{u}}=f^{eq}_{ai}(\boldsymbol{u}). \end{equation}

Some comments are in order.

1. (i) The $M$-component kinetic system satisfies $M+D$ conservation laws, where $D$ is the space dimension: the densities $\rho _1, \ldots , \rho _M$ and the vector of fluid momentum $\rho \boldsymbol {u}$ are locally conserved fields.

2. (ii) Thanks to the matching condition (2.19), the relaxation term on the right-hand side of (2.10) vanishes only at the equilibrium.

We now proceed with the identification of the relaxation parameters $\theta _{ab}$ in terms of the Stefan–Maxwell binary diffusion coefficients.

2.3. Hydrodynamic limit of kinetic equations for the species

Evaluation of the zeroth and first moments of the kinetic equation (2.10) results in the balance equations for the species densities and species velocities,

(2.20)\begin{gather} \partial_t\rho_a={-}\boldsymbol{\nabla}\cdot(\rho_a \boldsymbol{u}_a), \end{gather}

(2.21)\begin{gather}\rho_a\partial_t \boldsymbol{u}_a=\boldsymbol{u}_a\boldsymbol{\nabla}\cdot(\rho_a\boldsymbol{u}_a)- \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{P}_a+\sum_{b\!{\ne} a}^M \frac{1}{\theta_{ab}} \left( \boldsymbol{u}_{b} - \boldsymbol{u}_a \right)\!. \end{gather}

Here, $\boldsymbol {P}_a$ is the partial pressure tensor,

(2.22)\begin{equation} \boldsymbol{P}_a=\sum_{i=0}^{Q-1} f_{ai}\boldsymbol{c}_i\otimes \boldsymbol{c}_i. \end{equation}

Upon summation over the components in (2.20), we arrive at the continuity equation for the mixture density,

(2.23)\begin{equation} \partial_t\rho={-}\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u}), \end{equation}

while the summation over components in (2.21) results in the mixture momentum balance,

(2.24)\begin{equation} \partial_t(\rho\boldsymbol{u})={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{P}, \end{equation}

where $\boldsymbol {P}$ is the mixture pressure tensor,

(2.25)\begin{equation} \boldsymbol{P}=\sum_{a=1}^{M}\boldsymbol{P}_a. \end{equation}

The low-order closure relation for the species balance equation (2.20) is established by considering a perturbation around the equilibrium,

(2.26)\begin{equation} \boldsymbol{u}_a=\boldsymbol{u}+\delta\boldsymbol{u}_a,\end{equation}

where the perturbation $\delta \boldsymbol {u}_a$ satisfies the consistency condition,

(2.27)\begin{equation} \sum_{a=1}^{M}\rho_a\delta\boldsymbol{u}_a=0. \end{equation}

To first order, upon substitution into (2.21), we get the constitutive relation for the diffusion velocity ${\delta }\boldsymbol {u}_a$,

(2.28)\begin{equation} \rho_a\partial_t\boldsymbol{u}-\boldsymbol{u}\boldsymbol{\nabla}\cdot(\rho_a\boldsymbol{u})+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{P}_a^{eq}=\sum_{b\!{\ne} a}^M \frac{1}{\theta_{ab}} \left({\delta}\boldsymbol{u}_{b} - {\delta}\boldsymbol{u}_a \right)\!. \end{equation}

Upon summation over the species, and by taking into account Dalton's law in the equilibrium pressure tensor (2.18), the compressible Euler equation for the flow velocity is established,

(2.29)\begin{equation} \partial_t\boldsymbol{u}={-}\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}-\frac{1}{\rho}\boldsymbol{\nabla}{P}. \end{equation}

By elimination of the time derivative in (2.28), we get the constitutive relation as

(2.30)\begin{equation} P\boldsymbol{\nabla} X_a+(X_a-Y_a)\boldsymbol{\nabla} P=\sum_{b\!{\ne} a}^M \frac{1}{\theta_{ab}} ({\delta}\boldsymbol{u}_{b} - {\delta}\boldsymbol{u}_a). \end{equation}

Constitutive relation (2.30) becomes the Stefan–Maxwell diffusion equation once the relaxation parameters $\theta _{ab}$ are identified in terms of the binary diffusion coefficients $\mathcal {D}_{ab}$ as

(2.31)\begin{equation} \theta_{ab}=\frac{\mathcal{D}_{ab}}{PX_aX_b}.\end{equation}

Summarizing, kinetic equations for the species (2.10) recover the Stefan–Maxwell law of diffusion in the hydrodynamic limit, with both the diffusion due to non-uniformity of the species molar concentration and the barodiffusion taken into account. The present model does not include thermodiffusion as it should be expected by the simplicity of the relaxation term. We comment that the above derivation of (2.30) assumes the validity of the equation of state. The latter, in turn, depends on the temperature derived from the mixture energy equation, and which shall be introduced in § 3. We now proceed with finalizing the continuous time–space kinetic equations by identifying the equilibrium and the quasi-equilibrium populations.

2.4. Realization on the standard lattice

The above kinetic model is realized on the standard three-dimensional $D3Q27$ lattice, where $D=3$ stands for three dimensions and $Q=27$ is the number of discrete velocities:

(2.32a,b)\begin{equation} \boldsymbol{c}_i=(c_{ix},c_{iy},c_{iz}),\quad c_{i\alpha}\in\{{-}1,0,1\}. \end{equation}

Following Karlin & Asinari (Reference Karlin and Asinari2010) we define a triplet of functions in two variables, $\xi$ and $\zeta >0$,

(2.33a–c)\begin{align} \varPsi_{0}(\xi,\zeta) = 1 - (\xi^2 + \zeta), \quad \varPsi_{1}(\xi,\zeta) = \frac{\xi + (\xi^2 + \zeta)}{2},\quad \varPsi_{{-}1}(\xi,\zeta) = \frac{-\xi + (\xi^2 + \zeta)}{2}. \end{align}

For a vector $\boldsymbol {\xi }=(\xi _x,\xi _y,\xi _z)$, we consider a product form associated with the discrete velocities $\boldsymbol {c}_i$ (2.32a,b),

(2.34)\begin{equation} \varPsi_i(\boldsymbol{\xi},\zeta)= \varPsi_{c_{ix}}(\xi_x,\zeta) \varPsi_{c_{iy}}(\xi_y,\zeta) \varPsi_{c_{iz}}(\xi_z,\zeta). \end{equation}

The equilibrium $f_{ai}^{eq}$ and the quasi-equilibrium $f_{ai}^*$ are represented with the product form (2.34) by choosing $\boldsymbol {\xi }=\boldsymbol {u}$ or $\boldsymbol {\xi }=\boldsymbol {u}_a$, respectively, and by assigning $\zeta =R_aT$ in both cases:

(2.35)\begin{gather} f_{ai}^{eq}(\rho_a,\boldsymbol{u},R_aT)= \rho_a\varPsi_{c_{ix}} (u_x,R_aT ) \varPsi_{c_{iy}} (u_y,R_aT ) \varPsi_{c_{iz}} (u_z,R_aT ), \end{gather}

(2.36)\begin{gather}f_{ai}^{*}(\rho_a,\boldsymbol{u}_a,R_aT)= \rho_a\varPsi_{c_{ix}} (u_{ax},R_aT ) \varPsi_{c_{iy}} (u_{ay},R_aT ) \varPsi_{c_{iz}} (u_{az},R_aT ). \end{gather}

One can readily verify that, the equilibrium (2.35) and the quasi-equilibrium (2.36) satisfy all the constraints put forward in § 2.2. We now proceed with the lattice Boltzmann discretization of the kinetic equations (2.10).

2.5. Lattice Boltzmann equation for the species

2.5.1. Kinetic equations in the relaxation form

With the mass diffusivity (2.31), the kinetic equation (2.10) is written as

(2.37)\begin{equation} \partial_t\, f_{ai} + \boldsymbol{c}_{i}\boldsymbol{\cdot}\boldsymbol{\nabla} f_{ai} = \sum_{b\!{\ne} a}^M \frac{P X_a X_b}{\mathcal{D}_{ab}} \left[ \left( \frac{f_{ai}^{eq}-f_{ai}}{\rho_a} \right) - \left( \frac{f_{bi}^{eq}-f^*_{bi}}{\rho_b} \right) \right]. \end{equation}

It can readily be seen that, in its present form, (2.37) is not well suited for numerical implementation. Indeed, in an actual problem, the density of some species can be small or even vanishing if a particular gas component is absent at some location. This is an inconvenience rather than a failure since vanishing of the density is compensated by the simultaneously vanishing molar fraction in the product $X_aX_b$. Hence, we first transform (2.37) in order to eliminate this artifact. Substituting the equation of state (2.4) into (2.37), we get

(2.38)\begin{equation} \partial_t\, f_{ai} + \boldsymbol{c}_{i}\boldsymbol{\cdot}\boldsymbol{\nabla} f_{ai} = \sum_{b\!{\ne} a}^M \left(\frac{ m }{m_a m_b}\right)\left(\frac{ R_UT }{\mathcal{D}_{ab}}\right) [ Y_b (\,f_{ai}^{eq}-f_{ai} ) - Y_a (\,f_{bi}^{eq}-f^*_{bi} ) ]. \end{equation}

Equation (2.38) is equivalent to (2.37) and does not suffer from a spurious division by a vanishing density. Furthermore, it proves convenient to recast (2.38) in a relaxation form. To that end, let us define characteristic times $\tau _{ab}=\tau _{ba}$,

(2.39)\begin{equation} \frac{1}{\tau_{ab}} = \left(\frac{R_UT}{\mathcal{D}_{ab}}\right)\left(\frac{m}{m_a m_b}\right), \end{equation}

and let us introduce the relaxation times $\tau _a$,

(2.40)\begin{equation} \frac{1}{\tau_a} = \sum_{b\!{\ne} a}^M \frac{Y_b}{\tau_{ab}} =R_a T\left(\sum_{b\!{\ne} a}^M \frac{X_b}{\mathcal{D}_{ab}}\right).\end{equation}

Finally, let us introduce a shorthand notation,

(2.41)\begin{equation} F_{ai} = Y_a \sum_{b\!{\ne} a}^M \frac{1}{\tau_{ab}} (\,f_{bi}^{eq}-f_{bi}^* ). \end{equation}

With these definitions, the kinetic equation (2.38) can be rearranged as follows:

(2.42)\begin{equation} \partial_t\, f_{ai} + \boldsymbol{c}_{i}\boldsymbol{\cdot}\boldsymbol{\nabla} f_{ai} = \frac{1}{\tau_a}(\,f_{ai}^{eq}-f_{ai} ) - F_{ai}.\end{equation}

The species kinetic equations (2.42) are now cast into the relaxation form, familiar from previous lattice Boltzmann models. The right-hand side comprises the conventional relaxation term and a source term (2.41). The latter depends on the populations only through the local densities, momenta and the temperature, as prescribed by the local equilibrium and quasi-equilibrium populations (2.35) and (2.36). Hence, kinetic equation (2.42) is amenable to a lattice Boltzmann discretization in time and space.

2.5.2. Derivation of the lattice Boltzmann equation

Following a procedure first introduced by He et al. (Reference He, Chen and Doolen1998), we integrate (2.42) along the characteristics and apply the trapezoidal rule on the right-hand side to obtain

(2.43)\begin{align} &f_{ai}(\boldsymbol{x}+\boldsymbol{c}_i \delta t, t+ \delta t) - f_{ai}(\boldsymbol{x},t) = \frac{\delta t}{2 \tau_a}[f_{ai}^{eq}(\boldsymbol{x}+\boldsymbol{c}_i \delta t, t+ \delta t) - f_{ai}(\boldsymbol{x}+\boldsymbol{c}_i \delta t, t+ \delta t)] \nonumber\\ &\quad + \frac{\delta t}{2 \tau_a}[f_{ai}^{eq}(\boldsymbol{x},t) - f_{ai}(\boldsymbol{x},t)] - \frac{\delta t}{2} F_{ai}(\boldsymbol{x}+\boldsymbol{c}_i \delta t, t+ \delta t)- \frac{\delta t}{2} F_{ai}(\boldsymbol{x}, t). \end{align}

Next, we introduce transformed populations $k_{ai}$,

(2.44)\begin{equation} f_{ai} = k_{ai} + \frac{\delta t}{2 \tau_a} (\,f_{ai}^{eq}-f_{ai}) - \frac{\delta t}{2} F_{ai}. \end{equation}

Let us evaluate the pertinent moments of the transform (2.44). Summation over the discrete velocities gives

(2.45)\begin{equation} \rho_a(\,f) = \rho_a(k),\end{equation}

where we have specified that the density $\rho _a(\,f)$ on the left-hand side is defined using the original populations $f_{ai}$ while the density $\rho _a(k)$ is defined by the zeroth moment of the $k$-populations,

(2.46)\begin{equation} \rho_a(k)=\sum_{i=0}^{Q-1}k_{ai}. \end{equation}

Thus, the species densities do not alter under the populations transformation, and the specification can be dropped: $\rho _a=\rho _a(\,f)=\rho _a(k)$. On the other hand, evaluating the first moment of (2.44) gives

(2.47)\begin{equation} \rho_a \boldsymbol{u}_{a}(\,f) \left( 1+ \frac{\delta t}{2 \tau_a}\right) - \frac{\delta t}{2} Y_a \sum_{b\!{\ne} a}^{M} \frac{1}{\tau_{ab}} \rho_b\boldsymbol{u}_{b }(\,f)=\rho_a \boldsymbol{u}_{a}(k), \end{equation}

where the species velocity $\boldsymbol {u}_a(k)$ is defined by the $k$-populations in a conventional way,

(2.48)\begin{equation} \rho_a\boldsymbol{u}_a(k)=\sum_{i=0}^{Q-1}k_{ai}\boldsymbol{c}_i. \end{equation}

Summation over the species in (2.47) shows that the momentum $\rho \boldsymbol {u}$ is also an invariant of the transform (2.44):

(2.49)\begin{equation} \rho\boldsymbol{u}(\,f)=\rho\boldsymbol{u}(k). \end{equation}

Since the term $F_{ai}$ vanishes at equilibrium, and also thanks to the invariance of the local conservation (2.45) and (2.49), the equilibrium is the fixed point of the map (2.44):

(2.50)\begin{equation} f_{ai}^{eq}(\rho_a,\boldsymbol{u},R_aT)=k_{ai}^{eq}(\rho,\boldsymbol{u},R_aT). \end{equation}

Substituting (2.44) into (2.43), and introducing the parameters $\beta _a\in [0,1]$,

(2.51)\begin{equation} \beta_a=\frac{\delta t}{2 \tau_a + \delta t}, \end{equation}

we get

(2.52)\begin{equation} k_{ai}(\boldsymbol{x}+\boldsymbol{c}_i \delta t, t+ \delta t) = k_{ai}(\boldsymbol{x},t)+ 2 \beta_a [k_{ai}^{eq}(\boldsymbol{x},t) - k_{ai}(\boldsymbol{x},t)] + \delta t (\beta_a-1) F_{ai}(\boldsymbol{x}, t).\end{equation}

The last term in (2.52) is spelled out as follows. The quasi-equilibrium population $f_{ai}^*$ in the expression $F_{ai}$ (2.41) depends on the species velocity $\boldsymbol {u}_a(\,f)$. The latter, unlike the mixture velocity $\boldsymbol {u}(\,f)$ (2.49), is not invariant under the transform to the $k$-populations. Rather, $\boldsymbol {u}_a(\,f)$ has to be evaluated using the linear relation (2.47) in terms of $\boldsymbol {u}_b(k)$ by solving a $M\times M$ linear system for each of the spatial components. In our realization we used Gauss elimination.

2.5.3. Summary of the lattice Boltzmann equation for the Stefan–Maxwell diffusion

For convenience, we summarize the lattice Boltzmann equation for the species. We return to a more conventional notation and rename $k_{ai}$ to $f_{ai}$ in (2.52),

(2.53)\begin{equation} f_{ai}(\boldsymbol{x}+\boldsymbol{c}_i \delta t, t+ \delta t) - f_{ai}(\boldsymbol{x},t)= 2 \beta_a [f_{ai}^{eq} - f_{ai}] + \delta t (\beta_a-1) F_{ai}.\end{equation}

The last term is written as

(2.54)\begin{equation} F_{ai} = Y_a \sum_{b\!{\ne} a}^M \left(\frac{R_UT}{\mathcal{D}_{ab}}\right)\left(\frac{m}{m_a m_b}\right) [\, f_{bi}^{eq}(\rho_b,\boldsymbol{u},R_bT)-f_{bi}^*(\rho_b,\boldsymbol{u}+\boldsymbol{V}_b,R_bT) ],\end{equation}

where we have introduced the transformed diffusion velocities of the components, $\boldsymbol {V}_a$, $a=1,\ldots ,M$. The latter are defined by (2.47) which can be recast as follows:

(2.55)\begin{equation} \rho_a \boldsymbol{V}_{a}-\left(\frac{\delta t}{2}\right)\sum_{b\!{\ne} a}^{M}\frac{PX_aX_b}{\mathcal{D}_{ab}}[\boldsymbol{V}_{b }-\boldsymbol{V}_{a }]=\rho_a (\boldsymbol{u}_{a}-\boldsymbol{u}). \end{equation}

The field $\rho _a\boldsymbol {u}_a$ on the right-hand side of (2.55) is defined by the moment relation as before,

(2.56)\begin{equation} \rho_a\boldsymbol{u}_a=\sum_{i=0}^{Q-1}f_{ai}\boldsymbol{c}_i. \end{equation}

The $M+D$ independent conservation laws of the $M$-component system (2.53) correspond to the mass of each component and the momentum of the mixture. The corresponding locally conserved fields are the species densities $\rho _a$ and the momentum flux $\rho \boldsymbol {u}$,

(2.57)\begin{gather} \rho_a=\sum_{i=0}^{Q-1}f_{ai}, \end{gather}

(2.58)\begin{gather}\rho\boldsymbol{u}=\sum_{a=1}^M \rho_a\boldsymbol{u}_a=\sum_{a=1}^M\sum_{i=0}^{Q-1}f_{ai}\boldsymbol{c}_i. \end{gather}

The conservation law of the fluid mass is the implication of the mass conservation of each component. The density of the mixture $\rho$ is defined by the densities of the components,

(2.59)\begin{equation} \rho=\sum_{a=1}^M \rho_a=\sum_{a=1}^M\sum_{i=0}^{Q-1}f_{ai}, \end{equation}

while the flow velocity $\boldsymbol {u}$ is

(2.60)\begin{equation} \boldsymbol{u}=\frac{\rho\boldsymbol{u}}{\rho}=\frac{\displaystyle\sum_{a=1}^M\sum_{i=0}^{Q-1}f_{ai} \boldsymbol{c}_i}{\displaystyle\sum_{a=1}^M\sum_{i=0}^{Q-1}f_{ai}}. \end{equation}

Finally, we recall that the temperature dependence in the equilibrium and in the quasi-equilibrium populations has to be supplied by the energy equation to be discussed in the next section.

The Stefan–Maxwell diffusion relation is a source of a variety of approximate constitutive relations for the multicomponent diffusion where the former is substituted by a simpler, often explicit expression. Diffusion models such as the mixture averaged approximation are widely used; see, e.g., Giovangigli (Reference Giovangigli2012), Williams (Reference Williams1985), Poinsot & Veynante (Reference Poinsot and Veynante2005), Kee et al. (Reference Kee, Coltrin and Glarborg2003) and Bird et al. (Reference Bird, Stewart and Lightfoot2006). In appendix A we show how approximate constitutive relations are derived based on the kinetic system (2.10).

3.1. First law of thermodynamics for ideal gas mixture

For further references and notation, we open this section with a summary of the first law of thermodynamics for ideal gas mixtures. Since non-reactive mixtures are considered in the following, the energy of formation is not included. The caloric equation of state of a single-component ideal gas provides for the specific mole-based sensible internal energy of species $a$:

(3.1)\begin{equation} \bar{U}_a=\int_{T_0}^T\bar{C}_{a,v}(T)\, \textrm{d} T. \end{equation}

Here, $\bar {C}_{a,v}$ is the specific heat at constant volume. Thus, the sensible specific enthalpy reads as

(3.2)\begin{equation} \bar{H}_a=\int_{T_0}^T\bar{C}_{a,p}(T)\, \textrm{d} T, \end{equation}

where $\bar {C}_{a,p}$ is the specific heat at constant pressure defined by Mayer's relation,

(3.3)\begin{equation} \bar{C}_{a,p}-\bar{C}_{a,v}=R_{U}. \end{equation}

Proceeding from the mole basis onto the mass basis, the specific heats are defined relative to the molar mass,

(3.4)\begin{gather} {C}_{a,v}=\frac{\bar{C}_{a,v}}{m_a}, \end{gather}

(3.5)\begin{gather}{C}_{a,p}=\frac{\bar{C}_{a,p}}{m_a}, \end{gather}

while the mass-based specific sensible internal energy and enthalpy are

(3.6)\begin{gather} {U}_a=\int_{T_0}^T{C}_{a,v}(T)\, \textrm{d} T, \end{gather}

(3.7)\begin{gather}{H}_a=\int_{T_0}^T{C}_{a,p}(T)\, \textrm{d} T. \end{gather}

Finally, the Mayer relation in the mass basis reads as

(3.8)\begin{equation} {C}_{a,p}-{C}_{a,v}=R_{a}, \end{equation}

where the gas constant $R_a$ is defined by (2.9).

Switching to the case of a $M$-component mixture, the mixture internal energy $\rho U$ is defined on the mass basis as follows:

(3.9)\begin{equation} \rho U=\sum_{a=1}^M\rho_a U_a.\end{equation}

The specific mixture internal energy $U$ can be rewritten as

(3.10)\begin{equation} U=\sum_{a=1}^MY_a U_a=\sum_{a=1}^MY_a \int_{T_0}^T{C}_{a,v}\, \textrm{d} T=\int_{T_0}^T\left[\sum_{a=1}^MY_a{C}_{a,v}\right]\, \textrm{d} T=\int_{T_0}^T{C}_{v}\, \textrm{d} T, \end{equation}

where the specific heat at constant volume is the mass-averaged value over the composition,

(3.11)\begin{equation} {C}_{v}=\sum_{a=1}^MY_a{C}_{a,v}. \end{equation}

Similarly, the mixture enthalpy $\rho H$ is defined as

(3.12)\begin{equation} \rho H=\sum_{a=1}^M\rho_a H_a, \end{equation}

while the specific mixture enthalpy $H$ can be transformed in the manner of (3.10),

(3.13)\begin{equation} H=\sum_{a=1}^MY_a H_a=\sum_{a=1}^MY_a \int_{T_0}^T{C}_{a,p}\, \textrm{d} T=\int_{T_0}^T\left[\sum_{a=1}^MY_a{C}_{a,p}\right]\, \textrm{d} T=\int_{T_0}^T{C}_{p}\, \textrm{d} T. \end{equation}

The specific heat at constant pressure reads as

(3.14)\begin{equation} {C}_{p}=\sum_{a=1}^MY_a{C}_{a,p}, \end{equation}

while both the specific heats satisfy the Mayer relation,

(3.15)\begin{equation} {C}_{p}-{C}_{v}=R, \end{equation}

with the mixture gas constant $R$ defined by (2.5). In the following, we formulate the lattice Boltzmann equation for the mixture density, momentum and energy for a generic case of temperature-dependent specific heats of the components.

3.2. Two-population lattice Boltzmann equation for the mixture

Point of departure is a lattice Boltzmann model for a single-component ideal gas with variable Prandtl number and adiabatic exponent. To that end, several suitable single-component lattice Boltzmann models exist in the literature; here we mention compressible LBM by Frapolli et al. (Reference Frapolli, Chikatamarla and Karlin2015), Frapolli, Chikatamarla & Karlin (Reference Frapolli, Chikatamarla and Karlin2016a) and Saadat et al. (Reference Saadat, Bösch and Karlin2019). The common feature of these single-component models is the use of the double-population construction, the idea first introduced in the context of incompressible thermal convective LBM in the classical paper by He et al. (Reference He, Chen and Doolen1998) and further expanded in Guo et al. (Reference Guo, Zheng, Shi and Zhao2007), Karlin, Sichau & Chikatamarla (Reference Karlin, Sichau and Chikatamarla2013), Frapolli, Chikatamarla & Karlin (Reference Frapolli, Chikatamarla and Karlin2018). One set of populations, commonly quoted as $f$-populations, represents the density and momentum as the locally conserved fields of the corresponding $f$-LBM equation while another set, the $g$-populations, represents the energy as the local conservation of the $g$-LBM kinetics. The coupling between the $f$- and $g$-LBM equations is also well understood and enables the realization of an adjustable Prandtl number and adiabatic exponent. Various realizations differ by the choice of the discrete velocities of the $f$- and $g$-sets; in particular, the compressible LBM of Frapolli et al. (Reference Frapolli, Chikatamarla and Karlin2015) employs higher-order lattices with higher isotropy while the two-dimensional model developed in Saadat et al. (Reference Saadat, Bösch and Karlin2019) uses the standard lattice with correction terms to compensate for insufficient isotropy.

Whichever single-component double-population model is taken as the starting point for representing a multicomponent mixture, the central question is how to modify it. Note that, this question would not arise if one would follow the conventional approach by extending the already available $M$ species LBM equations of § 2 to represent the energy equation of the mixture. However, with the double-population approach, this would lead to $2\times M$ lattices since the lattice for each component would need to be doubled to represent the energy of that component. On the contrary, the mean-field approach pursuit here avoids the kinetic representation of partial energies, instead it addresses only the total energy of the mixture by a single $g$-set. This requires only $M+2$ lattices, $M$ for the species and two for the mixture momentum and energy.

In the following, we refer to the $f$-populations as the momentum lattice, and the $g$-populations as the energy lattice. For the momentum lattice, the locally conserved fields are the density and the momentum of the mixture,

(3.16)\begin{gather} \sum_{i=0}^{Q-1} f_i = \rho, \end{gather}

(3.17)\begin{gather}\sum_{i=0}^{Q-1} f_i \boldsymbol{c}_{i} = \rho \boldsymbol{u}. \end{gather}

For the energy lattice, the locally conserved field is the total energy of the mixture,

(3.18)\begin{equation} \sum_{i=0}^{Q-1} g_i = \rho E.\end{equation}

Here, the total energy $\rho E$ is the sum of the mixture internal energy $\rho U$ (3.10) and the kinetic energy $\rho u^2/2$,

(3.19)\begin{equation} \rho E=\rho U + {\frac{\rho u^2}{2}}. \end{equation}

Since the mixture internal energy (3.10) depends on the composition, it is the first instance where the species kinetic equations become coupled with the kinetic equations for the mixture. Conversely, the temperature of the mixture is computed from the integral equation,

(3.20)\begin{equation} \int_{T_0}^T C_{v}(T)\, \textrm{d} T=E-\frac{u^2}{2}. \end{equation}

The temperature evaluated from solving (3.20) is used as the input in the definition of the pressure $P$ in the equilibrium and quasi-equilibrium constraints of the species lattice Boltzmann system. This furnishes the input from the energy lattice into the species lattices.

We comment that, in the present section, the mixture density (3.16) and momentum (3.17) are defined self-consistently in the sense of $f$-populations of the momentum lattice. On the other hand, quantities carrying the same physical meaning were independently and also self-consistently defined earlier using the species populations, (2.59) and (2.58), respectively. Doubling of the conservation of the total mass and momentum is the feature of the intermediate steps of the construction during which the species subsystem and the mixture subsystem are set-up independently from one another. At the end of the construction, the doubling of the conservation shall be resolved through a coupling of both the species and the mixture subsystems in § 3.5.

The lattice Boltzmann equations for the momentum and for the energy lattice are patterned from the single-component developments,

(3.21)\begin{gather} f_i(\boldsymbol{x}+\boldsymbol{c}_i \delta t,t+\delta t)- f_i(\boldsymbol{x},t)= \omega (\,f_i^{eq} -f_i), \end{gather}

(3.22)\begin{gather}g_i(\boldsymbol{x}+\boldsymbol{c}_i \delta t,t+ \delta t) - g_i(\boldsymbol{x},t)= \omega_1 (g_i^{eq} -g_i) + (\omega - \omega_1) (g_i^{*} -g_i), \end{gather}

where relaxation parameters $\omega$ and $\omega _1$ shall be related to the viscosity and thermal conductivity in the following. We now proceed with specifying the constraints on the equilibrium populations $f_i^{eq}$ and $g_i^{eq}$, and the quasi-equilibrium $g_i^{*}$ in order that the system (3.21) and (3.22) recovers the momentum and energy equations of the mixture.

First, the equilibrium populations must satisfy the $D+2$ conservation laws,

(3.23)\begin{gather} \sum_{i=0}^{Q-1} f_i^{eq} = \rho, \end{gather}

(3.24)\begin{gather}\sum_{i=0}^{Q-1} f_i^{eq} \boldsymbol{c}_{i } = \rho \boldsymbol{u}, \end{gather}

(3.25)\begin{gather}\sum_{i=0}^{Q-1} g_i^{eq} = \rho E. \end{gather}

Second, the equilibrium pressure tensor $\boldsymbol {P}^{eq}$ and the tensor of equilibrium third-order moments $\boldsymbol {Q}^{eq}$ of the momentum lattice must verify the Maxwell–Boltzmann relations in order to recover the compressible flow momentum equation,

(3.26)\begin{gather} \boldsymbol{P}^{eq} =\sum_{i=0}^{Q-1} f_i^{eq} \boldsymbol{c}_i\otimes\boldsymbol{c}_i = P \boldsymbol{I}+\rho \boldsymbol{u}\otimes\boldsymbol{u}, \end{gather}

(3.27)\begin{gather}\boldsymbol{Q}^{eq} =\sum_{i=0}^{Q-1} f_i^{eq} \boldsymbol{c}_i\otimes\boldsymbol{c}_i\otimes\boldsymbol{c}_i = P\overline{\boldsymbol{u}\otimes\boldsymbol{I}} + \rho \boldsymbol{u}\otimes\boldsymbol{u}\otimes\boldsymbol{u}, \end{gather}

where the overline denotes symmetrization. Similarly, the equilibrium mixture energy flux $\boldsymbol {q}^{eq}$ and the second-order moment tensor $\boldsymbol {R}^{eq}$ pertinent to the energy lattice are

(3.28)\begin{gather} \boldsymbol{q}^{eq}= \sum_{i=0}^{Q-1} g_i^{eq} \boldsymbol{c}_{i} = \left(H+\frac{u^2}{2}\right)\rho\boldsymbol{u}, \end{gather}

(3.29)\begin{gather}\boldsymbol{R}^{eq}=\sum_{i=0}^{Q-1} g_i^{eq} \boldsymbol{c}_i\otimes\boldsymbol{c}_i = \left(H+\frac{u^2}{2}\right) \boldsymbol{P}^{eq} + P\boldsymbol{u}\otimes\boldsymbol{u}. \end{gather}

The mixture equation of state $P$ (2.4), the mixture gas constant $R$ (2.5) and the specific enthalpy of the mixture $H$ (3.13) entering the constraints (3.26)–(3.28), (3.28) and (3.29) depend linearly on the composition through the local mass fractions $Y_a$.

To that end, the constraints on the equilibrium populations of the mixture momentum and energy lattices is a straightforward extension of those of the single-component double-population LBM for compressible flows where the ideal gas equation of state, the internal energy and the enthalpy are merely replaced by their mixture-averaged counterparts. A major difference comes next with the constraints for the quasi-equilibrium. The zeroth-, first- and second-order moments of the quasi-equilibrium populations $g_i^{*}$, or the quasi-equilibrium energy $\rho E^{*}$, the energy flux $\boldsymbol {q}^{*}$ and the flux of the energy flux $\boldsymbol {R}^{*}$, respectively, have to satisfy the following relations:

(3.30)\begin{gather} \rho E^{*}=\sum_{i=0}^{Q-1} g_i^{*}=\rho E, \end{gather}

(3.31)\begin{gather}\boldsymbol{q}^{*} =\sum_{i=0}^{Q-1} g_i^{*} \boldsymbol{c}_{i} = \boldsymbol{q} - \boldsymbol{u} \cdot (\boldsymbol{P} - \boldsymbol{P}^{eq}) +\boldsymbol{q}^{diff}+\boldsymbol{q}^{corr}, \end{gather}

(3.32)\begin{gather}\boldsymbol{R}^{*}=\sum_{i=0}^{Q-1} g_i^{*} \boldsymbol{c}_{i}\otimes \boldsymbol{c}_i=\boldsymbol{R}^{eq}. \end{gather}

The first and third of these quasi-equilibrium constraints, (3.30) and (3.32), as well as the first and second terms in the quasi-equilibrium energy flux (3.31) are again the direct extension of the single-component LBM. Specifically, the two first terms in (3.31), comprising the energy flux $\boldsymbol {q}$ and the pressure tensor $\boldsymbol {P}$,

(3.33)\begin{gather} \boldsymbol{q}=\sum_{i=0}^{Q-1} g_i \boldsymbol{c}_{i}, \end{gather}

(3.34)\begin{gather}\boldsymbol{P}= \sum_{i=0}^{Q-1} f_i \boldsymbol{c}_{i}\otimes \boldsymbol{c}_{i}, \end{gather}

are needed to decouple the viscosity from thermal conductivity, and to maintain a variable Prandtl number, in both the single-component and multicomponent cases.

The remaining two terms in the quasi-equilibrium energy flux (3.31), $\boldsymbol {q}^{diff}$ and $\boldsymbol {q}^{corr}$, are specific to the multicomponent case and appear due to the mean-field approach to the energy representation. The interdiffusion energy flux $\boldsymbol {q}^{diff}$ reads as

(3.35)\begin{equation} \boldsymbol{q}^{diff}=\left(\frac{\omega_1}{\omega-\omega_1} \right) \rho\sum_{a=1}^{M}H_aY_a\boldsymbol{V}_a, \end{equation}

where the transformed diffusion velocities $\boldsymbol {V}_a$ are defined according to (2.55). The interdiffusion energy flux contributes the enthalpy transport due to diffusion and, hence, it vanishes in the single-component case. The effect of the interdiffusion energy flux is typically significant at the initial stages of the diffusion process and cannot be neglected.

Finally, the correction flux $\boldsymbol {q}^{corr}$ reads as

(3.36)\begin{equation} \boldsymbol{q}^{corr}=\frac{1}{2}\left(\frac{\omega_1-2}{\omega_1-\omega}\right) {\delta t}P \sum_{a=1}^M H_{a}\boldsymbol{\nabla} Y_a, \end{equation}

and is explained by the following consideration. The thermal flux is the mixture average of the component thermal fluxes, $\boldsymbol {q}^{th}=\sum _{a=1}^{M}Y_a \boldsymbol {q}_a^{th}$, where $\boldsymbol {q}_a^{th}= -\tau PC_{a,p}\boldsymbol {\nabla } T$ is the Fourier law for the component, $\tau$ is a parameter of no importance to the current consideration. On the other hand, in the single-component LBM, the thermal flux is $\boldsymbol {q}^{th}_{sc}= - \tau P\boldsymbol {\nabla } H_{sc}$, and with the single-component enthalpy $H_{sc}$ it returns the Fourier law in this case. However, the extension of the single-component to the multicomponent case so far invokes only the replacement of the single-component enthalpy with the ‘lumped’ mixture enthalpy and without any correction one gets

(3.37)\begin{equation} \boldsymbol{q}^{lump}={-}\tau P\boldsymbol{\nabla}\left(\sum_{a=1}^{M}Y_aH_a\right)=\boldsymbol{q}^{th}-\tau P\sum_{a=1}^{M}H_a\boldsymbol{\nabla} Y_a. \end{equation}

Thus, apart from the mixture-averaged Fourier law $\boldsymbol {q}^{th}$, the thermal flux also contains a spurious term. The spurious term is eliminated by the correction flux $\boldsymbol {q}^{corr}$ (3.36), where the prefactor is chosen by considering the hydrodynamic limit; see appendix B. The correction flux vanishes if all components are thermodynamically indistinguishable, that is, if all species have the same specific heat. In many cases, the correction flux contributes negligibly, for example, for air at moderate temperatures where the standard-air assumptions for diatomic molecules holds to a good approximation.

3.3. Hydrodynamic limit of the two-population lattice Boltzmann model for mixtures

Constraints on the pertinent equilibrium and quasi-equilibrium moments (3.26)–(3.31), (3.35), (3.36) and (3.32) are sufficient to study the hydrodynamic limit of the two-population lattice Boltzmann system (3.21) and (3.22) without a complete specification of the equilibrium and the quasi-equilibrium populations. The analysis follows the route of expanding the propagation to second order in the time step $\delta t$ and evaluating the moments of the resulting expansion. Details of the derivation are included in appendix B, here we present the final result.

The continuity equation:

(3.38)\begin{equation} \partial_t \rho + \boldsymbol{\nabla} \cdot(\rho \boldsymbol{u})=0. \end{equation}

The momentum equation:

(3.39)\begin{equation} \partial_t (\rho\boldsymbol{u}) + \boldsymbol{\nabla} \cdot({\rho\boldsymbol{u}\otimes\boldsymbol{u} })+ \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\rm \pi}=0. \end{equation}

Here, the pressure tensor $\boldsymbol {{\rm \pi} }$ reads as

(3.40)\begin{equation} \boldsymbol{\rm \pi}=P\boldsymbol{I} -\mu \left( \boldsymbol{S} -\frac{2}{D} (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u})\boldsymbol{I} \right) -\varsigma (\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}) \boldsymbol{I}, \end{equation}

where $\boldsymbol {S}$ is the strain rate,

(3.41)\begin{equation} \boldsymbol{S}=\boldsymbol{\nabla}\boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^{{\dagger}}. \end{equation}

The dynamic viscosity $\mu$ and the bulk viscosity $\varsigma$ are related to the relaxation parameter $\omega$,

(3.42)\begin{gather} \mu = \left( \frac{1}{\omega} - \frac{1}{2}\right) P{\delta t}, \end{gather}

(3.43)\begin{gather}\varsigma =\left( \frac{1}{\omega}-\frac{1}{2}\right)\left( \frac{2}{D} - \frac{R}{C_v} \right) P{\delta t}. \end{gather}

The energy equation:

(3.44)\begin{equation} \partial_t (\rho E)+\boldsymbol{\nabla}\cdot(\rho E\boldsymbol{u})+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}+\boldsymbol{\nabla}\cdot(\boldsymbol{\rm \pi}\boldsymbol{u})=0. \end{equation}

Here, the heat flux $\boldsymbol {q}$ reads as

(3.45)\begin{equation} \boldsymbol{q}={-}\lambda\boldsymbol{\nabla} T+\rho\sum_{a=1}^{M}H_aY_a\boldsymbol{V}_a. \end{equation}

The first term is the Fourier law of thermal conduction, with thermal conductivity $\lambda$ related to the relaxation parameter $\omega _1$,

(3.46)\begin{equation} \lambda= \left(\frac{1}{\omega_1} - \frac{1}{2}\right) P C_p{\delta t}. \end{equation}

The second term in (3.45) is the interdiffusion energy flux. Some comments are in order.

1. (i) The continuity, the momentum and the energy equations are the standard equations for multicomponent compressible mixtures (Williams Reference Williams1985; Bird et al. Reference Bird, Stewart and Lightfoot2006).

2. (ii) The bulk viscosity vanishes if all components are monatomic, $\bar {C}_{a,v}=DR_U/2$.

3. (iii) Introducing the thermal diffusivity $\alpha =\lambda /\rho C_p$ and the kinematic viscosity $\nu =\mu /\rho$, the Prandtl number becomes

   (3.47)\begin{equation} Pr = \frac{\nu}{\alpha} = \frac{\omega_1(2-\omega)}{\omega(2-\omega_1)}.\end{equation}

4. (iv) Using the equation of state of the mixture (2.4) in (3.42) and (3.46), relaxation parameters $\omega$ and $\omega _1$ are expressed in terms of dynamic viscosity and thermal conductivity,

   (3.48)\begin{gather} \omega = \frac{2 P\delta t}{P\delta t+2 \mu}, \end{gather}
   (3.49)\begin{gather}\omega_1 = \frac{2 PC_p\delta t}{PC_p\delta t+2 \lambda}. \end{gather}

Finally, the dynamic viscosity $\mu$ and the thermal conductivity $\lambda$ of the mixture at any point is evaluated as a function of the local composition by using the methods described in Wilke (Reference Wilke1950) and Mathur, Tondon & Saxena (Reference Mathur, Tondon and Saxena1967), respectively:

(3.50)\begin{equation} \mu = \sum_{a=1}^{M} \frac{\mu_a X_a}{\displaystyle\sum_{b=1}^{M} \phi_{ab} X_b}. \end{equation}

Here $\mu _a$ are the dynamic viscosity of the components while the dimensionless factor $\phi _{ab}$ is given by

(3.51)\begin{equation} \phi_{ab} = \frac{\displaystyle\left[1+\sqrt{\frac{\mu_a}{\mu_b} \sqrt{\frac{m_b}{m_a}}} \right]^2}{\displaystyle\sqrt{8} \sqrt{1+\frac{m_a}{m_b}}}. \end{equation}

The thermal conductivity of the mixture $\lambda$ is calculated from the thermal conductivity of the components $\lambda _a$,

(3.52)\begin{equation} \lambda = \frac{1}{2}\left(\sum_{a=1}^{M} X_a \lambda_a + \frac{1}{\displaystyle\sum_{a=1}^{M} \frac{X_a}{\lambda_a}} \right). \end{equation}

3.4. Realization on the standard lattice

3.4.1. Equilibrium and quasi-equilibrium

In order to finalize the construction of the lattice Boltzmann equations for the mixture, we need to specify the choice of the momentum and the energy lattices, and to provide the corresponding equilibrium and quasi-equilibrium populations. To that end, the single-component lattice Boltzmann models satisfying the moment constraints of § 3.2 are known in the literature. These employ higher-order lattices with a relatively large number of discrete velocities such as $D2Q49$ ($Q=49$ in two dimensions, Frapolli et al. Reference Frapolli, Chikatamarla and Karlin2015) or $D3Q39$ ($Q=39$ in three dimensions, Frapolli, Chikatamarla & Karlin Reference Frapolli, Chikatamarla and Karlin2020).

In this paper we develop the standard $D3Q27$ lattice realization as in the above case of the species LBM of § 2.4. We thus consider a two-dimensional compressible single-component lattice Boltzmann model by Saadat et al. (Reference Saadat, Bösch and Karlin2019) on the standard $D2Q9$ velocity set. Since the $D2Q9$ and $D3Q27$ belong to the same family of product lattices, cf. Karlin & Asinari (Reference Karlin and Asinari2010), it is natural to consider compressible LBM of Saadat et al. (Reference Saadat, Bösch and Karlin2019) for our purpose. In the following, we extend the model of Saadat et al. (Reference Saadat, Bösch and Karlin2019) to the three-dimensional $D3Q27$ discrete velocities set.

For the evaluation of the equilibrium $g_i^{eq}$ and of the quasi-equilibrium $g_i^*$ of the energy lattice, we proceed with the following ansatz $\mathcal {G}_i$, parameterized with a scalar $\theta \in [0,1]$, a scalar $M_0$, a vector $\boldsymbol {m}$ and a second-order tensor $\boldsymbol {M}$,

(3.53)\begin{gather} \mathcal{G}_i(\theta, M_0,\boldsymbol{m},\boldsymbol{M})= h_i(\theta, M_0,\boldsymbol{m},\boldsymbol{M}) + \boldsymbol{B}_{i}\boldsymbol{\cdot} \boldsymbol{Z}(\theta,M_0,\boldsymbol{M}), \end{gather}

(3.54)\begin{gather}h_i(\theta,M_0,\boldsymbol{m},\boldsymbol{M}) = w_i(\theta) \left(M_0 + \frac{ \boldsymbol{m} \boldsymbol{\cdot}\boldsymbol{c}_{i}}{\theta} + \frac{(\boldsymbol{M}-M_0 \theta \boldsymbol{I}): (\boldsymbol{c}_{i}\otimes\boldsymbol{c}_{i} - \theta \boldsymbol{I})}{2 \theta^2}\right). \end{gather}

Here, the weights $w_i$ are calculated in the product form as

(3.55)\begin{equation} w_i = w_{c_{ix}} w_{c_{iy}} w_{c_{iz}},\end{equation}

based on the fundamental triplet,

(3.56a–c)\begin{equation} w_{0} = 1 - \theta, \quad w_{1} = \frac{\theta}{2},\quad w_{{-}1} = \frac{\theta}{2}.\end{equation}

Furthermore, in (3.53), $\boldsymbol {Z}$ is a vector with the components

(3.57)\begin{equation} Z_{\alpha} = \frac{(1-3 \theta)}{2 \theta} (M_{\alpha \alpha} - \theta M_0). \end{equation}

Here, $M_{\alpha \alpha }$ is the diagonal component of the second-order tensor $\boldsymbol {M}$, while the components of vectors $\boldsymbol {B}_i$ are defined as follows:

(3.58)\begin{equation} \left.\begin{aligned} B_{i \alpha} &= 1 & \text{for } c_i^2=0, \\ B_{i \alpha} &={-}\tfrac{1}{2} \lvert c_{i \alpha} \rvert & \text{for } c_i^2=1, \\ B_{i \alpha} &= 0 & \text{otherwise}. \end{aligned}\right\} \end{equation}

Note that, when the parameter $\theta$ is set to the lattice reference temperature $\theta =1/3$, the term in (3.57) vanishes and the remaining term (3.54) becomes the familiar second-order Grad's approximation (Grad Reference Grad1949). By construction, the form (3.53) satisfies the moment relations, for any $\theta$:

(3.59)\begin{equation} \sum_{i=0}^{Q-1}\{1, \boldsymbol{c}_i, \boldsymbol{c}_i\otimes\boldsymbol{c}_i\}\mathcal{G}_i(\theta, M_0,\boldsymbol{m},\boldsymbol{M})=\{M_0,\boldsymbol{m},\boldsymbol{M}\}. \end{equation}

The equilibrium and the quasi-equilibrium populations $g_i^{eq}$ and $g_i^{*}$ are defined with the help of the form (3.53) by specifying the parameters $\theta =RT$, $M_0=\rho E$ and $\boldsymbol {M}=\boldsymbol {R}^{eq}$ (3.29) in both cases, and $\boldsymbol {m}=\boldsymbol {q}^{eq}$ (3.28) or $\boldsymbol {m}=\boldsymbol {q}^*$ (3.31) for the equilibrium or the quasi-equilibrium, respectively, following Saadat et al. (Reference Saadat, Bösch and Karlin2019):

(3.60)\begin{gather} g_i^{eq} = \mathcal{G}_i( RT,\rho E,\boldsymbol{q}^{eq},\boldsymbol{R}^{eq}), \end{gather}

(3.61)\begin{gather}g_i^{*} = \mathcal{G}_i( RT,\rho E,\boldsymbol{q}^{*} , \boldsymbol{R}^{eq}). \end{gather}

We shall now proceed with identifying the equilibrium of the momentum lattice and the modification of the lattice Boltzmann equation necessary for the $D3Q27$ model.

3.4.2. Augmented lattice Boltzmann equation for the momentum lattice

Equilibrium populations of the momentum lattice $f_i^{eq}$ are evaluated in the conventional way with the help of the product form (2.34) and using $\boldsymbol {\xi }=\boldsymbol {u}$ and $\zeta =RT$,

(3.62)\begin{equation} f_i^{eq} = \rho \varPsi_i(\boldsymbol{u},RT).\end{equation}

It is well known that the diagonal element of the equilibrium third-order moment of the momentum lattice $Q_{\alpha \alpha \alpha }^{eq}$ cannot satisfy the required moment relation (3.27). This happens due to the lattice constraint, $c_{i\alpha }^3=c_{i\alpha }$, which makes the diagonal third-order moments linearly dependent on the momentum, cf., e.g., Karlin & Asinari (Reference Karlin and Asinari2010). Following Saadat et al. (Reference Saadat, Bösch and Karlin2019), we consider the augmented lattice Boltzmann equation on the momentum lattice as

(3.63)\begin{equation} f_i(\boldsymbol{x}+\boldsymbol{c}_i,t+1) = f_i(\boldsymbol{x},t) + \omega (\,f_i^{eq} -f_i) + \boldsymbol{A}_{i}\boldsymbol{\cdot} \boldsymbol{X}, \end{equation}

where $\boldsymbol {X}$ is the vector with the components $\alpha =x,y,z$,

(3.64)\begin{equation} X_{\alpha} ={-}\partial_\alpha \left[ \left( \frac{1}{\omega}-\frac{1}{2} \right)\delta t \partial_\alpha (\rho u_\alpha (1 - 3 R T) - \rho u_\alpha^3) \right],\end{equation}

and where the components of vectors $\boldsymbol {A}_i$ are defined as

(3.65)\begin{equation} \left.\begin{aligned} A_{i \alpha} &= \tfrac{1}{2} c_{i \alpha} & \text{for } c_i^2=1, \\ A_{i \alpha} &= 0 & \text{otherwise}. \end{aligned}\right\} \end{equation}

This completes the realization of the mixture momentum and energy lattice Boltzmann equations on the standard $D3Q27$ lattice. In the next section we shall specify the coupling between the species and the mixture lattice Boltzmann subsystems.

3.5. Weak and strong coupling

3.5.2. Strong coupling: Matching of mixture density and momentum

With the two subsystems, the species and the mixture, first constructed independently from each other and after that being coupled weakly in the way described previously, we are left with two independent definitions of the mixture density and the mixture momentum. On the one hand, the mixture density $\rho$ (3.16) and the mixture momentum $\rho \boldsymbol {u}$ (3.17) are defined as the moments of the populations $f_i$. On the other hand, the same quantities are defined with the species populations as the sum of partial densities and partial momenta. The number of the conservation laws for the species subsystem is $M+D$, while for the mixture subsystem it is $D+2$. The total number of the conservation laws in the weakly coupled combined system is $M+2D+2$. Thus, the weakly coupled system is in excess of $D+1$ conservation laws as compared with the $M+D+1$ conservation laws of the mixture. This redundancy can be eliminated by removing one set of species populations (here, the $M$th) and writing

(3.66)\begin{equation} f_{Mi}=f_{i}-\sum_{a=1}^{M-1}f_{ai}. \end{equation}

As a consequence, the $M$th component is not an independent field anymore but is slaved to the remaining species and mixture populations.

While this is the method of choice to avoid over-determination of the system, various other, weaker coupling strategies exist. For example, considering the mixture density and mixture momentum defined by the momentum lattice as primary, the excess in conservation laws can be resolved by replacing the density of a selected component $\rho _M$ by the deficit of density once the $M-1$ other components are taken into account. Similarly, the momentum of the $M$th component $\rho _M\boldsymbol {u}_M$ accounts for the deficit of the mixture momentum once the momenta of the other species are counted,

(3.67)\begin{gather} \rho_M=\rho-\sum_{a=1}^{M-1}\rho_{a}=\sum_{i=1}^{Q-1}f_i- \sum_{a=1}^{M-1}\sum_{i=1}^{Q-1}f_{ai}, \end{gather}

(3.68)\begin{gather}\rho_M\boldsymbol{u}_M=\rho\boldsymbol{u}-\sum_{a=1}^{M-1}\rho_{a}\boldsymbol{u}_a= \sum_{i=1}^{Q-1}f_i\boldsymbol{c}_i-\sum_{a=1}^{M-1}\sum_{i=1}^{Q-1}f_{ai}\boldsymbol{c}_i. \end{gather}

In other words, only the density and momentum of the $M$th component are slaved by the corresponding mixture quantities and the rest of the mixture composition. Hence, the lattice Boltzmann equation for the $M$th component becomes purely relaxational, stripped of its conservation law. At the same time, the total momentum conservation is slaved by the momentum conservation of the momentum lattice. Since the relation (3.66) also implies (3.67) and (3.68), the number of independent conservation laws in both versions of the strongly coupled system is thus

(3.69)\begin{equation} (D+2)+ [(M+D)-1-D]=M+D+1, \end{equation}

and corresponds to the locally conserved fields, $\rho _1,\ldots ,\rho _{M-1}$, $\rho$, $\rho \boldsymbol {u}$ and $\rho E$. While the assignment of the slaved component $M$ is not unique, it is advisable to select the component which carries the majority of mass in the mixture.

In practice all couplings, the weak or the strong versions, yield identical results but the strongest coupling as in (3.66) is recommended as it reduces the number of lattices from $M+2$ in other coupling strategies to $M+1$.