1. Introduction
Reactive flows are a complex physicochemical phenomenon where different chemical species collide randomly and react violently, various interfacial and/or mechanical structures coexist, the chemical, hydrodynamic and thermodynamic nonequilibrium effects play significant roles [1–4]. Due to its practical importance in both nature and society, reacting flows have been widely studied in human history. With the rapid development of computer hardware and computational science in recent decades, numerical simulation has become indispensable for academic research. Roughly speaking, there are three levels of physical description of reactive fluids, i.e., the macroscopic, mesoscopic and microscopic models. The most commonly used method is the macroscopic description based upon the continuous models, such as the reactive Euler or Navier–Stokes (NS) equations, where hydrodynamic quantities such as the density, velocity, temperature, and pressure are utilized to characterize the reactive system [1–3]. Despite their great success, the traditional hydrodynamic governing equations ignore detailed thermodynamic nonequilibrium effects that often play essential roles, especially in a microscopic system or a local structure with sharp physical gradients [3, 4]. By contrast, the microscopic description is generally based on molecular dynamics (MD) where the interaction potential between molecules is considered [5]. Although the exact position and velocity of each molecule can be obtained dynamically, the MD is not capable of mimicking a relatively large system due to its excessive computational cost.
To solve the aforementioned issues, one way is to resort to the mesoscopic description that bridges the microscopic molecular and macroscopic continuous models. As a kinetic methodology, the discrete Boltzmann method (DBM) is regarded as a promising tool for mimicking multi-physics flow phenomena [6–18]. The DBM is a variant version of the standard lattice Boltzmann method (LBM) that has achieved great success in simulating complex reactive or nonreactive flows [19–37]. Early in 1997, Succi et al proposed the pioneering simple extension of the LBM for numerical combustion where reactive flow dynamics is handled in the limit of fast chemistry [21] . In 2010, Chiavazzo et al presented the modelling of reactive flows based upon a coupling between accurate reduced reaction mechanism and the LBM simulation of flow phenomena [22]. In 2014, Kang et al proposed a thermal multicomponent LBM for catalytic reactive flows where the Bhatnagar-Gross-Krook (BGK) relaxation process is split into two parts: the first part is characterized by the relaxation toward an auxiliary state and the second part describes the relaxation toward the thermodynamic equilibrium [23]. In 2020, Dugast et al proposed a topology optimization algorithm based on a multi-relaxation-time (MRT) LBM coupled with a level-set method for reactive fluid flows, allowing higher Reynolds numbers flow simulations compared to the ordinary single-relaxation-time model [28]. In 2021, Lei and Luo developed a sophisticated LBM for reactive flows in porous media, where separate equations describe the evolution of multicomponent flows and chemical species [33]. In 2022, Jiang et al proposed an immersed boundary LBM for particle combustion with varying thermodynamic and transport properties, and conducted the hydrodynamics-resolved simulation of a char particle combustion [34].
Although a series of LBM studies have been performed for reactive flows over the past two decades, they have not been developed to simulate complex compressible reacting flows where significant hydrodynamic, thermodynamic and chemical nonequilibrium effects coexist and interact with each other. In fact, both the DBM and basic LBM are intended as solvers for a truncated version of the Boltzmann equation, and they have the following same advantages: (i) Simple scheme. It is easy to code the algorithm for the discrete Boltzmann equation with a simplified collision operator. (ii) Parallel computation. It is convenient for massively parallel computing as all the information transfer is local in time and space [38]. (iii) Boundary condition. It is straightforward to deal with complex geometry by changing the discrete distribution functions [26, 39]. Moreover, the differences between the DBM and basic LBM are as follows: (i) Phase-space discretization. Classical LBM solvers rely on a specific stencil to discretize the phase space and a Lagrangian approach is used for the time-space discretization (integration along characteristic lines) resulting in the streaming-collision algorithm; The DBM is a classical Eulerian solver for the discrete Boltzmann equation where the discretization is performed using a moment-matching approach. (ii) Physical model. Traditional LBMs mainly serve as the solver of NS equations or other partial differential equations, while the DBM is fairly equivalent to a modified hydrodynamic model plus a coarse-grained model of thermodynamic nonequilibrium behaviors. Namely, the DBM is strictly inherited by the Boltzmann equation and can describe nonequilibrium system beyond macroscopic governing equations, while the standard LBM is a numerical scheme characterized by the streaming-collision steps and can be used to model a large family of advection equations beyond the topic of the Boltzmann equation.
It should be pointed out that the last aforementioned difference is the key reason why the DBM has been developed. The DBM is derived from the nonequilibrium statistical physics and has been successfully applied to investigate compressible reacting flows [7–18], multiphase [40–43], and fluid instabilities [44–52], etc. In 2013, the pioneering DBM for combustion was developed by using a hybrid scheme, and the hydrodynamic and thermodynamic nonequilibrium effects were investigated around the detonation wave [7]. In 2016, the DBM was employed to probe detonation with negative temperature coefficient from three aspects: hydrodynamic quantities, nonequilibrium quantities and entropy productions [10]. In the same year, a DBM was constructed where one (another) set of distribution function describes chemical reactant (product) [9]. In 2017, a BGK DBM for multicomponent reactive flows was presented, where the relaxation time has the same value for one species and the Prandtl number is fixed to Pr = 1 [11]. In 2019, an efficient MRT DBM was developed to tackle steady or unsteady supersonic reactive flows [13]. Since 2021, the DBM has been extended to three-dimensional steady and unsteady detonation [15, 16]. In 2022, the DBM was used to study the quantitative discrepancy between equilibrium and nonequilibrium distribution functions around the detonation wave [17]. In 2023, via the DBM and the fast Fourier transform, the deviations of the velocity distribution function from the equilibrium state have been investigated and the kinetic moments of reaction terms have been discussed in the evolution of unsteady detonation [18].
Base on previous works [11, 13, 14, 50], we develop an MRT DBM with a split collision term for reacting flows with both hydrodynamic and thermodynamic nonequilibrium effects. Compared with the BGK model [11], the MRT DBM has various relaxation times for different nonequilibrium processes and a flexible Pr. In contrast to the single-distribution-function DBM [13, 14], the current DBM takes account of collisions and reactions among different chemical species. Different from the MRT DBM for multicomponent mixtures [50], this model involves the effects of chemical reaction and external force, and a splitting technique is applied to the collision term. In historical context, the splitting technique originated within the realm of plasma physics for electron-ion systems, where the distinct impacts of self-collision and cross-collision on each species were recognized [53, 54]. Physically, complex systems often undergo multiple evolution stages characterized by varying time scales, and their components, such as electrons or ions, may exhibit divergent temperatures [53, 54]. Consequently, employing an appropriate methodology becomes imperative to investigate such intricate physical systems, especially those involving reactive flows with numerous chemical species. This serves as the primary motivation for this study. Specifically, the MRT DBM with split collision is developted for a mixture system to delineate the effects of self-collision and cross-collision on individual species. In contrast, the DBM without split collision represents a simplified version of the current model, applicable under the condition where relaxation frequencies during self-collision match those in cross-collision scenarios.
The rest of the paper is organized as follows. First, we introduce details of the MRT DBM with the split collision term for compressible reacting flows. The model is subsequently validated by six benchmarks, i.e., the multicomponent diffusion, homogeneous mixture in the force field, Kelvin–Helmholtz (KH) instability, flame at constant pressure, opposing chemical reaction, and steady detonation. Finally, conclusions are drawn.
2. Discrete Boltzmann method
The coarse-grained physical modeling from the Boltzmann equation to the DBM mainly involves three key steps [45, 55]. (i) Simplification of the collision term: it is difficult to solve the original Boltzmann equation where the collision term is too complex in the integral form. Hence, it is necessary to simplify the collision term in order to utilize the Boltzmann equation. In this paper, we employ the widely used MRT collision model [8, 13, 47]. (ii) Discretization of the particle velocity: for the purpose of physical accuracy and numerical efficiency, the particle velocity space is discretized with the matrix inversion method [6, 13]. (iii) Description of nonequilibrium effects: the main purpose of the DBM is to probe and extract essential nonequilibrium information beyond traditional hydrodynamic models [7, 9, 44]. Note that the first two steps are for coarse-grained modeling, and the third one is the core of DBM. The last step is not only an extension to the first two, but also puts forward stricter physical requirements for the construction. To be specific, all kinetic moment relations in the DBM should be consistent with those in nonequilibrium statistical physics, and multi-physics fields (including density, temperature, and velocity) should be coupled naturally as various physical quantities are obtained from the kinetic moments of the same distribution function.
In this work, the MRT discrete Boltzmann equations, which describe the spatio-temporal evolution of reacting flows, take the form,
$\begin{eqnarray}\displaystyle \frac{\partial {{\boldsymbol{f}}}^{\sigma }}{\partial t}+{{\boldsymbol{v}}}^{\sigma }\cdot {\rm{\nabla }}{{\boldsymbol{f}}}^{\sigma }={{\boldsymbol{\Omega }}}^{\sigma }+{{\boldsymbol{F}}}^{\sigma }+{{\boldsymbol{R}}}^{\sigma }.\end{eqnarray}$
Here, t represents the time. The superscript σ = 1, 2, ..., Ns indicates chemical species with the number Ns in total. The column matrices ${{\boldsymbol{f}}}^{\sigma }={\left(\begin{array}{l}{f}_{1}^{\sigma }\ {f}_{2}^{\sigma }\ \cdots \ {f}_{N}^{\sigma }\end{array}\right)}^{{\rm{T}}}$, ${{\boldsymbol{\Omega }}}^{\sigma }={\left(\begin{array}{l}{{\rm{\Omega }}}_{1}^{\sigma }\ {{\rm{\Omega }}}_{2}^{\sigma }\ \cdots \ {{\rm{\Omega }}}_{N}^{\sigma }\end{array}\right)}^{{\rm{T}}}$, ${{\boldsymbol{F}}}^{\sigma }={\left(\begin{array}{l}{F}_{1}^{\sigma }\ {F}_{2}^{\sigma }\ \cdots \ {F}_{N}^{\sigma }\end{array}\right)}^{{\rm{T}}}$, and ${{\boldsymbol{R}}}^{\sigma }={\left(\begin{array}{l}{R}_{1}^{\sigma }\ {R}_{2}^{\sigma }\ \cdots \ {R}_{N}^{\sigma }\end{array}\right)}^{{\rm{T}}}$ denote the discrete distribution functions, collision terms, force terms, and reaction terms, respectively. The diagonal matrix ${{\boldsymbol{v}}}^{\sigma }=\mathrm{diag}\left(\begin{array}{l}{{\boldsymbol{v}}}_{1}^{\sigma }\ {{\boldsymbol{v}}}_{2}^{\sigma }\ \cdots \ {{\boldsymbol{v}}}_{N}^{\sigma }\end{array}\right)$ stands for the discrete velocities, and the subscript i = 1, 2, ..., N is the index of the discrete velocities, with N = 16 in this paper. As shown in figure 1, a discrete velocity set reads $\begin{eqnarray}{{\boldsymbol{v}}}_{i}^{\sigma }=\left\{\begin{array}{ll}\mathrm{cyc}:{v}_{a}^{\sigma }\left(\pm 1,0\right), & 1\leqslant i\leqslant 4,\\ {v}_{b}^{\sigma }\left(\pm 1,\pm 1\right), & 5\leqslant i\leqslant 8,\\ \mathrm{cyc}:{v}_{c}^{\sigma }\left(\pm 1,0\right), & 9\leqslant i\leqslant 12,\\ {v}_{d}^{\sigma }\left(\pm 1,\pm 1\right), & 13\leqslant i\leqslant 16,\end{array}\right.\end{eqnarray}$
where cyc indicates the cyclic permutation and (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$) are the flexible parameters. Besides, to take account of the part of internal energies due to molecular rotation and/or vibration, we introduce the symbol ${\eta }_{i}^{\sigma }={\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, and ${\eta }_{d}^{\sigma }$ for 1 ≤ i ≤ 4, 5 ≤ i ≤ 8, 9 ≤ i ≤ 12, and 13 ≤ i ≤ 16, respectively, where the parameters (${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$) are tunable as well.
Figure 1. Sketch of the discrete velocities. |
It is worth emphasizing that the values of ${{\boldsymbol{v}}}_{i}^{\sigma }$ and ${\eta }_{i}^{\sigma }$ can be adjusted to optimize the DBM robustness and accuracy. On the one hand, the parameters (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$) can be chosen around the values of flow velocity uσ and sound speed ${v}_{s}^{\sigma }=\sqrt{{\gamma }^{\sigma }{T}^{\sigma }/{m}^{\sigma }}$, where Tσ stands for the temperature, mσ the molar mass, γσ = (D + Iσ + 2)/(D + Iσ) the specific heat ratio, D = 2 the spatial dimension in this paper, and Iσ extra degrees of freedom due to molecular rotation and/or vibration. On the other hand, among (${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$), one (another) parameter should be less (greater) than ${\bar{\eta }}^{\sigma }=\sqrt{{I}^{\sigma }{T}^{\sigma }/{m}^{\sigma }}$, because the extra internal energy is about $\displaystyle \frac{1}{2}{m}^{\sigma }{\bar{\eta }}^{2}=\displaystyle \frac{1}{2}{I}^{\sigma }{T}^{\sigma }$ according to the equipartition of energy theorem.
2.1. Macroscopic quantities
The individual molar concentration nσ, mass density ρσ, velocity uσ, energy Eσ, and temperature Tσ are given by,6 ) and E in equation (11 ) are not conserved during chemical reaction and may be called the ‘sensible' or ‘total nonchemical' energies.
$\begin{eqnarray}{n}^{\sigma }=\displaystyle \sum _{i}{f}_{i}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{\rho }^{\sigma }={m}^{\sigma }{n}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{n}^{\sigma }{{\boldsymbol{u}}}^{\sigma }=\displaystyle \sum _{i}{f}_{i}^{\sigma }{{\boldsymbol{v}}}_{i}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{E}^{\sigma }=\displaystyle \frac{{m}^{\sigma }}{2}\sum _{i}{f}_{i}^{\sigma }\left(\left|\displaystyle {{\boldsymbol{v}}}_{i}^{\sigma }{|}^{2}+{{\eta }_{i}^{\sigma }}^{2}\right.\right),\end{eqnarray}$
$\begin{eqnarray}{T}^{\sigma }=\displaystyle \frac{2{E}^{\sigma }-{\rho }^{\sigma }{{u}^{\sigma }}^{2}}{\left(D+{I}^{\sigma }\right){n}^{\sigma }},\end{eqnarray}$
respectively. The mixing number density n, mass density ρ, velocity u, energy E, internal energy Eint, and temperature T are obtained from $\begin{eqnarray}n=\displaystyle \sum _{\sigma }{n}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}\rho =\displaystyle \sum _{\sigma }{\rho }^{\sigma },\end{eqnarray}$
$\begin{eqnarray}\rho {\boldsymbol{u}}=\displaystyle \sum _{\sigma }{\rho }^{\sigma }{{\boldsymbol{u}}}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}E=\displaystyle \sum _{\sigma }{E}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{E}_{\mathrm{int}}=E-\displaystyle \frac{1}{2}\rho {\left|{\boldsymbol{u}}\right|}^{2},\end{eqnarray}$
$\begin{eqnarray}T=\displaystyle \frac{2{E}_{\mathrm{int}}}{{\sum }_{\sigma }\left(D+{I}^{\sigma }\right){n}^{\sigma }},\end{eqnarray}$
respectively. Actually, energies Eσ in equation (It should be mentioned that, with the substitution of the equilibrium discrete distribution functions fσseq for the discrete distribution functions fσ, the formulas (3 )–(6 ) still hold. In fact, the aforementioned physical quantities in equations (3 )–(13 ) are macroscopic parameters that are statistical results of particles with random motion, and can be conveniently measured by traditional numerical or experimental methods.
2.2. Split collision term
In fact, the discrete Boltzmann equation (1 ) is a simplified form of the original Boltzmann equation, and the collision term is a reduced expression of its original nonlinear integral term. To be specific, the collision term is composed of three parts, i.e.,A . The square matrix ${{\boldsymbol{M}}}^{\sigma }=\left({M}_{{il}}^{\sigma }\right)$ and its inverse ${({{\boldsymbol{M}}}^{\sigma })}^{-1}=\left({\left({M}_{{il}}^{\sigma }\right)}^{-1}\right)$, both of which have N × N elements, act as the link between the velocity and moment spaces. To be specific,18 ) is equivalent to the following relationship,21 ) is a necessary condition to recover the NS equations in the hydrodynamic limit.
$\begin{eqnarray}{{\boldsymbol{\Omega }}}^{\sigma }={{\boldsymbol{\Omega }}}^{1\sigma }+{{\boldsymbol{\Omega }}}^{2\sigma }+{{\boldsymbol{\Omega }}}^{3\sigma },\end{eqnarray}$
in terms of $\begin{eqnarray}{{\boldsymbol{\Omega }}}^{1\sigma }=-{\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}{{\boldsymbol{S}}}^{1\sigma }\left({\hat{{\boldsymbol{f}}}}^{\sigma }-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}\right),\end{eqnarray}$
$\begin{eqnarray}{{\boldsymbol{\Omega }}}^{2\sigma }=-{\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}{{\boldsymbol{S}}}^{2\sigma }\left({\hat{{\boldsymbol{f}}}}^{\sigma {\rm{seq}}}-{\hat{{\boldsymbol{f}}}}^{\sigma {\rm{eq}}}\right),\end{eqnarray}$
and $\begin{eqnarray}{{\boldsymbol{\Omega }}}^{3\sigma }={\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}{\hat{{\boldsymbol{A}}}}^{\sigma },\end{eqnarray}$
where the diagonal matrix ${{\boldsymbol{S}}}^{1\sigma }=\mathrm{diag}\left(\begin{array}{l}{S}_{1}^{1\sigma }\ {S}_{2}^{1\sigma }\ \cdots \ {S}_{N}^{1\sigma }\end{array}\right)$ indicates the relaxation frequencies that control the relaxation speed of kinetic moments ${\hat{{\boldsymbol{f}}}}^{\sigma }={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma }\ {\hat{f}}_{2}^{\sigma }\ \cdots \ {\hat{f}}_{N}^{\sigma }\end{array}\right)}^{{\rm{T}}}$ approaching their individual intermediate equilibrium counterparts ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma \mathrm{seq}}\ {\hat{f}}_{2}^{\sigma \mathrm{seq}}\ \cdots \ {\hat{f}}_{N}^{\sigma \mathrm{seq}}\end{array}\right)}^{{\rm{T}}}$, and ${{\boldsymbol{S}}}^{2\sigma }\,=\mathrm{diag}\left(\begin{array}{l}{S}_{1}^{2\sigma }\ {S}_{2}^{2\sigma }\ \cdots \ {S}_{N}^{2\sigma }\end{array}\right)$ denotes the relaxation frequencies which govern the relaxation speed of ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}\,={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma \mathrm{seq}}\ {\hat{f}}_{2}^{\sigma \mathrm{seq}}\ \cdots \ {\hat{f}}_{N}^{\sigma \mathrm{seq}}\end{array}\right)}^{{\rm{T}}}$ approaching the ultimate equilibrium counterparts ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma \mathrm{eq}}\ {\hat{f}}_{2}^{\sigma \mathrm{eq}}\ \cdots \ {\hat{f}}_{N}^{\sigma \mathrm{eq}}\end{array}\right)}^{{\rm{T}}}$. Here ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}$ is the function of (nσ, uα, T), and ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}$ is expressed by substituting (${u}_{\alpha }^{\sigma }$, Tσ) for (uα, T) in the formula of ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}$, see appendix $\begin{eqnarray}{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}={{\boldsymbol{M}}}^{\sigma }{{\boldsymbol{f}}}^{\sigma \mathrm{eq}},\end{eqnarray}$
$\begin{eqnarray}{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}={{\boldsymbol{M}}}^{\sigma }{{\boldsymbol{f}}}^{\sigma \mathrm{seq}},\end{eqnarray}$
$\begin{eqnarray}{\hat{{\boldsymbol{f}}}}^{\sigma }={{\boldsymbol{M}}}^{\sigma }{{\boldsymbol{f}}}^{\sigma },\end{eqnarray}$
with ${{\boldsymbol{f}}}^{\sigma \mathrm{eq}}={\left(\begin{array}{l}{f}_{1}^{\sigma \mathrm{eq}}\ {f}_{2}^{\sigma \mathrm{eq}}\ \cdots \ {f}_{N}^{\sigma \mathrm{eq}}\end{array}\right)}^{{\rm{T}}}$ and ${{\boldsymbol{f}}}^{\sigma \mathrm{seq}}\,={\left(\begin{array}{l}{f}_{1}^{\sigma \mathrm{seq}}\ {f}_{2}^{\sigma \mathrm{seq}}\ \cdots \ {f}_{N}^{\sigma \mathrm{seq}}\end{array}\right)}^{{\rm{T}}}$. Here (${\hat{f}}^{\sigma \mathrm{eq}}$, ${\hat{f}}^{\sigma \mathrm{seq}}$, ${\hat{f}}^{\sigma }$) and (fσeq, fσseq, fσ) correspond to the moment and velocity spaces, respectively. In fact, equation ( $\begin{eqnarray}\iint {f}^{\sigma \mathrm{eq}}{\rm{\Psi }}{\rm{d}}{\boldsymbol{v}}{\rm{d}}\eta =\displaystyle \sum _{i}{f}_{i}^{\sigma \mathrm{eq}}{{\rm{\Psi }}}_{i},\end{eqnarray}$
in terms of $\Psi$ = 1, v, $\left(| {\boldsymbol{v}}{| }^{2}+{\eta }^{2}\right)$, vv, $\left(| {\boldsymbol{v}}{| }^{2}+{\eta }^{2}\right){\boldsymbol{v}}$, vvv, $\left(\left|{\boldsymbol{v}}{|}^{2}+{\eta }^{2}\right.\right){\boldsymbol{vv}}$, and their corresponding discrete counterparts $\Psi$i = 1, ${{\boldsymbol{v}}}_{i}^{\sigma }$, $\left(| {{\boldsymbol{v}}}_{i}^{\sigma }{| }^{2}+{{\eta }_{i}^{\sigma }}^{2}\right)$, ${{\boldsymbol{v}}}_{i}^{\sigma }{{\boldsymbol{v}}}_{i}^{\sigma }$, $\left(| {{\boldsymbol{v}}}_{i}^{\sigma }{| }^{2}+{{\eta }_{i}^{\sigma }}^{2}\right){{\boldsymbol{v}}}_{i}^{\sigma }$, ${{\boldsymbol{v}}}_{i}^{\sigma }{{\boldsymbol{v}}}_{i}^{\sigma }{{\boldsymbol{v}}}_{i}^{\sigma }$, $\left(| {{\boldsymbol{v}}}_{i}^{\sigma }{| }^{2}+{{\eta }_{i}^{\sigma }}^{2}\right){{\boldsymbol{v}}}_{i}^{\sigma }{{\boldsymbol{v}}}_{i}^{\sigma }$. In this work, the theoretical equilibrium distribution function is expressed by [11] $\begin{eqnarray}\begin{array}{l}{f}^{\sigma \mathrm{eq}}={n}^{\sigma }{\left(\displaystyle \frac{{m}^{\sigma }}{2\pi T}\right)}^{D/2}{\left(\displaystyle \frac{{m}^{\sigma }}{2\pi {I}^{\sigma }T}\right)}^{1/2}\\ \,\times \,\exp \left[-\displaystyle \frac{{m}^{\sigma }{\left|{\boldsymbol{v}}-{\boldsymbol{u}}\right|}^{2}}{2T}-\displaystyle \frac{{m}^{\sigma }{\eta }^{2}}{2{I}^{\sigma }T}\right].\end{array}\end{eqnarray}$
Actually, the formula (Note that in the simplification process from the original Boltzmann equation to the DBM, the physical quantities (such as density, momentum, energy, and lower-order kinetic moments) under consideration remain unchanged, while some other physical information (such as higher-order kinetic moments and the interactions between them) may be lost. The loss of relevant information constrains the applications of the physical model, and may lead to inaccuracy for particular real situations. To rectify this, Chapman–Enskog (CE) multiscale analysis can be employed to identify and amend the deficiencies in the physical model. For this purpose, an additional term ${\hat{{\boldsymbol{A}}}}^{\sigma }$ is incorporated into the collision term to make up for the missing relation between the physical quantities ${\hat{f}}_{5}^{\sigma }$, ${\hat{f}}_{6}^{\sigma }$, ${\hat{f}}_{7}^{\sigma }$, ${\hat{f}}_{8}^{\sigma }$, and ${\hat{f}}_{9}^{\sigma }$. Specifically, the term ${\hat{{\boldsymbol{A}}}}^{\sigma }={\left(\begin{array}{l}0\ \cdots \ 0\ {\hat{A}}_{8}^{\sigma }\ {\hat{A}}_{9}^{\sigma }\ 0\ \cdots \ 0\end{array}\right)}^{{\rm{T}}}$ depends upon
$\begin{eqnarray}{\hat{A}}_{8}^{\sigma }=2\left({S}_{8}^{1\sigma }-{S}_{5}^{1\sigma }\right){u}_{x}^{\sigma }{{\rm{\Delta }}}_{5}^{\sigma }+2\left({S}_{8}^{1\sigma }-{S}_{6}^{1\sigma }\right){u}_{y}^{\sigma }{{\rm{\Delta }}}_{6}^{\sigma },\end{eqnarray}$
$\begin{eqnarray}{\hat{A}}_{9}^{\sigma }=2\left({S}_{9}^{1\sigma }-{S}_{7}^{1\sigma }\right){u}_{y}^{\sigma }{{\rm{\Delta }}}_{7}^{\sigma }+2\left({S}_{9}^{1\sigma }-{S}_{6}^{1\sigma }\right){u}_{x}^{\sigma }{{\rm{\Delta }}}_{6}^{\sigma },\end{eqnarray}$
with $\begin{eqnarray}\begin{array}{l}{{\rm{\Delta }}}_{5}^{\sigma }=\displaystyle \frac{2{n}^{\sigma }{T}^{\sigma }}{{S}_{5}^{1\sigma }{m}^{\sigma }}\left(\displaystyle \frac{1-D-{I}^{\sigma }}{D+{I}^{\sigma }}\displaystyle \frac{\partial {u}_{x}^{\sigma }}{\partial x}+\displaystyle \frac{1}{D+{I}^{\sigma }}\displaystyle \frac{\partial {u}_{y}^{\sigma }}{\partial y}\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{{\rm{\Delta }}}_{6}^{\sigma }=-\displaystyle \frac{{n}^{\sigma }{T}^{\sigma }}{{S}_{6}^{1\sigma }{m}^{\sigma }}\left(\displaystyle \frac{\partial {u}_{x}^{\sigma }}{\partial y}+\displaystyle \frac{\partial {u}_{y}^{\sigma }}{\partial x}\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{{\rm{\Delta }}}_{7}^{\sigma }=\displaystyle \frac{2{n}^{\sigma }{T}^{\sigma }}{{S}_{7}^{1\sigma }{m}^{\sigma }}\left(\displaystyle \frac{1}{D+{I}^{\sigma }}\displaystyle \frac{\partial {u}_{x}^{\sigma }}{\partial x}+\displaystyle \frac{1-D-{I}^{\sigma }}{D+{I}^{\sigma }}\displaystyle \frac{\partial {u}_{y}^{\sigma }}{\partial y}\right).\end{array}\end{eqnarray}$
It should be further explained that the physical meaning of the first two parts in equation (14 ) is as follows: there are two split steps during the thermodynamic relaxation process. The distribution function ${\hat{{\boldsymbol{f}}}}^{\sigma }$ firstly approaches its temporary equilibrium state ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}$ under the control of relaxation frequency S1σ, then tends toward the local ultimate equilibrium state ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}$ with relaxation frequency S2σ. It is worth mentioning that, there are more flexible parameters in the two-step-relaxation collision term, which is suitable for a wider application range of physical systems. Through the CE expansion, the relations can be determined between the relaxation parameters and other physical quantities, such as the nonequilibrium quantities (48 )–(52 ), diffusivity (55 ), thermal conductivity (B10 ), and dynamic viscosity (B11 ). Consequently, compared with the one-step-relaxation MRT or BGK model, the two-step-relaxation collision term presents a more detailed relationship between the thermodynamic relaxation process and nonequilibrium effects.
Moreover, substituting equations (15 )–(17 ) into (14 ) leads to the following expression28 ) (called two-step-relaxation collision model) reduces to the popular MRT model (named a one-step-relaxation collision model)23 ) and (24 ). Namely, the additional term ${\hat{{\boldsymbol{A}}}}^{\sigma }$ disappears in the single-relaxation case. A widely used single-relaxation model is the BGK model. In fact, both MRT and BGK models are also suitable for a fluid system if there is only one component or the average mixing effect is under consideration. Besides, in order to achieve local momentum conservation for a two-component system, the relaxation times of the two components should be equal: ${\tau }^{\sigma }={\tau }^{\bar{\sigma }}=\tau $, which is a constraint of the single-relaxation BGK model [54].
$\begin{eqnarray}\begin{array}{l}{{\boldsymbol{\Omega }}}^{\sigma }=-{\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}\left[{{\boldsymbol{S}}}^{1\sigma }\left({\hat{{\boldsymbol{f}}}}^{\sigma }-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}\right)\right.\\ \,+\,\left.{{\boldsymbol{S}}}^{2\sigma }\left({\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}\right)-{\hat{{\boldsymbol{A}}}}^{\sigma }\right].\end{array}\end{eqnarray}$
Clearly, in the case of S1σ = S2σ = Sσ, the split collision model ( $\begin{eqnarray}{{\boldsymbol{\Omega }}}^{\sigma }=-{\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}\left[{{\boldsymbol{S}}}^{\sigma }\left({\hat{{\boldsymbol{f}}}}^{\sigma }-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}\right)-{\hat{{\boldsymbol{A}}}}^{\sigma }\right],\end{eqnarray}$
which further reduces to the single-relaxation model $\begin{eqnarray}{{\boldsymbol{\Omega }}}^{\sigma }=-\displaystyle \frac{1}{{\tau }^{\sigma }}\left({{\boldsymbol{f}}}^{\sigma }-{{\boldsymbol{f}}}^{\sigma \mathrm{eq}}\right),\end{eqnarray}$
if Sσ = I/τσ, and I denotes the unit tensor. Clearly, ${\hat{A}}_{8}^{\sigma }={\hat{A}}_{9}^{\sigma }=0$ when ${S}_{5}^{1\sigma }={S}_{6}^{1\sigma }={S}_{7}^{1\sigma }={S}_{8}^{1\sigma }={S}_{9}^{1\sigma }$, as seen in equations (In addition, the collision term in equation (14 ) can be rewritten as,
$\begin{eqnarray}{{\boldsymbol{\Omega }}}^{\sigma }={{\boldsymbol{\Omega }}}^{1\sigma * }+{{\boldsymbol{\Omega }}}^{2\sigma * }+{{\boldsymbol{\Omega }}}^{3\sigma },\end{eqnarray}$
with $\begin{eqnarray}{{\boldsymbol{\Omega }}}^{1\sigma * }=-{\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}{{\boldsymbol{S}}}^{1\sigma * }\left({\hat{{\boldsymbol{f}}}}^{\sigma }-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}\right),\end{eqnarray}$
$\begin{eqnarray}{{\boldsymbol{\Omega }}}^{2\sigma * }=-{\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}{{\boldsymbol{S}}}^{2\sigma * }\left({\hat{{\boldsymbol{f}}}}^{\sigma }-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}\right),\end{eqnarray}$
where S1σ* = S1σ − S2σ and S2σ* = S2σ. The terms Ω1σ* and Ω2σ* are related to the self-collision and cross-collision among various particles, respectively. In fact, both self-collision and cross-collision affect the evolution of the discrete distribution functions. In other words, the effects of self-collision and cross-collision on the evolution of physical systems are taken into consideration.2.3. Force term
Physically, the force term denotes the change rate of distribution function due to the external force. How to calculate the force term is a key to an accurate physical model. In this part, we introduce three ways to obtain the mathematical expression of the force term. The first two methods, which were actually proposed for two-component fluids in [45], are extended to multicomponent systems in this work. The last method that is named the matrix inversion method [6, 13] is developed for multicomponent systems for the first time as well.
Method I
Via the Taylor expansion, it can be found that the main part of the distribution function is the equilibrium distribution function in a system not too far from equilibrium [45]. Theoretically, fσ is close to fσseq(nσ, uσ, Tσ) rather than fσeq(nσ, u, T), especially in a non-premixed or partially premixed system. Hence the approximation fσ ≈ fσseq(nσ, uσ, Tσ) can be used as follows,35 ) is a conventional way to calculate the force terms [46, 47].
$\begin{eqnarray}{F}^{\sigma }=-{{\boldsymbol{a}}}^{\sigma }\cdot \displaystyle \frac{\partial {f}^{\sigma }}{\partial {\boldsymbol{v}}}\approx -{{\boldsymbol{a}}}^{\sigma }\cdot \displaystyle \frac{\partial {f}^{\sigma \mathrm{seq}}}{\partial {\boldsymbol{v}}}={{\boldsymbol{a}}}^{\sigma }\cdot \left({\boldsymbol{v}}-{{\boldsymbol{u}}}^{\sigma }\right)\displaystyle \frac{{m}^{\sigma }}{{T}^{\sigma }}{f}^{\sigma \mathrm{seq}},\end{eqnarray}$
where ${{\boldsymbol{a}}}^{\sigma }={a}_{\alpha }^{\sigma }{{\boldsymbol{e}}}_{\alpha }$ stands for the body acceleration of species σ, and eα the unit vector in the α direction. Then, the force terms are obtained in the discretization form directly $\begin{eqnarray}{F}_{i}^{\sigma }={{\boldsymbol{a}}}^{\sigma }\cdot \left({{\boldsymbol{v}}}_{i}^{\sigma }-{{\boldsymbol{u}}}^{\sigma }\right)\displaystyle \frac{{m}^{\sigma }}{{T}^{\sigma }}{f}_{i}^{\sigma \mathrm{seq}}.\end{eqnarray}$
In fact, equation (Method II
The force terms are used to incorporate forcing effects into the Boltzmann equation. According to its physical meaning, the force terms can be expressed by the change of discrete distribution functions $\delta {f}_{i}^{\sigma }$ due to the external force over a small time interval δt, i.e.,
$\begin{eqnarray}{F}_{i}^{\sigma }={\left.\displaystyle \frac{\partial {f}_{i}^{\sigma }}{\partial t}\right|}_{\mathrm{Force}}=\mathop{\mathrm{lim}}\limits_{\delta t\to 0}\,{\left.\displaystyle \frac{\delta {f}_{i}^{\sigma }}{\delta t}\right|}_{\mathrm{Force}}\approx {\left.\displaystyle \frac{{\rm{\Delta }}{f}_{i}^{\sigma \mathrm{seq}}}{{\rm{\Delta }}t}\right|}_{\mathrm{Force}},\end{eqnarray}$
where ${\rm{\Delta }}{f}_{i}^{\sigma \mathrm{seq}}$ represents the corresponding change of the equilibrium distribution function within a time step Δt and is a function of the concentration, velocity and temperature.In classical physics, the impulse (work) done by an external force changes the momentum (kinetic energy) of a system directly, while the mass, internal energy, and temperature remain constant in the force field. In other words, the force changes the velocity and energy of the fluid components, but does not have a direct influence on the density or temperature. Consequently, ${\rm{\Delta }}{f}_{i}^{\sigma \mathrm{seq}}$ in equation (36 ) can be written as
$\begin{eqnarray}{\rm{\Delta }}{f}_{i}^{\sigma \mathrm{seq}}={f}_{i}^{\sigma \mathrm{seq}\dagger }-{f}_{i}^{\sigma \mathrm{seq}},\end{eqnarray}$
where the equilibrium distribution functions change from ${f}_{i}^{\sigma \mathrm{seq}}={f}_{i}^{\sigma \mathrm{seq}}\left({n}^{\sigma },{{\boldsymbol{u}}}^{\sigma },{T}^{\sigma }\right)$ to ${f}_{i}^{\sigma \mathrm{seq}\dagger }={f}_{i}^{\sigma \mathrm{seq}}\left({n}^{\sigma },{{\boldsymbol{u}}}^{\sigma \dagger },{T}^{\sigma }\right)$, and the flow velocity changes from uσ to uσ† = uσ + aσΔt within a time step.Theoretically, because of the external force, the energy of component σ change from Eσ into Eσ† = Eσ + ρσuσ · aσΔt, then it can be derived from equation (7 ) that the temperature of component σ is38 ) that the expression in equation (37 ) is of the second order accuracy.
$\begin{eqnarray}{T}^{\sigma \dagger }={T}^{\sigma }-\displaystyle \frac{{m}^{\sigma }{\left|{{\boldsymbol{a}}}^{\sigma }\right|}^{2}}{D+{I}^{\sigma }}{\left({\rm{\Delta }}t\right)}^{2}\approx {T}^{\sigma }.\end{eqnarray}$
Clearly, Tσ† equals Tσ as Δt approaches zero. In other words, the temperature is not changed by the external force. It can be found from equation (Method III
Let us consider the following relation [13], 21 ), and Fσ and ${F}_{i}^{\sigma }$ denote the force terms in the continuous and discrete velocity spaces, respectively.
$\begin{eqnarray}\iint {F}^{\sigma }{\rm{\Psi }}{\rm{d}}{\boldsymbol{v}}{\rm{d}}\eta =\displaystyle \sum _{i}{F}_{i}^{\sigma }{{\rm{\Psi }}}_{i},\end{eqnarray}$
where $\Psi$ and $\Psi$i are the same as those in equation (In fact, the formula (39 ) is equivalent to the following matrix formA . From equation (40 ), the force terms can be expressed by
$\begin{eqnarray}{\hat{{\boldsymbol{F}}}}^{\sigma }={{\boldsymbol{M}}}^{\sigma }{{\boldsymbol{F}}}^{\sigma },\end{eqnarray}$
where the elements of ${\hat{{\boldsymbol{F}}}}^{\sigma }$ are given in appendix $\begin{eqnarray}{{\boldsymbol{F}}}^{\sigma }={\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}{\hat{{\boldsymbol{F}}}}^{\sigma }.\end{eqnarray}$
It is noteworthy that there are two similarities among the above three methods. (I) Based upon the approximation fσ ≈ fσseq(nσ, uσ, Tσ), the force terms are expressed with the discrete equilibrium distribution functions; (II) the relations satisfied by the force terms are sufficient to recover the NS equations, see appendix B . Besides, the differences among above three methods are as follows. (I) As for Method I, there are nine relations satisfied by the force terms, which are the necessary and sufficient conditions to recover the NS equations. In contrast, besides the nine relations, another seven relations are satisfied by the force terms in the last two methods as well. That is to say, there are sixteen relationships in Methods II and III, respectively. (II) As shown in equations (35 ) and (37 ), the force terms are a function of the equilibrium discrete distribution functions. Consequently, it is necessary to calculate the equilibrium discrete distribution functions in the program for Method I or II. In contrast, as shown in equation (41 ), the force terms are computed by using an inverse matrix in the last method, which has a higher computational efficiency. In other words, the matrix inversion method has the characteristic of high physical accuracy and computational efficiency, see table 2. Consequently, this methodology is utilized to calculate the reaction terms in the next subsection as well.
2.4. Reaction term
The reaction term, which represents the change rate of distribution function because of the chemical reaction, satisfies the following relationship [13], 21 ), and Rσ and ${R}_{i}^{\sigma }$ are the reaction terms in the continuous and discrete velocity spaces, respectively. The expression of Rσ reads [8],45 ), ωov stands for the chemical reaction rate, and Q stands for the chemical heat release of reactant per unit mole.
$\begin{eqnarray}\iint {R}^{\sigma }{\rm{\Psi }}{\rm{d}}{\boldsymbol{v}}{\rm{d}}\eta =\displaystyle \sum _{i}{R}_{i}^{\sigma }{{\rm{\Psi }}}_{i},\end{eqnarray}$
where $\Psi$ and $\Psi$i are the same as those in equation ( $\begin{eqnarray}\begin{array}{l}{R}^{\sigma }={f}^{\sigma \mathrm{eq}}\displaystyle \frac{{{n}^{\sigma }}^{{\prime} }}{{n}^{\sigma }}\\ \,+\,{f}^{\sigma \mathrm{eq}}\displaystyle \frac{-\left(1+D\right){I}^{\sigma }T+{m}^{\sigma }{I}^{\sigma }{\left|{\boldsymbol{v}}-{\boldsymbol{u}}\right|}^{2}+{m}^{\sigma }{\eta }^{2}}{2{I}^{\sigma }{T}^{2}}T^{\prime} ,\end{array}\end{eqnarray}$
where ${n}^{\sigma ^{\prime} }$ stands for the concentration variation rate of species σ due to the chemical reaction, $\begin{eqnarray}T^{\prime} =\displaystyle \frac{2\left[E^{\prime} \cdot {\sum }_{\sigma }{n}^{\sigma }\left(D+{I}^{\sigma }\right)-{E}_{\mathrm{int}}\cdot {\sum }_{\sigma }{{n}^{\sigma }}^{{\prime} }\left(D+{I}^{\sigma }\right)\right]}{{\left[{\sum }_{\sigma }{n}^{\sigma }\left(D+{I}^{\sigma }\right)\right]}^{2}},\end{eqnarray}$
is the temperature variation rate, and $\begin{eqnarray}E^{\prime} ={\omega }_{\mathrm{ov}}Q,\end{eqnarray}$
indicates the release rate of chemical heat that equals the energy variation rate due to the chemical reaction. In equation (In addition, the formula (42 ) is equivalent to the following matrix formA . From equation (42 ), the expressions of reaction terms can be obtained
$\begin{eqnarray}{\hat{{\boldsymbol{R}}}}^{\sigma }={{\boldsymbol{M}}}^{\sigma }{{\boldsymbol{R}}}^{\sigma },\end{eqnarray}$
where the elements of Rσ are given in appendix $\begin{eqnarray}{{\boldsymbol{R}}}^{\sigma }={\left({{\boldsymbol{M}}}^{\sigma }\right)}^{-1}{\hat{{\boldsymbol{R}}}}^{\sigma }.\end{eqnarray}$
As for the description of the chemical reactions, we can adopt the one-step reaction, two-step reaction, detailed or reduced multi-step chemical kinetics. For example, without loss of generality, three simple chemical reaction models are adopted in this manuscript.
It should be emphasized that, via the reaction terms on the right-hand side of the discrete Boltzmann equation (1 ), the multi-physics and chemical reactions are naturally coupled. The matrix inversion method is a precise and efficient calculation approach to compute the reaction terms in equation (47 ), because the sixteen moment relations of reaction terms are satisfied in an elegant way.
2.5. Nonequilibrium effects
In fact, the difference between equations (18 ) and (20 ) indicates the nonequilibrium departure degree of the physical system,49 ) means the kinetic moment deviation of a fluid component from its individual temporary equilibrium state, and equation (48 ) denotes the departure of a chemical species from the local mixing eventual equilibrium.
$\begin{eqnarray}{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{neq}}={\hat{{\boldsymbol{f}}}}^{\sigma }-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}={{\boldsymbol{M}}}^{\sigma }\left({{\boldsymbol{f}}}^{\sigma }-{{\boldsymbol{f}}}^{\sigma \mathrm{eq}}\right)={{\boldsymbol{M}}}^{\sigma }{{\boldsymbol{f}}}^{\sigma \mathrm{neq}},\end{eqnarray}$
with ${{\boldsymbol{f}}}^{\sigma \mathrm{neq}}={\left(\begin{array}{l}{f}_{1}^{\sigma \mathrm{neq}}\ {f}_{2}^{\sigma \mathrm{neq}}\ \cdots \ {f}_{N}^{\sigma \mathrm{neq}}\end{array}\right)}^{{\rm{T}}}$ and ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{neq}}\,={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma \mathrm{neq}}\ {\hat{f}}_{2}^{\sigma \mathrm{neq}}\ \cdots \ {\hat{f}}_{N}^{\sigma \mathrm{neq}}\end{array}\right)}^{{\rm{T}}}$. Similarly, we can define $\begin{eqnarray}{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{sneq}}={\hat{{\boldsymbol{f}}}}^{\sigma }-{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}={{\boldsymbol{M}}}^{\sigma }\left({{\boldsymbol{f}}}^{\sigma }-{{\boldsymbol{f}}}^{\sigma \mathrm{seq}}\right)={{\boldsymbol{M}}}^{\sigma }{{\boldsymbol{f}}}^{\sigma \mathrm{sneq}},\end{eqnarray}$
with ${{\boldsymbol{f}}}^{\sigma \mathrm{sneq}}={\left(\begin{array}{l}{f}_{1}^{\sigma \mathrm{sneq}}\ {f}_{2}^{\sigma \mathrm{sneq}}\ \cdots \ {f}_{N}^{\sigma \mathrm{sneq}}\end{array}\right)}^{{\rm{T}}}$ and ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{sneq}}\,={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma \mathrm{sneq}}\ {\hat{f}}_{2}^{\sigma \mathrm{sneq}}\ \cdots \ {\hat{f}}_{N}^{\sigma \mathrm{sneq}}\end{array}\right)}^{{\rm{T}}}$. Physically, equation (It is noteworthy that the physical quantities ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{neq}}$ and ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{sneq}}$ represent the nonequilibrium effects from various aspects. To be specific, ${\hat{f}}_{1}^{\sigma \mathrm{neq}}={\hat{f}}_{1}^{\sigma \mathrm{sneq}}\,=\,0$ in line with the mass conservation; ${\hat{f}}_{2}^{\sigma \mathrm{sneq}}\,=\,0$ and ${\hat{f}}_{3}^{\sigma \mathrm{sneq}}\,=\,0$ due to the momentum conservation, ${m}^{\sigma }{\hat{f}}_{2}^{\sigma \mathrm{neq}}={\rho }^{\sigma }({u}_{x}^{\sigma }-{u}_{x})$ and ${m}^{\sigma }{\hat{f}}_{3}^{\sigma \mathrm{neq}}={\rho }^{\sigma }({u}_{y}^{\sigma }-{u}_{y})$ denote the individual mass diffusion fluxes in the x and y directions, respectively, ${\hat{f}}_{4}^{\sigma \mathrm{sneq}}\,=\,0$ on account of the energy conservation. Neither ${\hat{f}}_{i}^{\sigma \mathrm{sneq}}$ nor ${\hat{f}}_{i}^{\sigma \mathrm{neq}}$ may be zero for i ≥ 5 when a system is in the nonequilibrium state.
Furthermore, equations (25 )–(27 ) are parts of the following nonequilibrium quantities,
$\begin{eqnarray}\begin{array}{l}{\hat{f}}_{5}^{\sigma \mathrm{sneq}}={{\rm{\Delta }}}_{5}^{\sigma }-\displaystyle \frac{{\rho }^{\sigma ^{\prime} }}{{S}_{5}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{{\left({u}_{x}^{\sigma }-{u}_{x}\right)}^{2}+{\left({u}_{y}^{\sigma }-{u}_{y}\right)}^{2}}{D+{I}^{\sigma }}\\ \ \ \ +\,\displaystyle \frac{{\rho }^{\sigma ^{\prime} }}{{S}_{5}^{1\sigma }{m}^{\sigma }}{\left({u}_{x}^{\sigma }-{u}_{x}\right)}^{2}\\ +\,\displaystyle \frac{{S}_{4}^{2\sigma }{\rho }^{\sigma }}{{S}_{5}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{{{u}^{\sigma }}^{2}-{u}^{2}}{D+{I}^{\sigma }}-\displaystyle \frac{{S}_{5}^{2\sigma }{\rho }^{\sigma }}{{S}_{5}^{1\sigma }{m}^{\sigma }}\left({{u}_{x}^{\sigma }}^{2}-{u}_{x}^{2}\right)\\ \ \ +\,\displaystyle \frac{{S}_{4}^{2\sigma }-{S}_{5}^{2\sigma }}{{S}_{5}^{1\sigma }}\displaystyle \frac{{\rho }^{\sigma }}{{{m}^{\sigma }}^{2}}\left({T}^{\sigma }-T\right)\\ +\,\displaystyle \frac{{S}_{2}^{2\sigma }{\rho }^{\sigma }}{{S}_{5}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{D+{I}^{\sigma }-1}{D+{I}^{\sigma }}2{u}_{x}^{\sigma }\left({u}_{x}^{\sigma }-{u}_{x}\right)\\ \ \ -\,\displaystyle \frac{{S}_{3}^{2\sigma }{\rho }^{\sigma }}{{S}_{5}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{2{u}_{y}^{\sigma }\left({u}_{y}^{\sigma }-{u}_{y}\right)}{D+{I}^{\sigma }},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{\hat{f}}_{6}^{\sigma \mathrm{sneq}}={{\rm{\Delta }}}_{6}^{\sigma }+\displaystyle \frac{{S}_{2}^{2\sigma }{\rho }^{\sigma }}{{S}_{6}^{1\sigma }{m}^{\sigma }}\left({u}_{x}^{\sigma }-{u}_{x}\right){u}_{y}^{\sigma }\\ \ \ +\,\displaystyle \frac{{S}_{3}^{2\sigma }{\rho }^{\sigma }}{{S}_{6}^{1\sigma }{m}^{\sigma }}{u}_{x}^{\sigma }\left({u}_{y}^{\sigma }-{u}_{y}\right)\\ \ \ -\,\displaystyle \frac{{S}_{6}^{2\sigma }{\rho }^{\sigma }}{{S}_{6}^{1\sigma }{m}^{\sigma }}\left({u}_{x}^{\sigma }{u}_{y}^{\sigma }-{u}_{x}{u}_{y}\right)\\ \ \ +\,\displaystyle \frac{{\rho }^{\sigma ^{\prime} }}{{S}_{6}^{1\sigma }{m}^{\sigma }}\left({u}_{x}{u}_{y}+{u}_{x}^{\sigma }{u}_{y}^{\sigma }-{u}_{x}{u}_{y}^{\sigma }-{u}_{x}^{\sigma }{u}_{y}\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{\hat{f}}_{7}^{\sigma \mathrm{sneq}}={{\rm{\Delta }}}_{7}^{\sigma }-\displaystyle \frac{{\rho }^{\sigma ^{\prime} }}{{S}_{7}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{{{u}^{\sigma }}^{2}+{u}^{2}-2{u}_{x}^{\sigma }{u}_{x}-2{u}_{y}^{\sigma }{u}_{y}}{D+{I}^{\sigma }}\\ \ \ +\,\displaystyle \frac{{\rho }^{\sigma ^{\prime} }}{{S}_{7}^{1\sigma }{m}^{\sigma }}{\left({u}_{y}^{\sigma }-{u}_{y}\right)}^{2}\\ +\,\displaystyle \frac{{S}_{4}^{2\sigma }{\rho }^{\sigma }}{{S}_{7}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{{{u}^{\sigma }}^{2}-{u}^{2}}{D+{I}^{\sigma }}-\displaystyle \frac{{S}_{7}^{2\sigma }{\rho }^{\sigma }}{{S}_{7}^{1\sigma }{m}^{\sigma }}\left({{u}_{y}^{\sigma }}^{2}-{u}_{y}^{2}\right)\\ \ \ +\,\displaystyle \frac{\left({S}_{4}^{2\sigma }-{S}_{7}^{2\sigma }\right){\rho }^{\sigma }}{{S}_{7}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{{T}^{\sigma }-T}{{m}^{\sigma }}\\ -\,\displaystyle \frac{{S}_{2}^{2\sigma }{\rho }^{\sigma }}{{S}_{7}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{2{u}_{x}^{\sigma }\left({u}_{x}^{\sigma }-{u}_{x}\right)}{D+{I}^{\sigma }}\\ \ \ +\,\displaystyle \frac{{S}_{3}^{2\sigma }{\rho }^{\sigma }}{{S}_{7}^{1\sigma }{m}^{\sigma }}\displaystyle \frac{D+{I}^{\sigma }-1}{D+{I}^{\sigma }}2{u}_{y}^{\sigma }\left({u}_{y}^{\sigma }-{u}_{y}\right),\end{array}\end{eqnarray}$
at the NS level, which can be proved via the CE expansion.Physically, both (${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{neq}}$, ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{sneq}}$) and (S1σ, S2σ) play roles in the thermodynamic and hydrodynamic behaviors, and various nonequilibrium modes are coupled in the relaxation process. Moreover, the nonequilibrium effects are important and traditional hydrodynamic models are not accurate in cases with small characteristic scales or large Knudsen numbers, particularly for multicomponent flows where various complex material and/or mechanical interfaces exist. For those complex nonequilibrium problems, the DBM provides a convenient tool to probe and analyze the nonequilibrium state and process.
2.6. Nondimensionalization
For numerical simulations and investigations, it is helpful to perform nondimensionalization. In this paper, Φd and Φn are designated as dimensional and nondimensional variables, respectively, and their ratio is Φr = Φd/Φn. The number density n, length L, flow speed $u=\left|{\boldsymbol{u}}\right|$, temperature T, and universal gas constant R are adopted as references, see table 1. Obviously, from these references, we can derive other ratios of dimensional to nondimensional variables, see appendix C . With the ratios and dimensional quantities, the nondimensional values are obtained in a straightforward way.
Table 1. References for nondimensionalization. |
| Variable | Dimension | Nondimension | Ratio |
|---|---|---|---|
| Number density | nd | nn | nr |
| Length | Ld | Ln | Lr |
| Flow speed | ud | un | ur |
| Temperature | Td | Tn | Tr |
| universal gas constant | R | 1 | R |
Additionally, various numerical schemes can be adopted to solve the discrete Boltzmann equations (1 ). In this paper, we employ the third-order total variation diminishing Runge–Kutta [56] for handling the time derivative and the fifth-order weighted essentially non-oscillatory scheme for the space derivative [57]. Consequently, in order to achieve good numerical stability, both the temporal and spatial steps must adhere to the Courant–Friedrichs–Lewy condition, and the temporal step should be smaller than the relaxation time. In the DBM with split collision, there are two sets of relaxation frequencies represented as ${{\boldsymbol{S}}}^{1\sigma }=\mathrm{diag}\left(\begin{array}{l}{S}_{1}^{1\sigma }\ {S}_{2}^{1\sigma }\ \cdots \ {S}_{N}^{1\sigma }\end{array}\right)$ and ${{\boldsymbol{S}}}^{2\sigma }=\mathrm{diag}\left(\begin{array}{l}{S}_{1}^{2\sigma }\ {S}_{2}^{2\sigma }\ \cdots \ {S}_{N}^{2\sigma }\end{array}\right)$, hence the temporal step Δt should be smaller than the minimum of ${\left({S}_{i}^{1\sigma }\right)}^{-1}$ and ${\left({S}_{i}^{2\sigma }\right)}^{-1}$. As for the DBM without split collision, there exists only one set of relaxation frequencies ${{\boldsymbol{S}}}^{\sigma }=\mathrm{diag}\left(\begin{array}{l}{S}_{1}^{\sigma }\ {S}_{2}^{\sigma }\ \cdots \ {S}_{N}^{\sigma }\end{array}\right)$, with the restriction ${\rm{\Delta }}t\lt {\left({S}_{i}^{\sigma }\right)}^{-1}$.
3. Verification and validation
To verify and validate the current model, six benchmarks are under consideration. First of all, multicomponent diffusion is adopted to confirm that the MRT DBM with the split collision term could provide a detailed relationship between the thermodynamic relaxation frequencies and the diffusivity of chemical species. Second, homogeneous mixture in a force field is used to test that the DBM can describe both thermal and isothermal systems where the acceleration and relaxation frequencies are tunable. Third, the KH instability is simulated to demonstrate that the DBM could capture fluid flows with complex interfacial structures. Fourth, the laminar flame of a propane-air mixture is chosen to demonstrate that the DBM is capable of mimicking combustion. Fifth, simulations of the opposing reactions are performed to verify that the DBM can capture the chemical nonequilibrium process accurately. Finally, steady detonation is simulated to validate that the DBM has the ability to capture the supersonic reactive front with a strong compressible effect.
3.1. Multicomponent diffusion
Diffusion [58, 59] widely exists in fields of physics, chemistry, biology, etc. It plays a vital role in the combustion process, particularly where non-premixed or partially mixed fuel and oxidant contact [1]. To mimic the diffusion in an accurate way is a prerequisite to simulating the combustion precisely. Actually, the diffusion (as well as viscosity and heat conduction) is a fundamental physical phenomenon in the thermodynamic nonequilibrium process. As the first benchmark, multicomponent diffusion is considered to validate our DBM for this kind of nonequilbrium process.
Give four chemical species σ = A, B, C, and D with the same molar mass mσ = 1 in the initial field as below
$\begin{eqnarray}\left\{\begin{array}{l}{\left({n}^{A},{n}^{B},{n}^{C},{n}^{D}\right)}_{L}=\left(1,2,3,4\right),\\ {\left({n}^{A},{n}^{B},{n}^{C},{n}^{D}\right)}_{R}=\left(4,3,2,1\right),\end{array}\right.\end{eqnarray}$
where the subscript L indicates 0 ≤ x ≤ 0.05, and R indicates 0.05 < x ≤ 0.1. The inflow and outflow boundary conditions are adopted in the x direction, and the periodic boundary conditions in the y direction. Moreover, the physical system is at rest (u = 0) and can be regarded as isothermal (T = 1). Consequently, the exact solutions of the concentrations read [58, 59], $\begin{eqnarray}{n}^{\sigma }=\displaystyle \frac{{n}_{L}^{\sigma }+{n}_{R}^{\sigma }}{2}-\displaystyle \frac{{n}_{L}^{\sigma }-{n}_{R}^{\sigma }}{2}\mathrm{Erf}\left(\displaystyle \frac{x-{x}_{0}}{\sqrt{4{\zeta }^{\sigma }t}}\right),\end{eqnarray}$
where $\mathrm{Erf}$ denotes the complementary error function, x0 = 0.05 the interfacial position, and the diffusivity $\begin{eqnarray}{\zeta }^{\sigma }=\displaystyle \frac{T}{{m}^{\sigma }{S}_{J}^{\sigma }},\end{eqnarray}$
in terms of ${S}_{J}^{\sigma }={S}_{2}^{2\sigma }={S}_{3}^{2\sigma }$.To begin with, let us perform a grid convergence analysis, which is of great importance for numerical simulations. For this sake, four simulations are conducted under various mesh grids Nx × Ny. To be specific, the mesh number is given as Nx = 10, 20, 40, and 80 in the horizontal direction and fixed as Ny = 1 in the vertical direction. The relaxation frequencies are chosen as ${S}_{i}^{1\sigma }={S}_{i}^{2\sigma }=1250$, the temporal step Δt = 2 × 10−5, the extra degrees of freedom Iσ = 3, and the parameters (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$, ${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$) = (0.01, 0.01, 1.75, 1.3, 3, 0, 2.1, 1). Figure 2 delineates the concentration of species A. The long-dashed, short-dashed, dash-dotted, and short-dotted lines denote simulated results in the four cases, and the solid line stands for the analytical solution in equation (54 ). As shown in figure 2, with increasing resolution, the numerical results approach the exact solution. That is to say, the differences between these simulation results and theoretical solutions decrease as the spatial step reduces. In particular, numerical results for both Nx = 40 and Nx = 80 are quite close to the solution, which is satisfying.
Figure 2. Grid convergence analysis: the horizontal distribution of concentration nA at the time t = 0.15. The solid line represents the exact solution, and the other lines indicate the simulation results under different mesh grids. |
Next, a comparison is made between the simulation results and the exact solutions at various time instants during the evolution of multicomponent diffusion. In order to obtain simulation results as accurately as possible, the mesh grid is chosen as Nx × Ny = 200 × 1. Other parameters are the same as those used in figure 2. Figure 3 displays the concentrations of chemical species A, B, C, and D along the x direction in the diffusion process. Symbols denote the DBM results, and solid lines indicate the corresponding exact solutions of equation (54 ). It can be observed that the simulated results match the exact solutions. Therefore, it is confirmed that the DBM can accurately describe the multicomponent diffusion.
Figure 3. Evolution of concentrations of four chemical species at time instants t1 = 0.01, t2 = 0.05, and t3 = 0.15, respectively. The left panel is for components A and D, and the right panel for B and C. Squares, circles, and triangles represent the simulation results at time instants t1 = 0.01, t2 = 0.05, and t3 = 0.15, respectively. Lines denote the corresponding exact solutions. |
Furthermore, as shown in equation (55 ), the diffusivity is a function of ${S}_{2}^{2\sigma }$ and ${S}_{3}^{2\sigma }$, but neither ${S}_{2}^{1\sigma }$ nor ${S}_{3}^{1\sigma }$. To verify it, a series of simulations are performed with different values of those parameters, i.e., ${S}_{2}^{1\sigma }={S}_{3}^{1\sigma }={S}_{J}^{1\sigma }$, and ${S}_{2}^{2\sigma }={S}_{3}^{2\sigma }={S}_{J}^{2\sigma }$. Figure 4(a) illustrates concentrations of species A versus x for cases of a fixed ${S}_{J}^{2\sigma }=1250$ and various ${S}_{J}^{1\sigma }=10000$, 5000, 2500, and 1250, respectively. Figure 4 (b) is for cases of a fixed ${S}_{J}^{1\sigma }=1250$ and various ${S}_{J}^{2\sigma }=10000$, 5000, 2500, and 1250, respectively. Obviously, the parameters ${S}_{2}^{1\sigma }$ and ${S}_{3}^{1\sigma }$ have a negligible impact on the diffusion process. The diffusivity depends upon the values of ${S}_{2}^{2\sigma }$ and ${S}_{3}^{2\sigma }$. The numerical simulations are entirely consistent with the exact solutions.
Figure 4. Concentrations of species A along the x direction at a time instant t = 0.15 in the diffusion process: (a) for a fixed ${S}_{J}^{2\sigma }$ and various ${S}_{J}^{1\sigma };$ (b) for a fixed ${S}_{J}^{1\sigma }$ and various ${S}_{J}^{2\sigma }$. Symbols denote the DBM results, and lines represent the corresponding exact solutions. |
It should be mentioned that the two-step-relaxation DBM presents the same results as the one-step-relaxation model for ${S}_{i}^{1\sigma }={S}_{i}^{2\sigma }$. It is verified that numerical results of the one-step-relaxation model are identical to those of the two-step-relaxation model when ${S}_{i}^{1\sigma }={S}_{i}^{2\sigma }$ in figure 4(a) or (b). Meanwhile, the simulations of the two-step-relaxation model in cases of ${S}_{i}^{1\sigma }\ne {S}_{i}^{2\sigma }$ are beyond the one-step-relaxation model. Compared to the latter one, the former model could present more details of nonequilibrium relaxation processes.
In addition, to further validate that the DBM can be used to measure thermodynamic nonequilibrium manifestations, figure 5 displays nonequilibrium quantities ${\hat{f}}_{5}^{\sigma \mathrm{sneq}}$ for species σ = A and B in the multicomponent diffusion process. The symbols stand for the numerical results, and the lines for the corresponding analytical solutions in equation (50 ). It can be observed that the DBM results coincide with the analytical solutions. The profiles of species C and D are similar to those in figure 5 and are not plotted here for brevity. As a result, it is demonstrated that the DBM is capable of describing nonequilibrium behaviors.
Figure 5. Nonequilibrium quantities ${\hat{f}}_{5}^{\sigma \mathrm{sneq}}$ in the diffusion process: (a) for species σ = A; (b) for species σ = B. Squares, circles, and triangles denote DBM results at time instants t1 = 0.01, t2 = 0.05, and t3 = 0.15, respectively. Solid lines represent the corresponding analytical solutions. |
3.2. Mixture in the force field
To confirm that chemical species with various accelerations and relaxation frequencies can be well described by using the current DBM, let us consider two situations, i.e., isothermal and thermal systems in force fields. Both the isothermal and thermal systems are homogeneous nonreactive mixtures that contain three components σ = A, B, and C, respectively. We choose the molar mass $\left({m}^{A},{m}^{B},{m}^{C}\right)=\left(2,1.5,1\right)$, molar concentrations $\left({n}^{A},{n}^{B},{n}^{C}\right)=\left(1,4,2\right)$, accelerations $\left({{\boldsymbol{a}}}^{A},{{\boldsymbol{a}}}^{B},{{\boldsymbol{a}}}^{C}\right)=\left(-{a}_{0},0,{a}_{0}\right){{\boldsymbol{e}}}_{x}$, extra degrees of freedom Iσ = 3, parameters (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$, ${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$) = (0.5, 1.5, 2.2, 3.5, 0, 5.2, 3, 0). Because the physical field is uniformly distributed, only one mesh grid Nx × Ny = 1 × 1 is used to have a high computing efficiency, and the periodic boundary conditions are adopted.
As for the isothermal mixtures, the initial temperatures and velocities are given as Tσ = 1 and uσ = 0, respectively. Theoretically, the temperatures and velocities are expressed by equations (D2 ) and (D4 ), when the systems reach steady states. To compare the DBM results with the theoretical solutions, figure 6 displays the velocities and temperatures as the mixtures are imposed on various accelerations and the relaxation frequencies are chosen as ${S}_{2}^{1\sigma }={S}_{2}^{2\sigma }=2000$. The squares, circles, triangles, and diamonds indicate the simulation results of the mixing system, chemcial species A, B, and C, respectively. The solid lines stand for the corresponding exact results. Obviously, all DBM results coincide exactly with the theoretical solutions.
Figure 6. Horizontal velocities (a) and temperatures (b) with different accelerations when the isothermal systems reach steady states. Symbols denote the DBM results, and solid lines represent the exact solutions. |
Let us now examine the impact of varying relaxation frequencies. In figure 7, we observe the horizontal velocities and temperatures when employing variable ${S}_{i}^{1\sigma }$ alongside a fixed ${S}_{i}^{2\sigma }=2000$ within the stationary isothermal system. Meanwhile, figure 8 showcases scenarios involving a constant ${S}_{i}^{1\sigma }=2000$ alongside a variable ${S}_{i}^{2\sigma }$. Numerical results are denoted by symbols, while theoretical outcomes are represented by lines. It is apparent that the physical fields remain unaffected by ${S}_{i}^{1\sigma }$ and are solely dependent on ${S}_{i}^{2\sigma }$ within this homogeneous mixture under the specified body force. Notably, the DBM consistently yields exact results across all cases, which is quite satisfying.
Figure 7. Horizontal velocities (a) and temperatures (b) with variable ${S}_{i}^{1\sigma }$ and fixed ${S}_{i}^{2\sigma }$ in the stationary isothermal system. Symbols denote the DBM results, and solid lines represent the corresponding exact solutions. |
Figure 8. Horizontal velocities (a) and temperatures (b) with fixed ${S}_{i}^{1\sigma }$ and changeable ${S}_{i}^{2\sigma }$ in the stationary isothermal system. Symbols denote the DBM results, and solid lines represent the corresponding exact solutions. |
Next, the thermal mixtures are under consideration. The ultimate steady physical fields of the isothermal mixtures in figure 6 are set as the initial configurations of the thermal systems. Figure 9 plots the evolution of temperatures when the thermal mixtures are in force fields. Four values of accelerations a0 = 10, 30, 50, and 70 are under consideration. The squares, circles, triangles, and diamonds represent the simulated temperatures of the mixing system, chemcial species A, B, and C, respectively. The solid lines are for the exact results of the mixing temperatures. It is clear that the DBM results are in good agreement with the exact solutions.
Figure 9. Evolution of temperatures with different accelerations. Symbols denote the DBM results, and solid lines represent the exact solutions. |
Figure 10 delineates the velocities and temperatures at the time t = 0.15 in the evolution of thermal systems with various accelerations. Note that, the velocities of thermal mixtures in figure 10(a) are exactly the same as those in figure 6(a). Besides, the simulated mixing temperatures agree well with the exact solutions, and the individual and mixing temperatures are close to each other in the thermal systems, as shown in figure 10(b). In fact, there are only minor differences among the individual and mixing temperatures in the isothermal systems as well, see figure 6.
Figure 10. Horizontal velocities (a) and temperatures (b) with different accelerations when the thermal systems are at the moment t = 0.15. Symbols denote the DBM results, and solid lines represent the exact solutions. |
It should be mentioned that the simulation results in figures 6–10 are calculated with Method I. Actually, Methods II and III give similar results. To compare Methods I, II, and III, table 2 lists the data in the simulation case of a0 = 10 in figure 10. Simulations are performed on a personal computer with Intel(R) Core(TM) i9-9880H CPU @ 2.30 GHz, RAM 64.0 GB, and a 64-bit version system for double precision floating point operations. As shown in table 2, the simulated density and velocities equal the corresponding exact solutions. It means that the conservation of mass and momentum is obeyed for all methods. Besides, the mixing temperatures simulated with Methods I, II, and III are T=1.01714857, 1.02572000, and 1.01714857, respectively. Compared with the exact value 1.01714286, the relative errors are 5.6 × 10−6, 8.4 × 10−3, and 5.6 × 10−6, respectively. That is to say, Methods I and III give the same simulation results, both have a higher accuracy than Method II. Additionally, the computing time is 54, 60, and 50 seconds for the three methods as the simulation runs 2 × 105 time steps. In other words, Methods I and II require 8% and 20% more computing time than Method III. Because it needs to calculate the equilibrium distribution functions fσseq once per loop for Method I, twice per loop for Method II, while it does not demand fσseq for Method III. Consequently, it takes less computational cost and running time for Method III.
Table 2. Simulation data about the three methods to calculate force terms. |
| Methods | ρ | ux | $\ {u}_{x}^{A}\ $ | $\ {u}_{x}^{B}\ $ | $\ {u}_{x}^{C}\ $ | T | Computing time |
| Method I | 10 | 0 | −0.005 | 0 | 0.005 | 1.01714857 | 54 s |
| Method II | 10 | 0 | −0.005 | 0 | 0.005 | 1.02572000 | 60 s |
| Method III | 10 | 0 | −0.005 | 0 | 0.005 | 1.01714857 | 50 s |
| Exact | 10 | 0 | −0.005 | 0 | 0.005 | 1.01714286 | / |
3.3. Kelvin–Helmholtz instability
As a fundamental interfacial instability in fluid mechanics, the KH instability occurs when there is velocity shear across a wrinkled interface in a fluid system, and leads to the formation of vortices and turbulence [60]. The KH instability is ubiquitous in nature and of considerable interest in scientific and engineering fields [60–62]. To show the capacity of our DBM in dealing with complex fluid systems, here we simulate the KH instability with a complicated interfacial dynamics.
Figure 11 portrays the sketch of the initial configuration of a three-component fluid. The length and height of the calculation domain are Lx = Ly = 1. The domain is divided into three parts, i.e., 0 < x ≤ x1, x1 < x ≤ x2, and x2 < x ≤ Lx. Between the left and middle parts is the interface located at x1 = 0.3Lx, between the middle and right parts is the interface located at x2 = 0.7Lx. To trigger the KH instability rollup, both interfaces have an imposed cosinusoid perturbation, $w={w}_{0}\cos \left(4\pi y/{L}_{y}\right)$, with an amplitude w0 = Lx/200. Initially, the left (middle, right) part is occupied by species A (B, C) moving upwards (downwards, upwards) with velocity uL = u0ey (uM = − u0ey, uR = u0ey). Both concentrations and temperatures in the three parts are equal, i.e., nL = nM = nR and TL = TM = TR, hence the pressure is homogeneous across the two interfaces, pL = pM = pR, due to the constitutive relation pσ = nσTσ. To be smooth across the interface, the initial profiles of the concentrations and velocities are given by
$\begin{eqnarray*}{n}^{A}=\displaystyle \frac{{n}_{L}}{2}-\displaystyle \frac{{n}_{L}}{2}\tanh \left(\displaystyle \frac{x-{x}_{1}+w}{{W}_{n}}\right),\end{eqnarray*}$
$\begin{eqnarray*}{n}^{C}=\displaystyle \frac{{n}_{R}}{2}+\displaystyle \frac{{n}_{R}}{2}\tanh \left(\displaystyle \frac{x-{x}_{2}+w}{{W}_{n}}\right),\end{eqnarray*}$
$\begin{eqnarray*}{n}^{B}={n}_{M}-{n}^{A}-{n}^{C},\end{eqnarray*}$
$\begin{eqnarray*}{{\boldsymbol{u}}}^{\sigma }=\left\{\begin{array}{ll}\tfrac{{{\boldsymbol{u}}}_{L}+{{\boldsymbol{u}}}_{M}}{2}-\tfrac{{{\boldsymbol{u}}}_{L}-{{\boldsymbol{u}}}_{M}}{2}\tanh \left(\tfrac{x-{x}_{1}+w}{{W}_{{\boldsymbol{u}}}}\right), & \mathrm{for}\ 0\lt x\leqslant \tfrac{{L}_{x}}{2},\\ \tfrac{{{\boldsymbol{u}}}_{M}+{{\boldsymbol{u}}}_{R}}{2}-\tfrac{{{\boldsymbol{u}}}_{M}-{{\boldsymbol{u}}}_{R}}{2}\tanh \left(\tfrac{x-{x}_{2}+w}{{W}_{{\boldsymbol{u}}}}\right), & \mathrm{for}\,\tfrac{{L}_{x}}{2}\lt x\leqslant {L}_{x},\end{array}\right.\end{eqnarray*}$
where Wn (Wu) indicates the width of concentration (velocity) transition layer.
Figure 11. Initial configuration of the KH instability. |
The boundary conditions are as follows: specular reflection boundary condition in the x direction and periodic boundary condition in the y direction. A simulation is carried out on a uniform mesh Nx × Ny = 2000 × 2000 with Δx = Δy = 5 × 10−4. The time step is set to be as small as Δt = 2.5 × 10−5 to keep the numerical dissipation negligible. The other parameters are mA = 1, mB = 1.5, mC = 2, Iσ = 3, ${S}_{i}^{1\sigma }={S}_{i}^{2\sigma }=5\times {10}^{3}$, and (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$, ${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$) = (2, 1.414, 3.9, 2.758, 1.5, 0, 5.5, 0). In addition, as the number of mesh grids is large enough, the parallel programming with the message-passing interface is implemented in Fortran to improve the computing speed. Actually, because all information transfer is local in time and space in the evolution of the discrete Boltzmann equation, the DBM has natural parallelism with excellent scalability [38].
Figure 12 depicts contours of the molar fraction of species B at representative times in the evolution of the miscible KH instability. It can be observed in figures 12(a)–(b) that the amplitude of perturbation w increases due to the shear effect and the width of concentration transition layer Wn increases because of diffusion. The fluid interface begins wiggling and its shape changes from regular to irregular gradually. It can be found in figure 12(c) that, as time advances, several pairs of vortices appear, and the middle fluid penetrates into the left and right ones. Figures 12(d)–(e) show that the continuous growth of vortices leads to formation of billows, and nonregular interfaces become more complex at the later stage. As shown in figure 12(f), fluid structures are chaotic and the KH instability promotes the mixture between different fluid species. The above dynamic process of the KH instability obtained by our model is basically consistent with the scenarios in previous studies [61, 62].
Figure 12. Molar fraction XB at time instants t = 0.0, 0.5, 1.0, 1.5, 2.0, and 3.0 in the evolution of the KH instability. The color from blue to red corresponds to the value from 0 to 1. |
To further verify our model, we measure conserved quantities in the process of KH instability, see figure 13. Squares, pentagrams, circles, and triangles indicate numerical results of the mass $\iint $ρdxdy, momentum in the x direction $\iint $Jxdxdy, momentum in the y direction $\iint $Jydxdy, and energy $\iint $Edxdy, respectively. Solid lines are for the corresponding exact solutions $\iint $ρdxdy = 1.5, $\iint $Jxdxdy = 0, $\iint $Jydxdy = 0.15, and $\iint $Edxdy = 2.685, respectively. It is clear in figure 13 that our computed results agree well with these exact solutions. Moreover, it is found that the numerical results of $\iint $ρσdxdy are exactly equal to their corresponding theoretical solutions 0.3, 0.6, and 0.6 for species σ = A, B, and C, respectively. The results are very good and satisfactory.
Figure 13. Conserved quantities ($\iint $ρdxdy, $\iint $Jxdxdy, $\iint $Jydxdy, $\iint $Edxdy) in the evolution of the KH instability. Symbols and lines indicate numerical and exact results, respectively. |
3.4. Laminar flame
In this subsection, the objective is to demonstrate that the DBM is suitable for combustion. For this purpose, we simulate a laminar flame of propane-air mixture. The combustion is controlled by the one-step overall reaction,
$\begin{eqnarray}{{\rm{C}}}_{3}{{\rm{H}}}_{8}+5{{\rm{O}}}_{2}\to 3{\mathrm{CO}}_{2}+4{{\rm{H}}}_{2}{\rm{O}},\end{eqnarray}$
$\begin{eqnarray}{\omega }^{\sigma }={s}^{\sigma }\cdot {m}^{\sigma }\cdot {\omega }_{\mathrm{ov}},\end{eqnarray}$
$\begin{eqnarray}{\omega }_{\mathrm{ov}}={k}_{\mathrm{ov}}{n}^{{{\rm{C}}}_{3}{{\rm{H}}}_{8}}{n}^{{{\rm{O}}}_{2}}\exp \left(-{E}_{{\rm{a}}}/{RT}\right),\end{eqnarray}$
where C3H8, O2, CO2, and H2O stand for propane, oxygen, carbon dioxide, and water, respectively. Nitrogen N2 is assumed to be inert. The stoichiometric coefficients are $[{s}^{{{\rm{C}}}_{3}{{\rm{H}}}_{8}}$, ${s}^{{{\rm{O}}}_{2}}$, ${s}^{{\mathrm{CO}}_{2}}$, ${s}^{{{\rm{H}}}_{2}{\rm{O}}}$, ${s}^{{{\rm{N}}}_{2}}]$ = [−1, −5, 3, 4, 0], the molar mass $[{m}^{{{\rm{C}}}_{3}{{\rm{H}}}_{8}}$, ${m}^{{{\rm{O}}}_{2}}$, ${m}^{{\mathrm{CO}}_{2}}$, ${m}^{{{\rm{H}}}_{2}{\rm{O}}}$, ${m}^{{{\rm{N}}}_{2}}]$ = [4.4, 3.2, 4.4, 1.8, 2.8] ×10−2 [kg/mol], the reaction coefficient ${k}_{\mathrm{ov}}=9.9\,\times {10}^{7}\left[{{\rm{m}}}^{3}\cdot {\mathrm{mol}}^{-1}\cdot {{\rm{s}}}^{-1}\right]$, the universal gas constant $R\,=8.315\left[{\rm{J}}\cdot {\mathrm{mol}}^{-1}\cdot {{\rm{K}}}^{-1}\right]$, the effective activation energy ${E}_{{\rm{a}}}=1.26\times {10}^{5}\left[{\rm{J}}\cdot {\mathrm{mol}}^{-1}\right]$, the chemical heat of overall reaction $Q=2.05\times {10}^{6}\left[{\rm{J}}\cdot {\mathrm{mol}}^{-1}\right]$, the overall reaction rate ωov, the mass change rate of species ωσ, and the parameters (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$, ${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$) = (0.7, 0.7, 3.9, 2.758, 1.5, 0, 5.5, 0). Initially, a channel with length L = 4 cm is filled with the propane-air mixture with an equivalence ratio 0.6. The molar concentration is 44.6 mol · m−3, the temperature 300 K, the pressure 1 atm. After ignition in the left part of the channel 0 ≤ x ≤ L/64, the flame starts to move downstream. The periodic boundary conditions are used at the upper and lower walls, the specular reflection (outflow) boundary condition at the left (right) side. The grid is chosen as Nx × Ny = 1600 × 1, the spatial step Δx = Δy = 2.5 × 10−5 m, and the temporal step Δt = 1.25 × 10−10 s.Figure 14 delineates the overall reaction rate at various time instants in the evolution of the laminar flame. Clearly, the overall reaction rate first increases then decreases and forms a peak at the combustion front as the flame propagates forwards. The profile of the reaction rate gradually becomes steady as time goes on. Let us define the flame position as the location where the reaction takes its maximum. Then we can obtain the temporal evolution of the flame position and its velocity, as shown in figure 15. It is clear that the flame moves rightwards after the initial ignition stage, and its speed tends to be a constant in the later period. To be specific, the flame speed is around Uf = 0.522 m/s at t = 0.05 s. Meanwhile, the value of flow velocity is about Uo = 0.413 m/s ahead of the moving flame. Note that the burning velocity can be estimated by the relation, Uf = Uo + Ub, hence the resultant burning velocity is Ub = 0.109 m/s.
Figure 14. The overall reaction rate in the evolution of the laminar flame at time instants t1 = 3.125 × 10−2 s, t2 = 3.75 × 10−2 s, t3 = 4.375 × 10−2 s, and t4 = 5 × 10−2 s from left to right. |
Figure 15. Temporal evolution of the flame position and speed presented by the solid and dashed lines, respectively. |
Figure 16 exhibits the evolution of the burning speed. To give a clearer depiction, an insert is attached for the enlargement within the period 0.045 s ≤ t ≤ 0.05 s, during which the flame speed is almost steady with only small perturbations. It is easy to calculate the average burning speed, Ub = 0.10941 m/s, within the time range 0.045 s ≤ t ≤ 0.05 s. Clearly, the numerical result approaches the experimental datum 0.11 m/s in [63]. Furthermore, figure 17 displays the temperature profiles at t = 0.05 s. Specifically, the temperatures are 1704 K and 1705 K in the simulation and experiment [1], respectively. The satisfying result of the DBM is due to its outstanding advantages: (i) physical properties such as the extra degrees of freedom and the specific heat ratio are flexible in the DBM; (ii) the physical fields of concentration, velocity, temperature, and pressure are naturally coupled together in the DBM; (iii) the DBM is suitable for both low-speed incompressible and high-speed compressible fluid flows.
Figure 16. The burning speed in the evolution of the laminar flame. The insert shows the zoom-in view within the range 0.045 s ≤ t ≤ 0.05 s. The dashed and solid lines denote simulation results and their average, and the squares represent experimental data. |
Figure 17. Temperature profiles in the evolution of the laminar flame at time instant t = 0.05 s. The line and squares stand for the DBM results and experimental data [1], respectively. |
3.5. Opposing reaction
To demonstrate that our DBM is suitable for a chemical nonequilibrium system, the opposing reactions A ⇌ B are used as a typical benchmark. The forward and reverse reaction coefficients are k1 and k−1, respectively. Hence the overall reaction rate reads ωov = k1nA − k−1nB. Initially, the concentrations are (nA, nB)=(n0, 0) = (1, 0), the velocity u = 0, the temperature T = 1. The concentrations are
$\begin{eqnarray*}{n}^{A}={n}_{e}^{A}+{n}_{e}^{B}\exp \left(-\displaystyle \frac{{k}_{1}{n}_{0}}{{n}_{e}^{B}}t\right),\end{eqnarray*}$
$\begin{eqnarray*}{n}^{B}={n}_{e}^{B}-{n}_{e}^{B}\exp \left(-\displaystyle \frac{{k}_{1}{n}_{0}}{{n}_{e}^{B}}t\right),\end{eqnarray*}$
during the nonequilibrium reaction process, and $\begin{eqnarray*}{n}_{e}^{A}=\displaystyle \frac{{k}_{-1}{n}_{0}}{{k}_{1}+{k}_{-1}},\end{eqnarray*}$
$\begin{eqnarray*}{n}_{e}^{B}=\displaystyle \frac{{k}_{1}{n}_{0}}{{k}_{1}+{k}_{-1}},\end{eqnarray*}$
when the chemical reaction reaches equilibrium. Given the heat release Q = 10 and extra degrees of freedom IA = IB = I = 3, the temperature reads $\begin{eqnarray*}T={T}_{0}+\displaystyle \frac{2{n}^{B}Q}{{n}_{0}\left(D+I\right)},\end{eqnarray*}$
as the chemical reaction takes place.Figure 18 illustrates physical quantities (nA, nB, ρ, ux, T, ωov) versus time in the case of (k1, k−1) = (0.3, 0.7). The relaxation frequencies are choosen as ${S}_{i}^{1\sigma }={S}_{i}^{2\sigma }={10}^{4}$, and the parameters (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$, ${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$) = (0.7, 0.7, 2.7, 2.7, 0, 5.3, 5, 0). To be specific, the concentration of species A (B) decreases (increases) gradually as time goes on. The mixture density remains ρ = 1, which equals the exact solution ρ = mAnA + mBnB with mA = mB = 1. The system remains motionless, i.e., ux = uy = 0, in line with the momentum conservation. The temperature rises and the reaction rate reduces over time. On the whole, the DBM results coincide well with the exact solutions in the evolution of chemical reaction. All quantities become constant when the opposing reaction reaches equilibrium. At t = 20, the numerical results (nA, nB, ρ, ux, T, ωov) = (0.7, 0.3, 1, 0, 2.2, 0) match the exact solutions. Consequently, it is verified that the DBM is capable of both chemical equilibrium and nonequilibrium processes.
Figure 18. Physical quantities versus time in the evolution of opposing reaction. Symbols denote DBM results in the legend, and lines denote corresponding exact solutions. |
Let us consider more cases of opposing reactions with various forward reaction coefficients k1. For simplicity, the corresponding reverse reaction coefficients are set as k−1 = 1 − k1. Figure 19 displays the concentrations of species A and B after the chemical reaction reaches equilibrium. It is clear that the concentration of species A (B) reduces (increases) linearly with the increasing forward reaction coefficient. There is a satisfying agreement between the simulation results and exact solutions.
Figure 19. Concentrations of species A and B versus the reaction coefficient. Symbols represent DBM results in the legend, and lines represent corresponding exact solutions. |
3.6. Detonation wave
Detonation is a particular type of combustion with violent chemical heat release around a supersonic exothermic front accelerating through a medium. The physical fields have strong temporal and spatial changes near the detonation wave, which poses a great challenge to the numerical robustness and physical accuracy of computational fluid dynamics. In this subsection, we demonstrate that the DBM has the capability of capturing the detonation wave travelling at a supersonic speed. As the detonation wave passes in the x direction, the chemical reactant A changes into the product B, i.e., A → B, and the chemical energy is released. In theory, the speed of the steady detonation front is a function of the chemical heat release of reactant per unit mass, q, i.e.,
$\begin{eqnarray}{D}_{s}=\sqrt{\displaystyle \frac{\left({\gamma }^{2}-1)q\right.}{2}+\gamma {T}_{0}}+\sqrt{\displaystyle \frac{\left({\gamma }^{2}-1)q\right.}{2}},\end{eqnarray}$
where T0 denotes the temperature in front of the detonation wave and γ represents the specific heat ratio. Then, the Mach number Ma is calculated by $\begin{eqnarray}\mathrm{Ma}=\displaystyle \frac{{D}_{s}}{\sqrt{\gamma {T}_{0}}}.\end{eqnarray}$
Here, we focus on a specific case where q = 1, yielding a corresponding Mach number of Ma = 1.74436. It is worth noting that the DBM excels in simulating detonations under a high Mach number, employing a robust numerical scheme.The reaction rate is controlled by
$\begin{eqnarray}{\omega }_{\mathrm{ov}}={k}_{\mathrm{ov}}{n}^{A}\exp \left(-\displaystyle \frac{{E}_{{\rm{a}}}}{{RT}}\right),\end{eqnarray}$
in terms of kov = 5 × 105 and Ea = 10. The initial configuration is set as $\begin{eqnarray}\left\{\begin{array}{l}{\left({n}^{A},{n}^{B},{u}_{x},T\right)}_{L}=\left(0,1.38837,0.57735,1.57856\right),\\ {\left({n}^{A},{n}^{B},{u}_{x},T\right)}_{R}=\left(1,0,0,1\right),\end{array}\right.\end{eqnarray}$
where the subscript L indicates 0 ≤ x ≤ 0.04, and R indicates 0.04 < x ≤ 0.4. The quantities in the left and right parts satisfy the Hugoniot relationship for detonation wave. The parameters are mσ = 1, γσ = 1.5, ${S}_{i}^{1\sigma }={S}_{i}^{2\sigma }=2\times {10}^{4}$, and (${v}_{a}^{\sigma }$, ${v}_{b}^{\sigma }$, ${v}_{c}^{\sigma }$, ${v}_{d}^{\sigma }$, ${\eta }_{a}^{\sigma }$, ${\eta }_{b}^{\sigma }$, ${\eta }_{c}^{\sigma }$, ${\eta }_{d}^{\sigma }$) = (0.5, 1.5, 2.2, 3.5, 0, 5.2, 3, 0). In addition, inflow and outflow boundary conditions are adopted in the x direction, and the periodic boundary conditions in the y direction. Figure 20 depicts the evolution of pressure profiles around the detonation wave that propagates forwards. It is clear that the spatial distribution of the pressure field is quite similar to each other at the four time instants t = 0.075, 0.1, 0.125, and 0.15. It indicates that the detonation wave moves forwards in a steady state. Then, the speed of the steady detonation front can be calculated, Ds = 2.058. Compared with the analytic solution 2.06395, the relative error is about 0.0029. The simulation results are satisfactory.
Figure 20. Pressure profiles at time instants t1 = 0.075, t2 = 0.1, t3 = 0.125, and t4 = 0.15 in the evolution of detonation. |
Figure 21 gives the density (a), temperature (b), pressure (c), and horizontal velocity (d) at t = 0.15 in the detonation process. The squares represent the DBM results, the triangles denote the numerical outcomes obtained from a Euler solver, and the lines stand for the analytical solutions of the Zel'dovich–Neumann–Döring (ZND) results in theory [1]. Obviously, there is a satisfying agreement among the three models in regions distant from the detonation wave. To be specific, the DBM results are (ρ, T, p, ux) = (1.38748, 1.57724, 2.18838, 0.57892) after the detonation wave, resulting in relative errors of (0.0006, 0.0008, 0.0015, 0.0027) compared to the ZND results [1]. However, slight differences between them emerge at the von-Neumann-peak. This disparity arises from the fact that both the ZND theory and Euler solver neglect the viscosity and heat conduction, assuming a sharp discontinuity at the von-Neumann-peak. In contrast, the DBM takes into account the viscosity, heat conduction, and other thermodynamic nonequilibrium effects. As a result, the physical fields simulated by the DBM exhibit smoothness across the detonation wave, aligning more closely with real-world conditions.
Figure 21. Physical quantities around the detonation wave. The squares, triangles, and lines represent the results obatined from the DBM, Euler solver, and ZND theory [1], respectively. |
It should be mentioned that previous research on detonation mechanisms has primarily relied on traditional computational fluid dynamics methods, which differ substantially from the current DBM. Physically, this DBM can be likened to a modified continuous fluid model augmented with a coarse-grained representation of significant thermodynamic nonequilibrium effects. Consequently, the DBM exhibits the capability to accurately capture detonation phenomena with nonequilibrium effects, encompassing diffusion, viscosity, and thermal conduction. In fact, the DBM is suitable for both steady and unsteady detonation scenarios, although the latter is not explicitly demonstrated in this manuscript due to space constraints.
4. Conclusion
An MRT DBM with split collision is presented for both subsonic and supersonic reactive flows. The external forces, chemical reactions, and multi-physical fields are coupled naturally through the collision, force and reaction terms on the right-hand side of the discrete Boltzmann equations that describe the evolution of reactive mixture. Through the CE expansion, it can be proved that the DBM is consistent with the reactive NS equations with external forces, the Fick's law and Stefan–Maxwell diffusion equation in the hydrodynamic limit. Each chemical species owns individual adjustable molar mass, concentration, velocity, acceleration, temperature, pressure, diffusivity, dynamic viscosity, thermal conductivity, specific heat ratio, Prandtl number, Reynolds number, and Schmidt number, etc.
Compared to the one-step-relaxation MRT or BGK model [11, 13], the DBM with the splitting technique, a two-step-relaxation model, offers greater flexibility in parameters and is applicable to a broader range of physicochemical systems. (i) The relaxation frequencies S1σ and S2σ govern the thermodynamic nonequilibrium process, guiding the approach toward temporary individual equilibrium and ultimate mixing equilibrium. (ii) The relaxation frequencies in the split collision term influence both self- and cross-collisions, thereby impacting the evolution of the mixture. (iii) Specific relationships can be established between the relaxation frequencies and other physical quantities, such as thermodynamic nonequilibrium quantities, diffusivity, dynamic viscosity, and thermal conductivity. (iv) Classical dimensionless numbers in fluid mechanics, such as the Reynolds number, Prandtl number, and Schmidt number of each species, can be adjusted. Consequently, the DBM with the split collision term presents a more detailed relationship between the thermodynamic relaxation process and nonequilibrium effects. Moreover, the two-step-relaxation DBM can reduce to the one-step-relaxation or BGK model under special conditions.
It should be stressed that the hydrodynamic, thermodynamic, and chemical nonequilibrium effects can be captured and measured by the versatile kinetic DBM dynamically. Physically, the DBM is more general than traditional NS solvers since it contains more detailed thermodynamic nonequilibrium information. Mathematically, a set of uniform discrete Boltzmann equations is used to describe the reactive mixtures, and the algorithm is easy to code due to the linearization of evolution equations. Computationally, it can be implemented on very large parallel clusters with exceptional scalability because all information transfer in DBM is local in time and space, which is similar to other LBMs.
In addition, three methods to calculate the source terms (including the force and reaction terms) are introduced into the multicomponent DBM. As a traditional idea, Method I uses the discretization form of the formula of source terms in the velocity space directly. Method II expresses the source terms as the change of discrete distribution functions due to the source influences over a small time interval. Method III gives the expression of source terms by using the matrix inversion method. Methods I and III own higher accuracy than Method II which possesses only the first-order accuracy. Besides, Method III has the highest computational efficiency because it requires to calculate the equilibrium distribution functions zero, one and two times per loop for Methods III, I and II, respectively.
Finally, several canonical systems, including the multicomponent diffusion, mixture in the force field, KH instability, laminar flame of propane-air mixture, opposing reactions, and detonation wave are simulated to validate this model. The first three benchmarks are physical systems without chemical reaction, and the last three benchmarks have chemical reactions, with the last one containing rather violent chemical heat release. It is demonstrated that the current DBM is suitable for multicomponent mixtures with or without the chemical reaction. The interplay among different chemical species can be described accurately. The complex interfacial structures could be captured dynamically. The essential nonequilibrium and compressible effects can be quantified. In the near future, this DBM will be employed to investigate more practical combustion problems with significant nonequilibrium and compressible effects.
Appendix A. Matrices
The square matrix Mσ takes the form,
$\begin{eqnarray}{{\boldsymbol{M}}}^{\sigma }=\left(\begin{array}{cccc}{M}_{11}^{\sigma } & {M}_{12}^{\sigma } & \cdots & {M}_{1N}^{\sigma }\\ {M}_{21}^{\sigma } & {M}_{22}^{\sigma } & \cdots & {M}_{2N}^{\sigma }\\ \vdots & \vdots & \ddots & \vdots \\ {M}_{N1}^{\sigma } & {M}_{N2}^{\sigma } & \cdots & {M}_{{NN}}^{\sigma }\end{array}\right),\end{eqnarray}$
with ${M}_{1i}^{\sigma }=1$, ${M}_{2i}^{\sigma }={v}_{{ix}}^{\sigma }$, ${M}_{3i}^{\sigma }={v}_{{iy}}^{\sigma }$, ${M}_{4i}^{\sigma }={{v}_{i}^{\sigma }}^{2}+{{\eta }_{i}^{\sigma }}^{2}$, ${M}_{5i}^{\sigma }\,={{v}_{{ix}}^{\sigma }}^{2}$, ${M}_{6i}^{\sigma }={v}_{{ix}}^{\sigma }{v}_{{iy}}^{\sigma }$, ${M}_{7i}^{\sigma }={{v}_{{iy}}^{\sigma }}^{2}$, ${M}_{8i}^{\sigma }=\left({{v}_{i}^{\sigma }}^{2}+{{\eta }_{i}^{\sigma }}^{2}\right)\times \,{v}_{{ix}}^{\sigma }$, ${M}_{9i}^{\sigma }\,=\left({{v}_{i}^{\sigma }}^{2}+{{\eta }_{i}^{\sigma }}^{2}\right){v}_{{iy}}^{\sigma }$, ${M}_{10i}^{\sigma }={{v}_{{ix}}^{\sigma }}^{3}$, ${M}_{11i}^{\sigma }={{v}_{{ix}}^{\sigma }}^{2}{v}_{{iy}}^{\sigma }$, ${M}_{12i}^{\sigma }={v}_{{ix}}^{\sigma }{{v}_{{iy}}^{\sigma }}^{2}$, ${M}_{13i}^{\sigma }={{v}_{{iy}}^{\sigma }}^{3}$, ${M}_{14i}^{\sigma }=\left({{v}_{i}^{\sigma }}^{2}+{{\eta }_{i}^{\sigma }}^{2}\right){{v}_{{ix}}^{\sigma }}^{2}$, ${M}_{15i}^{\sigma }=\left({{v}_{i}^{\sigma }}^{2}+{{\eta }_{i}^{\sigma }}^{2}\right){v}_{{ix}}^{\sigma }{v}_{{iy}}^{\sigma }$, and ${M}_{16i}^{\sigma }=\left({{v}_{i}^{\sigma }}^{2}+{{\eta }_{i}^{\sigma }}^{2}\right){{v}_{{iy}}^{\sigma }}^{2}$. Its inverse ${({{\boldsymbol{M}}}^{\sigma })}^{-1}$ could be obtained by using a calculation software, such as Matlab or Mathematica.The column matrix ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}$ is given by
$\begin{eqnarray}{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{eq}}={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma \mathrm{eq}}\ {\hat{f}}_{2}^{\sigma \mathrm{eq}}\ \cdots \ {\hat{f}}_{N}^{\sigma \mathrm{eq}}\end{array}\right)}^{{\rm{T}}},\end{eqnarray}$
in terms of ${\hat{f}}_{1}^{\sigma \mathrm{eq}}={n}^{\sigma }$, ${\hat{f}}_{2}^{\sigma \mathrm{eq}}={n}^{\sigma }{u}_{x}$, ${\hat{f}}_{3}^{\sigma \mathrm{eq}}={n}^{\sigma }{u}_{y}$, ${\hat{f}}_{4}^{\sigma \mathrm{eq}}={n}^{\sigma }[(D+{I}^{\sigma })T/{m}^{\sigma }+{u}^{2}]$, ${\hat{f}}_{5}^{\sigma \mathrm{eq}}={n}^{\sigma }(T/{m}^{\sigma }+{u}_{x}^{2})$, ${\hat{f}}_{6}^{\sigma \mathrm{eq}}={n}^{\sigma }{u}_{x}{u}_{y}$, ${\hat{f}}_{7}^{\sigma \mathrm{eq}}={n}^{\sigma }(T/{m}^{\sigma }+{u}_{y}^{2})$, ${\hat{f}}_{8}^{\sigma \mathrm{eq}}={n}^{\sigma }{\xi }^{\sigma }{u}_{x}$, ${\hat{f}}_{9}^{\sigma \mathrm{eq}}={n}^{\sigma }{\xi }^{\sigma }{u}_{y}$, ${\hat{f}}_{10}^{\sigma \mathrm{eq}}\,=\,3{n}^{\sigma }{u}_{x}T/{m}^{\sigma }+{n}^{\sigma }{u}_{x}^{3}$, ${\hat{f}}_{11}^{\sigma \mathrm{eq}}\,={n}^{\sigma }{u}_{y}T/{m}^{\sigma }+{n}^{\sigma }{u}_{x}^{2}{u}_{y}$, ${\hat{f}}_{12}^{\sigma \mathrm{eq}}\,=\,{n}^{\sigma }{u}_{x}T/{m}^{\sigma }+{n}^{\sigma }{u}_{x}{u}_{y}^{2}$, ${\hat{f}}_{13}^{\sigma \mathrm{eq}}\,=\,3{n}^{\sigma }{u}_{y}T/{m}^{\sigma }+{n}^{\sigma }{u}_{y}^{3}$, ${\hat{f}}_{14}^{\sigma \mathrm{eq}}$ $={n}^{\sigma }{\xi }^{\sigma }T/{m}^{\sigma }\,+{n}^{\sigma }{u}_{x}^{2}({\xi }^{\sigma }+2T/{m}^{\sigma })$, ${\hat{f}}_{15}^{\sigma \mathrm{eq}}\,={n}^{\sigma }{u}_{x}{u}_{y}({\xi }^{\sigma }+2T/{m}^{\sigma })$, and ${\hat{f}}_{16}^{\sigma \mathrm{eq}}$ $={n}^{\sigma }{\xi }^{\sigma }T/{m}^{\sigma }$ $+\,{n}^{\sigma }{u}_{y}^{2}({\xi }^{\sigma }+2T/{m}^{\sigma })$, with ξσ = (D + Iσ + 2)T/mσ + u2.The column matrix ${\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}$ is expressed byA2 ). The elements ${\hat{f}}_{i}^{\sigma \mathrm{seq}}$ are functions of (nσ, ${u}_{\alpha }^{\sigma }$, Tσ), and ${\hat{f}}_{i}^{\sigma \mathrm{eq}}$ are functions of (nσ, uα, T). The former are given by substituting (nσ, ${u}_{\alpha }^{\sigma }$, Tσ) for (nσ, uα, T) in the latter, respectively.
$\begin{eqnarray}{\hat{{\boldsymbol{f}}}}^{\sigma \mathrm{seq}}={\left(\begin{array}{l}{\hat{f}}_{1}^{\sigma \mathrm{seq}}\ {\hat{f}}_{2}^{\sigma \mathrm{seq}}\ \cdots \ {\hat{f}}_{N}^{\sigma \mathrm{seq}}\end{array}\right)}^{{\rm{T}}},\end{eqnarray}$
whose elements are similar to those in equation (Moreover, the column matrix ${\hat{{\boldsymbol{F}}}}^{\sigma }$ takes the form
$\begin{eqnarray}{\hat{{\boldsymbol{F}}}}^{\sigma }={\left(\begin{array}{l}{\hat{F}}_{1}^{\sigma }\ {\hat{F}}_{2}^{\sigma }\ \cdots \ {\hat{F}}_{N}^{\sigma }\end{array}\right)}^{{\rm{T}}},\end{eqnarray}$
where ${\hat{F}}_{1}^{\sigma }\,=\,0$, ${\hat{F}}_{2}^{\sigma }\,=\,{n}^{\sigma }{a}_{x}^{\sigma }$, ${\hat{F}}_{3}^{\sigma }\,=\,{n}^{\sigma }{a}_{y}^{\sigma }$, ${\hat{F}}_{4}^{\sigma }\,=\,2{n}^{\sigma }\left({u}_{x}^{\sigma }{a}_{x}^{\sigma }+{u}_{y}^{\sigma }{a}_{y}^{\sigma }\right)$, ${\hat{F}}_{5}^{\sigma }=2{n}^{\sigma }{u}_{x}^{\sigma }{a}_{x}^{\sigma }$, ${\hat{F}}_{6}^{\sigma }={n}^{\sigma }\left({u}_{x}^{\sigma }{a}_{y}^{\sigma }+{u}_{y}^{\sigma }{a}_{x}^{\sigma }\right)$, ${\hat{F}}_{7}^{\sigma }=2{n}^{\sigma }{u}_{y}^{\sigma }{a}_{y}^{\sigma }$, ${\hat{F}}_{8}^{\sigma }=2{n}^{\sigma }{u}_{x}^{\sigma }\left({u}_{x}^{\sigma }{a}_{x}^{\sigma }+{u}_{y}^{\sigma }{a}_{y}^{\sigma }\right)+{n}^{\sigma }{a}_{x}^{\sigma }\left[{{u}^{\sigma }}^{2}+\left(D+{I}^{\sigma }+2\right){T}^{\sigma }/{m}^{\sigma }\ \right]$, ${\hat{F}}_{9}^{\sigma }=2{n}^{\sigma }{u}_{y}^{\sigma }\left({u}_{x}^{\sigma }{a}_{x}^{\sigma }+{u}_{y}^{\sigma }{a}_{y}^{\sigma }\right)+{n}^{\sigma }{a}_{y}^{\sigma }[{{u}^{\sigma }}^{2}+\left(D+{I}^{\sigma }+2\right){T}^{\sigma }/{m}^{\sigma }]$, ${\hat{F}}_{10}^{\sigma }\, =\,3{n}^{\sigma }{a}_{x}^{\sigma }\left({{u}_{x}^{\sigma }}^{2}\,+\,{T}^{\sigma }/{m}^{\sigma }\ \right)$, ${\hat{F}}_{11}^{\sigma }\,=\,2{n}^{\sigma }{a}_{x}^{\sigma }{u}_{x}^{\sigma }{u}_{y}^{\sigma }\, +\,{n}^{\sigma }{a}_{y}^{\sigma }({{u}_{x}^{\sigma }}^{2}\,+\,{T}^{\sigma }/{m}^{\sigma })$, ${\hat{F}}_{12}^{\sigma }={n}^{\sigma }{a}_{x}^{\sigma }\left({{u}_{y}^{\sigma }}^{2} +{T}^{\sigma }/{m}^{\sigma }\right)+2{n}^{\sigma }{a}_{y}^{\sigma }{u}_{x}^{\sigma }{u}_{y}^{\sigma }$, ${\hat{F}}_{13}^{\sigma }$ $=3{n}^{\sigma }{a}_{y}^{\sigma }\times \left({{u}_{y}^{\sigma }}^{2}+{T}^{\sigma }/{m}^{\sigma }\ \right)$, ${\hat{F}}_{14}^{\sigma } $ $=2{n}^{\sigma }{a}_{x}^{\sigma }{u}_{x}^{\sigma } $ $ \left[2{{u}_{x}^{\sigma }}^{2} +{{u}_{y}^{\sigma }}^{2}+\left(D+{I}^{\sigma } +5\right){T}^{\sigma }/{m}^{\sigma }\ \right] $ $+\,2{n}^{\sigma }{a}_{y}^{\sigma }{u}_{y}^{\sigma }\left({{u}_{x}^{\sigma }}^{2}+{T}^{\sigma }/{m}^{\sigma }\right)$, ${\hat{F}}_{15}={n}^{\sigma }{a}_{x}^{\sigma }{u}_{y}^{\sigma }[3{{u}_{x}^{\sigma }}^{2}$ $+{{u}_{y}^{\sigma }}^{2}+\left(D+{I}^{\sigma } +4\right){T}^{\sigma }/{m}^{\sigma }]\,+\,{n}^{\sigma }{a}_{y}^{\sigma }{u}_{x}^{\sigma }[{{u}_{x}^{\sigma }}^{2}+3{{u}_{y}^{\sigma }}^{2}$ $+(D+{I}^{\sigma }+4){T}^{\sigma }/{m}^{\sigma }]$, and ${\hat{F}}_{16}^{\sigma }=2{n}^{\sigma }{a}_{x}^{\sigma }{u}_{x}^{\sigma } $ $ \left({{u}_{y}^{\sigma }}^{2} +{T}^{\sigma }/{m}^{\sigma }\right) $ $+2{n}^{\sigma }{a}_{y}^{\sigma }{u}_{y}^{\sigma }\times [{{u}_{x}^{\sigma }}^{2} $ $+2{{u}_{y}^{\sigma }}^{2}+(D+{I}^{\sigma }+5){T}^{\sigma }/{m}^{\sigma }]$.Additionally, the column matrix ${\hat{{\boldsymbol{R}}}}^{\sigma }$ is calculated by
$\begin{eqnarray}{\hat{{\boldsymbol{R}}}}^{\sigma }={\left(\begin{array}{l}{\hat{R}}_{1}^{\sigma }\ {\hat{R}}_{2}^{\sigma }\ \cdots \ {\hat{R}}_{N}^{\sigma }\end{array}\right)}^{{\rm{T}}},\end{eqnarray}$
with ${\hat{R}}_{1}^{\sigma }={{n}^{\sigma }}^{{\prime} }$, ${\hat{R}}_{2}^{\sigma }={{n}^{\sigma }}^{{\prime} }{u}_{x}$, ${\hat{R}}_{3}^{\sigma }={{n}^{\sigma }}^{{\prime} }{u}_{y}$, ${\hat{R}}_{4}^{\sigma }={{n}^{\sigma }}^{{\prime} }[(D+{I}^{\sigma })T/{m}^{\sigma }+{u}^{2}]+(D+{I}^{\sigma }){n}^{\sigma }T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{5}^{\sigma }={{n}^{\sigma }}^{{\prime} }\times (T/{m}^{\sigma }+{u}_{x}^{2})+{n}^{\sigma }T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{6}^{\sigma }={{n}^{\sigma }}^{{\prime} }{u}_{x}{u}_{y}$, ${\hat{R}}_{7}^{\sigma }\,={{n}^{\sigma }}^{{\prime} }\times (T/{m}^{\sigma }$ $+{u}_{y}^{2})+{n}^{\sigma }T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{8}^{\sigma }\,=\,{{n}^{\sigma }}^{{\prime} }{u}_{x}{\xi }^{\sigma }\,$ $+\,(D+{I}^{\sigma }+2)\times {n}^{\sigma }{u}_{x}T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{9}^{\sigma }={{n}^{\sigma }}^{{\prime} }{u}_{y}{\xi }^{\sigma }$ $+(D+{I}^{\sigma }+2){n}^{\sigma }{u}_{y}T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{10}^{\sigma }=3{{n}^{\sigma }}^{{\prime} }{u}_{x}T/{m}^{\sigma }$ $+{{n}^{\sigma }}^{{\prime} }{u}_{x}^{3}+3{n}^{\sigma }{u}_{x}T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{11}^{\sigma }={{n}^{\sigma }}^{{\prime} }{u}_{y}T/{m}^{\sigma }$ $+{{n}^{\sigma }}^{{\prime} }{u}_{x}^{2}{u}_{y}+{n}^{\sigma }{u}_{y}T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{12}^{\sigma }={{n}^{\sigma }}^{{\prime} }{u}_{x}T/{m}^{\sigma }$ $+{{n}^{\sigma }}^{{\prime} }{u}_{x}{u}_{y}^{2}+\,{n}^{\sigma }{u}_{x}T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{13}^{\sigma }=3{{n}^{\sigma }}^{{\prime} }{u}_{y}T/{m}^{\sigma }$ $+{{n}^{\sigma }}^{{\prime} }{u}_{y}^{3}+3{n}^{\sigma }{u}_{y}T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{14}^{\sigma }={{n}^{\sigma }}^{{\prime} }{\xi }^{\sigma }T/{m}^{\sigma }+{{n}^{\sigma }}^{{\prime} }{u}_{x}^{2}({\xi }^{\sigma }+2T/{m}^{\sigma })+{n}^{\sigma }\times [2(D\,$ $+\,{I}^{\sigma }\,+\,2)T/{m}^{\sigma }+{u}^{2}+(D+{I}^{\sigma }+4){u}_{x}^{2}]T^{\prime} /{m}^{\sigma }$, ${\hat{R}}_{15}^{\sigma }={{n}^{\sigma }}^{{\prime} }{u}_{x}{u}_{y}({\xi }^{\sigma }+2T/{m}^{\sigma })$ $+{n}^{\sigma }[(D+{I}^{\sigma }+4){u}_{x}{u}_{y}]T^{\prime} /{m}^{\sigma }$, and ${\hat{R}}_{16}^{\sigma }={{n}^{\sigma }}^{{\prime} }{\xi }^{\sigma }T/{m}^{\sigma }$ $+{{n}^{\sigma }}^{{\prime} }{u}_{y}^{2}({\xi }^{\sigma }+2T/{m}^{\sigma })+{n}^{\sigma }\times [2(D+{I}^{\sigma }+2)T/{m}^{\sigma }+{u}^{2}+(D+{I}^{\sigma }+4){u}_{y}^{2}]T^{\prime} /{m}^{\sigma }$.Appendix B. Hydrodynamic equations
Via the CE expansion, it can be found that the current DBM is consistent with the reactive NS equations in the hydrodynamic limit. Based on the Einstein summation convention, the NS equations of individual species read,
$\begin{eqnarray}\displaystyle \frac{\partial {\rho }^{\sigma }}{\partial t}+\displaystyle \frac{\partial }{\partial \alpha }\left({\rho }^{\sigma }{u}_{\alpha }^{\sigma }\right)={{\rho }^{\sigma }}^{{\prime} },\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left({\rho }^{\sigma }{u}_{\alpha }^{\sigma }\right)+\displaystyle \frac{\partial }{\partial \beta }\left({\delta }_{\alpha \beta }{p}^{\sigma }+{\rho }^{\sigma }{u}_{\alpha }^{\sigma }{u}_{\beta }^{\sigma }+{P}_{\alpha \beta }^{\sigma }+{U}_{\alpha \beta }^{\sigma }+{V}_{\alpha \beta }^{\sigma }\right)\\ =\,{S}_{J\alpha }^{\sigma }{\rho }^{\sigma }\left({u}_{\alpha }-{u}_{\alpha }^{\sigma }\right)+{\rho }^{\sigma }{a}_{\alpha }^{\sigma }+{\rho }^{\sigma ^{\prime} }{u}_{\alpha },\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial {E}^{\sigma }}{\partial t}+\displaystyle \frac{\partial }{\partial \alpha }\left({E}^{\sigma }{u}_{\alpha }^{\sigma }+{p}^{\sigma }{u}_{\alpha }^{\sigma }-{\kappa }_{\alpha }^{\sigma }\displaystyle \frac{\partial {T}^{\sigma }}{\partial \alpha }+{u}_{\beta }^{\sigma }{P}_{\alpha \beta }^{\sigma }+{X}_{\alpha }^{\sigma }+{Y}_{\alpha }^{\sigma }\right)\\ =\displaystyle \frac{1}{2}{S}_{4}^{2\sigma }{\rho }^{\sigma }\left[\left(D+{I}^{\sigma }\right)\displaystyle \frac{T-{T}^{\sigma }}{{m}^{\sigma }}+{u}^{2}-{{u}^{\sigma }}^{2}\right]+{\rho }^{\sigma }{u}_{\alpha }^{\sigma }{a}_{\alpha }^{\sigma }+{{E}^{\sigma }}^{{\prime} },\end{array}\end{eqnarray}$
in terms of $\begin{eqnarray}{P}_{\alpha \beta }^{\sigma }={\mu }_{\alpha \beta }^{\sigma }\left(\displaystyle \frac{2{\delta }_{\alpha \beta }}{D+{I}^{\sigma }}\displaystyle \frac{\partial {u}_{\chi }^{\sigma }}{\partial \chi }-\displaystyle \frac{\partial {u}_{\alpha }^{\sigma }}{\partial \beta }-\displaystyle \frac{\partial {u}_{\beta }^{\sigma }}{\partial \alpha }\right),\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{U}_{\alpha \beta }^{\sigma }=\displaystyle \frac{{S}_{4}^{2\sigma }}{{S}_{\alpha \beta }^{1\sigma }}\displaystyle \frac{{\delta }_{\alpha \beta }{\rho }^{\sigma }}{D+{I}^{\sigma }}\left({{u}^{\sigma }}^{2}-{u}^{2}\right)-\displaystyle \frac{{S}_{J\chi }^{\sigma }}{{S}_{\alpha \beta }^{1\sigma }}\displaystyle \frac{2{\delta }_{\alpha \beta }{\rho }^{\sigma }}{D+{I}^{\sigma }}{u}_{\chi }^{\sigma }\left({u}_{\chi }^{\sigma }-{u}_{\chi }\right)\\ -\displaystyle \frac{{S}_{\alpha \beta }^{2\sigma }}{{S}_{\alpha \beta }^{1\sigma }}{\rho }^{\sigma }\left({u}_{\alpha }^{\sigma }{u}_{\beta }^{\sigma }-{u}_{\alpha }{u}_{\beta }\right)+\displaystyle \frac{{S}_{J\alpha }^{\sigma }}{{S}_{\alpha \beta }^{1\sigma }}{\rho }^{\sigma }\left({u}_{\alpha }^{\sigma }-{u}_{\alpha }\right){u}_{\beta }^{\sigma }\\ +\displaystyle \frac{{S}_{J\beta }^{\sigma }}{{S}_{\alpha \beta }^{1\sigma }}{\rho }^{\sigma }{u}_{\alpha }^{\sigma }\left({u}_{\beta }^{\sigma }-{u}_{\beta }\right)+{\delta }_{\alpha \beta }\displaystyle \frac{{S}_{4}^{2\sigma }-{S}_{\alpha \beta }^{2\sigma }}{{S}_{\alpha \beta }^{1\sigma }}{\rho }^{\sigma }\displaystyle \frac{{T}^{\sigma }-T}{{m}^{\sigma }},\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{V}_{\alpha \beta }^{\sigma }=\displaystyle \frac{{\rho }^{\sigma ^{\prime} }}{{S}_{\alpha \beta }^{1\sigma }}\left(\Space{0ex}{3.55ex}{0ex}{u}_{\alpha }{u}_{\beta }+{u}_{\alpha }^{\sigma }{u}_{\beta }^{\sigma }-{u}_{\alpha }{u}_{\beta }^{\sigma }-{u}_{\alpha }^{\sigma }{u}_{\beta }\right.\\ \ \ -\,\left.{\delta }_{\alpha \beta }\displaystyle \frac{{{u}^{\sigma }}^{2}+{u}^{2}-2{u}_{\chi }^{\sigma }{u}_{\chi }}{D+{I}^{\sigma }}\right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{X}_{\alpha }^{\sigma }=\displaystyle \frac{{S}_{4}^{2\sigma }}{{S}_{\kappa \alpha }^{1\sigma }}\displaystyle \frac{{\rho }^{\sigma }{u}_{\alpha }^{\sigma }}{D+{I}^{\sigma }}\left({{u}^{\sigma }}^{2}-{u}^{2}\right)\\ \ \ -\displaystyle \frac{2{S}_{J\beta }^{\sigma }}{{S}_{\kappa \alpha }^{1\sigma }}\displaystyle \frac{{\rho }^{\sigma }{u}_{\alpha }^{\sigma }}{D+{I}^{\sigma }}{u}_{\beta }^{\sigma }\left({u}_{\beta }^{\sigma }-{u}_{\beta }\right)\\ +\displaystyle \frac{{S}_{J\alpha }^{\sigma }-{S}_{J\beta }^{\sigma }}{{S}_{\kappa \alpha }^{1\sigma }}{\rho }^{\sigma }{u}_{\alpha }^{\sigma }{u}_{\beta }^{\sigma }\left({u}_{\beta }^{\sigma }-{u}_{\beta }\right)\\ +\displaystyle \frac{{S}_{4}^{2\sigma }{u}_{\alpha }^{\sigma }-{S}_{\kappa \alpha }^{2\sigma }{u}_{\alpha }}{2{S}_{\kappa \alpha }^{1\sigma }}{\rho }^{\sigma }\left[\left(D+{I}^{\sigma }+2\right)\displaystyle \frac{{T}^{\sigma }-T}{{m}^{\sigma }}+{{u}^{\sigma }}^{2}-{u}^{2}\right]\\ +\displaystyle \frac{{S}_{J\alpha }^{\sigma }-{S}_{\kappa \alpha }^{2\sigma }}{2{S}_{\kappa \alpha }^{1\sigma }}{\rho }^{\sigma }\left({u}_{\alpha }^{\sigma }-{u}_{\alpha }\right)\left[\left(D+{I}^{\sigma }+2\right)\displaystyle \frac{{T}^{\sigma }}{{m}^{\sigma }}+{{u}^{\sigma }}^{2}\right],\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}{Y}_{\alpha }^{\sigma }=\displaystyle \frac{{{\rho }^{\sigma }}^{{\prime} }}{2{S}_{\kappa \alpha }^{1\sigma }}\left({u}_{\alpha }^{\sigma }-{u}_{\alpha }\right)\left[\left(D+{I}^{\sigma }+2\right)\displaystyle \frac{{T}^{\sigma }-T}{{m}^{\sigma }}+{{u}^{\sigma }}^{2}-{u}^{2}\right]\\ -\displaystyle \frac{1}{{S}_{\kappa \alpha }^{1\sigma }}\displaystyle \frac{{\rho }^{\sigma ^{\prime} }{u}_{\alpha }^{\sigma }}{D+{I}^{\sigma }}{\left({u}_{\beta }^{\sigma }-{u}_{\beta }\right)}^{2}+\left(D+{I}^{\sigma }+2\right)\displaystyle \frac{{\rho }^{\sigma }T^{\prime} }{2{S}_{\kappa \alpha }^{1\sigma }{m}^{\sigma }}\left({u}_{\alpha }-{u}_{\alpha }^{\sigma }\right),\end{array}\end{eqnarray}$
with the change rate of individual energy due to the chemical reaction $\begin{eqnarray}{{E}^{\sigma }}^{{\prime} }={{\rho }^{\sigma }}^{{\prime} }\left(\displaystyle \frac{D+{I}^{\sigma }}{2}\displaystyle \frac{T}{{m}^{\sigma }}+\displaystyle \frac{{u}^{2}}{2}\right)+\displaystyle \frac{D+{I}^{\sigma }}{2}\displaystyle \frac{{\rho }^{\sigma }}{{m}^{\sigma }}T^{\prime} =\displaystyle \frac{{m}^{\sigma }{\hat{R}}_{4}}{2},\end{eqnarray}$
the thermal conductivity $\begin{eqnarray}{\kappa }_{\alpha }^{\sigma }=\displaystyle \frac{D+{I}^{\sigma }+2}{2{S}_{\kappa \alpha }^{1\sigma }}\displaystyle \frac{{\rho }^{\sigma }{T}^{\sigma }}{{{m}^{\sigma }}^{2}},\end{eqnarray}$
the dynamic viscosity $\begin{eqnarray}{\mu }_{\alpha \beta }^{\sigma }=\displaystyle \frac{{p}^{\sigma }}{{S}_{\alpha \beta }^{1\sigma }},\end{eqnarray}$
and the parameters (${S}_{{Jx}}^{\sigma }$, ${S}_{{Jy}}^{\sigma }$, ${S}_{{xx}}^{1\sigma }$, ${S}_{{xy}}^{1\sigma }$, ${S}_{{yy}}^{1\sigma }$, ${S}_{{xx}}^{2\sigma }$, ${S}_{{xy}}^{2\sigma }$, ${S}_{{yy}}^{2\sigma }$, ${S}_{\kappa x}^{1\sigma }$, ${S}_{\kappa y}^{1\sigma }$ ${S}_{\kappa x}^{2\sigma }$, ${S}_{\kappa y}^{2\sigma }$) = (${S}_{2}^{2\sigma }$, ${S}_{3}^{2\sigma }$, ${S}_{5}^{1\sigma }$, ${S}_{6}^{1\sigma }$, ${S}_{7}^{1\sigma }$, ${S}_{5}^{2\sigma }$, ${S}_{6}^{2\sigma }$, ${S}_{7}^{2\sigma }$, ${S}_{8}^{1\sigma }$, ${S}_{9}^{1\sigma }$, ${S}_{8}^{2\sigma }$, ${S}_{9}^{2\sigma }$).In the case of ${S}_{5}^{1\sigma }={S}_{6}^{1\sigma }={S}_{7}^{1\sigma }={S}_{\mu }^{\sigma }$ and ${S}_{8}^{1\sigma }={S}_{9}^{1\sigma }\,={S}_{\kappa }^{\sigma }$, the Prandtl number isB1 )–(B3 ) over all species σ, we get the following reactive NS equations,
$\begin{eqnarray}{\Pr }^{\sigma }=\displaystyle \frac{{S}_{\kappa }^{\sigma }}{{S}_{\mu }^{\sigma }}.\end{eqnarray}$
Moreover, the specific heat at constant pressure and volume are ${c}_{{\rm{p}}}^{\sigma }=(D+{I}^{\sigma }+2)/(2{m}^{\sigma })$ and ${c}_{{\rm{v}}}^{\sigma }=\left(D+{I}^{\sigma }\right)/\left(2{m}^{\sigma }\right)$, respectively. Consequently, the specific heat ratio is $\begin{eqnarray}{\gamma }^{\sigma }=\displaystyle \frac{{c}_{{\rm{p}}}^{\sigma }}{{c}_{{\rm{v}}}^{\sigma }}=\displaystyle \frac{D+{I}^{\sigma }+2}{D+{I}^{\sigma }}.\end{eqnarray}$
Additionally, summing equations ( $\begin{eqnarray}\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{\partial }{\partial \alpha }\left(\rho {u}_{\alpha }\right)=0,\end{eqnarray}$
$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}\left(\rho {u}_{\alpha }\right)+\displaystyle \frac{\partial }{\partial \beta }{\sum }_{\sigma }\left({\delta }_{\alpha \beta }{p}^{\sigma }+{\rho }^{\sigma }{u}_{\alpha }^{\sigma }{u}_{\beta }^{\sigma }+{P}_{\alpha \beta }^{\sigma }+{U}_{\alpha \beta }^{\sigma }+{V}_{\alpha \beta }^{\sigma }\right)=\rho {a}_{\alpha },\end{eqnarray}$
$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{\partial E}{\partial t}+\displaystyle \frac{\partial }{\partial \alpha }{\sum }_{\sigma }\left({E}^{\sigma }{u}_{\alpha }^{\sigma }+{p}^{\sigma }{u}_{\alpha }^{\sigma }-{\kappa }_{\alpha }^{\sigma }\displaystyle \frac{\partial {T}^{\sigma }}{\partial \alpha }\right.\\ \,+\left.\,{u}_{\beta }^{\sigma }{P}_{\alpha \beta }^{\sigma }+{X}_{\alpha }^{\sigma }+{Y}_{\alpha }^{\sigma }\Space{0ex}{2.8ex}{0ex}\right)=\rho {u}_{\alpha }{a}_{\alpha }+E^{\prime} ,\end{array}\end{eqnarray}$
in the case of ${a}_{\alpha }^{\sigma }={a}_{\alpha }$, where $E^{\prime} ={\sum }_{\sigma }{{E}^{\sigma }}^{{\prime} }$ denotes the change rate of mixing energy due to the chemical reaction.Actually, under some corresponding conditions, we could obtain the Fick's laws of diffusion and Maxwell-Stefan diffusion equation from equations (B1 ) and (B2 ) as well [9]. More discussion is beyond this work.
Appendix C. Nondimensional parameters
Now, let us demonstrate some important parameters in the current DBM.
(i) Time. The dimensional time is td = Ld/ud, the nondimensional time tn = Ln/un, and the time ratio
$\begin{eqnarray}{t}_{{\rm{r}}}=\displaystyle \frac{{L}_{{\rm{d}}}{u}_{{\rm{n}}}}{{L}_{{\rm{n}}}{u}_{{\rm{d}}}}.\end{eqnarray}$
(ii) Energy. The internal energies with dimension and nondimension are Eintd = (D + I)ndRTd/2 and Eintn = (D + I)nnTn/2, respectively, and the energy ratio reads
$\begin{eqnarray}{E}_{{\rm{r}}}=\displaystyle \frac{{E}_{\mathrm{intd}}}{{E}_{\mathrm{intn}}}={n}_{{\rm{r}}}{{RT}}_{{\rm{r}}}.\end{eqnarray}$
(iii) Mass. Given the dimensional and nondimensional molar mass md and mn, the mass densities are ρd = mdnd and ρn = mnnn in dimensional and nondimensional forms, the kinetic energies are ${E}_{\mathrm{kd}}={\rho }_{{\rm{d}}}{u}_{{\rm{d}}}^{2}/2$ and ${E}_{\mathrm{kn}}={\rho }_{{\rm{n}}}{u}_{{\rm{n}}}^{2}/2$ in the two forms, so the energy ratio isC2 ) and (C3 ), we obtain the mass ratio
$\begin{eqnarray}{E}_{{\rm{r}}}=\displaystyle \frac{{E}_{\mathrm{kd}}}{{E}_{\mathrm{kn}}}={m}_{{\rm{r}}}{n}_{{\rm{r}}}{u}_{{\rm{r}}}^{2}.\end{eqnarray}$
From equations ( $\begin{eqnarray}{m}_{{\rm{r}}}=\displaystyle \frac{{{RT}}_{{\rm{r}}}}{{u}_{{\rm{r}}}^{2}},\end{eqnarray}$
which leads to the mass density ratio ρr = mrnr.(iv) Viscosity. The Reynolds number is defined as
$\begin{eqnarray}\mathrm{Re}=\displaystyle \frac{{\rho }_{{\rm{d}}}{u}_{{\rm{d}}}{L}_{{\rm{d}}}}{{\mu }_{{\rm{d}}}}=\displaystyle \frac{{\rho }_{{\rm{n}}}{u}_{{\rm{n}}}{L}_{{\rm{n}}}}{{\mu }_{{\rm{n}}}},\end{eqnarray}$
where μd and μn denote the dynamic viscosity in dimensional and nondimensional forms, respectively. The ratio of dimensional to nondimensional dynamic viscosity takes the form $\begin{eqnarray}{\mu }_{{\rm{r}}}=\displaystyle \frac{{n}_{{\rm{r}}}{L}_{{\rm{r}}}{{RT}}_{{\rm{r}}}}{{u}_{{\rm{r}}}},\end{eqnarray}$
from which we could get the kinematic viscosity ratio νr = μr/ρr.(v) Thermal diffusivity. The Prandtl number is a dimensionless number defined as the ratio of momentum diffusivity (i.e., kinematic viscosity) to thermal diffusivity,
$\begin{eqnarray}\Pr =\displaystyle \frac{{\nu }_{{\rm{d}}}}{{K}_{{\rm{d}}}}=\displaystyle \frac{{\nu }_{{\rm{n}}}}{{K}_{{\rm{n}}}},\end{eqnarray}$
where (νd, νn) and (Kd, Kn) denote the momentum diffusivity and thermal diffusivity in dimensional and nondimensional forms, respectively. Consequently, the ratio of dimensional to nondimensional thermal diffusivity is Kr = νr.(vi) Mass diffusivity. The Schmidt number is a dimensionless number defined as the ratio of momentum diffusivity (kinematic viscosity) to mass diffusivity,
$\begin{eqnarray}\mathrm{Sc}=\displaystyle \frac{{\nu }_{{\rm{d}}}}{{\zeta }_{{\rm{d}}}}=\displaystyle \frac{{\nu }_{{\rm{n}}}}{{\zeta }_{{\rm{n}}}},\end{eqnarray}$
where ζd and ζn stand for the dimensional and nondimensional mass diffusivity, respectively. Hence, the ratio of dimensional to nondimensional mass diffusivity is ζr = νr.Appendix D. Homogeneous mixture with body force
When a homogeneous nonreactive mixture reaches its steady state in the force field, the time and space derivatives equal zero, hence equation (B2 ) leads to
$\begin{eqnarray}{S}_{J\alpha }^{\sigma }{\rho }^{\sigma }\left({u}_{\alpha }-{u}_{\alpha }^{\sigma }\right)+{\rho }^{\sigma }{a}_{\alpha }^{\sigma }=0.\end{eqnarray}$
Namely, the individual velocities read $\begin{eqnarray}{u}_{\alpha }^{\sigma }={u}_{\alpha }+{\left({S}_{J\alpha }^{\sigma }\right)}^{-1}{a}_{\alpha }^{\sigma }.\end{eqnarray}$
Moreover, if the homogeneous nonreactive mixture is an isothermal system in the steady state, equation (B3 ) gives
$\begin{eqnarray}\displaystyle \frac{1}{2}{S}_{4}^{2\sigma }{\rho }^{\sigma }\left[\left(D+{I}^{\sigma }\right)\displaystyle \frac{T-{T}^{\sigma }}{{m}^{\sigma }}+{u}^{2}-{{u}^{\sigma }}^{2}\right]+{\rho }^{\sigma }{u}_{\alpha }^{\sigma }{a}_{\alpha }^{\sigma }=0.\end{eqnarray}$
Consequently, the individual temperatures take the form $\begin{eqnarray}{T}^{\sigma }=T+\displaystyle \frac{{m}^{\sigma }}{D+{I}^{\sigma }}\left[{\left|{\boldsymbol{u}}\right|}^{2}-{\left|{{\boldsymbol{u}}}^{\sigma }\right|}^{2}+2{\left({S}_{4}^{2\sigma }\right)}^{-1}{{\boldsymbol{u}}}^{\sigma }\cdot {{\boldsymbol{a}}}^{\sigma }\right].\end{eqnarray}$
On the contrary, for a thermal homogeneous mixture in the force field, if the individual accelerations are not equal, the chemical species propagate collectively with different velocities and collide randomly with each other, which leads to the change of energies and temperatures. That is to say, the work done by external forces transforms into the kinetic and internal energies of the thermal system, i.e,B16 ). Let us consider a special case where the mixing velocity keeps constant u = 0 in the force field, then equation (D5 ) changes into13 ) and (D6 ), we get
$\begin{eqnarray}\displaystyle \frac{\partial E}{\partial t}=\displaystyle \sum _{\sigma }{\rho }^{\sigma }{{\boldsymbol{a}}}^{\sigma }\cdot {{\boldsymbol{u}}}^{\sigma },\end{eqnarray}$
which can be derived from equation ( $\begin{eqnarray}\displaystyle \frac{\partial \displaystyle {E}_{\mathrm{int}}}{\partial t}=\displaystyle \sum _{\sigma }{\left({S}_{J}^{\sigma }\right)}^{-1}{\rho }^{\sigma }{{\boldsymbol{a}}}^{\sigma }\cdot {{\boldsymbol{a}}}^{\sigma }.\end{eqnarray}$
From equations ( $\begin{eqnarray}{T}^{\sigma }\approx T=\displaystyle \frac{{\sum }_{\sigma }(D+{I}^{\sigma }){n}^{\sigma }{T}_{0}+2{\sum }_{\sigma }{\left({S}_{J}^{\sigma }\right)}^{-1}{\rho }^{\sigma }{\left|{{\boldsymbol{a}}}^{\sigma }\right|}^{2}t}{{\sum }_{\sigma }\left(D+{I}^{\sigma }\right){n}^{\sigma }},\end{eqnarray}$
where T0 is the initial mixing temperature.