Numerical simulation of reacting flow in the combustion chamber and study of the impact of turbulent diffusion coefficients

A methodology for combustion modeling with complex mixing and thermodynamic conditions, especially in thrusters, is still under development. The resulting flow and propulsion parameters strongly depend on the models used, especially on the turbulence model as it determines the mixing efficiency. In this paper, the effect of the sigma-type turbulent diffusion coefficients arriving in the diffusion term of the turbulence model is studied. This study was performed using complex modeling, considering the conjugate effect of several physical phenomena such as turbulence, chemical reactions, and radiation heat transfer. To consider the varying turbulent Prandtl, an algebraic model was implemented. An adiabatic steady diffusion Flamelet approach was used to model chemical reactions. The P1 differential model with a WSGG spectral model was used for radiation heat transfer. The gaseous oxygen (GOX) and methane (GCH4) operating thruster developed at the Chair of turbomachinery and Flight propulsion of the Technical University of Munich (TUM) is taken as a test case. The studies use the 3D RANS approach using the 60° sector as the modeling domain. The normalized and absolute pressures, the integral and segment averaged heat flux are compared to numerical results. The wall heat fluxes and pressure distributions show good agreement with the experimental data, while the turbulent diffusion coefficients mostly influence the heat flux.


Introduction
The complex procedure of nowadays design of propulsion devices now contains numerical simulation of different physical phenomena inside combustors as a necessary part of the process. This is because of the significant reduction of costs and time needed for the production and development cycle. The main requirements for rocket thrusters are precise prediction of combustion efficiency, wall heat transfer, specific impulse, and pressure. The modeling tool should both properly predict the parameter distribution and be robust enough and requiring only limited amount of resources to respond the engineering needs. Numerical codes, therefore, should be validated for this criteria and the most suitable approaches of physical phenomena modeling should be chosen. [1][2][3][4][5] For the several past years, at the Chair of Turbomachinery and Flight Propulsion (LTF) at the Technical University of Munich a lot of experimental rocket propulsion framework studies have been carried out.
Among these, a GCH4/GOx multi-element injector test case was chosen to validate the numerical calculations, as it is both well-described and suitable for numerical verification settings. A paper by Silvestri et al. gives detailed data concerning the test bench and experimental setting. [6][7][8] The discussed test case and the similar ones have already been numerically studied by several authors on the effects of turbulence and combustion modeling, which agreed on the primary influence of the turbulence model and then the combustion model on the calculated pressures and heat fluxes. Ivancic et al. 9 studied combustion process in a Penn State combustor and concluded that RANS approaches are still capable of predicting flow parameters in such test cases and showed that turbulence has the main effect on the efficiency and heat flux predictions. C. Roth et al. also discussed simulation results based on a single injector methane-oxygen case, where unless generally nice results the heat release was not computed sufficiently well, which was likely due to improper mixing modeling. 10 Chemnitz et al. 11 studying the same singleinjector case found that the differences in mixing behavior originate from the fields of turbulent viscosity, which means the key effect of the turbulence model behavior.
For the multi-injector methane-oxygen case, Perakis et al. 6,12,13 also showed the contrast in predicting integral parameters by different turbulence models probably originating from the different turbulent viscosity computation.
Currently the majority of RANS-based studies of reacting flow in combustion devices take the turbulence models with their empirical coefficients set by default in the (mostly commercial) codes. It is often assumed that the coefficients derived during the creation and validation process of the model would work smoothly with the applied model geometry and conditions. However, sometimes changes to at least some of the coefficients is needed and can lead to different results.
Studying reactive flow in a combustor, Poschner et al. applied variations of the multiplying coefficients in the turbulence model, which was addressed to obtaining more physical turbulence viscosity field. 14  Many studies of other combustion devices aimed at altering the turbulent diffusion coefficients in the energy (turbulent Prandtl number) and scalar transport (turbulent Schmidt number) equations. [19][20][21][22][23][24][25] All studies show high impact of these coefficients and the approaches to vary these values locally. Similarly, the turbulent closure diffusion coefficients (TCC) arriving in the diffusion term of the differential equations for turbulent scalars can be determining for the resulting eddy viscosity, which means resulting mixing intensity and finally the flow parameter distribution with the integral efficiency of the combustion process.
However, no studies discussed the impact of the TCCs. Therefore, in this paper, the main point is studying the influence of these coefficients based on the Shear-stress transport model. It must be mentioned, however, that changes in the turbulence model, in general, need to be first calibrated on special test cases, as any amplification in the model actually changes the model consistency and finally one gets a new model. Such calibrating studies are usually very consuming in time and sometimes omitted by the researchers of more complex processes such as combustion. The turbulence modeling community had considered this ambiguity and a new General equation k-omega (GEKO) model has been created, 26,27 which is only available in Ansys FLUENT now. This model applies to a larger amount of flows than any of the existing models by eliminating this consistency problem when needed. A user can calibrate the model himself for the type of flows needed by the values of ''representing'' coefficients. The study of the present test case with the use of the GEKO model is planned for the future. It should be mentioned that even the appearance of such turbulence model in a commercial code shows the current gap in engineeringapplied turbulence modeling methodology.
In this paper the TCC variation pursued two goals: (a) to observe the influence of the TCCs on the mean pressure and heat flux and, if large, that would mean a need to take a deeper insight into these values; (b) and/ or to show the need to use the GEKO model instead of the widely used SST k-omega model and others.
Moreover, one of the main concepts of this research was also using the mostly relatively robust and effortless models for the studied physical processes, which could be then used in the routine engineering design optimization process. Besides, most of the studies observed used the ''resolved to the wall'' approach with RANS. On the opposite, here an attempt to use a coarse near-wall refinement in such cases is presented. Also, a low-consuming algebraic model for the turbulent Prandtl variation was chosen, whereas a relatively robust and cheap weighted sum of gray gases (WSGG) approach was used to simulate the spectral radiative field coupled with the P1 radiation transfer model. Also, a robust and efficient adiabatic steady diffusion flamelet model was chosen to simulate the reacting mixture.

The test case
The view of the experimental combustor is given on Figure 1. One of the main features of the system is that it contains four water-cooled segments and one nozzle segment allowing to determine the wall heat flux. The total length of the combustor is 381 mm, which makes it possible for the injected jets to mix and react sufficiently.
An operating point according to the mixture ratio about 2.65 and the mean pressure of 18.3 bar was chosen for the numerical study. The mixture ratio is defined by the mass flows which are equal to 0.08 kg/s for CH 4 and 0.211 kg/s for O 2 . The mixing process of the components is induced by seven coaxial injectors, oxygen injected from in the inner part and methane from the outer. No swirling motion inside the injectors was reported, instead, the length of the injectors is supposed to make the gases sufficiently homogenously and smoothly distributed inside the injectors. The experimental data available for comparison includes the heat flux for each segment, the mean pressure the wall pressure distribution, the wall temperatures, and the methane/oxygen flow rates. In general, the test setup provides a wide area to amplify the numerical approaches to get best representation of both pressure and heat flux. Another interesting point is that it has the ability to use various mixture ratios, different fuels and a injector for the wall film cooling which makes it possible to study the influence of configurations and make the simulation approaches applicable to different modeling tasks.

Numerical setup
The commercial 3D Ansys CFX 28 coupled algebraic multigrid steady-state solver was used for the study. The Favre-averaged Navier-Stokes system of equations along with the equations for scalars was solved to get a solution.

Simulation domain
The 60°sector of the combustor was chosen for the presented simulations, which included one peripheral injector and 1/6 of the central injector. To consider the developed velocity profile, the injector pipes were also modeled. As depicted in Figure 2, the symmetry boundaries were applied to planes corresponding to 6 30°f rom the radial position of the outer injector. The application of symmetry boundaries has already been widely used in the combustor studies, both the same with the discussed and the others. 12,13,29 The mass flows and temperatures were defined for the fuel and oxidizer inlets and the outlet boundary had an ''opening'' typewith 1 atm pressure value. The no-slip condition with the experimental temperature profile as shown in Figure 3 was set at the thruster wall. The other wall faces were modeled with a no-slip adiabatic setting. Due to the previously discussed ambiguities in the nozzle heat flux estimations, 12,13 no comparison between the experimental and CFD data for the nozzle was produced, as the current research was focused on the fields and integral parameters in the cylindrical part of the combustor. The nozzle temperature was set equal to the last value of the temperature profile in the cylindrical part.
The approximate size of 3,4 M hexahedral cells was used for the numerical grid after a convergence study based on the average pressure and integral wall heat flux criteria, which can be seen in Table 1. Most of the research papers on the present and other typical combustors applied a resolving-to-the-wall grid to satisfy the y + ' 1 condition. 6,12,13,29 Some groups focus on the development of wall-functions, but now it is mainly coupled with LES based turbulence modeling. 10,11,29-35 Figure 1. View of the studied thrust chamber. 7,8 In this study, a default near-modeling approach being coupled with RANS was used, namely an ''Automatic wall treatment'', which CFX implements 28 and proposes a smooth shift from a resolved viscous sublayer to a coarser wall-function based method for logarithmic region, which can make the engineering optimization routine computations in industry more fast and easy. Thus, a mesh with y + ' 45 . . . 240 was used to assure that default wall function approaches were used. On Figure 4, a view on some of the mesh details is shown.

Numerical models
Turbulent diffusion closure coefficients. The turbulence model, as previously discussed, is determining for mixing of fuels, therefore coefficients which a model implies    may have an important role in the general operation of the combustor. As already mentioned, in this study, the changes to the TCC values of the SST model were made (1,2).
where s k1 , s k2 , s v1 , s v2 are sigma-type TCCs varied in this study. In the recent typical papers, 9,11 the variation of the turbulent Prandtl numbers was made in the region 0.6-1.1, which is approximately 67% or 120% of the default turbulent Prandtl number of 0.85-1 being usually used in the CFD codes. 27,28 As the TCCs and turbulent Prandtl number are both the turbulent diffusion-type coefficients, which determine the intensity of the turbulent diffusion of energy (in case of Prandtl) and turbulent scalar (in case of TCC), similar ranges were applied for TCCs. Four cases were studied, which is presented in Table 2.
Turbulent Prandtl number models. To consider for the variations of turbulent Prandtl, an algebraic model developed by Kays and Crawford (3) 36 was implemented. The model is strongly dependent on the eddy viscosity produced by the turbulence model, and this makes the study of the TCCs effect even more important as the turbulent Prandtl mainly determines the thermal diffusion, which can lead to different distribution of flow parameters inside and the integral heat fluxes. This model has recently been observed by D. Yoder 37 and showed satisfactory performance (compared to the differential models studied) for the near-wall, jet and pipe flow. Here, the main reason to use the algebraic approach was its robustness and simplicity.
Combustion model. An adiabatic steady diffusion flamelet approach was used to model combustion. The theory of this approach is well described in Ansys Theory guides. 28,38 Its main advantage is high computational performance while it still considers kinetics due to only two scalars to transport -the mixture fraction and mixture fraction variance. The turbulence-chemistry interactions are accounted for by the scalar dissipation rate. The C1 CH 4 /O 2 kinetic mechanism was used to create a flamelet database in CFX-RIF.
Radiation model. The P1 differential model was used to account for radiation due to its simplicity and wide applicability and a WSGG model was used to consider spectral absorption and emission. No scattering was included in this study as no particles were introduced; the studies of radiation heat transfer without scattering had been conducted before and provided decent results. 39 The absorption coefficients in the implemented WSGG model by Centeno et al. 40 are derived from the partial pressures of most radiating species, H 2 O and CO 2 . The radiation transfer equation would have a form 28 : dI y ðr;sÞ where yfrequency, r-position vector, s-direction vector, s -path length, K a -absorption coefficient, K ay À the blackbody absorption coefficient, I y -spectral radiation intensity, I b -blackbody emission intensity. The weighting factors were determined by the relation below: where b j, i are the polynomial coefficients taken from reference. 40 The absorption coefficients were introduced by the equation: where p CO 2 and p H 2 O are CO 2 and H 2 O partial pressures respectively, and k p, j is taken from reference. 40

Results and discussion
The normalized (Figures 5 and 6) and absolute pressure curves do not show any dependence on the TCC values. The absolute pressures are underestimated, as well as for the studies by other researchers. 10,12,13.29 It is seen that for the first 0.1 m of the combustor the normalized pressures also have an underestimation of the values, this area being the primary mixing area where the most of fuels are supposed to mix. However, the experimental peaks and cavity of the pressure in the vicinity of the injectors are reproduced. Here, the error for the pressures is less than 6.5%, which is also attributed to the underestimation of mixing intensity determined by the eddy viscosity distribution given by the SST model.
However, it can be still considered a satisfactory fit for the results.
Opposite to the pressures, the local heat flux and segment average heat flux values show (Figures 7 and  8) a slight influence for the varied range. In general, the TCCs can impact the observed heat flux in several ways: (a) the change of coefficients influence the eddy viscosity field, which is determining for variable turbulent Prandtl calculation by the KC model, which     Obviously, this TCC and turbulence model behavior is a reason for ambivalent values of heat fluxes near the injector and down the combustor. However, the distinctions of heat flux peaks in the injector vicinity are too small to influence the absolute and normalized pressure field. Interesting is also the total radiative heat flux presented in Table 3. The radiation heat flux fraction tends to rise while the TCCs decrease. This is explained by the decay of the convective and conductive heat transport due to the changed turbulent Prandtl, which results in higher rates of radiative flux. This outcome is an interesting feature to be studied in the upcoming research.
The observed effects show the need either to have an attentive look at the TCC values which would mean more studies with other turbulence models, operating conditions, etc.; or to use the new GEKO model and calibrate it for such case and/or develop similar models keeping the whole turbulence model consistent.

Conclusion
A variation in the coefficients of turbulent diffusion closure had a significant effect on the distribution of the heat flux of the wall, while no difference was observed in the calculated pressure and a slight difference for radiative to total heat flux ratio was noted. This is mainly attributed to the application of a variable turbulent Prandtl approach and therefore contradictory impact of eddy viscosity on the heat fluxes. However, the results for the default TCC values show a good match with experimental data, both on wall heat flux and pressure profiles, which also shows that even with the use of the default wall-function-based near-wall modeling approaches it is possible to get decent results in terms of pressure and heat fluxes at least for engineering design. It can be also mentioned that there is a need to use the GEKO model or similar approaches for future studies to get a more insightful view of the impact of other turbulence model parameters and handier calibration for combustor cases.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was funded by the Russian Ministry of Education and Science (project 13.7418.2017/8.9).