Numerical Simulation of Reactive Flows in Overexpanded Supersonic Nozzle with Film Cooling

Reignition phenomena occurring in a supersonic nozzle flow may present a crucial safety issue for rocket propulsion systems. These phenomena concern mainly rocket engines which use H 2 gas (GH 2 ) in the film cooling device, particularly when the nozzle operates under over expanded flow conditions at sea level or at low altitudes. Consequently, the induced wall thermal loads can lead to the nozzle geometry alteration, which in turn, leads to the appearance of strong side loads that may be detrimental to the rocket engine structural integrity. It is therefore necessary to understand both aerodynamic and chemical mechanisms that are at the origin of these processes. This paper is a numerical contribution which reports results from CFD analysis carried out for supersonic reactive flows in a planar nozzle cooled with GH 2 film. Like the experimental observations, CFD simulations showed their ability to highlight these phenomena for the same nozzle flow conditions. Induced thermal load are also analyzed in terms of cooling efficiency and the results already give an idea on their magnitude. It was also shown that slightly increasing the film injection pressure can avoid the reignition phenomena by moving the separation shock towards the nozzle exit section.


Introduction
One of the major challenges that the aerospace industry continues to face is the continued increase in launchers payload. For example, the latest version of the Vulcain-II engine of the European Ariane 5 ECA launcher is able to put into geostationary orbit a payload of about 10 tons. The goal is however to reach 12 tons in the near future. There are two challenges here: the increase in payload and the performance consolidation in terms of reliability.
These challenges promote development of nozzles with higher performances, which are substantially achieved by increasing the nozzle expansion area ratio or by developing new innovative nozzle concepts. The rocket engine nozzles, with high expansion area ratio, are generally optimized for operating at high altitudes. At sea level and at low altitudes, the nozzle operates in overexpanded flow conditions; that is, the ambient pressure is higher than the nozzle exit pressure. The resulting adaptation shock may lead to flow separation, unsteadiness, and shock interaction. The ensuing side loads may be detrimental for both nozzle and other engine components. In addition, these nozzles are designed to expand and accelerate combustion gases at high temperature. To avoid thermal loads, designers adopt several nozzle cooling methods. The most effective one uses the film cooling technique. For example, for Vulcain-II rocket engine, the cooling system is designed in two parts: a dump cooling for the first expansion part of the nozzle and GH 2 film cooling for the second part [1,2].
The present study is a contribution to the works initiated by the CNES during the last ten years in the field of innovative nozzle researches. The main goal is to develop new supersonic nozzle concepts for cryogenic rocket engines. One issue that attracted the interest of scientists is that related to the risk of reignition that can occur in the main flow, resulting from combustion of GH 2 used as film cooling. This phenomenon may happen when the GH 2 film is mixed with the air engulfed into the separated region along the mixing shear layer at high  temperature. In this occurrence, the film cooling efficiency is altered and the induced thermal loads may be critical for these nozzles. To understand these phenomena, a campaign of experimental measurements was conducted at the cryogenic LH 2 -LO 2 ONERA/Mascotte test bench facility [3]. The results yet obtained for low oxidizer/fuel ratios ( / ), on the basis of PLIF and Kulite pressure measurements, already constitute a very useful database for understanding the involved physicochemical phenomena. The first series of measurements were performed under / ratios of up to 3 and combustion chamber pressure up to 40 bars. The first results obtained during this measurements campaign highlighted a reignition phenomenon occurring in the nozzle main flow, inside the separation area.
Given the complexity of these phenomena, numerical simulations performed by different computing codes can give a better interpretation of the experimental results. The present numerical work is a further contribution in this regard and consists in simulating the turbulent reactive flows in the same operating conditions as those performed experimentally.
This work was carried out within the framework and with the support of the French research group ATAC (Aérodynamique des Tuyères et Arrières Corps (Nozzles and Afterbodies Aerodynamics)).

Objectives
The aim of this paper is to perform 2D RANS turbulent reacting flows simulations, for three test cases carried out experimentally in the ONERA/Mascotte test bench. For this purpose, two chemical kinetic schemes and two turbulence models were tested to investigate their relevance in the reignition process inside the nozzle.
Basically, the test bench consists of a subscale planar nozzle connected to a cylindrical combustion chamber fed with a cryogenic mixture LO x -LH 2 (Figure 1(a)). This hardware version has been developed especially to investigate the flow separation in overexpanded regimes. The combustion chamber operates under controlled total pressure and oxidizer/fuel ratios ( / ). The main geometrical characteristics of the nozzle are given in Figure 1(b). Pure hydrogen gas (GH 2 ), serving as film cooling, is injected tangentially into the upper wall of the nozzle through a slightly supersonic injector. Five Kulite pressure transducers are taped along this wall. A schematic longitudinal cross section of the nozzle is shown in Figure 1(b).
The first results [3] obtained by PLIF visualizations performed in this test bench show an increase in concentrations of OH radical near the upper wall of the nozzle's divergent, in the separation zone. This suggests a reactivation of the combustion inside this area ( Figure 2). During these International Journal of Aerospace Engineering 3 experiments, it was also revealed that there is unsteadiness in the radical OH emission inside the separation zone. The authors related this phenomenon to the pressure fluctuations observed in the combustion chamber. An example of PLIF images obtained by this technique is depicted in Figure 2 [3]. However, no detailed explanations have been given so far about the occurrence of this phenomenon.
Parallel to these experiments, steady RANS calculations have been carried out at ONERA/DEFA using the CEDRE Code. In these calculations the reduced chemical kinetic model of Eklund et al. [4] was used for chemical rates. However no process of reignition in the flow has been pointed out from the numerical results [3].
From these experimental and numerical conclusions further calculations using new chemical kinetics based on the well-known Evans and Schexnayder model [5] were recommended. This last model was slightly modified by adding two additional reactions, without excessively penalizing the computation time. Furthermore, two RANS turbulence models have also been tested in these calculations. The main objective of this study was to investigate relevance of chemical kinetics and turbulence models, in the context of reignition phenomenon.

Numerical Code and Equation Formulation
3.1. Numerical Code. The calculations presented in this paper were performed using the FASTRAN code. This code was specifically designed for compressible flow studies at high Mach numbers and based on solving the multispecies Reynolds-Averaged Navier-Stokes (RANS) equations with a finite volume formulation. The code offers two upwind differencing schemes with a variety of high order slope limiters to calculate the convective terms in the transport equations. Both explicit and implicit iterative and noniterative time integration schemes are available for steady-state and time accurate flow simulations. The convective fluxes are calculated by means of either flux difference splitting scheme (Roe) or flux vector splitting scheme (van Leer). The code also offers a choice of several turbulence models ( -, -, --SST-Menter, Spalart Allmaras, and Baldwin Lomax). The following sections recall the main fluid physical modelling and the flow field numerical formulations methodology used in this code for solving reactive Navier-Stokes equations.

Thermodynamics Gas
Properties. The thermal equation of state for a mixing or reacting gas is given by Dalton's Law of partial pressures such that where is the static pressure, is the species or mixture density, is the universal gas constant, is the static temperature, and is the molecular weight of the species .
The caloric equation of state relates the total energy to other gas dynamic variables and gas properties. For calorically perfect gas where is the total energy per unit volume and , V , , and are the gas density, the heat capacity of the gas mixture at constant volume, the static temperature, and the mass averaged velocity, respectively.
The form of the caloric equation of state for a mixing or reacting gas depends on the database used for describing the molecular properties of the species. Two databases are available: the first is based on molecular (or spectroscopic) data for chemical species and the second is based on fifthorder polynomial curve fits for each chemical species [6]. Using molecular properties, (2) can be written as where is the species density, V,tr, is the translationalrotational heat capacity for species at constant volume, Δℎ 0 , , is the heat of formation at reference temperature and pressure for species , and V is the molecular vibrational energy per volume.
Using of polynomial curve fits for properties gives two forms of caloric equation of state.
For thermal equilibrium For thermal nonequilibrium = ∑ [ V,tr, + Δℎ 0 , , − ( V,tr, − )] where ℎ is the sensible enthalpy per unit mass for species defined as where is calculated from fifth-order polynomial curve fits, for each chemical species, from Gordon database [6].
The viscosity of each fluid species , in the case of reactive flows, is calculated by the relationship of Bird et al. [7] = 2.66693⋅10 −6 √ , where is the characteristic molecular diameter and Ω is the viscosity collision integral. The characteristic molecular diameter is based on Lennard-Jones potentials [8]. For the mixture, the viscosity is calculated by the semiempirical relationship of Wilke [9] = ∑ Φ , where is the mole fraction of specie and Φ is given by For the diffusivity vector , the mass diffusivity can be represented by either Fick's law [10] or a binary diffusion model [11]. Consider where is the diffusion coefficient and is the species mass fraction. is given by where Sc is the Schmidt number.

Chemical Production.
The reactive flow calculation is obtained by solving the flow conservation equations in which one integrates a source term expressing the mixture chemical composition variation resulting from chemical reactions.
In the approach used for reacting flows, the general finite rate reaction is written as where ] and ] are the stoichiometric coefficients of the reaction and represents an arbitrary molecule in the reaction. According to Kuo [11], the source term for species is given by where is the coefficient of efficiency of the third body for the reaction , is the species concentration, and and are forward and backward reaction rates of a reaction , respectively. The concentration powers and are identical to ] and ] , respectively, for most applications, particularly for chemical kinetic reaction governed by Arrhenius rates of reaction where , , and / must be specified for each reaction under investigation.

Flow Field Numerical Method.
Basically, the conservation equations with appropriate closure models are expressed in vector form as In this expression, and represent the convective and diffusive fluxes, respectively, such as ,̇=  ] . (16) In this system of equations, , V, and are the velocity components ( = 0 for 2D calculations), is the mixture density, and is the shear stress tensor. and are, respectively, the species density and species velocity. One can note that, for calorically perfect gas, only one "species" is tracked such that = 1, = , and = 0. Note that, for thermal equilibrium calculation, all terms relating the contribution of vibrational internal energy ( int , int , V , . . .) are no longer required.

Chemical Kinetic Models
To investigate the pertinence of the chemical reactions associated with the reacting mixture issued from the combustion chamber, it was expedient to test the most suitable kinetic schemes for reactive H 2 -O 2 flow. Two kinetic schemes were selected for this study: the modified Evans-Schexnayder model and Eklund's kinetic model, commonly used by ONERA and CNES.
International Journal of Aerospace Engineering 5  [12]. In this system, N 2 operates as the third body and does not dissociate. This mechanism has been widely used for simulation in supersonic and hypersonic flows, particularly in the case of combustion initiation around obstacles or in scramjets [13][14][15]. This model was proved to be less expensive in terms of computation time; however it presents weakness in modeling the self-ignition delay (induction time) and in estimating the reaction heat release. This is mainly due to the absence of hydroperoxyl radical (HO 2 ) in this scheme. Indeed, studies have shown that fast three body recombination reactions, involving the radical HO 2 , have been identified as major contributor in the heat release process during the combustion of hydrogen with air [16].
To overcome this deficit, two reactions taken from the model of Rogers and Chinitz [14], involving this radical, have been added to the original Evans model, adding no substantial computation time Another insufficiency attributed to this model is its autoignition delay, which is relatively long, especially for reactions at low temperature (≈1000 K). This problem is related to the absence of hydrogen peroxide (H 2 O 2 ) in the model. Adding more reactions involving this species to correct this deficiency substantially complicates the model. An alternative solution would be to increase the production rate of Reaction 7 of the original model Indeed, this reaction has been identified as important in the case of inflammation at low temperatures [14]. Initially, the forward rate equation for this reaction is expressed as [5] = 2.2 ⋅ 10 14 exp (− 8455 ) .
This value was obtained with an accuracy of 50% for a temperature range of 300 to 2000 K. By multiplying the coefficient , by 3, the rate of hydroperoxyl radical production OH is increased which leads to reduction in the ignition delay [12]. The corresponding ignition delay becomes comparable to those obtained by more complex chemical kinetic models with more reactions. Finally, the modified Evans-Schexnayder kinetic model, with ten chemical reactions, is given in Table 1. Figure 3 depicts the results obtained for validation of this kinetic model, in the case of H 2 -O 2 combustion. The results are presented in terms of pressure and temperature rise, H 2 consumption, and OH and H 2 O formation, from onedimensional combustion simulation. In Figures 3(a) and 3(b), the initial and the modified Evans-Schexnayder model results are compared to those from more complex kinetic schemes of Rogers and Chinitz [17], Drummond [18], and Bitker and Scullin [19], respectively. The results clearly highlight the relevance of the added specific reactions on the ignition time delay [12].

Eklund
Model. This reaction simplified scheme, proposed by Eklund et al. [4] and implemented on both CEDRE and CPS codes, has been widely used by ONERA and CNES [3,20] for nozzle reactive flow studies. This scheme consists of 7 reversible reactions and 6 chemical species [O 2 , H 2 , OH, H 2 O, O, and H]. As can be seen in Table 2, this scheme does not involve any third body reaction which can present an advantage in terms of computing time.

Initial Conditions and Implementation of Calculations.
The calculations were performed over a computational domain which includes the nozzle, the injector, and the outside experimental environment. The calculation was initiated at the nozzle inlet using the same initial data as the experimental test conditions described below. No-slip conditions along the nozzle walls were assumed. For the outlet condition, 1 bar fixed pressure is applied at the downstream exit section. Adiabatic no-slip conditions are imposed for the rest of the surrounding block boundaries.
Combustion temperature and species mass fractions for the cryogenic LH 2 -LO 2 combustion products, at desired operating pressure chamber and / ratio, are obtained from separate calculations using the CEA thermochemical code [21].
In order to perform 2D CFD simulations of the nozzle's flow field, finite volume grids have been constructed using algebraic grid generator software. Multiblock structured grids have been used in this calculation. The computational domain includes the convergent-divergent parts of the nozzle, the injector, the zone downstream, and the area located on both sides of the nozzle. The dimensions of the downstream block are more than 70ℎ in both and directions, ℎ = 0.028 m being the nozzle height at its exit section. A part of the meshing used in the computational domain, representing the nozzle and part of the outer domain, is illustrated in Figure 4(a). The right part of Figure 4(b) depicts the mesh system used for the injection slot and the nearby region. It also depicts the mesh refinement near the wall. To achieve the CFD calculation, up to 40000 cells were required. The grid refinement near the wall, in the direction, such as the first cell-height leads to + = 1, where + is the dimensionless distance defined as with , , and ] being the shear stress at the wall, the density, and the kinematic viscosity, respectively. An average of 25 cells in the normal direction of the nozzle wall was used to resolve the boundary layer thickness. Furthermore, a preliminary study of the mesh refinement sensitivity was beforehand performed before achieving the final calculations. The fluid is considered to be an ideal gas mixture in thermal equilibrium. The mixture heat capacity at constant pressure ( ) is calculated as a function of temperature from the calorific capacity of each species, expressed in polynomial form. The Schmidt number Sc and the turbulent Prandtl number are set equal to 0.9 for all calculations.
Boundary conditions for turbulence model are also needed. For the -turbulence model, the turbulent kinetic energy and the turbulent dissipation rate required at the inlets for both nozzle and injector are calculated on the basis of an estimated turbulence intensity of 5% for this nozzle [3].
Moreover, the Spalart Allmaras turbulence model requires the kinematic turbulence viscosity ] of the mixture. For Reynolds numbers Re > 10 5 , which is the case here, a value of the turbulent viscosity ratio, namely, the ratio of turbulent viscosity to laminar viscosity ( = ] /]), up to 100 is a reasonable estimate.
In this study, the fluxes are evaluated at each time step using the Roe scheme associated with the second-order spatial accuracy MINMOD flux limiter. The time integration is performed by a fully implicit scheme with a local time step based on CFL number ramped from 0.1 to 0.5 within the first 1000 iteration cycles. Convergence to 10 −6 residuals was usually attained after 20,000 to 25,000 iterations.
Finally the code is based on the finite rate mechanism in calculating the chemical composition. Each chemical reaction obeys the law of mass action; hence the backward rate coefficient is calculated on the basis of the equilibrium constant for each reversible reaction.
The results obtained by the numerical simulations are presented in terms of flow Mach number field, temperature, and OH mass fraction fields. In addition to the PLIF pictures, five experimental wall pressure values are available. These pressures are measured along the cooled upper nozzle wall and will be compared to the calculated values in the same plots. Each test case is simulated using two chemical kinetic schemes, namely, "the modified Evans-Schexnaiyder model" and "the Eklund model. " Two turbulence models are also considered in this calculation, the standard two-equationand the one-equation Spalart-Allmaras model.
The operating conditions of the tested cases and the initial data for the numerical simulations inputs are summarized in Table 3.

Aerodynamic Field.
To understand the flow topology, we plotted the aerodynamic flow fields in terms of Mach number of the different tested cases. As depicted in Figure 5, for each test case, both turbulence models used in the computation predict a similar flow pattern. Indeed, in such nozzle regimes (overexpanded conditions), a separation shock is created inside the nozzle extension to adapt the flow to the outside pressure. This configuration is known as the Free Separation Shock (FSS). As sketched in Figure 5(a), in all configurations, the separation shock impacts the opposite wall and leads the boundary layer to separate similarly as in the so-called Shock Wave Boundary Layer Interaction (SWBLI). The other point the flow topology reveals is that relating to the shock position. As expected, this position depends on both nozzle pressure ratio NPR, for a given turbulence model, and turbulence model, for a given NPR, as illustrated in the comparisons given in Figure 5. This point will be emphasized with more details in wall pressure discussion below.
Moreover, it should be noted that an increase in the GH 2 injection pressure, used in the film cooling, as in the calculated Case 4 pushes back the separation shock downstream the injection slot towards the nozzle exit section. This result is depicted in Figure 6 in terms of calculated pressure profiles for Case 1 (NPR = 25.9, 0 = 3 bars) and for Case 4 (NPR = 25.9, 0 = 4.3 bars). As shown, the relative displacement of shock, determined with respect to the length of the nozzle part concerned with the fim cooling (Δ / ), is about 36.3%. A similar trend has been reported in previous x (m) Rogers     works [23], where similar film injection is used to control the flow separation in overexpanded supersonic nozzles.

Wall Pressure.
The calculated pressure profiles are compared to the experimental pressure measured by means of Kulite transducers. Figure 7 depicts the pressure distribution along the upper wall of the nozzle's divergent, downstream of the injection slot.
The results in Figures 7(a) and 7(b) show that the turbulence models used in this work were able to predict with good accuracy the pressure rise in the separation zone, for Test Case 1 (NPR = 25.9) and Test Case 2 (NPR = 22.2). The numerical results also show that the reignition phenomenon significantly affects the flow field in the separation zone. Due to the combustion process, the separation area becomes larger than the corresponding frozen cases and the separation shock moves upstream. This result highlights the chemical reactions influence on the separation shock displacement and therefore on the pressure rise within the recirculation zone. This result can be connected to the boundary layer, which becomes thicker in reactive flow and easier to detach in response to an adverse pressure gradient. This finding is quite consistent with previous works such as in [24]. Note that Test case 3 is not represented here in terms of comparison between the calculated and measured pressure profiles. Indeed, for this particular test case, large pressure fluctuations in the combustion chamber were recorded during the course of experiments. This led to great separation shock instabilities and large discrepancy in measurements.
Additional Test Case 4 was carried out only by numerical means, with an oxidizer/fuel ratio ( / = 6), close to that of real rocket engines, functioning with LO 2 -LH 2 . The nozzle pressure ratio NPR is fixed the same way as for Test Case 1; that is, NPR = 25.9, while the total pressure of the H 2 film cooling is set to 0 = 4.3 bars. As can be seen in Figure 7(c), augmenting the film cooling total pressure leads the separation shock to move downstream as well as the reignition attachment point consequently. This possibility can be used, as will be mentioned below, as means to avoid the reignition of the mixture inside the nozzle.

Chemical Production.
To highlight the reignition phenomenon in the separation zone, we plotted the flow OH mass fraction. The results obtained show that OH radical concentrations are at levels high enough to readily suggest  Figure 9. It is worth noting that these results were observed for both turbulence models, with modified Evans-Schexnayder chemical kinetics model but not with Eklund model. This can be explained by the temperature levels in the nozzle main flow, across the shear layer, which are not high enough to reactivate the chemical reactions with the above mentioned reaction model. On the other hand, when the nozzle operates at pressure ratios close to the adaptation, as in Case 3 (NPR = 36.5), no reignition process has been observed in the flow inside the nozzle. Indeed, as mentioned above, the separation zone size is too small that it does not allow the full completion of the combustion process within this area. In this case, the combustion takes place in the outer caisson, far away from the nozzle exit section. This situation is depicted in Figure 10 for both OH concentrations and temperature plots.

Thermal
Loads. The film cooling efficiency is evaluated by a widely used parameter , describing the difference between the film cooled wall temperature, at each point, and the combustion chamber temperature, relative to the difference between the coolant temperature and the combustion chamber temperature. The cooling efficiency is expressed by the nondimensional efficiency defined as where is the calculated wall temperature, the calculated combustion chamber temperature, and the injected film temperature.
The cooling efficiency along the wall depends, as before, on whether the nozzle flow is fully expanded or overexpanded. Figures 11(a) and 11(b) represent an example of wall temperature calculated in the Case 1 (NPR = 25.9) and Case 2 (NPR = 22.2), respectively, where the nozzle operates in overexpanded conditions. The combustion of injected GH 2 leads to an increase in wall temperature in the region where the flame is attached. The wall temperature plots are characterized by marked peaks up to 1600 K, compared to frozen flow calculations. These results represent a further confirmation of the hypotheses raised previously about the reignition phenomenon. As can be seen, these peaks are shifted according to the observations made above about the separation zone. In this case the cooling efficiency falls to values close to 0.05 and 0.4, respectively, for Case 1 and Case  2 (Figures 12(a) and 12(b)), where the temperature of the attached flame is close to the adiabatic temperature of the combustion chamber, particularly for Case 1. Meanwhile Case 3 shows a typical behavior of film cooling applications as it was observed in several works [2]. In this case the temperature of the film cooling grows gradually by mixing with the hot main flow gas (Figure 11(c)) and the corresponding efficiency ( > 0.8) is within the range of the expected values ( Figure 12(c)).

Conclusions
This numerical study was conducted with the goal of reproducing some experimental tests carried out previously. The main objective was to check to which extent the chemical kinetics and turbulence models may be relevant to highlight the reignition phenomenon, in the case of a subscale supersonic nozzle hot flow cooled by H 2 film cooling device. The main findings of this study can be summarized as follows.
Even though the real separated flow is naturally unsteady, the steady RANS approach, using -and Spalart-Allmaras models, combined with a modified kinetic model, was found to be able to reproduce the accurate averaged flow field and the reignition phenomenon observed experimentally.
Globally, the calculated wall pressures are in good agreement with the experimental data for the whole tested cases and the effect of chemical reactions on these quantities has been pointed out.
The computed steady reacting regions are slightly different from the instantaneous PLIF measurements which illustrate the complexity of this phenomenon. Furthermore, it is demonstrated that the steady attached flame occurs along the mixing shear layer and the combustion regime is mostly diffusion.
When reignition occurs, the temperature at the nozzle wall is such that the film cooling becomes locally inoperative and the induced thermal loads is so high that the nozzle's aerodynamic properties can be altered up to being detrimental for the nozzle structural integrity.
As expected, the choice of a turbulence model does not globally impact the phenomena that depend on the chemical kinetics, as the reignition process; however, it significantly affects the phenomena related to the boundary layer. This is how Spalart-Allmaras model gives a better prediction of the shock position and the pressure rise in the separation zone, compared to the experimental results.