1. Introduction
Detonation is a kind of supersonic combustion where a chemical reaction is coupled with shock dynamics [1, 2]. Detonation is the key factor causing catastrophic damage in various explosion accidents, such as gas explosions [3]. Besides, detonation has been widely used in aerospace, defense, mining, and demolition applications [4–6]. In view of its importance, detonation has been researched by experimental [7–10], theoretical [11,12], and numerical methods [13–17]. Especially, with the rapid development of computer technology and computational fluid dynamics, numerical simulation has become an indispensable method to study detonation [2, 18]. However, the evolutionary process of detonation involves complex physicochemical phenomena [11, 12]. During the evolution of the detonation, the chemical reaction and hydrodynamic processes interact with each other, which produces various interfacial and mechanical structures, and contains a wealth of hydrodynamic and thermodynamic nonequilibrium effects [19]. Therefore, the accurate prediction and effective control of detonation are still open issues.
In recent years, the kinetic methods based on the Boltzmann equation have become a popular mesoscopic simulation approach in various fields [20–23]. In particular, the discrete Boltzmann method (DBM) [24, 25] is suitable for supersonic waves with essential nonequilibrium effects, such as the steady or unsteady detonation. A finite number of velocities is employed to discretize the Boltzmann equation in the modeling of DBM. Thanks to its finite discrete velocities and distribution functions, the DBM owns a high computational efficiency [24, 26]. Moreover, the hydrodynamic and thermodynamic nonequilibrium effects beyond the macroscopic governing equations can be obtained and quantified as the higher-order kinetic moments in the DBM [27]. Therefore, as a mesoscopic kinetic model, the DBM not only inherits the function of the Boltzmann equation to recover the macroscopic transport equations in the continuous limit, but also has the ability to capture the thermodynamic nonequilibrium effects beyond the macroscopic equations. The DBM has been widely utilized to obtain some deeper insights into nonequilibrium behaviours in various complex fluids and has excellent applications in the analysis of nonequilibrium, nonlinear, transient complex flow systems [22–27, 29].
Recently, the DBM has achieved remarkable success in simulating combustion and detonation [24, 25, 30–38]. The pioneering DBM for detonation [30] was presented by Yan et al who adopted the Lee-Tarver model [39] to control the chemical reaction. In 2014, Lin et al [31] proposed a novel polar coordinate DBM and investigated typical implosion and explosion processes. In 2015, Xu et al presented a two-dimensional multiple-relaxation-time model [32] and studied the hydrodynamic and thermodynamic nonequilibrium in the process of steady detonation. In 2016, a double-distribution-function DBM for combustion and detonation was constructed by Lin et al [33], where the chemical reactant is described by one distribution function, and the chemical product by the other distribution function. Then, Lin et al reported a multi-component DBM for premixed, nonpremixed, or partially premixed nonequilibrium reactive flows [35]. In 2019, Lin et al developed a multiple-relaxation-time DBM for unsteady reactive flows [24], where the specific heat ratio and Prandtl number are adjustable. In 2020, Lin et al simulated unsteady detonation with essential hydrodynamic and thermodynamic nonequilibrium effects, and studied the impact of the perturbation amplitude, wavelength, and chemical heat release on the physical field of unsteady detonation with nonequilibrium effects [25]. Next, a three-dimensional DBM for the simulation of detonation was developed by Ji et al [36], and they captured the features of the three-dimensional detonation via the DBM [37]. Very recently, Su et al recovered the particle velocity distribution function of reactive flows from the gas kinetic theory and investigated the higher-order kinetic moments and nonequilibrium quantities around the detonation wave by using the DBM [38].
In the present study, we carry out more in-depth and comprehensive research on unsteady detonation via the DBM. In section 2 , the details of the modeling construction are presented. Then, the DBM is employed to simulate the unsteady detonation. The deviation of the velocity distribution function from the equilibrium state and the kinetic moments of the reaction term are studied in section 3 . Reaction manifestations in moment space are investigated from a kinetic perspective in section 4 . Finally, section 5 provides a summary.
2. Discrete Boltzmann model
The DBM is based on the following equation,
$\begin{eqnarray}\frac{\partial {f}_{i}}{\partial t}+{v}_{i\alpha }\frac{\partial {f}_{i}}{\partial {r}_{\alpha }}=-\frac{1}{\tau }\left({f}_{i}-{f}_{i}^{\mathrm{eq}}\right)+{R}_{i},\end{eqnarray}$
where τ denotes the relaxation time, t the time, fi and ${f}_{i}^{\mathrm{eq}}$ represent the discrete distribution function and its equilibrium counterpart, respectively. Moreover, vi denotes the discrete velocity with i = 1, 2, 3, ⋯, N, and N = 16 is the total number of discrete velocities. Here a two-dimensional sixteen-velocity model is utilized, see figure 1.
Figure 1. Sketch of the two-dimensional sixteen-velocity model. |
In especial, the original expression of equilibrium distribution function feq is as follows [24]
$\begin{eqnarray}{f}^{\mathrm{eq}}=n{\left(\displaystyle \frac{1}{2\pi T}\right)}^{D/2}{\left(\displaystyle \frac{1}{2\pi {IT}}\right)}^{1/2}\exp \left[-\displaystyle \frac{{\left|{\boldsymbol{v}}-{\boldsymbol{u}}\right|}^{2}}{2T}-\displaystyle \frac{{\eta }^{2}}{2{IT}}\right],\end{eqnarray}$
where n, u, and T represent the particle number density, flow velocity and temperature, respectively. Besides, D = 2 denotes the dimension, I counts extra degrees of freedom due to vibration and/or rotation, and η is used to describe the corresponding vibrational and/or rotational energies.On the right-hand side of equation (1 ), Ri indicates the reaction term describing the change rate of the distribution function due to chemical reactions. In this paper, a two-step reaction scheme is employed to control the chemical reaction process [40],3 ) denotes a thermally neutral induction zone or ignition process, and the second step in equation (4 ) describes the rapid energy release after the branched-chain thermal explosion. $\xi ^{\prime} $ denotes the change rate of the reaction progress variable in a thermally neutral induction period, while $\lambda ^{\prime} $ is the change rate of the mass fraction of the chemical product. Besides, Ts is the temperature after the preshock wave. EI and ER represent the activation energies. KI and KR indicate the rate constants.
$\begin{eqnarray}\xi ^{\prime} ={{HK}}_{I}\exp \left[{E}_{I}\left({T}_{s}^{-1}-{T}^{-1}\right)\right],\end{eqnarray}$
$\begin{eqnarray}\lambda ^{\prime} =\left(1-H\right){K}_{R}\left(1-\lambda \right)\exp \left(-{E}_{R}{T}^{-1}\right).\end{eqnarray}$
The first step in equation (Under the assumption that the time scale of molecular collisions < the time scale of chemical reactions < the characteristic time scale of flows, the original reaction term can be obtained [38],
$\begin{eqnarray}R={f}^{\mathrm{eq}}\times \left[-\displaystyle \frac{D+1}{2T}+\displaystyle \frac{{\left|{\boldsymbol{v}}-{\boldsymbol{u}}\right|}^{2}}{2{T}^{2}}+\displaystyle \frac{{\eta }^{2}}{2{{IT}}^{2}}\right]\displaystyle \frac{2Q\lambda ^{\prime} }{D+I},\end{eqnarray}$
where Q indicates the chemical heat release per unit mass of chemical reactant.In this work, in order to recover the Navier–Stokes (NS) equations, the following moment relations should be satisfied,2 ) into (6 ) leads to an explicit expression:A .
$\begin{eqnarray}\iint {f}^{\mathrm{eq}}{\rm{\Psi }}{\rm{d}}{\boldsymbol{v}}{\rm{d}}\eta ={\sum }_{i}{f}_{i}^{\mathrm{eq}}{{\rm{\Psi }}}_{i},\end{eqnarray}$
with the particle velocities ψ = 1, v, $\left({\boldsymbol{v}}\cdot {\boldsymbol{v}}+{\eta }^{2}\right)$, vv, $\left({\boldsymbol{v}}\cdot {\boldsymbol{v}}+{\eta }^{2}\right){\boldsymbol{v}}$, vvv, $\left({\boldsymbol{v}}\cdot {\boldsymbol{v}}+{\eta }^{2}\right){\boldsymbol{v}}{\boldsymbol{v}}$, and the corresponding discrete velocities, ψi = 1, vi, $\left({{\boldsymbol{v}}}_{i}\cdot {{\boldsymbol{v}}}_{i}+{\eta }_{i}^{2}\right)$, vivi, $\left({{\boldsymbol{v}}}_{i}\cdot {{\boldsymbol{v}}}_{i}+{\eta }_{i}^{2}\right){{\boldsymbol{v}}}_{i}$, vivivi, $\left({{\boldsymbol{v}}}_{i}\cdot {{\boldsymbol{v}}}_{i}+{\eta }_{i}^{2}\right){{\boldsymbol{v}}}_{i}{{\boldsymbol{v}}}_{i}$. Substituting equations ( $\begin{eqnarray}{{\boldsymbol{H}}}_{{{\boldsymbol{f}}}^{\mathrm{eq}}}={\boldsymbol{M}}{{\boldsymbol{f}}}^{\mathrm{eq}},\end{eqnarray}$
where M is a 16 × 16 square matrix bridging the velocity and moment space, consisting of the functions of the discrete velocities. In addition, feq and ${{\boldsymbol{H}}}_{{{\boldsymbol{f}}}^{\mathrm{eq}}}$ are matrix expressions of the equilibrium distribution function in velocity and moment space, respectively. The details are in appendix Similarly, the reaction terms follow the relations,6 ). Substituting equations (5 ) into (8 ) results in, HR = MR, where R and HR correspond to the matrices of reaction term in velocity and moment space, respectively. Hence,
$\begin{eqnarray}\iint R{\rm{\Psi }}{\rm{d}}{\boldsymbol{v}}{\rm{d}}\eta ={\sum }_{i}{R}_{i}{{\rm{\Psi }}}_{i}.\end{eqnarray}$
The elements of ψ and ψi are the same as those in equation ( $\begin{eqnarray}{\boldsymbol{R}}={{\boldsymbol{M}}}^{-1}{{\boldsymbol{H}}}_{{\boldsymbol{R}}},\end{eqnarray}$
where M−1 is the inverse of M. Mathematically, the matrix HR is expressed by ${{\boldsymbol{H}}}_{{\boldsymbol{R}}}={\left(\begin{array}{llll}{H}_{R1} & {H}_{R2} & \cdots & {H}_{R16}\end{array}\right)}^{{\rm{T}}}$. Physically, HR1 represents the change rate of mass density due to chemical reactions. HR2 and HR3 describe the corresponding change rate of momentum in the x and y directions, respectively. HR4 indicates the change rate of energy because of the chemical reactions. HR5, HR6 and HR7 describe the change rate of nonorganised momentum fluxes due to chemical reactions. HR8, HR9, HR10, HR11, HR12, and HR13 are related to the change rate of nonorganised energy fluxes due to chemical reactions. HR14, HR15, and HR16 correspond to the change rate of the fluxes of nonorganised energy fluxes due to chemical reactions.3. Investigation of unsteady detonation
In this section, the DBM is employed to simulate the one-dimensional unsteady detonation. The initial configuration, obtained from the Hugoniot relation, takes the form
$\begin{eqnarray*}\left\{\begin{array}{l}{\left(\rho ,{u}_{x},{u}_{y},T,\xi ,\lambda \right)}_{L}=\left(1.56244,\,1.1547,\,0,\,3.01066,\,1,1\right),\\ {\left(\rho ,{u}_{x},{u}_{y},T,\xi ,\lambda \right)}_{R}=\left(1,\,0,\,0,\,1,\,0,\,0\right),\end{array}\right.\end{eqnarray*}$
where the subscript L indicates 0 ≤ x ≤ 0.01, and R indicates 0.01 < x ≤ 1.6. The Mach number is 2.71. To ensure the resolution is high enough, the grid is chosen as Nx × Ny = 40 000 × 1, the space step Δx = Δy = 4 × 10−5, the time step Δt = 2 × 10−6, the specific heat ratio γ = 1.66, and the chemical heat release per unit mass is Q = 4. Furthermore, the inflow and/or outflow boundary conditions are employed in the x direction, and the periodic boundary condition is adopted in the y direction.The pressure, temperature, density and horizontal velocity are illustrated in figure 2(a). Obviously, as the detonation wave propagates from left to right, the macroscopic physical quantities in the medium first rise drastically and then fall sharply to form a high peak. Then there are regular fluctuations after the peak. The three dashed lines in figure 2(a) denote the locations x = 1.1724, x = 1.1799, and x = 1.2518, respectively.
Figure 2. Profiles of physical quantities (a) and nonorganized energy (b) of the unsteady detonation wave. |
Figure 2(b) exhibits the nonorganized energy in the x direction. The inset shows the corresponding profiles within the spatial range 1.245 < x < 1.255. With the detonation wave propagating forward, the nonorganized energy in the x direction first forms a high positive peak and then a negative trough. It can be seen that in the post-wave region, the nonorganized energy also changes in a form of small regular oscillations. Correspondingly, the dashed line in figure 2(b) represents the position x = 1.2518. Moreover, there is a comparison between these DBM results and CFD simulations, which is displayed in appendix B .
Next, let us investigate the velocity distribution function of unsteady detonation. In [38], via the Chapman-Enskog analysis, the local particle velocity distribution function of reactive flow is derived from the gas kinetic theory. In short, on the NS level, the local particle velocity distribution function reads10 ) and (11 ), it can be seen that both spatial gradients and reaction term make contributions to the deviation of the local velocity distribution function from the equilibrium state, while the velocity distribution function is composed of the equilibrium particle velocity distribution function and the corresponding deviation counterpart.
$\begin{eqnarray}f={f}^{\mathrm{eq}}+{f}^{\mathrm{neq}},\end{eqnarray}$
where feq is the equilibrium particle velocity distribution function, and fneq is the deviation of the velocity distribution function from the equilibrium state. Besides, fneq can be written as [38], $\begin{eqnarray}{f}^{\mathrm{neq}}=-\tau \left(\displaystyle \frac{\partial {f}^{\mathrm{eq}}}{\partial t}+{v}_{\alpha }\displaystyle \frac{\partial {f}^{\mathrm{eq}}}{\partial {r}_{\alpha }}\right)+\tau R.\end{eqnarray}$
From equations (Next, we conduct an in-depth study of fneq. Since the macroscopic quantities show regular oscillations in figure 2, we first investigate the deviations at a trough (x = 1.1724) and a peak (x = 1.1799) of the pressure in the region after the detonation front, and the position of the highest peak (x = 1.2518) of pressure around the detonation wave is also researched. To this end, figure 3 depicts the deviation of the particle velocity distribution function from the equilibrium state ${f}^{\mathrm{neq}}={f}^{\mathrm{neq}}\left({v}_{x},{v}_{y}\right)$ in the left column and corresponding contours in the velocity space $\left({v}_{x},{v}_{y}\right)$ in the right column. To be specific, figures 3(a) and (b) exhibit the deviation fneq and its corresponding contour at the valley of the pressure. Figures 3(c) and (d) display the deviation fneq and the relevant contour at the peak of the pressure. Besides, the deviation fneq and its contour at the position of the highest peak of pressure around the detonation wave are drawn in figures 3(e) and (f), respectively.
Figure 3. The first column shows the distribution of ${f}^{\mathrm{neq}}\left({v}_{x},{v}_{y}\right)$ from a three-dimensional perspective, and the second column plots the corresponding contours of ${f}^{\mathrm{neq}}\left({v}_{x},{v}_{y}\right)$ in the two dimensional velocity space. The three rows from top to bottom depict the deviations at locations x = 1.1724, x = 1.1799 and x = 1.2518, respectively. |
As shown in figures 3(a) and (c), there are both positive and negative deviations. Figure 3(a) depicts that along the vx direction, a deep negative valley first appears, and then there is a positive peak, afterwards it decreases to form a shallow negative valley. On the contrary, there is an opposite trend in figure 3(c). Firstly, a small positive peak appears, then drops to a negative trough, and later it increases to a high positive peak. Meanwhile, it is interesting to note that the curves of figures 3(b) and (d) look quite similar to each other. Both structures are composed of three parts, and the corresponding shapes of the three parts are the same. However, the difference between figures 3(b) and (d) mainly lies in the values of each contour line, especially whether the values are positive or negative. Specifically, from left to right, the values of three parts of figure 3(b) are first negative, then positive, and finally negative, which is consistent with the trend of figure 3(a). Similarly, the values of the three parts in figure 3(d) are consistent with the trend of figure 3(c), which contains positive, negative, and positive parts. Moreover, From the three-dimensional view, figures 3(e) looks like a bulging positive circle surrounding a negative valley. Besides, it can be found from the contour that figure 3(f) is symmetric in the vx and vy directions.
In sum, it can be found in figure 3 that the characteristics of the deviation near the wavefront of the detonation are significantly different from those in the post-wave region. Besides, the nonequilibrium intensity around the detonation wave is much greater than that in the post-wave region. Moreover, in the post-wave region of the unsteady detonation wave, the deviation of the particle velocity distribution function from the equilibrium state changes with the macroscopic quantities and nonorganized energy.
It should be pointed out that the DBM is a mesoscopic kinetic method, so the chemical reaction of unsteady detonation can be studied from the kinetic moment viewpoint. Figure 4 depicts the sixteen moments of reaction term HR.
Figure 4. Kinetic moments of reaction term around the detonation wave. (a) HR1 to HR4, (b) HR5 to HR7, (c) HR8 to HR13, and (d) HR14 to HR16. The symbols show the numerical results, and the solid lines represent the theoretical solutions. |
Figure 4(a) exhibits the profiles of moments of reaction term HR1 to HR4. The symbols denote the DBM results, and the solid lines represent the theoretical solutions. Clearly, the simulation results are in nice agreement with the theoretical solutions. Besides, it can be seen that during the evolution of unsteady detonation wave, HR1, HR2 and HR3 keep zero due to the conservation of mass and moment. HR4 first increases sharply and then decreases smoothly, and forms a towering asymmetric peak. The reason is that there is violent chemical heat released around the detonation wave.
Figure 4(b) displays the moments of reaction term HR5, HR6 and HR7. It is obvious that HR6 is zero, while HR5 and HR7 first grow and then decline, forming a positive crest. Besides, the trends of HR5 and HR7 are exactly the same, which is also consistent with the theoretical solutions ${H}_{R5}=\rho T^{\prime} $ and ${H}_{R7}=\rho T^{\prime} $.
The patterns of moments of reaction term HR8, HR9, HR10, HR11, HR12, and HR13 are shown in figure 4(c). Clearly, HR9, HR11, and HR13 keep zero, while HR8, HR10, and HR12 all increase first and then drop sharply. Moreover, the peak of HR8 is higher and wider than HR10, and the peak of HR12 is thinner and lower than HR10. The evolution in figure 4(c) is in accordance with the theoretical solutions ${H}_{R8}\,=\left(D+I+2\right)\rho {u}_{x}T^{\prime} $, ${H}_{R9}=\left(D+I+2\right)\rho {u}_{y}T^{\prime} $, ${H}_{R10}\,=3\rho {u}_{x}T^{\prime} $, ${H}_{R11}=\rho {u}_{y}T^{\prime} $, ${H}_{R12}=\rho {u}_{x}T^{\prime} $, and ${H}_{R13}=3\rho {u}_{y}T^{\prime} $, respectively.
The rest change rate of the fluxes of nonorganized energy flux due to chemical reactions HR14 to HR16 is plotted in figure 4(d). Similarly, HR15 keeps zero as the detonation wave propagates from left to right, while HR14 and HR16 show positive peaks. In addition, the peak of HR14 is taller and wider than HR16. The theoretical solutions ${H}_{R14}\,=\rho \left[2T\left(D+I+2\right)+\left(D+I+5\right){u}_{x}^{2}+{u}_{y}^{2}\right]T^{\prime} $, ${H}_{R15}\,=\,\rho {u}_{x}{u}_{y}\left(D\,+\,I\,+\,4\right)T^{\prime} $, and ${H}_{R16}\,=\,\rho [2T(D\,+\,I\,+\,2)\,+{u}_{x}^{2}\,+\left(D+I+5\right){u}_{y}^{2}]T^{\prime} $ explain the variation patterns of figure 4(d).
4. Reaction manifestation in moment space
In order to have a better understanding of chemical reactions of unsteady detonation from the perspective of kinetic moments, the reaction manifestation GR is defined to describe the global effects due to the chemical reaction in moment space12 ), which denotes the magnitude of the quantity in the moment space (HR1, HR2, ..., HR16). That is to say, the mathematical expression of equation (12 ) consists of 16 independent variables.
$\begin{eqnarray}{G}_{R}=\iint \sqrt{\sum {{H}_{{Ri}}}^{2}}{\rm{d}}x{\rm{d}}y,\end{eqnarray}$
where the integral is over the whole computational domain. It should be mentioned that it is a coarse-grained way to investigate physical characteristics and reaction manifestations by using equation (Figure 5 shows the evolution of the reaction manifestation in moment space. Figure 5(a) shows the whole process of reaction manifestation, it can be seen that when t < 0.1, the reaction manifestation changes very drastically, in especial, an extremely high peak was observed around t = 0.04. After t > 0.15, the reaction manifestation begins to show regular oscillations. In figure 5(b), the red dashed lines represent the envelope curves. From the upper envelop curve, it can be seen that the maximum of reaction manifestation displays a periodic small change when t > 0.15. One cycle of the first type of periodic oscillation is shown by the red rectangular frame, and the first cycle is about T1 = 0.040 28. Then, we further narrow the time scale of the observed reaction manifestation. Figure 5(c) illustrates the second type of cycle of the reaction manifestation around T2 = 0.006 69, i.e. the pattern of change from a valley to a peak, and one of the cycles is displayed by the green rectangular frame. Besides, after partially zooming in figure 5(c), we see that the data points which are plotted by hollow squares also seem to evolve with a much smaller period from figure 5(d), i.e. they alternate in groups of about seven points. Similarly, one cycle of the third period is about T3 = 0.000 012 signed by the black rectangular frame.
Figure 5. Evolution of the reaction manifestation in moment space. (a) 0 ≤ t ≤ 0.4, (b) 0.15 ≤ t ≤ 0.35, (c) 0.24 ≤ t ≤ 0.26, and (d) 0.25 ≤ t ≤ 0.2501. |
As shown in figure 5, there are three types of periodic changes in the reaction manifestation in the moment space of unsteady detonation during the post-stabilization evolution. Therefore, the fast Fourier transform is introduced to deal with the reaction manifestation in figure 5(b). Figure 6 illustrates the amplitude versus frequency f(HZ). It is apparent that there are many peaks in various frequencies. As the frequency increases, the amplitude decreases gradually. The first eight peaks are more distinct, which denotes that the reaction manifestation is composed mainly of the first several groups' signal of different frequencies.
Figure 6. The amplitude versus frequency from the fast Fourier transform performed on the reaction manifestation. |
In addition, the specific frequency and corresponding amplitude are shown in table 1. The first line denotes the number of the first eight frequencies, and the second and third lines give the values of the corresponding frequency and amplitude, respectively. It can be found from the table that the frequencies are all roughly multiples of f = 150 (HZ).
Table 1. The frequencies and corresponding amplitudes of the reaction manifestation. |
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |
|---|---|---|---|---|---|---|---|---|
| Frequency | 150.003 | 300.006 | 450.009 | 590.012 | 740.015 | 890.018 | 1040.02 | 1190.024 |
| Amplitude | 0.008 72 | 0.004 63 | 0.002 27 | 0.001 87 | 0.001 45 | 0.001 05 | 0.000 69 | 0.000 44 |
The data in table 1 are fitted to further explore the relations among the number, frequency, and amplitude. Figure 7(a) depicts the values of frequency versus the number of frequency. The squares denote the numerical data, and the solid line represents the fitting function f = 150n, where f and n stand for the frequency and number, respectively. It is evident that the original data agree well with the fitting solutions in figure 7(a), which means that the reaction manifestation may consist of a superposition of frequencies of different multiples of f = 150 (HZ). Figure 7(b) displays the amplitude (A) versus frequency. In the same way, the circles represent the numerical data, and the line is for the fitting function $A=0.01637\times \exp \left(-f/210.9455\right)$. From the fitting result, it can be seen that there is an exponentially decreasing relationship between the frequency and amplitude. The amplitude decreases exponentially as the frequency increases.
Figure 7. The fitting results of the reaction manifestation. (a) the relation between the number and frequency, and (b) amplitude versus frequency. |
Next, let us study the effects of different physical factors on the reaction manifestation in the detonation process. We focus on the changes in the maximum and amplitude of the reaction manifestation. In fact, the bifurcation boundary between steady and unsteady detonations depends on the ratio of the length of the heat release layer to that of the induction zone layer [40, 41]. Generally, the length of the induction zone can be adjusted by adjusting KI, meanwhile, the length of the heat release layer can be determined by KR. The two rate constants play an important role in the evolution of the detonation. Therefore, we first consider unsteady detonations with the rate constant KI tunable, and the other parameters fixed the same as the above simulation.
Figures 8(a) and (b) exhibit the maximum and amplitude of the reaction manifestation versus the rate constant KI, respectively. The symbols denote the simulation results and the lines are for the fitting results. The fitting equations in figures 8(a) and (b) are ${Y}_{\mathrm{Max}}=0.11305\,-2.2505\times {10}^{-4}\times {K}_{I}$ and A = 0.10767 − 2.41021 × 10−4 × KI, respectively. Obviously, with the increasing rate constant KI, the maximum and amplitude of the reaction manifestation decrease. Figure 9 illustrates the physical quantities around the detonation wave with two different rate constants KI. As shown in figures 9(a)–(d), with the increasing KI, the amplitudes of the oscillation of density, pressure, horizontal velocity and temperature around the detonation wave become smaller, which means that the chemical reactions become less violent. In other words, the reaction manifestation becomes weaker.
Figure 8. The maximum (a) and amplitude (b) of the reaction manifestation versus rate constant KI. |
Figure 9. The physical quantities with different rate constants KI, (a) density, (b) pressure, (c) horizontal velocity, and (d) temperature. |
Next, we simulate the unsteady detonation wave with various values of rate constant KR from KR = 7000 to KR = 13000. Figure 10 illustrates the simulation results and fitting solutions which are displayed by the symbols and solid lines, respectively. The fitting equations are ${Y}_{\mathrm{Max}}=0.02175\,+2.3068\,\times \,{10}^{-6}\,\times \,{K}_{R}$ and A = 0.00818 + 2.64065 × 10−6 × KR, respectively. Obviously, the values of maximum and amplitude grow linearly with the increasing KR. Figures 11(a)–(d) display the evolutions of density, pressure, horizontal velocity and temperature with different rate constants KR, respectively. It is clear in figure 11 that with the increase of KR, the macroscopic physical quantities change more violently, which means that the intensity of chemical reactions is enhanced. Consequently, the global effect caused by chemical reactions on the moment space is enhanced, and the reaction manifestation becomes stronger.
Figure 10. The maximum (a) and amplitude (b) of the reaction manifestation versus rate constant KR. |
Figure 11. The physical quantities with different rate constants KR, (a) density, (b) pressure, (c) horizontal velocity, and (d) temperature. |
Finally, the simulations of unsteady detonation waves are carried out for various values of chemical heat release per unit mass from Q = 3 to 8. Figures 12(a) and (b) show the simulation and fitting results. Symbols represent the simulation results, and the solid lines denote the fitting equations ${Y}_{\mathrm{Max}}$ = −0.12108 + 0.04272 × Q and A = −0.10301 + 0.03554 × Q, respectively. It can be found that as the chemical heat increases, the values of maximum and amplitude increase in a linear form. Physically, for larger chemical heat, the chemical reaction occurs more violently, the physical quantities change more sharply around the detonation wave, and the reaction manifestation becomes larger.
Figure 12. The maximum (a) and amplitude (b) of the reaction manifestation versus chemical heat. |
5. Conclusions
In this work, we have investigated the characteristics of unsteady detonation from the kinetic viewpoint. The unsteady detonations are simulated via the DBM, which has the capability of capturing detonation waves with nonequilibrium effects. As the detonation wave propagates forwards, the macroscopic physical quantities in the medium first rise drastically and then fall sharply to form a peak. There are regular oscillations after the detonation peak. Besides, the deviation of the velocity distribution function from the local equilibrium state of unsteady detonation oscillates with the macroscopic quantities. Especially, the deviations show opposite patterns at the peak and trough of pressure after the detonation wave. Moreover, the characteristics of the deviation at the peak of pressure are significantly different from those in the post-wave region. In addition, the nonequilibrium intensity around the detonation wave is much greater than that in the post-wave region. Furthermore, the DBM is adopted to investigate the chemical reaction of unsteady detonation, and the change rate of the kinetic moments due to chemical reactions is consistent with their theoretical solution.
In order to have a deeper understanding of chemical reactions from the kinetic perspective, the reaction manifestation is defined to describe the global effects. There are three types of periodic oscillations in the reaction manifestation in moment space during the evolution of the unsteady detonation, and the corresponding cycles are about T1 = 0.04 028, T2 = 0.006 69 and T3 = 0.000 012, respectively. Via the fast Fourier transform, it can be seen that the reaction manifestation is mainly composed of several signal frequencies. Besides, it has been demonstrated that the reaction manifestation consists of a superposition of frequencies of different multiples of f = 150 (HZ), and the corresponding amplitude decreases exponentially as the frequency increases. Finally, the influences of rate constants and the chemical heat on the reaction manifestation have been investigated. From the numerical simulations and corresponding fitting equations, the following points can be obtained: (I) with the increasing rate constant KI, the maximum and amplitude of the reaction manifestation decrease. (II) The values of maximum and amplitude of the reaction manifestation grow linearly with the increasing rate constant KR. (III) As the chemical heat increases, the values of maximum and amplitude increase in a linear form.
