Cite this Article
Zhang Yu-Dong, Xu Ai-Guo, Zhang Guang-Cai, Chen Zhi-Hua. A One-Dimensional Discrete Boltzmann Model for Detonation and an Abnormal Detonation Phenomenon. Communications in Theoretical Physics, 2019, 71(1): 117  
Permissions
A One-Dimensional Discrete Boltzmann Model for Detonation and an Abnormal Detonation Phenomenon
Zhang Yu-Dong1, 2, Xu Ai-Guo2, 3, †, Zhang Guang-Cai2, Chen Zhi-Hua1
1Key Laboratory of Transient Physics, Nanjing University of Science and Technology, Nanjing 210094, China
2Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100088, China
3Center for Applied Physics and Technology, MOE Key Center for High Energy Density Physics Simulations, College of Engineering, Peking University, Beijing 100871, China

 

† Corresponding author. E-mail: Xu_Aiguo@iapcm.ac.cn

Abstract
Abstract

A one-dimensional discrete Boltzmann model for detonation simulation is presented. Instead of numerical solving Navier-Stokes equations, this model obtains the information of flow field through numerical solving specially discretized Boltzmann equation. Several classical benchmarks including Sod shock wave tube, Colella explosion problem, and one-dimensional self-sustainable stable detonation are simulated to validate the new model. Based on the new model, the influence of negative temperature coefficient of reaction rate on detonation is further investigated. It is found that an abnormal detonation with two wave heads periodically appears under negative temperature coefficient condition. The causes of the abnormal detonation are analyzed. One typical cycle of the periodic abnormal detonation and its development process are discussed.

Keyword:detonation;discrete Boltzmann model;negative temperature coefficient;abnormal detonation
1 Introduction

Detonation is one kind violent combustion mode accompanied with a large amount of heat release within a short time.[1] It can be treated as a shock wave driven by chemical reaction and propagates with a supersonic speed.[2] Detonation is closely related to the energy use and production safety.

In some cases, it is necessary to generate detonation waves to improve the utilization efficiency of fuels. Because detonation possesses an approximate isovolumetric characteristic during chemical reaction, it has a higher mechanical efficiency than the general combustion mode.[3] Based on the detonation mechanism, several kinds of aero-engines conception including pulse detonation engine,[4] rotating detonation engine,[5] oblique detonation ramjet-in-tube,[6] etc. have been presented and well investigated recently. While in other cases, the formation of detonation should be avoided as far as possible such as in coal mines.[7] Detonation is closely related to both the industrial production and our daily life. However, there still exists much unknown for its deep formation and propagation mechanisms.[3,8]

It has been well known that combustion and detonation are complex chemical reaction processes with various non-equilibrium behaviors including Hydrodynamic Non-Equilibrium (HNE), Thermodynamic Non-Equilibrium (TNE) and chemical reaction non-equilibrium.[9–10] For detonation research, traditional methods are mainly by using Navier-Stokes (NS) equations to describe the flow behaviour and using phenomenological reaction rate formula to describe the reaction process.[11] In addition, high resolution difference schemes are often needed to track the detonation interface and improve the numerical accuracy.[12–13] Of course, great progress has been made on the studies of detonation by the traditional method, especially in recent years.[3,9] However, NS equations themselves are not sufficient in describing the non-equilibrium effects in the reactive flow. The coefficients of viscosity and heat conduction in NS equations are generally calculated by empirical formula, such as Sutherland equation, or measured by experiments.[14–15] This method is not accurate enough when simulating the flow phenomena with strong non-equilibrium characteristics. Compared with NS equations, Boltzmann equation is more fundamental to describe the flow process. Rooted from the non-equilibrium statistic mechanics, Boltzmann equation is a mesoscale model and contains more kinetic information. By means of the Chapman-Enskog analysis,[16] a well-known multi-scale asymptotic expansion, the Euler equations can be obtained from the Boltzmann equation when the system is exactly in its local thermodynamic equilibrium state, and the NS equations can be obtained when the system linearly, in the Knudsen number, deviates from its local thermodynamic equilibrium state. However, when the system deviates much farther from its local thermodynamic equilibrium state, NS equations will not be accurate enough and fail to capture many important flow behaviors, whereas the Boltzmann equation naturally adapts to all of the above situations.

Unfortunately, the original Boltzmann equation is too complex to be solved directly, and some attempts have been made to simplify this model, among which the Lattice Boltzmann Model (LBM)[17–20] is a typical one. The first LBM for combustion simulation was presented by Succi in 1997.[21] Subsequently, Filippova,[22] Yamamoto,[23] Chiavazzo,[24] Chen[25] and other researchers further developed the application of LBM to combustion simulation. However, all of those previous works aim only at simulating incompressible combustion and cannot satisfy the requirements of detonation simulation which shows pronounced compressible behaviors. Recently, Xuʼs group made some attempts in simulating high speed compressible flows using LBM and developed it into a kinetic modeling method to investigate non-equilibrium characteristics.[26–30] To distinguish from the LBM aiming at numerical solving partial differential equations, the kinetic modeling LBM is referred to simply as discrete Boltzmann method/model (DBM). Discrete Boltzmann method now has been well developed and widely used in various complex flow simulations including compressible flows,[31] multiphase flows,[32] Rayleigh-Taylor instability,[33–34] combustion and detonation.[35–38] The new observations brought by DBM are helpful to understand the mechanisms for formation and effects of shock wave, phase transition, energy transformation, and entropy increase in various complex flows. Besides by theory, results of DBM have been confirmed and supplemented by results of molecular dynamics,[39–41] direct simulation Monte Carlo[42–43] and experiment.[44] In the system containing both material interface and mechanical interface (such as shock wave and rarefaction wave), non-equilibrium characteristics have been used to recover the main feature of real distribution function, and to provide physical criteria for discriminating various interfaces. The latter has been used to design appropriate tracking schemes for various interfaces.[28,33,45–46] In recent studies,[34,47] the correlations between the various macroscopical nonuniformity and the relevant non-equilibrium strengths in systems with Rayleigh-Taylor instability and/or Richtmyer-Meshkov instability are used to probe the material mixing processes. In 2013, the first DBM for detonation was presented[35] then a series of extensions have been made. For example, the Multiple-relaxation-time DBM[38,48] and double-distribution-function DBM[36] have been developed. However, those works are all two-dimensional model and at least 16 discrete velocities are needed.[37] The calculation efficiency is much low for some cases where the main behaviors can be described by one-dimensional model. In those cases, a one-dimensional DBM model for detonation is in demand.

Generally speaking, during a combustion process, many species of reactants and a large number of reactions are involved. For example, the combustion of CH4 in air involves 53 species and 325 reactions.[49] The reaction rate varies with the special reaction paths. A detonation process may guide the reactions into different reaction chains because of the type of fuel, shock strength, premixing homogeneity, etc. Consequently, the global reaction rates show different behaviors for different conditions and may possess a non-monotonic dependence on the temperature. Although most of the reaction rates have an exponential temperature dependence and the Arrhenius model is commonly be adopted to describe the reaction rate, the phenomenon of Negative Temperature Coefficient (NTC) in reaction rate has been observed in combustion process and has drawn great attention in recent years.[9,50–51] The existence of NTC may also cause significantly special phenomena during the detonation process. However, to the authors’ knowledge, possible effects of NTC on detonation have not been well studied. In 2016, we conducted a preliminary study on the effects of NTC during detonation.[37] In that work, we found that the effect of NTC during detonation is to lower the reaction rate and delay the formation of detonation wave. In this paper, we further study the characteristics of detonation under NTC based on the one-dimensional DBM. An abnormal detonation phenomenon with two wave heads is observed and carefully studied.

2 Model and Verification
2.1 Discrete Boltzmann Model and Chemical Reaction Rate Model

For the one-dimensional discrete Boltzamnn model, the evolution of the distribution function fi for the discrete velocity vi is governed by Eq. (1),
The subscript i in fi indicates the index of discrete velocity vi. The variables, t and x, are the time and spatial coordinates, respectively, and τ is the relaxation parameter. is the local equilibrium distribution function taking into account the effects of chemical reaction and can be calculated by
where ρ, u, and T are the density, velocity, and temperature, respectively. γ is the ratio of special heats, ω is the chemical reaction rate, and Q is the amount of heat release per unit mass of reactant. The on the right-hand of Eq. (2) can be solved by the following seven equations:






where ηi is a free parameter introduced to describe the n extra degrees of freedom corresponding to molecular rotation and/or vibration.[27] In this paper for i = 1, 2, 3 and ηi = 0 for i = 4, 5, 6, 7.

From multi-scale asymptotic expansion, we know that NS equations with chemical reactions, Eqs. (10)–(12), can be deduced from Eq. (1) under the conditions of Eqs. (3)–(9).


where P = ρ T is the pressure, e = (n + 1)T/2 is the internal energy density, and are the coefficients of viscosity and heat conduction, respectively, the constant-pressure specific heat is cp = (n + 3)/2 and the specific heat ratio is γ = (n + 3)/(n + 1).

In order to solve the above seven equations, at least seven velocities are needed. The discrete velocities model adopted in this paper is shown schematically in Fig. 1. The corresponding values of discrete velocities are listed in Table 1 and c0 is a free parameter depends on the numerical stability.

Fig. 1. Discrete velocity model for one-dimensional DBM.
Table 1.

Values of discrete velocities.

.

The chemical reaction rate ω is described by
where k is reaction rate constant and λ indicates reaction process. It has λ =0 at the beginning of the reaction and λ =1 at the end of the reaction. Tth is defined as ignition temperature and there is no reaction happens below this threshold value.

In addition to being able to recover to NS equations with chemical reactions, DBM provides a set of non-equilibrium measurements, which is also its advantage over the traditional hydrodynamic model. The non-equilibrium quantities are represented by the difference between the moments of distribution function and its corresponding local equilibrium distribution function at a certain time. Those non-equilibrium quantities are defined as



Another set of non-equilibrium quantities can also be defined by the kinetic center moment in a similar way. Each of those quantities reflects the non-equilibrium characteristic of system from a different perspective which cannot be provided by the traditional hydrodynamic model. Those non-equilibrium quantities are significant around the detonation wave front because of the effects of impact compression and chemical reactions.

2.2 Model Verification

In order to validate the new model, several typical benchmarks are carried out. Firstly, two shock wave problems, including the Sod shock tube and Colella explosion wave problem, are simulated and compared with the Riemann solutions. Then a one-dimensional self-sustainable stable detonation is simulated and compared with CJ theoretical solutions.

3 An Abnormal Detonation Induced by Negative Temperature Coefficient
3.1 Simulation of the Abnormal Detonation

In order to investigate the NTC of chemical reaction rate, we adopt the following formula to describe the temperature dependence of rate constant k in Eq. (13) which has been presented in Ref. [37].
with

where h1 and h2 are the peak and valley value of k, respectively. T1 and T2 are values of temperature corresponding to h1 and h2, respectively. In this work, we set h1 = 2000, h2 = 10, T1 = 1.14, T2 = 1.45. Besides, the ignition temperature in Eq. (13) is set as Tth = 1.1. The rate constant-temperature curve is shown in Fig. 5. It can be clearly seen that there is an NTC interval where the reaction rate constant decreases with the increase of temperature.

Fig. 5. Relation curve of rate constant with temperature.

Except the chemical reaction rate and Nx = 50000, other simulation conditions are all the same with those in Fig. 4. Under this chemical reaction condition, an abnormal detonation phenomenon is obtained, which is shown in Fig. 6(b). The normal stable detonation processes with a constant speed and a stable waveform are also shown in Fig. 6(a) for comparison. For the abnormal detonation, there is no constant wave speed and the waveform changes with time periodically.

Fig. 6. Comparisons between (a) the normal detonation wave and (b) the abnormal detonation wave.

Figure 6(b) gives the evolution of the waveform in one cycle. Firstly, at a certain time, a local hotspot is generated behind the wave head. Then a local detonation wave appears and develops with a speed faster than its downstream wave front. So the new formed detonation wave chases the front wave and finally the two waves merge into an over-driven detonation wave. However, the over-driven detonation wave cannot self-sustain, so the wave velocity gradually decreases until it reaches the CJ detonation wave velocity. After that, a local hotspot generates again and a new cycle begins.

The corresponding non-equilibrium quantities Δ2, Δ3, Δ3,1, and Δ4,2 defined in Eqs. (14)–(17) are plotted in Fig. 7. For the normal detonation shown in Fig. 6(a), the shape of the wave front does not depend on time so the non-equilibrium characteristics around the wave front keep unchange with time. The profiles of Δ2, Δ3, Δ3,1, and Δ4,2 for the normal detonation are given in Fig. 7(a). It can be seen the system deviates from the equilibrium in the opposite direction at two sides of the wave front. In fact, it has been known that the non-equilibrium quantities are close to zero near the von Neumann peak point.[28,35] However, for the abnormal detonation, the profiles of non-equilibrium quantities vary with time. Figure 7(b), 7(c), and 7(d) give the profiles of non-equilibrium quantities for the abnormal detonation at three typical moments. From Fig. 6, we can see the wave front of the abnormal detonation is very similar to the normal detonation at t = 0.135 but the profiles of non-equilibrium quantities in Fig. 7(b) and those in Fig. 7(a) are much different. The non-equilibrium deviations in Fig. 7(b) are always in the same direction. From Figs 7(c) and 7(d) which correspond to t = 0.145 and t = 0.155, respectively, double wave fronts can be observed. The non-equilibrium characteristics around the back detonation wave (the left wave) is very similar to those in Fig. 7(a) and the non-equilibrium characteristics around the front detonation wave (the right wave) is similar to those in Fig. 7(b). In addition, the non-equilibrium strength is more significant for the abnormal detonation than the normal detonation when two wave fronts exist.

Fig. 7. Profiles of non-equilibrium quantities for (a) normal detonation, (b) abnormal detonation at t = 0.135, (c) abnormal detonation at t = 0.145, and (d) abnormal detonation at t = 0.155.
3.2 Analysis and Discussion about the Abnormal Detonation

In this section, we will discuss the causes of the abnormal detonation wave. For convenience, we roughly divide the reaction into three stages according to the temperature. Those three stages are denoted as S1, S2, and S3, respectively, as shown in Fig. 8. The first stage (S1) is in the low temperature region but has a fast reaction rate because of NTC. The second stage (S2) has a slower reaction rate at a specific temperature range. The third stage (S3) has a fast reaction rate at high temperature and the reaction rate increases dramatically with the increment of temperature.

Fig. 8. Three stages of chemical reaction rate.

The development of the abnormal detonation is shown in Fig. 9. At t = 0.12, as shown in Fig. 9(a), the detonation wave is still a normal detonation followed by a long rarefaction wave region. The temperature behind the wave head has not reached the third stages (S3). So it has only the first two stages, S1 and S2, of the chemical reaction. At t = 0.13, as shown in Fig. 9(b), a local hotspot appears in the rarefaction wave region and the temperature in this hotspot achieves to S3, then the third stage of chemical reaction appears. Because it has a rapid reaction rate in S3, a local detonation wave is formed and developed quickly. The local detonation wave moves forward and obtains more fuel, so it continues being strengthened and has an increasing wave speed. From Fig. 9(c) we can see S3 gradually widens while S2 becomes narrower. At t = 0.15, the new formed detonation wave almost catches up with the front detonation wave. At this time, S2 is nearly disappeared, which can be seen from Fig. 9(d). After that, two waves merge and the overdriven detonation occurs.

Fig. 9. Detonation waveforms at (a) t = 0.155, (b) t = 0.165, (c) t = 0.175, and (d) t = 0.18.

The evolution of the process from overdriven detonation to the normal detonation is shown in Fig. 10. Fig. 10(a) shows the overdriven detonation wave at t = 0.155, which has a wave speed faster than the CJ detonation. At this time, almost all of the chemical reactions occur in S3. However, overdriven detonation cannot self-sustain and rarefaction waves would gradually form behind the wave head. At t = 0.165, as shown in Fig. 10(b), the temperature behind the detonation wave front begins to go down due to the effect of the rarefaction waves. When the temperature goes down to the region of S2, the second stage of the chemical reaction occurs. Then the rarefaction waves behind the wave head continue growing and temperature behind the wave keeps going down. When the temperature declines to the region of S1, the first stage of reaction occurs. As more and more fuel reacts in S1 and S2 reactions occur in S3 gradually decrease, which can be seen from Figs. 10(c) and 10(d). As a result, the overall reaction rate slows down. With the chemical reaction rate slows down, the detonation wave speed gradually decreases to CJ detonation wave speed.

Fig. 10. Detonation waveforms at (a) t = 0.155, (b) t = 0.165, (c) t = 0.175, and (d) t = 0.18.

After that, a local hotspot reappears and a new local detonation wave is developed again, which means the above process is repeated. The whole process of the development of the abnormal detonation within a period can be summarized by Fig. 11.

Fig. 11. Schematic diagram of development process of the abnormal detonation.
4 Conclusion

In this work, we present a one-dimensional discrete Boltzmann model for detonation. The validity of the new model is verified by three tests. The new detonation model possesses both high computational efficiency and numerical accuracy. Based on the new model, the effects of negative temperature coefficient of reaction rate in a detonation are further investigated. An abnormal detonation phenomenon is presented and its development process is analyzed. It is found that the main reason for the abnormal detonation is that the chemical reaction has three stages, namely S1, S2, and S3. The first stage, S1, is in the low temperature region but has a fast reaction rate, the second stage, S2, has a slower reaction rate at a specific temperature range, and the third stage, S3, has a fast reaction rate at high temperature and the reaction rate increases dramatically with the increase of temperature. For a normal detonation, the chemical reactions occur mainly in S1 and S2. For the abnormal detonation, however, at a certain time a local hotspot is formed as a consequence of the S3. Then a new detonation with a more violent chemical reaction appears behind the old detonation wave front. The new detonation wave has a faster speed than the wave ahead, it catches up with the front wave and finally two waves merge. Then the speed of the detonation wave begins to slow down until it reaches a CJ detonation value. After that, the local hotspot appears again and the previous process reappears.

Reference
[1] Mader C. L. Numerical Modeling of Detonations University of California Press California 1979
[2] Sun J. Zhu J. Theory of Detonation Physics National Defence Industry Press Beijing 1995 (in Chinese)
[3] Jiang Z. Teng H. Liu Y. Adv. Mech. 42 2012 129 (in Chinese)
[4] Kailasanath K. AIAA J. 41 2003 145
[5] Shao Y. Liu M. Wang J. Combust. Sci. Technol. 182 2010 1586
[6] Brackett D. Bogdanoff D. J. Propul. Power 5 1989 276
[7] Yu H. Zhang S. Min. Saf. Environ. Prot. 43 2016 108 (in Chinese)
[8] Shepherd J. E. Proc. Combust. Inst. 32 2009 83
[9] Ju Y. Adv. Mech. 44 2014 26
[10] Xu A. Zhang G. Ying Y. Acta Phys. Sin. 64 2015 184701 (in Chinese)
[11] Xu A. Zhang G. Ying Y. Wang C. Sci. China: Phys. Mech. Astron. 59 2016 650501
[12] Wang C. Shu C. W. Han W. Ning J. Combust. Flame 160 2013 447
[13] Li P. Wang C. Transactions of Beijing Institute of Technology 37 2017 1211 (in Chinese)
[14] Fu D. Wang Y. Computational Aerodynamics Aerospace press Beijing 1994 (in Chinese)
[15] Anderson J. D. Wendt J. Computational Fluid Dynamics McGraw-Hill New York 1995
[16] Chapman S. Cowling T. The Mathematical Theory of Nonuniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Cambridge university press Cambridge 1970
[17] Qian Y. H. D’ Humières D. Lallemand P. Europhys. Lett. 17 1992 479
[18] He X. Luo L. S. Phys. Rev. E 56 1997 6811
[19] Chen S. Doolen G. D. Ann. Rev. Fluid Mech. 30 1998 329
[20] Succi S. The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond Oxford university press New York 2001
[21] Succi S. Bella G. Papetti F. J. Sci. Comput. 12 1997 395
[22] Filippova O. Haenel D. Comput. Phys. Commun. 129 2000 267
[23] Yamamoto K. Takada N. Misawa M. Proc. Combust. Inst. 30 2005 1509
[24] Chiavazzo E. Karlin I. V. Gorban A. N. Boulouchos K. Combust. Flame 157 2010 1833
[25] Chen S. Mi J. Liu H. Zheng C. Int. J. Hydrogen Energy 37 2012 5234
[26] Xu A. Zhang G. Gan Y. et al. Front. Phys. 7 2012 582
[27] Gan Y. Xu A. Zhang G. Yang Y. Europhys. Lett. 103 2013 24003
[28] Lin C. Xu A. Zhang G. et al. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 89 2014 013307
[29] Xu A. Zhang G. Li Y. Hua L. Prog. Phys. 34 2014 136 (in chinese)
[30] Xu A. Zhang G. Gan Y. Mech. Eng. 38 2016 361 (in chinese)
[31] Xu A. Zhang G. Zhang Y. Chap. 2 Kinetic Theory Kyzas G. InTech Rijeka 2018
[32] Gan Y. Xu A. Zhang G. Succi S. Soft Matter 11 2015 5336
[33] Lai H. Xu A. Zhang G. et al. Phys. Rev. E 94 2016 023106
[34] Chen F. Xu A. Zhang G. Front. Phys. 11 2016 114703
[35] Yan B. Xu A. Zhang G. et al. Front. Phys. 8 2013 94
[36] Lin C. Xu A. Zhang G. Li Y. Combust. Flame 164 2016 137
[37] Zhang Y. Xu A. Zhang G. et al. Combust. Flame 173 2016 483
[38] Lin C. Luo K. H. Comput. Fluids 166 2018 176
[39] Liu H. Kang W. Zhang Q. et al. Front. Phys. 11 2016 1
[40] Liu H. Zhang Y. Wei K. et al. Phys. Rev. E 95 2017 023201
[41] Liu H. Kang W. Duan H. et al. Sci. China: Phys. Mech. Astron. 47 2017 070003 (in Chinese)
[42] Zhang Y. Xu A. Zhang G. et al. Submitted to Comput. Phys. Commun. arXiv:1801.02649 2018
[43] Meng J. Zhang Y. Hadjiconstantinou N. G. et al. J. Fluid Mech. 718 2013 347
[44] Lin C. Luo K. H. Fei L. Succi S. Sci. Rep. 7 2017 14580
[45] Chen F. Xu A. Zhang G. Wang Y. Front. Phys. 9 2014 246
[46] Lin C. Xu A. Zhang G. et al. Phys. Rev. E 96 2017 053305
[47] Chen F. Xu A. Zhang G. Phys. Fluids 30 2018 102105
[48] Xu A. Lin C. Zhang G. Li Y. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 91 2015 043306
[49] Faghih M. Gou X. Chen Z. J. Loss Prev. Process Ind. 40 2016 131
[50] Ksandopulo G. I. Russ. J. Phys. Chem. B 5 2011 701
[51] Pan J. Wei H. Shu G. et al. Combust. Flame 174 2016 179
[52] Zhang H. Acta Aerodynamica Sinica 6 1988 143 (in Chinese)