Extended Eddy-Dissipation Model for Modeling Hydrogen Rocket Combustors

ABSTRACT In the work, an extension of the eddy-dissipation model (EDM) is developed to simulate turbulent combustion of hydrogen in undiluted oxygen in rocket combustion chambers. The modification of the eddy-dissipation model allows eliminating of main demerits of the original EDM model. This is achieved by introducing additional parameters into the model, which limit the reaction rate and depend on the local stoichiometry and temperature. The main such parameter is “Maximum flame temperature,” which depends on local stoichiometry and takes into account the dissociation of combustion products. The extension of the EDM model is based on the framework provided by ANSYS CFX. The new turbulent combustion model is validated against experimental data from three different subscale rocket combustors. The validation of the model is carried out against data on pressure and wall heat flux, which are the main targets of simulations of rocket combustion chambers.


Introduction
Most of us associate rocket combustion chambers with high turbulence, pressure, and temperature, which is true. At high pressure and temperature, chemical reactions are so fast in such a way that they allow to use the assumption of thin flame, that is, of infinitely fast chemistry. This assumption holds very well for hydrogen in rocket combustion chambers (Ivancic and Mayer, 2002). The thin-flame assumption greatly simplifies the modeling of turbulent combustion as now there is no need to solve kinetic equations, as they are infinitely fast. However, this does not mean infinitely fast combustion. Now, the combustion rate is limited by other processes: turbulent mixing or fuel evaporation. In the most turbulent combustion models, which are used in the thin-flame assumption, the combustion rate is proportional to the rate of turbulent mixing. The most popular model of turbulent combustion in computational fluid dynamics simulations is the eddy-dissipation model by EDM. The model is attractive not due to the accuracy but simplicity. The combustion is described as a single-stage process. The present work deals with the case of the combustion of pure hydrogen and oxygen, so The reaction rate is given by the following expression: where A is a model constant, which usually equals to 4, and k are turbulent eddy dissipation and turbulent kinetic energy, respectively. The model is so crude that the authors of the EDM model tried to solve some of the problems of the model already in their original work (Magnussen and Hjertager, 1977). They introduced a product limiter, which makes reaction rate dependent on the concentration of reaction products when their concentration is small enough. However, this limiter is not used here.
The EDM model has a clear physical interpretation: the reaction rate is proportional to the turbulent mixing timescale and to the average concentration of a deficient reactant. In contrast to other combustion models and other similar models like the eddy breakup model by Spalding (1971) and the eddy-dissipation-concept model by Magnussen (1981), there is no complicated kinetic mechanism with intermediates in the model, and there is no splitting in different turbulent scales. There is only one global reaction step and one general turbulent timescale.
The main disadvantage of the EDM model in rocket combustion chamber is that it gives a very high flame temperature. Actually, the prediction of gas temperature is the main goal of computational fluid dynamics (CFD) simulations of rocket combustion chambers. Nowadays, it is also possible using CFD simulations to predict the dynamics of combustion processes in rocket combustion chambers: engine firing and onset of combustion instability. However, the main thing, which is required from CFD simulations for the design of rocket combustion chambers, is an accurate prediction of thermal loads. The accurate prediction of heat fluxes requires the accurate simulation of temperature field in combustion chamber. Thus, the main requirement for combustion model is the accurate prediction of the temperature of burned gases. The direct use of Equation (1) gives a flame temperature near 5000 K while the flame temperature in rocket combustion chambers amounts to around 3500 K. The temperature of H 2 /O 2 flame is significantly lower than 5000 K due to the fact that the significant part of H 2 O dissociates at a temperature above 3500 K. At T ¼ 3600 K and p ¼ 60 bar, the equilibrium composition of water vapor in mole fractions is following, according to Vargaftik (1975): The dissociation of water vapor can be taken into account through a kinetic mechanism with several reactions and intermediates: OH, H, O 2 , etc. However, the EDM model does not assume multistep reaction or mechanism with many reactions. The direct use of the EDM model with a kinetic mechanism of many reactions results in nonphysical results because Equation (2) results in a wrong balance between rates of different reactions and, consequently, in a wrong chemical equilibrium state. To overcome this and other demerits of the original EDM model, the following extension of the model is proposed.

Extension of the eddy-dissipation model
All simulation results presented in this work were obtained using the commercial computational fluid dynamics code ANSYS CFX (2012b). The proposed extension of the EDM model is based on the framework provided by CFX for this model. CFX also allows users to define their own constants, expressions, functions, routines, etc. To solve the problem of too high flame temperature, CFX offers a simple solution that is to specify explicitly a maximum flame temperature. When the temperature of gas exceeds this parameter, the reaction rate is set to zero. However, this does not solve the problem completely. The temperature of burnt gases will exceed the adiabatic flame temperature in areas where the mixture equivalence ratio is not the same as the average equivalence ratio. To overcome this problem, it is necessary to set "Maximum flame temperature" dependent on mixture composition. In the present case, a new parameter that reflects a local equivalence ratio was defined. This is the mass ratio of oxygen to hydrogen in the mixture (or abbreviated MROH): The additional term "0.01" is needed to avoid the situation when MROH is equal to zero or infinity. The adiabatic flame temperature is calculated using CEA (McBride and Gordon, 1996) or other similar code, see Table 1. (The focus of the present work is hydrogen rocket combustion chambers, which are usually cryogenic; thus, the initial temperature of the propellants is 160 K in Table 1.) However, the effect of the variable "Maximum flame temperature" instead of a single value (e.g., 3660 K) is not large and is shown in the work by Zhukov (2015). The difference between the single value and the variable "Maximum flame temperature" is visible in very lean or reach mixtures. In the EDM model, the reaction rate does not depend on temperature. This causes two problems: combustion of propellants even at cryogenic temperatures and a sharp drop of the reaction rate behind the flame front when the gas temperature exceeds the "Maximum flame temperature." The first problem is solved by a "standard" approach: the introduction of an auxiliary parameter called "Extinction temperature." When the fuel mixture has a temperature below the "Extinction temperature," the reaction rate is set to zero. This parameter should reflect the flammability limits of H 2 /O 2 mixture at low temperatures. There is no literature data on the flammability limits at cryogenic temperatures. However, an idea about the flammability limits can be obtained by modeling a laminar free premixed flame. In this study, we used the kinetic model by Burke et al. (2012), which is the closest to experimental data at the pressures above 10 atm, and the computer code PREMIX (Kee et al., 1985), which is the part of the software package CHEMKIN. As a result of the analysis of H 2 /O 2 flames, the following expression has been used for the "Extinction temperature": T ext ¼ 170 ½K Á ððlnðMROHÞ À 1:7Þ 2 þ 1Þ: This expression is approximate and is only a convenient function. In reality, the value of the "Extinction temperature" can be set by any other method or function.  Explosion limits at 100 bar according to Schröder and Holtappels (2005).
The sharp drop of the reaction rate behind the flame front does not adversely affect the simulation results itself, but the presence of such singular points decreases the convergence of the solver. This problem was solved by setting model constant A dependent on temperature and the mass fraction of water vapor Now, the reaction rate smoothly decreases at temperatures near T ext and when the mass fraction of water (the product of the reaction) approaches 100%. This modification of the model results in a smoother spatial distribution of reaction rate across the chamber and allowed a tenfold decrease of mean residuals.
The present extended EDM model has two additional parameters for modeling the interaction between flame and turbulence. The turbulent mixing rate =k becomes large close to walls due to the drop of k. Therefore, the value of =k in Equation (2) is limited to a value of 5 Á 10 3 s À1 , which is set by a parameter called "Mixing Rate Limit" in (CFX, 2012a); otherwise, the reaction rate goes up unnaturally near walls. At a certain level of turbulence, the dissipation of heat and radicals leads to flame quenching. In regions of high turbulence, when the turbulence time scale is smaller than a chemical time scale, local extinction occurs in the present model, that is, the reaction rate is set to zero. The chemical time scale is defined in the present model as a ratio of laminar flame thickness δ to laminar flame velocity S u . The flame thickness is evaluated using Blint's correlation (Blint, 1986): where indexes u and b denote unburnt and burnt states, respectively. The laminar flame velocity is calculated using the already mentioned mechanism of Burke et al. (2012) and flame code PREMIX (Kee et al., 1985). In our earlier studies (Zhukov, 2015;Zhukov and Suslov, 2016), turbulence mixing time scale k= was used for comparison with the chemical timescale; however, later the Kolmogorov time scale has been used instead of the mixing time scale (τ mix ¼ k=). This mechanism provides flame extinction near the injector tip. The weaker the flame near injector, the more oxygen crosses the flame. This slightly changes the distribution of oxygen over the cross-section of the combustor. The present simulation results obtained using the Kolmogorov time scale are marginally better than the earlier results (Zhukov, 2015;Zhukov and Suslov, 2016).

Numerical modeling
It necessary to say a few words on numerical models and setups used here. As it was mentioned earlier, the modeling was done using ANSYS CFX (2012b). The flow is modeled using the Favre averaged Navier-Stokes equations. Turbulence has been modeled using the SST k À ω turbulence model using the standard values of the coefficients and the "automatic" wall treatment (CFX, 2012b). The transport has been modeled with a turbulent Schmidt number of 0.7 (the value of 0.7 is recommended for high-Reynolds-number jet flows by Yimer et al. (2002)). The turbulent Prandtl number was set to the default value in CFX that equals 0.9. All components of the gas mixture (H 2 , O 2 , and H 2 O) have significant distinctions from ideal gas under conditions typical for rocket engines; therefore, an accurate modeling of thermodynamic properties is required. Cryogenic O 2 and H 2 are modeled using the Peng-Robinson real gas equation of state. The enthalpy and the entropy of the individual components have been defined using NASA polynomials (Burcat, 2001). The dynamic viscosity and the thermal conductivity of H 2 , O 2 , and H 2 O have been defined using Sutherland's law with coefficients recommended by White (1991). The diffusion coefficients have estimated using the data from Kikoin (1976). The mixture molecular transport properties are calculated from the properties of the individual components using the equation by Mathur et al. (1967) for the thermal conductivity and Wilke's formula for the viscosity (Wilke, 1950), that is, using the same equations as Chemkin (Kee et al., 1986). The accuracy of both mixing rules (Mathur et al. and Wilke) in undiluted H 2 /O 2 mixtures is good enough according to Zhukov and Pätz (2017).

Penn State test case
A test case, which became the first well-known test case for the validation of CFD models for rocket combustion chambers, got the unofficial name "Penn State test case" (also known as RCM-1) (Pal et al., 2006). This test case is a starting point for CFD modeling of rocket combustion chambers due to its relative simplicity. By this reason, this test case was modeled by many researchers: Tucker et al. (2008) Zhukov (2015). Two works (Ivancic et al., 2013;Tucker et al., 2008) are collaborative works where authors compared different CFD approaches.
The configuration of the experimental setup corresponds to a staged combustion cycle operating with gaseous oxygen and hydrogen propellants (Pal et al., 2006). The experimental setup consisted of two preburners and a main combustion chamber. The main combustion chamber has a single co-axial injector and is fueled by fuel-rich and oxygenrich preburner gases. Within this test case, the results of wall heat flux measurements are supposed to be a target of simulations.
In Figure 1, one can see a comparison of different CFD models with the experiment. The results that are obtained using the extended EDM model are named as CFX, DLR-LA (the used code, and the German abbreviated name of our institute). The selection of other simulations is determined by the similarity with our numerical model. All simulation results shown in Figure 1 were obtained in quasi-two-dimensional axisymmetric domains by solving the Reynolds-averaged Navier-Stokes equations (RANS). The numerical setup used by us is described in detail in the earlier work by Zhukov (2015). The results presented here are slightly better than those presented earlier due to the use of the Kolmogorov time scale instead of the mixing time scale for the modeling of extinction at high turbulence. Figure 2 shows simulated flow and temperature fields in the Penn State combustor. The simulation results look typical for RANS simulations. Many researchers were not very successful in modeling of heat transfer in this test case: Karl and Hannemann (2006); Tucker et al. (2008); Lian and Merkle (2011);Lempke et al. (2015). The main critical point in this test case is the modeling of the recirculation zone in the corner between the sidewall and the front wall because the maximum of heat flux is located near the flow attachment point. If a model correctly predicts the size of the recirculation zone and the flame temperature, then the simulation results automatically lie near the experimental points. That is why, the accuracy of the simulation results given by our model is not surprising because the SST turbulence model used here shows very good results in the prediction of the recirculation zone behind a backward facing step (Ul Haque et al., 2007). Indeed, simulation results in the Penn State test case are not very sensitive to used combustion model. The simulations presented in Figure 1 use different combustion models ("Rocflam3, Astrium 1 "-equilibrium based PPDF, "CFX, Astrium"-flamelet, "TAU, DLR-GÖ 2 "-finite rate chemistry without turbulencechemistry interaction, and "CFX, DLR-LA"-the extended EDM model), but the results are close to each other. This is due to the fact that the Penn State combustor was fueled by hot  partially burned gases, which immediately react with each other, so the duty of the combustion model is only to predict accurately the temperature of burnt gases.

Combustion chamber with porous injector
The second test case is much more complex than the Penn State test case. The details of this case are published by Zhukov andSuslov (2016, 2015). The test case has a feature that is a porous injector head. In this test case, liquid oxygen (LOx) is fed through many distributed single injectors while cold hydrogen (100 K) is fed into a combustion chamber through a porous plate. This injector is a new concept under development at the Institute of Space Propulsion. The concept has some advantages over coaxial injectors, which are conventional for rocketry. According to the hot-tests at the Institute of Space Propulsion of the German Aerospace Center (DLR-Lampoldshausen), the porous injector head (Figure 3) allows to maintain the high combustion efficiency over the wide throttling range from 40% to 130% (Deeken et al., 2011). Besides the low manufacture cost and the throttling capability, porous injector head has two additional advantages over conventional coaxial injectors. Porous injector head operates at a smaller pressure drop than injector heads with coaxial injectors. Second, the small diameter of the injectors in a porous head results in a small jet break-up distance that allows reducing the combustor length.
The numerical setup used in this test case is described in detail in the earlier work by Zhukov and Suslov (2016). However, the present results are obtained on a much refined mesh with 12.5 Â 10 6 nodes and using the Kolmogorov time scale for modeling extinction at high turbulence. Liquid oxygen was treated in the simulations as a real gas that obeys the Peng-Robinson equation of state. In the earlier study, Zhukov et al. (2011) showed that simulations of the combustion chamber with the porous injector head should be carried out in a three-dimensional (3D) formulation. To have a correct representation of the arrangement of oxygen injectors, the simulations were carried out in the eighth part of The transition from 2D to 3D problem formulation significantly increases the computational cost. Certainly, using a heavy combustion model with many transport and differential equations would make it impossible to simulate the combustion chamber in 3D geometry using a single workstation, in this case Dell T7500 with two Intel Xeon E5645 processors. Other options for limited computational resources can be the use of the flamelet approach or the equilibrium chemistry model (Ivancic et al., 2013). Both these combustion models also use the thin-flame assumption; however, they both also assume unlimited mixing rate or infinitely fast chemistry at "mesoscopic" level. Moreover, the adiabatic flamelet model assumes no further reaction behind the flame front. Zhukov and Suslov (2016) simulated this test case (chamber "B" with the porous injector) also using the flamelet approach, in this case the Extended Coherent Flame Model (ECFM) by Colin and Benkenida (2004), which includes also  turbulence-chemistry interaction. The comparison of the models (extended EDM and ECFM) shows that the ECFM model predicts pressure in the combustion chamber lower by 1-1.5 bar while the prediction of the EDM model is within the error margins, see Figure 5. The low pressure in the combustion chamber with the ECFM model is explained by the absence of chemical reactions in burnt gases, namely in the nozzle. From the equations given in Rocket Propulsion by Barrère et al. (1960), it is possible to connect combustion chamber pressure p c with the speed of sound c th and temperature T th in the throat where _ m is the total propellant mass flow rate, A th is the throat area, and μ is molar mass. The extended EDM model predicts higher temperature in the throat, because the mixture retains the reactivity in the nozzle because the temperature there is below the "Maximum flame temperature." In contrast, the flamelet combustion model means a chemically frozen flow in the nozzle. The lower temperature in the nozzle given by the ECFM model also results in the too low wall heat flux in the nozzle (Zhukov and Suslov, 2016) while the results of the extended EDM model agree with the experimental data, see Figure 6.
The model is also validated at different oxidizer-to-fuel ratios (ROF). The results are presented in Figure 6. The difference between simulation and experiment exceeds the experimental error, which is not less than AE 3 MW/m 2 , but is still in an acceptable range. The accuracy of wall heat flux predictions depends not only on the accuracy of combustion model but also on the accuracy of the boundary layer modeling, which is very difficult in this particular case where flames of some injectors penetrate into the boundary layer (Zhukov and Heinrich, 2018).

Single coaxial injector combustion chamber
The third test case simulated here was obtained at test facility P8 at the German Aerospace Center (DLR-Lampoldshausen). The test case was presented on conferences by Suslov et al. (2015Suslov et al. ( , 2016), and it was already simulated by Seidl et al. (2017) and by Fechter et al. (2017). Although the combustion chamber has only one coaxial injector, this test case is not easier for CFD modeling than the previous one. Now, both propellants are injected at cryogenic temperatures, and a significant part of the fuel is used for a film cooling and injected at a very high speed.
The used combustion chamber was specifically developed and manufactured for intrachamber studies of injection and combustion by providing optical access for the application of optical diagnostics. A sketch of the chamber is shown in Figure 7. The combustion chamber is segmented into four interchangeable water-cooled sections and a nozzle. Thus, the instrumented section can be placed at various axial locations. This feature has been used to achieve optical access along the full combustion chamber length from x ¼ 0 mm to x ¼ 370 mm. This distance is the same as the length of real rocket engines up to nozzle throat. The inner diameter of the presented subscale combustion chamber is 50 mm. The detailed descriptions of the combustion chamber can be found in the original works by Suslov et al. (2015Suslov et al. ( , 2016 and by Mayer et al. (2001).
The hot-fire tests were carried out at pressures of 40, 50, and 60 bar for three different ratios of oxidizer to fuel at the injector (ROF inj ): 4, 5, and 6. Rearranging the sections of the combustion chamber, pressure and temperature measurements were obtained at different locations along the chamber axis. To carry out the measurements over the whole length of the chamber for the nine different load points, the hot-fire tests were repeated a considerable number of times. The reproducibility of the test parameters for different chamber configurations amounts to AE 2%.
The simulations have been performed in a 2D axisymmetric domain using the same numerical setup as before; however, all components (O 2 , H 2 , H 2 O, and N 2 (,1%)) are modeled using the Peng-Robinson real gas equation  Figure 7. DLR subscale combustion chamber model "C" .
has supercritical pressure and ROF inj = 6, has been simulated. The simulation results are presented in Figures 8 and 9.
The model shows a good agreement with the experimental data by Suslov et al. (2015). The simulation is complicated by the presence of a massive film cooling and by an uncertainty in the wall heat flux, which was not measured in the experiment. The cooling film injected at high speed and the flame of the coaxial injector generate an unstable system of vortices between theses two flows near the front wall. Small changes in the flame or in the film result in a significant change of the arrangement of the vortices, and this significantly slows down the solver convergence. Nevertheless, the numerical model has delivered good results.

Discussion
The extended EDM model was successfully validated against the experimental data from the three different rocket combustors. The simulation results show no problems or  nonphysical behavior associated with the low injection temperature of propellants. In other words, the model shows the good performance at the different injection conditions. Thus, the developed model can operate in the wide ranges of fuel-to-oxidizer ratios and of injection temperatures starting from cryogenic. Another advantage of the new model is that it can be easily adapted for methane, which is considered as a future rocket fuel by ESA, SpaceX, Blue Origin, and Roskosmos.
The new model was compared with the flamelet based model, but there is also a need for comparison with the equilibrium chemistry model by Preclik et al. (2004). It also utilizes single-step global reaction scheme, but the reaction rate is not limited and equals to infinity. Using the equilibrium chemistry model, Preclik et al. (2004) simulated fullscale rocket engines and obtained a good agreement with experimental data on wall heat fluxes in subscale rocket combustors. However, the assumption about the infinite reaction rate may lead some problems. This assumption is very strong. It is certainly not valid close to walls, in regions with low temperatures or where the mixture is very lean or rich, and it is not valid for hydrocarbons.
The EDM model can be adapted for the use with scale-resolved turbulence models: SAS, DES, and LES. In this case, the simulated structures will transport reactants toward the flame front while the modeled structures will be responsible for the limiting mixing within the flame front. The mixing rate can be defined in the same way as it is defined by ANSYS in the realization of the EDM model for LES in Fluent (2009): The developed method also can be used for other fuels: CH 4 , C 3 H 8 , etc. However, this requires, foremost, that the assumption of the thin flame holds. The parameters of a new extended EDM model can be calculated with a kinetic mechanism, which is the most suitable for a chosen fuel and conditions. The extended EDM model has shortcomings as well. First, it inherits all weaknesses of the thin-flame assumption and depends on the validity of this assumption. The validity of the model is questionable outside rocket engine conditions, which mean a non-premixed combustion of high energetic propellants at high pressures and temperatures. The model considers the combustion of propellants as a single-stage process without the formation of intermediates. There is no self acceleration of reaction due to the formation of radicals in the model. In the absence of OH and CH in the model, the comparison of simulation results with experimental data on flame chemiluminescence requires additional modeling.
Temperature makes the most impact on the reaction rate; however, the reaction rate (Equation (2)) in the EDM model, as in many other turbulent combustion models, does not depend on temperature. The reaction rate should be low at low temperatures and high at high temperatures; close to the chemical equilibrium state, the reaction should slow down. The proposed extension of the EDM model partially resolves this problem; however, it cannot offer the same accuracy as a kinetic model with many reactions.

On the role of detailed kinetic mechanisms
The eddy-dissipation model exploits the assumption of the thin flame. Indeed, this assumption is valid in hydrogen rocket combustors (Ivancic and Mayer, 2002;Tucker et al., 2008). One may think it means that the detailed chemical kinetics plays no role, as it is infinitely fast. However, it is not true. The extension of the EDM model is required precisely because the original EDM model disregards the detailed kinetics of the process. Despite the apparent simplicity of the H 2 /O 2 system, it is not simple. First, to accurately estimate the temperature of hydrogen-oxygen flame, it is necessary to know the intermediate species. As we saw in the Introduction, this requires taking into account the formation of at least O, OH, and H.
Second, both parameters connected with the extinction ("Extinction temperature" and "Chemical timescale") are calculated with the use of the detailed kinetic mechanism, in this case by Burke et al. (2012). The extinction cannot be taken into account without the knowledge of the kinetics. The extinction temperature is approximated by the "convenient" function; however, it is based on the kinetic modeling. It has a bell shape on a logarithmic scale, and its ends are where the speed of the premixed H 2 /O 2 flame goes down below 0.2 m/s.

Conclusions
The extension of the EDM model has been developed for the application in hydrogen rocket combustors. The extension resolves the major shortcomings of the EDM model. The extended EDM model was successfully validated against experimental data on wall heat fluxes and pressures in the three different combustion chambers in the wide ranges of oxidizer-to-fuel ratios and of injection temperatures.
The provided coefficients and data allow to use the extended EDM model for modeling hydrogen rocket combustors at pressures about 80 bar directly or after corresponding correction of the "Maximum flame temperature" at any other pressures. The given approximations for "Extinction temperature" T ext and for model constant A are "convenient" functions; thus, they can be defined by other ways. The chemical time scale can be defined by an alternative way as well.
The advantages and disadvantages of the new model in comparison to other models are described in this paper. The proposed methodology also allows extending the EDM model for methane. In the paper, it is shown how the new model can be also used with scaleresolved turbulence models: SAS, DES, and LES.