NUMERICAL INVESTIGATION OF AN UNSTEADY INJECTION ADAPTED TO THE CONTINUOUS DETONATION WAVE ROCKET ENGINE OPERATION

Detonation applied to propulsion could result in a promising increase of the thermodynamic e©ciency of the engine cycle. Numerical simulations of the detonation propagating in the Continuous Detonation Wave Rocket Engine (CDWRE) are currently performed but still do not account for realistic injection process. The assumption of an ideal injected premix is generally chosen for convenience to obtain theoretical results. Comparison of the numerical results with experiments is di©cult because of the clear di¨erence of the injection con¦gurations. Some physical as-pects of the separate injection of the components used in experiments are not clearly assessed. This study is included in a wider numerical project aimed at designing and optimizing a realistic CDWRE. The optimization process is presently focused on the injector. One element of the injection hole pattern is considered assuming that this element is periodically repeated over the injector head. The aim of the work presented here is to model and analyze the re¦ll process of the components in the combustion chamber behind the rotating detonation. The simulation starts just after the passage of the detonation over the considered injection element. This simulation gives information on the way the injected propellants recreate the reactive mixture for the next detonation. In the ¦rst step, two-dimensional (2D) computations helped us to set up the methodology and to study the dynamic response of the fresh components injected. A comparison between 2D homogeneous and separate injections is provided. In the second step,


NUMERICAL INVESTIGATION OF AN UNSTEADY INJECTION ADAPTED TO THE CONTINUOUS DETONATION WAVE ROCKET ENGINE OPERATION
T. Gaillard, D. Davidenko, and F. Dupoirieux ONERA Chemin de la Hunière, BP 80100, Palaiseau CEDEX 91123, France Detonation applied to propulsion could result in a promising increase of the thermodynamic e©ciency of the engine cycle. Numerical simulations of the detonation propagating in the Continuous Detonation Wave Rocket Engine (CDWRE) are currently performed but still do not account for realistic injection process. The assumption of an ideal injected premix is generally chosen for convenience to obtain theoretical results. Comparison of the numerical results with experiments is di©cult because of the clear di¨erence of the injection con¦gurations. Some physical aspects of the separate injection of the components used in experiments are not clearly assessed. This study is included in a wider numerical project aimed at designing and optimizing a realistic CDWRE. The optimization process is presently focused on the injector. One element of the injection hole pattern is considered assuming that this element is periodically repeated over the injector head. The aim of the work presented here is to model and analyze the re¦ll process of the components in the combustion chamber behind the rotating detonation. The simulation starts just after the passage of the detonation over the considered injection element. This simulation gives information on the way the injected propellants recreate the reactive mixture for the next detonation. In the ¦rst step, two-dimensional (2D) computations helped us to set up the methodology and to study the dynamic response of the fresh components injected. A comparison between 2D homogeneous and separate injections is provided. In the second step, three-dimensional (3D) computations have been performed with a separate injection suitable for the CDWRE operation. Some performance parameters are evaluated such as mixing e©ciency or ¦lling of the domain.

INTRODUCTION
It is widely accepted that detonation engines represent an alternative to the conventional ones, whether it is applied to aircraft or rocket propulsion. Its bene¦cial e¨ect on the thermodynamic cycle has been theoretically proved by Zel£dovich [1] and Wintenberger and Shepherd [2]. The additional pressure gain provided by detonation makes it a more interesting combustion regime than classical de §agration.
Over the past decade, Rotating Detonation Engines (RDE) have aroused scientists£ interest in both numerical and experimental ¦elds. Since the ¦rst numerical simulations of a continuous detonation under 2D hypothesis [3], some important achievements have greatly helped to simulate RDE operation in three dimensions. Eude et al. [4] have succeeded in computing a 3D rotating detonation in an RDE chamber with a uniformly distributed injection using adaptive mesh re¦nement. Three-dimensional propagation of a periodic rotating detonation has been studied by Schwer and Kailasanath [5] over a row of several holes for a premix injection. They have observed the feedback pressure e¨ect on the continuous injection and the consequences on detonation stability. Recently, Stoddard and Gutmark [6] have tried to simulate a 3D rotating detonation in a centerbodiless RDE fed with a separate injection. Only one detonation has been propagated in the chamber, indicating the need of additional improvements. This numerical study focuses on analyzing an injector applicable to a rocket engine operating in the continuous detonation mode, which is also called CDWRE. In the CDWRE combustor, stable propagation of detonation waves is ensured by a layer of fresh mixture continuously formed close to the injector head. The present authors consider that detonation waves travel along the azimuth of an annular combustor. The feeding process involves gaseous injection of two propellants, H 2 and O 2 , at a stoichiometric ratio. As detonation speed exceeds 2 km/s, the re¦ll process is a critical phase of the CDWRE operation and must not exceed a few tens of microseconds. In this short amount of time, high uniformity of the fresh mixture is required as well as low losses of the injection total pressure to keep the theoretical bene¦t.
The propellants can be introduced separately or as a premix in the combustor. Premixed injection represents the ideal injection case for which the highest possible compression from the detonation can be obtained. Although it is massively used to simplify the computation of rotating detonation, it is hardly applicable to experiments. Detonation is likely to propagate back into the feeding system through the injection holes if their diameter is not inferior to a critical tube diameter dependent on the initial pressure of the fresh mixture [7]. Some experimental attempts have been made by Andrus et al. [8] to propagate two consecutive detonations in a premixed layer re¦lled as if it were in a RDE. The appearance of §ashbacks and blowo¨s con¦rms that it is not a safe operation to use a premix.
In case of transmission into the feeding collector and mixing device, a brutal energy release in a con¦ned volume can cause an important deterioration of the combustor, as reported by Thomas et al. [9]. To avoid the aforementioned risk, separate injection is commonly used in experiments. The advantage of the perfect mixing is now lost. Thus, the issue is to obtain fast and e©cient mixing with turbulent structures while limiting the total pressure losses to preserve the potential of detonation. Some experimental studies tend to show that mixing is essential for detonation stability but still little information is available to characterize mixing state in the fresh layer. What is known is that a de¦cit of 10% to 40% of the detonation speed from the theoretical ChapmanJouguet (CJ) velocity is generally observed [10,11]. This di¨erence may be due to inhomogeneities of composition and pressure losses. That is why, this numerical study is dedicated to the characterization and optimization of separate injection. Recent numerical works have been undertaken to propose new con¦gurations for mixing propellants. Stoddard et al. [12,13] tried di¨erent geometries for separate air and hydrogen injection tubes with a single detonation propagation over a linear series of 10 injection elements. They wanted to study the expansion of burnt gases and the ability to re¦ll the fresh mixture by propagating a real detonation. Such kind of simulation must be developed to ¦nd out more about unsteady process occurring during RDE operation.
The present study is included in a wider project aimed at computing an entire CDWRE chamber with a suitable feeding system. Besides mixing issues, we are also interested in evaluating the propulsive performance of a CDWRE. Future studies may involve full-scale engine design with a thrust nozzle. A speci¦c modeling approach for the injector might be needed to have such a simulation at a¨ordable cost. In a previous analysis [14], 3D Large Eddy Simulation (LES) was performed to compute the mixing of propellant jets and to evaluate H 2 /O 2 mixing e©ciency and pressure losses after a su©cient time to establish the §ow. Since then, we are progressing step by step to understand the mixing interactions at low scales and to choose the best way to inject the components. The focus is on understanding the mixing interactions at the scale of the injected jets. The turbulence integral scale is considered to be of the order of the tube diameters. Then, a single element simulation on a ¦ne mesh enables us to obtain a resolved scale about one tenth of the integral scale. Even if the simulated conditions did not correspond to the CDWRE chamber environment, it was necessary to choose the injection element design for the following study. This paper presents the latest new simulations in which unsteady re¦lling process after detonation passage is taken into account. Even if the injection blocking and jet development phases are reproduced, chemical reactions in the injected propellants due to partial mixing with the burnt products are not modeled. Detonation propagation is not directly simulated but the expansion of burnt gases past the detonation wave is modeled. This new methodology is presented hereafter and together with results for premixed and separate injection in 2D and 3D simulations.

MODELING STRATEGY FOR THE CONTINUOUS DETONATION WAVE ROCKET ENGINE INJECTOR
Before computing a complete  Two points are of particular interest when designing and analyzing an injection process: (i) the way reactants mix; and (ii) the total pressure losses involved. After optimization, the injection device will be studied under conditions of the 3D propagation of a rotating detonation.
The present study is based on a particular injector design. It is assumed that the CDWRE has an annular cylindrical combustor delimited by inner and outer walls and the injector face. The injector is designed as a repetition of one element both along the azimuth and radius. An example of such a geometry is depicted in Fig. 1.
From the whole injector, one single element is extracted (see Fig. 1) and becomes the reduced computational mixing domain shown in Fig. 2. Flow interactions between the adjacent elements can be taken into account with either periodic or symmetric boundary conditions (3 and 4) according to the adopted repetition principle. In this study, 3 and 4 are the periodic boundaries. The inlet conditions (1a and 1b) can be de¦ned by imposing a mass §ow rate or a total pressure together with a total temperature and injected gas composition. For the outlet 5, only pressure is prescribed. The boundary layer is not taken into account. Hence, the mesh is not re¦ned close to the tube walls and the injection plane. As a consequence, slip conditions are applied on tube walls (2a and 2b). This condition does not feature reality but helps to better respect the mass §ow rates of each propellant that leads to stoichiometric proportions, which is more convenient for the analysis.
Geometric design for the injection element is presented in Fig. 3 with all the parameters of interest as de¦ned for the best con¦guration in the previous study (referred to as 3a ′ ) [14]. This con¦guration is a semi-impingement of propellant jets with suitable distribution. It is our 3D case for this study.

NUMERICAL PROCEDURES
Computations are done with the CEDRE code [15,16], which is a multiphysics software developed at ONERA to solve energetic and propulsion related problems. Only the gaseous solver CHARME has been used to perform the simulations.
The compressible NavierStokes equations are solved for a nonreacting §ow on a regular mesh with a 50-micron cell size in order to capture a wide range of vortex scales. As the mixing process is strongly unsteady due to hydrodynamic instabilities, the LES approach is chosen to simulate the energetic turbulent scales. The subgrid viscosity is modeled by the Smagorinsky model for its simplicity and as a ¦rst hypothesis to account for the nonresolved turbulent scales. A ¦nite-volume method for general unstructured meshes is adopted for the spatial discretization. Second-order accuracy is obtained with a MUSCL (Monotonic Upstream Scheme for Conservation Laws) interpolation scheme coupled with the Van Leer slope limiter for the convective §uxes. The viscous §uxes are computed using an adapted second-order centered scheme. An implicit time integration is performed based on the Generalized Minimal RESidual (GMRES) method to solve the linearized equations system.
Hereafter, two di¨erent numerical procedures are proposed to study, on the one hand, an established §ow with prescribed Mach injection and outlet pressure; and on the other hand, a transitory phase corresponding to the re¦ll of the chamber linked to the expansion of the burnt gases produced after detonation passage over the injection element. While the ¦rst procedure has already been developed [14], the second one is explained and tested ¦rst on 2D cases and then applied to the 3D con¦guration.

Particular Procedure for Established Injection Regimes
In [14], three di¨erent jet interactions (sheared, impingement and semiimpingement) with two arrangements (periodic or symmetric repetition) along the x-axis were compared. The aim of this ¦rst analysis was to determine the best injection element based on mixing e©ciency and total pressure recovery after a su©cient time to establish the §ow. The present authors selected the semiimpingement con¦guration with periodic arrangement and this selected con¦guration is used for the calculations presented in sections 4 and 5. A mass §ow rate is imposed at the inlets in order to obtain a stoichiometric composition with the following §uxes related to the injection sections of each component: j m,O2 = 844 kg/(m 2 s) and j m,H2 = 213.5 kg/(m 2 s). The corresponding injection Mach is equal to 0.8 with a prescribed pressure of 0.25 MPa at the initial time and on the outlet boundary 5 (see Fig. 2). Total temperature is set to 300 K. Tube lengths are set to 4d H2 . Mixing domain length L is equal to 15d H2 . A ¦rst computation is made with a large time step of 1 µs with a ¦rst-order implicit Euler method to evacuate the initial conditions for a physical time of 1 ms. Then, a second-order implicit RungeKutta method is applied to simulate 3 to 4 residence times with a 0.01-microsecond time step. The §ow¦eld of interest is obtained with another 0.01-microsecond time step computation during a physical time of 100 µs.  from the very beginning of the process. Before starting studies on the reactant mixing in the reduced domain, 2D computations of a rotating detonation were performed with the CEDRE code to observe the physical and representative content of the unsteady re¦ll. Stoichiometric H 2 /O 2 premix was injected through the lower boundary of a rectangular periodic domain with a mass §ux of 100 kg/(m 2 s) in the chamber and an established §ow¦eld was obtained with a detonation speed D = 2670 m/s. An instantaneous temperature ¦eld from this simulation is presented in Fig. 4. The feeding slot corresponds to the area of negative y-value where chemical reactions are not activated. From the slot to the chamber, the depth of the domain cross section is increased by a coe©cient 10/3. As the detonation pressure at y = 0 is stronger than the injection pressure, a shockwave propagates in the feeding slot but it does not a¨ect the main physical phenomena. The higher burnt gases pressure leads to the stop of the injection for a few microseconds.
The evolution of the vertical Mach number M y and the hydrogen mass fraction Y H2 are extracted along the level y = 0.5 mm and shown in Fig. 5. The coordinate x in Fig. 4 has been transformed in the time variable in Fig. 5 according to: where x front is about 0.034 m. In Fig. 5, one can see that it takes 6 µs to reinject the right Y H2 value. This blocking time represents up to 30% of the re¦ll period, the later evaluated to 18.7 µs. When it is re¦lled, the Mach number of the fresh mixture continuously increases to reach a value higher than 1. This transition leads to some pressure losses in the propellants §ow before detonation burns the mixture. The re- Figure 5 Results from a 2D CDWRE case: vertical Mach number (1) and hydrogen mass fraction (2) pro¦les at the injection plane y = 0.5 mm during a re¦ll period ¦ll time is very short but as we inject ideal premix, we only have to take care about dilution of the reactant mixture by the burnt gases. Thanks to these intermediate results, we can now derive the developed methodology to compute the mixing and re¦ll process. Unlike in [12,13], we do not study multiple injection elements aligned side by side in a row but only one element. We cannot propagate a detonation over one single injection element with periodic boundary conditions because the re¦ll time would be too short. In the present calculations, we only intend to reproduce the physical effect of the burnt gases expansion. To determine the initial conditions of such a calculation, it is useful to examine the vertical pressure pro¦les after the detonation propagation (x < 34 mm) over the injection plane obtained from the aforementioned calculation (Fig. 6). These pro¦les account for pressure evolution across the wave structure composed of the rarefaction waves, the slip line and the shockwave. Only the rarefaction wave is important for the injection as the others never reach the injection plane.
With our methodology, we want to produce initial conditions and control the expansion process of the burnt gases that are local in x but su©ciently representative for our mixing study. The wave structures displayed in Fig. 6 resemble to the solutions of a Riemann problem extracted at particular times in an (yt) diagram. To model this Riemann problem, let us ¦rst assume that the detonation propagates with a plane front whose height is h (Fig. 7). The initial discontinuity occurs at the top of the detonation front at y(t 0 ) = h. The discontinuity divides the t 0 state in two homogeneous states, referred to as the right ¤rh¥ (y < h) and left ¤lh¥ (y > h) states. For convenience, the (yt) diagram is displayed with a representation adapted to the detonation front. Also, the t-axis is turned into a x-axis to be related to Fig. 4.
The lh state is de¦ned by the CJ conditions, representing the burnt gases state just behind the detonation front. It is calculated by extracting the mean properties of the fresh gases (m) just in front of the detonation wave from the CDWRE case (see Fig. 4), which are the following: The rh state de¦nes some arti¦cial state to control the expansion process of burnt gases. To obtain this state, it is assumed that the CJ burnt gases from the previous detonation have undergone an isentropic expansion, from the CJ state with vertical Mach (M CJ ) equal to 0 to a prescribed vertical Mach (M rh ) in the rh state. The mass fractions are constant across the discontinuity since in the present study, chemical reactions have been neglected.
The value of M rh in §uences on the opening of the rarefaction wave. For example, for a supersonic M rh , the area that in §uences the injection plane is   Fig. 7. Indeed, the sonic Mach limit is located at y = h. All the waves above this limit will not reach the injection plane.
From the initial discontinuity between the CJ and the rh states, the exact Riemann problem has been solved with an in-house code. The self-similar solution can be obtained from t 0 (x 0 ) to t 1 (x 1 ) where t 1 is the instant at which the expansion wave reaches the injection plane. After this time, the solution is not considered because the re §ection of the expansion fan is out of the scope of the Riemann problem. We obtain a multitude of possible pro¦les between x 0 and x 1 to be tested as initial conditions. The solutions obtained in x = x 0 and x = x 1 are given in Fig. 8 for M rh = 1.5.
Note that with a high vertical velocity in the chamber, the expansion process leads to an expansion pressure that is lower than P rh . Hence, the evolution of P , T , and M is reversed from y/h = 3 to 3.5, compared to what the pro¦les in Fig. 6 suggest. Nevertheless, with the proposed technique, a suitable P (t) evolution of the expansion process has been obtained along y in the zone of in §uence (y < h). Now, let us make the following assumptions to apply these pro¦les as initial conditions for the injection and mixing simulation: (1) the injection element is su©ciently small to neglect §ow parameter variation along x and z; and (2) x-and z-wise velocities are initially set to 0. By analyzing numerical results, it was found that the x-wise velocity is important right behind the detonation but then it rapidly decreases and is relatively low during the reinjection phase.

Two-dimensional test cases
For the ¦rst time, the methodology is tested on simple 2D cases for premixed and separate injection. The corresponding 2D domains are presented in Fig. 9. For the separate injection, the widths of the H 2 and O 2 injection pipes are equal to d/3 and 2d/3, respectively, with d = 1 mm. The mesh is re¦ned in area 1 in Fig. 9 (−15d < y < 15d) to obtain cells of 50-micron side, whereas it is coarsened in areas 2. The upper initial conditions (y > 0) are determined by the present authors£ previous methodology whereas the lower initial conditions (y < 0) correspond to a uniform §ow separated by a discontinuity at y = 0. By doing this, it is assumed that before the detonation arrival, the injection was completely established; then, the fresh mixture above the injection plane was instantly burnt by detonation without perturbing the initial state in the pipes. Two validation cases, referred to as 1.1 and 1.2, are introduced in Table 1 taking the solutions of the Riemann problem (see Fig. 8) in x = x 0 and x = x 1 as the initial upper conditions, with M rh = 1.5. They correspond to premixed injection with an initial Mach M ini = 1 in the injection pipes. The inlet conditions set at y = −50d for premixed injection are de¦ned by total pressure (P t,inj ) and temperature (T t,inj ). From a given mass §ux j m = 100 kg/(m 2 s) related to the chamber cross section, the total injection pressure was derived through the following relationship: where W is the gas molar mass in kg/mol; A %,inj is the injection area ratio equal to 0.2; and R = 8.314 J/mol/K and γ = 1.4 are the fresh mixture parameters. The initial static conditions in the tubes (P s,ini and T s,ini ) are calculated to comply with the initial injection Mach number M ini = 1:    The physical e¨ect of the burnt gases expansion is better featured in case 1.2. Despite a slower pressure decrease, it can be considered that the time evolution obtained with the conditions of case 1.2 correctly reproduces the evolution given by the CDWRE calculation. The time at which low-temperature propellants replace the burnt gases is quite satisfactory in case 1.2, whereas this time is much too long in case 1.1. Therefore, the initial conditions deduced from the solution of the Riemann problem at x = x 1 are considered appropriate for the following calculations. Now, let us analyze the e¨ect of M rh on the re¦ll process. Cases 2.1 to 2.4 ( Table 2) are the premix injection cases with subsonic M ini and M rh varying from 0.8 to 2.
As expected for cases 2.1 to 2.4, the stabilized injection Mach (M inj ) after burnt gases expansion during 100 µs is di¨erent from the initial prescribed M ini because the ¦nal pressure level (P ch ) in the chamber is conditioned by M rh . These dependences are shown in Table 3. With M rh greater than 1.5 (cases 2.3 and 2.4), the acceleration induced by the burnt gases movement in the chamber produces a transition from the subsonic injection regime (Mach of 0.8) to a sonic injection regime. Conversely, M rh = 0.8 (case 2.1) restrains the re¦ll process to M inj = 0.3.
The actual injected mass §ow rate in cases 2.1 to 2.4 is then di¨erent from the initial value since the injected §ow adapts to the stabilized chamber conditions. The nondimensional mass of the fresh mixture m % injected in the chamber is evaluated as follows: where S det = 5dh is the imposed detonation front area; S ch is the domain area above the injection plane; and ρ m is the mean density of the fresh mixture  Fig. 11 for the four cases. When m % reaches 1, it means that the initial mass of fresh gases is fully restored in the chamber. Fig.12 to show the di¨erent injection phases. Figure 12a is the initial pressure ¦eld. In Fig. 12b, the burnt gases have still a high pressure leading to a push back e¨ect on the fresh mixture in the feeding pipes. Pressure in the chamber continues decreasing and the re¦ll phase begins (Fig. 12c). Finally, the jet develops in the chamber with typical vortex structures shown in Fig. 12d. Cases of separate injection 3.1 to 3.4 are de¦ned in Table 4 similarly to cases 2.1 to 2.4 but for the computational domain shown in Fig. 9b. Compared to the previous cases, the re¦ll process is di¨erent because the separate injection uses an H 2 mass §ux j m,H2 = 11.11 kg/(m 2 s) and an O 2 mass §ux j m,O2 = 88.89 kg/(m 2 s). Both §uxes are related to the chamber cross section given that A %,inj,H2 = 1/15 and A %,inj,O2 = 2/15. For both propellants, the total injection pressure is now P t,inj = 0.274 MPa to keep a total mass §ux of 100 kg/(m 2 s) with M ini = 0.8. This di¨erence can be seen in Fig. 13 where m % is plotted vs. time. Contrary to Fig. 11, m % does not reach 1 at t = 50 µs whatever the case.

Case 2.2 is illustrated in
What is more interesting in the case of separate injection is to analyze the equivalence ratio of injected propellants vs. time. In Fig. 13b, one can see that for case 3.2, it takes 50 µs to inject the required global ER. It shows the delay that exists between the re¦ll of H 2 and O 2 . The time necessary to obtain the aimed value of global ER is an important issue in the case of separate injection.
With respect to the CDWRE case, the di¨erent phases of the re¦ll process are well reproduced. The time dedicated to the blocking phase represents 20% of the re¦ll period (de¦ned by the instant when m % reaches 1) and lasts about 10 µs for the premix and separate cases with maximum mass §ow rates. Compared to the CDWRE case, the process is longer and this will certainly result in a global e¨ect on chamber operation.

Preparation of three-dimensional unsteady computations
To adapt the 2D procedure for the 3D computations, some modi¦cations have to be done to limit the computational cost. Until now, the feeding pipes were excessively long to prevent the initial shock wave from re §ecting at the tube inlets and perturb the injection. Nonre §ective boundary conditions have been tested and validated to absorb the shockwave. The tube length can be reduced to 10 diameters. This length is su©cient to compute the burnt gases penetration in the feeding tubes.
Another change consists in unifying the burnt gases composed of six species in one equivalent gas referred to as BG to reduce the number of resolved equations. Let us evaluate the thermodynamic properties of BG: molar mass W ; enthalpy of formation h f ; heat capacity at constant pressure c p ; dynamic viscosity µ; Prandtl and Schmidt numbers. Molar mass W is calculated by harmonic mass weighted average; and h f and c p,i (coe©cients of c p ) are calculated by arithmetic molar weighted average. Parameters of the Sutherland law for µ BG are obtained with a least squares method to ¦t the viscosity of the burnt gases evaluated from the Wilke formula. The Prandtl and Schmidt numbers are assumed to be constant for the computations. The Prandtl number is evaluated from the thermal conductivity pro¦le and chosen for a temperature of 2000 K. A constant Schmidt number is evaluated to correspond to a stoichiometric mixture of H 2 and O 2 . These assumptions for the Prandtl and Schmidt numbers might be strong but turbulent transport is supposed to be much stronger than the molecular one in the present simulations.
Finally, to improve results visualization, let us introduce a new variable Z, by analogy with the well known mixture fraction, based on molar fractions of the three species (H 2 , O 2 , and BG): The Z parameter is designed for proper visualization of computational results for the mixing zone. Using mass or mole fractions of the gases is not convenient. Because of very di¨erent molar masses, the graphical representation gives the impression that one of the components is dominating. In a stoichiometric mixture without burnt gases, mole fractions of H 2 and O 2 are related as X H2 /X O2 = 2. From formulas (2) and (3), one obtains for the same mixture Z H2 = Z O2 = 0.5 that is more convenient for graphical representation of the mixing zone. Now, consider that the stoichiometric mixture is burnt by 50%. By using the global reaction expression 2H 2 + O 2 = 2H 2 O and substituting BG for H 2 O, one ¦nds Z H2 = Z O2 = 0.25 and Z BG = 0.5 from formula (4). The same values for Z H2 , Z O2 , and Z BG are obtained in the case of dilution by the burnt gases by 40%. Therefore, Z BG represents a convenient parameter to estimate the degree of dilution or burning of the fresh mixture. We can go further by introducing three more parameters de¦ned by formulas: where Z H2,ex and Z O2,ex represent excess of the respective components from the stoichiometric proportion whereas Z st is the fraction of stoichiometric mixture. By comparing ¦elds for Z H2,ex , Z O2,ex , and Z st , one can easily differentiate the zones of good and poor mixing. In case of a §ow of pure components, the equivalence ratio is a convenient parameter for visualization. However, it is less useful if the mixture quality is a¨ected by dilution. In the latter case, Z st will help to identify zones where the stoichiometric mixture is nondiluted.

RESULTS FOR THREE-DIMENSIONAL ESTABLISHED INJECTION CASES
For this computation, the inlet mass §uxes related to the injection sections are the following: j m,O2 = 844 kg/(m 2 s); j m,H2 = 213.5 kg/(m 2 s); and T t,inj = 300 K. The established injection results give a reference §ow¦eld representing what mixing is expected to be after a time long enough for the vortex structures to develop. The Q-criterion is used to visualize the instantaneous eddies in the turbulent §ow (as shown in Fig. 14). In Fig. 14a, isosurface of the Q-criterion (10 8 ) is colored by the value of ER to give an idea about local distribution of ER when the stoichiometric proportions of the propellants are globally respected. The ER color scale is logarithmic from 0.1 to 10 whereas the overall ER variation is from 0 in pure O 2 to in¦nity in pure H 2 . One can easily deduce from Fig. 14a that the components are not well-mixed in the lower half of the domain, as intermittent zones of unmixed propellants are observed whereas in the upper half, the mixture is almost homogeneous. This representation does not consider the dilution of the fresh mixture by the burnt gases. As the established results are compared with the unsteady ones, Z st is more appropriate as it takes into account the burnt gases fraction; hence, it is possible to distinct the zones ¦lled with well-mixed propellants from the zones with strong dilution by burnt gases.
One can also determine the zones where H 2 or O 2 are in excess with the parameters Z H2,ex or Z O2,ex . In Fig. 14, it is demonstrated how the di¨erent parameters characterize the mixture quality. From a quantitative point of view, it is worth analyzing the evolution of mixing according to the nondimensional vertical coordinate y adim = y/d H2 . Instantaneous mixing e©ciency can be evaluated from two formulas, with either mass or mass- §ow weighting: η mix,ρ = Sy ρY H2 dS/max(ϕ, 1) Sy ρY H2 dS/ min(ϕ, 1) ; (5) Figure 15 Comparison of time-averaged mixing e©ciencies evaluated according to Eqs. (5) (1) and (6) (2) for 3D injection where ϕ = W O2 Y H2 /(2W H2 Y O2 ) is the local ER, with W O2 and W H2 being the molar mass of O 2 and H 2 , respectively. η mix is varying in the range [0, 1], the value of 1 being obtained when ϕ is equal to 1 everywhere in the considered S y section, i. e., when the mixing is perfect in this section. From the instantaneous pro¦les, time-averaged e©ciencies η mix,ρ and η mix, ' m are calculated by simple arithmetic average and displayed in Fig. 15. The two e©ciencies show similar behavior after y adim = 6. Close to the injection plane, some stagnation and recirculation zones exist. Propellants in these low-velocity zones are partially premixed resulting in η mix,ρ > 0 at y = 0. The mass- §ow-weighted e©ciency (see Eq. (6)) does not take into account stagnant zones but only the zones of high velocity that are fully unmixed at y = 0. That is why η mix, ' m is 0 at the beginning and rapidly increases thanks to the jet interactions. Mixing by jet interactions becomes e¨ective from y adim = 4 only where η mix,ρ joins the trend of η mix, ' m . Equation (6) is well adapted to established computations whereas it is not well suited for transitory re¦ll because mass- §ow rate is not constant and stabilized. Equation (5) seems better ¦tted for transitory re¦ll simulations. Total pressure losses are not considered in this present study. In [14], it was shown that a strong correlation exists between mixing e©ciency and total pressure recovery based on mass- §ow weighted average for di¨erent injector con¦gurations.

RESULTS FOR THREE-DIMENSIONAL TRANSITORY INJECTION
The aim of these simulations is to identify the temporal evolution of the injection and mixing from the instant just after the passage of detonation over the injection element to the achievement of a mixing state comparable with the established situation. For the initial conditions, it was adopted: M ini = 0.8 and M rh = 1. The inlets are de¦ned by the following total conditions: P t,inj = 0.381 MPa and T t,inj = 300 K. Same mass-weighted mixing e©ciency calculations can be done from the instantaneous §ow¦eld. Figure 16 shows instantaneous pro¦les of mixing e©ciency extracted from the transitory computation every 20 µs. They are compared with the established pro¦le as a reference state. The re¦ll process is supposed to converge toward the stabilized §ow¦eld but at 100 µs, it has still not completely reached the ¦nal state. The interaction between burnt and fresh gases can be analyzed using parameters Z st and Z BG . In Fig. 17, upper row, the progression of mixing in the computational domain is illustrated by the isosurfaces of the Q-criterion, colored by Z st , every 20 µs. The Q-criterion isosurfaces are also shown in Fig. 17, lower row, but colored by Z BG . In the ¦rst 20 µs (Fig. 17a), burnt gases are pushed back from the feeding tubes and propellants start to be reinjected in the domain. There are still much burnt gases in the mixing zone (Fig. 17a, lower row) and almost no good mixture (Fig. 17a, upper row).
In the next 20 µs (Fig. 17b), the turbulent §ow ¦lls the domain and interactions with neighboring injection elements appear as the turbulent structures cross the periodic boundaries. However, mixing is developing in the middle of the domain but is not yet high enough because the §ow¦eld is intermittent.  Some burnt gases remain near the injector face. Between 60 and 100 µs, turbulent structures become ¦ner and improve the mixing of the propellants. Burnt gases have been rejected from the lower half of the domain but remain in the upper half, as mixing e©ciency continues to improve in this part of the domain.
At the end (Fig. 17e), the ¦nal mixing state is similar to the one obtained for the established case with a di¨erence from 7% to 30%, depending on the y adim coordinate. This di¨erence seems to be due only to the imperfect mixing as burnt gases have nearly left the domain. The mean Z BG level in Fig. 17e (lower row) is 4%.
These results provide an evaluation of the propellant dilution in a nonreacting context. Dilution is related to the mixing of H 2 or O 2 or both with the surrounding burnt gases and is responsible of two important loss factors. First, a diluted mixture is less energetic; hence, the detonation strength in such a mixture is reduced. Second, diluted mixture can self-ignite before the next detonation passage, resulting in a loss of available mass for detonation. Both factors conduct to a loss of pressure gain in the chamber and by consequence to a loss of RDE performance. With the present 3D methodology, the authors propose the ¦rst step in understanding the transient phenomena of propellant injection and mixing occurring between two detonations in a nonreacting environment. A numerical study accounting for chemical reactions and the coupling between the injector response and mixture quality, on the one hand, and the detonation strength and burnt gases expansion, on the other hand, is under way.

CONCLUDING REMARKS
In this paper, numerical studies were presented in the framework of a global project to simulate CDWRE operation with a realistic injector for separate injection of the propellants. Two numerical procedures were presented. The ¦rst procedure is proposed to compute an established injection of the two propellant jets, giving the best mixing state that is expected to be after a long time. The second procedure is aimed at computing the transitory re¦ll phase from the moment the detonation has passed over the injection element. Test cases associated with this procedure have shown that the expansion e¨ect of the burnt gases on the fresh mixture is well reproduced. The results analysis points out that the stabilized conditions are greatly in §uenced by the Mach number spec-i¦ed for the right state of the Riemann problem. The Mach number in the injection pipe deviates from its initial value during the injection process, thus leading to a mass §ow rate change during injection. A 3D transitory simulation was ¦nally performed to study the re¦ll process up to the almost established §ow.