Multidimensional Numerical Modeling of Combustion Dynamics in a Non-Premixed Rotating Detonation Engine With Adaptive Mesh Refinement

Conventional gas turbine engines widely deployed in both aerospace propulsion and stationary power generation applications operate on the Brayton cycle, wherein combustion occurs in the constant-pressure deflagrative mode. In recent years, detonative combustion has been explored as an alternative approach to significantly enhance the efficiency of modern gas turbine engines [1]. A detonation cycle leads to nearly constant-volume (as opposed to constant-pressure) heat release, thereby resulting in a corresponding pressure gain. The thermal efficiency from a detonation cycle could be 20% higher than that of conventional isobaric combustion [2]. Moreover, extremely high thrust can be generated from a relatively small combustor by leveraging the more intense and rapid energy release from detonation. The prospect of these substantial benefits has spurred a number of research efforts worldwide over the last few decades pertaining to the development of advanced detonation engines [1].

One such novel combustor concept that has received wide attention over the last few years is known as a rotating detonation engine (RDE). The operation of an RDE can be summarized as follows. In an RDE, fuel and air are injected at one end of an annular combustion chamber, and a continuous detonation, once initiated, propagates in the azimuthal direction within the annulus near the inlet, while the detonation products are expanded and purged from the chamber further axially downstream [1,3,4]. RDEs can provide a nearly steady source of thrust at very high frequencies without having to initiate a detonation repeatedly. Moreover, their annular geometry and compact design with no moving parts make them very adaptable to existing aircraft or rocket engine hardware.

Many experimental studies have been performed to gain an understanding of RDE operation and performance and visualize the associated combustion phenomena, employing techniques such as static pressure and thrust measurements [5], OH* chemiluminescence imaging [6,7], tunable diode laser absorption spectroscopy [8], and more recently, mid-infrared imaging [9] and megahertz-rate imaging [10]. However, rotating detonation engines still remain challenging to study using only experiments, owing to the harsh environment in the combustor chamber comprised supersonic flows, high-frequency oscillations, high pressures, and temperatures. In this context, multidimensional modeling and computational fluid dynamics (CFD) can complement engine experiments and serve as a valuable tool to provide more insights into the RDE combustion dynamics and enable optimization of critical design parameters that influence RDE operation and performance.

Several numerical investigations have been conducted to understand the fundamental physics of rotating detonation engines. Most of the earliest CFD studies [11–15] were focused on the overall description of the flow-field within the RDE combustion chamber and typically employed “unrolled” two-dimensional (2D) geometries with premixed fuel/oxidizer mixture injected into the annulus, although a few three-dimensional (3D) simulation studies are also available in the literature [16,17]. Schwer and Kailasanath [13,14] identified seven major regimes within an RDE. In order to capture detonation, they utilized the induction parameter model (IPM) [18] wherein an induction variable (τ) was convected, with its source-term dependent on the induction time at local composition and thermodynamic conditions. This induction time was computed from either experimental correlations or lookup tables generated a priori using detailed kinetic models. In regions where τ > 1, the reactants were converted to products, and the heat was released over a characteristic reaction time scale. To ensure the stability of the numerical scheme, the heat release was spread out over several time-steps. IPM was extensively used for computational analysis of the effects of operating conditions, combustor geometric parameters as well as fuel chemistry on RDE performance [14,15,19–21]. Paxson [22] and Paxson and Naples [23] studied RDE combustion using a high-resolution CFD algorithm that integrated the quasi-2D, single-species, reactive Euler equations with source terms in the detonation frame of reference. The reaction rate was considered to be proportional to the product of a non-dimensional rate constant, density, and reactant mass fraction. Moreover, a user-defined threshold temperature was prescribed below which no reactions were allowed to occur. Recently, 2D CFD simulations of an idealized RDE geometry were performed using this code [24], and it was shown that turbulence could cause reaction zone enlargement resulting in the reduction of detonation wave speed.

One of the limitations of these previous studies was the assumption of reactants being perfectly premixed. Hence, the effect of fuel-air mixing on RDE combustion could not be properly incorporated, which is a critical governing factor for RDE wave dynamics [25]. In addition, the combustion models employed in these works were not capable of properly capturing deflagrative burning and spatially varying local mixture composition which is encountered in more realistic non-premixed configurations having separate fuel and air injection schemes. A few attempts were made to account for low-pressure deflagrative combustion, for example, by either utilizing a non-detonation model (based on global Arrhenius reaction rate) below a certain threshold value of local pressure [26] or artificially lowering the reaction rate constant in specified regions of the computational domain [22]. The effect of fuel-air mixing was accounted for in some premixed RDE simulations [23] by specifying a certain number of axial computational cell rows near the inlet where no reaction would occur even if the temperature exceeds the threshold value, thereby mimicking the delay due to finite mixing time before combustion starts. However, these approaches were based on ad hoc user-defined criteria which would need to be varied case-by-case depending on the operating condition and configuration.

More recent numerical studies have explored RDE combustion phenomena in full-scale 3D configurations incorporating the fuel and oxidizer injection plenums. Strakey et al. [27] simulated the non-premixed RDE configuration at US Air Force Research Laboratory (AFRL) [6] using an unsteady Reynolds-averaged Navier–Stokes (URANS) turbulence model and a tuned 1-step kinetic mechanism for hydrogen/air mixtures. Although the experimental trends [6] were well captured, there were large errors (∼15–20%) in some of the model predictions. Cocks et al. [28] performed 3D improved delayed detached eddy simulations (IDDES) [29] of the AFRL RDE using a detailed hydrogen/air mechanism [30] and finite-rate chemistry-based combustion model wherein the chemical source terms were calculated in each computational cell assuming simple laminar chemistry. Although a good agreement was found with experimental data for mean axial static pressures, the simulations over-predicted wave speeds by 10–30%. Moreover, these simulations were computationally expensive with a total cell count of approximately 93 million. In another study, Yellapantula et al. [31] simulated the same RDE using URANS simulations with a cell count of 24 million per simulation. High-fidelity large eddy simulations (LESs) of RDEs have also been performed by various groups [32,33], but remain computationally expensive.

It has been well established that a multitude of factors affects the flow field and operation of RDEs, such as fuel and oxidizer compositions, global equivalence ratio, mass flow rates, stagnation and back pressures, injector geometry, and detonation channel geometry [6,15]. Optimizing these multiple design parameters simultaneously to maximize RDE performance requires CFD tools that are: (1) sufficiently predictive in capturing the essential combustion dynamics and parametric trends in realistic 3D configurations and (2) computationally affordable to explore a large design space. The primary objective of the present work is to develop a well-validated CFD model for combustion in non-premixed RDEs that meets the above two requirements. A unique feature of the modeling strategy employed in this study, which hasn’t been explored in previous works [27,28,31–33], is the application of adaptive mesh refinement (AMR) to achieve a trade-off between simulation accuracy and computational expense.

The remainder of the paper is organized as follows. The next section provides a brief description of the RDE configuration and experimental conditions studied. Then, the details of the numerical model are presented in the following section. Subsequently, results from the numerical simulations are discussed and compared with experiments. Finally, the paper concludes with a summary of the main findings and future work.

In this study, numerical simulations were performed for the representative non-premixed RDE at US AFRL [6] as shown in Fig. 1. In the experiments, the air was injected from a plenum into an annular detonation channel through a circumferential slot with a height of 1.78 mm and a diameter of 123 mm. Fuel (hydrogen in this case) was injected from a separate plenum through 120 holes, each with a diameter of 0.89 mm. The fuel injection holes were evenly spaced on a circle with a diameter of 134 mm located inside the inner edge of the detonation channel. The width of the channel was 7.6 mm, with inner and outer diameters of 138.7 mm and 153.9 mm, respectively. The height of the annular channel was 101.6 mm. The fuel and air mass flow rates were metered upstream of the respective plenums using sonic nozzles. The RDE operation sequence involved establishing the air and fuel flow followed by detonation initiation. Detonation in the channel was initiated using a small tube (6.35 mm diameter, 63.5 mm long). The small tube was filled with a stoichiometric hydrogen-oxygen mixture, which was spark ignited. Deflagration-to-detonation transition (DDT) occurred in the small tube. The detonation then entered the annular channel and initiated the RDE operation.

Time-averaged static pressure measurements were acquired in the air plenum, fuel plenum, and between 6.35–70 mm downstream of the fuel injection surface at spatial increments of 6.35 mm. Experimental data were reported after steady-state conditions had been achieved in the plenums, which was defined to occur when the detonation limit-cycle had been reached and the phase-averaged static pressure in the detonation channel remained independent of time [6].

In this work, three operating conditions corresponding to different air/fuel mass flow rates but the same global equivalence ratio of around unity were chosen for the numerical study. Cases 1 and 2 exhibited one detonation wave each and case 3 exhibited two detonation waves in experiments. These cases are listed in Table 1. The ambient pressure and temperature for both cases were 101.1 kPa and 295 K, respectively.

A commercial CFD code, converge (version 2.4) [34], was employed to perform the numerical study. The computational domain for the RDE geometry is shown in Fig. 2. A detonation tube of 6.35 mm diameter was added to the combustor in order to replicate the detonation initiation in experiments (described in the previous section) as closely as possible. On comparing Figs. 1 and 2, it can be observed that the fuel and air plenums are not included in the simulation setup.

converge uses a modified, cut-cell Cartesian method for grid generation directly during runtime. In addition, it has the capability to include fixed embedding of cells i.e., increasing the grid resolution with respect to the base grid size a priori and AMR to locally refine regions with the steepest gradients in the domain [34]. In the present work, a base mesh size of 2 mm was used. Three levels of fixed boundary embedding (cell size of 0.25 mm) were employed near the fuel/air inlets and fuel ports, while two levels of fixed boundary embedding (cell size of 0.5 mm) were specified near all the remaining boundaries of the combustor. Two levels of fixed region embedding (cell size of 0.5 mm) were prescribed within the detonation tube. In addition, four levels of spherical embedding (cell size of 0.125 mm) with a diameter of 6.35 mm were imposed at the far end of the detonation tube to resolve the spark ignition kernel. Two levels of AMR were employed based on velocity, temperature, and density sub-grid values of 0.5 m/s, 1.5 K, and 0.01 kg/m³, respectively. This resulted in the minimum grid size of 0.5 mm and peak cell count of around 4.5 million per simulation, which is much lower than previous studies [28,31]. Grid convergence tests were performed by (a) reducing the base grid size up to 1 mm and (b) increasing the AMR level to 3. Negligible variation in the detonation wave frequency (within 1%) with grid refinement was observed.

Governing transport equations for mass, momentum, species, and total energy was solved in the URANS framework. A variety of turbulence models (standard k–ɛ, realizable k–ɛ, k–ω) [35] were tested with no significant variation in detonation wave frequency. Eventually, the realizable k–ɛ model with wall functions was used for the final simulation study to provide closure for the unsteady turbulent flow equations. The model proposed by O’Rourke and Amsden [36] was employed to account for wall heat transfer. An ideal gas equation of state was used as the constitutive relation. As the CFD solver is based on collocating all computed values at the center of the cell, the Rhie-Chow algorithm [37] was utilized to prevent checkerboarding. A second-order-accurate spatial discretization scheme (with step flux limiter) was used for the governing equations with a fully-implicit first-order-accurate time integration scheme. The transport equations were solved using the Pressure Implicit with Splitting of Operators (PISO) method of Issa [38]. A geometric multigrid solver was used for the pressure solution. In addition, a variable time-stepping algorithm was employed in the current study with a minimum time step of 1 ns. The time step during the simulation was automatically calculated by the solver during each computational cycle based on the maximum allowed Courant-Friedrichs-Levy (CFL) numbers of 0.5, 2.0, and 1.0, corresponding to convection, diffusion, and speed of sound, respectively.

A finite-rate chemistry model [39] was employed to solve for the reacting flow in the combustor, thereby not accounting for turbulence-chemistry interaction [40–43]. A kinetic mechanism containing ten species and 21 reactions [30] was employed to describe detailed hydrogen-air chemistry. This mechanism has been shown to capture the detonation wave speed accurately under a wide range of conditions in one-dimensional detonation tube simulations [28]. Constant mass flowrates from experiments were imposed together with the reflective condition at the fuel and air inlets, while a transonic boundary condition was specified at the outlet. The law-of-the-wall boundary condition for velocity and isothermal boundary conditions with a temperature of 300 K was applied to the solid surfaces. In light of the absence of any wall temperature measurements from experiments, the prescribed value in simulations was considered to be a reasonable approximation considering that the experimental run times were quite short (1 s). Lastly, several monitor points were set up at the fuel/air inlets as well as multiple locations along the axis of the detonation channel to record local static pressures, which were then used to compute the detonation wave frequency/velocity and time-averaged static pressure values for comparison against experimental data.

Each simulation was initialized with both the RDE combustor and detonation tube at constant temperature and pressure of 295 K and 101.1 kPa, respectively. The combustor was filled with air whereas the detonation tube was filled with a stoichiometric hydrogen-oxygen mixture. In addition, the interface between the detonation channel and detonation tube was kept closed. Then, a non-reacting flow simulation was performed until it reached a steady-state. Thereafter, this steady flow solution was used as the initial condition for the reacting flow simulation. The detonation tube was the spark ignited at one end, leading to DDT within the tube. The interface between the channel and the tube was then opened for a short duration to let the detonation wave enter the annulus and initiate detonation in the chamber. Afterward, the interface was closed for the remainder of the simulation. It took approximately 5–6 ms of physical time for the simulations to reach the detonation limit cycle. Completion of each simulation took ∼60–72 h on 144 processors.

Numerical simulations were performed for the three operating conditions listed in Table 1. In this section, the simulation results are discussed and compared with available experimental data.

Figure 3 shows the instantaneous temperature and density contours on the mid-plane of the detonation channel (r = 73.15 cm) for case 1. These contour plots depict key features of the flow that also have been recently observed in experiments by Rankin et al. [9] under the same conditions: (a) detonation wavefront propagating from right to left, (b) reactant fill region ahead of the detonation wave where fuel and air mix, (c) products from the current detonation cycle that expand behind the detonation, (d) oblique shock from the downstream edge of the detonation wavefront to the combustor outlet, (e) products from the previous detonation cycle downstream of the reactant fill region, (f) shear layer between products from the current and previous detonation cycles, and (g) reflected shock arising from the interaction of the oblique shock with channel exit plane. However, it must be noted that the present simulations do not resolve certain flow characteristics, such as Kelvin–Helmholtz instabilities at the reactant–product interface and vortical structures in the shear layer, which were observed in previous studies with detailed numerical simulations [14,44]. This is primarily attributed to the URANS approach and relatively coarse grid sizes employed in this study.

The detonation wave heights for cases 1 and 2 are qualitatively compared in Fig. 4, which shows pressure contours on the channel mid-plane under both operating conditions. It is evident that case 1 exhibits a taller detonation wave than case 2. This is due to the fact that higher mass flowrates in case 1 result in a taller reactant fill region ahead of the detonation wave. On the other hand, Fig. 5 shows the temperature contour plots on the channel mid-plane for case 3, which has the highest air/fuel mass inflow rates. It can be seen that two detonation waves are spawned under these conditions. A quantitative comparison between the numerical results and experimental data with respect to detonation wave dynamics is provided in Table 2. A clear trend of the detonation wave height initially increasing (from case 2 to case 1) and then decreasing (for case 3) with mass flowrate can be observed, similar to experiments [6,9]. Note that the wave frequency of case 3 is nearly double that of case 1 due to the presence of two waves. The predictions of detonation wave frequency/height, reactant fill region height, and oblique shock angle are evidently in good agreement with experiments. Overall, the CFD model is able to capture the essential characteristics of detonation wave structure and speed in the RDE reasonably well.

Figure 6 depicts the initial transient dynamics leading to the formation of two detonation waves in case 3. Initially, only one wave is observed. However, subsequent mixing of fresh reactants and products leads to a high-temperature zone in the wake region behind the original detonation wave. This zone later results in the spawning of a second wave which starts to propagate in the left direction due to the presence of a fresh mixture. After a while, the two waves ultimately get stabilized leading to the quasi-steady-state behavior as shown in Fig. 5, with identical wave heights.

A comparison of time-averaged pressures in the fuel/air plenums and detonation channel was also carried out between experiments and simulations as shown in Table 3. The third column in the table corresponds to the location in the RDE annulus at an axial distance of 2.54 cm from the fuel injection plane. Interestingly, simulations overpredict the plenum pressures while showing good agreement with experiments as far as channel pressure is concerned. This is likely due to the fact that the plenums were not included in the computational domain; the reported values are, therefore, calculated at the locations of the fuel and air inlets instead of the actual air/fuel plenums shown in Fig. 1. Moreover, reflective boundary conditions were employed at the fuel and air inlets in simulations. This caused the pressure waves to bounce around after upstream propagation from the passing detonation. As shown in Fig. 7, these pressure oscillations never really dissipate, thereby increasing the local pressure. The experimental measurements, on the other hand, were acquired comparatively much further upstream and were likely less affected by pressure waves. On the other hand, Fig. 8 shows the evolution of static pressure with time within the annulus at the axial distance of 2.54 cm from the fuel injection plane for case 1. Clearly, pressure recovery after the passing of the detonation wave occurs rapidly at this location unlike the fuel/air inlets.

Figure 9 shows a comparison of the experimental mean static pressure distribution to the corresponding numerical predictions at different axial locations on the outer wall of the detonation channel. Good agreement was obtained between the two for both cases, with simulations correctly predicting the trend of highest pressure at the most upstream location.

Apart from model validation, qualitative and quantitative analysis of detonative versus deflagrative burning was also carried out for the three cases. Figure 10 shows the azimuthal distribution of static pressure for cases 1 and 2 on a horizontal cut plane passing through the respective detonation waves. The corresponding pressure distribution for case 3 was similar to that of case 1 (except for the presence of two waves instead of one), hence is not shown here for the sake of brevity. It can be observed that for case 1, the detonation wave nearly spans the full channel width. However, for case 2, the high-pressure wave only seems to span roughly half of the annulus width. The corresponding temperature contour plots for the two cases on the same cut plane are shown in Fig. 11. Clearly, near the inner diameter of the annulus, high-temperature regions can be observed for case 2 ahead of the detonation wave, indicating deflagrative burning. This explains the apparent propagation of detonation waves only through the outer half of the channel width.

Table 4 shows the fractions of overall heat release occurring under lean and rich conditions for the three cases. It can be seen that the majority of the heat release occurs in fuel-rich regions for cases 1 and 3, whereas heat release corresponding to case 2 is biased towards lean conditions. Table 5 shows the heat release as a function of pressure for each case. Here, 5 bar has been chosen as the threshold beyond which heat release is considered to occur at high pressures. Evidently, heat release at high pressures is significantly greater for cases 1 and 3, indicating higher effective pressure gain than case 2. Interestingly, detonative heat release shows a non-monotonic trend with air/fuel mass flowrate, in particular case 1 exhibiting greater high-pressure heat release compared to case 1. This can be attributed to the presence of two waves, which leads to less time available for recovery of low-temperature fresh reactants and fuel-air mixing in the fill region before being consumed by the detonation wave. On the other hand, case 2 exhibits higher combustion efficiency (96.6%) than case 1 (90.4%) and case 3 (88.3%), where combustion efficiency is defined as the fraction of fuel consumed before exiting the RDE annulus. This trend in combustion efficiency is directly linked to the relative lean bias of the three cases as shown in Table 4.

In summary, a computationally efficient and predictive CFD model for realistic full-scale non-premixed RDEs was developed in the present work. The URANS approach along with AMR was employed. Such a simulation tool can be readily employed for rapid design space exploration and optimization of RDE performance. Through extensive verification and validation against experimental data, it was shown that the model could capture the essential detonation wave dynamics with good accuracy. In addition, heat release analysis based on simulation results provided a good understanding of the detonative burning occurring under different operating conditions. Recently, the CFD modeling framework with AMR presented in this paper has also been demonstrated for high-fidelity LES of hydrogen-air [45] and ethylene-air [46] non-premixed RDEs.

A novel CFD modeling approach was developed to simulate complex combustion phenomena in full-scale non-premixed RDEs. To achieve a good trade-off between simulation accuracy and computational expense, the URANS turbulence modeling framework, and AMR were employed. Simulations were performed for three operating conditions corresponding to the same global equivalence ratio but different mass flowrates. The simulations were able to capture key features of the reacting flow within the RDE reasonably well, including multiple waves. Model validation was carried out against available experimental data on a number of waves, wave frequency/height, reactant fill height, oblique shock angle, fuel/air plenum pressure, and axial pressure distribution in the channel. The CFD model was shown to be able to capture the sensitivity of all these wave characteristics to the operating conditions with good accuracy, both qualitatively and quantitatively. In addition, a heat release analysis was performed to compare the combustion behaviors observed under the simulated conditions. The main advantage of the present CFD modeling framework is higher computational efficiency compared to previous studies, and hence, the potential capability to perform rapid design space exploration and/or performance optimization studies for RDEs. Future work will involve the incorporation of full fuel and air plenums into the computational domain, and modeling of hydrocarbon combustion and back-pressurized RDEs with the present numerical approach.

The authors would like to thank Dr. Peter Cocks from the United Technologies Research Center (UTRC) for many helpful discussions during the course of this study. In addition, thanks to the visiting Ph.D. student, Krishna Kalvakala, from the University of Illinois at Chicago (UIC), for assistance in postprocessing some of the simulation data. Lastly, the authors would like to acknowledge the computing core hours on the Bebop cluster provided by the Laboratory Computing Resource Center (LCRC) at Argonne National Laboratory for this research. The submitted paper has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (Argonne). Argonne, a U.S. Department of Energy Office of Science Laboratory, is operated under Contract No. DEAC02-06CH11357. The research work was funded by Argonne’s Laboratory Directed Research and Development (LDRD) Swift Project No. 2018-108-N0.

There are no conflicts of interest.