Three-dimensional multiple-relaxation-time discrete Boltzmann model of compressible reactive flows with nonequilibrium effects

Based on the kinetic theory, a three-dimensional multiple-relaxation-time discrete Boltzmann model (DBM) is proposed for nonequilibrium compressible reactive flows where both the Prandtl number and specific heat ratio are freely adjustable. There are 30 kinetic moments of the discrete distribution functions, and an efficient three-dimensional thirty-velocity model is utilized. Through the Chapman–Enskog analysis, the reactive Navier–Stokes equations can be recovered from the DBM. Unlike existing lattice Boltzmann models for reactive flows, the hydrodynamic and thermodynamic fields are fully coupled in the DBM to simulate combustion in subsonic, supersonic, and potentially hypersonic flows. In addition, both hydrodynamic and thermodynamic nonequilibrium effects can be obtained and quantified handily in the evolution of the discrete Boltzmann equation. Several well-known benchmarks are adopted to validate the model, including chemical reactions in the free falling process, thermal Couette flow, one-dimensional steady or unsteady detonation, and a three-dimensional spherical explosion in an enclosed cube. It is shown that the proposed DBM has the capability to simulate both subsonic and supersonic fluid flows with or without chemical reactions. © 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0047480


I. INTRODUCTION
Reactive flows encompassing a wide variety of nonlinear, unsteady, and nonequilibrium processes are common in nature and industry. 1 In fact, more than four-fifths of mankind's utilized energy is generated from the exothermic reactive flows. 2 In the past few decades, considerable research has been devoted to reactive flows, including supersonic and hypersonic flows related to the supersonic aircraft, rocket engine, detonation engine, supersonic combustion ramjet, etc. However, there are still many problems that have not been solved due to the complex physical and chemical processes involved, such as high compressibility, strong flow discontinuity, and combustion instability. In addition, these processes are usually accompanied by the hydrodynamic and thermodynamic nonequilibrium effects as well as the complex interplay between the chemical reaction and fluid flow. Moreover, these processes cover a wide range of spatial and temporal scales. [3][4][5][6] In order to study the reactive flow in detail, experimental techniques have been widely used. 7,8 However, experiments for supersonic and hypersonic reactive flows are difficult and expensive to conduct; and the measurements are usually global quantities only. On the other hand, numerical simulations provide a convenient tool for relevant research and become more and more reliable and cost-effective.
The supersonic reactive flow modeling approaches can generally be classified into three groups. The first group is the conventional simulation methods based on the continuum assumption, such as direct numerical simulation, 9 large eddy simulation, 10 and Reynolds-averaged Navier-Stokes (NS). 11 These conventional methods are good at capturing the main hydrodynamic characteristics. However, they lack the ability to describe thermodynamic nonequilibrium effects under highly nonequilibrated conditions, such as the shock or detonation front. 12,13 The second group is the microscopic methods, like molecular dynamics. 14 The molecular dynamics gets rid of the local equilibrium assumption and thus can be used to study the detailed hydrodynamic and thermodynamic properties. [15][16][17] On the other hand, the computational cost of molecular dynamics is usually prohibitive and the computational domain is always limited. The third group is the kinetic method 14 based on the Boltzmann equation, which removes the limit of the local equilibrium assumption. A kinetic method, the lattice Boltzmann method (LBM), simulates the evolution of probability distribution functions in a discretized phase-space, which can be related to macroscopic quantities by moment relationships. The LBM has been successfully applied to simulate a variety of complex flows, 18 including reactive flows. [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] Merits of the LBM include the algorithm simplicity and locality, which lead to excellent performance on parallel clusters.
The key to dealing with reactive flows in the LBM is how to describe the energy and species mass fractions by LB formulation and couple them consistently. Current LBM formulations for reactive flows can be divided into two categories. The first category is the hybrid method, in which the flow simulation is accomplished through the weakly compressible LBM solver, while the transport equations for energy and species are solved by a finite difference scheme. In 2000, Filippova and Häenel proposed an LBM for reacting flows at a low Mach number with variable density. 20,21 Later, Hosseini et al. modified the LBM solver and successfully simulated premixed and non-premixed flames with a detailed thermo-chemical model. 22,23 In 2020, Shu et al. developed a simplified sphere function-based gas-kinetic flux solver for compressible viscous reacting flows. 37 This model applies the finite volume method to discretize the multi-component NS equations and computes numerical flux at the cell interface by using local solution of the Boltzmann equation. Although these models can handle the variable density, they lose the simplicity and the parallel efficiency of the pure LBM scheme. 24 The other category is the pure LB formulation. In 1997, with the assumption of irreversible infinite fast chemistry reactions, Succi et al. adopted the conserved scalar approach to describe the temperature and concentration fields. 19 Similarly, in 2002, Yamamoto et al. assumed that the temperature field does not affect the flow field and simulated diffusion flames with a double-distribution-function LBM where the flow, temperature, and species fields are represented by two sets of distribution functions. 25,26 Later, Lee et al. 27 and Chen et al. 28,29 also utilized the double-distribution-function LBM to solve the low Mach number flows. Instead of using two sets of distribution functions, in 2012, Prasinari et al. extended a consistent LBM 38 by introducing correction terms, recovering the third-and fourth-order moments, and describing the temperature field. 30 In conclusion, these pure models are all limited to low Mach number flows. However, for supersonic and hypersonic reactive flows, high compressibility is an important factor.
Briefly, the aforementioned LBMs mainly focus on low Mach number reactive flows, which cannot make the best of the kinetic theory. To establish an LBM for reactive flows possessing more kinetic information beyond the NS equations, the LBM should be compressible, thermal, and, at the same time, coupling the chemical reaction naturally. In addition, the model should take hydrodynamic and thermodynamic nonequilibrium effects into consideration.
In recent years, the discrete Boltzmann method (DBM) has been constructed to model and simulate nonequilibrium systems, with various velocity and time-space discretization schemes. [31][32][33][34][35]39 The DBM and LBM can be viewed as two distinctive classes of discrete numerical methods based on the Boltzmann equation. Despite sharing the same origin and some similarities, the objectives, numerical implementation, and capabilities of the LBM and DBM are different. Basically, the LBM aims to recover the macroscopic behaviors of the flow system and acts as an alternative method to the continuum-based Navier-Stokes solver. On the other hand, the DBM is a discrete Boltzmann equation solver that aims to keep some kinetic features that go beyond the macroscopic behaviors. Numerically, discretization in space, time, and particle velocities is inter-dependent in the LBM algorithm. In the DBM, however, discretization in space, time, and particle velocities is decoupled, which allows a variety of numerical methods to be applied. As a result, the LBM requires different models for different types of flows (e.g., incompressible, compressible, thermal, or reactive), while a uniformed DBM framework can simulate all types of flows. Moreover, the DBM can capture both hydrodynamic and thermodynamic nonequilibrium effects explicitly.
To some extent, the mathematical framework of the DBM is similar to the rational extended thermodynamics (RET), which is also capable of explaining nonequilibrium phenomena. The RET proposed by Ruggeri and Sugiyama 40,41 is motivated by the Boltzmann equation. Different from the kinetic theory, the RET focuses on the relations of the moments and forms a hierarchy of the balance laws. To achieve the closure of the moments, the RET truncates the hierarchy by adding restrictions from the universal principle. In this respect, the RET theory still belongs to the continuum approach but is applicable to the nonequilibrium state. On the other hand, the DBM is a kinetic approach and possesses kinetic features beyond the macroscopic equations.
Since 2013, several DBMs for reactive flows have been proposed. [31][32][33][34][35]39 In 2013, Yan et al. proposed the first DBM for combustion. 31 Very recently, Lin et al. proposed a two-dimensional model for detonations and investigated the main features of the hydrodynamic and thermodynamic nonequilibrium effects. 36 However, these formulations are two-dimensional but realistic reactive flows are a three-dimensional (3D) phenomenon and some important patterns and structures cannot be obtained from twodimensional models. In 2017, Gan et al. proposed a 3D DBM for compressible flows without reaction based on the single-relaxationtime Boltzmann equation, fixing the Prandtl number Pr = 1. 42 In this work, we extend the model to 3D reactive flows and present a multiple-relaxation-time (MRT) method to make the Prandtl number adjustable. Moreover, the hydrodynamic and thermodynamic nonequilibrium effects can be investigated to study the nonequilibrium behaviors.
The rest of this paper is organized as follows: In Sec. II, we formulate a 3D MRT DBM for reactive flows based on a three-dimensional thirty-velocity (D3V30) model. Through the Chapman-Enskog multiscale analysis, the reactive NS equations are recovered and the nonequilibrium effects are derived in Sec. III.

ARTICLE scitation.org/journal/adv
Section IV presents numerical tests, and Sec. V provides a summary and conclusions.

II. DISCRETE BOLTZMANN METHOD
In the Bhatnagar-Gross-Krook (BGK) model, 43,44 a single relaxation time τ determines all discrete distribution functions approaching their equilibria at the same speed, and the Prandtl number is fixed at Pr = 1. To overcome this shortage, we construct an MRT DBM to make the Prandtl number adjustable, which takes the form ∂f ∂t where f = ( f 1 , f 2 , . . . , f N ) ⊺ and f eq = ( f ⊺ denote kinetic moments of discrete distribution function and their equilibrium counterparts, respectively. Here, the subscript N = 30 is the total number of discrete velocities. In fact, the terms in moment space are transformed from those in velocity space by a transformation matrix M, which is a square matrix with N × N elements in terms of discrete velocities. The elements of the transformation matrix are given in the Appendix. Specifically, (2) Similarly, R = (R 1 , R 2 , . . . , RN) ⊺ and F = (F 1 , F 2 , . . . , FN) ⊺ stand for the reaction and force terms in velocity space, respectively.R = MR andF = MF stand for discrete reaction and force terms in moment space, respectively. S = diag (S 1 , S 2 , . . . , SN) ⊺ is a diagonal matrix that consists of relaxation coefficients Si determining the speed off i approachingf eq i . The additional termÂ = MA = (0, . . . , 0,Â 12 ,Â 13 ,Â 14 , 0, . . . , 0) ⊺ is used to modify the collision operator so that the discrete Boltzmann equations could recover the correct reactive NS equations via the Chapman-Enskog analysis in terms of where ρ, T, p(= ρT), and uα are the density, temperature, pressure, and velocity, respectively. Here, D = 3 stands for the number of dimensions and I represents the extra degrees of freedom. The discrete equilibrium distribution function satisfies the following moment relations: Here, α, β, and χ denote the direction that can be x, y, or z. Based on physical considerations and following the model proposed by Bourgat et al., 45 we introduce a single additional variable I to represent non-translational degrees of freedom and utilize a free parameter η i to describe molecular rotation and/or internal vibration energies. The specific heat ratio is defined as γ = (D + I + 2)/(D + I).
According to Eqs. (3) and (7)- (13), the discrete equilibrium distribution functions can be expressed as To ensure the matrix M invertible and get better stability, we improve the discrete velocity model D3V30 proposed by Gan et al. 42 Instead of using one parameter, we adopt four parameters va, v b , vc, and v d to determine the magnitude of four sets of discrete velocities, respectively, in terms of More forms of the antisymmetric part can be found in Ref. 42. In Eqs. (7)- (9), we see f eq i can be replaced by f i according to the conservation laws. However, in Eqs. (10)-(13), replacing f eq i with f i may lead to the imbalance between left and right sides. The differences are actually departures of high order kinetic moments from their equilibria and can be utilized to investigate the nonequilibrium effects. We define the following nonequilibrium quantities: Moreover, we introduce the symbols linked with the viscous stress tensor, related to the nonorganized energy fluxes, and associated with the flux of nonorganized energy flux. Physically, iα is defined as the translational energy in the α direction, and 1 2 ∑ i f eq i v 2 iα is its equilibrium counterpart. The nonequilibrium part 1 2 Δv α v α is the nonorganized energy in the α direction, which corresponds to the molecular individualism on top of the collective motion. Clearly, the DBM can provide the above nonequilibrium information beyond traditional NS equations.
The force and reaction terms are the variation rates of the distribution functions resulting from the external force and chemical reaction, respectively. The two terms are derived based on the following assumptions: (1) Over a small time interval, the change of the distribution function due to the external force and the chemical reaction can be treated as the change of the equilibrium distribution function, which is the leading part of the distribution function when the system is not too far from the equilibrium state. 46 (2) The effect of the external force is to change the hydrodynamic velocity u with the acceleration a. In a small time interval, the velocity becomes u † = u + aτ. (3) The temporal scale of the chemical reaction is much smaller than that of fluid flow, and the chemical reaction leads to the change of energy with the varying rate, where Q indicates the chemical heat release per unit mass of fuel and λ indicates the mass fraction of chemical product. From Eq. (9), we obtain the temperature after the chemical reaction T ♢ = T + τT ′ with the varying rate of temperature Now, we can derive the force term, and the reaction term, The discrete forms Fi and Ri satisfy the relation with (23) and (24) results inF = MF andR = MR, respectively, and the elements of the force and reaction terms are given in the Appendix.
To describe the dynamics of detonations and compare with the previous study on the instability of detonations, we consider a simple model of a chemical reaction between two perfect gases, assuming irreversible, one-step Arrhenius kinetics, where K is the rate constant and Ea is the activation energy. The DBM employs larger velocity stencil and introduces higher-order moments, by which the hydrodynamic and thermodynamic fields are fully coupled. It is also worth mentioning that we utilize the matrix inversion method here instead of the commonly adopted polynomial approach where equilibrium distribution functions are expanded upon the terms of the macroscopic quantities due to the following merits. The number of discrete velocities of the DBM equals exactly the number of required kinetic moments of equilibrium discrete distribution functions, while the polynomial method always requires more discrete velocities because the discrete velocity sets of the latter should have enough isotropy to recover the hydrodynamic equations correctly. 47,48 Consequently, the presented method is more efficient. Additionally, in the DBM, each kinetic moment computed by summation of the discrete equilibrium distribution functions is exactly equal to the one calculated by integral of the Maxwellian function. Furthermore, the discrete Boltzmann equations are expressed in a uniform form, and thus, the DBM is easier to code. The DBM code is parallelized on the UK National Supercomputing Service ARCHER and runs efficiently.

III. THE CHAPMAN-ENSKOG ANALYSIS
In this section, we show the main procedure of the Chapman-Enskog analysis. With respect to an expansion parameter ε, which is a quantity of the order of the Knudsen number, 49 we introduce the following expansions: where the part of distribution function f 1, 2, . . .). By substituting Eqs. (28)-(33) into Eq. (1), we can obtain the following equations in consecutive order of the parameter ε: where j α = ρu α is the momentum and ξ = (D + I)ρT + (j 2 x + j 2 y + j 2 z )/ρ is twice the total energy. Combining Eqs. (37)- (39) and the sixth to the fourteenth relations of Eq. (35), we can derive the following quantities: (1)

S7f
(1) S 8f (1) S 9f (1) S 10f (1) S 11f (1) S 12f S 13f (1) S 14f (1) The additional terms in Eqs. (46)- (48) are determined in Eqs. (4)- (6), which are used to modify the collision terms so that the NS equations can be correctly recovered. In the case that the system is isothermal, the additional terms can be eliminated. The above quantitiesf i are the exact solutions of nonequilibrium quantities on the level of the first order accuracy. Similar to the ARTICLE scitation.org/journal/adv above derivation of t 1 = εt scale, the macroscopic equations on the t 2 = ε 2 t time scale are derived by the first five relations in Eq. (36), With the aid of Eqs. (37)-(48), the final equations can be written as in terms of where Sxx = S 6 , Syy = S7, Szz = S 8 , Sxy = S 9 , Sxz = S 10 , and Syz = S 11 . Here, we explain the reason for the additional term. Mathematically, we can adjust the relaxation coefficients independently. From the point view of physics, the coefficients are not completely independent for the system with isotropy constraints. 50 From the derivation, we can see that S 1 , S 2 , S 3 , S 4 , and S5 have no influence on the equation due to the conservation laws. However, in order to recover the NS equations in the continuity limit, the terms referring to the stress tensor and the energy flux should be coupled, respectively, which leads to related to viscosity, and related to heat conductivity, and then Eqs. (54)-(56) reduce to where and the dynamic viscosity μ, the bulk viscosity μ B , and the thermal conductivity κ are defined as and κ = ( respectively. Furthermore, the Prandtl number, is adjustable in the model. In conclusion, we modify the collision operator by the additional term so that the NS equations can be recovered correctly. Specifically, we can find that the nonequilibrium quantity is just the viscous stress tensor, and is related to the thermal conductivity (see Ref. 51). As we can see, the diagonal matrix S controls the speed off i approaching its equilibrium counterpartf eq i and the elements of the matrix Si are related to the thermodynamic nonequilibrium manifestations. S 1 , S 2 -S 4 , and S5 are related to mass, momentum, and energy conservation, respectively. Actually, the values of the above coefficients have no effect on the flow system due to the conservation laws. S 6 -S 11 are related to dynamic viscosity and S 12 -S 14 are linked with heat conductivity in the NS equations. Finally, S 15 -S 24 are associated with nonorganized energy fluxes and S 25 -S 30 correspond to the flux of nonorganized energy flux.

IV. NUMERICAL TESTS
To validate the proposed model and showcase its performance, we show simulation results of some classical benchmarks. First, we simulate the free falling process with a chemical reaction to investigate the model in the flow field coupled with a chemical reaction and external force. Then, we simulate the thermal Couette flow to verify that both the specific heat ratio and Prandtl number are adjustable. Furthermore, we simulate a 1D detonation to show that the model is capable of describing supersonic flows and measuring the nonequilibrium effects. Finally, a spherical explosion is modeled to demonstrate the ability of the proposed model to deal with a 3D configuration. In this work, we utilize the third-order Runge-Kutta scheme for the time derivative in Eq. (1) and the second-order nonoscillatory and nonfree-parameter dissipation difference scheme for the space derivative. 52

A. Chemical reaction in the free falling process
The exothermic chemical reaction in a free falling box is simulated as the first numerical test. At the beginning, the box is filled with the chemical reactant and the initial physical quantities are density ρ 0 = 1, velocity u 0 = 0, and temperature T 0 = 1.0. Since the field is uniform, we use only one mesh to simulate the process, and thus, the computational domain is Nx × Ny × Nz = 1 × 1 × 1, where the periodic boundary is used in each direction. We choose Δx = Δy = Δz = 10 −3 as the spatial step and Δt = 10 −4 as the temporal step. Reaction parameters are K = 5 × 10 2 and Ea = 8, and relaxation coefficients are Si = 10 4 . Figure 1(a) plots the vertical velocity uz vs various acceleration a = azez with fixed chemical reaction heat release Q = 10 at time t = 0.5. The theoretical vertical velocity in the external force field is uz = azt. Figure 1(b) illustrates the temperature T vs chemical heat release with fixed acceleration az = 10. The theoretical solution for temperature after the chemical reaction is T = T 0 + (γ − 1)Q. From the above pictures, we can find the simulation results match well with the analytical ones, which shows good performance of the model to describe the effects of the external force and chemical reaction.

B. Thermal Couette flow
To validate that the model can be used with the adjustable specific heat ratio and Prandtl number, the thermal Couette flow is simulated here. The initial configuration is a viscous fluid flow between two infinite parallel flat plates, and the physical quantity of the flow is ρ 0 = 1, u 0 = 0, and T 0 = 1. The plate below the flow is

ARTICLE
scitation.org/journal/adv stationary with temperature TL = 1, and the top plate moves along the horizontal direction with a constant speed ux = 1 and temperature TU = 1.1. The distance between the plates is H = 0.1. We carry out the simulation with Δx = Δy = Δz = 10 −3 , Δt = 10 −4 , va = v b = vc = v d = 1.6, and η 0 = 2.1, and the grid mesh is Nx = Ny = Nz = 1 × 1 × 100. Periodic boundary conditions are utilized in the x and y directions, and the nonequilibrium extrapolation method is applied in the z direction. The analytical solution to the temperature profile along the z direction is Figure 2(a) shows the comparisons of the DBM simulation results and exact solutions of the temperature along the z direction with different specific heat ratios γ = 3/2, 7/5, and 4/3 and fixed Prandtl number Pr = 1.0, where collision parameters Si = 2 × 10 3 . It is evident that the simulation results agree excellently well with the analytical ones. Figure 2(b) plots the temperature distribution with various Prandtl numbers Pr = 0.5, 1.0, and 2.0, where collision parameters Si = 2 × 10 3 and Sμ = 4 × 10 3 , 2 × 10 3 , and 1 × 10 3 , respectively. The specific heat ratio is 7/5 for the three cases. The simulation results also match well with the exact solutions for various Prandtl numbers.

C. One-dimensional detonation
In order to make the calculation simple and explicit, all the quantities have been made dimensionless by reference to the unburned state ahead of the detonation front. The reference length scale x ref is chosen such that the half of 1D Zel'dovich, von Neumann, and Döring (ZND) induction length is unity, Adopting the notations in Ref. 53, the symbol ∼ stands for dimensional quantities and subscript 0 denotes the unburnt state. Previous research indicates that the overdrive factor f (the square of the ratio of imposed detonation front velocity to the Chapman-Jouguet velocity), the specific heat ratio γ, the heat release Q, and the activation energy Ea are all related to the stability of the detonations. [54][55][56][57] In this part, we first simulate an unstable 1D detonation with the following parameters to verify the effectiveness of the model: f = 1.6, γ = 1.2, Q = 50, and Ea = 50. Given the above dimensionless values, the pre-exponential parameter equals K = 230.75. The initial condition is given by the theory developed independently by ZND. The physical domain is set to be 0 ≤ x ≤ 800 with an inflow at the left boundary and an outflow at the right boundary. The periodic boundary condition is imposed on the upper and lower boundaries. The initial location of the ZND detonation front is set at x = 8.
We carry out the simulation with va = 1.0, v b = vc = 4.5, v d = 0.1, η 0 = 2.7, and Si = 5 × 10 3 . Different grid resolution δn (denotes the number of grid points per half-reaction length) is employed to find an appropriate resolution for computation effectiveness.
As common practice, we use the peak pressure, which is the maximum pressure at the precursor explosion, to validate the performance of the numerical schemes. 54,[56][57][58] Figure 3 shows temporal histories of the peak pressure under different resolutions. As one can see, all the results oscillate periodically for t > 20. The solution with the lower resolution δ 50 yields a poor result. With the increase in the resolution, the higher resolutions δ75 and δ 100 behave similarly. To further validate the accuracy of the DBM, we compare the maximum, minimum, and the period of the oscillations of the DBM (Pmax, P min , T period ) = (99.05, 57.5, 7.25) with the predicted peak pressure (Pmax = 98.6) in the literature 54   In addition to the unstable detonation, we simulate a stable detonation with the following parameters: f = 1.0, γ = 1.4, Q = 1.0, and Ea = 8.0. The pre-exponential parameter equals K = 122.77. The physical domain is set to be 0 ≤ x ≤ 500, and the initial location of the ZND detonation front is set at x = 5. For stable detonations, δ 100 is used. The parameters adopted are va = v b = vc = v d = 1.8, η 0 = 2.0, and Si = 2 × 10 2 . respectively. The results of the DBM have a satisfying coincidence with the ZND theory. Figure 6 displays the kinetic moments and nonequilibrium quantity around the detonation wave at time t = 160 in order to demonstrate the capability of the DBM of capturing the nonequilibrium effects. This result is physically reasonable. At the front shock, the chemical energy is continuously released and due to the rapid change of physical quantities, the system departs from its thermodynamic equilibrium state.
To demonstrate the evolution process of the spherical explosion, we choose several typical snapshots of the pressure isosurfaces at various times in Fig. 7. Since the computational system is symmetrical, only half of the system is depicted for clear illustration. Figure 7(a) shows the initial configuration, and (b) and (c) display the evolution process at time t = 0.0025 and 0.005, respectively. First, the spherical shock wave expands in the enclosed box and contacts with the wall when t = 0.005. Afterward, the shock wave is reflected and propagates inwards with the increase in time. As we can see, the shocks are well captured in the DBM. Figure 8 illustrates the temporal histories of the mass and energies in the explosion process. The total reactant mass decreases and the total product mass increases gradually due to the chemical reaction in the initial stage (0 ≤ t ≤ 0.011 54). After the chemical reaction ends (t > 0.011 54), the reactant density equals zero. The total mass remains constant in the whole evolution. Figure 8(b) shows in the evolution, the chemical energy decreases gradually and is transformed into the internal and kinetic energies, while the sum of all energies remains unchanged. Figure 8 indicates that the mass and energy of the system are conserved. It is shown that the proposed DBM has a satisfying performance for simulations of 3D reacting flows.

V. CONCLUSION
In this paper, a 3D MRT DBM is presented for reactive flows where both the Prandtl number and specific heat ratio are freely adjustable. There are 30 kinetic moments of discrete distribution functions, and an efficient discrete velocity model, D3V30, is utilized. The NS equations can be recovered at macroscales from the DBM through the Chapman-Enskog analysis. Unlike existing LBMs for reactive flows, the hydrodynamic and thermodynamic fields are fully coupled in the DBM to simulate combustion in subsonic, supersonic, and potentially hypersonic flows. In addition, the nonequilibrium effects of the system can be probed and quantified. In this model, the reaction and force terms added into the discrete Boltzmann equation describe the change rates of discrete functions due to the chemical heat release and external force, respectively.
The DBM has been validated for several typical applications. The capability of the DBM to simulate the flow field fully coupled with a chemical reaction and external force is verified using the free falling reacting box. The case of the thermal Couette flow demonstrates that both the specific heat ratio and Prandtl number are adjustable. Furthermore, the DBM is shown to accurately simulate supersonic flows and quantify the nonequilibrium effects in 1D detonation. Finally, the main features of a spherical explosion in an enclosed box are successfully captured to showcase the ability of the DBM to deal with a 3D configuration. In conclusion, a new 3D DBM for reactive flows featuring full coupling of flow and combustion fields has been developed and successfully validated. This opens up possibilities for exploring a variety of reactive flows at subsonic, supersonic, and hypersonic speeds with both hydrodynamic and thermodynamic nonequilibria.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.