Numerical Investigation of the Effects of the Diffusion Time on the Mechanisms of Transition from a Turbulent Jet Flame to Detonation in a H2-Air Mixture

1. Introduction

In the event of a hypothetical accident in a pressurized water reactor (PWR), there is the potential for a substantial release of hydrogen gas. This hydrogen can result from the oxidation of metals, either present in the reactor containment's basement during the phase of molten corium interacting with concrete or within the corium recovery pool. The H₂-air mixture can create a flammable mixture, leading to an explosion hazard. Local ignition of this flammable mixture can give rise to slowly propagating flames. The occurrence of the acceleration of a subsonic flame or deflagration transitioning to detonation depends on the turbulence level, mixture composition, and geometry. Hydrogen combustion, including detonation, could pose a significant threat to the reactor building and the integrity of the containment vessel in the event of an accident at a nuclear power station. Detonation events generate temperature peaks, shock waves, large pressure gradients, and high-pressure pulses, all of which have the potential to cause severe damage to specific equipment, internal walls, and containment components [1].

Some research in the field of detonation is focused on the fundamental problem of the deflagration-to-detonation transition (DDT) [2,3]. Nguyen et al. [4] conducted a numerical examination of the impacts of the equivalence ratio of ethylene fuel-air on the onset of detonation. They proposed a relationship between the equivalence ratio and the flame speed as a function of mixture fraction, pressure, and temperature.

In recent years, numerous researchers have conducted studies focused on detonation, with hydrogen serving as a fuel source [5,6,7,8]. Dorofeev et al. [9] analyzed the data from FA experiments in-depth. The data cover a wide range of tube scales and obstacle configurations, as well as mixture compositions. Temperatures and pressures were at normal and elevated levels. The suggested critical conditions for effective FA were presented in the form of correlations between the critical expansion ratio (${\mathsf{\sigma }}^{\mathrm{*}}$) and the dimensionless effective activation energy.

Considerable research effort has been dedicated to expediting the DDT process [10,11]. The placement of solid obstacles in the channel is a recognized method for reducing both the DDT distance and time. This is achieved by enhancing turbulence within the flow and increasing the interaction between the shock wave and the flame. Consequently, this approach enlarges the flame area and accelerates the release of the reaction heat rate. Ogawa et al. [12] conducted a comprehensive study on FA and DDT within a stoichiometric H₂-air mixture using a two-dimensional array of obstacles and a flame is ignited at the center of the obstacle array and propagates outward in all directions. Notably, their research highlighted the significant influence of flame propagation direction on both FA rate and DDT run-up distance. Breitung et al. [13] investigated the behavior of turbulent flames in H₂-air mixtures through a series of experimental studies. It was observed that within tubes featuring various obstacle arrangements, a clear distinction in flame behavior can be noted between slow, subsonic flames and fast-moving flames. Ciccarelli and Dorofeev [14] provided an overview of the FA process and discussed the mechanism of flame propagation in obstructed and unobstructed tubes with uniform H₂-air mixtures. Porowski et al. [15] employed the ddtFoam solver running on OpenFOAM and examined the impacts of obstacles on the onset of detonation in tubes with a stoichiometric H₂-air mixture. They noted that the ddtFoam solver had industrial safety applications. Their results showed that to obtain the shortest run-up distance, the flame needs a high blockage ratio (BR) with broad obstacle spacing or a lower BR with narrow gaps between subsequent rows of obstacles. Pinos and Ciccarelli [16] found that staggered obstacles noticeably enhance the average flame speed in the quasi-detonation regime in comparison with inline obstacles. Xiao and Oran [17] used numerical simulations to reveal the mechanism behind detonation transition in staggered obstacle arrays, demonstrating that it involves shock focusing at the flame. They also explored various obstacle shapes, such as circular, square, left triangular, and right triangular, and their distinct effects on FA and DDT [18]. Additionally, Han et al. conducted experimental examinations of DDT in obstructed channels [19]. They investigated mechanisms of FA in syngas-air in a closed tube. They also studied the effects of the distance of obstacles from the ignition source and hydrogen volume fractions on the pressure dynamics and FA. Their results showed that reducing the distance of obstacles from the ignition source affects the velocity and acceleration of the flame and maximizes overpressure. They described two processes contributing to FA: 1) the flame front is pushed by the unburned gas because its burns are delayed behind the flame front, and 2) the obstacles transition the laminar flame front to turbulent flame or intensify the turbulence rate. Baiwei et al. [20] conducted a numerical study of the transition from a turbulent jet flame to detonation in an H₂-air pre-mixed gas within a pipeline. The results revealed that the proportion of heat losses within the total heat release continuously increased before the flame crossed obstacles. However, as the flame advanced through the obstacles, it steadily decreased.

When hydrogen is released into an air-filled duct, specific behaviors occur. Hydrogen tends to rise to the top of the duct, creating an initially stratified mixture. This inhomogeneity results from the interplay between diffusion forces and buoyancy. Consequently, the mixture initially establishes a vertical concentration gradient, with hydrogen concentration being higher toward the top of the duct. Over time, this concentration gradient evens out, and the mixture becomes more uniform on a global scale. In the case of a large volume, like a reactor containment, this process can extend over several hours before achieving a fully uniform state. Despite this realistic scenario, research on the combustion of mixtures with non-uniform concentrations has been relatively limited. Given that an initially non-uniform fuel distribution more accurately represents the initial state of the H₂-air mixture in an accident scenario [21], it is crucial to consider this factor in safety-related combustion studies. Ishii and Kojima [22] focused on detonation propagation in mixtures with vertical concentration gradients in a horizontal tube. Their results demonstrated that the local deflection angle increases with an increase in the local concentration gradient. Bleyer et al. [23] investigated flame propagation with a vertical hydrogen gradient in a vertical tube. Their results demonstrated that the flame propagates in the upward direction, away from the igniter, into a gradually richer or leaner mixture depending on the concentration gradient. Also, they showed that the concentration at the ignition source has a significant influence on the flame development. Song et al. [24] conducted numerical investigations into self-sustaining detonation modes within heterogeneous hydrogen-oxygen mixtures. They observed that concentration gradients played a pivotal role in alternating between multi-head and single-head detonation modes. Saied et al. [25] investigated DDT in uniform and stratified mixtures with 15% and 30% hydrogen concentrations. They found that in the lean mixture (15% hydrogen), DDT occurred only in stratified conditions. They also showed that reducing obstacle spacing caused DDT to happen earlier in time and space. Luan et al. [26] performed numerical simulations of 2D rotating detonations. They discovered that the wave number progressively increased with higher equivalent ratios in stratified environments. Jiang et al. [27] examined the impact of concentration gradients on detonation and re-initiation behavior within oxygen-containing mixed gases in bifurcation ducts. Their study revealed that in mixtures with concentration gradients, the second reflection plays a crucial role in successfully re-initiating detonation. This is in contrast to the uniform mixture. Saied et al. [28] conducted a numerical investigation into the effect of hydrogen mixture inhomogeneity on the mechanisms underlying the DDT. Their study identified three distinct mechanisms of DDT within a rich stratified mixture at a 35% concentration. In the first regime, the creation of a robust Mach stem among the obstacles results in the direct initiation of detonation in the vicinity of the triple point. In the second regime, reflected shock waves from the obstacles serve as the necessary conditions for DDT. In the third regime, neither the Mach stem nor the reflected shock waves play a substantial role in the onset of detonation.

In recent years, there has been extensive experimental research on FA and DDT. However, detonation is a very complex phenomenon with a microsecond time-scale, and investigating the details of its occurrence requires advanced and costly techniques. Therefore, numerical simulation of this phenomenon, as a suitable and accurate tool, can minimize the need for experimental tests. On the other hand, many of the related studies have focused on the structure of gas detonation, often overlooking the concentration on FA and DDT phenomena. In the last few years, much research has been dedicated to understanding the factors affecting DDT and how to control it. Accurate calculations necessitate more realistic models and appropriate numerical methods. As such, this research has focused on stratified mixtures for more realistic modeling and more detailed investigations. The low density of hydrogen makes it prone to forming non-uniform concentrations when mixed with air in a combustion chamber. This non-uniformity significantly impacts both FA and DDT processes. Therefore, it is essential to conduct research aimed at understanding how jet flames influence FA and DDT processes within stratified H₂-air premixed gases.

This research has focused on the numerical 2D investigation of the propagation, acceleration, and the transition from a turbulent jet flame to detonation. The flame acceleration process comprises both an initial and a final stage. In the first stage, the flame begins accelerating after the reactions commence, reaching the sound velocity (in reactive materials). In the final stage, the process initiates with the formation of strong pressure waves in front of the flame. As the flame continues to accelerate, its speed approaches the velocity of sound in the combustion products, resulting in a choked flame. Several lines of evidence showed that this choked flame speed can be reached in a turbulent flame where hot spots are formed within the channel, ultimately leading to DDT. This research has employed numerical simulations to explore the mechanisms influencing flame acceleration in both the initial and final stages and has also discussed the interaction between combustion and turbulence. Also, this study focuses on evaluating the influence of diffusion time, delving into factors that have received less attention in previous research on the transition from a turbulent jet flame to detonation. Through numerical simulations of gas detonation, we have examined how diffusion time affects the fundamental processes responsible for the transition from a turbulent jet flame to detonation.

2. Governing Equation

The governing equations for the numerical simulations are expressed as follows:

Continuity equation

$$\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial }{{\partial x}_j}\left(\overline{\rho }{\tilde{u}}_j\right)=0$$

(1)

Momentum equation

$$\frac{\partial }{\partial \mathrm{t}}\left(\overline{\rho }{\tilde{u}}_i\right)+\frac{\partial }{{\partial \mathrm{x}}_{\mathrm{j}}}\left(\overline{\rho }{\tilde{u}}_i{\tilde{u}}_j+\frac{2}{3}{\delta }_{ij}\overline{\rho }\mathrm{k}\right)=-\frac{\partial \overline{p}}{{\partial x}_i}+\frac{\partial }{{\partial x}_j}\left({\mu }_{eff}\left(\frac{\partial {\tilde{u}}_i}{{\partial x}_j}+\frac{\partial {\tilde{u}}_j}{{\partial x}_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial {\tilde{u}}_m}{{\partial x}_m}\right)\right)+\overline{ \rho }{\mathrm{g}}_{\mathrm{i}}$$

(2)

Energy equation

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{e}}_t\right)+\frac{\partial }{\partial xj}\left(\left(\overline{\rho }{\tilde{e}}_t+\overline{p}\right){\tilde{u}}_j\right)=\frac{\partial }{{\partial x}_j}\left(\overline{\rho }{a}_{eff}\frac{\partial }{{\partial x}_j}\left({\tilde{e}}_t+\frac{\overline{p }}{\overline{\rho }}\right)+\mu \left(\frac{\partial {\tilde{u}}_i}{{\partial x}_j}+\frac{\partial {\tilde{u}}_j}{{\partial x}_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial {\tilde{u}}_m}{{\partial x}_m}\right){\tilde{u}}_i\right)$$

(3)

Reaction progression variable equation

$$\frac{\partial }{\partial t}\left(\overline{\rho }\tilde{\mathrm{c}}\right)+\frac{\partial }{{\partial x}_j}\left(\overline{\rho }{\tilde{u}}_j\tilde{\mathrm{c}}\right)=\frac{\partial }{{\partial x}_j}\left(\overline{\rho }{D}_{eff}\frac{\partial \tilde{\mathrm{c}}}{\partial x_j}\right)+{\overline{\mathsf{\omega } }}_{\mathrm{c}.\mathrm{d}\mathrm{e}\mathrm{f}}+{\overline{\mathsf{\omega } }}_{\mathrm{c}.\mathrm{i}\mathrm{g}\mathrm{n}}$$

(4)

Mixture fraction equation

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{f}}_H\right)+\frac{\partial }{\partial x_J}\left(\overline{\rho }{\tilde{u}}_j{\tilde{f}}_H\right)=\frac{\partial }{\partial x_J}\left(\overline{\rho }{D}_{eff}\frac{\partial {\tilde{f}}_H}{\partial x_j}\right)$$

(5)

Additionally, the chemical kinetic model used is based on the detailed reaction mechanism of O'Conaire et al. [29], which includes 9 species and 21 reversible chemical reactions and the transport model for species is based on the Sutherland model.

2.1. The Turbulence Models

The mechanism of the explosion is described as a loop with positive feedback. In this process, the flame acceleration, and as a result, the increase of turbulence rate, cause the flame to propagate faster. When the ignition occurs in an explosive mixture of fuel and air, the hot gases from combustion exert a force on the forefront of the unburned gases, pushing them into the path of the flame propagation. Due to the pressure gradient, the flow accelerates and turns to turbulent regime once it collides with the obstacles. This collision of flow with obstacles increases turbulence mixing. As a result, the reaction rate grows as the combustion process accelerates, where the flow velocity and hence the turbulence rises. The whole process repeats itself in a periodic manner. Thus, the role of turbulence in modeling the DDT is vital. The turbulence model used in this study was the two-equation model of 𝑘−𝜔 𝑆𝑆𝑇. It is an improved form of the 𝑘−𝜔 model to achieve more accurate results near all walls in the combustion chamber. In the boundary layer, the 𝑘−𝜔 𝑆𝑆𝑇 turbulence model outperforms 𝑘−𝜔 in the viscous sub-layer (near the wall). This model can be utilized as a low Reynolds turbulence model without any additional damping function. In the 𝑘−𝜔 𝑆𝑆𝑇 model, the turbulence relations in the free-flow are turned into relations of the 𝑘-ε model; however, the 𝑘−𝜔 𝑆𝑆𝑇 model does not have 𝑘-ε problems, such as sensitivity to the properties of input free-flow, lack of sensitivity to adverse pressure-gradients, and numerical stiffness of the equations when integrated through the viscous sub-layer. The 𝑘−𝜔 𝑆𝑆𝑇 turbulence model based on RANS is written as follows:

$$\frac{\partial }{\partial t}\left(\overline{\rho }k\right)+\frac{\partial }{{\partial x}_i}\left(\overline{\rho }k{\tilde{u}}_i\right)=\frac{\partial }{{\partial x}_j}\left((\mu +\frac{{\mu }_T}{\partial k})\right)+2{{\mu }_T{\mathrm{S}}_{\mathrm{i}\mathrm{j} }{\mathrm{S}}_{\mathrm{i}\mathrm{j} }}_{ }-\frac{2}{3}\overline{\rho }k\frac{\partial {\tilde{u}}_i}{{\partial x}_i} {\delta }_{ij}$$

(6)

$$\frac{\partial }{\partial t}\left(\overline{\rho }\mathsf{\omega }\right)+\frac{\partial }{{\partial x}_i}\left(\overline{\rho }\mathsf{\omega }{\tilde{u}}_i\right)=\frac{\partial }{{\partial x}_i}\left((\mu +\frac{{\mu }_T}{\partial \mathsf{\omega },1})\frac{\partial \mathsf{\omega }}{{\partial x}_i}\right)+{\gamma }_2\left(2\overline{\rho }{\mathrm{S}}_{\mathrm{i}\mathrm{j} }{\mathrm{S}}_{\mathrm{i}\mathrm{j} }-\frac{2}{3}\overline{\rho }\mathsf{\omega }\frac{\partial {\tilde{u}}_i}{{\partial x}_i} {\delta }_{ij}\right)-{\beta }_2\overline{\rho }{\omega }^2+2\frac{\overline{\rho }}{{\sigma }_{\mathsf{\omega },2}\mathsf{\omega }}\frac{\partial \mathrm{k}}{{\partial x}_j}\frac{\partial \mathsf{\omega }}{\partial \mathrm{k}}$$

(7)

In the above equations, the turbulent kinetic energy $\mathrm{k}=\frac{1}{2}\widetilde{\stackrel{´}{u_i}\stackrel{´}{u_j}}$, specific turbulent dissipation rate $\mathsf{\omega }=\frac{\epsilon }{k}{\beta }^{\mathrm{*}}$, and model coefficients ${\gamma }_2=0.44,{ \beta }_2=0.083,$ and ${\beta }^{\mathrm{*}}=0.09$ [30].

3. Numerical Method

In the present numerical simulation, the ddtFoam solver [31] was used to simulate FA and DDT in the H₂-air mixture. The primary objective was to compute the macroscopic properties of shock propagation and the characteristics of the flame. DDT is a complex gas dynamic phenomenon that involves flame-turbulence interaction and flame instabilities. It is virtually impossible to resolve all microscopic details of the flow in industry-scale geometries with current resources and capabilities. Moreover, not all microscopic-scale phenomena are necessary for simulating macroscopic DDT events. In the current study, a coarse grid size was employed. Using under-resolved grids serves two purposes: it enables computational scalability for larger domains and facilitates the study of FA and DDT within a limited timeframe. To simulate flame propagation in large three-dimensional geometries without resolving the microstructure of the flow within the Computational Fluid Dynamics (CFD) cell, a combination of the PISO solver and the ddtFoam solver [32] was employed. The PISO solver is well-suited for Mach numbers less than 0.3, while the ddtFoam solver is activated when the Mach number exceeds 0.3. This combined approach serves as a foundation for studying FA and the initiation of detonation processes.

In the current study, the Unsteady Reynolds-Averaged Navier-Stokes (URANS) modeling approach was chosen. This selection was used by its computational efficiency, which permits the use of a comparatively larger mesh size. Additionally, URANS modeling offers reasonably accurate estimations of macroscopic parameters, making it a suitable option for the research objectives. The ddtFoam is a solver to simulate the partial-premixed combustion flows with the two-equation model of 𝑘−𝜔 𝑆𝑆𝑇. The combustion model in this solver is the two-equation model of c-Ξ Weller, which is one type of Flamelet model. Here, c indicates the mass fraction of products (the indication of the reaction progress), and Ξ represents the flame surface wrinkle. Instead of solving separate transport equations for individual chemical species reactions, a unified transport equation was created and resolved for the variable c. The auto-ignition delay time (τ_ign) was devoted using Arrhenius equations, relying on the comprehensive reaction scheme described by O'Conaire et al. [29]. The mixture composition is explained by the mixture fraction (𝑓_H), which indicates the hydrogen amount that would be found if the cell were completely unburned. To avoid repetitive recalculations of local ignition delay times, a table of τ_ign based on temperature (T), pressure (P), and mixture fraction (𝑓_H) was generated using Cantera [33]. Throughout the numerical solutions, for each computational cell, the local ignition delay time was determined by searching this tabulated data. In this solver, an additional energy equation was solved for the unburned mixture. By solving this equation and obtaining the enthalpy of other unburned gases, the thermodynamic properties of these gases, such as temperature, density, and dynamic viscosity were calculable.

In the ddtFoam solver, the detonation process involving DDT and the effects of auto-ignition are described by the detonative source term (${\overline{\mathsf{\omega }}}_{\mathrm{c}.\mathrm{i}\mathrm{g}\mathrm{n}}$) in the transport equation of the reaction progress variable. If a gas mixture ignites spontaneously without any ignition source, auto-ignition has occurred. The mixture composition, pressure (P), and local temperature (T) control the auto-ignition delay. In this solver, a sub-model is used, which raises the accuracy of the auto-ignition model in the coarse grids. This sub-model is presented, because the auto-ignition can occur due to shock-induced heating. Providing that the model predicts the auto-ignition based on average pressure and temperature, incorrect results can be obtained. The transition to detonation was modeled using the auto-ignition source term as follows:

$${\overline{\omega }}_{c.ign }=\overline{\rho }\frac{1-\tilde{c}}{\Delta t} H (\tilde{\tau }-1)$$

(8)

In the above equation, ∆𝑡 indicates the present time step, H indicates the Heaviside function, and 𝜏 is a dimensionless variable expressing the auto-ignition process [32]:

$$\mathsf{\tau }=\frac{\mathrm{t}}{t_{\mathrm{i}\mathrm{g}\mathrm{n}}}=\frac{{\mathrm{y}}_R}{y_{\mathrm{R},\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}}$$

(9)

The mixture is ignited when the ignition variable 𝜏 reaches unity. Nevertheless, when the ignition delay time has not yet passed (𝑡 < 𝑡_ign), there is no effect on the properties of flow. The local ignition delay time is a function of local pressure 𝑝, local temperature 𝑇, and mixture composition. This ignition delay time of a mixture is calculated using lookup tables, obtained from isochoric explosion calculations, which are calculated using a detailed reaction mechanism and Cantera.

The average temperature may be adequate to initiate auto-ignition, but the shock responsible for temperature increase might not have fully traversed the entire computational cell. As a result, the computational cell can split into two segments, one with a volume fraction α and the other with a volume fraction of 1-α (as illustrated in Figure 1). In one segment, the temperature and pressure are denoted as $T_{high}$ and $P_{high}$, while in the other segment, they are represented as $T_{low}$ and $P_{low}$, respectively. These values are calculated based on the properties of the surrounding computational cells. The mean pressure within the computational cell is denoted as $p̀$, and the volume fraction is determined using the following equation.

$$\frac{{p̀}_{-}{p}_{low}}{p_{high-}{p}_{low}}=\mathsf{\alpha }$$

(10)

The values $T_{low}$ and ${ T}_{high}$ can be computed using dynamic shock relations for an ideal gas with a specific heat ratio represented by 𝛾 [34]. This model enables the distinct calculation of the ignition delay time across the shock within the computational cell and is applicable to the entire computational domain, even in regions where no shock is present.

In every grid cell, the process of auto-ignition is assessed individually on both sides of the discontinuity. By solving transport equations for 𝜏_high and 𝜏_low, mixing and transport consequences are accounted.

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{\mathsf{\tau }}}_{high}\right)+\frac{\partial }{{\partial x}_j}\left(\overline{\rho }{\tilde{u}}_j{\tilde{\mathsf{\tau }}}_{high}\right)=\frac{\partial }{{\partial x}_j}\left(\overline{\rho }\left(D{+D}_T\right)\frac{\partial {\tilde{\mathsf{\tau }}}_{high}}{\partial x_j}\right)+\frac{\overline{\rho }}{t_{ign,high}}$$

(11)

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{\mathsf{\tau }}}_{low}\right)+\frac{\partial }{{\partial x}_j}\left(\overline{\rho }{\tilde{u}}_j{\tilde{\mathsf{\tau }}}_{low}\right)=\frac{\partial }{{\partial x}_j}\left(\overline{\rho }\left(D{+D}_T\right)\frac{\partial {\tilde{\mathsf{\tau }}}_{low}}{\partial x_j}\right)+\frac{\overline{\rho }}{t_{ign,low}}$$

(12)

If 𝜏 = 1 (the critical value) is reached on one side of the discontinuity (i.e., 𝜏_high = 1 or 𝜏_low = 1), only the corresponding volume fraction of a computational cell is ignited. Therefore, the source term of auto-ignition will be calculated as shown [35]:

$${\overline{\mathsf{\omega }}}_{\mathrm{c}.\mathrm{i}\mathrm{g}\mathrm{n}}=\mathsf{\alpha }\frac{1-\tilde{\mathrm{c}}}{\Delta t} H \left( {\tilde{\mathsf{\tau }}}_{high}-1\right)+(1-\mathsf{\alpha })\frac{1-\tilde{\mathrm{c}}}{\Delta t} H ( {\tilde{\mathsf{\tau }}}_{low}-1)$$

(13)

In the above equation, ∆𝑡 indicates the current time step and H indicates the Heaviside function.

$$H \left(x\right)=\left\{\begin{array}{c}0, x<0\\ 1, x \ge 0\end{array}\right.$$

(14)

The Heaviside function triggers solely the part of the computational cell in which the local ignition delay time has passed, gained by weighting the volumetric source term by the volume fraction of either 𝛼 or 1−𝛼.

In computational simulations, it's crucial to use a grid resolution of 30-50 grid points within the flame thickness region to accurately capture the propagation of the flame. However, this finer resolution can significantly increase computational costs. To address this challenge, the ddtFoam solver has been specifically designed for simulating accidental explosions in nuclear plants. It employs a volume fraction method within each cell to determine the auto-ignition delay time, delivering satisfactorily accurate solutions even on coarser grids. The ddtFoam solver is versatile and can be effectively applied to problems involving detonation diffraction and various time-dependent extreme combustion scenarios. It offers the advantage of providing relatively quick results while maintaining an acceptable level of accuracy, making it a valuable tool for analysis and design purposes. The material properties were chosen from the Chemkin database through look-up tables. The Sutherland correlation was utilized to evaluate the molecular transport coefficients. The problem was a transient mode and was solved with the explicit Euler scheme for time discretization. In spatial discretization for the diffusion terms, second-order schemes were used to obtain higher accuracy. Tolerance of 10^-6 was set for all solution algorithms, the value set for the time step size was 2*10^-5, which was also the initial size of the time step.

3.1. Characteristics of Combustion Geometry

The numerical investigation of FA in stratified mixtures was carried out within a closed duct measuring 5.1 meters in length (Figure 2a). The duct had a rectangular cross-section with a height of 60 mm and a width of 300 mm. Inside the duct, a turbulence-inducing obstacle was installed. In this particular configuration, the blockage ratio was set at 0.9, and the obstacle was positioned at a distance of 1 meter from the ignition source, which was mounted on one of the end plates of the tube. Notably, the dimensions of this duct are identical to those of the Gravent facility [36]; however, the Gravent facility features seven obstacles. All cases were conducted with an H₂-air mixture at an initial pressure of 1 atm and an initial temperature of 293 K. The vertical concentration gradient achieved depends on the diffusion time, which is defined as the time between the end of hydrogen injection and ignition. Gradients with a predetermined slope can be created by regulating the diffusion time (t_d) between H₂ injection and ignition. Figure 2b provides an initial visualization of concentration gradients.

3.2. Mesh Sensitivity Analysis

To assess the sensitivity of the results to grid size, we generated three different computational grids in the GraVent facility [36] (Figure 3) with cell sizes of 4 mm (mesh-1), 1 mm (mesh-2), and 0.5 mm (mesh-3). Figure 4.a demonstrates the effects of different grids on the position of the flame front at different times. In Figure 4b, the channel pressure history recorded at x = 2.4 m over time for these three grids are shown. These simulations were conducted using a stratified mixture of 20% hydrogen with diffusion time of 3 seconds. Notably, our analysis revealed no significant differences between cell sizes of 1 mm (mesh-2) and 0.5 mm (mesh-3).

To validate our current simulations, we compared our numerical results with the experimental observations conducted by Boeck et al. [36] and those obtained by Karanam et al. [37] in the GraVent facility. The flame velocity, featuring a diffusion time of 3 seconds and an average hydrogen concentration of 20%, is depicted in Figure 5. This graph demonstrates a strong agreement between the experimental observations and our simulation results, confirming the capability of our solver to accurately reproduce properties such as turbulent and laminar flame velocities, as well as the speed of unsteady detonation propagation in stratified mixtures.

3.3. Numerical Validation

In our simulations, we used a uniform structured grid within a 2D computational domain. The smallest cell size was set to 1 mm, following a similar methodology as the simulations carried out in the previous section and as described in the work by Saeid et al. [28]. No-slip and adiabatic boundary conditions were imposed on the walls of the channel, and the spark zone was represented as a hot spherical region with a radius (R) of 25.4 mm, where the reaction progress variable reached 1, and the temperature reached 2440 K. Due to the absence of any existing simulation or experimental results for the channel shown in Figure 2, we examined the channel depicted in Figure 2 in its unobstructed state to validate the accuracy of our current simulations. Figure 6 presents a comparison between our simulation results and the experimental data provided by Boeck et al. [38]. Clearly, our numerical simulations closely match the overall trend observed in the experimental data. Furthermore, when we compare our results to those presented by Karanam et al. [39], it becomes evident that our calculated results are in strong agreement with their outcomes. This collective evidence underscores the satisfactory performance of our model, yielding qualitatively acceptable results for simulating both deflagration and DDT phenomena.

4. Results

4.1. Effects of Inhomogeneity on the FA and DDT in Lean Mixtures

In Figure 7, which shows the effects of inhomogeneity on the flame speed in the channel depicted in Figure 2 for the H₂-air mixture with an average concentration of 22.5%, when the mixture is uniform, DDT occurs at the channel's end. However, when its inhomogeneity increases (diffusion time reduces to 3 seconds), the run-up distance of DDT occurrence decreases.

Figure 8 and Figure 9 illustrate temperature contours within a channel featuring a 0.9 blockage ratio, filled with a stratified and uniform H₂-air mixture with an average hydrogen concentration of 22.5%. These figures provide insights into the two distinctive stages of FA:

Initial flame propagation stage:

In Figure 8, we observe the initial stages of ignition and flame propagation. In the stratified mixture (Figure 8a), during the first 4.0 ms, the ignition zone gradually expands. The flame front undergoes wrinkle formation due to natural flame instability in lean mixtures, leading to an increase in surface area around t = 5.0 ms. As time progresses, the flame propagates upward, driven by the higher hydrogen concentration present there compared to the duct floor. In a uniform mixture (as shown in Figure 8.b), it is evident that the flame shape continuously changes during its propagation. Initially, at t = 5.0 ms, the flame exhibits a laminar behavior, expanding in a symmetric shape. However, at t = 7.0 ms, the flame undergoes significant wrinkling and propagates towards an obstacle. In the stratified mixture, as the flame surface becomes increasingly elongated, it is anticipated that flame acceleration will be enhanced in the direction of propagation. This elongation directly influences flame speed by increasing the flame's surface area and, consequently, boosting the rate of combustion. As a result, FA in stratified mixtures is expected to be higher than in uniform mixtures.

Jet flame passing through obstacle stage:

Figure 9 captures the jet flame stage for stratified and uniform mixtures. In Figure 9a, at t = 11.7 ms, the flame's behavior can be observed, showing a preference for propagating towards the upper regions of the domain, where the hydrogen concentration is higher. At t = 11.8 ms, the flame successfully passes through the obstacle, giving rise to a jet-like flame on the right side. At t = 11.9 ms, the effects of vortices create an unburned zone in front of the obstacle. The significant difference in flame propagation between the stratified and uniform configurations can be attributed to the flame's shape and surface area. In Figure 9b, the uniform configuration presents an almost symmetrical flame front situated at a constant axial position at t = 17.8 ms, prompting the flame to ascend due to the buoyancy force in the upper part of the duct. As the flame continues, it expands and burns, leading to the formation of a relatively not-so-strong compression wave at t = 18.1 ms, just ahead of the advancing flame front.

Now, the jet flame to detonation transition mechanism is addressed in uniform and stratified H₂-air mixtures with diffusion times of 3 and 5 seconds inside a duct with a blockage ratio of 0.9. Contours of the temperature and pressure in Figure 10a belong to a 3 second-diffusion time stratified mixture. In this figure, a hot region has been formed upstream of the obstacle at 12.355 ms, causing an explosion center to form. Subsequently, due to the concentration of pressure waves, a detonation wave is formed, gradually occupying the channel width. In Figure 10b, which shows the temperature/pressure contours for a 5 second-diffusion time, the preheated area, in the temperature contour at 13.375 ms, lies ahead of the flame front. Here, the upper-wall shock reflection creates a hot spot at 13.385 ms, which, over time, causes a detonation to begin. The new, detonation-induced combustion wave moves fast and, finally, the detonation wave due to local explosion and the leading shock wave merge. In Figure 10c that shows the temperature/pressure contours at the channel end for a fully uniform mixture, its flame front is symmetrical before reaching the obstacle, but becomes asymmetric after passing it. This flame front asymmetry in uniform mixture has also been observed in the experimental work of Gaathaug et al. [40]. An explosive center is formed, at about 22.980 ms, behind the leading shock at the lower wall of the channel and moves towards unburned materials as a detonation wave, which gradually covers the entire channel width.

4.2. Effects of Inhomogeneity on FA and DDT in Stoichiometric Mixtures

Figure 11 illustrates the impact of mixture inhomogeneity on flame speed within the channel, represented in Figure 2. This channel contains an H₂-air mixture with an average concentration of 30%. As the flame advances towards an obstacle positioned 1 meter from the channel's entrance, it accelerates and attains a higher velocity as it passes through the obstacle, resulting in the formation of a supersonic jet. Later, the interaction of the flame with this jet flow will create a high-speed turbulent flame. In fact, when the flame enters this turbulent jet, its surface area experiences a sudden sharp increase, causing the energy release rate and, hence, the flame propagation speed to increase suddenly. In this scenario, the flame's surface area increases significantly compared to previous conditions, leading to the formation of a robust shock wave in front of it. This shock wave causes the flame to undergo a transition to detonation in less than 0.5 meters of travel. As depicted in this figure, in contrast to lean mixtures, for stoichiometric mixtures, decreasing the diffusion time (or increasing the mixture's non-uniformity) diminishes the FA and postpones the transition to detonation. In the uniform mixture, the DDT event has occurred in the shortest time and at the closest distance from the channel head.

Changing the concentration and diffusion time will change the DDT-process governing-mechanisms. Figure 12a shows the temperature/pressure contours for a stratified mixture with an average concentration of 30%, a diffusion time of 3 seconds and a blockage ratio of 0.9. Here, at 11.590 ms, an explosion center occurs in an area containing a lean mixture, which is the place where the reflected shock wave comes from Mach reflection and the flame front collide. This hot spot is the main source of the detonation onset, which gradually occupies the channel width. Figure 12b demonstrates the final steps of the DDT process in a stratified H₂-air mixture with a 5 second-diffusion time. At 11.135 ms, a high-temperature zone or a hot spot appears in the temperature contour near the upper wall of the channel, which leads, over time, to a local explosion and transition to detonation. This overdriven detonation wave propagates rapidly and merges with the leading shock wave. In short, the reason for the DDT occurrence is the interaction between the upper wall, leading shocks, and the area with higher hydrogen concentration. In Figure 12.c, which shows the temperature/pressure contours of the DDT process for a uniform mixture, a relatively strong curved shock wave is propagating in front of the flame. At 11.120 ms, a high-temperature area (hot spot) appears in the temperature/pressure contour near the lower wall of the channel due to the concentration of the pressure waves near the flame front. With the passage of time, this hot spot leads to a local explosion causing an overdriven detonation wave to form.

5. Conclusions

This study examined the effects of the fuel diffusion time on the FA and detonation initiation in uniform and stratified H₂-air mixtures with different concentrations in a turbulent jet flame. According to the results, the diffusion time and concentration of the H₂-air mixture are influential factors in the FA and the DDT location.

• In the average 22.5% concentration H₂-air mixture, detonation began in all the cases, and comparing the flame speed graphs showed that a decrease in the diffusion time and an increase in the mixture inhomogeneity sped up the flame propagation and the DDT occurrence in the channel.
• In the average 30% concentration H₂-air mixture too, deflagration to detonation transition occurred in all the cases, and the mixture inhomogeneity weakened the FA and delayed the DDT.
• This research showed that the inhomogeneity behavior differed between lean and stoichiometric mixtures; in the lean mixture, increasing the inhomogeneity reduced the run-up distance to detonation while in the stoichiometric mixture, increasing the diffusion time (homogeneity) did so.
• A feature of the turbulent jet flame to detonation transition study was the formation of a supersonic turbulent jet in the downstream part of the obstacle. At a certain stage of the FA, the speed at the obstacle bottleneck reached that of sound, causing the flame to chocked.

In general, since this study considers different physical phenomena, it enables an accurate evaluation of the present code. An accurate modeling of the energy released during the initial spark has great effects not only on the early stages of the flame expansion, but also on the shape of the core of the expanding flame. Therefore, an accurate modeling of the spark to achieve the times reported in the experimental results is suggested to continue this research in the future.