The first case studied is the annular RDC designed and operated at Nagoya University by Ishihara et al. (2017). The RDC inner and outer diameters are 62 mm and 78 mm respectively, resulting in an 8 mm annulus width, while the length of the chamber is 70 mm. A converging-diverging nozzle is attached at the exit of the combustion chamber. The contraction ratio at the throat is 2.5. In the present simulation, the nozzle is not accounted for since its effect is also modeled by the initial condition. On the other hand, the length of the combustion chamber is greatly increased to reduce the potential reflective effect from the outlet boundary to the injection plane.

C₂H₄ and O₂ are injected through 120 injection elements to obtain a global equivalence ratio of 0.9, corresponding to the test of Ishihara et al. (2017). The chamber mass flux is 120 kg/m²/s. One injection element is composed of two perpendicular injection tubes of 1 mm in diameter. The two tubes are inclined at 45° with respect to the chamber axis and they are radially aligned to provide a direct impact of the unlike jets.

In addition to the detonation speed, experimental data on the mean pressure at various locations in the RDC are available in the study of Ishihara et al. (2017) and will be used as reference values to correctly set the initial condition. In the following, a parametric study is conducted to determine the effect of the Riemann problem parameters on the simulation results.

The effect of the two states of the Riemann problem is investigated in the Nagoya RDC case. In this RDC, detonation does not propagate at the theoretical CJ speed and it does not produce the theoretical CJ pressure, mainly because of the heterogeneities in the fresh mixture. Therefore, two reinjection cases are compared in Figure 4A: Case 1, for which the states of the Riemann problem are determined from the exact CJ conditions, and Case 2, for which the states are determined from a deteriorated CJ condition ( ( . ) ̃ conditions) accounting for certain unmixedness of fresh propellants, obtained from the in-house code DetonHeter (see Supplementary Material). The corresponding upper and lower initial states of the Riemann problem are summarized in Table 2 for the two cases. Since the initial condition of Case 2 is obtained from a deteriorated CJ condition, the initial pressure is logically below the pressure of Case 1. In the following, h _D is arbitrarily fixed at 10 mm.

Figure 4A presents pressure evolution at the injection plane from the reinjection simulations. Pressure in Case 2 decreases slower than in Case 1 and reaches its final value around t = 100 µs, while pressure in Case 1 drops faster until it reaches a steady value around t = 25 µs. Thus, the conditions of the lower and upper states have a strong impact on the injection dynamics. The mean pressure in these two simulations is either too high (5.4 bar in Case 1) or too low (1.2 bar in Case 2) compared to the experimental pressure measured at the injection plane (approximately 3.7 bar).

For the next cases, the initial condition obtained for Case 2 (deteriorated CJ condition for the lower state values) is selected because it seems to produce a pressure decrease better suited to represent the mean experimental pressure.

The second parameter studied here is the upper state Mach number (M _y,up ), which is linked to the expansion of the burnt gases produced by the previous detonation wave. The pressure decrease at the injection plane in the simulations is plotted in Figure 4B for three different values of M _y,up , the other parameters from Case 2 remaining unchanged.

From the results shown in Figure 4B, M _y,up has an effect on the final pressure: when M _y,up is reduced, burnt gases expansion is weaker and the pressure at the injection plane has a higher final value. In addition, the pressure decrease shown here is quite different from the classical exponential decay illustrated in (Section 3.2.2). Consequently, M _y,up is not a suitable parameter to control the rate of the pressure decrease in spite of its effect on the final pressure at the injection plane.

The effect of M _y,up is also shown in Figure 5 with the field of Z _excess . If Z _excess is positive (respectively negative), it corresponds to the volume fraction of C₂H₄ (respectively O₂) in excess relative to stoichiometry. The two fields obtained at the end of the reinjection period τ are very different from each other. Since the final pressure in the chamber increases when M _y,up is reduced, the velocity of the propellant jets in the chamber is also reduced for the same injected mass flow rates. Thus, the fresh gases have more time to mix close to the injection plane. This leads to a better mixing in the case of M _y,up = 1 as shown in Figure 5, on the right. Also in the same reinjection simulation, the O₂ jet impinges the inner wall of the engine, creating a lean mixture close to it. This effect is due to the fact that the radial component of the O₂ jet momentum is about four times higher than that of the C₂H₄ jet (the estimation is based on the mass fluxes of the injected propellant flows, which are 470 kg/m²/s for C₂H₄ and 1800 kg/m²/s for O₂). This behavior is consistent with the observation of Ishihara et al. (2017), that a part of the wall coating in the experimental chamber was oxidized during the RDC operation.

The second geometry studied is an annular RDC (named as ONERA RDC in the following), in which the injector is composed of elements described in the patent of Davidenko and Gaillard. (2022). H₂ and O₂ are used as propellants for this RDC. Without any preliminary information on the operational conditions in this RDC, they are evaluated from a 2D simulation as a first approximation.

The numerical approach used for the 2D RDC simulations has been detailed by Gaillard et al. (2017), and an example of the simulated flowfield is shown in Figure 6. Premixed propellants are uniformly injected through a slot at the bottom of the chamber. This kind of simulation serves to evaluate a reference pressure decrease on the injection plane, for ideal injection process and detonation propagation.

Figure 7 shows a pressure decrease along the injection plane obtained from three 2D simulations, with the same mass flux of propellants as in Figure 6. The dimensionless abscissa x* is the circumferential position x divided by the spatial period between detonation waves, noted L _x in the following. The detonation wave is located at x * = x L x = 0 .

The only parameter that differs in these simulations is the number of detonation waves considered in the RDC and hence, the spatial period L _x . Nevertheless, the three pressure profiles are almost perfectly superimposed, meaning that for 2D premixed simulations, the period L _x does not have an impact on the pressure profile in the engine if the dimensionless coordinate x* is used. This similarity was also mentioned by Davidenko et al. (2007) concerning simulations using Euler equations.

According to this result, 2D simulations can be performed on a small computational domain (i.e., a small spatial period), leading to lower computational cost. The pressure profile describing the expansion process for any selected number of waves can then be derived from such a 2D simulation.

In determining the initial condition for a reinjection simulation, the position of the Riemann problem discontinuity, h _D , is a free parameter unless one knows the height of the detonation front. The method described in this section allows finding the right h _D for a given pressure variation. The 2D simulation in which 5 detonations were assumed is used to define the reference pressure profile. The pressure profiles obtained from reinjection simulations with h _D = 3 and 10 mm are compared with the pressure profile of the 2D case in Figure 8A.

It can be seen in Figure 8A that h _D has an impact only on the rate of the pressure decay. The initial pressure is the same, but increasing h _D slows down the expansion of burnt gases at the injection plane. Therefore, h _D is a key parameter to set a proper initial condition that can correctly model the effect of burnt gas expansion.

Now, a strategy to determine the correct value of h _D is proposed. As for Figure 7, the idea is to find a dimensionless coordinate to compare the pressure profiles obtained from reinjection simulations. The Riemann problem is supposed to model the discontinuity located at h _D , separating the burnt gases right behind a detonation, and the expanded burnt gases produced by the previous detonation. This discontinuity position is on the top of the fresh mixture layer in front of the detonation wave. If the period between two detonation waves increases, there is more time to inject propellants and the thickness of the fresh layer grows correspondingly. Therefore, it can be assumed that, like the thickness of the fresh layer, h _D is proportional to the time period between detonations (τ). Also, the temporal and spatial periods are linked as L _x = τ × V _D , where V _D is the detonation speed in the laboratory reference frame. Hence, the thickness of the fresh mixture layer is proportional to the spatial period (Equation 2), as experimentally shown by Bykovskii et al. (2006) and Rankin et al. (2017). h D ∝ L x (2) It was shown previously that the pressure profiles in 2D simulations with different L _x are similar in dimensionless coordinate x * = x L x . Thus, by using Equation 2 it comes: x * = x L x ∝ x h D (3) Then, the coordinate change t = x V D in Equation 3 leads to the following equation: x * ∝ t × V D h D = x + (4)

Pressure at the injection plane of the RDC is plotted against the dimensionless variable x ⁺ for the two reinjection simulations represented in Figure 8B. The two curves are perfectly superimposed. Thus, it is possible to predict the pressure profile for a reinjection simulation with a different h _D . In fact, one reinjection simulation (obtained with an assumed discontinuity height h _D,1) allows to plot the pressure profile P ( x 1 + = t V D h D , 1 ) . Then, it is possible to reconstruct the pressure profile P(t) that will be obtained for a different value of the discontinuity height h _D,2 with Equation 5: t = x 1 + h D , 2 V D (5)

To illustrate this method, the optimum height h _D = 6.5 mm is used to obtain a pressure profile, close to the one of the 2D simulation, plotted with green dashed-dotted line in Figure 8A.

The geometry studied hereafter is the annular chamber of the experimental RDC designed and operated at the Technische Universität of Berlin (TU Berlin) by Bach et al. (2020). The annulus inner and outer diameters are 74.8 mm and 90 mm respectively, resulting in a 7.6 mm radial gap, whereas the length of the chamber is 110 mm. Air is injected radially from the outer wall through a 1 mm slot and H₂ is injected through 100 evenly spaced 0.5 mm holes. The conditions in the chamber are taken or evaluated from the experiment. The global equivalence ratio is 1 and the mass flux related to the chamber cross-section is 100 kg/m²/s. The selected chamber geometry has an outlet restriction of 50%, creating a sonic throat for the selected operating point.

Firstly, a non-reactive reinjection simulation is carried out to make a comparison with an established injection. The lower state of the Riemann problem corresponds to the CJ condition for a stoichiometric H₂/Air mixture at 1 bar (approximately 50% of the mean experimental pressure) and ambient temperature (300 K). The discontinuity is set at h _D = 20 mm and the upper state Mach number is M _y,up = 1.5. The initial condition has been validated by comparing the experimental and numerical mean pressure measured 10 mm upstream of the throat (2.3 bar in both cases).

To study the mixing process in the chamber, the Z _i variables (proposed by Gaillard et al. (2019)) will be used below. Z _i are computed in every mesh cell and have a particular meaning: Z _st corresponds to the mixture volume fraction at stoichiometry (i.e., the volume fraction of propellants that can react completely), Z H 2 , e x (respectively Z O 2 , e x ) is the H₂ (respectively O₂) volume fraction in excess compared to stoichiometry, while Z _bg is the burnt gases volume fraction. Hence the sum of Z _i is equal to unity.

Figure 9 shows the Z H 2 , e x field at t = 5, 10 and 15 µs. Injection blocking and backflow of burnt gases in the fuel tube can be seen at 5 µs. The same phases happen for the Air slot (not shown here). At 10 µs, the chamber pressure falls below the H₂ supply pressure, allowing the injection of H₂. On the other hand, the Air slot is still blocked. This phenomenon creates an axial stratification of H₂ up to t = 15 µs.

The final state (see the top of Figure 10), corresponds to the mixture the detonation will consume, meaning that a good mixing is obtained above y = 5 mm and slightly closer to the inner wall of the engine (see Figure 10A). The rich zone created close to the outer wall at the beginning of the simulation can still be observed (see Figure 10B). It is located in a recirculation zone, in which some burnt gases from the previous detonation is stuck (see Figure 10D). Meanwhile, a lean mixture is formed in the recirculation zone located at the bottom of the chamber, and along the inner wall (see Figure 10C).

Mixing efficiency from the reinjection and established flow simulations can be compared as a function of axial coordinate. Mixing efficiency is defined by Equation 6. It corresponds to the mass fraction of fresh gases which are in proportions that respect the global equivalence ratio of the chosen operating point (ER _glob ). ER is the local equivalence ratio in a cell. η mix y , t = ∬ S y ρ Y H 2 max E R , E R glob d x d z ∬ S y ρ Y H 2 min E R , E R glob d x d z (6)

Instantaneous profiles of mixing efficiency at different instants from the reinjection simulation and the time-averaged mixing efficiency of the established injection are compared in Figure 11. Mixing efficiency increases with time and gets to a maximum at t = 120 µs. At the same time, the fresh mixture layer reaches y = 20 mm (the maximum axial distance shown in Figure 11). The mixing efficiency profiles of the reinjection simulation are far below the one of the established simulation. In fact, the established injection simulation was performed at a pressure of 2.1 bar, corresponding to the mean experimental pressure at the chamber outlet obtained by Bach et al. (2020). On the contrary, the chamber pressure is not constant in the reinjection simulation, and it falls below the mean outlet pressure after approximately 50 µs. This lower pressure increases the flow velocity in the chamber, thus the reactants mix further away from the bottom of the chamber (above 20 mm) compared to the continuous injection.

The Z _i and temperature fields obtained at the end of the reactive simulation are shown at the bottom of Figure 10. In comparison with the non-reactive simulation, the fraction of burnt gases indicated by Z _bg (Figure 10D) and the temperature (Figure 10E) are significantly higher within the recirculation zone above the Air slot where combustion is stabilized. At the same time, the flow dynamics in the recirculation zone change in the way that a larger area is filled with the stoichiometric mixture (Figure 10A).

The reactive simulation makes it possible to compute the deflagration losses in the combustion chamber. Multiple strategies can be used to determine the mass of fresh gases consumed by deflagration. Here, it was decided to use the mass balance of fresh gases between the inlet and outlet of the combustion chamber. The CEDRE software allows to save the instantaneous mass flow rates at the inlet and outlet boundaries ( m ̇ i n and m ̇ out ) and the mass integral in the computational domain (m _cmb ). Hence for each reactant, the rate of consumption by deflagration ( m ̇ def ) can be computed for each timestep with Equation 7. m cmb t + Δ t = m cmb t + ∫ t t + Δ t m ̇ i n − m ̇ out − m ̇ def d t (7)

The height of the combustion chamber is high enough to prevent fresh gases from leaving the domain. Therefore, the mass of each propellant in the RDC corresponds to the difference between the injected mass and the mass consumed by deflagration. It is then divided by the total mass injected during the simulation and shown in Figure 12.

The time delay between the start of H₂ and Air reinjections is visible. The deflagration losses are of the same order of magnitude for the two reactants: only 5% are consumed by deflagration. Therefore, the early consumption of fresh gases during the time period between detonations is small under the conditions of the present study.

The final flowfield from the reactive simulation is used as input for the in-house code DetonHeter, which can account for mixture heterogeneities to evaluate the detonation propagation speed (see Supplementary Material). The simulated flowfield is divided, along the distance from the injection plane, into 10 slices of 4 mm in height, and the computation is performed with the heterogeneous state of each slice.

The obtained detonation propagation speed is displayed in Figure 13. The abscissa represents the center of each slice. The experimental value of Bach et al. (2020) and the numerical speed range of Nassini. (2022) are shown for comparison. The simulation of Nassini. (2022) was performed on the TU Berlin RDC but at a different operating point. The computed detonation speed in Figure 13 presents a bell shape. In fact, fresh gases are not well mixed in the bottom of the chamber, inducing an important speed deficit. As shown before, mixing improves along the chamber height, thus increasing the detonation speed, which reaches a maximum at 26 mm above the injection wall. The following decrease is caused by the growth of the burnt gases mass fraction that remains on the top of the fresh mixture layer. The reduced model results are also coherent with the numerical work of Nassini (2022).

The experimental detonation propagation speed lies within the predicted range. Nevertheless, the experimental speed is about 400 m/s lower than the maximum estimated detonation velocity. RDC simulations often result in such discrepancies with experimental observations. Different reasons can be given, such as neglected viscous interactions on the walls or too simplified kinetic mechanism. Apart from the modeling aspects, we also propose a discussion in relation to the detonation cell size for the present RDC. Figure 14 compares the combustor annular gap with the experimental cell width for a similar H₂/Air mixture at different pressures. The cell size is close to the annular gap. In this particular case, the detonation cellular structure tends to change from 3D to 2D (Ishii et al. (2002)) and the effect of the boundary layer (friction and heat transfer) is non-negligible. Multiple experiments such as Ishii et al. (2002) or Xiao et al. (2021) showed that the speed deficit depends on the gap-to-cell width ratio and that it can reach 30% before detonation failure. The present detonation speed calculation does not account for these effects, which are probably important in the TU Berlin RDC and can be responsible for the observed detonation speed deficit. From an experimental point of view, it could be interesting to see if the detonation speed can be increased with a larger annular gap in the TU Berlin RDC.