Mesoscopic simulation of nonequilibrium detonation with discrete Boltzmann method

Thanks to its mesoscopic nature, the recently developed discrete Boltzmann method (DBM) has the capability of providing deeper insight into nonequilibrium reactive flows accurately and efficiently. In this work, we employ the DBM to investigate the hydrodynamic and thermodynamic nonequilibrium (HTNE) effects around the detonation wave. The individual HTNE manifestations of the chemical reactant and product are probed, and the main features of their velocity distributions are analyzed. Both global and local HTNE effects of the chemical reactant and product increase approximately as a power of the chemical heat release that promotes the chemical reaction rate and sharpens the detonation front. With increasing relaxation time, the global HTNE effects of the chemical reactant and product are enhanced by power laws, while their local HTNE effects show changing trends. The physical gradients are smoothed and the nonequilibrium area is enlarged as the relaxation time increases. Finally, to estimate the relative height of detonation peak, we define the peak height as H(q ) = (q max − q s ) / (q v on − q s ) , where q max is the maximum of q around a detonation wave, q s is the CJ solution and q von is the ZND solution at the von-Neumann-peak. With increasing relaxation time, the peak height decreases, because the nonequilibrium effects attenuate and widen the detonation wave. The peak height is an exponential function of the


Introduction
Detonation is a type of violent combustion that propagates through an exothermic supersonic wave [1,2] . In practice, detonation is extremely complex as it encompasses various interactive physicochemical phenomena, covers a wide range of spatiotemporal scales, and incorporates hydrodynamic, thermodynamic and chemical nonequilibrium effects [3][4][5][6] . Specifically, a flow is in hydrodynamic equilibrium if the net force on each of its individual parts is zero, otherwise, it is hydrodynamic nonequilibrium (HNE). A system is thermodynamic equilibrium if the velocity distribution function is the Maxwellian distribution function, otherwise, it is thermodynamic nonequilibrium (TNE). In fact, both hydrodynamic and thermodynamic nonequilibrium (HTNE) effects often exert significant influences on fluid systems, especially on detonation. The nonequilibrium detonation phenomenon has long been studied [7,8] . Nichols examined the influence of nonequilibrium diffusional flow upon detonation velocities in composite ex-E-mail addresses: chuandonglin@163.com (C. Lin), K.Luo@ucl.ac.uk (K.H. Luo). plosives [7] . Romick et al. investigated the dynamics of a onedimensional detonation using a one-step irreversible Arrhenius kinetic model, which considered partly nonequilibrium features of mass, momentum and energy diffusion [8] . Previous works, however, were mainly focused on macroscopic manifestations by using traditional hydrodynamic methods, usually ignoring more abundant, complex and essential thermodynamic nonequilibrium effects caused by interactions at microscopic scales.
To predict fluid behaviors accurately, it is necessary to take HTNE into consideration. A possible method is to use molecular dynamics (MD) [9][10][11] or direct simulation Monte Carlo (DSMC) [12,13] . Kawakatsu and Ueda employed MD to study the energy relaxation processes around a detonation wave [9] . Liu et al. presented an MD simulation of shock waves, and probed the ion velocity distributions and nonequilibrium characteristics near the shock front [11] . With the DSMC method, Bruno et al. showed that the nonequilibrium feature of the velocity distribution function at the shock front changes significantly the kinetics of chemical reactions and the global profile of the reaction zone [12] . But both MD and DSMC encounter the issue that the spatio-temporal scales amenable are rather limited due to the high computation cost. To address this problem, a promising way is to resort to a kinetic methodology based on the Boltzmann equation. Up to today, great efforts have been made and various kinetic models have been developed, such as the fast spectral method [14] , unified gas kinetic scheme [15] , discrete unified gas kinetic scheme [16] , discrete velocity method [17] , Lattice Boltzmann model (LBM) [18] , and discrete Boltzmann method (DBM) [19] , etc.
In fact, as a mesoscopic method, the DBM has the capability of modeling and simulating nonequilibrium systems with relatively high efficiency [19][20][21][22][23][24][25] . The idea of investigating nonequilibrium behavior with the DBM was explored by Xu et al. [19] . Specifically, kinetic moments of the departure of the distribution function from its equilibrium counterpart were regarded as the nonequilibrium manifestations. Furthermore, Yan et al. proposed the first DBM for combustion and detonation in 2013 [20] . Later, Lin et al. presented a polar-coordinate DBM for explosion or implosion, and investigated the main features of the velocity distribution function by analyzing nonequilibrium manifestations [21] . Zhang et al. then demonstrated that nonequilibrium characteristics are related to the entropy increasing rate [26] . More recently, Zhang et al. constructed a discrete ellipsoidal statistical BGK model and gave a scheme to recover the actual velocity distribution function quantitatively [27] . The DBM results are consistent with the MD results [11,27] . In the following, we use a recently developed two-component DBM [23] to study the individual HTNE effects of chemical reactant and product around the detonation wave.

Discrete Bolzmann method
The DBM is a variant of the LBM that has been employed as an effective tool for complex systems [28][29][30][31][32][33][34][35][36][37][38][39] . The LBM has been widely adopted as a solver for partial differential equations, such as incompressible Navier-Stokes (NS) equations. Several sets of discrete distribution functions are used to describe physical quantities, such as the flow velocity, temperature, and species concentration. Consequently, the flow field and chemical reaction are not coupled naturally in traditional LBMs. In contrast, the hydrodynamic and thermodynamic quantities interact with each other in the DBM, which is consistent with the Boltzmann equation. All physical quantities are described by the same set of discrete distribution functions. In addition, the DBM provides more nonequilibrium information beyond the macroscopic governing equations or the traditional LBM solvers.
The DBM is based on the following equation, where the superscript σ denotes the fluid component, t the time, τ the relaxation time, f σ i the discrete distribution function. In Eq. (1) , the time derivative is solved in analytic form [40] , while the space derivatives are treated with the second order nonoscillatory and nonfree-parameter dissipation difference scheme [41] . The discrete velocities take the form, Figure 1 illustrates the sketch of the discrete velocities. It is noteworthy that various species have independent discrete velocities. We use the same lattice for different species in this work. The discrete velocities, spatial step, and temporal step are decoupled in the DBM, which is different from the standard LBM. This feature makes the DBM possess better numerical robustness than the standard LBM. The molar concentration n σ , density ρ σ , hydrodynamic velocity u σ , and temperature T σ are calculated with the following formula, where m σ = 1 represents the molar mass, D = 2 stands for the space dimension, I σ indicates extra degrees of freedom corresponding to molecular rotation and/or vibration. η σ i = η σ α for And the mixing internal energy is E = σ n σ ( D + I σ ) T / 2 . The discrete equilibrium distribution function, f isfies the following seven kinetic moments, n σ = i with I = δ αβ e α e β , ξ σ = (D + I σ + 2) T m σ + u · u , = ( δ βχ u α + δ αχ u β + δ αβ u χ ) e α e β e χ , and e α the unit vector in the α direction.
The above seven relations can be uniformly written as where the column matrix f eq = ( f σ eq , and M is a square matrix with 16 × 16 elements (see the appendix in Ref. [23] for details). Then, the discrete equilibrium distribution function is expressed as The reaction term R σ i is the change rate of the discrete distribution function due to the chemical reaction. The treatment of the chemical reaction is based on the following considerations. The chemical reaction is irreversible and exothermic. The electronic excitation, ionization and radiation are negligible. The time scale of the chemical reaction is much longer than that of the kinetic process but shorter than that of the hydrodynamic flow behavior. Therefore, during a relatively short period of the local chemical reaction, both mass density and hydrodynamic velocity of the mixture remain unaffected while the temperature changes with the chemical heat released. The thermodynamic nonequilibrium behaviors take place on the relaxation time scale, which is much shorter than that of the chemical reaction. The reaction term is calculated as below, where f bution function after a time step, t , during the chemical reaction [23] . The physical quantities n σ * and T * are computed as below, . (23) In this work, σ = A and B stand for the chemical reactant and product, respectively. The chemical reactant changes into the prod- The Arrhenius function is adopted to describe the reaction, where k ov is the reaction coefficient, n A the molar concentration of chemical reactant, E a the effective activation energy, R = 1 the universal gas constant. Note that variables and parameters used in this paper are expressed in nondimensional forms, i.e., the widely accepted LB units. The results are then expressed in the physical units based on the similarity. It has been proved that the DBM can recover the reactive NS equations [23,24] and capture the following nonequilibrium quantities [23,24] , Physically,  of species σ in the α direction, i.e., the departure of the energy flux of species σ in the α direction from its equilibrium counterpart. Clearly, the above nonequilibrium information provided by the DBM is beyond traditional NS models [23,24] .

Nonequilibrium detonation
Let us investigate the nonequilibrium detonation. First of all, we study the main HTNE characteristics of the chemical reactant and product, from which we obtain the fundamental features of the velocity distributions near the detonation front. Next, we investigate the impact of the chemical heat and relaxation time on the HTNE. For the sake of convenience and simplicity, only one-dimensional steady detonation is considered. Moreover, the inflow and periodic boundary conditions are adopted in the x and y directions, respectively. To ensure the resolution is high enough, we adopt the mesh N x × N y = 50 0 0 × 1 , the lattice width x = y = 2 × 10 −4 , and the time step t = 5 × 10 −6 .

Nonequilibrium characteristics of detonation
The first simulation is conducted with the relaxation time τ = 4 × 10 −5 , the chemical heat per unit mass Q = 2 , and the specific heat ratio γ = 1 . 4 . The initial field is, where the subscripts L and R designate 0 ≤ x < 0.1 and 0.1 ≤ x < 1.0, respectively. The left part is occupied by the chemical product, and the right part is filled with the chemical reactant. The quantities of the two parts satisfy the Rankine-Hugoniot relations. Other pa- ical detonation speed is v s = 2 . 516 , whose relative error is only 0.001% compared to the solution 2.51603. It is satisfying. In addition, there are a few differences between the numerical and analytical results near the detonation front. The ZND theory assumes a sharp discontinuity at the detonation wave and ignores all thermodynamic nonequilibrium effects [1] , while we take into consideration the viscosity, heat conduction and other nonequilibrium effects. Figure 2 (b) exhibits the mole fraction Y σ . As the detonation wave travels from left to right, the chemical reactant changes into the product, and the chemical heat converts to kinetic energy and internal energy. Obviously, Y A changes from 1 into 0, Y B changes from 0 into 1, and Y A + Y B = 1 . Figure 2 02486 at the guideline. It can be found that the total nonequilibrium effects of the fluid mixture are weak due to the opposite nonequilibrium effects of the chemical reactant and product near the von-Neumann-peak, while their individual nonequilibrium effects are relatively strong.
In fact, from above nonequilibrium manifestations, we can obtain the fundamental characteristics of the distribution functions f σ , see Fig. 3 . To be specific, the equilibrium distribution function takes the form [42] f M = n σ m σ 2 π kT m σ 2 π I σ kT where k = 1 is the Boltzmann constant. Mathematically, f M is a normal distribution function whose peak is located at u . At the pressure maximum, A 2 ,xx < 0 implies that, for some v y , the distribution function of the chemical reactant f A ( v x ) is "thinner" and "higher" than f M ( v x ). Meanwhile, A 3 , 1 ,x < 0 indicates that, for some v y , f A ( v x ) is asymmetric, i.e., the portion of f A for v x > u x is smaller than that for v x < u x . Similarly, B 2 ,xx > 0 means that, for some v y , the distribution function of the product f B ( v x ) is "fatter" and "lower" than f M ( v x ). B 3 , 1 ,x > 0 results from that, for some v y , the portion of f B ( v x ) for v x > u x is larger than that for v x < u x . Moreover, with the nonequilibrium quantities A 2 ,yy = 0 . 02046 , B 2 ,yy = 0 . 00197 , A 3 , 1 ,y = B 3 , 1 ,y = 0 at the guideline, we obtain that f A ( v y ) and f B ( v y ) are "fatter" and "lower" than f M , and both are symmetrical about the axis v y = u y , see Fig. 3 (b). Combining the sketches of f σ versus v x and v y , it is easy to obtain the main features of f σ in the velocity space ( v x , v y ), see Fig. 3 (c) and (d). Clearly, the peak of f B is shifted towards the burnt gas, whereas that of f A is leaning towards the fresh gas.
It is worth noting that the source of the nonequilibrium effect is the physical gradients and chemical reaction [23,43] . In turn, the nonequilibrium affects the physical gradients and chemical reaction. The physical gradients and chemical reaction interact with each other as well. Next, let us study the interplay between them.

Effects of chemical heat
It is clear that the departure of f σ from f σ eq is far with increasing nonequilibrium physical quantities. To gain a deeper understanding of nonequilibrium effects on detonation systems, we probe them with various chemical heat. Figure 4 exhibits physical quantities around the detonation wave with chemical heat Q = 64 , which is compared to the case with Q = 2 in Fig. 2 . We can find that the length of the reaction zone (or the nonequilibrium zone) is about 0.01 in Fig. 2 , while it is around 0.02 in Fig. 4 . That is to say, the nonequilibrium area becomes wider with increasing chemical heat. Meanwhile, the physical quantities ( p, T, ρ, u x ) increase faster and the physical gradients become steeper for increasing chemical heat. Physically, increasing chemical heat promotes the chemical reaction rate and results in a sharp increase of physical quantities and their gradients near the detonation front. The HTNE effects are strong if the chemical reaction is violent, and/or the physical quantities as well as their gradients are large [23] . Figure 5 (a) demonstrates the amplitude of σ 2 ,xx versus the chemical heat. Its amplitude is defined as the half distance between the peak and trough of σ 2 ,xx around the detonation wave.
note the chemical reactant and product, and lines represent the fit- The relationships between these nonequilibrium effects and the chemical heat roughly obey a power law. Mathematically, the nonorganized energy depends on the chemical heat, reaction rate, relaxation time, specific heat ratio, and physical field around the detonation wave [23] . Physically, the chemical heat released plays a direct role in the mixing temperature (an equilibrium physical quantity) and the nonorganized energy (a nonequilibrium physical quantity). The chemical heat released, which results in the temperature gradient and causes the heat flux, also exerts an indirect influence on the nonorganized energy. Similarly, the chemical heat affects the nonorganized heat flux directly and indirectly, as well.

Effects of relaxation time
The discrete Boltzmann equation is a special discretization and simplification of the Boltzmann equation, which is a fundamental equation in nonequilibrium statistical physics. In Eq. (1) , the relaxation time τ describes the speed of a fluid system relaxing to its local thermodynamic equilibrium state. In fact, ν = 1 /τ is the molecular collision frequency. The effect of molecular collisions is to force the nonequilibrium distribution function f σ i back to the most probable distribution f σ eq i [42] , and the rate is proportional to the molecular collision frequency. Moreover, via the Chapman-Enskog expansion, it is found that the dynamic viscosity coefficient μ σ and heat conductivity κ σ of species σ are related to the relaxation time, i.e., μ σ = p σ τ and κ σ = γ σ μ σ , with the specific heat ratio γ σ = (D + I σ + 2) / (D + I σ ) [24] .
In this subsection, we investigate the nonequilibrium effect on the detonation. Because the relaxation time is a key parameter relative to the nonequilibrium effect. To investigate the nonequilibrium influence upon the detonation, we simulate detonations with various relaxation time. Now, let us compare the case with  Fig. 6 to the case with τ = 4 × 10 −5 in Fig. 2 .
The length of the reaction zone is about 0.02 in Fig. 6 , which is wider than that in Fig. 2 . The length of the nonequilibrium zone is about 0.04 in Fig. 6 , which is also wider than that in Fig. 2 . Namely, both the reaction and nonequilibrium areas are enlarged with increasing relaxation time, and the nonequilibrium area increases faster than the reaction area. In addition, the physical quantities ( p, T, ρ, u x ) increase more slowly and the physical gradients are smaller for increasing relaxation time. Physically, the local HTNE becomes stronger with the increase in either the relaxation time or physical gradients [23] , while the global HTNE is proportional to both the local nonequilibrium amplitude and the entire nonequilibrium area [22] . The physical gradients decrease and the nonequilibrium area enlarges when the relaxation time increases. Figure 7 exhibits the nonequilibrium effects versus the relaxation time. Squares and triangles denote the chemical reactant and product, respectively. The lines stand for the fit- fects of the chemical reactant are greater (less) than those of the product for smaller (larger) relaxation time. Physically, the local HTNE becomes stronger with the increasing relaxation time, chemical heat, reaction rate, or physical gradients [23] . The global HTNE depends upon both the local nonequilibrium amplitude and the entire nonequilibrium area [22] . Specifically, the physical gradients reduce and the nonequilibrium area enlarges when the relaxation time increases [22] . The increasing relaxation time and decreasing physical gradients exert opposite influences on the local HTNEs. Consequently, the local HTNEs increase when the impact of increasing relaxation time dominates, otherwise, they reduce. Meanwhile, increasing relaxation time, reducing physical gradients and enlarging nonequilibrium area play competitive roles in the global HTNEs. The global HTNEs increase as the impact of increasing relaxation time and enlarging nonequilibrium area dominates. Finally, it is evident in Figs. 2 (a) and 6 (a) that the peaks of pressure, density, and horizontal velocity decrease with increasing relaxation time. For the sake of quantitative study, we introduce a dimensionless parameter, H , named peak height. The peak height of a physical quantity q is defined as where q max is the maximum of q around a detonation wave, q s the value of q behind the detonation wave, and q von the maximum of q in the case without nonequilibrium effects. Actually, q s is the Chapmann-Jouguet (CJ) solution in the steady production zone, q von is the ZND solution of q at the von-Neumannpeak. For one-dimensional steady detonation, q s ≤ q max ≤ q von , i.e., 0 ≤ L ( q ) ≤ 1. q max approaches q s and H ( q ) becomes close to 0 when the relaxation time is large. Otherwise, q max approaches q von and H ( q ) becomes close to 1. As shown in Figs. 2 and 6 , with increasing relaxation time, the global nonequilibrium effects become stronger, the detonation peak becomes lower, and the detonation wave becomes wider. It indicates that the nonequilibrium effects attenuate and widen the detonation wave. To further demonstrate this conclusion, we conduct a simulation with all nonequilibrium terms switched off, i.e., f σ i − f σ eq i = 0 , in each iterative step. The configuration is the same as the one in Fig. 6 . The results are plotted in Fig. 9 . Figure 9 (b) and (c) shows that the nonequilibrium effects vanish. The peak heights of pressure, density, and horizontal velocity in Fig. 9 (a) are 0.984, 0.987, 0.968, respectively. While the peak heights of pressure, density, and horizontal velocity in Fig. 6 (a) are 0.0 07, 0.0 01, 0.0 05, respectively. The differences between the peak heights in Figs. 6 (a) and 9 (a) mainly result from the significant nonequilibrium effects. Comparisons between Figs. 6 (a) and 9 (a) lead to three points. (I) The detonation peak under nonequilibrium effects is lower than that without nonequilibrium effects; (II) The detonation wave under nonequilibrium effects is wider than that without nonequilibrium effects; (III) The physical gradients under nonequilibrium effects are sharper than those without nonequilibrium effects.

Conclusions
It is a challenging issue to investigate the complex detonation system which encompasses various interactive physicochemical phenomena, covers a wide range of spatio-temporal scales, and incorporates significant HTNE behaviors. Based on the kinetic theory, the DBM has the capability of providing a deeper insight into nonequilibrium reactive flows in an accurate, efficient and robust way. In the present study, the DBM is employed to investigate the HTNE around the detonation wave. We probe the individual HTNE manifestations of the chemical reactant and product, and analyze the fundamental features of their velocity distribu-tions in velocity space. The velocity distribution functions depart far from their equilibrium counterparts when the local nonequilibrium effects are strong. Moreover, we investigate the impact of the chemical heat and relaxation time on local and global HTNE effects, respectively. The local HTNE effects become strong with increasing relaxation time, physical gradients, and/or chemical reaction, while the global HTNE effects depend on the local HTNE effects and entire nonequilibrium area. It is found that both local and global HTNE effects of the chemical reactant and product increase approximately as a power of the chemical heat which promotes the chemical reaction rate and sharpens the detonation front. With increasing relaxation time, the global HTNE effects of the chemical reactant and product are enhanced by power laws, while their local HTNE effects show changing trends. The physical gradients are smoothed and the nonequilibrium area is enlarged as the relaxation time increases. Finally, to estimate the relative height of detonation peak, we introduce and define the peak height as H(q ) = (q max − q s ) / (q v on − q s ) , where q max is the maximum of q around a detonation wave, q s is the CJ solution and q von is the ZND solution at the von-Neumann-peak. With increasing relaxation time, the height of detonation peak reduces, because the nonequilibrium effects attenuate and widen the detonation wave.
We find that the peak heights of physical quantities ( p, ρ, u x ) are exponential functions of the relaxation time.