A multi-component discrete Boltzmann model for nonequilibrium reactive flows
AA multi-component discrete Boltzmann model for nonequilibriumreactive flows
Chuandong Lin ∗1 , Kai Hong Luo †1,2 , Linlin Fei , and Sauro Succi Center for Combustion Energy, Key Laboratory for Thermal Science and Power Engineeringof Ministry of Education, Department of Thermal Engineering, Tsinghua University, Beijing100084, China Department of Mechanical Engineering, University College London, Torrington Place,London WC1E 7JE, UK Istituto Applicazioni Calcolo, CNR, Via dei Taurini 19, 00185 Rome, ItalyJune 30, 2020
Abstract
We propose a multi-component discrete Boltzmann model (DBM) for premixed, nonpremixed, or partiallypremixed nonequilibrium reactive flows. This model is suitable for both subsonic and supersonic flows withor without chemical reaction and/or external force. A two-dimensional sixteen-velocity model is constructedfor the DBM. In the hydrodynamic limit, the DBM recovers the modified Navier-Stokes equations forreacting species in a force field. Compared to standard lattice Boltzmann models, the DBM presents notonly more accurate hydrodynamic quantities, but also detailed nonequilibrium effects that are essential yetlong-neglected by traditional fluid dynamics. Apart from nonequilibrium terms (viscous stress and heatflux) in conventional models, specific hydrodynamic and thermodynamic nonequilibrium quantities (highorder kinetic moments and their departure from equilibrium) are dynamically obtained from the DBM ina straightforward way. Due to its generality, the developed methodology is applicable to a wide range ofphenomena across many energy technologies, emissions reduction, environmental protection, mining accidentprevention, chemical and process industry.
Introduction
Reactive flows are ubiquitous in nature and paramount to the sustainable development of society and ecologicalenvironment all over the world. For example, chemical energy released from fossil fuel in combustion comprisesover 80% the world’s energy utilization [1]. As the main human’s threats, atmospheric pollution, climatechange and global warming are directly relevant to harmful emissions from reactive flows, which involve abroad range of physicochemical phenomena, interacting over various spatial and temporal scales [2]. Besides,understanding of reactive flows is helpful to prevent fires in buildings, gas explosion in mines, burst in chemicalfactories, etc. Due to their significant importance to human society, reactive flows have attracted considerableattention in experimental, theoretical, and numerical fields. Actually, it is a challenging issue for traditionalmacroscopic or microscopic models to efficiently and accurately describe combustion phenomena where thespan of spatial-temporal scales is relatively large and nonequilibrium phenomena play essential roles [3, 4, 5]. Infact, the nonequilibrium effects always change physical quantities (such as density, velocity, temperature, etc)in the evolution of fluid systems away from equilibrium, especially in transient and/or extreme conditions. Apromising way to address this issue is to employ a mesoscopic kinetic model, lattice Boltzmann model (LBM),based on suitably simplified versions of the Boltzmann equation [6, 7, 8, 9, 10, 11, 12, 13, 14].Recently, LBM has emerged as a versatile tool to simulate various complex systems, including reactiveflows [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. Previous LBMs were successfullyemployed as solvers of macroscopic governing equations, such as hydro-chemical equations for incompressiblesystems with low Mach number [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. Physicalquantities (such as pressure, velocity, temperature) can be described separately by several distribution functionsin traditional LBMs, which are different from the Boltzmann equation where a single distribution functioncontains all information. For traditional LBMs, only a few low order kinetic moments of discrete equilibrium ∗ Electronic address: [email protected] ; Corresponding author † Electronic address:
[email protected] ; Corresponding author a r X i v : . [ phy s i c s . c o m p - ph ] J un igure 1: Sketch of the two-dimensional sixteen-velocity model.function are used, and the high order moments are not correctly reproduced [32]. This limitation resultsin the failure of recovering the complete Navier-Stokes (NS) equations and providing more information onnonequilibrium behaviours. To overcome those problems, one promising method is to resort to a variant oftraditional LBM, discrete Boltzmann model (DBM), where a required number of high order moments are satisfied[33, 34, 35, 36, 37, 38]. Different from traditional LBMs, DBM contains both equilibrium and nonequilibriumphysical quantities that stem from the same discrete distribution function [33, 34, 35, 36, 37, 38].Over the past years, the versatile DBM has been effectively applied to thermal phase separation, fluidinstabilities, reactive flows, etc. [33, 34, 35, 36, 37, 38] The DBM for reactive flows was firstly presented by Yanet 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 Prandtlnumber are adjustable [37]. The next year, a DBM is formulated for reactive flows where chemical reactantand product are described by two coupled distribution functions [38]. However, previous DBMs are suitable forpremixed reactive flows, but not for nonpremixed or partially premixed reactive flows [35, 36, 37, 38]. For thesake of simulating both subsonic and supersonic nonequilibrium reactive flows with premixed, nonpremixed, orpartially premixed reactants, we propose a multi-component DBM in this work. The DBM presents two ways toaccess the thermodynamic nonequilibrium behaviours. One is to measure the viscous stress and heat flux thatare described by traditional NS models; The other is to calculate the kinetic moments of the difference betweenequilibrium and nonequilibrium discrete distribution functions, which is beyond conventional hydrodynamicmodels. 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,C H + 5O → + 4H O , (1)where C H , O , CO , and H O denote propane, oxygen, carbon dioxide, and water, respectively. The stoi-chiometric coefficients for them are [ a C H , a O , a CO , a H O ] = [ − , − , , ω ov = k ov n C H n O exp ( − E a /RT ) , (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, ∂f σi ∂t + v σi · ∇ f σi = Ω σi + R σi + G σi , (3)2ith Ω σi = − τ σ [ f σi − f σeqi ( n σ , u , T )] , (4) R σi = 1 τ σ [ f σeqi ( n σ ∗ , u , T ∗ ) − f σeqi ( n σ , u , T )] , (5) G σi = 1 τ σ (cid:2) f σeqi (cid:0) n σ , u † σ , T † σ (cid:1) − f σeqi ( n σ , u σ , T σ ) (cid:3) . (6)Here f σi ( f σeqi ) 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 molecularcollision, chemical reaction and external force, respectively. The collision term in Eq. (4) obeys the conversionof 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 thephysical quantities of the species (mixture). In Eq. (5), n σ ∗ and T ∗ ( n σ and T ) denote the molar concentrationand temperature after (before) chemical reaction within time step τ σ . Similarly, in Eq. (6), the hydrodynamicvelocity changes from u σ to u † σ within time τ σ due to external force, meanwhile the temperature changes from T σ to T † σ . The discrete equilibrium function f σeqi is linked with the formula, f σeq = n σ (cid:18) m σ πT (cid:19) D/ (cid:18) m σ πI σ T (cid:19) / exp " − m σ ( v − u ) T − m σ η I σ T , (7)in the way that a required number of kinetic moments calculated by the integral of f σeq are equivalent tothose by the summation of f σeqi . In Eq. (7), m σ stands for molar mass, D = 2 space dimension, I σ extradegrees of freedom corresponding to molecular rotation or vibration. There are 7 moments satisfied by f σeqi = f σeqi ( n σ , u , T ) in this work. Specifically, n σ = X i f σeqi , (8) n σ u = X i f σeqi v σi , (9) n σ (cid:20) ( D + I σ ) Tm σ + u · u (cid:21) = X i f σeqi (cid:0) v σi · v σi + η σ i (cid:1) , (10) n σ (cid:18) δ αβ Tm σ e α e β + uu (cid:19) = X i f σeqi v σi v σi , (11) n σ u (cid:20) ( D + I σ + 2) Tm σ + u · u (cid:21) = X i f σeqi (cid:0) v σi · v σi + η σ i (cid:1) v σi , (12) n σ ( u α e β e χ δ βχ + e α u β e χ δ αχ + e α e β u χ δ αβ ) T σ m σ + n σ uuu = X i f σeqi v σi v σi v σi , (13) n σ δ αβ e α e β (cid:20) ( D + I σ + 2) Tm σ + u · u (cid:21) Tm σ + n σ uu (cid:20) ( D + I σ + 4) Tm σ + u · u (cid:21) = X i f σeqi (cid:0) v σi · v σi + η σ i (cid:1) v σi v σi , (14)which can be expressed in an uniform form ˆf eq = Mf eq , where f eq = (cid:0) f σeq , f σeq , · · · , f σeqN i (cid:1) T , ˆf eq = (cid:16) ˆ f σeq , ˆ f σeq , · · · , ˆ f σeqN i (cid:17) T , M = M σ M σ · · · M σ N i M σ M σ · · · M σ N i ... ... . . . ... M σN i M σN i · · · M σN i N i , with N i = 16. The elements of ˆf eq are ˆ f σeq = n σ , ˆ f σeq = n σ u x , ˆ f σeq = n σ u y , ˆ f σeq = n σ [( D + I σ ) T /m σ + u ],ˆ f σeq = n σ ( T /m σ + u x ), ˆ f σeq = n σ u x u y , ˆ f σeq = n σ ( T /m σ + u y ), ˆ f σeq = n σ u x [( D + I σ + 2) T /m σ + u ], ˆ f σeq = n σ u y [( D + I σ + 2) T /m σ + u ], ˆ f σeq = 3 n σ u x T /m σ + n σ u x , ˆ f σeq = n σ u y T /m σ + n σ u x u y , ˆ f σeq = n σ u x T /m σ + n σ u x u y , ˆ f σeq = 3 n σ u y T /m σ + n σ u y , ˆ f σeq = n σ [( D + I σ + 2) T /m σ + u ] T /m σ + n σ u x [( D + I σ + 4) T /m σ + u ],ˆ f σeq = n σ u x u y [( D + I σ +4) T /m σ + u ], ˆ f σeq = n σ [( D + I σ +2) T /m σ + u ] T /m σ + n σ u y [( D + I σ +4) T /m σ + u ],3nd those of M are M σ i = 1, M σ i = v σix , M σ i = v σiy , M σ i = v σ i + η σ i , M i = v σ ix , M i = v σix v σiy , M i = v σ iy , M i = ( v σ i + η σ i ) v σix , M i = ( v σ i + η σ i ) v σix , M i = v σ ix , M i = v σ ix v σiy , M i = v σix v σ iy , M i = v σ iy , M i = ( v σ i + η σ i ) v σ ix , M i = ( v σ i + η σ i ) v σix v σiy , M i = ( v σ i + η σ i ) v σ iy .The discrete velocities v σi and η σi are (see Fig. 1),[ v σi , η σi ] = [cyc : v σa ( ± , , η σa ] 1 ≤ i ≤ , [cyc : v σb ( ± , ± , η σb ] 5 ≤ i ≤ , [cyc : v σc ( ± , , η σc ] 9 ≤ i ≤ , [cyc : v σd ( ± , ± , η σd ] 13 ≤ i ≤ , (15)where cyc indicates the cyclic permutation, ( v σa , v σb , v σc , v σd , η σa , η σb , η σc , η σd ) are adjustable parameters. And thediscrete equilibrium distribution is expressed by f eq = M − ˆf eq . (16)One significant capability of the DBM is to investigate nonequilibrium manifestations by measuring thefollowing physical variables, ∆ σ = X i m σ ( f σi − f σeqi ) v σi v σi , (17)∆ σ , = X i m σ ( f σi − f σeqi ) (cid:0) v σi · v σi + η σ i (cid:1) v σi , (18)∆ σ = X i m σ ( f σi − f σeqi ) v σi v σi v σi , (19)∆ σ , = X i m σ ( f σi − f σeqi ) (cid:0) v σi · v σi + η σ i (cid:1) v σi v σi . (20)Mathematically, ∆ σ = ∆ σ αβ e α e β is a second order tensor with four components, among which only three (∆ σ xx ,∆ σ xy , ∆ σ yy ) are independent; ∆ σ = ∆ σ αβγ e α e β e γ is a third order tensor with eight components where only four(∆ σ xxx , ∆ σ xxy , ∆ σ xyy , ∆ σ yyy ) are independent; ∆ σ , = ∆ σ , α e α is the first order tensor (i.e. vector) contractedfrom a third order tensor and have two independent components; Similarly for ∆ σ , = ∆ σ , αβ e α e β .It is easy to prove that, via the Chapman-Enskog multiscale analysis, the DBM is in line with the followingmodified NS equations, ∂ρ σ ∂t + ∂∂r α ( ρ σ u σα ) = ρ σ , (21) ∂∂t ( ρ σ u σα ) + ∂∂r β (cid:0) δ αβ p σ + ρ σ u σα u σβ + ∆ σ αβ (cid:1) = ρ σ a α − ρ σ τ σ ( u σα − u α ) + ρ σ u α , (22) ∂∂t ρ σ (cid:18) e σ + 12 u σ (cid:19) + ∂∂r α (cid:20) ρ σ u σα (cid:18) e σ + 12 u σ (cid:19) + p σ u σα + ∆ σ , α (cid:21) = ρ σ u σα a α + ρ σ (cid:18) D + I σ Tm σ + 12 u (cid:19) + D + I σ m σ ρ σ T − ρ σ τ σ (cid:18) D + I σ T σ − Tm σ + u σ − u (cid:19) , (23)in the hydrodynamic limit, with ∆ σ αβ = P σαβ + U σαβ , ∆ σ , α = − κ σ ∂∂r α ( D + I σ T σ m σ ) + u σβ P σαβ − X σα − Z σα , P σαβ = − µ σ ( ∂u σα ∂r β + ∂u σβ ∂r α − δ αβ D + I σ ∂u σχ ∂r χ ), U σαβ = − ( ρ σ + ρ σ τ σ )( δ αβ u σ + u − u σχ u χ D + I σ + u α u σβ + u σα u β − u σα u σβ − u α u β ), X σα = τ σ ( D + I σ + 2) ρ σ ( u σα − u α ) m σ T , Z σα = ( ρ σ + τ σ ρ σ ) { u σα D + I σ ( u σβ − u β ) − u σα − u α [ D + I σ +2 m σ ( T σ − T ) + u σ − u ] } , where p σ = n σ T σ indicates pressure, e σ = ( D + I σ ) T σ / (2 m σ ) internal energy per unit mass, µ σ = p σ τ σ dynamicviscosity coefficient, κ σ = γ σ µ σ heat conductivity, and γ σ = ( D + I σ + 2) / ( D + I σ ) specific heat ratio. Thesuperscript “ " represents the change rate of physical quantities due to the chemical reaction. In fact, applyingthe operator P σ to both sides of Eqs. (21) − (23) gives NS equations for the whole system, which reducesto conventional NS equations when u σ = u and T σ = T . Obviously, Eqs. (21) − (23) gives a more detaileddescription 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 importantthermodynamic nonequilibrium manifestations or physical effects on fluid flows. In fact, a more detailed wayto study the nonequilibrium behaviours is to investigate the departure of high order velocity moments fromtheir local equilibrium counterparts, as shown in Eqs. (17)-(20). Those kinetic moments of the differencebetween nonequilibrium and equilibrium distribution functions have significant physical meanings. In particular,∆ σ is associated with viscous stress tensor and nonorganised momentum fluxes, ∆ σ , and ∆ σ are related tononorganised energy (heat) fluxes, ∆ σ , corresponds to the flux of nonorganised energy (heat) flux [36, 37]. Theterminology “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]. Moreover, M αα ( f σ ) = P i m σ f σi v σ iα is defined as the translational energy of species σ in α direction, M αα ( f σeq ) is4igure 2: The average physical quantities of the reactive mixture in the free-falling box versus time t : (a) massdensity ρ , (b) vertical speed u y , (c) the sum of internal energy and chemical heat E . The symbols (squares,circles, and triangles) denote simulation results (D2V16, D2V24, and D2V65), the solid lines exact solutions.its equilibrium counterpart, and ∆ σ αα = M αα ( f σ ) − M αα ( f σeq ) the nonorgnised energy of species σ in α direction; M , α ( f σeq ) = P i m σ f σeqi ( v σ i + η σ i ) v σiα refers to the organised flux of energy in α direction, ∆ σ , α the nonorganised flux of energy in α direction, and M , α ( f σ ) = M , α ( f σeq ) + ∆ σ , α the total fluxof energy in α direction. Obviously, DBM provides more nonequilibrium information on various species in fluidflows, which is an essential advantage over traditional models. 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 theexternal force effects are demonstrated. The second test is a subsonic flame at constant pressure. In the thirdpart, 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 simulatedin 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 thepurpose of accuracy and robustness, it is preferable to set the values of discrete velocities ( v σa , v σb , v σc , v σd ) aroundthe values of | u | and p DT /m σ , and choose ( η σa , η σb , η σc , η σd ) around the value p I σ T /m σ , which is reasonableon account of Eqs. (8)-(10). Combustion in constant volume
First of all, we simulate the combustion of propane-air filled in a free-falling box, which consists of threeparts with volumes V , V and V , respectively. The fixed volume of the box is V = V + V + V , and V : V : V = 3 : 119 : 78. Initially, the left part is filled with propane, the middle part is full of air, and theright part is occupied by the propane-air mixture with equivalence ratio 0 .
6. In each part, the particle numberdensity is 40 . · m − , temperature 300 K, and pressure 1 atm. Premixed, nonpremixed and partially premixedcombustion phenomena take place simultaneously in this box after ignition. Specifically, the nonpremixedcombustion takes place between the left and middle parts, the partially premixed combustion occurs betweenthe middle and right parts with a changing equivalence ratio, and the premixed combustion is in the rightmostpart 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 × − m, temporalstep ∆ x = ∆ y = 1 . × − s.Figure 2 illustrates the simulation results and exact solutions during the chemical reaction in the free-fallingbox. Theoretically, the density remains constant, ρ = 1 . · m − , the velocity changes as u y = gt , with g = − . · s − , and the sum of internal energy and chemical heat remains constant, E = 2 . × J · m − . Asfor the simulation, each model (D2V16, D2V24 [37], and D2V65 [40]) gives the density ρ = 1 . · m − andthe energy E = 2 . × J · m − in the whole process, which coincide with the exact solutions. There are tinydifferences between the simulation results and exact solutions of the velocity. For example, at time t = 3 × − s, the three models (D2V16, D2V24 [37], and D2V65 [40]) give simulation results u y = − . × − , − . × − , and − . × − m · s − , respectively. Compared to the exact value u y = − . × − m · s − ,5able 1: Computing times for simulation of combustion in constant volume with various models.Model Number of discrete velocities Number of moment relations Computing timeD2V16 16 16 1560 sD2V24 24 24 3960 sD2V65 65 16 4980 sFigure 3: The flame profiles at constant pressure at times t = 0 . s , t = 0 . s , and t = 0 . s . The leftpanel is for the overall reaction rate ω ov , and the right for ∆ O xx which is two times the departure of translationalenergy of O in x direction from its equilibrium counterpart.their relative errors are 0 . . . K calculated by the three DBMs, while it is 2614 K obtained by the standard LBM [12, 41]. The parameters for theLBM in this work are the same as those in Ref. [41]. Compared with the experimental datum 2080 K [42], therelative differences are 0 .
1% for the DBM and 25 .
7% for the standard LBM, respectively. Physically, the DBMis suitable for compressible systems with adjustable ratio of specific heats, while the LBM in Refs. [12, 41] canonly 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 timesrequired by the aforementioned simulation in Table 1. The computational facility is a personal computer withIntel(R) Core(TM) i7-6700K CPU @ 4.00GHz 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 computingtimes are 1560 s, 3960 s, and 4980 s for the three models, respectively. Obviously, D2V24 and D2V65 modelsneed 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. Infront of the flame is the propane-air mixture with equivalence ratio 0 .
6, particle number density 44 . · m − ,temperature 300 K, and pressure 1 atm. The grid is N x × N y = 2500 ×
1, spatial step ∆ x = ∆ y = 2 × − m,temporal step ∆ t = 1 . × − s.Figure 3 shows the evolution of ω ov (left) and ∆ O xx (right) versus x . The peak of ω ov corresponding to themost active chemical reaction is ahead of the trough of ∆ O xx where nonequilibrium manifestations are intenseand physical gradients are sharp. Note that the nonequilibrium manifestations can be employed to capture theflame or other interfaces [34]. The flame speed, 0 . m / s , can be obtained from the profiles of either ω ov or ∆ O xx .And the flow velocity is 0 . m / s in front of the flame. Hence, the burning velocity is (0 . − . m / s = . / s ,which equals the experimental result 0 . m / s [43]. While the standard LBM [41] gives a relative error, 9 . atm around the flame, and the temperature is1705 K behind the flame, which is consistent with the experiment [2], while the temperature is 2028 K in thetraditional 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 chemicalreaction and fluid flow are naturally coupled in our DBM.6igure 4: The flame temperature at constant pressure. The solid (dashed) line denotes DBM (LBM) results,the squares experimental data.Figure 5: The profiles of a shock wave: (a) mass density ρ , (b) velocity u x , (c) temperature T , (d) pressure p .The squares represent DBM results, the lines exact solutions. Shock wave
A shock wave is a type of disturbance that propagates faster than the local speed of sound in a fluid withsignificant 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 celltransformation, preservative impregnation in bamboo, sandal oil extraction, and removal of micron size dustfrom silicon wafer surfaces [44]. To validate the DBM for high-speed compressible systems, we conduct thesimulation of a shock wave. The wave propagates in the air from left L to right R . The initial field is, (cid:26) ( ρ, u x , u y , T ) L = (cid:0) . · m − , . · s − , · s − , . (cid:1) ( ρ, u x , u y , T ) R = (cid:0) . · m − , · s − , · s − , (cid:1) . (24)The grid is N x × N y = 40000 ×
1, spatial step ∆ x = ∆ y = 1 × − m, temporal step ∆ t = 1 . × − s. Figure 5 plots the profiles of the shock: (a) ρ , (b) u x , (c) T , (d) p . The squares denote DBM results, thelines exact solutions. The DBM results behind the shock wave are ( ρ, u x , u y , T ) = (1 . · m − , . · s − , · s − , . x direction M O xx ( f ), its equilibrium counterpart M O xx ( f eq ), and the exact solution n O T + ρ O u x . Figure 6 (b) illustrates7igure 6: Nonequilibrium manifestations around the shock wave: (a) translational energy of oxygen in x direction, (b) departure of translational energy of oxygen in x direction from equilibrium state, (c) translationalenergy of nitrogen in x direction, (d) departure of translational energy of nitrogen in x direction from equilibriumstate.the departure of translational energy of oxygen in x degree of freedom from its equilibrium state ∆ O xx . Figures6 (c) and (d) are for nitrogen. It is clear that M σ xx ( f eq ) coincides with the solution n σ T + ρ σ u x in panels (a)and (c), respectively. Physically, the translational energy of oxygen (or nitrogen) in x degree of freedom travelsfaster than its equilibrium counterpart. Consequently, its departure from the equilibrium state is greater thanzero around the shock wave. Furthermore, there are few differences between the DBM and the approximatesolution [38] in panels (b) and (d), respectively. Because the approximate solution is obtained by the first-ordertruncation 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, minefor minerals, weld metals, etc. Numerical research of supersonic reacting wave has practical significance for theprevention of gas explosion in mining, flammable dust fires, and furnace burst in industry, etc. For the sake ofverifying its suitability for supersonic reactive flows, the model is used to simulate a reacting wave. The initialfield is divided into two parts. The right part is occupied by the premixed propane-air with equivalence ratio0 . (cid:26) ( ρ, u x , u y , T ) L = (cid:0) . · m − , . · s − , · s − , . (cid:1) ( ρ, u x , u y , T ) R = (cid:0) . · m − , · s − , · s − , (cid:1) . (25)The grid is N x × N y = 8000 ×
1, spatial step ∆ x = ∆ y = 5 × − m, temporal step ∆ t = 6 . × − s.Figure 7 displays the wave profiles: (a) ρ , (b) u x , (c) T , (d) p . The squares indicate DBM results, thelines analytic solutions of Zeldovich-Neumann-Doering (ZND) theory [2]. The DBM results behind the waveare ( ρ, u x , u y , T ) = (2 . · m − , . · s − , · s − , . . / s, and the analytic solution is 1631 . / s. The relative difference between them is 0 . ρ , (b) velocity u x , (c) temperature T , (d)pressure p . The squares stand for DBM results, the lines exact solutions.Figure 8: The molar fraction of propane in the evolution of KH instability at various times: (a) 0 s, (b) 5 × − s, (c) 10 − s, (d) 2 × − s.Figure 9: Pressure and velocity fields in the evolution of KH instability at time 10 − s.9igure 10: Evolution of the velocity perturbation, ln | u x − min | . The squares stand for DBM results, the solid(dashed) line is for the fitting function (analytic solution). Kelvin-Helmholtz instability
To demonstrate that the DBM has a good ability of capturing interface deformation, we simulate a typical 2-Dphenomenon, KH instability. The initial field, with area 0 . m × . , consists of two parts. The left part is fullof propane with vertical speed 200 m · s − , while the right part is filled with air with − m · s − . Consideringthe transition layer between the two parts, the field jump at the interface is smoothed by a tanh profile withwidth 0 . m . The temperature is 300 K, and pressure 1 atm. Between the propane and air is a sinusoidalinterface with amplitude 0 . m , which is used to promote the KH instability. Moreover, the outflow and periodicboundary conditions are adopted in the x and y directions, respectively. The grid is N x × N y = 3000 × x = ∆ y = 2 × − m, temporal step ∆ t = 2 . × − s.Figure 8 shows the molar fraction of propane at four different times. Initially, the interface starts to wrinkledue to the initial perturbation and velocity shear. A rolled-up vortex emerges after the initial linear growthstage. Then there is a large vortex around the interface. The evolution of the field is qualitatively similar toprevious studies [38, 45]. Moreover, Fig. 9 delineates the contour of pressure with velocity field, correspondingto Fig. 8 (c). Compared Fig. 9 with Fig. 8 (c), we can find that the minimum pressure, p = 0 . atm , is locatedat, (0 . m , 0 . m ), the center of the vortex. While the maximum, p = 1 . atm , takes place at, (0 . m ,0 . 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 perturbedhorizontal velocity, ln | u x − min | , versus time, see Fig. 10. The squares are for DBM results, the solid linerepresent the fitting function F ( t ) = − . . t , and the dashed line stands for the analytic solution F ( t ) = − . At with the growth rate ˙ A = 12 . . . Conclusions
We present a reactive multi-component DBM in combination with a one-step overall chemical reaction. Theeffects of chemical reaction and external force are considered. A two-dimensional sixteen-velocity model D2V16is 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 sinceextra degrees of freedom are taken into account. This model is suitable for premixed, nonpremixed or partiallypremixed 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 externalforce effects in the hydrodynamic limit. In addition to the usual nonequilibrium terms (viscous stress andheat flux) in NS models, more detailed hydrodynamic and thermodynamic nonequilibrium quantities (highorder kinetic moments and their departure from equilibrium) can be calculated in the DBM dynamically andconveniently. Since the DBM can provide detailed distributions of nonequilibrium quantities, it permits to assessthe 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 transient10nd/or extreme conditions away from equilibrium. Furthermore, thanks to its mesoscopic nature, the DBMcould 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 ofphenomena across many energy technologies, emissions reduction, environmental protection, mining accidentprevention, chemical and process industry. References [1] S. Chu, A. Majumdar, Opportunities and challenges for a sustainable energy future, Nature (London) 488(2012) 294–303.[2] C. K. Law, Combustion physics, Cambridge University Press, Cambridge, 2006.[3] H. Liu, Y. Zhang, W. Kang, P. Zhang, H. Duan, X. He, Molecular dynamics simulation of strong shock wavespropagating in dense deuterium, taking into consideration effects of excited electrons, Physical Review E95 (2) (2017) 023201.[4] Y. Ju, Recent progress and challenges in fundamental combustion research, Adv. Mech. 44 (2014) 1–72.[5] E. Nagnibeda, E. Kustova, Non-equilibrium reacting gas flows: kinetic theory of transport and relaxationprocesses, Springer, Berlin, 2009.[6] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, NewYork, 2001.[7] V. Sofonea, R. F. Sekerka, Bgk models for diffusion in isothermal binary fluid systems, Physica A: StatisticalMechanics and its Applications 299 (3) (2001) 494–520.[8] A. Montessori, P. Prestininzi, M. La Rocca, G. Falcucci, S. Succi, Lattice kinetic approach to non-equilibrium flows, in: AIP Conference Proceedings, Vol. 1738, AIP Publishing, 2016, p. 090005.[9] Y. Zhang, R. Qin, D. R. Emerson, Lattice boltzmann simulation of rarefied gas flows in microchannels,Physical Review E 71 (4) (2005) 047702.[10] G. H. Tang, Y. H. Zhang, D. R. Emerson, Lattice boltzmann models for nonequilibrium gas flows, PhysicalReview E 77 (4) (2008) 046701.[11] J. Meng, Y. Zhang, X. Shan, Multiscale lattice boltzmann approach to modeling gas flows, Physical ReviewE 83 (4) (2011) 046701.[12] X. He, L.-S. Luo, Lattice boltzmann model for the incompressible navier–stokes equation, Journal ofstatistical Physics 88 (3) (1997) 927–944.[13] R. Benzi, M. Sbragaglia, S. Succi, M. Bernaschi, S. Chibbaro, Mesoscopic lattice boltzmann modeling ofsoft-glassy systems: theory and simulations, The Journal of Chemical Physics 131 (10) (2009) 104903.[14] Q. Li, K. Luo, Q. Kang, Y. He, Q. Chen, Q. Liu, Lattice boltzmann methods for multiphase flow andphase-change heat transfer, Progress in Energy and Combustion Science 52 (2016) 62–105.[15] S. Ponce Dawson, S. Chen, G. D. Doolen, Lattice boltzmann computations for reaction-diffusion equations,The Journal of chemical physics 98 (2) (1993) 1514–1523.[16] D. H. Zanette, Interplay of reaction and transport in a perfect fluid, Physical Review E 50 (2) (1994) 1171.[17] Y. Qian, S. Orszag, Scalings in diffusion-driven reaction a+b → c: Numerical simulations by lattice bgkmodels, Journal of Statistical Physics 81 (1) (1995) 237–253.[18] J. R. Weimar, J. P. Boon, Nonlinear reactions advected by a flow, Physica A: Statistical Mechanics andits Applications 224 (1-2) (1996) 207–215.[19] Z. Tian, H. Xing, Y. Tan, S. Gu, S. D. Golding, Reactive transport lbm model for co 2 injection in fracturedreservoirs, Computers & Geosciences 86 (2016) 15–22.[20] S. Succi, G. Bella, F. Papetti, Lattice kinetic theory for numerical combustion, J. Sci. Comput. 12 (1997)395–408.[21] O. Filippova, D. Hänel, A novel numerical scheme for reactive flows at low mach numbers, Comput. Phys.Commun. 129 (2000) 267–274. 1122] H. Yu, L. S. Luo, S. S. Girimaji, Scalar mixing and chemical reaction simulations using lattice boltzmannmethod, Int. J. Comput. Eng. Sci. 3 (2002) 73–87.[23] K. Yamamoto, N. Takada, M. Misawa, Combustion simulation with lattice boltzmann method in a three-dimensional porous structure, Proc. Comb. Inst. 30 (2005) 1509–1515.[24] T. Lee, C. Lin, L. D. Chen, A lattice boltzmann algorithm for calculation of the laminar jet diffusion flame,J. Comput. Phys. 215 (2006) 133–152.[25] E. Chiavazzo, I. V. Karlin, A. N. Gorban, K. Boulouchos, Efficient simulations of detailed combustion fieldsvia the lattice boltzmann method, Int. J. Numer. Methods Heat Fluid Flow 21 (2011) 494–517.[26] S. Chen, J. Mi, H. Liu, C. Zheng, First and second thermodynamic-law analyses of hydrogen-air counter-flow diffusion combustion in various combustion modes, Int. J. of Hydrogen Energ. 37 (2012) 5234–5245.[27] S. Succi, O. Filippova, G. Smith, E. Kaxiras, Applying the lattice boltzmann equation to multiscale fluidproblems, Computing in Science & Engineering 3 (6) (2001) 26–37.[28] K. Furtado, J. Yeomans, Lattice boltzmann simulations of phase separation in chemically reactive binaryfluids, Physical Review E 73 (6) (2006) 066124.[29] M. Ashna, M. H. Rahimian, A. Fakhari, Extended lattice boltzmann scheme for droplet combustion,Physical Review E 95 (5) (2017) 053301.[30] G. Falcucci, S. Succi, A. Montessori, S. Melchionna, P. Prestininzi, C. Barroo, D. C. Bell, M. M. Biener,J. Biener, B. Zugic, et al., Mapping reactive flow patterns in monolithic nanoporous catalysts, Microfluidicsand Nanofluidics 20 (7) (2016) 1–13.[31] G. Falcucci, G. Amati, V. K. Krastev, A. Montessori, G. S. Yablonsky, S. Succi, Heterogeneous catalysisin pulsed-flow reactors with nanoporous gold hollow spheres, Chemical Engineering Science 166 (2017)274–282.[32] J. Kang, N. I. Prasianakis, J. Mantzaras, Thermal multicomponent lattice boltzmann model for catalyticreactive flows, Phys. Rev. E 89 (2014) 063310.[33] Y. Gan, A. Xu, G. Zhang, S. Succi, Discrete boltzmann modeling of multiphase flows: hydrodynamic andthermodynamic non-equilibrium effects, Soft Matter 11 (26) (2015) 5336–5345.[34] H. Lai, A. Xu, G. Zhang, Y. Gan, Y. Ying, S. Succi, Nonequilibrium thermohydrodynamic effects on therayleigh-taylor instability in compressible flows, Physical Review E 94 (2) (2016) 023106.[35] B. Yan, A. Xu, G. Zhang, Y. Ying, H. Li, Lattice boltzmann model for combustion and detonation, Frontiersof Physics 8 (1) (2013) 94–110.[36] C. Lin, A. Xu, G. Zhang, Y. Li, Polar coordinate lattice boltzmann kinetic modeling of detonation phe-nomena, Communications in Theoretical Physics 62 (5) (2014) 737.[37] A. Xu, C. Lin, G. Zhang, Y. Li, Multiple-relaxation-time lattice boltzmann kinetic model for combustion,Physical Review E 91 (4) (2015) 043306.[38] C. Lin, A. Xu, G. Zhang, Y. Li, Double-distribution-function discrete boltzmann model for combustion,Combustion and Flame 164 (2016) 137–151.[39] H. Zhang, F. Zhuang, Nnd schemes and their applications to numerical simulation of two-and three-dimensional flows, Advances in Applied Mechanics 29 (1991) 193–256.[40] M. Watari, Finite difference lattice boltzmann method with arbitrary specific heat ratio applicable tosupersonic flow simulations, Physica A: Statistical Mechanics and its Applications 382 (2) (2007) 502–522.[41] K. Yamamoto, X. He, G. D. Doolen, Simulation of combustion field with lattice boltzmann method, J.Stat. Phys. 107 (2002) 367–383.[42] R. Zalosh, Flammable Gas and Vapor Explosions[M]//SFPE Handbook of Fire Protection Engineering,Springer, New York, 2016.[43] I. Yamaoka, H. Tsuji, Determination of burning velocity using counterflow flames, in: Symposium (Inter-national) on Combustion, Vol. 20, Elsevier, 1985, pp. 1883–1892.1244] G. Jagadeesh, Industrial applications of shock waves, Proceedings of the Institution of Mechanical Engi-neers, Part G: Journal of Aerospace Engineering 222 (5) (2008) 575–583.[45] Y. Gan, A. Xu, G. Zhang, Y. Li, Lattice boltzmann study on kelvin-helmholtz instability: Roles of velocityand density gradients, Physical Review E 83 (5) (2011) 056704.[46] L. Wang, W. Ye, Y. Li, Combined effect of the density and velocity gradients in the combination ofkelvin-helmholtz and rayleigh-taylor instabilities, Phys. Plasmas 17 (2010) 042103. Acknowledgements
The authors thank Profs. Aiguo Xu and Guangcai Zhang for their helpful suggestions. This work is supportedby the Natural Science Foundation of China (NSFC) under Grant No. 91441120 and the Center for CombustionEnergy at Tsinghua University. Support from the UK Engineering and Physical Sciences Research Council underthe project “UK Consortium on Mesoscale Engineering Sciences (UKCOMES)” (Grant No. EP/L00030X/1) isgratefully acknowledged.
Author contributions
S.S. contributed to the analysis of the theory and results of the DBM; L.F. provided the code of standandLBMs and assisted in performing the simulations; K.L. initiated and supervised the study; C.L. developed andvalidated the DBM as well as drafted the manuscript. All authors modified and approved the manuscript.
Additional information