Numerical Study of the Effects of Injection Conditions on Rotating Detonation Engine Propulsive Performance

: A three-dimensional upwind conservation element and solution element method (CESE) in cylindrical coordinates is ﬁrst developed to effectively solve the unsteady reactive Euler equations governing a hydrogen–air rotating detonation engine (RDE) with coaxial structures. The effects of the annular width on the structure of the detonation front and the relationship between the thrust and mass ﬂow rate are then investigated. Additionally, RDEs with various injection conditions are systematically analyzed regarding ﬂow patterns and propulsion performance. The results reveal a positive correlation between the speciﬁc impulse and the area ratio of the injection slot to the head-end wall. Nevertheless, the speciﬁc impulse shows minimal dependence on the injector slot’s location when the area ratio is constant. Ultimately, it is concluded that the area ratio between the injector and the head-end wall is critical in determining the loss of speciﬁc impulse.


Introduction
The rotating detonation engine (RDE) has attracted significant interest due to its potentially high propulsive efficiency.Since its early conceptual development, significant experimental [1,2] and numerical [3][4][5] advances have been made in recent years to enrich the fundamental understanding of the stability and performance of RDEs.
Most RDEs typically employ a structure consisting of coaxial cylinders.It is also important to note that the concept of non-circular RDEs has recently been proposed and successfully demonstrated in experimental studies [6,7].Numerically, if the combustor width is small compared to the size of the RDE, it can be unwrapped as a two-dimensional (2D) problem.Based on this assumption, many numerical studies have been carried out.For instance, Schwer and Kailasanath [8] demonstrated the viability of RDEs for both hydrogen and hydrocarbon fuels.In addition, Tsuboi et al. [9] studied the performance of RDEs for a stoichiometric hydrogen-oxygen gas mixture; the specific impulse (I sp ) and thrust were correlated with the mass flow rate.Zhao and Zhang [10] focused on the chaotic propagation of detonation with discrete reactant inlets, where explosive spots were found arising from the interaction between the shock waves and the deflagrative fronts.Furthermore, solvers for two-phase RDEs have also been developed to study the effects of the fuel droplet diameter and pre-vaporization [11], the effects of the pre-vaporized gas temperature and equivalence ratio [12], and the limit of operational stability [13].If the combustor width is not negligible, or the geometry of the RDE is a hollow cylinder [14][15][16], the threedimensional (3D) effect should be modeled.For RDEs with a wide combustor chamber, Katta et al. [17] reported that compression emerges near the outer wall, and expansion appears at the inner wall, making the detonation wave sensitive to the combustor channel width.Nevertheless, the influence of the combustor width on the propulsion performance was not reported.In a series of experimental studies conducted by Kawasaki et al. [18], the radius of the outer wall was fixed, and the injectors were located near the outer wall.The measured specific impulse decreased when reducing the radius of the inner wall.This deficit was concluded as a loss through deflagration.Despite significant advancements, our understanding of rotating detonation waves still needs to be improved for a comprehensive analysis.As observed in these studies, the wave structures and engine performance can be affected by the geometry of the combustor and the positions of the injectors.Towards this end, one objective of the current study is to identify the wave dynamics and engine performance of different injector positions at varied injector to head-end wall area ratios of the RDE.
According to the characteristics of the coaxial cylinders of the RDE geometry under consideration in the present study, the space-time conservation element and solution element method (CESE) adapted to cylindrical coordinates are developed.As a scheme with a compact stencil, the CESE method with an arbitrary order of accuracy only requires neighboring cell information for the update of a cell.This characteristic greatly benefits parallel computation.In a study by Jiang et al. [19], the computational efficiency of a 2D secondorder CESE scheme was compared to that of a finite volume method (FVM) that utilized the second-order monotonic upstream-centered scheme for conservation law (MUSCL) reconstruction, the van Leer limiter, and the HLLC Riemann solver.Results showed that, under the same resolution, the CESE scheme exhibited superior computational efficiency compared to the FVM method mentioned.Additionally, the CESE scheme demonstrated better capturing of narrower contact discontinuities and small flow structures, which is crucial for simulations involving highly compressible reactive flows.The original CESE method was proposed by Chang [20], and it unifies the treatments of space and time with a staggered advancing strategy, where both the conserved variables and their spatial derivatives are computed simultaneously.Recently, a characteristic (upwind) CESE method was developed by Shen et al. [21,22] to overcome the excessive dissipation of the original central CESE methods.In the upwind versions, the marching scheme for the conserved variables is the same as in the original CESE scheme.When computing the spatial derivatives, the equation of the non-dissipative CESE scheme is utilized.Meanwhile, Riemann problems are solved to determine the inner fluxes [21][22][23].These updated CESE schemes have shown small numerical dissipation even with minimal CFL numbers.Numerical tests have shown their satisfying conservativeness and accuracy.This method has been widely extended to solve hypersonic viscous flows [23], detonation waves [24][25][26], magnetohydrodynamics [27], shock-bubble interactions [28], and two-phase flow problems [29].The CESE schemes have been applied to solve the above flow dynamics mainly using either rectangular or unstructured meshes.However, for the simulation of rotating detonation waves in coaxial cylinders, it would be preferable to employ an extended version of the upwind CESE scheme that incorporates coordinate transformation.
The rest of the manuscript is organized as follows.The governing equations and physical models for a 3D RDE are first introduced in Section 2. A meticulous extension of the numerical scheme is presented in Section 3, where the 3D upwind CESE method is applied to Euler equations in cylindrical coordinates.In Section 4, the effects of the geometrical size of the RDE and the locations of the injection slots are discussed both in regard to wave structures and propulsion performance.Conclusions are summarized in Section 5.

Physical Model
Following previous studies [30][31][32], neglecting the effects of viscosity, thermal conduction, and molecular diffusion, the flow field of an RDE is governed by the three-dimensional, reactive Euler equations: Here, U = [ρ, ρu, ρv, ρw, E, ρλ] T is the vector of conservative variables, and ρ, u, v, w, p, E, λ are the density, velocities, pressure, total energy, and mass fraction of the reactant, respectively.The vectors of flux and the vector of the source are expressed as where ω = −Kρλe − Ea RT denotes the mass production rate of the reactant governed by the one-step reaction model for hydrogen-air detonation [4,33].K, E a , R, and T are the preexponential factor, activation energy, gas constant, and temperature, respectively.The equation of state for a perfect gas p = ρRT is adopted here, and the energy per unit volume is defined as where γ is the specific heat ratio.In the present study, the thermodynamic and chemical parameters for the reaction from Ma et al. [33] are used as follows: γ = 1.29,K = 7.5 × 10 9 /s, E a = 4.794 × 10 6 J/kg, R = 368.9J/(kg • K), and heat release Q = 2.72 × 10 6 J/kg.A Strang operator splitting method is used to split the computation of flow and reactions.The integration of the reaction source term is achieved implicitly to avoid the stiffness problem.
For continuous rotating detonation waves that propagate within an annulus combustor, four types of injection arrangements are considered, as schematized in Figure 1.For ideal (full) injection, a combustible gas can be injected through any part of the head-end wall.For non-ideal (inner, outer, and middle) injection, a combustible gas can be only injected through a slot on the head-end wall.This is to replicate certain experimental configurations where injectors are positioned exclusively near the outer wall [18].Although the injectors in experiments are commonly angled about the sidewall, and the non-premixed phenomenon is crucial, these features are not within the scope of this study.Here, the RDE is assumed to operate without any exit nozzle effects [30,34] and directly exhaust the burnt gas with a back pressure of 1 atm.The combustible gas is assumed to be injected into the manifold through Laval micronozzles, and the boundary condition at the combustor inlet is specified according to the local pressure near the wall p w and the inlet stagnation pressure p 0 : (a) if p w ≥ p 0 , there is no injection, and the slip wall boundary will be implemented; (b) if p cr1 ≤ p w < p 0 , the inflow velocity is locally subsonic; (c) if p cr2 ≤ p w < p cr1 , the throat of the nozzle maintains choking conditions, and the injection is subsonic; (d) if p w < p cr2 , the injection is supersonic and not affected by p w .
The detailed calculation of the critical pressures, p cr1 and p cr2 , can be found in the study by Zhdan et al. [35].The non-reflective boundary condition is imposed at the exhaust.The adiabatic slip wall condition is implemented on the inner and outer chamber walls.

Extension of 3D CESE Method to Cylindrical Coordinate System
Here, we develop the second-order 3D upwind CESE scheme in cylindrical coordinates.For a more comprehensive review of the development of other CESE schemes, the readers can refer to Jiang et al. [19] and Wen et al. [36].
By disregarding the source term due to reaction, the governing Equation (1) in the physical domain is transformed into the computational domain with uniform meshes (Figure 2) with the mapping relation The quantities that are labeled with a tilde denote the variables in the computational domain, and U = |J|U, where J ≡ ∂(x,y,z) ∂(ξ,η,ζ) is the transformation Jacobian.For a typical cell in 3D cylindrical coordinates with the resolution of ∆r, r∆θ, and ∆z along the radial, azimuthal, and axial directions, x = ξ cos η, y = ξ sin η, z = ζ.Thus J can be analytically computed.As a result, the above equations can be further deduced as By implementing the product rule, the governing equation can be rearranged as where S g = − F/ξ is the source term due to geometrical transformation.Equation ( 7) is the exact equation that needs to be solved.In the CESE scheme, the conservative variables and their spatial derivatives are required and must be updated.If the spatial derivatives of U, defined as U m ≡ U m in the mth direction (m = ξ, η, ζ), are known, the derivatives of fluxes, i.e., F m , G m , and H m in the physical domain, as well as F m , G m , H m , and S g,m in the computational domain, can be straightforwardly computed by the chain rule.The temporal derivative U t is then obtained using the Cauchy-Lowalevski procedure as Then, F t , G t , H t , and S g,t are well defined in a similar manner as when computing the spatial derivatives.The detailed descriptions and formulations of the 3D LLF CESE can be found in Appendix A. The Sedov blast wave problem and 1D detonation are used as benchmark tests to validate the present 3D upwind CESE scheme in cylindrical coordinates.Details of the tests are provided in Appendix B.

Results and Discussion
The length of the RDE model in the present study is held constant at 40 mm, the outer radius (R o ) is 20 mm, the stagnation temperature is fixed at 300 K, and the ambient pressure is 1 atm, while the stagnation pressure (p 0 ), the area ratio of the injector throat to the injector exit of the Laval micronozzles (A * /A), the annulus width (w), and the area ratio of the injection slot to the head-end wall (α) are adjusted accordingly.Tables 1-4 summarize the conditions for different cases, and these details will be further discussed in the following sections.The inner radius is limited to not less than 11 mm to ensure that the resolution near the inner wall does not change much compared to that near the outer wall.All cases are run with 15∼23 operation cycles depending on different setups to ensure that the data are analyzed when the RDEs are statistically stable.The annulus width varies within 1∼9 mm, and the area ratio of the injection slot to the head-end wall changes from 0.33 to 1 with different positions of the injector slot.It is noteworthy that Sow et al. [37] attributed the deficit in detonation velocity observed in an extremely narrow 2D channel to the viscous boundary layer.Since most cases involve a relatively wide channel, to maintain consistency, the viscous effects were not considered in the present study.To preface, all of the tested conditions in the present study demonstrate the presence of continuous rotating detonation waves.

Grid Sensitivity Analysis
To ensure the accuracy and reliability of the numerical results, a convergence study is performed for a combustor of 1 mm width with three sets of mesh grids (Table 1), M1: ∆r = 0.5 mm, R o ∆θ = 0.4 mm, ∆z = 0.4 mm; M2: ∆r = 0.2 mm, R o ∆θ = 0.2 mm, ∆z = 0.2 mm; and M3: ∆r = 0.1 mm, R o ∆θ = 0.1 mm, ∆z = 0.1 mm.Ignition is attained by a segment of a detonation wave that propagates circumferentially near the inlet.Usually, stable propagation can be established after a few cycles of propagation.Figure 3 shows the temporal traces of pressure at a point located at the outer wall and 1 mm away from the inlet (noting that no attempt has been made to align the pressure profiles in these cases).It is shown that the pressure periodically oscillates after about four cycles of propagation.The detonation velocity averaged from five stable cycles at the outer wall for M2 case is 1976 m/s, which is close to the Chapman-Jouguet (CJ) value [4].Overall, no significant difference is observed for these three resolutions, except for the peak pressure values.Furthermore, all results can produce similar macroscopic structures, as depicted in Figure 4.For the consideration of both the accuracy and computational efficiency, M2 is adopted as the resolution for the remaining simulations.This resolution is on a similar scale to previous studies on hydrogen-air RDEs [17,38].It is indeed important to note that the resolution used in our study is still insufficient to fully capture the detailed cellular structures present in RDEs.To accurately capture these structures, a resolution of at least one order smaller than the reaction zone is typically required, as suggested by previous studies [24].Achieving such a fine resolution for a wide channel would demand computational resources surpassing the current capabilities.Therefore, in the present study, our focus was primarily on investigating the macro flow structures and overall performance of the RDE.Future research with more advanced computational resources could aim to explore and analyze these structures in greater detail.In the present study, for a typical case with a channel width of 9 mm, approximately 4.89 million computational cells were used.A total of 25 cases were conducted, with each case involving 15∼23 cycles of detona-tion propagation to ensure reliability.For a typical case with a 9 mm channel width and 4.89 million computational cells, approximately 211 CPU-hours were required per one cycle of detonation wave propagation.All cases were performed in parallel with MPI standards on the Tianhe supercomputer in the Tianjin Supercomputer Center, China.

Effect of Channel Width with Full Injection
To understand the effects of the increasing combustor width on the wave structures and performance of the RDE, the temperature contours of four channel widths (1 mm, 3 mm, 6 mm, and 9 mm) are compared in Figure 5.The test condition is as shown in the first row of Table 2.The injector fully occupies the head-end wall in these cases.For channel widths of 3 mm, the flow fields are close to that of 1 mm.When the channel width is increased to 6 mm, the detonation front is inclined to the outer wall and reflected shocks are generated.The strong expansion near the inner chamber wall leads to the significant curvature of the detonation wave front.The velocity component normal to the wave front is lower than the CJ value near the inner wall but approaches CJ as it approaches the outer wall [39,40].As a result, the velocity measured along the outer wall may be more appropriately referred to as the velocity of the wave signal rather than the detonation velocity.This curved detonation wave induces a series of shock waves that bounce back and forth, similar to the phenomena observed in the study conducted by Katta et al. [17].These shock trains can be barely observed in the cases with narrower channels.
By further increasing the channel width to 9 mm, the distance between the leading detonation wave and the first reflected shock near the inner wall increases.This increased distance allows for greater pressure relaxation below the critical pressure, resulting in the faster injection of fresh combustible gas near the inner wall compared to the outer wall due to compression by the reflected shock.Furthermore, it can be observed that wider chambers result in shorter propagating detonation waves.The recorded heights of the detonation waves, ranging from narrow to wide chambers, are measured as follows: 12.2 mm, 11.7 mm, 10.4 mm, and 8.8 mm, respectively.Per the authors' knowledge, there is limited numerical research available on the impact of various annulus sizes on performance.The comparison of the propulsive performance is achieved in terms of the thrust (F) and mixture-based specific impulse (I sp ), defined by where v z is the flow velocity normal to the RDE exit, A is the area of the exit, ṁ is the mass flow rate of the combustible gas, and g is the acceleration of gravity.Table 5 lists the propulsive performance for different combustor widths.In a highly intuitive manner, it is noted that increasing the channel width leads to an observed increase in thrust.However, the specific impulse remains relatively unchanged.This observation suggests that the thrust exhibits a linear relationship with the mass flow rate.Nevertheless, there is a slight variation in the standard deviation (SD) of F and I sp , suggesting that the RDEs maintain stable operation across all channel widths.Note that when increasing the channel width, the detonation velocity at the inner wall (V i ) drops due to the vigorous expansion caused by the convex surface.In contrast, the velocity of the wave signal at the outer wall (V o ) increases.Overall, these simulations suggest that although the compression and expansion near the inner and outer walls significantly modify the wave structures, their influence on performance is negligible.

Effect of Inner Radius with Fixed Injection Area
In the experimental apparatus, the injectors are sometimes located only near the outer wall [18,41].In this section, we investigate the impact of the inner radius on the system while maintaining a constant injection area near the outer wall.We vary the inner radius within a certain range to analyze its effects on the overall system performance.By keeping the injection area constant (the width of the injection slot is w inj = 3 mm near the outer wall), we can isolate the influence of the inner radius and gain insights into its role in the system dynamics.The test conditions are summarized in Table 2.A typical flow field with the inner radius of 11 mm is shown in Figure 6.In this case, as the injector slot is near the outer surface, the gas near the inner wall is filled with the burned gas from previous cycles, forming a stratification structure.Similar to the findings in gaseous layered detonations bounded by an inert gas with a small inert to reactant acoustic impedance ratio [42], the difference in acoustic impedance results in a detached shock (DS) in the burned gas that propagates ahead of the detonation wave.A transmitted oblique shock (TS) links the detached shock with the detonation front.This transmitted shock seems inadequate to directly ignite the gas.Additionally, from the pressure history at the outer surfaces of the three cases (Figure 7a), the pressure peak is weaker for a smaller inner radius, and the period of two successive peaks decreases (indicating higher rotating frequency).The specific impulse also decreases for a smaller inner radius.It is further noted that for cases with w = 6 mm and w = 9 mm, the specific impulse strongly fluctuates (Figure 7b), with standard deviations of 4.4 s and 5.0 s, respectively.

Effect of Injector Location and Size
To further analyze the effect of the injector location, the injection patterns depicted in Figure 1 are examined.The test conditions are listed in Table 3.While changing the location of the injection slots, the widths of the injection slots are adjusted accordingly to keep the area ratio of the injector to the head-end wall fixed.We resort to the cases with an area ratio of 0.5, of which the pressure contours with density lines are plotted in Figure 8.When the injector slot is positioned close to the inner wall, there are notable differences compared to the case with an outer injector.The detonation wavefront appears visually straight (Figure 8c), and there is a reduction in the distance between the detached shock and the detonation front.From close observation of the case in which the injector slot is located in the middle, DS-TS structures exist on both sides of the detonation wave (Figure 8d).All three cases demonstrate relatively stronger fluctuations in specific impulses compared to the scenario with full injection.The standard deviations for the inner, middle, and outer injector cases are 5.6 s, 2.1 s, and 4.1 s, respectively.In addition, we utilize the total pressure as a measure of the system's capacity to perform work.This is calculated using the equation where M is the Mach number.However, it is crucial to accurately define the total pressure for an unsteady system [30].To achieve this, we adopt the mass-weighted average total pressure (MAP) definition proposed by Wen et al. [43].This Eulerian-based definition has been shown to be equivalent to the Lagrangian-based definition in stable RDEs.This approach enables us to accurately provide a reliable metric for comparison between different injection patterns.The MAP is obtained through integration over a specified surface: P t (z) = S ρ(x, y, z)ω(x, y, z)P t (x, y, z)dA S ρ(x, y, z)ω(x, y, z)dA (11) As depicted in Figure 9, all four cases, whether full or partial injections, exhibit a higher total pressure at the exit (z = 40 mm) of the combustor compared to the inlet (z = 0 mm), a testament to the pressure gain combustion.In the case of full injection, there is approximately a 67% pressure gain in comparison to the total pressure near the inlet.The total pressure increases from the inlet until reaching a peak, which roughly corresponds to the triple point of the detonation wave with the contact discontinuity.Following the peak, the total pressure gradually decreases due to nonisentropic expansion and interactions with shocks.This variation in total pressure aligns well with the trends observed in the theoretical work by Wen et al. [43].However, when comparing the total pressure to the stagnation pressure of the gas reservoir, there is still approximately a 29% pressure loss, even in the case of full injection.This can be attributed to the characteristics of the injectors used.Only injectors with relatively large area ratios can provide significant pressure gains [44], but this also increases the risk of poor isolation of the mixture in the plenum from the gas and waves originating from the combustor in practical applications.Careful design is necessary to achieve a total pressure gain while preventing backflow from the hot gas to the supply manifold [45][46][47].To demonstrate, in the case of full injection using a micro-Laval nozzle with an area ratio of A * /A = 0.5, the pressure loss across the inlet boundary is significantly reduced, resulting in a net pressure gain of 32% compared to the stagnation pressure in the reservoir.Additionally, it should be noted that the full injection case with a 6 mm channel width exhibits a similar total pressure distribution to that of the 9 mm case.This further supports the previous finding that the channel width alone does not have a significant impact on the performance of the RDE.The effect of the injection location on the total pressure can be explained by analyzing time-averaged contours across the combustor channel and variable distributions near the head-end wall.Figure 10 represents the time-averaged contour, which is obtained by averaging the data azimuthally over the entire RDE.In the case of full injection, the injected reactant and the heat release are nearly uniform along the radial direction.However, the partial injection cases show varying degrees of inhomogeneities.This indicates that, apart from the lack of confinement in the axial direction, the detonation waves in these partial injection cases also suffer additional losses due to the presence of bounded inert gas in the radial direction.This observation is further supported by the measured detonation velocity, where the case with middle injection demonstrates approximately 91% of the velocity achieved in the full injection case.Figure 11 plots the total pressure and injection Mach number distributions along the centerline at the bottom of the channel, with the detonation propagating from left to right (azimuthally).The weaker detonation wave in the middle injection case leads to lower pressure compared to the others.In the studied cases, the micro-Laval nozzles at the inlet boundary operate at the choking condition.However, with the collective impact of the velocity deficit, the extra expansion caused by the bounded inert layers after passing through the detonation wave, and the absence of strong reflected shock trains as in the full injection case, the pressure rapidly drops to a lower value (Figure 11b).Consequently, the partial injection cases experience a stronger total pressure loss in the inlet gas.In summary, it can be stated that the confinement provided by hot burnt gas, as opposed to a rigid wall, in partial injection scenarios leads to a reduction in post-detonation pressure and an elevation in total pressure loss.The performance of the RDEs when the area ratio of the injector to the head-end wall ranged from 0.33 to 0.75 and with different injector locations is summarized in Figure 12, along with the data of the full injection case.Once again, the RDEs with non-ideal injection conditions are prone to large fluctuations in specific impulse.It is also evident that the area ratio between the injection slot and the head-end wall is the most critical factor in determining the specific impulse.In contrast, the injector locations appear to cause minor changes.

Effect of Stagnation Pressure and Area Ratio of the Injector Nozzle
In this section, RDEs with a stagnation pressure of 20 atm or with a nozzle area ratio of 0.1 are simulated (Table 4).Figure 13 displays the specific impulses corresponding to different area ratios of the injector to the head-end wall.The observation of similar trends in all three scenarios suggests a conjectured universality regarding the positive correlation between the area ratio of the injection slot to the head-end wall and the loss in specific impulse, as discussed in the previous section.It is worth noting that for higher stagnation pressure, the percentage of specific impulse loss compared to the full injection scenario is smaller.To elucidate, at an area ratio α = 0.33, the specific impulse for the "20 atm − A * /A = 0.2" scenario is about 79% compared to that with the ideal injector, while it is 65% for the "10 atm − A * /A = 0.2" scenario and 54% for the "10 atm − A * /A = 0.1" scenario.

Conclusions
The present paper extended the three-dimensional upwind CESE method to solve the RDE problems in cylindrical coordinates.In particular, the effects of the combustor width and area ratio of the injection slot to the head-end wall on the dynamic features and propulsion performance were quantitatively determined.For the cases with ideal injection, by increasing the combustor width, the velocity measured on the outer wall is increased.It should be noted that this velocity does not correspond to the detonation velocity as the wavefront is inclined to the wall.On the other hand, the wavefront near the inner wall appears to be nearly perpendicular to it.Moreover, the measured detonation velocity in this region is observed to be smaller than the CJ value.This velocity deficit is more significant when decreasing the radius of the inner wall.However, there seems to be no significant influence on the propulsion performance.For the cases where the injector occupies only part of the head-end wall, due to the difference in the impedance of the burned gas and fresh reactant, the wave structures are different from the cases with ideal injectors, and the stratification of burned and fresh gas can be observed near the head-end wall.The total pressure loss for the cases with partial injection can be primarily attributed to the presence of bound inert gas.Furthermore, the specific impulse decreases compared to the ideal cases.With the geometry of the RDE fixed, the specific impulse is measured to be mainly dependent on the area ratio between the injector slot and the head-end wall, regardless of whether the injector slot is located at the inner, outer, or middle head-end wall.The drop in specific impulse manifests less severely at higher stagnation pressure.These numerically gathered data are anticipated to provide useful insights for the formulation of guidelines that can aid in the design of practical RDEs.

Appendix A.1. Upwind CESE Method on 3D Meshes
The CESE unifies the treatment of space and time and uses a staggered time marching strategy.To be specific, each timestep is divided into two half-steps, i.e., the 1st and 2nd half-steps.Initially, the values (both conservative variables and their spatial derivatives) are known at the vertices of each computational cell, e.g., A1∼A8 (at step n) in Figure A1a.For the 1st half-step, the values at the center of the cell (P * ) need to be computed.For the 2nd half-step, the values are updated at the vertices (A i ) based on the values at the centers (P i ).The approaches in advancing these two half-steps are in a similar manner.Without loss of generality, here, we illustrate the procedure for advancing the 1st half-step.For easy description, B1∼B12 are the midpoints of corresponding grid edges, and C1∼C6 are the centroids of corresponding surfaces.We use superscript * to denote the grid point in the next half-step (step n + 1/2).The conservation element (CE) of the 3D CESE is defined as a Euclidean four-dimensional (4D) space-time region (Figure A1b).For example, the CE corresponding to P * (denoted as CE(P Consequently, the term surface of a CE in the context mentioned below refers to a 3D hexahedron, e.g., as the right-side surface of CE(P * ).Furthermore, it is worth noting that every CE can be subdivided into eight subordinate CEs (sub-CE).Each sub-CE is associated with one vertex A i and encompasses eight surfaces within the four-dimensional space (Figure A1c).For instance, the sub-CE associated with and inner surfaces that are separated from the neighbor sub-CEs, including ( 6 , where superscript * * denotes the point at step n + 1.Within each SE, the conserved variables and fluxes are assumed to be linearly distributed.where H = ( F, G, H, U), and ds = dδ • n, with dδ being an infinitesimal area and n the corresponding unit outward normal vector.The integration over the surfaces of each sub-CE can be achieved with the multiplication of the surface area and the average flux on that surface.For the second-order scheme, the average flux is the value at the centroid of the surface and is computed by the first-order Taylor expansion.For the purpose of clarity and simplicity in our illustrations, we will omit the use of hat symbols.The dimensions in the subsequent paragraphs will also be denoted as x, y, and z.Therefore, the integration of the surfaces of the sub-CE associated with A 2 can be written as approximated by Taylor expansions from A 1 and A 2 , respectively.In order to suppress spurious oscillations, U 1,x and U 2,x are reconstructed with the WBAP limiter [22,48].It should be noted that the values of derivatives are only stored in the computational domain, and it is not necessary to transform back into the physical domain.With the estimations of U l and U r , F l , F r , and F L,2 can be easily computed.The computation procedures of G B,2 and H U,2 bear resemblance to those of F L,2 and are omitted here.
Likewise, the discrete form for the ith sub-CE (i = 1∼8) is summarized here: Moreover, the inner fluxes across two neighboring sub-CEs have the same magnitude, e.g., F R,1 = F L,2 .By adding the discrete Equation (A8) of all the sub-CEs, the inner flux terms are balanced, x (P * ) = Using the same approach, one can derive U (2) x (P * ) by U * 4 and U * 3 , U x (P * ) by U * 5 and U * 6 , and U x (P * ) by U * 8 and U * 7 .Then, U x (P * ) is finally computed with a limited average:

Figure 1 .
(a) Schematic model of an RDE with an annulus combustion chamber, and injection patterns on the head-end wall, (b) ideal, (c) inner, (d) outer, and (e) middle injector slots.

Figure 2 .
Figure 2. Mesh in (a) physical domain and (b) computational domain.

Figure 3 .Figure 4 .
Figure 3. Temporal traces of pressure at a point 1 mm away from the inlet.

Figure 9 .
Figure9.The distribution of total pressure along the axial direction.Typical pressure gain relative to the total pressure near the inlet is provided, together with that relative to the stagnation pressure in the reservoir in the brackets.

Figure 10 .Figure 11 .
Figure 10.(a) Averaged density profile along the radial-axial plane and (b) averaged heat release rate for channel width of 9 mm with different injection patterns.The dashed lines indicate the top of the detonation wave.

Figure 12 .
Figure 12. Specific impulses as a function of area ratio of the injector to the head-end wall with different injector locations.P 0 = 10 atm, A * /A = 0.2.Error bars indicate the standard deviations.

Figure 13 .
Figure 13.Influence of stagnation pressure and area ratio of the injector nozzle on specific impulse.For all the cases, the injector slot is located near the outer wall.Error bars indicate the standard deviations.
(a) Mesh arrangement in three-dimensional CESE method.A denotes the vertex; B and C denote the centers of line segments and surfaces, respectively; P is the staggered point that needs to be solved in the 1st half-step.(b) Definitions of CE.(c) Definition of sub-CE, with one outer flux F R,2 and one inner flux F L,2 .Below, we impose the space-time flux balance over the mth sub-CE,

Figure A3 .
Figure A3.Pressure history of the 1D detonation wave.

Table 1 .
Simulation conditions for the effects of grid resolution.

Table 2 .
Simulation conditions for the effects of annulus width and inner radius.

Table 3 .
Simulation conditions for the effects of injector location and size.

Table 4 .
Simulation conditions for the effects of stagnation pressure and area ratio of nozzle.

Table 5 .
Performance of RDE with different combustor widths.