MRT discrete Boltzmann method for compressible exothermic reactive ﬂows

An eﬃcient, accurate and robust multiple-relaxation-time (MRT) discrete Boltzmann method (DBM) is proposed for compressible exothermic reactive ﬂows, with both speciﬁc heat ratio and Prandtl number being ﬂexible. The chemical reaction is coupled with the ﬂow ﬁeld naturally and the external force is also incorporated. An eﬃcient discrete velocity model which has sixteen discrete velocities (and kinetic moments) is introduced into the DBM. With both hydrodynamic and thermodynamic nonequilibrium ef- fects under consideration, the DBM provides more detailed and accurate information than traditional Navier–Stokes equations. This method is suitable for ﬂuid ﬂows ranging from subsonic, to supersonic and hypersonic ranges. It is validated by various benchmarks. http://creativecommons.org/licenses/by/4.0/


Introduction
Exothermic reactive flows are commonplace in nature and industry which play significant roles in economic and social development all over the world. In fact, more than 80% utilizable energy is transformed through exothermic reactive phenomena in the world [1] . On the other hand, they are associated with environmental problems, accidents or even disasters. For example, atmospheric pollution, global warming and climate change are closely linked to harmful emissions from reactive flows. In particular, fire hazards, which often induce explosion and shock, may cause huge danger and damage to human life, property and environment. Although considerable researches have been devoted to these fields, there are still many open issues due to their complexity. To be specific, they have a wide span of physicochemical phenomena, interact over various spatio-temporal scales, and involve various hydrodynamic and thermodynamic nonequilibrium behaviours [2][3][4] . Especially, for a spacecraft flying from the earth surface to outer space, where the chemical reaction and gravity exist, it covers a wide range of Knudsen numbers and various essential nonequilibrium phenomena. To describe such complex systems, traditional macroscopic models have the benefit of high computing efficiency, but could not capture detailed information accurately. While microscopic models have the merit of an accurate and full description, they encounter spatio-temporal constraints because of their high computing costs.
At the mesoscopic level, the lattice Boltzmann method (LBM) may overcome aforementioned problems [5][6][7][8][9][10][11][12][13][14][15][16] . In the past three decades, the LBM has achieved significant success in the simulation of complex systems, including reactive flows [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] . The traditional LBM usually works as an alternative tool to solve macroscopic equations, such as incompressible Navier-Stokes (NS) equations. Various physical quantities, such as flow velocity and temperature, may be described by different sets of the discrete distribution function. Recently, a novel variant of LBM, discrete Boltzmann method (DBM), has emerged as an efficient kinetic model to capture both hydrodynamic and thermodynamic nonequilibrium effects in fluid flows [36,37] . Different from traditional LBMs, the DBM employs only one set of discrete distribution function to describe various physical quantities, including the density, temperature, velocity, and other high order kinetic moments, which is in line with the Boltzmann equation. Since 2013, several Single-Relaxation-Time DBMs have been formulated for exothermic reactive flows [38][39][40] . Yet, the Prandtl number in those proposed model is fixed at Pr = 1 . To overcome this, a multiple-relaxation-time (MRT) DBM was presented [41] . There are 24 independent kinetic moments satisfied by 24 discrete equilibrium distribution functions in this work [41] . These kinetic moments are necessary for the DBM to recover the reactive NS equations in the hydrodynamic limit [41] . Besides, the effects of external force are neglected in this model [41] . However, external forces (such as gravity) often have essential influences upon reactive flows. In the present work, we introduce a new form of reaction and force terms, and reduce the 24 kinetic moments (and discrete equilibrium distribution functions) to only 16 while the recovery of the NS equations is made as well. Besides its practical value as an efficient computational tool for the traditional dynamics of complex systems, this model also provides details of nonequilibrium behaviours dynamically and conveniently. We describe the DBM in Section 2 , validate it in Section 3 , and finally summarize this work in Section 4 .

Discrete Boltzmann method
The DBE takes the form, Here · · · f eq N ) T denote discrete distribution functions and their equilibrium counterparts, re- represent kinetic moments of discrete distribution function and their equilibrium counterparts, respectively. M −1 is the inverse matrix of M , and M is a square matrix, see Appendix A .
, N and N = 16 . As shown in Fig. 1 with tunable parameters v a and v b controlling the value of v i .
The artificial term ˆ A = (0 · · · 0 ˆ A 8 ˆ A 9 0 · · · 0) T is used to modify the collision operator The reason for this modification is as follows. Although the tunable relaxation coefficients ˆ S i seem mathematically independent of each other, coupling may exist among the relaxation processes of various kinetic modes ( ˆ f ne i = ˆ f i −ˆ f eq i ) from the physical point of view. For the sake of correct description of macroscopic behaviours, we should perform the Chapman-Enskog expansion, analyze the consistency of nonequilibrium transportation terms in the recovered hydrodynamic equations, and find a solution for the modification to the collision term. In short, this modification is incorporated in the DBM to recover the consistent NS equations in the hydrodynamic limit, see Appendix A . The artificial term is the function of the velocities ( u x , u y ) and the first-order partial derivatives of them with respect to x or y . These derivatives can be solved by various finite difference schemes. In this work, the central difference scheme is adopted. For example, at the node ( i x , i y ). Numerical tests demonstrate that the artificial term does not induce significant numerical problems. Furthermore, the artificial term can be removed for the case ˆ S 5 = ˆ S 8 and ˆ S 7 = ˆ S 9 .
The force and reaction terms, Mathematically, the difference of the equilibrium distribution functions over a small time interval is an approximation to the change rate of distribution functions, based on the assumption f i ≈ f eq i . The physical reason for Eq. (6) is as follows. It is regarded that neither external force nor chemical reaction changes the density ρ.
The external force affects the hydrodynamic velocity u with acceleration a . Consequently, the velocity changes from u into u + a τ within a small time interval τ due to the external force. Meanwhile, the temperature changes into T + τ T on account of the chemical reaction. Specifically, the change rate of energy is because of the external force and chemical reaction. From Eq. (7) and the definition E = D + I 2 ρT + 1 2 ρu · u , we obtain the change rate of temperature where D = 2 stands for the number of dimensions, I the number of extra degrees of freedom corresponding to molecular rotation and/or internal vibration. The reaction process λ is defined as the mass ratio of the chemical product to mixture. The chemical reaction is controlled by the Cochran's rate function which depends upon the pressure, p = ρT , in terms of adjustable parameters ω 1 , ω 2 , m and n [42] . Here λ is defined as the local mass fraction of the reaction product. Without loss of generality, we choose ( ω 1 , ω 2 , m, n ) = (2, 100, 2, 2.5), and employ the ignition temperature T ig = 1 . 1 in this work. Only when T > T ig can the chemical reaction take place.
For the sake of recovering the NS equations, the discrete equilibrium distribution function should satisfy the following relations f eq i = ρ, (10) f eq f eq In contrast, all the amplification factors are identical in the SRT model, i.e., Pr = 1 , which is only a special case of the MRT model.
It can be found that discrete Boltzmann equation is in a simple form and its algorithm is easy to code. In contrast, the NS equations depend upon both the first-order and second-order partial derivatives of velocities ( u x , u y ) with respect to x or y , which are nonlinear terms relatively difficult to be treated with [40] . Moreover, it often needs to solve the Poisson equation based on global data transfer in NS method, while all spatio-temporal information communication is local in DBM that is suitable for massively parallel computing. In addition, the DBM provides an efficient tool to study detailed nonequilibrium effects and/or rarefied effects of gas flows beyond NS equations by capturing the departures of kinetic moments from their equilibrium counterparts [40,43] . Finally, it is easy to have a proper kinetic boundary condition for DBM to describe the velocity slip and the flow characteristics in the Knudsen layer that cannot be well described by traditional hydrodynamic models [43] .

Validation and verification
For validation and verification purposes, four benchmark tests are performed. (i) The chemical reaction in a free falling box is simulated to verify the effects of external force and chemical reaction. (ii) The simulation of a detonation wave is carried out to demonstrate the DBM in the case with violent chemical heat release. Additionally, we assess the spatial and temporal convergence of the numerical results. (iii) To verify the DBM for adjustable specific heat ratios and Prandtl numbers, we simulate Couette flow. Moreover, it is demonstrated that the nonequilibrium information provided by the DBM coincides with its analytical solution. (iv) Finally, a typical two-dimensional benchmark, shock reflection, is simulated successfully. Besides, it is demonstrated in the first two tests that the discrete velocity model D2V16 has higher efficiency and better robustness than D2V24 [41] . Note that the second order Runge-Kutta scheme is adopted for the time derivative, while the second order nonoscillatory and nonfree-parameter dissipation difference scheme [44] is employed for the space derivative in Eq. (1) . It is preferable to set t ≤ 1 / Max ( ˆ S i ) due to the explicit scheme for the time derivative, where Max ( ˆ S i ) denotes the maximum among ˆ S i . The relation between the time step t and space step x = y should satisfy convergence conditions. Additionally, variables and parameters used in this paper are expressed in nondimensional forms, i.e., the widely accepted LB units [45,46] .

Detonation wave
In order to test the present DBM under the condition with violet chemical heat release, we target the detonation wave. The ini-tial configuration is (ρ, T , u x , u y , λ) L = (1 . 38837 , 1 . 57856 , 0 . 57735 , 0 , 1) (ρ, T , u x , u y , λ) R = (1 , 1 , 0 , 0 , 0) (17) where the suffix L indexes the left part, 0 ≤ x ≤ 0.05, and R the right part 0.05 < x ≤ 1, see Fig. 3 . The inflow or outflow condition is adopted in the x direction, the period condition is employed in the y direction. The parameters are I = 3 , Q = 1 , (1 . 7 , 3 . 7 , 3 . 3) , t = 10 −5 , x = y = 10 −4 , and N x × N y = 10 , 0 0 0 × 1 . The collision parameters are ˆ S i = 10 5 except ˆ S μ (i.e., ˆ S 5 , ˆ S 6 , ˆ S 7 ) = 2 × 10 4 . The detonation wave travels from left to right with speed v s . The chemical reactant is in front of the detonation wave with λ = 0 , and it changes into the product after the wave with λ = 1 . Fig. 4 illustrates the propagation of pressure at time instants, ferences are (0.05%, 0.04%, 0.08%, 0.02%), respectively. Obviously, the numerical and analytical results coincide well in Fig. 5 . The tiny differences between them are due to the fact that the ZND theory ignores the viscosity and heat conduction, and the von Neumann peak is assumed as a strong discontinuity which is not a truth. The DBM considers the viscosity, heat conduction as well as other nonequilibrium effects. Note that, with the decrease of collision parameters, the nonequilibrium effects are enhanced, and the differences between the DBM and analytical solutions become large [41] .
To compare the numerical robustness of D2V16 and D2V24 [41] , the aforementioned detonation wave is simulated by us-  Obviously, D2V16 gives a smooth profile around the detonation front, while D2V24 gives an oscillating profile. This nonphysical oscillation is soon amplified and results in the stop of the simulation program. Moreover, further tests demonstrate that D2V16 is capable of simulating the detonation wave for Mach number Ma > 100. However, it is difficult and even impossible to use D2V24 to simulate such high-Mach systems.
Next, let us assess the spatial and temporal convergence of the DBM results. The spatial convergence is proved considering several values of the space step, x = y = 5 × 10 −6 , 1 × 10 −5 , 2 × 10 −5 , 4 × 10 −5 , 8 × 10 −5 , 1 . 6 × 10 −4 , with fixed time step t = 1 × 10 −6 . The relative difference of the minimum value of ˆ f ne 5 around the detonation wave is chosen as the numerical error. Fig. 7 (a)   When the field reaches steady, the temperature is different for various γ or Pr , see Fig. 9 . The space step is x = y = 10 −3 , the time step t = 5 × 10 − 5 , and the parameters ( v a , v b , η a ) = (1.1, 1.7,   2.3). Periodic boundary conditions are employed for the left and right boundaries, and the nonequilibrium extrapolation method is applied to the top and bottom boundaries. The sketch of the initial configuration for Couette flow is shown in Fig. 8 . Fig. 9 illustrates the temperature T versus y when the Couette flow reaches equilibrium. Fig. 9 (a) shows the cases with γ = 1 . 3 , 1.5, 1.8, and fixed Pr = 1 . 0 ; Fig. 9 (b) shows the cases with Pr = 0 . 5 , 1.0, 2.0, and fixed γ = 1 . 5 . The collision parameter ˆ S μ is 2 × 10 3 for Pr = 0 . 5 , 1 × 10 3 for Pr = 1 . 0 , and 5 × 10 2 for Pr = 2 . 0 , the other collision parameters ˆ S i are 1 × 10 3 . The symbols represent DBM results, the lines denote the corresponding analytical solutions, Clearly, the numerical and analytical results coincide well with each other. Hence, the DBM has the capability of capturing the flow field in the dynamic process of the Couette Flow. In panel (b), the lines stand for the analytical solutions It can be found that the DBM results are in good agreement with the analytical values. That is to say, the DBM could describe the nonequilibrium behaviours accurately.

Shock reflection
For the purpose of verifying the model for two dimensional systems, we use a typical benchmark: regular shock reflection. The computational domain is a rectangle. The reflecting surface is imposed on the bottom, the supersonic outflow is adopted for the right boundary, and the Dirichlet conditions are utilized on the top and left boundaries, i.e., The interesting readers refer to Ref. [41] for more details of the initial configuration. The parameters are N x × N y = 300 × 100 , x = y = 10 −3 , t = 5 × 10 −6 , I = 2 , ( v a , v b , η a ) = (1.7, 2.9, 3.0). The collision parameters are ˆ S μ = 1 . 8 × 10 5 , and 2 × 10 5 for the others. Fig. 11 exhibits the density contour of the steady regular shock reflection. Theoretically, the angle between the incident shock wave and the wall is φ = π / 6 while the DBM gives the angle φ = ArcTan (0 . 1 / 0 . 173) . The relative difference between them is only 0.1%, which is satisfying.

Conclusions
We present an MRT DBM for compressible flows, taking both chemical reaction and external force into account. The specific heat ratio as well as the Prandtl number are flexible. This model recovers the reactive NS equations in the hydrodynamic limit.
Meanwhile, thermodynamic nonequilibrium effects are dynamically taken into account through considering the departures of kinetic moments from their equilibrium counterparts. In fact, the nonequilibrium effects together with their relaxation parameters play a crucial role in fluid systems.
Compared with a previous MRT DBM where 24 discrete velocities (and kinetic moments) are employed to couple the chemical reaction with fluid flows [41] , our model requires only 16 discrete velocities (and kinetic moments) and thus less computing efforts. Compared to another MRT DBM with the incorporation of only a conventional force term [37] , our model introduces a new form for both force and reaction terms, which are physically more general. In this paper, we also demonstrate that the present model provides high computational efficiency, physical fidelity, and numerical robustness iy , Then the first nine elements of ˆ F and ˆ R are obtained, i.e., ˆ ˆ F 6 = ρu x a y + ρu y a x , ˆ F 7 = 2 ρu y a y , ˆ F 8 = 2 ρu x u x a x + u y a y + ρa x u 2 + ρa x ( D + I + 2 ) T , ˆ F 9 = 2 ρu y u x a x + u y a y + ρa y u 2 + ρa y ( D + I + 2 ) T , Substituting the variables' expansion, (A.14) From Eqs. (A.3) to (A.5) , we get  ∂ξ ∂t 25) where j α = ρu α is the momentum in α direction, and ξ = (D + I) ρT + ρu 2 is twice the total energy, with