A multi-component discrete Boltzmann model for nonequilibrium reactive flows

We propose a multi-component discrete Boltzmann model (DBM) for premixed, nonpremixed, or partially premixed nonequilibrium reactive flows. This model is suitable for both subsonic and supersonic flows with or without chemical reaction and/or external force. A two-dimensional sixteen-velocity model is constructed for the DBM. In the hydrodynamic limit, the DBM recovers the modified Navier-Stokes equations for reacting species in a force field. Compared to standard lattice Boltzmann models, the DBM presents not only more accurate hydrodynamic quantities, but also detailed nonequilibrium effects that are essential yet long-neglected by traditional fluid dynamics. Apart from nonequilibrium terms (viscous stress and heat flux) in conventional models, specific hydrodynamic and thermodynamic nonequilibrium quantities (high order kinetic moments and their departure from equilibrium) are dynamically obtained from the DBM in a straightforward way. Due to its generality, the developed methodology is applicable to a wide range of phenomena across many energy technologies, emissions reduction, environmental protection, mining accident prevention, chemical and process industry.

Over the past years, the versatile DBM has been effectively applied to thermal phase separation, fluid instabilities, reactive flows, etc [33][34][35][36][37][38] . The DBM for reactive flows was firstly presented by Yan et al. in 2013 35 . Then, Lin et al. extended the DBM to reactive flows in a polar coordinate 36 . In 2015, Xu et al. proposed a multiple-relaxation-time DBM for reactive flows where the specific heat ratio and Prandtl number are adjustable 37 . The next year, a DBM is formulated for reactive flows where chemical reactant and product are described by two coupled distribution functions 38 . However, previous DBMs are suitable for premixed reactive flows, but not for nonpremixed or partially premixed reactive flows [35][36][37][38] . For the sake of simulating both subsonic and supersonic nonequilibrium reactive flows with premixed, nonpremixed, or partially premixed reactants, we propose a multi-component DBM in this work. The DBM presents two ways to access the thermodynamic nonequilibrium behaviours. One is to measure the viscous stress and heat flux that are described by traditional NS models; The other is to calculate the kinetic moments of the difference between equilibrium and nonequilibrium discrete distribution functions, which is beyond conventional hydrodynamic models. Such capability is the main object of the present work.

Discrete Boltzmann model
Without loss of generality, we consider the oxidation of propane in air using the one-step overall reaction, = [−1, −5, 3,4]. Nitrogen is assumed to be inert. The overall reaction rate reads a ov ov C H O 3 8 2 with k ov the reaction coefficient, n σ molar concentration, E a effective activation energy, R universal gas constant, T temperature. The mass change rate of species σ is ω σ = a σ ⋅ M σ ⋅ ω ov . In addition to the one-step reaction, detailed or reduced multi-step chemical kinetics can also be employed. Let us introduce the discrete Boltzmann equation, Here σ f i ( σ f i eq ) indicates the discrete (equilibrium) distribution function, σ v i the discrete velocity, t (τ σ ) the (relaxation) time. Ω σ i , σ R i , and σ G i are the collision, reaction and force terms accounting for the molecular collision, chemical reaction and external force, respectively. The collision term in Eq. (4) obeys the conservation of mass, momentum, and energy, from which the relations between the physical quantities ( σ n , n, σ u , u, σ T , T) and the distribution function σ f i are obtained 38 . The symbols with (without) superscript σ denote the physical quantities of the species (mixture). In Eq. (5), σ⁎ n and ⁎ T ( σ n and T) denote the molar concentration and temperature after (before) chemical reaction within time step τ σ . Similarly, in Eq. (6), the hydrodynamic velocity changes from σ u to σ † u within time τ σ due to external force, meanwhile the temperature changes from σ T to σ † T . The discrete equilibrium function σ f i eq is linked with the formula, eq D/2 1/2 2 2 in the way that a required number of kinetic moments calculated by the integral of f σeq are equal to those by the summation of f i σeq . In Eq. (7), m σ stands for molar mass, D = 2 space dimension, I σ extra degrees of freedom corresponding to molecular rotation or vibration. There are 7 moments satisfied by f i σeq = f i σeq (n σ ,u,T) in this work. Specifically, i i 2 . The discrete velocities v i σ and η i σ are (see Fig. 1), And the discrete equilibrium distribution is expressed by = .
−f M f (16) eq eq 1 One significant capability of the DBM is to investigate nonequilibrium manifestations by measuring the following physical variables, is a second order tensor with four components, among which only three 3 3 is a third order tensor with eight components where only is the first order tensor (i.e. vector) contracted from a third order tensor and have two independent components; Similarly for Δ = Δ . It is easy to prove that, via the Chapman-Enskog multiscale analysis, the DBM is in line with the following modified NS equations, 2 )/( ) specific heat ratio. The superscript "′" represents the change rate of physical quantities due to the chemical reaction. In fact, applying the operator ∑ σ to both sides of Eqs (21)- (23) gives NS equations for the whole system, which reduces to conventional NS equations when = σ u u and = σ T T . Obviously, Eqs (21)-(23) gives a more detailed description than the conventional NS equations. The latter is just a special case of the former.
Furthermore, dynamic viscosity and heat conductivity in the NS equations are regarded as two important thermodynamic nonequilibrium manifestations or physical effects on fluid flows. In fact, a more detailed way to study the nonequilibrium behaviours is to investigate the departure of high order velocity moments from their local equilibrium counterparts, as shown in Eqs (17)- (20). Those kinetic moments of the difference between nonequilibrium and equilibrium distribution functions have significant physical meanings. In particular, Δ σ 2 is associated with viscous stress tensor and nonorganised momentum fluxes, Δ σ 3,1 and Δ σ 3 are related to nonorganised energy (heat) fluxes, Δ σ 4,2 corresponds to the flux of nonorganised energy (heat) flux 36,37 , The terminology "nonorganised" is relative to "organised". The latter refers to the collective motion of a fluid flow, while the former corresponds to the molecular individualism on top of the collective motion 34

Numerical simulation
To validate this DBM, we conduct five simulation tests. Test one is the combustion of (premixed, nonpremixed, and partially premixed) propane-air filled in a free-falling box. The released heat in constant volume and the external force effects are demonstrated. The second test is a subsonic flame at constant pressure. In the third part, to show its suitability for high speed compressible systems, the DBM is used to simulate a shock wave. Its capability to investigate nonequilibrium effects is verified as well. A supersonic reacting wave is simulated in the fourth part. The first four tests are 1-dimensional (1-D) cases. The last is for a typical 2-D case, Kelvin-Helmholtz (KH) instability.
Moreover, the first order Euler forward time discretization and the second order nonoscillatory and nonfree-parameters dissipative finite difference scheme 39 are adopted for the temporal and spatial derivatives in Eq. (3) in this section. Hence, the discrete velocities v i are independent of the grid mesh Δx and Δy. For the purpose of accuracy and robustness, it is preferable to set the values of discrete velocities Combustion in constant volume. First of all, we simulate the combustion of propane-air filled in a free-falling box, which consists of three parts with volumes V 1 , V 2 and V 3 , respectively. The fixed volume of the box is V 0 = V 1 + V 2 + V 3 , and V 1 :V 2 :V 3 = 3:119:78. Initially, the left part is filled with propane, the middle part is full of air, and the right part is occupied by the propane-air mixture with equivalence ratio 0.6. In each part, the particle number density is 40.6mol ⋅ m −3 , temperature 300 K, and pressure 1 atm. Premixed, nonpremixed and partially premixed combustion phenomena take place simultaneously in this box after ignition. Specifically, the nonpremixed combustion takes place between the left and middle parts, the partially premixed combustion occurs between the middle and right parts with a changing equivalence ratio, and the premixed combustion is in the rightmost part with a constant equivalence ratio. Three discrete velocity models (D2V16, D2V24 37 , and D2V65 40 ) are employed for this simulation. The grid is N x × N y = 200 × 1, spatial step Δx = Δy = 5 × 10 −7 m, temporal step Δx = Δy = 1.25 × 10 −10 s. Figure 2 illustrates the simulation results and exact solutions during the chemical reaction in the free-falling box. Theoretically, the density remains constant, ρ = 1.30290kg ⋅ m −3 , the velocity changes as u y = gt, with g = −9.8m ⋅ s −2 , and the sum of internal energy and chemical heat remains constant, E = 2.59050 × 10 6 J ⋅ m −3 . As for the simulation, each model (D2V16, D2V24 37 , and D2V65 40 ) gives the density ρ = 1.30290kg ⋅ m −3 and the energy E = 2.59050 × 10 6 J ⋅ m −3 in the whole process, which coincide with the exact solutions. There are tiny differences between the simulation results and exact solutions of the velocity.  Furthermore, after the chemical reaction is completed, the adiabatic constant volume temperature is 2078K calculated by the three DBMs, while it is 2614K obtained by the standard LBM 12,41 , The parameters for the LBM in this work are the same as those in ref. 41 . Compared with the experimental datum 2080K 42 , the relative differences are 0.1% for the DBM and 25.7% for the standard LBM, respectively. Physically, the DBM is suitable for compressible systems with adjustable ratio of specific heats, while the LBM in refs 12,41 , can only be used for the case with constant pressure and fixed ratio of specific heats.
To discuss computational costs of various discrete velocity models, we keep a record of computing times required by the aforementioned simulation in Table 1. The computational facility is a personal computer with Intel(R) Core(TM) i7-6700K CPU @ 4.00 GHz and RAM 32.00 GB. There are 16, 24, and 65 (16, 24, and 16) discrete velocities (moment relations) in D2V16, D2V24 37 , and D2V65 40 , respectively. And the computing times are 1560 s, 3960 s, and 4980 s for the three models, respectively. Obviously, D2V24 and D2V65 models need larger RAM and longer time than D2V16 model. Flame at constant pressure. Let us simulate a flame at constant pressure. It travels with subsonic speed in a channel from left to right. In front of the flame is the propane-air mixture with equivalence ratio 0.6, particle number density 44.6 mol ⋅ m −3 , temperature 300 K, and pressure 1 atm. The grid is N x × N y = 2500 × 1, spatial step Δx = Δy = 2 × 10 −5 m, temporal step Δt = 1.25 × 10 −10 s. Figure 3 shows the evolution of ω ov (left) and Δ xx Moreover, in the DBM simulation, the pressure is close to 1atm around the flame, and the temperature is 1705K behind the flame, which is consistent with the experiment 2 , while the temperature is 2028K in the traditional LBM 12,41 , (see Fig. 4). The latter's relative error is 18.9% compared with the experimental result 2 . Physically, the ratio of specific heats in the DBM is tunable, while the one in the LBM in refs 12,41 , is fixed at 2. Besides, the chemical reaction does not affect the flow field in this LBM 41 , while the chemical reaction and fluid flow are naturally coupled in our DBM.

Shock wave.
A shock wave is a type of disturbance that propagates faster than the local speed of sound in a fluid with significant compressible effects. Its applications cover the fields of medicine, astrophysics, industrial engineering, etc. For example, it becomes effective medical treatment for kidney and ureteral stones. It can be used for cell transformation, preservative impregnation in bamboo, sandal oil extraction, and removal of micron size dust from silicon wafer surfaces 44 . To validate the DBM for high-speed compressible systems, we conduct the simulation of a shock wave. The wave propagates in the air from left L to right R. The initial field is,  The grid is N x × N y = 40000 × 1, spatial step Δx = Δy = 1 × 10 −8 m, temporal step Δt = 1.25 × 10 −12 s. To exhibit the capability of the DBM to study nonequilibrium behaviours, Fig. 6 shows the nonequilibrium manifestations around the shock wave. Figure 6 n T u x 2 in panels (a) and (c), respectively. Physically, the translational energy of oxygen (or nitrogen) in x degree of freedom travels faster than its equilibrium counterpart. Consequently, its departure from the equilibrium state is greater than zero around  the shock wave. Furthermore, there are few differences between the DBM and the approximate solution 38 in panels (b) and (d), respectively. Because the approximate solution is obtained by the first-order truncation of distribution function 38 . The simulation results are satisfactory.
Supersonic reacting wave. Supersonic reactive flows have been successfully used to deposite coating to a surface, clean equipment, mine for minerals, weld metals, etc. Numerical research of supersonic reacting wave has practical significance for the prevention of gas explosion in mining, flammable dust fires, and furnace burst in industry, etc. For the sake of verifying its suitability for supersonic reactive flows, the model is used to simulate a reacting wave. The initial field is divided into two parts. The right part is occupied by the premixed propane-air with equivalence ratio 0.524865, the left part by the chemical products. The reacting wave travels from left to right. And physical quantities satisfy the Hugoniot relations, i.e., The grid is N x × N y = 8000 × 1, spatial step Δx = Δy = 5 × 10 −7 m, temporal step Δt = 6.25 × 10 −11 s. Figure 7 displays the wave profiles: (a) ρ, The squares indicate DBM results, the lines analytic solutions of Zeldovich-Neumann-Doering (ZND) theory 2 . The DBM results behind the wave are (ρ,u x ,u y ,T) = (2.00166kg ⋅ m −3 ,666.356m ⋅ s −1 ,0m ⋅ s −1 ,2383.86K). Compared with the first row in Eq. (25), the relative differences are 0%, 0%, 0%, and 0.8%, respectively. Moreover, DBM gives the wave speed 1632 m/s, and the analytic solution is 1631.6 m/s. The relative difference between them is 0.02%. Additionally, there are slight differences between the theoretical and numerical results around the wave peak. Physically, the ZND theory assumes a sharp discontinuity at the wave peak and ignores the viscosity, heat conduction and other nonequilibrium effects 2 . On the other hand, the DBM takes into account the viscosity, heat conduction and other transport processes. Thus, the DBM is more reliable than the simple ZND theory.

Kelvin-Helmholtz instability.
To demonstrate that the DBM has a good ability of capturing interface deformation, we simulate a typical 2-D phenomenon, KH instability. The initial field, with area 0.6 m × 0.2 m, consists of two parts. The left part is full of propane with vertical speed 200 m · s −1 , while the right part is filled with air with −200 m · s −1 . Considering the transition layer between the two parts, the field jump at the interface is smoothed by a tanh profile with width 0.002 m. The temperature is 300 K, and pressure 1 atm. Between the propane and air is a sinusoidal interface with amplitude 0.003 m, which is used to promote the KH instability. Moreover, the outflow and periodic boundary conditions are adopted in the x and y directions, respectively. The grid is N x × N y = 3000 × 1000, spatial step Δx = Δy = 2 × 10 −4 m, temporal step Δt = 2.5 × 10 −8 s.  Figure 8 shows the molar fraction of propane at four different times. Initially, the interface starts to wrinkle due to the initial perturbation and velocity shear. A rolled-up vortex emerges after the initial linear growth stage. Then there is a large vortex around the interface. The evolution of the field is qualitatively similar to previous studies 38,45 , Moreover, Fig. 9 delineates the contour of pressure with velocity field, corresponding to Fig. 8(c). Compared Fig. 9 with Fig. 8(c), we can find that the minimum pressure, p = 0.257atm, is located at, (0.3004 m, 0.0750 m), the center of the vortex. While the maximum, p = 1.24atm, takes place at, (0.3008 m, 0.1708 m), where counterflows above the vortex encounter each other and the horizontal velocity is close to zero. Physically, the pressure gradient around the vortex provides the centripetal force for the rotating flow.
To quantitatively validate the results, we plot the logarithm of absolute value of the minimum perturbed horizontal velocity, ln|u x−min |, versus time, see Fig. 10. The squares are for DBM results, the solid line represents the fitting function F(t) =−3.89713 + 11.3302t, and the dashed line stands for the analytic solution = − . +  F t At ( ) 3 95868 with the growth rate = .  A 12 0995 46 . Here the nondimensionalization is used the same as ref. 46 . The relative difference of the growth rate between the fitting function and the analytic solution is 6.4%. Furthermore, we compare the simulation frequency 1256 Hz with the analytic solution 1248 Hz 38,46 . The relative difference is 0.6%. The difference mainly comes from the fact that the effects of compressibility, viscosity, and heat conduction are considered by the DBM, but ignored by the analytic theory 46 .

Conclusions
We present a reactive multi-component DBM in combination with a one-step overall chemical reaction. The effects of chemical reaction and external force are considered. A two-dimensional sixteen-velocity model D2V16  is proposed with adjustable parameters ( σ v a , σ v b , σ v c , σ v d ) controlling discrete velocities and (η σ a , η σ b , η σ c , η σ d ) for internal energies in extra degrees of freedom. The specific heat ratio of each species σ is flexible since extra degrees of freedom are taken into account. This model is suitable for premixed, nonpremixed or partially premixed combustion, from subsonic to supersonic fluid flows, in or out of equilibrium. Through the Chapman-Enskog multiscale analysis, the DBM recovers the modified NS equations for reactive species with external force effects in the hydrodynamic limit. In addition to the usual nonequilibrium terms (viscous stress and heat flux) in NS models, more detailed hydrodynamic and thermodynamic nonequilibrium quantities (high order kinetic moments and their departure from equilibrium) can be calculated in the DBM dynamically and conveniently. Since the DBM can provide detailed distributions of nonequilibrium quantities, it permits to assess the corresponding numerical predictions of NS models without considering the nonequilibrium effects. Hence, the DBM has the potential to offer more accurate information to help design devices operating in transient and/or extreme conditions away from equilibrium. Furthermore, thanks to its mesoscopic nature, the DBM could provide deeper insight into ubiquitous reactive or nonreactive fluid flows with a large span of spatial-temporal scales. Finally, due to its generality, the developed methodology is applicable to a wide range of phenomena across many energy technologies, emissions reduction, environmental protection, mining accident prevention, chemical and process industry.