Abstract
In this study the hazardous potential of flammable hydrogen-air mixtures with vertical concentration gradients is investigated numerically. The computational model is based on the formulation of a reaction progress variable and accounts for both deflagrative flame propagation and autoignition. The model is able to simulate the deflagration-to-detonation transition (DDT) without resolving all microscopic details of the flow. It works on relatively coarse grids and shows good agreement with experiments. It is found that a mixture with a vertical concentration gradient can have a much higher tendency to undergo DDT than a homogeneous mixture of the same hydrogen content. In addition, the pressure loads occurring can be much higher. However, the opposite effect can also be observed, with the decisive factor being the geometric boundary conditions. The model gives insight into different modes of DDT. Detonations occurring soon after ignition do not necessarily cause the highest pressure loads. In mixtures with concentration gradient, the highest loads can occur in regions of very low hydrogen content. These new findings should be considered in future safety studies.
1. Introduction
Due to its potentially catastrophic consequences, the accidental release of hydrogen is a major concern in process engineering [1, 2], power generation [3–5], and future automotive concepts [6, 7]. A small heat source can be sufficient to initiate a deflagration that causes a sudden temperature and pressure rise. The severity of the consequences depends on the propagation speed of the deflagration. In general it can be expected that the flame accelerates due to a positive feedback loop: the expansion of the combustion products induces motion in the unburned gas which generates turbulence and thereby increases the burning speed of the flame. This process is further complicated through instabilities [8–10], formation of shocks, and interaction with the surrounding structure. Under certain circumstances the deflagration can undergo a transition to a detonation. The hazardous potential of a detonation (and especially that of the unsteady transition process) is considerably higher than that of a deflagration.
Therefore, the probability of the occurrence of a deflagration-to-detonation transition (DDT) needs to be considered in safety studies involving scenarios with accidental hydrogen release. In the past, many researchers have contributed impressive DDT simulation studies by applying a very high grid resolution to observe all microscopic flow details [11–15]. This approach has helped a lot to understand the formation of DDT but is only applicable to very small domains. When using simple reaction models like one-step Arrhenius kinetics, many publications can be found using a grid resolution of 8 or less computational cells per half-reaction length. However, it is generally agreed that when employing one-step Arrhenius kinetics, a resolution of at least 20 computational cells per half-reaction length is required in inviscid flow to sufficiently resolve the flow structure and obtain correct results for heat release profile, flame-shock interaction, detonation cell size, and so forth [16–18]. This is already challenging to achieve in large computational domains. Even though it is often argued that diffusive effects are negligible when simulating established detonations, they are of importance in fast deflagrations and the onset of detonations. Mazaheri et al. [19] showed that the inclusion of diffusive terms increases the required spatial resolution to approximately 25 to 300 computational cells (depending on the activation energy of the mixture) per half-reaction length. Powers and Paolucci [20] demonstrated that when using a detailed chemical mechanism instead of simplified kinetics, the demand for spatial resolution increases to an order of magnitude of approximately cells per half-reaction length. Due to the enormous computational costs, such “direct” simulations of DDT (even when using only simplified kinetics at 20 cells per half-reaction length) cannot be applied to large scale engineering simulations and are expected to remain impossible for years to come. Thus, the current state of the art of large scale accident simulations (e.g., in nuclear safety studies) is as follows [3, 4].(1)Determination of spatial hydrogen distribution.(2)Usage of empirical criteria to determine whether DDT can occur or not.(3)Simulation of the combustion with either a deflagration or a detonation code.
Step 1 can be performed with a lumped parameter code [5] or a CFD code [21, 22]. Step 2 is performed by applying empirical criteria [23, 24] which have been obtained from experiments in explosion tubes [25–29] to the geometry enclosing the hydrogen cloud. This step can only be performed with considerable uncertainty regarding not only the different scale and geometry (explosion tube with regular, periodic obstacles versus intricate three-dimensional building structure), but also the fact that virtually all criteria have been gained from experiments with perfectly homogeneous mixtures. In accident scenarios, however, mixtures are likely to include strong concentration gradients [30, 31]. Step 3 can be performed with a variety of CFD codes. While slower deflagrations can be simulated satisfactorily with commercial codes, simulations of fast deflagrations (where gas dynamics have a major influence) and especially simulations of detonations are usually performed more accurately with in-house codes that are generally not made available to the public.
In an OECD report prepared by a group of experts this “significant lack of numerical tools available to safety analysts” [32] has been criticized and the following requirements have been identified.(1)Development of approximate but reliable methods for simulating both flame acceleration and detonation in such a fashion that the simulation can be run within a single software framework.(2)Development of reliable combustion models that can be used to model DDT without the judgment and intervention of the simulator.
As the current situation is very unsatisfying, the present study was aimed at developing a CFD solver capable of simulating both deflagrations and detonations and especially the transition between both regimes. While the application of 3D computations at full reactor scale remains a long-term objective, this project goes a first step into this direction: the development of a solver to show the technical feasibility of simulating DDT experiments in explosion tubes in 2D without resolving the microstructure of the flow in the CFD grid. This forms an important prerequisite for the future simulation of flame acceleration and DDT in large, three-dimensional domains which will necessarily be performed on underresolved grids.
In a second step, this new solver is used to investigate the influence of concentration gradients on DDT formation. While many earlier studies considered a homogeneous mixture a worst-case scenario and consequently neglected the role of inhomogeneities, the influence of concentration gradients has only recently come into the focus of research [29, 33, 34].
The results of this project are shown in this paper.
2. Model Description
As the direct initiation [35] of a detonation is very unlikely, a detonation usually occurs after a turbulent flame (deflagration) has accelerated to a sufficiently high velocity [32]. This deflagration-to-detonation transition is a complex phenomenon including flame instabilities, interaction with turbulence, and gas dynamic phenomena. With all these effects occurring at very high Reynolds numbers, it is virtually impossible to resolve them completely in numerical computations. However, not all phenomena occurring on microscopic scale are necessarily relevant for an accurate simulation of macroscopic DDT events. Thomas [36] concluded from a comparison of experimental and numerical results that a reliable DDT model does not have to resolve all details of the flow. Instead, only correct turbulent burning rates, local density increase due to shocks and the capability of giving rise to detonations as a result of blast waves have been identified as necessary features of an accurate DDT model.
For the present work, the CFD code OpenFOAM [37] has been chosen as a basis for model development. One reason for this choice was that (contrary to many proprietary detonation codes) OpenFOAM has the built-in capability of dealing with unstructured grids—a clear advantage in view of the intended future application to sophisticated geometries. Another reason for the usage of OpenFOAM was that the code is free and open source and thus also the model developed in this project can be made available to the scientific community at no cost.
The density-based code developed under OpenFOAM solves the unsteady, compressible Navier-Stokes equations. All convective fluxes are determined using the HLLC scheme [38] with multidimensional slope limiters (“cellMDLimited” [37, 39]). This scheme is very suitable for the simulation of high Mach number compressible flow as it leads to much better shock capturing than the standard schemes used in most pressure-based codes like the PISO scheme [40]. As an example, the results of a 1D shock-tube calculation are compared in Figures 1 and 2. The initial data for this shock-tube problem is an ideal gas with molecular weight and specific heat ratio at , , for and , , for . The results shown have been obtained on an equidistant grid with spacing. For comparison, the analytical solution is displayed as a light gray line. It can be seen from Figure 1 that the PISO scheme not only predicts a wrong shock location (i.e., wrong propagation speed), but also is in general very dissipative and displays overshoots at discontinuities. This can be attributed to the nonconservative formulation of the Navier-Stokes equations which is inherent to the pressure-based scheme. With regard to the intended application to transonic reactive flow (including autoignition caused by shocks) this has to be regarded as a very critical issue. It can be concluded that standard pressure-based solvers are not suitable for the simulation of fast deflagrations and detonations.
Density-based solvers, especially Riemann solvers like the HLLC scheme employed in Figure 2, show a much better performance, producing accurate shock propagation speeds and far less dissipation at discontinuities. Thus it was decided to use the HLLC scheme as a basis for the DDT solver developed in this work. The only disadvantage of the scheme is that it does not work in very low Mach number flow. Therefore, the PISO scheme is also implemented and can be used to start computations in stagnant flow and then switch to the HLLC scheme once a combustion-induced flow has developed (see Section 3).
Realistic material properties for the reacting hydrogen-air mixture are obtained from the Chemkin database [41] and molecular transport coefficients are determined using the Sutherland correlation [42].
Combustion is described via a reaction progress variable [43]. corresponds to an unburned mixture, to a completely burned mixture. Within the context of Favre-averaging [44, 45], can be interpreted as the density-weighted probability of encountering burned gas at a particular instance of space and time. The transport equation of the reaction progress variable reads where the overbar denotes Reynolds-averaging and the tilde denotes Favre-averaging. The equation contains two source terms that account for deflagrative and detonative combustion, respectively. This concept is related to the simulation of autoignition in gas turbines or internal combustion engines as proposed by [46–48].
The deflagrative source term is modelled using the RANS version of the Weller combustion model [49], extended by a factor [50, 51] which takes quenching of turbulent flames into account [52]: is the density of the unburned mixture and the turbulent burning speed which is modelled as the product of the laminar burning speed and a flame wrinkling factor . The latter is obtained from a transport equation: Details about the expressions and describing the generation and destruction of flame surface can be found in [37, 49]. Due to the gradient ansatz, (2) leads to correct consumption rates even on underresolved grids. This means, even if the flame thickness is smeared out over several computational cells and thus thicker than the physical flame brush, the overall consumption rate is not affected [49, 52]. An adjustment of model constants to grid size or domain geometry is not required.
Experimental values for the laminar burning speed of hydrogen-air flames at standard conditions (temperature and pressure ) have been published by Konnov [53]. The dependence of flame speed on molar hydrogen fraction can be approximated as a polynomial [52]: In inhomogeneous mixtures it is essential to compute a correct flame speed in partially burned cells. Therefore, in the computation of the laminar flame speed, is not to be based on the actual hydrogen content, but on the mixture fraction , that is, the amount of hydrogen that would be present if the cell was completely unburned [54]. This is achieved by evaluating the hydrogen content (molar fraction of the hydrogen molecule) based on the mixture fraction (mass fraction of the hydrogen atom). The spatial distribution of the mixture fraction is described by a transport equation [54]:
The dependence of laminar flame speed on pressure and temperature can be approximated as [55]: Constant values of and have been used in this study.
The detonative source term in (1) accounts for autoignition effects. Autoignition occurs after the expiry of the local autoignition delay time . The autoignition delay time of a gaseous mixture is a function of local temperature , pressure , and mixture composition. In pure hydrogen-air mixtures, the mixture composition can be described using the mixture fraction (6). In order to avoid frequent recomputation of the local ignition delay time, a table of as a function of , , and is generated using Cantera [56] and the Ó Conaire mechanism [57]. This mechanism is valid for pressures up to 87 atm and temperatures up to 2700 K. The flow solver can access the table and query the local autoignition delay time in each computational cell.
An alternative means of modelling autoignition is the concept of a virtual radical species preferred by some authors [46, 47]. While in this study the concept of a tabulated ignition delay time is preferred, it can be shown that both models equivalently lead to the definition of a dimensionless variable describing the autoignition process [52]: The ignition variable reaches unity when the ignition delay time has passed (equivalently, it can be thought of the virtual ignition radical reaching a critical concentration). As long as the ignition delay time has not been reached yet (), there is no impact on the flow properties. Only if is reached, the mixture is ignited.
Due to the fact that autoignition in a DDT context can be triggered by shock-induced heating, a submodel is introduced that increases the accuracy of autoignition modelling on coarse grids. Tosatto and Vigevano [58] demonstrated that, while the average temperature in a computational cell can be high enough to trigger autoignition, it is important for detonation simulations to account for the fact that the shock causing the temperature rise might not have traversed the entire computational cell yet. Consequently a model that predicts autoignition of a computational cell based on average temperature and average pressure leads to incorrect results. Due to the large disparity of scales, that is, the shock being far too thin to be resolved on the computational grid, each computational cell is divided into two parts: in one part (volume fraction ) temperature and pressure are elevated to and ; in the other part (volume fraction ) the values remain at and . is defined by the highest value and by the lowest value that can be found in the surrounding computational cells. From consistency with the average pressure in the computational cell, the volume fraction can be determined: and are computed subsequently from gas dynamic shock relations for an ideal gas with heat capacity ratio [59]: When and are known, (10) can be solved for the shock Mach number and (11) yields the according temperature ratio . The temperatures can be determined by requiring consistency of and with the cell-averaged temperature :
The resulting representation of a computational cell is given in Figure 3. While the presence of a shock alters temperature and pressure, the mixture fraction is not affected. The ignition delay time is evaluated separately on each side of the shock: It is important to note that this model does not only work in computational cells where a shock is present but can also be applied to the entire computational domain. The reason is that in the case of low pressure differences, , the temperature rise computed from (10) and (11) is quasi-identical to the temperature rise gained from an isentropic compression:
In each grid cell, the process of autoignition is evaluated separately on both sides of the discontinuity. Transport and mixing effects are accounted for by solving transport equations for and :
If the critical value of is reached on one side of the discontinuity (i.e., or ), only the corresponding volume fraction of a computational cell is ignited. Consequently the autoignition source term in (1) can be formulated as
Here, represents the current time step and represents the Heaviside function:
The Heaviside function activates only the part of the computational cell in which the local ignition delay time has expired. This is achieved by weighting the volumetric source term by the volume fraction of either or .
3. Experimental and Numerical Setup
The model has been tested against experimental results gained in a closed rectangular channel of length , height , and width . The channel is equipped with flat plate obstacles (thickness ) of height spaced at a distance of from each other. The first obstacle is placed at from the front plate where a spark plug ignites the mixture. The last obstacle is placed at and the remaining part of the channel is unobstructed (see Figure 4). The obstacle blockage ratio BR is determined by the obstacle height :
The top wall of the channel is equipped with 42 UV sensitive photodiodes and 6 pressure transducers operated at a sampling rate of . Another pressure transducer is mounted head-on in the center of the end wall (). A flame position versus time ( versus ) correlation is obtained from the photodiode measurements. The flame velocity between two subsequent photodiodes is calculated by applying a first order derivative: Here, and represent the time at which the flame passes the photodiodes located at and , respectively. The same procedure is applied for evaluating the flame velocity in the numerical simulations.
Within the channel defined vertical concentration gradients can be generated. The overall amount of hydrogen is controlled via the partial pressure method. First, the air-filled channel is partially evacuated. Then, hydrogen is injected through several nozzles located at the top wall. The injection velocity is constant due to a choked nozzle upstream of the point of injection. The injection time defines the amount of hydrogen injected. Subsequently there is a defined time interval (waiting time ) during which diffusion takes place. Due to the strong density difference between hydrogen and air, a defined vertical concentration gradient is achieved while horizontal concentration gradients remain negligible. Finally the mixture is ignited by a spark plug. For a more detailed description of the experimental setup, the procedure of hydrogen injection, and mixture generation it is referred to the publications of Vollmer et al. [29, 60].
In order to determine the hydrogen distribution before ignition for many different hydrogen/air ratios and different waiting times, numerical simulations of the injection process have been conducted. Exemplary results for local hydrogen mole fraction over channel height at a waiting time (the strongest gradient under investigation) are shown in Figure 5. For waiting times the mixture can be considered as homogeneous. The results of the injection simulations are stored as polynomials (hydrogen content versus channel height) which are used as initial conditions for the combustion simulations.
In the two-dimensional combustion simulations presented in the following section the channel is discretized with a uniform, rectangular grid of 2 mm grid spacing. Test runs showed that this resolution is the minimum resolution required to achieve grid independence with respect to the location of DDT. On coarser grids, DDT occurred mostly later or not at all. On finer grids, the location of DDT did not vary any more. However, at higher grid resolution, pressure peaks still got a little sharper. This should be kept in mind for the interpretation of the pressure plots shown in this paper.
Initially the fluid is at rest at a temperature of 293 K and a pressure of 1.01 bar. The boundary conditions are defined as adiabatic no-slip walls. Turbulence is modelled using the --SST model which is known for its good performance for both free-stream jets and wall-bounded flow [61, 62]. The initial hydrogen distribution either is homogeneous or corresponds to a concentration gradient of waiting time (see Figure 5).
Ignition is modelled by patching the site of ignition at with a burned mixture (, see Figure 4). The initial turbulence is vanishingly small and consequently equals unity so that follows from (3). This means, the flame starts to propagate at laminar flame speed. However, turbulence is quickly generated by the flow itself so that the flame starts to accelerate. Test runs showed that the actual choice of initial turbulence variables is insignificant as long as is ensured. As the HLLC scheme gets unstable in the incompressible limit where no coupling between pressure and density exists, the first few time steps are calculated with a pressure-based solver [37] using the PISO scheme [40]. Before the flame reaches the first obstacle, the combustion-driven flow is usually strong enough to switch to the HLLC scheme that enables better shock capturing. Test runs showed that the transition between both schemes is smooth if it occurs while the maximum Mach number in the flow is in the range of .
4. Results and Discussion
Experimental and numerical results for a homogeneous case with 15% hydrogen (volumetric) and blockage ratio are shown in Figure 6. It can be seen that the agreement between experiment and simulation is very good. The flame velocity rises continually in the obstructed part of the channel (). This can be attributed to the mutual amplification of combustion-induced expansion and turbulence generation due to interaction with obstacles. Shortly after passing the final obstacle the flame speed reaches a maximum and then decreases slowly. At the flame comes to a nearly complete rest before it accelerates again. This can be explained as follows: after passing the final obstacle, turbulence generation is diminished so that decelerating effects like friction outweigh the accelerating ones. The flame continuously gets slower. Simultaneously, while the flame has been consuming fresh gas, it generated pressure waves and displaced the unburned gas into the positive direction. Shocks were generated that propagated towards the end wall from where they are being reflected. These reflected shocks now generate fluid flow in negative direction. When the leading, backwards-running shock reaches the flame (this happens at ), negative flow velocity and positive burning speed nearly cancel out so that the resulting net propagation velocity approaches zero. However, as there is still unburned gas in front of the flame, it recovers and accelerates again. The maximum flame propagation velocity of approximately indicates that no DDT occurred and the combustion process remained entirely deflagrative.
Figure 7 shows the results for a mixture that contains an average hydrogen content of 15% as well, but with a vertical concentration gradient as shown in Figure 5. All other parameters are kept identical. In the early acceleration phase the flame velocity increases continually. This is in good agreement with the experiment. Then the flame is decelerated for the first time, due to a first shock front reflected from the end wall. The difference in the experiment can be attributed to the different ignition process: the spark generated by the spark plug in the experiment is considerably smaller than the ignition patch used in the simulation which is limited by the grid resolution. Thus the initial pressure wave generation in the experiment might be a little different from the initial pressure rise caused by the ignition in the numerical simulation. At the end of the obstacle region the flame speed peaks and then loses some driving force but eventually recovers. Although there is a considerable velocity difference between experiment and simulation in the unobstructed part, the final velocity is nearly the same. The pressures recorded by the sensor in the end wall reach extremely high values close to 120 bar (see Figure 8). In the homogeneous case, for comparison, the maximum pressure is in the range of 10 bar. The reason for the extreme pressure rise in the inhomogeneous case is revealed in Figure 9 where the temperature and pressure distribution in the rear part of the explosion channel () is displayed.
At the flame approaches the end wall. Due to the inhomogeneous fuel distribution the flame is highly asymmetric and propagates mainly in the upper part of the channel. A leading shock has already been reflected from the end wall and moves towards the propagating flame. At it reaches the flame. From this point onwards the flame burns into a precompressed mixture where the heat release rate is increased due to the increased density and increased laminar burning velocity (see (2) and (7)). The increased reaction rate leads to a strong pressure rise and causes an explosion at . A radial detonation wave emanates from the explosion center and ignites the gas over the whole channel height. The newly formed detonation front runs towards the end wall where it causes an enormous pressure rise. This DDT mechanism has been suggested as one possible explanation for the high pressure loads observed in the experimental work of Eder [63]. In Eder’s work, high pressure loads on the end wall of an explosion channel have been observed, but the flame velocity measurements indicated only a fast deflagration, not a detonation. As in the present simulation, the DDT in Eder’s experiments obviously occurred so late (behind the final photo diode) that the DDT was not identified as one; only the high pressures on the end wall gave rise to speculation. Recent experimental investigations of Boeck et al. [64] support the conclusion that a DDT mechanism as identified in Figure 9 is responsible for the high pressure peak.
Due to the limited spatial resolution the present simulation does not resolve the interaction with the boundary layer. Moreover, it does not capture the shock-flame interaction in such a detailed manner as previous numerical studies on highly resolved grids (e.g., [12, 13]). Nevertheless the model is able to correctly predict the consequence of the backwards-running shock hitting the flame: an increased reaction rate due to precompression and intensified mixing which consequently triggers DDT.
From the pressure records in Figure 8 it can be concluded that there is a slight difference between experiment and simulation: the initial pressure rise in the simulation at (caused by the reflection of the leading shock) does not appear in the experimental record. Due to the highly nonlinear dependence of ignition delay time on temperature and pressure, the higher propagation velocity in the experiment (Figure 7) is obviously sufficient to cause a strong autoignition quasi-instantaneously when the leading shock reaches the end wall. The resulting pressure load on the end wall, however, is nearly the same in experiment and simulation.
Increasing the hydrogen content leads to an earlier occurrence of DDT. At a hydrogen content of 25% (again with a concentration gradient as described in Figure 5) it can be seen from Figure 10 that the flame velocity rises continually to approximately 1000 m/s in the obstructed part of the channel and then suddenly jumps to 2500 m/s and finally relaxes to approximately 2000 m/s.
This is a clear indication for the occurrence of a DDT with an initially overdriven detonation decaying to a Chapman-Jouguet detonation. The large fluctuations in the experimental velocity after the onset of DDT can be explained by small measurement errors in flame arrival time having a relatively large effect when the derivative (19) is applied to the data. Using only the - diagram (Figure 11) as it is common in most publications does not reveal this difference.
The DDT process occurring in this case is visualized in Figure 12. At , the flame approaches the final obstacle. The curved shock in front of the flame is reflected from the bottom wall by forming a Mach stem. At , autoignition occurs behind the Mach stem. At , the oblique shock hits the upper obstacle which initiates a second autoignition event. From there a circular detonation emanates and unites with the autoignition front from the lower part of the channel. While the detonation front moves through the gap between the obstacles into the unburned gas (), the opposite front of the reaction wave (“retonation wave” [65]) runs backwards and consumes the remaining fresh gas in the lower part of the channel. It is important to note that the two autoignition kernels in Figure 12 are both well ahead of the flame but occur due to different reasons: the one on the bottom wall is due to shock compression ahead of the flame while the one on the upper wall occurs only due to reflection of the shock from the upper obstacle.
Another simulation with only six obstacles showed that the final obstacle was not necessary to achieve DDT. Instead, the autoignition occurring behind the Mach stem at is sufficient to trigger DDT and is only amplified by the second autoignition event occurring on the upper obstacle. At lower fuel content (20% H2) however, the seventh obstacle is required to obtain a DDT.
It is interesting to note that the first autoignition in Figure 12 occurs at the bottom wall where the mixture is leanest. This phenomenon can be explained by taking a closer look at the shock propagation: the leading shock approaches the final obstacle at a constant speed of . Near the bottom wall the hydrogen content is 7% (see Figure 5) which results in a local speed of sound of . Thus, in the near vicinity of the bottom wall, the Mach stem can be seen as a normal shock propagating at Mach number . For fresh gas properties and the normal shock relations (10) and (11) yield the postshock state and . Near the top wall where the mixture is rich (45%, see Figure 5) the local sound speed has a value of . This corresponds to a Mach number of and yields a postshock state of and . The corresponding ignition delay time is by orders of magnitude higher than on the bottom wall. Thus, autoignition near the top wall can only be achieved via shock reflection from the upper obstacle.
Another important finding is that although DDT occurs much earlier in the case of 25% hydrogen, the pressure loads on the wall are higher in the case of 15% hydrogen. This is due to the different mode of DDT. At 25% hydrogen (Figure 12), the detonation has sufficient space to expand. At 15% hydrogen (Figure 9), shock and flame interact very close to the end wall and the overdriven detonation emerges from a preshocked mixture. As there is no space to expand, the overdriven detonation hits the end wall without losing much of its strength.
While both simulations and experiments in this configuration confirm the trend that (at equal average hydrogen content) a concentration gradient increases the DDT tendency, this must not be understood as a general conclusion. In other configurations the opposite effect can occur. For example, if the blockage ratio is increased from 30% to 60% it has been observed that the probability of DDT decreases. The reason for this phenomenon can be seen in Figure 13. At , the flame approaches the obstacle. Some pressure waves have already been generated and passed the obstacle. The constriction formed by the obstacle acts like a Laval nozzle behind which an area of supersonic flow (visualized by the pink line representing the contour) is generated. The obstacles cause the pressure waves to be reflected on 60% and to pass only on 40% of the cross-sectional area. Due to the gas dynamic constraint formed by the line, the total mass flow through the orifice between the obstacles is limited. At , a shock front is reflected from the obstacle which leads to a shock propagation into negative direction. This process is repeated with an even stronger shock between and . As the parts of the shock that are reflected from the upper and lower obstacle unite, a backwards-running shock over the whole channel height is formed. In the trailing flow behind this shock an area of negative flow velocity causes the approaching flame to decelerate, while the part of the shock that passed in between the obstacles to the right continues to propagate at nearly unaffected speed. Thus a high blockage ratio destroys the coupling between shock and flame which existed before approaching the obstacle. At , the flame finally manages to enter the region of high flow speed that prevails between the obstacles. It is convected into the next cavity where strong shocks are generated and the process described repeats.
These numerical results explain the experimental observation by Vollmer et al. [29] that a concentration gradient can either increase or decrease the tendency towards DDT with the decisive factor being the obstacle geometry. If the obstacles are too large, they can lessen the DDT tendency as the majority of the strong pressure waves causing DDT are blocked. This is especially valid for an inhomogeneous mixture as shown in Figure 13, where the flame mainly burns in the upper part of the channel. In this case the obstacles are more obstructive than in a case with homogeneous mixture where the flame can be expected to propagate through the center of the channel where no obstacles are present.
Another question that has been addressed with the newly developed solver concerns the pressure loads that are caused by a steadily propagating detonation front, that is, after the occurrence of DDT. Therefore a look is taken at a detonation propagating in an unobstructed channel. First, this is demonstrated for a homogeneous mixture with 25% hydrogen. Pressure records are taken from the top and the bottom wall of the channel while a detonation passes. Test runs showed that the axial location of the pressure sensors did not influence the result any more as soon as a steadily propagating detonation was achieved. The result is shown in Figure 14. As the detonation front is nearly planar, the pressure records from the bottom and the top wall are virtually simultaneous. Upon arrival of the detonation front the pressure jumps to approximately 18 bar. The following expansion lets the pressure decrease slowly.
A completely different picture is found for a case with the same average hydrogen content, but a strong concentration gradient (Figure 15). As the leading shock is curved, it reaches the pressure sensors on the top wall earlier. They show maximum values of approximately 15 bar. On the bottom wall, however, a pressure of nearly 38 bar is reached. This is especially striking as the hydrogen content on the lower wall is only 7% and a homogeneous mixture with 7% hydrogen is basically nondetonable. Here, however, the lack of fuel does not lead to lower, but to higher, pressure loads. Again, the reason for this seeming paradox can be found in the particular structure of the leading shock front: on the bottom wall it is reflected via a Mach stem. Due to the lower speed of sound this causes a higher pressure rise on the bottom wall than on the top wall. After a short decline of the pressure, a second pressure rise is observed on both walls. This is due to secondary reflections of the leading shock that can be seen in the pressure field in Figure 16. Behind the secondary reflections the pressure equalizes so far that it drops simultaneously on the bottom and the top wall.
The results demonstrate that the pressure loads caused by a detonation in an inhomogeneous mixture can be considerably higher than in a homogeneous mixture of the same hydrogen content. Moreover, the location of the highest impact can be in fuel-lean regions. Further calculations showed that even if the hydrogen content on the bottom wall is reduced to zero, the maximum pressures observed there can still exceed those of the homogeneous mixture: the concentration gradient only needs to be strong enough to form a Mach stem. A simple method for predicting whether a detonation front in an inhomogeneous mixture develops a Mach stem can be found in [66].
5. Summary, Critical Analysis, and Outlook
Motivated by the current lack of suitable tools for DDT-related safety studies [32], this paper presented a newly developed solver able to simulate flame acceleration, deflagration-to-detonation transition, and detonation propagation within a single run. The target was not to obtain detailed insight and maximum accuracy of the complex interaction between flow and reaction on microscopic scale, but to obtain a tool for engineering purposes that works on comparatively coarse grids and enables numerical safety studies at acceptable computational costs. The applicability to coarse grids is achieved by the inclusion of subgrid models. The agreement with experimental results is very good and the simulation gives additional insight into phenomena which cannot be easily observed in experiments. Although the simulations presented do not resolve all details of the flow, they are able to capture fundamental phenomena known from highly resolved simulations and experiments (e.g., [12, 13, 36]): DDT due to shock compression/Mach stem formation ahead of the flame, DDT due to shock reflection from obstacles, and DDT due to shock-flame interaction.
It has been found that concentration gradients, which are likely to occur in accident scenarios, can have a considerable effect on the nature of flame propagation. Depending on the enclosing geometry, the presence of a concentration gradient can decrease or increase flame propagation velocities, the probability of DDT, and the pressure loads associated with it. Thus, existing safety criteria developed for homogeneous mixtures can be inaccurate and nonconservative. Neither does a homogeneous mixture pose the highest threat regarding the probability of DDT nor does it cause the highest pressure loads. Due to gas dynamic phenomena within inhomogeneous mixtures, fuel-lean regions can be more DDT-prone than stoichiometric or rich regions. This should be taken into account in future safety studies.
However, although the general agreement with experiments is good, it has to be kept in mind that all results have been gained on relatively coarse grids, without resolving the induction distance between shock and reaction in a detonation front. Boundary layers are not resolved either and they are known to have an effect on the onset of detonations. Moreover, all the results presented in this study have been obtained on 2D grids. Therefore, the authors have deliberately chosen a geometry which is relatively wide (300 mm) compared to its height (60 mm). Nevertheless transversal waves are expected to play a role in flame acceleration and the onset of DDT. While this project has already finished, a follow-up project has started where 3D simulations are conducted [67] and also simulations in large, complex geometries with the aim of reproducing realistic accident scenarios. Approaches are being developed to use even coarser grids by applying subgrid models not only to shock propagation as shown in this paper, but also to deflagrative flame propagation.
The advantage of the solver developed and its implementation in OpenFOAM is that it is not limited to structured grids and thus can be applied to intricate geometries using unstructured grids as well. The solver and its source code are made freely available to the public [68].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This project is funded by the German Federal Ministry of Economics and Technology on the basis of a decision of the German Bundestag (Project no. 1501338) which is gratefully acknowledged. The authors would like to thank Oliver Borm for sharing a Riemann solver code within the framework of OpenFOAM which provided a valuable basis for this study.