Experimental and Numerical Investigation of Gaseous Detonation in a Narrow Channel with Obstacles

: Gaseous detonation propagation in a thin channel with regularly spaced cylindrical obstacles was investigated experimentally and numerically. The wave propagation with substantial velocity deﬁcits is observed and the details of its propagation mechanism are described based on experimental measurements of the luminosity and pressure and on three-dimensional ﬂow ﬁelds obtained by numerical simulations. Both experiments and simulations indicate a signiﬁcant role of shock–shock and shock–obstacle interactions in providing high-temperature conditions necessary to sustain the reaction wave propagation.


Introduction
Quasi-detonation propagation and flame acceleration in ducts with porous media or obstructions has been studied over many years. The main velocity propagation modes have been identified and classified depending on the gas mixture composition, stoichiometric ratio, initial pressure, porosity, pore size, and blockage ratio. Low-speed deflagration, highspeed deflagration, quasi-detonation, and Chapman-Jouguet (CJ) detonation regimes have been found and described, and empirical criteria for flame acceleration and transition from one regime to another have been proposed. In particular, Refs. [1][2][3][4][5][6] report on experimental observations of the multitude of propagation regimes in reacting flow in porous media and the role played by the material of the porous medium, of the type of reacting gas, and of the mixture initial pressure in the observed regimes. Gaseous detonation propagation in channels with various types of obstacles and with different levels of the channel obstruction was investigated in [7][8][9][10][11][12][13][14][15]. Despite the large amount of experimental and theoretical data, the physical mechanisms of wave propagation for different regimes remain unclear and would benefit from further research.
In a broad sense, high-speed reaction waves in obstructed channels are driven by local shock wave compression and adiabatic compression resulting in spontaneous ignitions of unburnt products at elevated thermodynamic conditions. In the case of a shock compression, the hot spots in the porous medium may be formed by single or multiple reflections of the leading shock wave, shock diffractions, and focusing effects within the porous matrix as well as interactions of shock-headed transient gas jets from the adjacent pores. It was shown in prior work that the normal shock reflection mechanism by itself cannot produce the auto-ignition of the mixture at realistic length scales in a porous bed and cannot ensure the self-sustained propagation of quasi-detonation near the limits [15,16]. From experimental and theoretical points of view, it is interesting to elucidate the dominating ignition mechanisms of quasi-detonation at various conditions and to identify critical conditions for transition between different mechanisms. At present, there is a lack of theoretical understanding of the relative importance of various physics at play for given experimental conditions.
The boundary layers, shock reflections, and diffraction that arise at the interfaces between the gas and obstacles or walls lead to momentum losses and thermal exchanges between the gas and the obstructions. Originally, these effects have been discussed in some detail in the seminal work by Zel'dovich [17], in [18], and in later works (e.g., [19][20][21]). The basic outcome of such analyses is that the obstacles can produce a significant decrease in the propagation velocity (called "velocity deficit") of the detonation compared to the ideal case without losses. An important feature of such non-ideal detonations is their ability to propagate with velocities ranging from near sonic value in the ambient gas to the ideal CJ value [3,4,6,8,22,23], the actual velocity depending on the mixture composition as well as on the nature of losses.
Three-dimensional numerical simulations carried out in this work consider the reaction wave in the same channel with a matrix of cylindrical obstacles as in the experiment and with the same mixture of acetylene and oxygen diluted with argon. In addition to the structure of the wave, the heat and momentum losses to the walls and obstacles are analyzed. Closely related simulations in a two-dimensional setting have been carried out in [23,24], however neglecting such factors as heat losses to the walls.
The remainder of the paper is organized as follows. In Section 2, we describe the experimental setting and measurement results. In Section 3, we explain the main modeling elements and simulation methods. Simulation results are then compared with the experimental data and we provide further details on the observed dynamics of detonation. Discussion of the results and concluding remarks are given in Section 4.

Experiments
The experiments were carried out in a planar channel with a matrix of cylinders fixed between the side walls as shown in Figure 1. The wave front in such a channel can be approximated as planar. It allows one to observe ignition phenomena inside a single pore and to explore the structure of the detonation and deflagration waves in the system as they interact with the obstacles.  Measurements were carried out in a test volume with dimensions 8 × 107 × 500 mm in which a matrix of cylinders was installed with diameter 8 mm, arranged in a checkerboard order. A spark plug initiated a flame in a small pre-chamber. It subsequently transmitted to a quasi-detonation mode in the porous layer. The propagation was observed through the frontal transparent wall along the entire length of the test volume using two high-speed cameras PHOTRON FASTCAM SA-X2 and PCO DICAM PRO equipped with dual narrowband filters with λ max = 430.4 nm, ∆λ 0.5 = 2.6 nm, and λ max = 430.6 nm, ∆λ 0.5 = 2.6 nm (photoluminescence of CH radicals). Various mixtures were used, including stoichiometric acetylene-oxygen diluted with 70% argon. Initial pressure was varied from 0.01 to 0.04 MPa. Commercial grade acetylene, oxygen, and argon of 99.9% purity were used for mixture preparations. A pressure meter controlled the initial pressure of the mixture with an accuracy of ±0.4 mm Hg. Sixteen ion gauges and three pressure transducers measured the velocity and pressure of the shock and reaction fronts. Two lines of ion sensors (eight in each line) provided measurements of a quasi-detonation velocity in different parts of the test volume along the direction of the wave propagation.
The propagating wave was recorded with a frame rate of 324,000 frames per second. In Figure 2, we display the video frames of a small zone in which the wave velocity is already established. The arrows indicate the direction and distance traveled by the bright region in the flame front. Based on the recordings, the local reaction-zone velocities around the cylinders were determined. The velocity varied from 400 m/s to 1500 m/s. If we define the average quasi-detonation velocity as the velocity of the leading shock wave heading a quasi-detonation structure, it is found to be V = 675 m/s. Note that the pressure and temperature behind the normally reflected shock wave at initial pressure 0.02 MPa are 0.35 MPa and 907 K, respectively. Based on these values, one can determine the ignition time and induction-zone length using a detailed kinetic mechanism from [25]: τ ind = 1.17 ms, L ind = τ ind V = 790 mm. Consequently, in the time between the initial shock compression and auto-ignition of the mixture in the porous layer, the leading shock wave would have traveled the distance approximately equal to 158 pore sizes. In contrast, the pressure and ion current measurements show that the induction-zone length is less than one pore size. This result indicates that the self-ignition of the mixture due to normal shock reflection is insufficient to maintain the propagation of quasi-detonation in a porous medium in agreement with similar conclusions reached in earlier works [16].
On the other hand, if we assume that the lead shock velocity is close to the maximum value V = 1516 m/s measured in the experiments, then pressure and temperature behind the lead shock would be 0.56 MPa and 1770 K, and the induction time and distance would be τ ind = 0.5 µs and L ind = 0.76 mm, respectively. This value of the induction-zone length is close to the pore size ( Figure 2 , frame 1) in agreement with both experimental observations and numerical simulations. This is an indication that the detonation propagates locally with its normal velocity while the average value across the channel is smaller due to both geometric factors and the extra deficit coming from thermal and momentum losses.  The arrows indicate the propagation direction. The angles are used to identify the reaction wave position and subsequent calculation of the local wave velocity. The luminosity is seen to be large behind the reflected wave as it collides with the obstacle (e.g., the middle obstacle in Frame 3) and also in the region behind the obstacle where the two shock waves diffracted around the obstacle collide and give rise to a new shock wave that moves in between the subsequent obstacles (Frames 8 and 9).

Simulation Results
Numerical simulations were carried out using the commercial solver ANSYS Fluent ® [26] and use the same 3D geometry with the same dimensions and material parameters as in the experiments. The governing compressible and reactive Reynolds-averaged Navier-Stokes (RANS) equations account for: viscous effects, turbulence, species and heat transport, and chemical reactions based on a global scheme of conversion [27]. Additionally, the wall heat transfer is accounted for in order to estimate its role in the observed dynamics. This is accomplished by setting the side-wall boundary condition at a fixed temperature while the obstacles have been kept adiabatic.
Specifically, the governing equations used are given as follows (see, e.g., [27] for general discussion of equations of reacting flows). The continuity equation is given by ∂ρ ∂t + ∇ · ρu = 0 with ρ the mixture density and u the mean velocity vector. The momentum equation is given by where p is pressure, µ is the mixture dynamic viscosity, the Reynolds stresses R ij = −ρu i u j are given in terms of the velocity fluctuations u i and are closed as with turbulent viscosity µ t and turbulence kinetic energy k, which are calculated based on equations below for k and turbulence dissipation rate .
The energy equation is written as ∂ρE ∂t where T T re f c p,j dT is the total enthalpy, T re f = 298.15 K, and c p,j is the specific heat capacity of j-th species, Y j is the mass fraction of species j, and τ e f f = µ e f f δ ij is the deviatoric stress tensor causing viscous heating in (3), where κ is the molecular thermal conductivity and κ t = c p µ t /Pr t is the turbulent thermal conductivity in which we take Pr t = 0.85. Also, j the enthalpy of formation and R j the rate of creation of species j. The molecular weight of species is M j .
The standard k − turbulence model is given by the equations for turbulence energy k and turbulence dissipation rate [28]: Here, is the generation of turbulence kinetic energy by the mean-velocity gradients. The influence of compressibility on turbulence is accounted for by [29]. The turbulent viscosity is given by µ t = ρC µ k 2 / , and the model constants are taken as C 1 = 1.44, C 2 = 1.92, C µ = 0.09, σ k = 1.0, and σ = 1.3.
The chemical reactions are described globally by a one-step forward reaction C 2 H 2 + 2.5O 2 + xAr → 2CO 2 + H 2 O + xAr, thus the gaseous mixture consists of five components, all considered as ideal gases. In all the simulations, this stoichiometric mixture of acetylene and oxygen is used diluted with 70% argon (x = 8.17). The Arrhenius rate expression is used with pre-exponential factor A = 3.655 × 10 10 (units consistent with the rate equation below) and activation energy E = 1.256 × 10 8 J/(kmol). The rate of reaction is then given by where the diffusion flux is calculated as ∇T T , D i,m is the mass diffusion coefficient of species i in the mixture, Sc t (taken equal to 0.7) is the turbulent Schmidt number, µ t is the turbulent viscosity. Generally, turbulent diffusion dominates over the molecular diffusion in our system. In the last term, D T,i is the thermal (Soret) diffusion coefficient.
As initial conditions, we have the driving section from x = 0 to x = 6 mm (almost touching the leftmost cylinders, see Figure 3(a1)) which was filled with pure CO 2 at temperature T d = 3000 K, pressure p d = 20 atm, and zero velocity. The rest of the domain was filled with the stoichiometric mixture of acetylene and oxygen diluted with argon at pressure p 0 = 0.2 atm and temperature T 0 = 300 K, and the gas was also at rest. As the computations begin, the high-pressure gas from the driving section moves towards the reactive mixture driving a shock wave into it that is capable of initiating chemical reactions. The chemical energy released in these reactions can sustain the lead shock in which case we have a self-sustained reaction wave of a detonation or quasi-detonation type. We note that this mode of initiation differs from the one used in experiments and was chosen to save computational time, which is substantial (taking about a week on a modern 16-core workstation) considering the three-dimensionality and complexity of the problem. As boundary conditions, we impose the no-slip conditions at all solid boundaries, including all the side walls and the surfaces of the cylindrical obstacles. For turbulence computations, the standard wall functions [30] are used in order to match the boundary layer region with the fully turbulent outer region, thus requiring no numerical resolution of the boundary layers. Wall y * is defined in terms of turbulence properties in the near-wall region as y * = ρC 1/4 µ k 1/2 p y p /µ, where y p is the distance from the centroid of the walladjacent cell to the wall, and variables with subscript p indicate the corresponding values at y p . If y * < 11.2 at the wall-adjacent cells, the laminar stress-strain relationship is used, which can be written as u * = y * . At y * > 11.2, the log-law u * = ln(Ey * )/κ is employed with κ = 0.4187 the von Kármán constant and E = 9.793 an empirical constant.
The wall heat transfer is accounted for by setting the side walls at a fixed temperature of T w = 300 K thus allowing for significant heat losses from the reaction zone and burnt products to the walls. The Reynolds analogy between momentum and energy transport allows for a similar logarithmic law for mean temperature [31]. It employs the linear law for the thermal conduction sublayer and logarithmic law for the turbulent region where turbulent effects dominate over conduction. In the near-wall region both heat conduction and heating due to viscosity can be important in compressible flows and are included in the boundary treatment. The species transport is treated in a manner similar to heat transfer using the wall functions in the boundary layers [31,32].
In addition, boundary conditions for k and are required. The turbulence kinetic energy satisfies the homogeneous Neumann condition, ∂k/∂n = 0 at all walls. The value of is not computed at the walls, instead its value at the wall-adjacent cell is used in terms of k p as p = C 3/4 µ k 3/2 p /(κy p ). The cylindrical obstacles were considered thermally insulated, which is justified to some extent by them not being exposed to the outside low-temperature conditions, unlike the side walls. Obviously, some energy would be lost through the cylinders, however, proper accounting for such losses would require solving a heat conduction equation in the interior of the cylinders, which is neither feasible nor necessary for the purposes of the present study.
The equations were integrated using an explicit first order upwind algorithm with AUSM flux type with a Courant number 0.5. Even though this algorithm clearly lacks in accuracy and is rather dissipative, it was chosen because of its better numerical stability compared to second-order methods available in Ansys Fluent. More refined computations will be needed for improved predictions of the flow, however for this first study the present choice is adequate. The computational mesh was refined near the obstacles and the walls in order to better resolve near-wall regions to the greatest extent possible (see Figure 5). Furthermore, grid adaptation based on pressure gradients was included in order to improve shock resolution. The agreement with experiment that we have seen in the simulations is reassuring, in that at least some of the global characteristics of the wave are calculated reasonably well. Thus, the present simulations can be considered as a step toward realistic simulations of this complex three-dimensional reacting flow field. It is evident from both experiments and simulations that the most important aspect of the process of detonation propagation in the porous space that is at the heart of the propagation mechanism of the wave is the presence of multiple collisions of the shock waves with each other and with the obstacles. These collisions are seen to create localized hot regions in which the reactions proceed very rapidly, consuming the reactants and strengthening the reflected shock waves. The reaction zone is seen to be highly nonuniform and rather far from the standard detonation structure. The incident shock is too weak to sustain chemical reactions, but the reflected waves are capable of burning the mixture efficiently and thus help sustain the overall wave propagation. Since these reflected waves propagate in transverse as well as in the axial direction locally, on the average we obtain a wave that propagates in the x direction with a velocity substantially lower than the CJ velocity for a given mixture.
Analysis of the simulation results indicates that the structure and dynamics of the propagating wave is substantially in agreement with experimental observations in the present work as well as in prior research on this problem. Namely, the wave propagation is controlled by the geometry of the obstructed space and is driven by shock collisions leading to substantial temperature increase locally that accelerates the heat release process and thus serves to strengthen the shock waves.
In particular, we see in Figures 3 and 6 that upon the diffraction around the cylinder, the shock waves collide, their reflection slowing down the flow locally, but giving rise to a strong shock wave that moves in the direction of the free space in between the two following obstacles (Figure 3(a3)). This strong shock then collides with the next obstacle in its path and generates a reflected shock wave and two shocks that diffract around the obstacle. The process repeats (see also the supplementary videos). Thus, the reaction zone is distributed over the space of a couple of cylinders with a highly non-uniform heat release that is essentially controlled by the shock-shock and shock-wall interactions. We also point out that the wave front undergoes some transverse deformations with the upper part of the front lagging initially behind the lower part and then catching up. This may indicate the development of a transverse structure akin to a transverse shock wave in cellular detonations. It may also indicate that the wave may be undergoing a transition to a different regime with a lower velocity. Revealing the actual final stage will require more refined and longer-time simulations. Figure 6. Simulated profiles of mass fraction of acetylene and temperature in the center-plane between the side walls at times: t 1 = 1.81 × 10 −5 s, t 2 = 5.18 × 10 −5 s, and t 3 = 7.70 × 10 −5 s corresponding to the three rows: (a1-a3)-mass fraction of C 2 H 2 ; (b1-b3)-temperature.
One may expect that some pore-space configurations act more favorably in terms of the efficiency of the wave propagation. To the best of our knowledge, such a topology optimization study in the realistic three-dimensional setting is lacking in the literature, and it would be an interesting problem to explore. That said, we mention the recent work [24], which considers obstacles of various shapes in two-dimensional simulations of detonation and deflagration-to-detonation transition (DDT) in a channel with obstacles. The authors find noticeable effects of the obstacle shape on the DDT and detonation propagation.
In Figure 7, we show the computed values of the heat flux to the walls at two different times, corresponding to the second and third rows in Figures 3 and 6. One can see that inside and nearby the reaction zone, the magnitude of the flux is on the order of q w ∼ 5-10 MW/m 2 . This compares well with the fluxes in rotating detonation engines reported in [33]. Somewhat lower values around 1 MW/m 2 were reported in the original measurements of the wall heat fluxes in detonation engines in [34]. The latter authors also estimated the relative significance of these losses compared to the power density of chemical energy release and concluded that the losses can often be neglected. However, more recent research such as [33] indicates that these losses and heating of the channel walls may be important and deserve more detailed and careful study. Fluxes on the order of 10 MW/m 2 or more were reported in [33]. In Figure 4, we show the experimental (a) and numerical (c) records of pressure at several different locations in the channel. We can see the characteristic spikes due to the multiple reflected and diffracted waves passing over the gauge locations. The maximum pressure recorded in the experiment in this figure is just under 3 atm. The numerical result in (c), even though qualitatively in agreement, shows larger values of pressure around 10 atm early on at position 25.3 mm from the left end wall and about 5 atm at 63.4 mm distance. The early high pressures are likely due to the slight overdrive still present from the initiating section. More generally, the reason for the larger value of the average pressure may be the fact that at initial pressures in the mixture near p 0 = 0.2 atm, the wave can propagate either as a low-velocity detonation (LVD) with a velocity about 600-700 m/s or as a quasi-detonation with a velocity about 1000-1100 m/s. In nearly identical settings in experiments, one may observe either one or the other of these regimes. In the particular simulation reported here, the wave propagates at the velocity about 1050 m/s as shown in Figure 4d. Figure 4b shows the experimental wave velocity at various initial pressures. In the case of p 0 = 9.2 kPa, we have a combustion wave that propagates at 400 m/s and that fails upon the exit of the obstructed region (shown as a vertical line at x ≈ 320 mm). The case with p 0 = 38 kPa shows that within the porous space, the wave propagates at the quasi-detonation speed of about 1050 m/s. However, it fails upon the exit into the unobstructed space while the slightly larger p 0 = 42 kPa results in detonation transmission upon the exit even though in the obstructed region the speed is nearly the same. The range of initial pressures around 20-40 kPa is near critical, and one can observe different outcomes both within the pore space and on the exit from it.
The reaction wave propagation in the obstructed media such as the one at hand is subject to multiple phenomena that affect its propagation mechanism. One observes wave velocities that range from those typical to flames (on the order of centimeters to tens of meters per second), high-speed turbulent deflagrations (near but below the speed of sound in the fresh mixture), and quasi-detonations (supersonic, but with substantial velocity deficits compared to the ideal case). In experiments, one observes transitions between these regimes that are typically abrupt, but can also be continuous [6]. Understanding the interplay of the multitude of physical and chemical processes in the dynamics of the combustion wave in such obstructed media is complicated and requires further and more detailed experimental, computational, and theoretical analysis.

Discussion and Conclusions
The propagation of detonation in a layer comprised of a matrix of 8-mm steel cylinders was studied experimentally using high-speed self-luminosity observations and computationally using three-dimensional numerical simulations. It was shown that in the pore space, the lead shock wave velocity varied between 500 and 1600 m/s. The induction zone estimates based on the maximum velocity of the lead shock were found to be in agreement with observations that indicate that the reaction-zone length is on the order of the pore size. Numerical simulations of the same mixture under conditions matching the experiments reproduced this and a number of other features observed in the experiments. Furthermore, the simulations indicate that the heat transfer to the walls and the cylindrical obstacles may play an important role in the structure and the dynamics of the detonation waves.
Another important observation is that the simulations indicate the formation of extremely high-temperature localized regions that form due to the interactions of shock waves and solid boundaries and caused by them rapid chemical reactions. The presence of such high temperatures implies that one must be careful with the nature of the chemical reactions taking place to make sure that the relevant detailed reactions are properly included in the reaction mechanism. Treatment of turbulence in this compressible reacting flow is another issue that requires further consideration, especially concerning its influence on the chemical reactions. Another consequence, this time concerning numerical integration of the governing equations, is that the rapid heat release will impose strict limitations on the computational time step that may require using more efficient numerical methods that utilize implicit algorithms.
Overall, it is evident from the present analysis that despite the reasonably good predictions obtained with the present numerical simulations, in order to fully characterize the experimental observations and elucidate the physical nature of the wave-wave and wave-obstacle interactions in the system and to understand various regimes of reactionwave propagation in this complex reactive system, further research is required with more resolved experiments as well as improved simulations.

Conflicts of Interest:
The authors declare no conflict of interest.