Numerical Investigation of the Detonation Cell Bifurcation with Decomposition Technique

Abstract

Bifurcation of the characteristic cells into multiple smaller cells and decay of those cells into single large characteristic cell is observed frequently. In the present study the bifurcation phenomenon of the detonation front is investigated for marginally unstable detonations using decomposition technique. Numerical analysis is carried out with detailed chemical kinetics for detonation propagation in H₂/O₂ mixtures at 10 kPa. The dynamic characteristics of the instability at the detonation front, such as the local oscillation frequency and the coherent spatial structure of the oscillation are also studied with dynamic mode decomposition (DMD) technique. The coherent structures of the primary and secondary detonation cells are analyzed during the cell bifurcation process and the mechanism in which the secondary cells are formed is investigated. It is demonstrated that the modal analysis categorizes the instability phenomena clearly and can be effectively utilized to identify the origin and source of the instability.

1. Introduction

Detonation remains a critical issue for industrial safety in chemical processing plants, refineries, and in nuclear power plants. Fundamental knowledge on the basic process such deflagrations, detonations, deflagration to detonation transition (popularly known as DDT process) in fuel–air mixtures have enhanced through various numerical and experimental research in the last few decades. Still there are various unresolved issues related to the detonation physics due to the expensive nature of the detonation research. All gaseous detonations are unstable and produce a characteristics cell structure depending on the initial properties of the combustible mixture such as initial pressure, temperature, and stoichiometric condition etc. The unstable detonation front locally decays and gets reaccelerated at the intersection of the oppositely moving transverse waves, called as triple point, and by tracing the triple point location, the characteristics cells are obtained. The resulting characteristic cellular pattern is the characteristics scale of the given combustible mixture from which the dynamic characteristics of the detonation are studied. Detonation propagation mechanism and the prediction of the dynamic parameters can be resolved only from the detailed understanding of the physical and chemical process involved in the non-equilibrium structure of the detonation wave. All detonations are unstable, and the instability at the detonation front has been experimentally observed as early as 1961 [1], subsequently followed by several researchers [2,3,4,5]. All experimentally observed detonation waves exhibit such unstable behavior [6] and was widely studied, both numerically and experimentally [7,8,9]. Austin [7] conducted a wide range of experiments in various detonation mixtures and reported wide range of instability modes in detonation.

The structure of the detonation cells is usually complex which ranges from regular to irregular structures and sometimes with fine structures within a detonation cell. The intermediate mode of detonation is the one in which the detonation cell has two different characteristics length scales. In such cases, the primary cells undergo cell bifurcation and produces multiple secondary detonation cells which again decays. The bifurcation of cells is also observed in moderately and highly unstable detonation cases in which the large characteristics cells are bifurcated into multiple small-scale cells. Possibility of shock bifurcation in the detonation front has been reported earlier [9,10,11]. Mach et al. [12] studied the bifurcation of the detonation front from the reflecting shocks and proposed the mechanism for the generation of new detonation cell. It is also reported that bifurcations can contribute to the irregularity of the detonation front as they are a source of new triple points and they are seen occurring only in highly unstable detonations. Jiang et al. [13] studied the characteristics of cylindrical detonation from direct initiation of the detonation and observed similar characteristics. Vasil’ev at al. [14] studied the bifurcation process; however, his works were primarily focused on the two-scale cellular structures in a blend of mixtures. Although few studies on the bifurcation of the shock in the detonation front are reported, the literature available on this phenomenon is very limited.

Bifurcations caused by the instability at the detonation front as a source of new triple points or transverse waves are not much explored. Moreover, previous numerical simulations involved a simple single-step Arrhenius model. While the single-step kinetics can capture the essential features in the detonation, using detailed kinetics in the simulation can provide crucial detail of the reaction zone dynamics, especially in the case of marginal and highly unstable cases. Moreover, one-step kinetic models are known to underpredict the induction zone as it involves only exothermic reactions. With the advancement of computational power, it is possible for high-fidelity simulations of such complex flows using detailed reaction mechanisms. Moreover, multi-dimensional analyzing techniques have been developed in the past decade by which the instability in the multi-dimensional problem can be studied in detail. The instability at the detonation front is a multi-dimensional phenomenon, and understanding the critical characteristics involves solving multi-dimensional simulation. The scope of the present study is to numerically simulate and analyze the multi-dimensional detonation structure in the marginal detonation case with cell bifurcations and understand the associated instability at the detonation front using modal decomposition techniques.

2. Numerical Description and Implementation

2.1. Governing Equations

The propagation velocity behind a detonation wave is so high such that the effect of diffusion is negligible in the flow. Hence, for this study, a fully coupled implicit solver is used to solve the two dimensional, unsteady, compressible, inviscid, reacting flow in the rectangular computational domain. Coupled conservation equations are given below:

$$\frac{\partial Q}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}=W$$

(1)

$$Q=\left[\begin{array}{c}{\rho }_k\\ \rho u\\ \rho v\\ \rho e\end{array}\right],F=\left[\begin{array}{c}{\rho }_{k}u\\ \rho u^2+p\\ \rho uv\\ \left(\rho e+p\right)u\end{array}\right],F=\left[\begin{array}{c}{\rho }_{k}u\\ \rho u^2+p\\ \rho uv\\ \left(\rho e+p\right)u\end{array}\right],G=\left[\begin{array}{c}{\rho }_{k}v\\ \rho uv\\ \rho v^2+p\\ \left(\rho e+p\right)v\end{array}\right],W=\left[\begin{array}{c}{\omega }_k\\ 0\\ 0\\ 0\end{array}\right]$$

(2)

where $u$ and $v$ are the velocity components in the $x$ and $y$ direction. Total density $\rho$ is expressed as the sum of the partial density ${\rho }_k$ of each species while $k=1,2,\dots ,N_s$ where $N_s$ is the total number of species involved in the reaction.

$$\rho ={\sum }_{i=1}^{N_s}{\rho }_k$$

(3)

Pressure $p$ is evaluated from the ideal gas law for a thermally perfect mixture pressure and total energy per unit volume, $\rho e$ is calculated by adding the kinetic energy and internal energy.

$$p={\sum }_{i=1}^{N_s}\frac{{\rho }_k}{M_k}RT$$

(4)

$$\rho e=\frac{\rho }{2}\left(u^2+v^2\right)+{\sum }_{i=1}^{N_s}{\rho }_k\left({\int }^{{}^T}\frac{c_{pk}-R}{M_k}dT+h_k^0\right)$$

(5)

where $M_k$ is the molecular weight for the $k^{th}$ species, $R$ is the universal gas constant. The specific heats of each species are obtained as functions of temperature from NASA thermochemical polynomial data.

$$c_{pk}/R=a_{1k}+a_{2k}T+a_{3k}{T}^2+a_{4k}{T}^3+a_{5k}{T}^4$$

(6)

The mass production rate of each species is given as

$${\omega }_k=M_k{\sum }_{i=1}^{N_s}\left(v_{k,r}^{″}-v_{k,r}^{\prime }\right)\left[k_{fr}{\displaystyle {\prod }_{i=1}^{N_s}{\left(\frac{{\rho }_k}{M_k}\right)}^{v_{k,r}^{\prime }}}-k_{br}{\displaystyle {\prod }_{i=1}^{N_s}{\left(\frac{{\rho }_k}{M_k}\right)}^{v_{k,r}^{″}}}\right]$$

(7)

where $k_{fr}$ is the forward reaction rate constants calculated using Arrhenius equation from the reaction coefficients $A,b,E_a$ in the reaction mechanisms expressed as:

$$k_{fr}=A{T}^b\exp \left(-E_a/RT\right)$$

(8)

The backward reaction rate constant is calculated using the forward reaction rate and equilibrium constant which is calculated from the Gibbs free-energy minimum condition as follows:

$$k_{eq}^r=k_{fr}/k_{br}$$

(9)

$$k_{eq}^r={\left(\frac{1atm}{RT}\right)}^{{\sum }_{i=1}^{N_s}\left(v_{k,r}^{″}-v_{k,r}^{\prime }\right)}\exp \left[{\sum }_{i=1}^{N_s}\left(v_{k,r}^{″}-v_{k,r}^{\prime }\right)\left(\frac{S_k^0}{R}-\frac{H_k}{RT}\right)\right]$$

(10)

Here, entropy $S_k^o$ and enthalpy $H_k$ are obtained as a function of temperature from the specific heat data at standard rate.

2.2. Numerical Methods

In the present study, RPL2D in-house code for compressible reactive flow is used. The code has been developed and validated through a variety of applications for shock-induced combustion and supersonic combustion problems [15,16,17,18,19,20,21,22,23,24,25,26]. Though the RPL2D is a comprehensive code encompassing a number of mathematical models and numerical schemes mathematical models and schemes used in the present study are as follows. A coupled form of the species, momentum, and energy conservation equations are solved for the two-dimensional inviscid flow in rectangular geometries using fully implicit formulation. The finite volume cell vertex scheme is used for spatial discretization of the governing equations. The convective terms are computed using advection upwind splitting method difference vector (AUSMDV) scheme [27]. A third-order weighted essentially non-oscillatory (WENO) scheme is implemented to extrapolate the primitive variables at the cell interface to ensure accuracy around the discontinuity regions. A fully implicit second-order time integration method is used for the analysis. A Newton sub-iteration method is used to minimize the temporal discretization errors and to ensure the second-order time accuracy. Detailed hydrogen combustion model [28] from University of California, San Diego is used to calculate the chemical source terms which consists of eight reaction species (H, H₂, O, O₂, H₂O, OH, HO₂, and H₂O₂). A structured grid system is used for this analysis as shown in Figure 1. Equally spaced dense regular grids followed by elongated grids are employed for this use. Thus, the detonation front will be always within the dense regular cells. The exponentially elongated cells in the downstream of the tube helps the flow to reach CJ state within the domain and if not, CJ condition is enforced at the exit. In the Figure 1, equally spaced grids are used until 3 cm in the axial direction and beyond that the grids are exponentially expanded. For efficient computation in the multi-core SMP (share memory processors) machines, the code is parallelized by OpenMP method.

2.3. Detonation Initiation

For initiation, CJ detonation properties were calculated first using NASA CEA code [29] for the combustible mixture by fixing the initial temperature, pressure, and mole fraction. From the calculated CJ detonation velocity, steady-state ZND solution is obtained using CALTECH FORTRAN-based ZND code [30]. The steady-state solution is spatially non-linear, and hence an interpolation code was applied to map the non-linear solution over equally spaced grid points depending on the computational domain. The interpolated one-dimensional (1D) steady-state ZND solution is stacked along the two-dimensional (2D) computational domain in which an inclination to create an initial perturbation that is necessary for the detonation front to oscillate. The incoming boundary condition is fixed to that of the CJ detonation velocity so that the fully developed detonation front oscillates within the equally spaced grid points in the computational domain. Top and bottom walls are given periodic boundary conditions. A characteristic boundary condition, as described in [31] is used as an exit boundary condition. Exit boundary condition is necessary for detonation initiation until the detonation front is fully developed after which the detonation front is self-sustained and does not depend on the exit boundary condition. Smoked foil record is obtained by extracting the peak pressure behind the shock wave across the width of the computational domain and transformed from the wave attached computational coordinates to the laboratory coordinates as detailed in [31].

For the grid sensitivity of the solution, the spatial resolution depends on the detonative capability of the mixture. Ohira et al. [32] used 10–15 points in the half-reaction zone (HRZ) for analyzing highly unstable detonation case with reduced hydrocarbon reaction mechanisms. Inaba et al. [33] used a grid resolution of around 20–40 grid points in the HRZ. Tsuboi et al. [34] used 32 points in HRZ for 3D simulation. Ando et al. [35] used 50 computational points in the HRZ but used 32 points in the transverse axis. Detonation propagates in both the direction and predicting the transverse wave propagation is as important as predicting the frontal detonation oscillation in the highly unstable detonation. Moreover, the reactivity of the detonation depends on the combustible mixture property. Hence, a grid dependence study in a narrow channel for each combustible mixture in

Table 1. Summarizing all the simulations with initial conditions.

Detonation Stability	Fuel Composition	Initial Pressure (kPa)	Initial Temperature (K)	CJ Velocity (m/s)	Domain Size –Dense Regular Cells (L × W)	Mesh Resolution
Weakly unstable	2H2 + O2 + 17Ar	20	298	1415.0	3 × 6	601 + 100 × 1201
Moderately unstable	H2 + O2	10	298	2710.9	2 × 2	1001 + 100 × 1001

is conducted first to ensure that the detonation frontal structure is resolved properly.

2.4. Verification of the Approach with Weakly Unstable Detonation

Detonation propagation in a weakly unstable detonation, for which the detonation front oscillates regularly, is first analyzed to validate the implementation of the numerical scheme and the initiation approach for detonation studies. For this analysis, a weakly unstable detonation case with hydrogen and oxygen mixture in the stoichiometric mix ratio is considered. The mixture is diluted with 85% of argon at a pressure of about 20 kPa. The CJ velocity for this case, calculated from the NASA CEA is about 1413.5 m/s. The induction zone length for this case is around 0.1655 cm calculated from ZND steady-state solution using UCSD reaction mechanism. It was found that 15 points in the induction zone length calculated from the steady ZND solution is sufficient to calculate the dynamics properties of the detonation front. Figure 2 shows that with further increase in the spatial resolution, the Von Neumann pressure behind the detonation front does not change much. Thus, increasing the grid resolution does not affect the solution. The width of the channel is increased to 6 cm with the same grid resolution, which is around 0.01 cm. The instantaneous plot of the frontal detonation structure for this case is shown in Figure 3. From the density plots, the typical keystone features as observed at the detonation front can be clearly seen. The induction zone behind the weak incident shock can also be seen from the results. Shear layers formed behind the leading shock can also be clearly seen. The transverse waves moving perpendicular to the flow direction into the burned gas is recorded very well. Post shock pressure is extracted along the centerline of the computational domain and plotted in Figure 4. The regularity of the detonation front oscillation can be clearly seen from the graph. Detonation front is fully stabilized and steadily propagates in the standing detonation case. A sequence of snapshots tracking the detonation front oscillation is shown in Figure 5. Maximum pressure behind the shock wave is extracted from each iteration, and then plotted as numerical smoked foil record as shown in Figure 6.

3. Modal Decomposition

The dynamic mode decomposition (DMD) is a data processing algorithm that extracts coherent structures with a single temporal frequency from a numerical or experimental data-sequence. DMD is a modern and more useful tool to calculate cohesive structure in the flowfield. It enables us to observe the structures of the modes and to evaluate acoustic damping as a single representative value over the whole domain [36]. The data are flowfields saved over a time and for a fixed time interval (snapshots). These snapshots contain information on both linear and nonlinear structures in the flowfield. DMD enables the extraction of the main linear structures existing in a flowfield. The DMD also has the advantage to be field selective. The field selection can be, for instance, the choice of specific 2-dimensional planes in a 3-dimensional space where snapshots are recorded. Moreover, in each snapshot, a physical variable can be selected, for instance, pressure or temperature. The processed decomposed spatial modes refer to a pattern in the flowfield that may have significant impact on the overall flow and is referred to as “coherent structure”. Each coherent structure is mapped to a particular frequency and the temporal characteristics of the decomposed modes are analyzed from damping/growth factor [36]. This means that unstable modes can be predicted and extracted by the DMD method.

• A negative value of the damping coefficient indicates that the oscillation is dampening over a long period.
• A positive value indicates that the oscillation is growing over time and is unstable.
• For a stable oscillation, the damping coefficient must be near one.

The algorithm for the decomposition technique is clearly explained in [36]. The characteristics of each mode can be described by damping coefficient or growth rate, $\xi$ and coherence, $E$ and is given by:

$${\xi }_i=\frac{\ln \left|{\lambda }_i\right|}{N_S\Delta t}(s^{-1})$$

(11)

where ${\lambda }_i$ denotes the eigenvector of the projected matrix [A] after SVD operation and $N_s\Delta t$ represents the physical time between two consecutive snapshots.

Coherence can be used as a parameter to rank the modes as considered by Schmid [37] in his original DMD formulation. This is done by first calculating matrix $G$

$$G=V_1^N{\sum }^{-1}Y$$

(12)

where $V_1^{N-1}$ represents the first $N$ snapshots ${\sum }^{-1}$ is calculated from the SVD operation of the projected matrix $\left[A\right]$ and $Y$ represents the eigenmodes of the projected matrix. Finally, coherence is calculated from each column $g_i$ of the $G$

$$E_i={\Vert g_i\Vert }^{-1}$$

(13)

The value of coherence of each mode is normalized with the mean mode to analyze the contribution of each mode to the flowfield while the damping coefficient gives information on whether the mode is stable or not.

4. Detonation Cell Bifurcations

The instability at the detonation front leads to frequent decay of the weaker transverse waves. At some point during the propagation, the single large characteristics cells disintegrate into multiple small characteristics cells, and that process is referred to as the “bifurcation phenomena” in this study. At the characteristics scale, the cell disintegrates into only two cells, and hence it is referred to as bifurcation phenomena. Experimental and numerical results have shown these fine structures within the characteristics of the detonation cell for highly unstable detonation cases. For experimentally observed irregular cellular structures in hydrocarbon mixtures [7,8,9] the characteristic cells undergo cell bifurcation multiple times during the propagation. The weaker cells decay into the larger characteristics again and this problem continues further. A similar trend was seen with hydrogen oxygen mixture with nitrogen diluents. These substructures are observed in various specific mixtures such as in lean mixtures of hydrogen, methane, as well as in N₂O diluted mixtures. The process of how this bifurcation occurs is under-reported as most of the detonation occurs in hydrocarbon mixtures. For a detailed analysis of the bifurcation process, detonation front self-propagating in a hydrogen-oxygen mixture at low pressure is considered.

For the analysis, a detonation propagation in hydrogen oxygen mixture at stoichiometric ration at 10 kPa is considered which is reported to have fine structures within the characteristic detonation cells in the experiments [38]. In Figure 7, a grid refinement study is done with a narrow channel of 1 cm with four different resolutions of 100 (case 01), 50 (case 02), 25 (case 03), and 20 μm (case 04). The post shock VN pressure is quite low at low resolution and produces regular oscillations. However, at 20 and 25 μm, the peak VN pressure reaches maximum at the strong triple points compared to the low grid resolutions cases. The computational domain used for the further analysis is 2 cm for which 20 μm sized grid points were used, and the number of grid points in the steady-state induction zone is about 31. Time step for each iteration is fixed based on the CFL number 0.2, which is about 2.64 ns for this case. A total of 1.1 million grid points were used for this simulation in which around 1 million grid points are used in the equally spaced domain where the detonation front is confined.

Initially, as the flow is developing, only uniform cells are formed with two full characteristics cells and exhibit different strengths in the oppositely moving transverse waves. However, as the detonation front reaches around x = −25, the flow starts to undergo the bifurcation process. Numerical smoked foil record from the simulations is shown in Figure 8. Bifurcation phenomena is observed regularly during the propagation such as at around, x = −26, −45, −60, −90 in Figure 8. In order to avoid analyzing the bifurcations during the detonation initiation, third regular bifurcation during the propagation, at x = −60, is explored in this study, and further simulation reveals there exist two more bifurcations after bifurcations at x = −90. One mode of detonation cell bifurcation is considered for detailed analysis and the smoked foil record for the extracted window is shown in Figure 9. From the preliminary observation, it can be observed that the weak triple points, travelling from top to bottom of the channel, starts to decay into the strong triple points and further strengthens it. Even though it is weak, weaker triple points still exist at until around x = −55.5, around which all the weaker triple points completely merge with the primary triple point structure and become one full characteristics cell. One full mode detonation happens from this triple point and later at around x = −58.75, the strong transverse wave of this full detonation cell interacts. Even at this time, it can be clearly seen that at the time of interaction, the strength of the transverse waves of the opposite family is different as observed after x = −58. After this strong interaction, bifurcation begins to occur at around x = −59.5, x = −61 and forms new additional cells which are stronger than the initially observed weak triple points. Detonation front oscillation is confined around 0.5 cm from the inlet boundary, and hence the induction zone and reaction zone properties are well resolved within the uniform space grid points.

From the preliminary investigation of the numerical smoked foil record, it can be observed that all the weaker triple points accelerate into the stronger transverse waves. The acceleration of the weaker triple points such as at x = −38, −52, and −81 results in a decaying of the weaker triple points into the stronger triple points and finally only single mode of detonation spanning the width is observed. They are highly unstable, and further propagation of the detonation front results in the formation of new additional detonation cells with additional triple points which are stronger than the previous weaker triple points. Moreover, in the single mode of detonation, during the decaying process, strong triple points are observed in all instances which interact and then break into multiple cells. Moreover, before bifurcation, the strength of the transverse waves of the opposite family are different and strong triple points can be seen moving in one direction while the other family has weaker triple points during the propagation. The direction of the movement of the strong triple points also changes during the propagation, that is, at x = −35 and −66 stronger triple points can be seen moving from the bottom to the top. At x = −50 and −95, the stronger triple points are seen moving from top to bottom. This shows that the direction of the transverse waves is defined after the breaking of the single-mode of detonation. Finally, the unstable single mode of detonation breaks into multiple cells (mostly two cells except at around x = −48). Moreover, the unstable mode of detonation runs for two complete cycles before breaking into multiple cells. Thus, the bifurcation process can be divided into three modes as:

• Decaying of the weaker triple points;
• Formation and decaying of a new single-mode of detonation;
• Bifurcation of the single-mode of detonation.

The mechanism in which this process is studied from the detailed analysis of the frontal detonation structure is explained as follows.

4.1. Decay of the Weaker Triple points

Evolution of the detonation front for x = −49.5 to x = −52.875 is shown in Figure 10. From the density gradient graph, it can be clearly seen that the shock wave and the reaction zone are decoupled for the transverse wave family moving from the bottom of the channel to the top of the channel. OH mass fraction figures clearly show the slow reaction behind the decoupled shock and fast reactions behind the coupled shock and reaction zones. These slow reactions interact with the transverse waves of the opposite family, which moves fast and produces additional weaker triple points. These additional triple points are seen in the numerical smoke foil record in Figure 9. These slow reactions create hot pockets of unburned gas and the transverse waves react with those pockets and gets strengthened again. The reaction front is coupled with the oppositely moving transverse waves of this detonation at first, however, because of the weaker transverse waves, they slowly decouple and weaken the triple point which moves from bottom to top. This decoupling of the shock wave and the reaction front result in weaker incident shock which starts to decay; whereas, on the other side, the transverse wave propagates much strongly by the coupling of shock wave and reaction zone.

4.2. Single-Mode Detonation

Transverse waves moving in the perpendicular direction to the flow direction from top to bottom are strong enough to consume any unburnt fuel left as observed from the radical mass fraction profiles. So, in the further propagation of the detonation front, as shown in Figure 11, the detonation results in the formation of a single cell. At somewhere in the evolution, the shock wave and the reaction zone are coupled again as can be observed in the same figure. This is because of the fast-moving transverse waves which react with the decoupled reaction zone. When all the weaker triple points are decayed, only one full characteristic detonation cell is formed in the system which produces enormously long induction zone behind the weak Mach stem. When the transverse wave from this single mode of detonation interacts, it results in a mild explosion as seen in the figure. Additional triple points are observed, but they are decoupled with the reaction zone and do not sustain. The fast reactions occur near the triple point interactions as observed from the HO₂ mass fraction as found in the figures in the next Section 4.3. Moreover, the oppositely moving waves exhibit different characteristics, even at this point of time, as seen from the density gradient plots.

This single mode of detonation continues again, but the cell length becomes shorter compared to the last cell which clearly shows that the detonation front is accelerating in the second single mode of detonation. Compared to the first single-mode, the triple point is very closely attached to the detonation front as the detonation front is accelerating. When these triple points interact, they produce strong reactions which push the detonation front further. The oppositely moving transverse shows different characteristics again even after the single mode of detonation, which results in acceleration of the slow reactions in the direction of the strong transverse waves. This acceleration moves closer to the detonation front and creates a new weak triple point. This decoupled new weaker triple point produces a new unburned pocket of the combustible mixture which reacts with the transverse waves.

4.3. Cell Bifurcation

Finally, the oppositely moving transverse waves with different strengths result in the combustion of the unburnt pocket of reactants, as observed in Figure 12. Moreover, the transverse waves moving from the top to the bottom of the channel get stronger because of these unburnt fuels. At some point, this results in a micro explosion, as seen in the early stages of Figure 13, which is very close to the Mach stem. The hot unburned pocket of gas explodes instantaneously as it interacts with the transverse waves and creates a new triple point configuration. This new triple point gets strengthened by the micro explosion that happens from the hot unburned reactions and forms a new cell during the propagation. It should be noted that there are other unburnt pockets of fuel observed before this stage, but they get burned in the fast-moving transverse waves.

This micro explosion creates shock bifurcation at the detonation front and produces a new set of triple points. This new set of triple points slowly evolves by interacting with the other set of triple points which already exists and finally stabilizes after producing new characteristics cells. Finally, stronger triple points of two different families are observed at around x = −61 in the numerical smoked foil record as shown in Figure 9.

For understanding the reaction progress at the frontal detonation structure, during this propagation, mass fraction distribution of the radicals is analyzed. As seen in Figure 14, Figure 15 and Figure 16, initially, multiple keystones such as feature are seen behind the leading shock wave where the reaction zone is decoupled from the shock wave. On the other side, shock wave and the reaction zone are attached in the other leading shock wave which is moving toward this detached induction zone. HO₂ mass fraction plots clearly show that these oppositely moving transverse waves have different strengths as they produce stronger reactions in one set of the triple point of interactions. Here, it can be seen that slow reactions behind the decoupled reaction wave accelerate and move into the incoming transverse waves. The decoupled and slow reactions can be seen clearly from the temperature and H mass fraction plots.

This strengthens the transverse wave which is moving toward the detached zone as seen in snap (b) where the transverse wave, moving toward the bottom wall is almost straight normal to the direction of the propagation of the detonation front. When this strengthened transverse wave crosses the hot compressed gas, it reacts and forms a new detonation front travelling in the transverse direction but quenches as soon as the reactants are all consumed. Since this transverse wave is moving faster, this could not support the newly formed leading shock seen on the top side of the snap (b). Ultimately it results in a detached shock wave and begins to quench. Another fast-moving transverse wave is strengthened because of this detachment as seen in snap (c). This process can be clearly understood by seeing the species profile concentration, as shown in Figure 14, Figure 15 and Figure 16. From the H radical which signifies the initiation of the chain branching reactions, the radicals are produced and consumed within a short distance. This is again confirmed from the OH profile in Figure 15. where the species concentration is high during the energy release. The leading shock wave formed from this detonation reacts very fast, as seen in snap (c), where the other detonation cell in which the leading shock wave detaches from the reaction wave again. This results in the formation of a new stronger detonation front and a new weaker detonation front, as seen in snaps (c) and (e). Over a period of time, the weaker detonation front completely detaches itself from the incident shock wave, which results in a large induction zone as seen in snap (h).

5. Modal Decomposition of the Bifurcation Process

This oscillation at the detonation front is a multi-dimensional instability problem. Hence the modal decomposition technique was applied to investigate this instability phenomenon. First, the decomposition analysis was applied to different variables such as temperature, pressure, density, velocities to identify an ideal variable, and density is rightly suited for explaining the instability behavior. Decomposed results of the weakly unstable case predict a regular fluctuation of the detonation front with a primary mode oscillating at a frequency of around 28 kHz. Primary structure of the decomposed mode is shown in Figure 17. It can be clearly seen that the extracted coherent structure clearly shows the interaction of the transverse waves, and the frequency for this mode is around 28 kHz.

Moderately unstable detonation case oscillates at a higher frequency and has more irregularity. At first, the POD and DMD spatial structures were compared as this process is not an oscillating phenomenon as observed with earlier shock-induced combustion case or regularly oscillating weakly unstable detonation case. Figure 18 shows the spatial structures of the coherent modes of POD and DMD, and they resemble the same with little differences in the secondary mode of the spatial structure. The difference is because the DMD method extracts the decomposed mode and maps it based on their temporal characteristics, the POD method extracts the coherent structure based on its eigenvalue. The primary (top left) and secondary POD (top right) modes with energy content, E of around 2.86 and 2.24 were mapped to the corresponding DMD modes which have a frequency of about 90.6 kHz (bottom left) and a low-frequency instability which is approximately 23.4 kHz (bottom right), respectively. The ensemble POD mode has E of around 22, which indicates that the flowfield is highly unstable. From the decomposed results, it becomes evident that the compression occurs at two different cycles in the dominant mode. The first cycle is the compression in the detonation front, and another cycle is the compression happening near the triple point region. The fluctuation becomes maximum in the region where the interaction of the oppositely moving transverse waves is strong. Similarly, the structure of the secondary decomposed mode has two different cycles. The low-frequency mode shows that the triple point where the Mach stem interacts with the incident wave oscillates along the plane normal to the detonation propagation.

The cause for this instability cannot be analyzed clearly as the waves are decaying and forming again at high frequency, and the decomposed result of high-frequency modes shows low temporal characteristics. Analyzing the full domain is possible; however, analyzing the results is challenging. Hence, the domain is divided into four different parts of the analysis. The first part consists of the decaying of the subscale detonations into single large-scale detonation. The second and third part contains the dynamics of one full single detonation cell in each part, and the fourth part includes the dynamic information of the creation of subscale detonations which results in the shock bifurcation in the detonation front.

The primary DMD modes of the divided parts, as shown in Figure 19, have the frequency of 73.39 kHz, 89.28 kHz, 71.64 kHz, and 90.54 kHz for parts 1–4 respectively. It is to be noted that the oscillating frequency increases initially when the weak triple point in the detonation front decays but decreases from part 2 to 3 within the single mode of detonation. It again increases from part 3 to 4, where multiple transverse waves are created too. Moreover, the location of the triple point during the interaction of the Mach stem and the incident shock varies. This is because the transverse waves properties propagating in the opposite direction are different, which influences the detonation cell structure. The structure of the primary mode in part 2 shows a strong fluctuation along the plane normal to the detonation front than in part 1. This indicates that the weak transverse waves are getting stronger, and the interaction of the strengthened transverse wave creates an impact in that region. As observed from the dominant mode in Part 4, there exist a secondary triple behind the primary triple point with a longer reflected shock along with the triple point in one cycle, and the other oscillating process has weaker triple.

The secondary DMD mode of the analysis shows a different trend in the variation of oscillation frequency, where the operating frequency is increased three times compared to the primary mode. The frequencies of the secondary modes are 150.94 kHz, 164.74 kHz, 208.93 kHz, and 151.64 kHz in the order. By analyzing part 1 result from Figure 20, the weaker triple points in the leading detonation front moving into the stronger transverse are recorded as one cycle, and the interaction of the transverse waves are recorded as the other cycle. The coherent structure of part 2 and 3, for which one full detonation cell is observed in each part, shows different features. The fluctuation in the cohesive structure of part 3 shows that the triple points oscillate primarily in the transverse direction and it becomes weak in the axial direction when compared to part 2. By observing part 4, the fluctuations in the triple point again start oscillating in both axial and transverse directions.

By analyzing the temporal characteristics of these modes in Figure 21, coherence of primary mode for parts 1–3 is high, but their growth rate indicates they are not stable. Part 4 is comparatively stable within the analyzed window compared to the other parts, and the coherence value is low compared to the other three parts. Low coherence value of the primary mode indicates that the flowfield is dynamically stable. For part 3, the secondary mode is operating at a higher frequency than the other parts. As mentioned before, parts 2 and 3 contain single-mode detonation each within the domain. So, the secondary mode oscillation in the single-mode detonation grows during this process and leads to the formation of new weaker detonation cells in part 4. The detonation stays relatively stable when the new cells are formed because of the shock bifurcation. However, the secondary mode is highly unstable and affects the stable propagation and causes instability again.

From this analysis, it can be concurred that the weaker triple points decay into the stronger triple points at a higher frequency and increases the fluctuation of the frontal detonation structure from part 1 to 2. During this process, the transverse wave properties are strengthened and create a stronger triple point during their interaction. Because of the stronger triple point, the frontal detonation structure is pushed farther, and this causes the frequency of the single-mode of detonation to decay in part 3. At this juncture, the interaction of triple point and the heated combustible mixture behind the weak incident shock increases the chemical sensitivity and propagates at high frequency along the transverse direction around the triple point region. This leaves a small pocket of the unreacted heated combustible mixture which interacts with the reflected shock waves and detonates. The shock wave from this secondary detonation interacts with the leading detonation front and initiates the bifurcation process, which results in multiple weak triple points in the detonation front. Those triple points in the detonation front slowly develop and increase the oscillation frequency again in part 4. Thus, multiple weaker triple points were formed and decayed in the moderately unstable detonation with shock bifurcation, oscillates the frontal detonation structure, and cause instability phenomena. The rate at which this process occurs determines the propagation characteristics of the detonation.

6. Conclusions

The dynamics of the bifurcation phenomena was studied for a marginal detonation case with detailed kinetics. For marginal detonation, the local oscillation of the detonation front is high resulting in frequent formation and decay of additional characteristic cells during the propagation. Combustion of the secondary explosion in the unreacted pockets causes the leading shock front to bifurcate, thereby creating new set of triple points and the characteristics cells. The newly formed cells gradually weaken and results in a single mode of detonation. The local oscillation in 2D detonation front is studied, and the bifurcation process is analyzed further using dynamic mode decomposition technique and it clearly shows the location of the secondary explosion. From the temporal characteristics of the decomposed modes, it is evident that the bifurcation phenomena are triggered by the secondary mode of oscillation that gets strengthened during the single mode of detonation. The computational cost associated with the current numerical simulation hampers the numerical simulation of such cases with detailed kinetics in 3D flowfield. For this case, 16,500 CPU hours are spent distributed on a 24 core SMP machine for 1.1 million grid points. Real detonations have three dimensional effects which may affect the frequency of the instability modes. The authors are confident that while the solution may be affected with 3D solution, such as the frequency of instability mechanism differing from the 2D solution, the overall bifurcation mechanism will not be affected.