# Role of the turbulent burning in the quenching limits of gaseous detonation confined by an inert compressible layer

Hiroaki Watanabe, Vincent Robin, Ashwin Chinnayya

► To cite this version:

Hiroaki Watanabe, Vincent Robin, Ashwin Chinnayya. Role of the turbulent burning in the quenching limits of gaseous detonation confined by an inert compressible layer. 30th International Colloquium on the Dynamics of Explosions and Reactive Systems (ICDERS 2025), Jul 2025, Ottawa, Canada. hal-05410515

HAL Id: hal-05410515

https://hal.science/hal-05410515v1

Submitted on 17 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Role of the turbulent burning in the quenching limits of gaseous detonation confined by an inert compressible layer

Hiroaki Watanabe, Vincent Robin, Ashwin Chinnayya

Institut Pprime UPR 3346 CNRS, ISAE-ENSMA, Université de Poitiers

1 avenue Clément Ader, BP 40109, 86961 Futuroscope-Chasseneuil CEDEX, France

## 1 Introduction

In the classical model of Zel'dovich, von Neumman and Döring (ZND model), a detonation has a steady and one-dimensional (1D) structure [1] and adiabatic shock compression is the mechanism for ignition. On the other hand, detonation exhibits unsteady multi-dimensional cellular structure. Nevertheless, previous studies revealed that the adiabatic shock compression of leading shock and transverse wave are the main factors in ignition in a weakly unstable mixture with a regular cellular structure [2–4]. However, in an unstable mixture with an irregular cellular structure or especially near limits, the adiabatic shock compression is not sufficient to ignite the unburned gas mixture, with the formation of unburned gas pockets behind the leading fronts [2, 5–7]. The latter are consumed by turbulent burning, which is not accounted for in the 1D steady ZND model and is not captured in inviscid simulations. Therefore, simulations without proper modeling of turbulent burning failed to reproduce the dynamic behavior of detonation near limits and predicted a different trend from experiments. Indeed, numerical cellular detonation is more difficult to propagate than ZND prediction or 1D simulations, and a weakly unstable mixture is easier to detonate than an unstable mixture, which is highlighted as a detonation paradox [6].

In response to this situation, Maxwell [8] formulated the compressible linear eddy model for Large Eddy Simulations (CLEM-LES) as a potential solution for modeling compressible reactive flows. Indeed, CLEM-LES has been successfully applied to transition of a turbulent shock-flame complex to detonation [10], detonation propagation in a narrow channel [9, 11], and detonation propagation confined by an inert layer [12]. While the formulation used in CLEM-LES offered a sophisticated approach, its complexity posed challenges in isolating and quantifying the role of turbulent burning in detonation. This underscores the need for a simple and predictive model to enhance simulation accuracy and facilitate its broader adoption.

The semi-confined detonation gathered attention from point of view of mimicking condensed detonation [13], rotating detonation engine [14], and safety engineering [15]. In this configuration, the expansion of the burned gas causes the streamline to diverge with velocity deficit, and quenching occurs when the height of the reactive layer becomes less than the critical value. Therefore, the measurement and prediction of critical conditions were of paramount importance. Nevertheless, Taileb et al. [16] have shown that chemical modeling gave a significant influence on the critical height from simulations in 2H₂ – O₂ mixture. Although critical height using simplified one-step and three-step models was about 4-5 times higher than the experimental value, simulations with detailed chemistry predicted a closer value to experiment [13] (6 mm in simulation and 4.6 mm in experiment). Thanks to recent studies [17, 18], predictive ability of critical height from simulations using simplified chemistry was improved.

Although Euler simulations with a proper chemical modeling gave closer results for critical height [16, 18], there is still a gap between experimental and numerical critical height. One of the factors most likely
contributing to the observed differences is the inability to accurately capture reactive fronts whose propagation is driven by thermal/molecular diffusion, as well as turbulent diffusion, as observed in turbulent flames. This mechanism is expected to play a crucial role in the ignition of unburned mixtures near detonation limits or in unstable mixtures [6].

The present study aims to propose a simple approach to model a deflagrative burning behind detonation front toward predictive simulations and discuss the relevance of deflagrative burning in detonation near limiting conditions. Simulations were conducted for a stoichiometric hydrogen-oxygen mixture, using an additional modeled source term to account for the consumption of unburned gas pockets, activated only below the crossover temperature. This source term was inspired by the propagative form of the mass fraction transport equation, which is similar to a level-set formulation [19, 20].

## 2 Numerical setup and problem statement

The chemical reaction was modeled by the modified three-step chemical model for ideal detonation (3SMI) proposed by Watanabe et al. [18]. Four chemical species are taken into account: reactants F, radicals R, products P, and inert gas I. The conservation law for the species mass fraction $Y_k$, $k \in \{F, R, P, I\}$ is

$$ \frac{\partial \rho Y_k}{\partial t} + \frac{\partial \rho u_i Y_k}{\partial x_i} = S_k. \qquad (1) $$

To account for flame propagation in unburned pockets, the right-hand side (RHS) term $S_k$ must include contributions from both diffusion and chemical reaction processes:

$$ S_k = \mathcal{F}_k + \dot{\omega}_k, \qquad (2) $$

where $\mathcal{F}_k$ and $\dot{\omega}_k$ represent the diffusion flux and the chemical reaction rate of species k, respectively. Nevertheless, to simplify the modeling, we introduce a filter $H$ to distinguish between regions where auto-ignition mechanisms predominantly occur ($H = 0$) and regions where flame constant-pressure propagation mechanisms predominantly occur ($H = 1$). Thus, the RHS can be expressed as:

$$ S_k = (1 - H)S_k|_{H=0} + H S_k|_{H=1}, \qquad (3) $$

where $S_k|_{H=0}$ and $S_k|_{H=1}$ represent the RHS terms conditioned for the auto-ignition combustion zone and the constant-pressure combustion zone, respectively.

In the auto-ignition combustion zone, diffusion fluxes can be neglected, so the RHS reduces to:

$$ S_k|_{H=0} = \dot{\omega}_k. \qquad (4) $$

However, in the constant-pressure combustion zone, diffusion processes cannot be neglected, and the flame structure induced may be too thin to be accurately captured on the mesh. Thus, in these regions, the computed variables represent filtered averages, similar to the filtering used in Large Eddy Simulations. For simplicity, we use the same notation as before, without explicitly indicating spatial filtering. Thus, in these zones, the diffusive fluxes include turbulent diffusion, and the chemical reaction rate requires a modeling approach. To preserve the Euler formulation of the two-dimensional (2D) compressible reactive equations considered, the RHS is expressed in its propagative form, which incorporates both diffusive and reactive effects. Finally, in this zone, the RHS for each species is written as:

$$ S_F|_{H=1} = -\rho_u S_u \Xi |\nabla Y_F|, \quad S_R|_{H=1} = 0, \qquad (5) $$

$$ S_P|_{H=1} = \rho_u S_u \Xi |\nabla Y_F|, \quad S_I|_{H=1} = 0. \qquad (6) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Mass flux for 1D steady laminar flame based on initial condition (a) along Hugoniot curve, and (b) behind one transverse wave.</figcaption></figure>

Here, $\rho$, $\rho_u$, $S_u$, and $\Xi$ are density, density for unburned gas, burning velocity for unburned gas, and the flame wrinkling factor, respectively. The mass flux through the laminar flame structure $\rho_u S_u$ and the flame wrinkling factor are input parameters and treated as constants. The filter $H$ is defined as a Heaviside step function based on the cell temperature $T$, the crossover temperature for chain-branching reactions $T_B$, the pressure $p$, and the initial pressure $p_0$, as follows:

$$ \begin{cases} H = 0 & (T \ge T_B \text{ or } p < 2p_0) \\ H = 1 & (T < T_B \text{ and } p \ge 2p_0) \end{cases} \qquad (7) $$

In addition, the threshold using pressure and mass fraction of reactant was included to prevent the deflagrative burning ahead of the detonation front. Near the edges of the unburned gas pockets, there can be some overlap of both mechanisms during the course of the overall burning process, specifically from deflagration to auto-ignition. The extension of modeling to other mixture without crossover temperature is a potential avenue for future research. The governing equations were solved using our parallel in-house code RESIDENT (REcycling mesh Simulation of DEtoNaTions), with the details given in Refs. [16, 18, 21]. The characteristic variables on cell boundaries were reconstructed using a ninth order MP-R interpolation [22]. The numerical flux were computed by the HLLCM solver [23]. The second order numerical Hamiltonian was evaluated [24] and the source term for deflagrative burning was integrated by second order two stage Strong Stability Preserving Runge-Kutta method. The model parameters in 3SMI were the same as those used in Watanabe et al. [18].

The target mixture was 2 H₂ – O₂ at ambient condition (101 kPa, 295 K). The detonation propagation in reactive mixture bounded by N₂ was simulated. The simulation procedure and conditions were the same as in previous studies [16, 18]. The grid resolution was 10 points per ZND induction length ($x_{ind}$), which was sufficient to capture the critical height and the flow fields [16, 25]. The height of the reactive layer $h$ and the product of mass flux and flame wrinkling factor $\rho_u S_u \Xi$ are parameters in the parametric studies. Each case was run only one time. Supercritical condition was defined as the successful transmission and propagation of detonation over 3000 $x_{ind}$.

The reference data from 1D steady laminar flame were obtained using the detailed chemistry (Mevel2017), with 9 species and 21 elementary reactions [26].

## 3 Results and discussions

In the source term for the deflagrative burning, the mass flux $\rho_u S_u \Xi$ was considered as a constant parameter in the present modelling (Eq. 5). The value of mass flux was estimated from 1D steady laminar flame using relevant initial conditions for detonation, i.e, along Hugoniot curve or behind a transverse wave in Fig. 1. Note that the unburned gas pockets were formed from the gas shocked by weak incident shock around the end of the cell cycle whose induction time was much longer than the
<figure><img src="image_3.png" /></figure>

Figure 2: Maximum pressure history for (a) $h = 5$ mm and $\rho_u S_u \Xi = 156.3 \text{ kg/(m}^2\text{s)}$, (b) $h = 5$ mm and $\rho_u S_u \Xi = 109.4 \text{ kg/(m}^2\text{s)}$, and (c) $h = 4$ mm and $\rho_u S_u \Xi = 156.3 \text{ kg/(m}^2\text{s)}$. White lines are shock front shape and the interval between the two successive lines are 1.0 $\tau_{cell}$. Here, the characteristic time for detonation cell $\tau_{cell}$ is defined as $\tau_{cell} = L_{cell}/D_{CJ}$, where $L_{cell}$ and $D_{CJ}$ are average cell length for ideal detonation in simulation and CJ velocity, respectively. Red dotted line is the initial position of the interface between the inert gas and reactive mixture. Note that the horizontal and vertical scales are different, to ease visualization of the cellular structure.

time for cellular cycle. The mass flux at the post-shock condition did not change by order of magnitude, and the change in the value between the post shock state with 0.60 $D_{CJ}$ and 0.80 $D_{CJ}$ was within the factor of 2.5 (Fig. 1(a)). In the second part of the cell, the unburned gas experienced a transverse wave after the passage of the leading shock [4]. The value of mass flux for various second shock Mach number is shown in Fig. 1(b). As the second shock Mach number increased, the mass flux also increased. The ratio of the increase in the mass flux with a second shock Mach number of 1.4 in the post shock condition with 0.60 $D_{CJ}$ was about 2.3. In addition, the increase ratio became lower as the leading shock velocity was lower. Therefore, the present modeling of constant mass flux for the deflagrative burning can be acceptable as an approximation especially for lower leading shock velocity near the limit.

The maximum pressure history for the selected cases is shown in Fig. 2. The case with $h = 5$ mm and $\rho_u S_u \Xi = 156.3 \text{ kg/(m}^2\text{s)}$ was a supercritical condition and the average propagation velocity was 0.870 $D_{CJ}$ (Fig. 2(a)). The failure and reinitiation by transverse detonation were repeated and the transverse detonation generated new cells enough. Nevertheless, the number of the triple points was relatively constant and the propagation was stable. Note that without deflagrative burning model, detonation can not keep propagation at this reactive layer height. With the decrease in the mass flux in $h = 5$ mm and $\rho_u S_u \Xi = 109.4 \text{ kg/(m}^2\text{s)}$, the consumption speed for unburned gas pockets became lower and the period of failure became longer (Fig. 2(b)). The divergence of streamlines was so strong and the deflagrative burning with this value of mass flux was not sufficient to burn the unburned gas pockets and sustain the detonation propagation. The cell width drastically increased and quenching was observed about 50 mm after interaction with an inert layer. With the shorter reactive layer height and the same mass flux as used in the supercritical condition in Fig. 2(a), the subcritical condition was observed (Fig. 2(c)). The reinitiation by transverse detonation mainly originated around the interface between the inert and reactive layer, which was in line with the observation on reinitiation in transient phase by Murray and Lee [27]. Nevertheless, the increased lateral expansion caused by a shorter reactive layer height led to the quenching.
<figure><img src="image_4.png" /><figcaption>Figure 3: Instantaneous flow fields for (a) temperature, (b) mass fraction of reactant, with $h = 5$ mm and $\rho_u S_u \Xi = 156.3 \text{ kg/(m}^2\text{s)}$, and (c) summary of results in the coordinate system of mass flux and height of reactive layer.</figcaption></figure>

Figure 3(a,b) depicts the 2D instantaneous flow fields in the supercritical condition for the case with $h = 5$ mm and $\rho_u S_u \Xi = 156.3 \text{ kg/(m}^2\text{s)}$. The classical detonation-shock combined wave structure consisting of detonation front and oblique shock was observed (Fig. 3(a)). The leading shock front was curved due to the divergence of streamlines. With the inclusion of the source term for deflagrative burning, the unburned gas pockets were burned from their edges and the energy release from the consumption of unburned gas pockets within the detonation driving zone further supported the detonation propagation (Fig. 3(b)).

The simulation results are summarized in the coordinates of mass flux $\rho_u S_u \Xi$ and reactive layer height $h$ in Fig. 3(c). As the mass flux increased, the average propagation velocity increased and the propagation became more stable. Moreover, the order of the velocity deficit in simulations (~ about 13%) agreed with the experimental observation [28] (~ about 16%), whereas simulations without deflagrative burning model predicted the order of about 5% as velocity deficit in the same level of unstable mixture [21]. In addition, the velocity deficit in simulations was higher than the velocity deficit at maximum curvature from steady 1D ZND model (~ 9%) [18], which highlighted that the deflagration extended the propagation limit compared to steady 1D ZND model and was in line with experimental trend [7]. Present simulations with a source term for deflagrative burning were able to capture the experimental critical height in the range of the mass flux used.

The ratio $\beta = (\rho_u S_u \Xi)/(\rho_u S_l)_{\text{ref}}$, which could reproduce the experimental quenching limit was between 1.09 and 1.86. The value $(\rho_u S_l)_{\text{ref}}$ was estimated in this study at post-shock state with state with a velocity of $0.6 D_{\text{CJ}}$, $S_l$ being the laminar flame speed. The value of $0.6 D_{\text{CJ}}$ was taken as representative of shock velocities in marginal conditions, around the end of the cell. This range of values was found to be lower to that reported by Maxwell et al. [9] ($\beta \approx 6.6-7.3$ from experiments, $\beta \approx 3.63-3.74$ from simulations) and Xiao et al. [7] ($\beta \approx 2.4$ to $6.5$ from experiments), probably due to their use of unstable mixtures, sensitivity in the selection of the reference value for normalization, and also in the determination of the deflagration surface. In the present study, this value was found to be closer to the laminar value, as suggested in ref [29]. This observation is consistent with the fact that developed turbulence occurred primarily near the end of the hydrodynamic thickness [30]. However, this ratio relies on several modeling assumptions and should be interpreted with care.
## 4 Conclusions

The inclusion of the source term for deflagrative burning significantly extended the propagation limit and successfully captured marginal detonation propagation. Furthermore, the experimental critical height was accurately reproduced. The mass flux for deflagration to recover the experimental critical height ranged from 109.4 kg/(m²s) to 187.5 kg/(m²s), which corresponds to 1.09 to 1.86 times the laminar value based on the post-shock condition with 0.60 $D_{cJ}$. The present results emphasized the pivotal role of deflagration in reproducing the dynamics of detonation near the limits.

## Acknowledgment

The computations were performed using HPC resources from GENCI-CINES (Grant A0172B07735). This work pertains to the French Government program "Investissements d'Avenir" (EUR INTREE, reference ANR-18-EURE-0010). Part of this research was supported by the French National Research Agency (ANR) under the project TRACKDEMO ANR-23-CE05-0021-010.

## References

[1] Lee JHS. (2008) The detonation phenomenon. Cambridge University Press.

[2] Radulescu MI, Lee JHS. (2002) The failure mechanism of gaseous detonations: experiments in porous wall tubes. Combust. Flame 131: 29.

[3] Xiao Q, Radulescu MI. (2020a) Dynamics of hydrogen-oxygen-argon cellular detonations with a constant mean lateral strain rate. Combust. Flame 215: 437.

[4] Watanabe H, Matsuo A, Chinnayya A, Itouyama N, Matsuoka K, Kasahara J. (2025) Lagrangian characterization of induction and reaction timescales in a cellular gaseous detonation. Phys. Fluids. 36: 026101.

[5] Radulescu MI, Sharpe GJ, Lee JHS, Kiyanda CB, Higgins AJ, Hanson RK. (2005) The ignition mechanism in irregular structure gaseous detonations. Proc. Combust. Inst. 30: 1859.

[6] Radulescu MI. (2018) A detonation paradox: why inviscid detonation simulations predict the incorrect trend for the role of instability in gaseous cellular detonations?. Combust. Flame 195: 151.

[7] Xiao Q, Radulescu MI. (2020b) Role of instability on the limits of laterally strained detonation waves. Combust. Flame 220: 410.

[8] Maxwell BM. (2016) Turbulent combustion modelling of fast-flames and detonations using compressible LEM-LES. PhD thesis, University of Ottawa.

[9] Maxwell BM, Bhattacharjee RR, Lau-Chapdelaine SSM, Falle SAEG, Sharpe GJ, Radulescu MI. (2017) Influence of turbulent fluctuations on detonation propagation. J. Fluid Mech. 818: 646.

[10] Maxwell BM, Pekalski A. Radulescu MI. (2018) Modelling of the transition of a turbulent shock-flame complex to detonation using the linear eddy model. Combust. Flame 192: 340.

[11] Maxwell BM, Wang WH. (2023) The influence of boundary conditions on three-dimensional Large Eddy Simulations of calorically perfect gas detonations. Flow Turbul. Combust. 111: 1279.

[12] Maxwell BM, Melguizo-Gavilanes J. (2022) Origins of instabilities in turbulent mixing layers behind detonation propagation into reactive-inert gas interfaces. Phys. Fluids 34: 106107.

[13] Dabora EK, Nicholis JA, Morrison RB. (1965) The influence of a compressible boundary on the propagation of gaseous detonations. Proc. Combust. Inst. 10: 817.

[14] Bykovskii FA, Zhdan SA, Vedernikov EF. (2006) Continuous spin detonations. J. Propul. Power 22: 1204.

[15] Rudy W, Kuznetsov M, Porowski R, Teodorczyk A, Grune J, Sempert K. (2013) Critical conditions of hydrogen-air detonation in partially confined geometry. Proc. Combust. Inst. 34: 1965.

[16] Taileb S, Melguizo-Gavilanes J, Chinnayya A. (2020) Influence of the chemical modeling on the quenching limits of gaseous detonation waves confined by an inert layer. Combust. Flame 218: 247.

[17] Veiga-López F, Taileb S, Chinnayya A, Melguizo-Gavilanes J. (2024) Towards predictive simplified chemical kinetics for hydrogen detonation. Combust. Flame 269: 113710.

[18] Watanabe H, Taileb S, Veiga-López F, Melguizo-Gavilanes J, Chinnayya A. (2024) A simple predictive three-step chemical model for gaseous stoichiometric hydrogen-oxygen detonation quenching. Combust. Flame 268: 113609.

[19] Pitsch H. (2005) A consistent level set formulation for Large-Eddy Simulation of premixed turbulent combustion. Combust. Flame 143: 587.

[20] Moureau V, Fiorina B, Pitsch H. (2009) A level set formulation for premixed combustion LES considering the turbulent flame structure. Combust. Flame 156: 801.

[21] Reynaud M, Virot F, Chinnayya A. (2017) A computational study of the interaction of gaseous detonations with a compressible layer. Phys. Fluids 29: 056101.

[22] He Z, Zhang Y, Gao F, Li X, Tian B. (2016) An improved accurate monotonicity-preserving scheme for the Euler equations. Comput. Fluids 140: 1.

[23] Shen Z, Yan W, Yuan G. (2016) A robust HLLC-type Riemann solver for strong shock. J. Comput. Phys. 309: 185.

[24] Osher S, Sethian JA. (1988) Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys. 79: 12.

[25] Mi XC, Higgins AJ, Kiyanda CB, Ng HD, Nikiforakis N. (2018) Effect of spatial inhomogeneities on detonation propagation with yielding confinement. Shock Waves 28: 993.

[26] Mével R, Javoy S, Lafosse F, Chaumeix N, Dupré G, Paillard CE. (2009) Hydrogen-nitrous oxide delay times: shock tube experimental study and kinetic modelling. Proc. Combust. Inst. 32: 359.

[27] Murray SB, Lee JHS. (1984) The influence of yielding confinement on large-scale ethylene-air detonations. Prog. Astronaut. Aeronaut. 94: 80.

[28] Grune J, Sempert K, Friedrich A, Kuzunetsov M, Jordan T. (2017) Detonation wave propagation in semi-confined layers of hydrogen-air and hydrogen-oxygen mixtures. Int. J. Hydro. Energy 42: 7589.

[29] Arienti M, Shepherd J. E. (2005) The role of diffusion at shear layer in irregular detonations. Joint Meeting of the US section of the Combustion Institute.

[30] Watanabe H, Matsuo A, Chinnayya A, Itouyama N, Kawasaki A, Matsuoka K, Kasahara J. (2023) Lagrangian dispersion and averaging behind a two-dimensional gaseous detonation front. J. Fluid Mech. 968: A28.