Multiplicity for idealized rotational detonation waves

Detonation-based pressure gain combustion (PGC) has the potential to provide higher propulsive efficiency than conventional deflagration combustion.^1–3 Besides the widely investigated and developed pulse detonation engine (PDE), an alternative system referred to as the rotating detonation engine (RDE) has attracted worldwide attention in recent years.^2–5 RDEs present an elegant solution to the problems of complex injection/ignition sequencing and discontinuous thrust of PDEs by harnessing a self-regulating detonation wave propagating continuously in an annular combustor.^3–5 The axisymmetric geometry of most RDEs is both favorable for sustaining the detonation wave and more readily combined with designs typical of existing stationary and transportation gas turbines.^1,6 These benefits do have a cost, and the exact flow dynamics in an RDE is an area of continuing research need.³ Recent experimental, theoretical, and numerical studies on RDEs have been comprehensively reviewed in the recent past.^1–3,7

RDEs are associated with technical challenges arising from stability and controllability, which often result from the bifurcation of detonation waves within the combustor.^8–15 Such a transition can be affected by injector dynamics, nonuniform mixing, parasitic contact surface burning, shear layer instability, wall curvature, and combustor geometries.^1,3,4 A common finding in experimental studies^9,12,16–18 is that increasing the mass flow rate seems to increment wave numbers.² Understanding the dynamics of rotating detonation waves, particularly the bifurcation and multiplicity of waves, remains a challenge due to the complex internal flow field in RDEs and its close interactions with the injection process. Single/dual-wave hybrid mode can be observed within one single experimental test,¹⁸ where the spike of detonation pressure in the dual-wave mode is approximately half of that in the one-wave mode. The individual velocity of each wave drops by almost 15% after each bifurcation as observed by Bykovskii et al.⁹ A larger wave speed deficit of approximately 50%² is reported when there are eight¹³ or nine waves.¹⁰ 

Wave multiplicity in annular RDEs is also reported in numerical simulations.^19–30 Schwer and Kailasanath¹⁹ conducted two-dimensional (2D) and three-dimension (3D) simulations of premixed H₂/air³¹ to assess the effect of different injector configurations on the stability and performance. In their 2D simulations with a slot injection system, they found that the lower and higher mixture plenum pressure cases are more unstable than the case with an ideal injection model,³² owing to the additional unsteadiness in the fill region due to micro-injection jets. At higher pressures, shock waves propagating in the opposite direction of the main detonation wave can accelerate and form a new detonation wave. Subsequently, the two waves collide and extinguish each other repeatedly. In the 3D simulations from Yao et al.,²⁰ the injection intervals are adjusted along the headwall to account for nonuniform injection. The spontaneous formation of new detonation wave is successfully captured, which is induced by collisions of wave fronts. The gas intake configuration on the headwall, the stagnation of fresh mixture, and the number of initiation positions²⁰ are reported to affect the number of detonation wave fronts in the stable state. Sun et al.²² performed 2D simulations of RDEs operated with premixed H₂/air mixture through injection nozzles. They observed that the new detonation wave can be induced by the heat release at the interface between the fresh premixed gas layer and detonation products. A series of 2D simulations with a nine-species hydrogen chemical mechanism³³ of RDE with separate injection of fuel and oxidizer are conducted by Zhao et al.²⁹ Spontaneous wave formation is observed due to mutual enhancement of an explosive hot spot and a traveling shock wave. Influences of various fuel and oxidizer compositions on the combustion mode and detonation wave bifurcation are also studied by Zhao et al.²⁹ It is found that the velocity deficit is reduced with increasing mixedness of the reactants. Lietz et al.³⁰ simulated a 3D RDE with discrete fuel and oxidizer injector orifices. Effects of varying equivalence ratios, mass flow rates and the addition of a throat at the end of the combustor are studied. Multiple waves are present consistent with the experimental measurement.

As observed in these studies, wave multiplicity can be caused by many factors, both deterministic and stochastic, including the stagnation pressure of injection system,¹⁹ the characteristics of injectors and subsequent mixing,²⁹ the geometry of the combustor, the number of initiations,²⁰ deflagration in the fill region and along the contact line,²² and interactions between shock waves and local hot spots.²⁹ Efforts have been made to predict the number of waves in the literature,^1,4,14 although conclusive arguments are difficult to be drawn due to the convoluted factors mentioned above.

To understand impact of individual physical processes on wave multiplicity, in particular the effect of increasing mass flow rates and wave numbers, a numerical experiment with an ideal premixed injection system is designed in this study to avoid complexities arising from interactions between the mixing processes and the detonation. The initiation of detonation is carefully implemented in a mode-locked fashion to avoid spontaneous development of local hot spots and wave-wave interactions. The reacting Euler equations are solved to minimize parasitic deflagration along the contact surface. Perfectly premixed fuel/air streams are employed to avoid altering the mixture reactivity arising from inhomogeneous mixing. The objective is to identify whether wave bifurcation is a fundamental feature with increasing mass flow rate and wave number, and to identify the limit wave numbers and fill height under these idealized conditions.

The rest of the manuscript is organized as follows. The computational methods and numerical solvers are first introduced in Sec. II, followed by the configuration of the model rotating detonation combustor (RDC) in Sec. III. The local detonation wave structures, impact of mass flow rates and initiations of ignition zones are discussed in Sec. IV, where a theoretical relation on predicting wave numbers is discussed. Conclusions are finally drawn in Sec. V.

The reactive Euler equations for mass, momentum, species and total energy are solved in the Eulerian framework, as shown in Eqs. (1)–(4), respectively,

where ρ is the mass density, $u$ is the fluid velocity vector, p is the pressure and Y_k is the mass fraction for species k. The non-chemical total energy E is defined by $es+1/2|u|2$ using the specific sensible internal energy e_s. $q̇c$ is the heat release rate (HRR) that is defined by $∑k=1Nsω̇khf,k0$⁠, where $ω̇k$ is the reaction rate of the species k and $hf,k0$ is its enthalpy of formation.

The transport equations of all species are solved except for the diluent inert species (e.g., N₂ or Ar). The mass fraction of the diluent inert species is obtained by $1−∑k=1Ns−1Yk$ to enforce mass conservation where N_s indicates the number of species. A second-order central-upwind scheme³⁴ with vanLeer limiter is adopted to discretize convective terms, and an explicit second-order multi-stage Runge-Kutta scheme is employed in the temporal discretization.

An in-house chemistry solver has been developed to efficiently compute the chemical source terms using an analytical Jacobian technique.³⁵ The chemistry solver is coupled with the hydrodynamic solver using operator splitting. Fixed time steps are used to maintain the maximum Courant-Friedrichs-Lewy (CFL) number below 0.25. When running in parallel, a round-robin load-balancing technique is applied to redistribute the computational load of chemical source terms evenly to all processors. One-dimensional and two-dimensional verifications of the solvers are detailed in  Appendix A. In this study, an 11 species and 34 reaction hydrogen sub-model of the Foundational Fuel Chemistry Model Version 1.0 (FFCM-1)³⁶ is used for all simulations. FFCM-1 is designed as an advanced reaction model for combustion of small hydrocarbon fuels using up-to-date kinetic knowledge. It has been widely validated and verified and shipped with well-defined predictive uncertainties.³⁷ The hydrogen sub-model used here has been employed in a recent publication (Ref. 38) to analyze detonation cellular stabilities, where good agreement of cell width between the modeling results and experimental measurements is obtained.

An idealized injection surface is designed in this study. A stoichiometric mixture of hydrogen and air at a temperature T_in is injected into the annular combustor through the entire head-end area with a constant mass flow rate of $ṁin$⁠. Zero gradient pressure boundary condition is applied at the inlet, where the inlet patch maintains identical pressure as the first-level cell center. The density $ρin,i$ of the inflow gas mixture for patch i on the inlet surface is determined by the ideal gas law using local conditions. The inlet surface has a uniform mass flux and the inlet velocity on the local patch i can be evaluated by $Uin,i=ṁin/(ρin,iAin)$ where A_in is the area of the entire inlet surface.

This boundary condition is a further simplification of the widely adopted total pressure boundary condition from the literature,³² as compared in Table I. Numerically, the current constant mass flux boundary condition can be considered to be almost equivalent to the total pressure boundary condition, except that an isentropic condition is not enforced with the constant mass flux boundary. The parametric increase/decrease in mass flow rate in this study can be viewed as increasing/decreasing the total pressure. The proportionality between the mass flow rate and pre-shock pressure is shown in Sec. IV C.

The constant mass flux boundary condition is purposefully employed in this study to attempt to decouple internal flow features from the boundary conditions, because such interactions contribute significantly to the creation/destruction of waves. Following highly idealized studies of chaotic behavior of detonations in 1D,³⁹ and multiplicity of stable solutions in 2D,⁴⁰ the current study is intended to identify if the multiplicity of stable wave solutions is a fundamental feature of the annular detonation combustor.

An annular confinement with an inlet and outlet is constructed as a simplified representation of a rotating detonation combustor.⁴¹ A schematic of the combustor geometry is presented in Fig. 1, superimposed by the initial conditions of hydrogen mass fraction $YH2$ and pressure p. The outer-wall diameter d_o and channel width h are 100 and 8 mm, respectively. The channel length L_z is 100 mm. A similar configuration can be found in Ref. 41 where h = 10 mm. Three numerical pressure probes, P_i, P_m, and P_o, are arranged at three different radial locations, r_i, r_m, and r_o, respectively, close to the inlet surface (z = 2 mm) to record pressure traces and to monitor detonation wave velocity u_D. r_i is the inner radius, r_o is the outer radius, and $rm=(ri+ro)/2$ is highlighted by the white line on the initial pressure contour in Fig. 1(b). Note that the configuration in this study represents a thin-channel combustor (⁠$rm/h=5.75$⁠), and observations here are specific to this configuration. Effects of wider channel width and stronger curvature effects will be investigated in future studies.

The combustor is initialized with cold air at the inlet temperature T_in and atmospheric pressure. To ignite the mixture, a one-dimensional detonation wave structure with a height of 5 mm is superimposed onto a subsection of the domain near the inlet surface, denoted as “IZ” in Fig. 1, similar to a few previous studies.^6,42 A triangular fuel (H₂/air) zone is added in front of “IZ” to support continuous fuel intake to the initial detonation wave [see Fig. 1(a)]. High pressure and temperature in the post-shock region are filtered to stabilize the tail of the ignition zone “IZ.” Both inner and outer walls are treated as adiabatic and slip surfaces. A zero gradient boundary condition is specified for velocity, species, temperature, and pressure at the outlet. Influence of the outlet boundary condition and initial conditions of ignition zone height and pressure are examined in  Appendixes B and  C, respectively. Negligible differences are observed in the predicted pressure traces when an exhaust plenum is added to the domain or a different ignition zone height/pressure is employed.

Structured non-orthogonal hexahedral meshes are employed throughout the study. Uniform grids are specified in the radial (r) and azimuthal (θ) directions, while stretched grids are assigned in the streamwise direction (z) to ensure better resolution near the inlet region while maintaining a manageable mesh size. Table II summarizes all the meshes employed in this study. The coarest mesh M1 contains 1.15 million (M1) cells, and three refined meshes, M2 to M4, are employed to investigate the grid dependence of the results.

Table III summarizes all the parametric cases simulated in this study. The inlet static temperature T_in is 250 K, which is chosen based on a commonly employed total temperature of 300 K, for example, in Refs. 5, 32, and 43. In those simulations, the static temperature at the inlet is around 250 K. The baseline case, referred to as “$m/iz1/M3$⁠,” adopts a medium mass flow rate (“m,” 1.16 kg/s) and one ignition zone (“iz1”) with the mesh M3. Additional cases with inlet mass flow rates ranging from “$0.05m$” to “5m” are listed in the second block of Table III. To study the impact of initial conditions on wave numbers, in the third block, two “iz2” to twenty-four “iz24” ignition zones are investigated. Each simulation is first run for 5–7 cycles to ensure statistical stationarity after which data are collected for another 10–14 cycles that contain approximately 130 snapshots. Animations of temperature isocontours from cases “$0.05m/iz1/M3$⁠,” “$3m/iz1/M3$⁠,” “$m/iz24/M3$⁠,” and “$m/iz36/M3$” are shown in Fig. 2 (Multimedia view), Fig. 3 (Multimedia view), Fig. 4 (Multimedia view), and Fig. 5 (Multimedia view), respectively. Case “$m/iz36/M3$” is not listed in Table III because all detonation waves fail, and deflagration combustion is observed when 36 ignition zones are initiated.

The impact of the grid resolution on wave dynamics is first investigated using the baseline mass flow rate of 1.16 kg/s “m” and one ignition zone “iz1” with different grid resolutions (M1 to M4). The isocontours of temperature at successive time instances are shown in Fig. 6. All cases are initiated with identical scalar and velocity fields. For case “$m/iz1/M1$” where the coarsest mesh M1 is applied, a small reversed detonation wave arises from the first collision between the main detonation wave and the tail of the ignition zone, as shown at $t=$ 0.18 ms. Its height quickly grows before the reversed wave collides with the main detonation wave at $t=$ 0.21 ms. After multiple collisions of the rotating waves subsequently, two co-rotating waves with different fill heights are sustained in the combustor annulus. Refinement in the radial and azimuthal directions in case “$m/iz1/M2$” leads to a similar wave bifurcation process. When the mesh M1 is refined in the streamwise direction, as in case “$m/iz1/M3$⁠,” only a single rotating wave is sustained throughout the simulation. The one-wave mode remains the same by further refining the mesh M3 in the radial and azimuthal directions (mesh M4).

Figure 7 presents the cell-centered values of temperature and pressure within the first layer of computational cells along the outer wall in cases “$m/iz1/M1$⁠,” “$m/iz1/M2$⁠,” “$m/iz1/M3$⁠,” and “$m/iz1/M4$” at $t=$ 0.12 ms before the first collision. An apparent difference in temperature can be found between M1 (or M2) and M3 (or M4) due to the differences in the cell height along the streamwise direction. Behind the detonation wave at $θ=100°$⁠, the predicted temperature in M1 and M2 is higher than that from M3 and M4 until the temperature drops to the unburnt value at $θ=$ 140°. The maximum difference in temperature is approximately 830 K. On the contrary, the mesh resolution in the radial and azimuthal directions has negligible effects on the cell temperature when we compare the results from M1 and M2 or from M3 and M4.

The over-prediction of temperature is introduced by numerically mixing hot products with fresh mixture when the computational cell is large in the streamwise direction. In meshes M1 and M2, such a numerical hot spot is formed along the contact surface after the first collision between the main detonation and its tail, which subsequently grows and transitions to form another detonation wave propagating in the opposite direction of the main detonation wave. Subsequently, the two waves collide with one another, and wave elimination and regeneration occur frequently. Eventually, two stable waves are sustained in the chamber. However, this numerically (as opposed to physically existing) hot spot can be avoided by further refining the grid size in the streamwise direction (e.g., in meshes M3 and M4) because relatively less amount of the product gas is numerically mixed with the fresh mixture near the inlet, which indicates that the bifurcation is an artifact resulting from the coarse mesh. For the 3D RDC simulation, it is critical to have sufficient streamwise resolution to capture the complex flow-chemistry coupling near the inlet, such that hot spots and wave bifurcation are not numerically introduced.

The radial and azimuthal resolutions have less impact on the temperature near the inlet because mixing with hot products mainly occurs in the streamwise direction. Comparing case “$m/iz1/M3$” and case “$m/iz1/M4$⁠,” the relative difference of detonation wave speed u_D and fill height l_cr are 1.5% and 2%, respectively. However, finer radial and azimuthal resolutions are necessary to capture the details of the cellular structure, as reported in Ref. 44. Figure 8 compares the pressure contours from the first-layer cells obtained from case “$m/iz1/M3$” and case “$m/iz1/M4$” at t = 1.98 ms. In case “$m/iz1/M4$⁠,” the detonation front moves further ahead of that on the outer wall due to a stronger expansion and compression of detonation waves captured by the finer resolution. Behind the detonation wave front, a stronger oblique shock is captured owing to the compression of the detonation wave near the outer wall and subsequently, shock trains that are established by the reflected oblique shock waves bouncing between the inner and outer walls⁶ are better captured in Fig. 8(b).

Based on the resolution studies in this section, the “M3” mesh is chosen as the baseline and employed in all subsequent simulations. It should be noted that the grid resolution used in this study is of similar order of magnitude to prior simulations^6,45 with similar geometries and fuels, which is not intended for capturing cellular structures or reaction zone details for hydrogen/air mixtures. However, as demonstrated in  Appendix A, the wave speed can be well captured with the current resolution, which is the target of the subsequent simulations.

The structure of the detonation wave is briefly discussed in this section. Figure 9(a) shows instantaneous streamlines colored by the frozen Mach number for case “$m/iz1/M3$⁠.” The majority of the outlet is supersonic for the baseline case, as seen from the color contour of the Mach number. The outlet is entirely subsonic when the mass flow rates are below “$0.5m$⁠,” while it is entirely supersonic when the mass flow rates exceed “1m.”

For better visualization, the 3D computational domain is unwrapped into a two-dimensional (2D) plane in Figs. 9(b) and 9(c). The temperature iso-contours on the outer wall are employed to show the main detonation features in the combustor. Behind the detonation wave, the fresh mixture penetrates into the annular chamber, expands almost linearly, and forms a triangular unburned fill region with a defection angle of contact surface α. The fill height l_cr is defined as the maximum height of the fill region, which is quantified using 10% of the maximum heat release rate in this study. In front of the detonation wave, the contact surface separates the products from the fresh mixture of different velocities. Minimal contact line burning is observed in the simulation, as evidenced by the heat release rate contour in Fig. 9(c). Above the detonation wave, a slip line originates at the intersection of the detonation wave, and an oblique shock wave exists between the fresh detonation products and products from the previous revolution in Fig. 9(b). The deflection angle of the slip line is referred to as δ. The back dashed line CD indicates the top side of the detonation wave.

The fill height l_cr has been a quantity of interest in many experimental^8,9,12,14,17 and numerical studies.⁵ It has been hypothesized in Ref. 14 that given the chamber width h and cell width λ, multiple waves are spawned when the upstream fresh mixture exceeds a critical fill height. l_cr depends on the stagnation pressure, temperature, mixture composition, its uniformity, and the channel geometry.⁴ In this study, we are particularly interested in whether a critical value of l_cr exists to lead to wave bifurcation, which is discussed next in Sec. IV C.

Simulations with a range of inlet mass flow rates are conducted using the mesh M3, as listed in the second block of Table III. All simulations observe one stable detonating wave. The fill height l_cr is plotted in Fig. 10(a) as a function of inlet mass flow rates. The fill height increases with increasing mass flow rates and plateaus around 23∼24 mm when the mass flow rate exceeds 1.16 kg/s (case “$m/z1/M3$”). The error bars indicate the fluctuations of l_cr. The fluctuations are largest at the largest mass flow rates due to a self-adjustment mechanism²¹ and appear to attain smallest values when mass flow rates of 1.16 kg/s (case “$m/z1/M3$”) and 1.74 kg/s (case “$1.5m/z1/M3$”) are specified. At lower mass flow rates (cases “$0.05m/z1/M3$⁠,” “$0.1m/z1/M3$⁠,” and “$0.2m/z1/M3$”), the fluctuations are partly caused by the wrinkled contact surface due to the Kelvin-Helmholtz instability as shown in Fig. 2 (Multimedia view). Figure 10(b) compares the detonation wave speed u_D and the Chapman–Jouguet (CJ) detonation velocity u_cj as a function of inflow mass flow rates $ṁin$⁠. u_D is computed by dividing $2πrm$ by time intervals between peak pressure that is recorded by the pressure probes. u_cj is calculated based on the pre-shock pressure p_pre and temperature T_pre extracted from the locations immediately in front of the shock near the inlet patch (first-layer cell). Both velocities drop as $ṁin$ becomes lower, and the velocity deficit (1−$uD/ucj$⁠) is maximum at lowest $ṁin$⁠. The reduction of u_cj is due to a decrease in pre-shock pressure p_pre. The larger disparity between u_cj and u_D at lower $ṁin$ is further discussed in Sec. IV E.

Figure 11 shows mean pressure profiles as a function of the azimuthal angle, normalized by the pre-shock pressure for each case. The profiles almost collapse into one line when the mass flow rate exceeds 0.58 kg/s (case “$0.5m/iz1/M3$”), which corresponds to outlet conditions becoming at least partially supersonic. A modeled pressure trace proposed by Shepard and Kasahara⁴⁶ is also plotted in Fig. 11, which is described by

where $ξ=x/lcr$ is the distance x from the shock front normalized by the fill height l_cr. The parameter β is a fitting parameter, which is equivalent to the inverse of the integrated normalized pressure over ξ and $β=0.96$⁠. From fitting the pressure trace, the fill height l_cr is determined to be 23 mm, which is consistent with the measured values using the heat release rate.

For mass flow rates larger than 0.58 kg/s (case “$0.5m/iz1/M3$”), $p¯/p¯pre$ decays toward unity right before reaching the shock front at $θ=180°$⁠. This indicates that the pressure expansion process takes the entire circumference to occur, which is the maximum distance x that can be supported by the current geometry. In contrast, the expansion processes are complete well ahead of the shock front for the smaller mass flow rate cases. These observations suggest that the averaged circumference introduces a geometrical constraint to the maximum fill height that can be supported by a given combustor. Pressure is adjusted to this constraint by maintaining similar $p¯/p¯pre$ profiles when the mass flow rate exceeds “0.5m” although the pre-shock pressure increases with increasing mass flow rates, as shown in Fig. 12(a). Note that no wave bifurcation is observed in these simulations even when the pre-shock pressure becomes extreme (e.g., 5–15 bars), which shows that wave bifurcation does not necessarily result from increasing stagnation pressure or inlet mass flow rates; when no local hot spots or parasitic burning are present, as carefully designed and implemented in this study, one-wave mode can be sustained across a wide range of inlet conditions.

The smallest mass flow rate in this study is 0.06 kg/s (case “$0.05m/iz1/M3$”), below which detonation cannot be sustained. At this condition, the fill height l_cr is approximately 8 mm [see Fig. 10(a)], and the wave speed is approximately 85% of the corresponding CJ velocity u_cj [see Fig. 10(b)]. Figure 12(b) shows the pressure expansion distance dis_exp, which is defined as the distance where post-shock pressure falls to 0.5% of the pre-shock pressure. The hydrodynamic thickness h_t, which is measured between the shock front and the Ma = 0.995 location from the Zeldovich-von Neumann-Doering (ZND) solution, is also plotted. With decreasing mass flow rates, the pre-shock pressure drops, and the hydrodynamic thickness h_t increases. Interestingly, the expansion distance at the smallest mass flow rate corresponds to the hydrodynamic thickness at corresponding conditions. Further reduction of the mass flow rate would induce even smaller pre-shock pressure and larger hydrodynamic thickness h_t, while the fill height l_cr and pressure expansion distance dis_exp continue to drop. Consequently, the expansion distance dis_exp falls below h_t, which can be interpreted as a violation of the necessary choking condition to sustain detonation waves. The hydrodynamic thickness h_t has also been employed to conduct scaling studies of the reaction layer thickness in a series of 2D simulations of detonation wave bounded by inert gases.^47,48 They show that the hydrodynamic thickness appears to allow a better scaling analysis of the results compared to the choice of cell sizes and can be a useful characteristic length scale for the detonation.^47,48

Wave bifurcation in the presence of multiple waves is investigated in this section through initiating multiple ignition regions that are evenly distributed, as shown in Fig. 13. All thermo-physical parameters are maintained identical for each ignition location, which corresponds to the mode-locked condition that is described in Ref. 49.

As an example, Fig. 13 shows the evolution of detonation waves in Case “$m/iz3/M3$⁠,” which is initialized with three evenly distributed ignition zones (see t = 0 ms). Three detonation fronts are formed (see $t=0.05∼0.21$ ms) and are sustained in the annular chamber (see t = 1.9 ms). The number of stable detonation waves is identical with the number of ignition regions. An obvious reduction in the fill height l_cr is observed when the number of waves increases. Similar observations are made for all other cases up to 24 ignition zones as shown in Fig. 4 (Multimedia view).

Figure 14(a) shows the fill height as a function of wave numbers. The fill height is inversely proportional to the number of waves. Simulations with 36 [see Fig. 5 (Multimedia view)] or 48 initial waves both fail to detonate, where deflagration combustion is advected downstream. Up to eight or nine waves have been reported in the literature,^10,13 where wave-wave interactions or deflagration to detonation transition (DDT) from local hot spots often contribute to wave bifurcation or coalition. The current study indicates that the presence of multiple existing waves in the system is not a sufficient condition for wave-wave interactions and/or wave bifurcation when waves are perfectly mode-locked and no additional stochastic hot spots are available. Figure 14(b) shows that the wave speed u_D almost drops linearly with increasing number of waves, while the CJ velocity u_cj is almost constant as there is little change in the pre-shock pressure p_pre as shown in Fig. 15(a). Decreasing wave speed has been reported when new waves are generated in experiments.^10,13,49 Intrinsically, variations of the velocity deficit $(ucj−uD)/ucj$ as a function of wave number are caused by the lack of confinement on the top side [back dash line CD in Fig. 9(b)] of detonation wave. More analysis is shown in Sec. IV E to explain the observation.

The averaged pre-shock pressures p_pre and pressure expansion distance dis_exp are plotted in Fig. 15. It can be seen that the pre-shock pressure slightly increases with increasing number of waves. For the 24 wave case, the pressure expansion distance is approximately twice the hydrodynamic thickness, and the distance between two neighboring waves is approximately three times the hydrodynamic thickness. It is expected that the pressure expansion distance is further reduced to be close to the hydrodynamic thickness with 36 waves, which eventually leads to failure.

In addition to viscous and heat loss, the lack of confinement on the top side of the traveling waves is a third mechanism leading to a velocity deficit in rotating detonation combustors.⁴⁹ The loss from pressure expansion in RDCs can be explained by a pioneer study by Fay,⁵⁰ where the effect of boundary layer on gaseous detonations in solids tubes was studied. Suppose the shock-front area of each streamtube is A₁ and the CJ-plane (or sonic-plane) area is A₂, the average fractional change in the area of each stream tube ξ is given by $(A2/A1−1)$⁠, which can be employed to represent the expansion loss. By considering a 2D channel of width b bounded by inert gases on one side normal to b, Dabora et al.⁵¹ showed that $ξ=(x¯· tan δ)/b$ where δ is the deflection angle of the interface between the inert layer and hot product and $x¯$ is the reaction length separating the CJ plane from the shock front. For a given fresh mixture, a larger velocity deficit is observed in the experiments when b is thinner and ξ is larger. Similar trends are also captured in the recent 2D simulations.^47,48 Following their analysis and replacing $x¯$ and b by h_t and l_cr, respectively, the fractional area increase ξ for the current RDC is defined as

where the deflection of contact surface $(ht· tan α)$ is subtracted from the deflection of the flow. Figure 16 shows the variations of velocity deficit $(ucj−uD)/ucj$ with fractional area increase $[ht·(tan δ−tan α)]/lcr$ from the cases with different inlet mass flow rates $ṁin$ and wave numbers W. For both cases, the velocity deficit is positively correlated with the fractional area. Figure 17 compares the fractions of $ht/lcr, (tan δ−tan α)$ and $[ht·(tan δ−tan α)]/lcr$ from all cases. It can be seen that the change in $[ht·(tan δ−tan α)]/lcr$ is mainly contributed by the variation of $ht/lcr$ as the change in $(tan δ−tan α)$ is relatively insignificant compared to $ht/lcr$⁠.

When the mass flow rate is decreased, h_t becomes longer with decreasing pre-shock pressure, l_cr becomes shorter, and hence $ht/lcr$ gets larger. When wave number is increased, h_t is almost fixed [see Fig. 15(b)] because p_pre barely changes as shown in Fig. 15(a). l_cr decreases with wave number W [see Fig. 14(a)]; therefore, $ht/lcr$ increases with W. In this study, we observed that detonation quenches when the pressure expansion distance is comparable to the hydrodynamic thickness. From a physical point of view, when the velocity deficit increases to certain point, the shock strength drops so much that the post-shock temperature and pressure cannot ignite the mixture within a reasonable amount of time, which eventually leads to failure.

Results in Sec. IV C demonstrate that the large mass flow rate is not a sufficient condition for wave bifurcation. Furthermore, results in Sec. IV D show that presence of multiple waves is also not a sufficient condition for wave bifurcation. In addition, the number of waves cannot be uniquely determined simply based on the inlet boundary condition and the combustor geometry because the number of initiation locations, or stochastic hot spots, can determine the number of waves present in a combustor. In either case, the fill height and wave number are closely related to the pressure expansion process along the azimuthal direction of the combustor, which seems to play a significant role here.

Attempt has been made in the literature to theoretically to predict the stability of detonation front and detonation wave number W.^1,4 For example, Eq. (7) from Ref. 1 predicts the wave number W by comparing the timescale of one revolution t_r to the timescale of consuming the fresh mixture ahead of the detonation wave t_mf,

Here, $tr=πdm/uD$⁠, and the detonation velocity u_D is also evaluated at the averaged diameter location. The critical volume V_cr is calculated through

where V_cr indicates the consumed volume that is swept by the detonation wave within one revolution. The volumetric supply rate of the fresh mixture can be computed as

The detonation wave number W equal to 1, 2, 3, or n suggests that there are 1, 2, 3, or n detonation heads propagating in one direction in the combustor. Re-arranging Eq. (7), we obtain Eq. (10), which essentially states that the number of detonation waves is determined by the ratio of the inlet mass flow rate $ṁin$ and the mass consumption rate by each detonation wave $ρmix(uDhlcr)$⁠,

Further re-arranging Eq. (10) and substituting $ṁin$ by $ρmixuinAin$⁠, l_cr can be expressed by

where $dm=(do+di)/2$⁠. Equation (11) shows that the fill height l_cr scales with the average circumference $πdm$ divided by the number of waves W. Consistent with the observations in Secs. IV C and IV D, $πdm/W$ is the maximum azimuthal expansion distance for each wave. The non-dimensionalized l_cr is determined by the velocity ratio u_in/u_D. The proportionality between l_cr and d_m has been discussed in previous 2D simulations.⁵² This equation quantitatively explains the change of the fill height due to varied u_in and d_m predicted by a fast 2D model developed by Sousa et al.⁵³ based on mass conservation.

Figure 18 shows the left-hand side vs the right-hand side of Eq. (11) using all the simulations data in this study. A linear relation is observed, which is a natural result of mass conservation in all the simulations. However, Eqs. (10) or (11) cannot be expected to “predict” the number of waves in a specific RDC, given specified operating conditions. To “predict” wave numbers, it is important to accurately account for all the deterministic and stochastic events that can lead to wave bifurcation, e.g., local hot spots and DDT. Therefore, it would be more meaningful to study the limit conditions and critical physico-chemical conditions that lead to those, rather than predicting the exact number of waves in a system. Results from the current study attempt to provide the upper and lower limits of l_cr when increasing mass flow rates and increasing wave initiations, under the assumptions of reactive Euler governing equations, perfectly mode-locked initial conditions, perfectly premixed inlet boundary conditions, and constant inlet mass flux. In addition, the various lengths and velocities are measured based on certain arbitrarily chosen thresholds. Therefore, the criteria of the lower and upper limit conditions in this study should be interpreted as an order of magnitude indication. Future studies will focus on exploring these limit conditions while considering viscous effects and inhomogeneous mixtures where more stringent resolution requirements apply.^54,55

Numerical simulations of a three-dimensional model RDC are conducted using a recently developed compressible reacting flow solver to study mechanisms of wave bifurcation. The solver has been verified using 1D hydrogen/oxygen/argon detonation and a 2D hydrogen/air configuration, where the detonation wave speed is well captured by a large range of grid resolutions. Subsequently, the wave dynamics in a model RDC with premixed hydrogen/air is studied. An ideal injection system with a fixed mass flow rate boundary condition is proposed to avoid complexities that arise from complex injector dynamics. Influences of grid resolutions, inflow mass flow rate, and initiations of ignition spots are investigated and main findings are summarized below.

The first-level cell height near the inlet surface should be sufficiently refined because it directly controls the mixing between hot products and the fresh mixture. Large cell height increases the cell-centered temperature near the inlet, which can serve as an igniter to numerically introduce multiple waves. One detonation wave is maintained throughout the simulation once the mesh is sufficiently refined in the streamwise direction.

As inflow mass flow rate increases, the peak pressure of the detonation wave becomes proportionally higher. The dynamics of fill height seems to be controlled by the coupling between the internal pressure field and the inlet boundary condition. The circumference of the combustor constrains the maximum expansion distance for the post-shock pressure, and a maximum value of fill height is achieved when the downstream outlet boundary condition becomes supersonic. Increasing mass flow rates beyond this point does not lead to wave bifurcation or further increase in the fill height although the fluctuation of the fill height seems to grow with extreme mass flow rates. The minimal fill height observed in these numerical experiments corresponds to the smallest mass flow rate at which the pressure expansion distance is approximately the same as the hydrodynamic thickness that is defined based on ZND solutions.

The number of waves is not uniquely determined by the inlet boundary conditions; multiple initial ignition regions can lead to multiple stably rotating waves. Consequently, the fill height cannot be solely determined based on the inlet condition as it is closely related to the number of existing waves in the combustor. The existence of multiple waves is also not a sufficient condition for wave-wave interactions or wave bifurcation; they can stably co-exist in a mode-locked manner when no hot spots or disturbance exist. The maximum number of waves, hence the minimum fill height, is again constrained by the pressure expansion distance, which is measured to be approximately two hydrodynamic thickness for 24 waves.

Either lower mass flow rates or larger wave numbers can increase the velocity deficit, which results from lack of confinement on the top side of the traveling waves in rotating detonation combustors. The concept of a fractional area increase in the CJ (or sonic) plane compared with the area of the shock front is employed to explain the loss. For both cases, the ratio of the hydrodynamic distance and fill height is the critical parameter that determines the level of velocity deficit.

A mass-conservation based interpretation of Wolański's formula^1,4 is derived, where the fill height scales with the circumference divided by the wave number, and can be determined by the ratio of inlet velocity and wave speed. Data from the simulations show consistency with this relation, simply indicating that all the fuel has been converted to products within the combustor. Synthesizing all the results in this study and those observed in the literature, it is shown that it is unlikely that the number of waves, hence the fill height, can be “predicted” based on the inlet condition and the geometry alone; other factors, such as stochastic local hot spots, can lead to wave bifurcation. Future studies will focus on understanding the limit conditions for the fill height when more realistic physics is considered, including the detailed detonation structure, the viscous effects and heat loss through walls.

The authors acknowledge the support of AFOSR under Grant No. FA9550-18-1-0173 under the supervision of project monitor Dr. Chiping Li.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Two mixtures are employed in the verification tests, including a stable mixture (H₂/O₂/Ar) and an unstable mixture (H₂/air). The 11 species, 34 reaction hydrogen sub-model of the Foundational Fuel Chemistry Model Version 1.0 (FFCM-1)³⁶ is employed in the verification tests. The solver is verified using one-dimensional stable detonation, 1D pulsating detonation and two-dimensional (2D) cellular detonations. The influence of the grid resolution is subsequently investigated using the three configurations.

A stoichiometric mixture of H₂/O₂ diluted with 70% argon (2H₂:O₂:7Ar) is employed with an unburnt pre-shock temperature T_pre of 300 K and pre-shock pressure p_pre of 0.2 atm. The detonation is initialized through a ZND solution obtained from the Shock and Detonation (SD) Toolbox.^56,57 The total domain length is 0.8 m and the length of the fresh mixture in front of the leading shock is 0.6 m. The detonation wave is expected to propagate with a Chapman-Jouguet (CJ) velocity u_cj= 1647.1 m/s into the quiescent unburnt gas. Zero gradient boundary condition is specified at the outlet boundary. Uniform grids of $Δx=10$⁠, 20, 50, 100, 200, and $500 μ$ m are employed, which correspond to 50, 25, 10, 5, 2.5, and 1 grids per induction length (0.5 mm), respectively.

Figure 19(a) shows the history of the peak pressure p_max obtained using different grid resolutions. The black dashed line indicates the reference p_max of the ZND structure evaluated in SD Toolbox.^56,57 Stable and accurate p_max can be maintained when the fine meshes of 10 and 20 $μm$ are applied. Once the resolution is reduced to 50 $μm$⁠, slight oscillations of the peak pressure occur. The p_max starts to drop as the grid size is further increased, and the oscillation of the peak pressure slightly grows. Figure 19(b) compares the history of shock locations identified by the peak pressure p_max obtained from different grid resolutions. The black dashed line indicates the reference shock location that shifts with the CJ speed u_cj at 1647.1 m/s. The trajectory of p_max from all simulations agrees with the reference line very well, which implies a good prediction of the detonation wave speed from all resolutions. Using linear regression of the shock location, the detonation speeds from all resolutions match u_cj within 0.1%, which is shown in Fig. 19(c) (red axis). Figure 19(c) provides the relative error of the mean peak pressure (black axis) and its corresponding standard deviation denoted by the error bars. The relative error of the peak pressure is bounded between 4.7% and 6.4%.

The detonation wave structures (at t = 0.3 ms) are further compared with the ZND solution in Fig. 20. As the grid becomes coarser, the discontinuity induced by the leading shock is smoother because the cell-centered value is obtained by volumetric averaging in finite-volume methods. The post-detonation states are well-predicted by all resolutions, which is reasonable since they are determined by thermodynamics only.

The 1D validation tests show that the detonation speed can be reasonably captured even when the grid resolution is very low. The details of the ZND structure will be compromised with very coarse meshes. This observation supports the usage of coarser meshes for the 3D simulation in Sec. IV because the detonation speed is the primary concern there.

Validation is further conducted for the less stable hydrogen/air mixture in this section. A one-dimensional detonation tube (0.6 m) is considered, which is filled with stoichiometric hydrogen/air at pre-shock temperature T_pre of 298 K and pressure p_pre of 42 kPa. A high-pressure p_driver, high-temperature T_driver driver gas consisting of hot product is used in a small region (10 mm) near the closed end of the tube. The driver pressure ratio $pdriver/ppre$ is 120 and the driver gas temperature is 2000 K. Zero gradient boundary condition is specified at the outlet boundary. The Chapman–Jouguet (CJ) velocity u_cj for the mixture is 1955.9 m/s. Uniform grids of $Δx=2.5$⁠, 5, 10, 20, 50, 100, 200, and $500 μ$m are employed, which correspond to 200, 100, 50, 25, 10, 5, 2.5, and 1 grids per induction length of 0.5 mm, respectively.

Figure 21 shows the peak pressure p_max as a function of time for grid resolutions of $Δx=2.5$⁠, 5, and $10 μ$m. A high-frequency mode of oscillation is observed, which is different from the stable p_max of hydrogen/oxygen/argon shown in Fig. 19(a). The frequencies are 1.10, 1.05, and 1.05 MHz for $Δx=2.5$⁠, 5, and $10 μ$m, respectively. The corresponding detonation velocity u_D and degree of overdrive $f=(uD/ucj)2$ are approximately 2107 m/s and 1.16, respectively, for the three grid resolutions. The computed frequency and degree of overdrive f are in close agreement with Lehr's experimental measurement,⁵⁸ where a oscillation frequency of 1.04 MHz is obtained for the same mixture at p_pre = 0.421 atm and f is approximately 1.10 estimated in Ref. 59. The current results are also similar to the simulation results reported by Yungster and Radhakrishnan,⁶⁰ where the oscillation frequency is 1.06 MHz for an overdrive of f = 1.09 at a pressure of p_pre = 42 kPa and T_pre = 298 K. Additionally, Romick et al.⁵⁹ reported an oscillation frequency of 0.97 MHz for an overdrive of f = 1.10 at a pressure of p_pre = 0.421 atm and T_pre = 293.15 K.

Once the resolution is reduced to 20 $μm$⁠, lower-frequency oscillations that have higher amplitude (not shown) appear. Further reducing the grid resolution damps out the pulsation. Although the frequency of pulsation is sensitive to grid resolutions, the detonation speed u_D from all resolutions remain comparable.

To further assess the performance of the solver on reproducing the intrinsic instability of 1D detonations at CJ state, a 1D detonation near CJ state (degree of overdrive f = 1.03) is obtained by carefully tuning the driver pressure ratio $pdriver/ppre$ to 80 and the driver gas temperature to 2500 K. The uniform grid size is $5 μ$ m as this resolution has been found to capture the regular high-frequency oscillations of p_max when f = 1.16. Figure 22 shows the peak pressure p_max as a function of time. The detonation wave first propagates in the high-frequency and low-amplitude oscillation mode, and subsequently the instability transitions to a low frequency (0.1 MHz) and high amplitude mode. The same trend is also observed in Ref. 60.

A 2D detonation with a stoichiometic mixture of hydrogen/air propagating in a narrow channel is employed for multi-dimensional verification of detonation structures under different grid resolutions. The quiescent unburnt gas is at a temperature T_pre of 300 K and a pressure p_pre of 1 atm. The channel width is 2 mm. The same conditions were adopted by Sato et al.⁶¹ for solver validation and the domain length is extended to 150 mm in this study. The detonation is initialized through a ZND solution sinusoidally perturbed in the transverse direction.⁶² Top and bottom walls are treated as adiabatic and slip. Zero gradient boundary condition is specified on both ends of the 2D channel. Three grid resolutions of uniform $Δx=40$⁠, 20, and $10 μ$m are employed, which correspond to 5, 10, and 20 grids per induction length (0.2 mm), respectively.

The numerical soot foils are employed to visualize the cellular patterns produced by the detonation, which were generated by storing the maximum pressure in each grid cell during the simulation. Figure 23 shows the smoke foils in the sub-domain within the first 100 mm. As the grid resolution is increased, more detailed cellular structures near triple points are captured. The maximum pressure increases with finer resolutions, which is consistent with the observation in previous 1D tests in Sec. A 1. All predicted cells have a uniform width of 1 mm, which agrees with the simulation results of Sato et al.⁶¹ 

As the detonation wave propagates further downstream beyond 100 mm as shown in Fig. 24, the cellular structure obtained using $Δx=10 μ$m becomes irregular both in shape and in size. At t = 7.2 × 10^–5s, the front of smoke foils are highly distorted by triple points and transverse waves. However, regular cellular patterns are maintained in the other two coarser-mesh cases and the smoke foil fronts tend to be flat. For the unstable mixture of H₂/air, a tendency to develop more irregular cell sizes has also been observed by Taylor et al.⁶³ as the numerical resolution increases.

To evaluate averaged shock locations and detonation wave speeds, the first grid cell with pressure larger than p_pre along the streamwise direction is first recorded at various wall-normal locations, and the averaged shock location is then determined by averaging the shock front location in the wall-normal direction. The history of shock locations from all resolutions overlaps with one another (not shown). The overall detonation wave speeds evaluated from the slope of the transient shock locations are 1990.5 m/s (+0.8 $%ucj$⁠), 1992.3 m/s (+0.9 $%ucj$⁠), and 1996.5 m/s (+1.1 $%ucj$⁠), respectively, for resolutions $Δx=$ 40, 20, and 10 μm.

Although different in structures, we found that the locations of the detonation front predicted by the three resolutions are almost identical, which indicates that the detonation wave speed can be captured by all three resolutions in this 2D test. The observations in 1D and 2D tests justify the usage of relatively coarse resolution if only the detonation speed is of primary interest.

An exhaust plenum with a diameter of 10d_o (1000 mm) and a length of 10 L_z (1000 mm) is added to the original domain to move the exit boundary farther away from the combustor^41,45,64 as shown in Fig. 25. The head-end and the sides of the exhaust plenum are all adiabatic slip walls.⁶⁴ A far-field boundary condition with a back pressure of 1 atm is applied to the outlet of the exhaust plenum. Figure 26 compares the pressure history of case “$m/iz1/M3$” with and without exhaust plenum at the probe location P_m. No obvious differences are observed throughout the simulations. Therefore, we consider the current outlet boundary condition without the exhaust plenum adequate for the purpose of this study.

The impact of initial conditions is investigated through case “$m/iz1/M3$” by changing the ignition zone height h_iz and initial pressure, and restarting the simulation from the solution of case “$0.5m/iz1/M3$⁠.” Figure 27 shows the pressure traces of the three cases with the ignition zone height h_iz of 1, 5, and 25 mm at the probe P_m location and no significant changes are observed. The strength of the oblique shock waves following the detonation wave is not affected by h_iz. Figure 28 shows the pressure traces of the three cases with the initial pressure of 0.5, 1, and 2 atm at the probe P_m location. After the transient stage (seven revolutions), the pressure immediately in front of the wave is independent of the initial pressure.