Numerical Study of Homogenous/Inhomogeneous Hydrogen–Air Explosion in a Long Closed Channel

Abstract

Hydrogen is regarded as a promising energy source for the future due to its clean combustion products, remarkable efficiency and renewability. However, its characteristics of low-ignition energy, a wide flammable range from 4% to 75%, and a rapid flame speed may bring significant explosion risks. Typically, accidental release of hydrogen into confined enclosures can result in a flammable hydrogen–air mixture with concentration gradients, possibly leading to flame acceleration (FA) and deflagration-to-detonation transition (DDT). The current study focused on the evolutions of the FA and DDT of homogenous/inhomogeneous hydrogen–air mixtures, based on the open-source computational fluid dynamics (CFD) platform OpenFOAM and the modified Weller et al.’s combustion model, taking into account the Darrieus–Landau (DL) and Rayleigh–Taylor (RT) instabilities, the turbulence and the non-unity Lewis number. Numerical simulations were carried out for both homogeneous and inhomogeneous mixtures in an enclosed channel 5.4 m in length and 0.06 m in height. The predictions demonstrate good quantitative agreement with the experimental measurements in flame-tip position, speed and pressure profiles by Boeck et al. The characteristics of flame structure, wave evolution and vortex were also discussed.

1. Introduction

The product of hydrogen combustion is only water and does not produce greenhouse gases or other harmful substances, making it an ideal clean energy source with a good prospect of development [1,2]. By electrolyzing water, excess electricity is converted into hydrogen, which can be used as a backup power source in the case of power shortages, as well as in transportation, industry, or residential areas. In addition, hydrogen can serve as a carbon-free long-term seasonal storage medium, solving the problem of seasonal supply–demand imbalance. Nowadays, hydrogen energy is driving the transformation and upgrading of the global energy structure. Various countries have incorporated hydrogen energy into their national energy strategies, promoting the development of the hydrogen energy industry through policy support and financial investment.

A booming demand for hydrogen energy requires a high degree of attention to the safety of generation, transportation, storage and usage of hydrogen, due to the lower minimum ignition energy and higher diffusion and combustion rates. The accidental release of hydrogen into confined spaces poses a significant safety concern, since it results in the formation of a flammable hydrogen–air mixture. If accidentally ignited, it can undergo rapid FA and potentially transition from deflagration to detonation, posing a high risk of explosion accidents. It is necessary to investigate the evolutions of possible flame acceleration and DDT in closed spaces.

Numerous investigations have been conducted to study hydrogen explosion, including flame propagation, deflagration-to-detonation transition and related parameters in the past few years [3,4,5,6,7]. Boeck et al. [8,9] focused on the effects of transverse concentration gradients on detonation propagation in hydrogen–air explosion tests in a closed rectangular channel. It was found that the mixture with transverse concentration gradients was accompanied by a lower detonation propagation speed compared with the identical average concentration mixture in the homogeneous case. The potential hazard of hydrogen–air mixtures has been reported in experiments by Vollmer et al. [10]. The tests take into account the effects of different hydrogen concentration, concentration gradient, and obstacles with differing blockage ratio and spacing. The rationality of 7λ criterion for DDT prediction in the homogeneous mixture was validated. Heidari et al. [11] developed two detonation solvers, one based on the solution of the reactive Euler equations and another based on the programmed CJ burn method, to investigate the mechanism of transition from deflagration to detonation. A further numerical study involving ignition strength was conducted [12]. Wang et al. [13] developed a density-based solver with single-step chemistry for hydrogen–air reaction to simulate flame acceleration and transition from deflagration to detonation. Azadboni et al. [14,15] conducted a numerical study within the OpenFOAM-7 platform and the effect of concentration gradients on DDT was investigated.

The physical understanding and quantitative description of FA and DDT are still critically needed, and moreover, numerical simulation would be an economical and efficient approach. Combustion in configuration with a large aspect ratio, such as with channels, tunnels and pipes, especially in obstructed conditions, is expected to provide more potential for spontaneous flame acceleration and transition from deflagration to detonation [16,17]. A numerical code for the 2D simulation of flame acceleration (FA) and deflagration-to-detonation transition (DDT) in homogenous/inhomogeneous hydrogen–air mixtures in a long-closed channel is developed in this study. In order to verify the reliability and accuracy of the numerical approach, the results of the simulation are compared against the available experimental data. Additionally, the evolutions of flame acceleration (FA) and deflagration-to-detonation transition (DDT) in homogenous/inhomogeneous hydrogen–air mixtures are also clarified, based on the contours of the pressure, temperature and hydrogen mass fraction. The purpose of this study aims at affording some theoretical support and fundamental data for hydrogen safety protection.

2. Numerical Models

2.1. Governing Equations

Both reactants and products are assumed to behave as ideal gases. The flow is governed by the compressible reactive Navier–Stokes equations and it can be written as follows:

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \stackrel{⃑}{U}\right)=0$$

(1)

$${\displaystyle \frac{\partial \rho \stackrel{⃑}{U}}{\partial t}}+\nabla ·\left(\rho \stackrel{\rightharpoonup }{UU}\right)+\nabla ·p=\nabla ·\tau$$

(2)

$$\begin{gathered} {\displaystyle \frac{\partial \rho E}{\partial t}}+\nabla ·\left(\rho \stackrel{⃑}{U}E\right)+\nabla ·\left(\stackrel{⃑}{U}p\right)=\nabla ·\left(\tau ·\stackrel{⃑}{U}\right)+\nabla ·\left({\lambda }_{eff}\nabla T\right)+\nabla · \\ \left({\sum }_{k=1}^{NS}\rho D_{eff}{h}_k\nabla Y_k\right) \end{gathered}$$

(3)

where $\rho$, $p$, $\stackrel{⃑}{\mathrm{U}}$, $\mathrm{E}$ and $Y_k$ are the density, pressure, velocity, total internal energy and mass fraction of the kth species, respectively. t is time, $\mathrm{\tau }$ is stress tensor, $h_k$ is the specific enthalpy of kth species, and NS is the species number. $D_{eff}$ and ${\lambda }_{eff}$ are the effective mass diffusion coefficient and the conductivity coefficient, and both involve the turbulent and laminar parts.

$$\rho E=\rho e+\rho {\displaystyle \frac{u^2+v^2+w^2}{2}}$$

(4)

$$\rho e={\sum }_{i=1}^{ns}{\rho }_i{h}_i-p$$

(5)

$$h_i={\int }_{T_{ref}}^T{C}_{pi}dT+h_{fi}^0$$

(6)

where $e$, $h_i$, $h_{fi}^0$ and $C_{pi}$ are the specific internal energy, enthalpy, the standard enthalpy and specific heat capacity at constant pressure of the ith species $u$, $v$ and $w$ denotes the velocities in $x$, $y$ and $z$ directions.

The combustion process is described by the reaction process variable c. c = 0 corresponds to the unburned gas, while c = 1 corresponds to the fully burned gas. The transport equation of the reaction process variable c can be expressed as:

$${\displaystyle \frac{\partial \left(\rho c\right)}{\partial t}}+\nabla ·\left(\rho \stackrel{⃑}{U}\mathrm{c}\right)-\nabla ·\left(\rho D\nabla c\right)={\mathsf{\omega }}_c$$

(7)

In RANS, the basic “laminar flamelet” models have been extensively extended. The flame front propagates locally as a laminar flame, but at the same time it is being wrinkled due to interactions with the turbulence. The flame propagation speed can be modelled using the laminar flame speed and flame wrinkling factor $\Sigma$. An alternative RANS model, Weller’s combustion model [18], represents the geometric properties of the flame front in terms of $\Sigma$, which is the flame area per unit area resolved in the mean direction of flame propagation. This choice makes this modelling easier than the equivalent equation for $\mathsf{\Sigma }$. The Weller’s combustion model [18] was employed to account for the deflagrative terms, together with a quenching factor of turbulent flames, $0\le G\le 1$ [19].

$${\omega }_c={\rho }_u\Sigma s_L\left|\nabla c\right|G$$

(8)

where $S_L$ and ${\rho }_u$ denote laminar burning speed and the density of the unburnt mixture.

As Bauwens et al. [20] suggested, the Darrieus–Landau (DL) and Rayleigh–Taylor (RT) instabilities and turbulence play an important role in hydrogen–air flame acceleration. Therefore, the flame wrinkling factor $\Sigma$ should be calculated as

$$\Sigma ={\Sigma }_t\times {\Sigma }_{DL}\times {\Sigma }_{RT}$$

(9)

where ${\Sigma }_t$, ${\Sigma }_{DL}$ and ${\Sigma }_{RT}$ denote the flame wrinkle factors caused by turbulence, Darrieus–Landau (DL) and Rayleigh–Taylor (RT) instabilities, respectively.

The transport equation model closure for the flame wrinkle factor ${\Sigma }_t$ is given as

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho {\Sigma }_t\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left(\rho {\Sigma }_t{u}_j\right)={\displaystyle \frac{\partial }{\partial x_j}}\left(\rho D_{eff}{\displaystyle \frac{\partial {\Sigma }_t}{\partial x_j}}\right)+\rho P_{{\Sigma }_t}{\Sigma }_t-\rho R_{{\Sigma }_t}{\Sigma }_t^2$$

(10)

where

$$P_{{\Sigma }_t}=0.28\sqrt{{\mathrm{C}}_{u^{\prime }}^3{\displaystyle \frac{\epsilon }{{\nu }_u}}}$$

(11)

$$R_{{\Sigma }_t}={\displaystyle \frac{P_{{\Sigma }_t}}{{\Sigma }_{eq}}}$$

(12)

$${\Sigma }_{eq}=1+C_{\xi }\left(1+2C_S\left(c-0.5\right)\right)\sqrt{{\displaystyle \frac{u^{\prime }}{s_L}}}{\mathrm{Re}}_{\eta }$$

(13)

$${\mathrm{Re}}_{\eta }={\displaystyle \frac{u^{\prime }}{{\left(C_{u^{\prime }}{\nu }_u\epsilon \right)}^{1/4}}}$$

(14)

where $C_{u^{\prime }}=1$, $C_{\xi }=0.62$ and $C_S=1$.${\nu }_u$, $\epsilon$ and $u^{\prime }$ are kinematic viscosity, the turbulence dissipation rate and the root-mean-square velocity, respectively.

The Darrieus–Landau and thermo–diffusive instabilities are modelled considering the expression as [20]

$${\Sigma }_{DL}=\mathrm{m}\mathrm{a}\mathrm{x}[1,{\alpha }_1{\left({\displaystyle \frac{\mathsf{\Delta }}{{\lambda }_c}}\right)}^{1/3}]$$

(15)

where, ${\mathsf{\lambda }}_c$, ${\mathsf{\alpha }}_1$ and Δ are the cutoff wavelength of unstable scales, a coefficient to uncertainty in ${\lambda }_c$ and the LES filter size, respectively.

The Rayleigh–Taylor instability is modelled as

$${\displaystyle \frac{\partial \overline{\rho }{\Sigma }_{RT}}{\partial t}}+{\widehat{U}}_s.\nabla {\Sigma }_{RT}=\overline{\rho }{G}_{RT}\left({\Sigma }_{RT}-1\right)-\overline{\rho }{R}_{RT}\left({\Sigma }_{RT}-1\right)$$

(16)

where, $G_{RT}\left({\Sigma }_{RT}-1\right)$ and $R_{RT}\left({\Sigma }_{RT}-1\right)$ are the rate of generation and the removal of sub-grid wrinkling owing to RT-instability.

The coefficients $G_{RT}$ and $R_{RT}$ are modelled as

$$G_{RT}=2{\left(k_{RT}{\displaystyle \frac{\sigma -1}{\sigma +1}}\stackrel{⃑}{a}.{\stackrel{⃑}{n}}_f\right)}^{1/2} and R_{RT}={\displaystyle \frac{8\sigma S_L{k}_{RT}}{\pi }}$$

(17)

where $\stackrel{⃑}{a}$ is the flame acceleration, $\sigma$ is the flame expansion ratio, $k_{RT}$ is the unstable wavenumber associated with the RT-instability.

2.2. Numerical Methods

In the Reynolds-averaged Navier–Stokes (RANS) framework, the convection term can be solved by the Harten-Lax-van Leer-Contact (HLLC) scheme [21,22,23,24,25]. The viscous terms are evaluated with second-order central-differencing discretization. The Crank–Nicholson second-order scheme was employed for time terms while the second-order central-difference scheme was used for the diffusion term. The Courant number $C_0=(U+a){\displaystyle \frac{∆t}{∆x}}$ is defined to calculate the time step. In the current study, the Courant number is set as 0.25, and therefore, the time step is adjustable since the flow speed and local acoustic speed change with time and space. The current case involves deflagration and detonation, so the time step ranges from 1 × 10⁻⁶ to 1 × 10⁻⁹.

3. Numerical Setup

In the current study, the numerical simulations were conducted in a long closed channel 5.4 m in length and 0.06 m in height for both obstructed and smooth channels full of homogeneous and inhomogeneous mixtures, which mimicked Boeck et al.’s [8]. The 30% blockage ratio (Br = 2h/H) was considered in the obstructed channel. Simultaneously, seven groups of obstacles with height h were installed on the top and bottom walls in the channel, as shown in Figure 1a. The obstacle space was 300 mm. The first obstacle was located at x = 250 mm, while the last obstacle was located at x = 2.05 m. The smooth channel is presented in Figure 1b, with the same height and length as the obstructed one. The average hydrogen mole fraction for both homogenous and inhomogeneous mixtures was 35%. The hydrogen concentration distribution for the inhomogeneous mixture is shown in Figure 2. The initial condition for the mixture was 1 atm and 293 K. In order to initiate the hydrogen–air mixture, a patch with a radius of 5 mm at the ignition point (x = 0, y = 0 m) was set with an adiabatic flame temperature of 2448 K and one atmospheric pressure.

4. Results and Discussion

Figure 3 presents the comparison between the predicted and measured flame-tip position vs. the time for the 35% homogenous mixture or the inhomogeneous mixture in the smooth or obstructed channel. It shows a good agreement between experimental and predicted data. The presence of obstacles has a significant effect on the flame propagation in the channel, since the time it takes for the flame to propagate from the ignition position to the end of the channel in a smooth pipe is about twice that in the obstructed channel.

Figure 4 and Figure 5 show the relationship between the flame velocity and the flame-tip position for the 35% homogenous/inhomogeneous hydrogen–air mixture in the smooth channel or the obstructed channel. It is easy to find a reasonable agreement between the predicted and measured values. At the initial stage, the flame propagates with a slow speed of around several meters per second. A great amount of heat is produced, accompanied with the combustion of the hydrogen–air mixture, which results in a volume expansion and a flame acceleration. Then, the flame propagates with a bigger flame speed and even reaches a high speed of more than 2000 m/s.

The flame speeds of the homogenous/inhomogeneous mixture in the experiment and the numerical simulation are presented in Figure 4 and Figure 5. Figure 4 is the one in the smooth channel while Figure 5 is the one in the obstructed channel. The predicted values are in reasonably good agreement with the experimental data. It can be found that the concentration gradient has a significant effect on flame speed in the smooth channel, and as a result, the flame speed increases faster in the inhomogeneous mixture, and the flame undergoes the transition to detonation at an earlier time than that of the homogenous mixture. The effect of the concentration gradient on the flame speed in the obstructed channel is not obvious.

In the smooth channel, with the homogeneous or inhomogeneous mixture, the flame-tip speed initially increases, which results in the enhancement of turbulence and DL and RT instabilities. These exert positive feedback to the flame. Resultantly, the flame speed increases continuously and is accelerated sharply at about the location of x = 4.0 m.

In the obstructed channel with the homogeneous or inhomogeneous mixture, the flame propagates with a bigger acceleration compared to that in the smooth channel. Flame speed increases fast and reaches the value of around 2000 m/s in the obstructed part (0 m < x < 2.05 m), indicating that the transition to detonation occurs here. With the presence of obstacles, the acceleration is extremely strong compared to that in the smooth channel. On the one hand, the obstacles induce the turbulence and the DL and RT instabilities, consequently increasing the flame’s surface area which promotes the burning rate. On the other hand, the obstacles provide a newly physical mechanism of flame acceleration, which is based on delayed burning between the obstacles, creating a strong jet flow and therefore causing the acceleration [26]. The flame propagates fast in the smooth part of the channel, leaving behind some unburned mixture between the obstacles which are burned later. The delayed burning in the empty space between the obstacles produces a strong jet flow in the unobstructed section of the channel. This jet flow causes a faster propagation of the flame tip, which creates new delayed burning in the front section and thereby affords positive feedback between the flame and the flow, driving the flame to be accelerated. This flame-acceleration mechanism can be powerful, causing supersonic flame propagation and then the transition from deflagration to detonation.

Figure 6 presents the predicted and measured pressure histories at specified probe locations (x = 0.4 m, x = 2.3 m, x = 3.2 m, x = 4.1 m) for the 35% homogenous hydrogen–air mixture in the obstructed channel (Br = 30%). The predicted pressure profiles are consistent with the experimental data. Moreover, the predicted first pressure peak times at x = 0.4 m, 2.3 m and 3.2 m are agreeable with the experimental data, while the predicted time for the first peak pressure at the position of 4.1 m is slightly earlier than the measured one. The predicted peak pressures are a little higher than in the experimental data, since the predicted detonation speed is about 10% faster than in the experimental data, and the time for deflagration to detonation is a little earlier than in the experimental value.

Figure 7 shows the predicted contours of temperature during the initial flame propagation for the 35% homogenous hydrogen–air mixture in the smooth and obstructed channels. At 4.1 ms, the flame front keeps the shape of a semi-circle and propagates at a slow laminar speed. Owing to the wave reflection at the obstacles, the flame propagation in the obstructed channel is slower than that in the smooth channel in front of the first group of obstacles. At 4.5 ms, the flame front is wrinkled and twisted. The flame surface area is increased, accompanied by a higher flame speed and a stronger heat release, which is caused by intrinsic flow instabilities. In the presence of obstacles, the flame is accelerated faster in the obstructed channel compared to the smooth channel.

The predicted contours of temperature, pressure and H₂ mass fraction during FA and DDT are presented in Figure 8, Figure 9, Figure 10 and Figure 11 for the 35% inhomogeneous mixture, respectively. For the convenience of observation, the parts of the smooth channel (from x = 0.85 m to x = 2.2 m) and the obstructed channel (from x = 2.95 m to x = 4.3 m) are presented. In Figure 8, before the transition from deflagration to detonation, the preheated areas are visible ahead of the flame front both in the smooth and obstructed channels. It indicates that the compressed wave induced by the burning of the mixture is ahead of the flame front and preheats the unburned mixture. At 16.5 ms in the smooth channel, a hot spot occurs near the location of 3.25 m on the bottom wall, which results in a subsequent explosion. A new flame produced by the explosion propagates with a faster speed than the flame front. Then it crosses the flame front and triggers the transition from deflagration to detonation at 16.8 ms.

In the obstructed channel, the transition from deflagration to detonation occurs much earlier, at 8.9 ms. More details can be observed in Figure 9. At 8.86 ms, the interaction between the wave and the obstacles induces a local explosion at the top of the sixth group of obstacles. At 8.87 ms, a subsequent one occurs at the bottom of the sixth group of obstacles. Resultantly, the transition from deflagration to detonation occurs.

In Figure 10, an obvious high-pressure region occurs in the bottom wall of the smooth channel at 16.5 ms, corresponding to the hot spot in temperature contour. It is caused by the reflection of the shockwave off the wall, which is accompanied by an increase in the temperature and the pressure near the wall. The reflection of shockwaves on the upper and bottom walls induces a powerful acceleration of the new flame, and a wider high-pressure region can be found in the smooth channel. The transition from deflagration to detonation is completed at 16.8 ms. In the obstructed channel, high-pressure regions appear near the obstacles, which is caused by the interaction between the shockwaves and the obstacles. It plays an important role in driving the transition from deflagration to detonation.

Figure 11 presents the predicted contours of H₂ mass fractions. There are pairs of vortexes ahead of the flame to the right of the obstacles, which are caused by the interaction between the obstacles and gas flow. Once the flame passes through these regions, the vortex pairs disappear, but the relatively high vertical hydrogen concentrations can still be found in the burned regions. However, no vortex can be seen ahead of the flame front in the smooth channel. It indicates that the presence of the vortex is caused by the obstacles. The hydrogen concentration gradient can still be found in both the smooth and obstructed channels in the burned regions. Since the process from ignition to detonation is fast and and about tens of ms, excess hydrogen remains in the upper part of the channel. This shows that a risk can still exist in the burned area, which may result in a re-ignition once enough fresh air is available in an accident scenario, inducing secondary explosion accidents.

5. Conclusions

Darrieus–Landau (DL) and Rayleigh–Taylor (RT) instabilities were taken into account together with the turbulence to modify the flame wrinkle factor in the Weller combustion model [18]. The flame acceleration (FA) and deflagration-to-detonation transition (DDT) were simulated in the smooth/obstructed tubes with the homogenous/inhomogeneous hydrogen–air mixtures. The predictions demonstrate good quantitative agreement with the experimental measurements in flame-tip position, speed and pressure profiles by Boeck et al. [8].

The analysis shows that a concentration gradient promotes flame acceleration in the smooth channel while the effect of the concentration gradient on flame acceleration and transition from deflagration to detonation is not obvious in the obstructed channel. DL and RT instabilities and turbulence play an important role in flame acceleration. The physical mechanism of flame acceleration by the obstacles in the obstructed channel is discussed in the current study, taking into account the wrinkle and jet flow which are caused by the obstacles. The details of the flame acceleration and the transition from deflagration to detonation have been discussed. Obstacles have an active influence on FA and DDT. DDT in the obstructed channel occurs much earlier than that in the smooth channel, due to the effect of the obstacles. The analysis on the distribution of inhomogeneous hydrogen concentration shows that the risk of secondary explosions still exists in the burned area.

In future work, more simulations will be carried out to further validate the current code, especially large-scale explosion scenarios. Moreover, more discussions will be focused on whether the above unburnt hydrogen left downstream of the flame will reburn or not if additional oxygen is involved.