Numerical Investigation of Contact Burning in an Air-Breathing Continuous Rotating Detonation Engine

Three-dimensional (3D) numerical simulations of a continuous rotating detonation engine are carried out with an unsteady Reynolds-averaged Navier-Stokes solver. The second-order upwind advection upstream splitting method and second-order Runge-Kutta method are used to discretize space and time terms, and detailed 9-species 19-step hydrogen-oxygen reactions are applied in this study. Nonpremixed rotating detonation is successfully realized numerically, and the characteristics of the detonation wave are revealed. The expanding angle of the combustor has a great impact on the shape of the detonation wave but has little influence on the propagation velocity. The evolution of combustion on the contact region is analyzed in detail; a more accurate schematic of non-premixed air-breathing rotating detonation engines is given in this paper. A rough analysis of the heat performance of the contact region shows that the heat release of the contact region is approximately one-third of the total heat release and the configurations of the combustors do not affect the proportion.


Introduction
Continuous rotating detonation engines (CRDEs) have attracted increasing attention due to their high efficiency and simple structure [1,2]. The CRDE usually adopts an annular combustion chamber. There are one or more detonation waves propagating continuously in the circumferential direction in the combustion chamber. The highpressure products expand rapidly and are spewed out at a high speed from the end of the combustion chamber [3][4][5][6][7][8][9].
However, in a real CRDE, complex secondary combustion (deflagration) effects will occur due to the nonuniformity of mixing of the fuel and air. Even in the ideal condition of two-dimensional premixing, there is still 23.6% deflagration in the detonation combustion chamber [30].
Numerous experimental studies have also shown that complex secondary combustion is present in the detonation combustion chamber. There are many manifestations of secondary combustion in the CRDE. Bykovskii et al. [7,8] found that mixed gas in front of the detonation wave is continuously ignited by the high-temperature product. In experiments by Frolov et al. [31], it was determined that there are some reaction fronts (one, two, or more) behind the detonation wave, which should be the premixed gas being ignited by the high-temperature products. The propagation speed of these deflagration points is still significantly small, so there is not enough time to induce a new detonation wave. Rankin et al. [12] also experimentally revealed the highbrightness hydroxyl (OH) region in the mixture before the detonation wave, referred to as the autoignition region. Chacon et al. [32] observed secondary combustion in the OH planar laser-induced fluorescence (OH-PLIF) test, including contact burning, autoignition near the injectors, and continuous deflagration due to inadequate mixing after the detonation wave. Peng et al. [13] studied the competitive relationship between detonation and deflagration through experiments. The deflagration before the continuous rotating detonation wave was parasitic combustion, which destroyed the accumulation of the combustible mixture. An increase in the intensity and longitudinal range of parasitic combustion reduces the operating range, propagation frequency, and stability of the CRDE. There are three types of secondary combustion: contact combustion at the front of the wave (termed parasitic combustion), continuous deflagration combustion behind the wave (termed commensal combustion) [32], and self-ignition.
There are many factors that affect secondary combustion including but not limited to the mixing process, injection conditions, and configuration of the combustion chamber. Gaillard et al. [33] studied various cold flow injection schemes and analyzed their mixing efficiency by large eddy simulations (LESs). Kawasaki et al. [34] and Jodele et al. [35] experimentally investigated the combustion efficiency of a CRDE with different inner diameters. A change in the inner diameter significantly affected the mixing environment. It was stated that if the inner diameter was less than a critical value, the velocity of the detonation wave would be greatly reduced. We can speculate that a small inner diameter makes the space near the nozzle too large, which causes the hot products to not be quickly pushed downstream, resulting in more serious secondary combustion, thereby destroying the conditions for the propagation of the detonation. Sun et al. [36] studied the effect of the nozzle width on the detonation process by three-dimensional (3D) numerical simulation. The interaction between the jet process and the propagation of the rotating detonation wave was very strong when the jet throat width was too large, and the detonation wave was quenched due to insufficient injection of reactants.
Many OH-PLIF tests [12,35,37] and simulations have shown that there is severe contact burning in the CRDE. Obviously, contact combustion is an important combustion component in the detonation combustion chamber. However, there is still a lack of sufficient research on this topic. Up till now, it is only known that 23.6% of heat release comes from contact burning in the detonation combustion chamber (all the deflagration comes from contact burning in this case) [30]. There is a lack of quantitative research on contact combustion in three-dimensional non-premixed detonation combustors because it is difficult to strictly distinguish the region of the contact combustion in threedimensional non-premixed detonation combustors whether experimental or numerical research.
The CRDE has been investigated by Euler equations [18,38], RANS equations [39][40][41][42], and large eddy simulations (LES) [43,44]. Edward et al. [45] compared the CRDE simulation results of Euler equations, unstable RANS equations, and LES equations. The study found that all methods can capture unstable detonation waves, while there are obvious differences in the FFT analysis of highly frequency bands, which is obvious because the unstable RANS cannot capture the fine detonation wave structures. It is found that in all methods with different viscosity models, the dominant frequencies and magnitudes are slightly different, which indicated that the heat releases, driving force of detonation waves, are only slightly different in different viscosity models. So unstable RANS is sufficient as the heat release of detonation wave is mainly concerned in this paper.  International Journal of Aerospace Engineering In this paper, 3D numerical simulations of a continuous rotating detonation engine are carried out with a RANS solver. Graphics processing units (GPUs) are applied to accelerate the calculations. The main objective of this paper is to develop further insight into the formation process of contact combustion and the proportion of heat release from contact combustion to the total heat release. Two configurations are considered to analyze the effects of the combustion chamber configuration on contact burning. This paper is organized as follows. Simulation details, such as numerical methods, physical models, inflow boundary conditions, free stream flow parameters, and code validation, are introduced in Section 2. In Section 3, the formation process and heat release of the contact combustion are analyzed. Moreover, the effect of combustor configurations on contact combustion is also analyzed by comparison. Lastly, the conclusions are given in Section 4.   3 International Journal of Aerospace Engineering isolator and the combustor is straight with a radius of 40 mm. The outer radius of the isolator is 43.31 mm, and the outer radius of the exhaust is 54 mm. The combustor of case 1 has an expansion angle of 3.4833°, while that of case 2 has an expansion angle of 6.6544°. The fuel (hydrogen) is ejected through 90 uniformly distributed injectors with a diameter of 0.6 mm in the inner cylinder.

Physical Model and Numerical Method
The simulations are run within an in-house 3D unsteady compressible RANS equation solver based on a GPU. More details about the code can be found in Ref. [46]. This code has been applied and verified in a variety of supersonic flow and combustion investigations [47][48][49]. The second-order upwind advection upstream splitting method (AUSM +-UP) scheme is used for inviscid fluxes, and the second-order central difference scheme is used for the viscous fluxes. The time advancement is performed by a four-step secondorder Runge-Kutta method, and the shear stress transport (SST) k − w model is used for turbulence modeling. A time-splitting method is used to uncouple the flow and reaction as follows.
where L f and L c denote the solver for flow and reaction, respectively, and n is the parameter that adjusts the numerical rigidity. In our previous study [49,50], it was found that n = 5 which makes the time step of the chemical reaction one order of magnitude smaller than that of flow and is enough to avoid numerical rigidity. A turbulence-chemistry interaction model with finiterate chemistry and compressibility correction finite-rate model is used for chemical reactions [51]. The model is modified from the Partially Stirred Reactor (PaSR) model [52] and was used in RANS simulations successfully. The reaction source term can be given as where _ ω is the source term of the laminar flow finite-rate model and γ is the fine-scale structure volume fraction, generally adopting the formula proposed by Here, the mixing time-scale τ m = ffiffiffiffiffiffiffiffiffi τ k τ Δ p , where Kolmogorov time-scale τ K = ffiffiffiffiffiffi ffi ν/ε p , ν is the viscosity coefficient, and ε is the dissipation rate and can be estimated by ε = 0:09kw in RANS study with the SST k − w turbulent model.

International Journal of Aerospace Engineering
The grid time-scale τ Δ can be estimated by τ Δ = Δ/u ′ , where velocity fluctuation u ′ = ffiffiffiffiffiffiffiffi ffi 2k/3 p and the grid characteristic length Δ = V 1/3 , where V is the local volume of the grid. And the chemical reaction time-scale τ ch can be deduced from the approximation of a laminar premixed flame, that is, τ ch ≈ ν/S 2 u , where S u is the propagation velocity. The flame propagation velocity is based on a one-dimensional conservation law and can be estimated by is the generation rate of the fuel, ρ u is the density of the unburned gas, the thermal diffusivity α = K T /ρ u c p , K T is the thermal conductivity, and c p is the specific heat capacity at constant pressure.
The ignition delay is significantly important for nonpremixed simulation since deflagration plays an important role in the CRDE. Therefore, a detailed classic 9-species 19-step reaction mechanism [53] is adopted in this study. This mechanism has been verified in detail in our previous supersonic combustion simulations [47,48], and it can accurately reflect the ignition delay of H 2 gas.
The flow rate is applied for both H 2 and air inlet boundaries, which is commonly used in the simulation calculation of the CRDE [36,43,54]. The flow rate of air is 600 g/s, and that of H 2 is 16.9 g/s. The global equivalence ratio is set as 0.98. The total temperature of H 2 is 300 K, and that of the air is 606 K (corresponding to the flight condition Ma = 3 at an altitude of 10 km). The nonslip and adiabatic walls are considered in all the wall boundaries. Supersonic extrapolation is applied for the outflow boundary condition.  Figure 2 shows the OH contours of three different grid sizes, and the corresponding detonation propagation velocities and detonation wave heights are listed in Table 1. From Figure 3, it can be seen that the three different grid sizes result in similar flow field structures, and every case has only one detonation wave in the combustion chamber. The propagation velocities and the detonation wave heights are very close. The similar results may be due to that the three meshes are all refined with grids of 0.1 mm near the nozzle, so the main jet mixing process of the three meshes is the same. Considering the amount of calculation, the grid accuracy of 0.2 mm is finer enough for the present study.
For the mixing efficiency, it can be obtained from the following equation: where where φ is the mixing efficiency, α is the mass fraction of the jet component, α mixed is the mixed fraction of the jet participating in the reaction, _ m fuel,mixed is the mass flow rate of the jet participating in the reaction, _ m fuel,total is the total jet flow rate, ρ and u are the local density and velocity, respectively, A is the cross-sectional area of the measuring position of the flow direction of the flow field, and α stoic is the mass fraction of the jet equivalent ratio. For hydrogen and air, α stoic = 0:0283 [55].  It can be seen from Figure 4 that the mixing distance between hydrogen and air is far less than the detonation wave height. Therefore, the fuel in front of the detonation wave can be consumed by the detonation wave, which shows that the combustion chamber configuration studied in this paper is reasonable and can provide the necessary mixture for the detonation wave under the working conditions.
The reliability of the simulation results is affected by the total number of computational time steps. Too many time steps will make the integration error too large [56]. In general, the range of the total error S max allowed for calculation is usually 1-5%. The integration error of the calculation domain can be obtained by S err ≈ ∑ 3 i=1 ðΔL i /L i Þ k+1 , where Δ L i and L i are the average size of the grid and the size of the computational domain, respectively, and k is the order of the numerical scheme. Therefore, the maximum number of computational steps allowed is n max = ðS max /S err Þ 2 . According to the calculation of the parameters in this paper, the integration error is S err = 1:14 × 10 −6 , and if S max = 1%, the maximum allowable number of calculation steps is n max = 7:7 × 10 7 . One cycle of the detonation wave takes approximately 17,000 time steps in this study, and the whole calculation period is far less than the maximum number of allowed calculation steps. The accumulation of the stochastic errors in numerical simulations of the unsteady flow was evaluated based on the error estimation method of Smirnov et al. [56,57]; thus, the simulation results have high reliability.

Results and Analysis
Two configurations with different expansion angles are investigated. The code is running on 6 V100 GPUs of a GPU workstation with a fixed step time of 5 * 9 −10 s. The detonation wave takes about 14 hours to propagate a cycle of the combustors. The contour results shown below are the tenth cycle of the detonation wave after the detonation reaches a steady state, while the statistical results are the average values of 10-15 cycles. Figure 3 shows the 3D instantaneous isosurfaces of pressure colored by temperature and isosurfaces of H 2 . There is a single detonation wave in both configurations. To analyze the propagation of DWs, the pressure history (15 mm downstream of the nozzle near the inner wall) is recorded in Figure 5. The main parameters are shown in Table 2, which are the average values of ten detonation wave cycles.

Characteristics of the Detonation Wave.
By comparing the parameters in Table 1, it seems that there is a large difference in the detonation height. This may be attributed to the flow path of case 1 which is narrower than that of case 2, which makes the flow velocity of case 1 larger than that of case 2. However, the difference in the wave height seems to have little impact on the propagation frequency and detonation velocity. The Chapman-Jouguet (CJ) speed is 1691.5 m/s calculated by using the NASA Chemical Equilibrium with Applications (CEA) code. The velocity loss is less than 10% for both cases and even less for case 2. Figure 6 shows a comparison of the current calculation results with the LES results and experimental results in the literature [32,58]. The current calculation results are in agreement with the results of LESs and experiments, and the main difference is that the results obtained by the unsteady RANS calculation in this paper are more regular than those of the LESs and experiments. Parasitic combustion in front of the wave is demonstrated, including two different combustion regions, namely, the first contact burning (CB1) and the second contact burning (CB2). The region between these two contact burning regions (called the buffer region) seems to be a fuel.

Evolution of Contact Burning.
The unwrapped equivalence ratio contour of case 1 shown in Figure 7 further supports that the buffer region is fuel-rich. Another difference between the present work and LESs and experiments is that there are regions of combustion below CB2 that are caused by autoignition. The total temperature of the incoming air is only 606 K in this study, which will not seriously cause autoignition. The last difference is the temperature of the buffer region, which will be explained below.
The fresh mixed gas affected by the combustion products is shown in Figure 8. In terms of scale, the contact burning region accounts for a large proportion of the fresh mixture.
Chacon et al. conducted a detailed experimental study on the parasitic combustion of CRDE [37], as represented in Figure 9. The meaning of the label can be seen below. The formation of CB1 (label 4) and CB2 (label 7) was described in detail. CB1 was ignited by the products of the previous detonation cycle, and the second contact combustion was caused by an additional ignition mechanism.
As seen from the schematic, CB1 and CB2 seem to be two parallel lines starting from two different positions, likely  7 International Journal of Aerospace Engineering because they are based on different causes described by Chacon et al. [37].
The slices in Figure 10 show the evolution of the OH contours in a circle time of the detonation wave. The newly generated OH region can be observed in Figure 10(a), after which the detonation wave has traveled some distance, which is because the pressure of the detonation wave may be higher than the injection pressure, and it may take some time for hydrogen to jet out. It is obvious that the inner part of the contact surface is the fuel-rich region, and the outer   International Journal of Aerospace Engineering part is the oxygen-rich area; the equivalence ratio demonstrates this. Then, as shown in Figures 10(b)-10(h), the OH region is constantly stretched and deformed streamwisely. The continuous production of OH means that continuous burning happens, which formed the CB1 and the CB2, and there is a fuel-rich region between them.
When the next detonation wave passes by, the mixing of hydrogen wrapped in the OH region and the oxygen outside the OH region will be strengthened due to the shock wave, resulting in a continuous OH region between the slip line and the oblique shock wave, which can also be found in Figures 6(b)-6(c). 9 International Journal of Aerospace Engineering Therefore, based on the above analysis, the flow field described in Figure 9 can be updated, as shown in Figure 11. It is mainly for non-premixed air-breathing CRDEs and should be more practical for the actual physical process of CRDEs.
Labels 1 through 5 in Figure 11 are consistent with Figure 9. They are the detonation wave, backflow, inflow mixture, first contact region, and hot product in order. These are the common characteristics of all CRDEs. The main differences are described as follows: (a) CB1 (label 4) and CB2 (label 7) start very close to each other and become further apart as time passes, so the buffer region between them (label 6) looks similar to a triangle, while CB1 and CB2 in Figure 7 are approximately parallel (b) Figure 11 highlights the continuous burning area (label 11) between the oblique shock wave and before the slip line. OH and heat release are continuously produced here. This is also shown in Figure 5(f) and other experimental studies [43,59]. Based on the present work, the continuous burning area is closely related to the buffer region (label 6). The fuel in the buffer region is intensively mixed with the air outside this region due to the detonation wave. Then, continuous deflagration occurs, and the OH region evolves (c) They have different temperatures of the buffer region and autoignition of the fresh reactant mixture of the fill region (label 8). The two differences are mainly due to the calculated condition, which has been explained above and will not be detailed here. It seems that the buffer region could be larger if the pressure of the jets is much larger than that in the present work 3.3. Heat Release of the Contact Burning. The percentage of deflagration combustion in the combustion chamber is very important but difficult to count precisely. The only method shown in the literature appears to analyze the flow particle path used in 2D premixed conditions [30], but this is unrealistic for 3D analysis. A careful observation of Figure 5(b) shows that the hightemperature region is concentrated in three regions: detonation and deflagration following the detonation wave; the contact region; and the continuous burning area between the slip line and oblique shock wave. Therefore, it can be assumed that heat release in the detonation combustion chamber is also concentrated in these three regions.
The enthalpy change is used to represent the heat release in this paper. Isosurfaces of different enthalpy changes are shown in Figure 12. The enthalpy change is very concentrated.
To quantify the degree of the concentration of heat release, the local enthalpy change ratio Γ is defined to represent the sum of the enthalpy changes that are greater than a certain heat release value.
where H is the enthalpy change. H +∞ and H -∞ are the maximum values and minimum values of the enthalpy change, respectively. Figure 13 shows that from 1E10 to 1E12, as the enthalpy change increases, the change of Γ is not very significant. There is a significant decline after 1E13. This can also show that the heat release in the combustion chamber is intensive. Since it is difficult to strictly distinguish the heat releases of detonation combustion and deflagration combustion in the CRDE, if we simply separate these three regions by their position, the error caused will be very small and will not affect our conclusions. Therefore, the heat release in the 10 International Journal of Aerospace Engineering combustion chamber is divided into three parts by position in this paper, as shown in Figure 14. Figure 14 shows the division of the burning regions. The combustion in region A includes almost all detonation combustion heat release and a large amount of deflagration combustion heat release. It is assumed that the heat release in region B is caused by contact combustion. In addition, region C covers the remainder of the burning regions. Of course, such a division is very rough, but we can obtain an approximate proportion of the heat release of the contact surface. Table 3 compares the heat release of case 1 and case 2 in the three regions. The combustion near the detonation wave (area A, which includes almost all detonation combustion and deflagration following the detonation wave) has a heat release ratio of approximately 50%. The percentage of heat release in the contact combustion region is slightly more than 30%, which is somewhat higher than Zhou and Wang's calculated value of 23.6% in the two-dimensional premixed case [30]. The possible cause should be because the contact region is larger in the non-premixed case than in the premixed case. Finally, it can be qualitatively seen from Figure 13 that the second contact burning (CB2) accounts for the major proportion of the heat release of contact combustion.

Conclusions
In this paper, the unsteady RANS solver is used to calculate the operating conditions of an air-breathing rotating detonation engine on a flight of Ma = 3 at an altitude of 10 km. Two configurations are considered. The characteristics of the DWs and contact burning are revealed in this paper. The main conclusions are described as follows.
A single-wave mode detonation wave is obtained in both cases with different isolator expansion angles. The expansion angle has a great influence on the shape of the detonation wave but has little influence on the detonation velocity.
In addition, the formation process of the buffer region was revealed, inside of which is a fuel-rich region. The fuel inside of the OH region continuously mixes and burns with the oxygen outside of the OH region, which maintains the existence of the OH region, and finally, the OH is exhausted in the region between the oblique shock wave and before the slip line.
Furthermore, a more accurate schematic of a nonpremixed air-breathing RDE is given in this paper. The CB1 (label 4 in Figure 11) and the CB2 (label 7 in Figure 11) start very close to each other and become further apart as time passes. The continuous burning area (label 11 in Figure 11) between the oblique shock wave and the slip line is highlighted, where the OH and the heat release are continuously produced.
Last, the heat release of contact burning accounts for approximately one-third of the total heat release, and the main heat release comes from CB2 because the fresh mixture is continually ignited by the hot products. Conversely, about one-third of the fuel is consumed by the hot products in the detonation combustor. The configuration has little effect on the heat release of contact burning.

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

Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.