Rocket Combustion Chamber Simulations Using High-Order Methods

High-order spatial discretizations significantly improve the accuracy of flow simulations. In this work, a multi-dimensional limiting process with low diffusion (MLPld) and up to fifth order accuracy is employed. The advantage ofMLP is that all surrounding volumes of a specific volumemay be used to obtain cell interface values. This prevents oscillations at oblique discontinuities and improves convergence. This numerical scheme is utilized to investigate three different rocket combustors, namely a seven injector methane/oxygen combustion chamber, the widely simulated PennState preburner combustor and a single injector chamber called BKC, where pressure oscillations are important.


Introduction
Accurate and reliable predictions of quantities such as the wall heat flux, the pressure fluctuations or the flow field are essential in the design process of rocket thrust chambers. In high-pressure environments, experiments, although being the most credible approach, are often restricted to single measurable quantities [16]. In contrast, numerical simulations provide extensive data sets.
With growing computational resources, three-dimensional and time-resolved combustion chamber simulations become more and more attractive. For this, however, numerical schemes are required that are capable of capturing the occurring phenomena. Amongst others, one challenge for the numerical scheme that strongly affects the result is the discretization of the inviscid fluxes and therefore the cell interface value reconstruction. High-order approaches improve the accuracy and are therefore demanded. This becomes even more relevant when performing large-eddy simulations (LES). Higher-order methods may provide results that can be achieved with lower-order methods only on significantly finer grids. In addition, at supersonic speed, it needs to be taken into account that no new extrema must occur. This is ensured if the method fulfills the total variation diminishing (TVD) criterion [8]. Classical approaches as the van Leer or the minmod limiter are applied separately in every coordinate direction. Hence, the cell interface values do not take information from the values of diagonally located cells. Therefore, Kim and Kim [11] introduced a multi-dimensional limiting process (MLP) that considers the latter. In addition, convergence is improved. Gerlinger [5] extended the scheme to obtain results with as little diffusion as possible. This approach is called MLP with low diffusion (MLP ld ). In this work, MLP ld with up to fifth order accuracy is used. However, MLP ld is not restricted to a specific accuracy order. Here, steady and unsteady calculations of rocket combustion chambers are carried out with MLP ld .
The remainder of this work is structured as follows. In the next section, the applied numerical method is presented. This includes a detailed description of the cell interface reconstruction using MLP ld . Combustion chamber simulations are performed and examined in Sect. 3. Section 4 gives a short summary.

Numerical Method
Over more than twenty years, the in-house code TASCOM3D (Turbulent All Speed Combustion Multigrid 3D) has been developed and successfully applied to a wide range of reacting as well as non-reacting subsonic and supersonic flows, e.g. [4,7,25]. The following subsections shortly describe the underlying equations and numerical methods.

Governing Equations
Turbulent flow and combustion processes are described by the compressible Navier-Stokes equations that, in addition, include equations for turbulence modeling and species transport. The set of governing equations is given by Here, Q denotes the vector of conservative quantities consisting of the density ρ, the velocities u, v, w in each direction, the total specific energy E, the turbulence quantities k and ω, the temperature variance σ T , the sum of the variances of all species mass fractions σ Y as well as the N k − 1 independent species mass fractions Y i . N k is the number of considered species. F, G and H indicate the inviscid fluxes in x-, y-and z-direction, respectively. F v , G v and H v are the corresponding viscous fluxes. The source vector includes terms resulting from turbulence and chemistry. The species source terms include k f r and k b r which denote the forward and backward reaction rates of reaction r . Furthermore, ν i,r and ν i,r represent the stoichiometric coefficients of species i for the forward and backward reaction, respectively. N r stands for the number of reactions, M i is the molar mass of the specified species, and c l the concentration of species l. A virtual species N k + 1 is introduced to account for three-body reactions.
To close this set of equations, an equation of state is required. This is either the ideal gas law or, for real gas effects, the Soave-Redlich-Kwong (SRK) equation.

Numerical Solver
The set of governing equations (1) is solved by an implicit lower-upper symmetric Gauss-Seidel algorithm [19,20] using a finite-volume approach which works on block-structured meshes. In steady-state simulations, the solution is advanced in time with a first order temporal discretization until convergence is achieved. Timeaccurate simulations utilize a second or third order BDF (backward differentiation formula) scheme [3] with a number of inner Newton iterations.
The inviscid fluxes are calculated with the AUSM+-up flux vector splitting method of Liou [14]. This or any other flux-vector splitting approach requires values for the variables at both sides of a cell interface. These values are determined with a highorder scheme as described in Sect. 2.3. The viscous fluxes are calculated by central differences in a cell-oriented coordinate system.
Throughout the course of this paper, various rocket thrust chambers are simulated with different levels of turbulence modeling. While Reynolds-averaged Navier-Stokes (RANS) simulations offer great simplicity, LES are more accurate, but, accordingly, also more tedious. Moreover, the impact of high-order discretizations is small for smooth steady-state RANS simulations, while its impact is large for time-resolved calculations. The applied RANS model is the q-ω model of Coakley and Huang [2]. Note that the variable vector (2) in this case contains the turbulent velocity scale q = √ k instead of the turbulent kinetic energy k. As a wall-resolved LES for rocket thrust chambers is extremely tedious and costly, the improved delayed detached-eddy simulation (iDDES) of Shur et al. [21], which belongs to the class of hybrid RANS-LES models, is used. The near-wall region is treated with an underlying unsteady RANS k-ω model, whereas the rest of the computational domain operates in LES mode.
To model combustion processes, finite-rate chemistry is employed. The corresponding reaction mechanisms are given in the respective sections. In addition, an assumed probability density function (APDF) approach [4] is used to account for turbulence-chemistry-interaction. As the used concept assumes statistical independence of temperature and species fluctuations, the joint pdf of temperature and species composition can be simplified into a product of the individual pdfs of temperature and composition. For the first one, a clipped Gaussian distribution which is defined by T and the temperature variance σ T is assumed. The joint pdf of all species concentrations is described by a multi-variate β-distribution using the mean mass fractions and the sum of all species mass fraction variances σ Y .

Cell Interface Value Reconstruction
AUSM+-up requires the values of the primitive variables at the cell interfaces. Using polynomial reconstruction, higher order schemes can be obtained. Here, the MLP technique is used [5,11]. This approach is an extension of the conventional second order limiters to higher order schemes while considering all surrounding neighbor cells.
Without loss of generality, let the interface (i + 1/2, j, k) be the interface of interest. Then, the left, q L i+1/2, j,k , and right, q R i+1/2, j,k , interface values need to be calculated. The term q stands for any of the primitive variables present in the vector Q in Eq. (2) as well as the the total enthalpy H and the integral specific heat ratio γ . In order to improve the accuracy of the flux calculation, a high-order approach is used. The interface values are reconstructed by Depending on the preferred order, the functions β L and β R in (5) contain information from various neighbor cells. For example, the discretization stencil of the fifth order upwind biased scheme used in this work contains three upwind and two downwind cells. In case of equidistant grid spacing, these functions can be derived by a polynomial reconstruction Here, r L i, j,k = q i+1/2, j,k / q i−1/2, j,k and r R i, j,k = 1/r L i, j,k denote the required slope ratios. If non-equidistant grids are used, the coefficients in (6) become grid size dependent. However, as these coefficients only depend on the grid, they can be calculated in advance.
Equation (5) may lead to oscillations at discontinuities and the creation of new maxima or minima at the cell corners. Therefore, limiters are required in order to disable these effects. For third and higher order reconstructions, Kim and Kim [11] introduce a TVD limiter instead of β L i, j,k in the first equation of (5). Calculating the parameter α i (i ∈ {x, y, z}) for each coordinate direction is the basis for the MLP concept. MLP ensures that no new extrema can occur. This is done by checking all neighbor cells, including those located diagonally [5,13].
In order to validate the proposed method, a simulation of Schardin's problem, which consists of a planar shock that impinges on a wedge, is conducted using the proposed fifth order spatial discretization in combination with a third order temporal discretization. The results, depicted in Fig. 1, show excellent agreement with experimental data [1].

Combustion Chamber Simulations
This section presents one RANS and two iDDES rocket combustion chamber simulations. The first one, a seven injector combustion chamber is designed to improve the understanding of injector-injector and wall-injector interactions. The next two test cases exhibit instationary behaviors and are therefore tackled with iDDES. One of them uses the ideal gas law, the other one the SRK equation of state. In those cases, high-order spatial discretization has a large impact [13]. Some exemplary results are prescribed in order to demonstrate the applicability of the proposed high-order method for steady and unsteady combustion. All cases represent laboratory-scale chambers. However, in the course of this project, a full-scale combustion chamber with 90 injectors has also been simulated (not shown here).

Multi-injector Combustion Chamber
This combustion chamber was experimentally investigated at the Technical University of Munich [22]. Seven coaxial injectors supply gaseous oxygen and gaseous methane at an oxidizer-to-fuel ratio (O/F) of 2.65. One injector is at the middle of the faceplate, the others surround it at constant distance. The nominal chamber  [1] pressure is 18.5 bar. The inner diameter of the chamber is 30 mm, its nozzle diameter 19 mm and the length 341 mm. More details concerning the experimental setup can be found in [22].
As only steady-state RANS simulations are performed, only a 90 • segment is simulated, which includes 1.75 injectors. The numerical grid that covers the injectors, the chamber and the convergent-divergent nozzle consists of 82 blocks with approximately 7.9 million volumes. Symmetry boundary conditions are applied at both symmetry planes. The wall of the combustion chamber is assumed to be isothermal with wall temperature values stemming from the experiment [22]. At the inlet, a mass flow rate is prescribed and the outlet uses a supersonic outflow condition. The 21 species, 98 reactions methane reaction mechanism of Slavinsakaya et al. [23] describes reaction kinetics.
The wall pressure distribution is depicted on the left-hand side of Fig. 2. All facets of the experimental values are reproduced correctly. This includes the sharp rise at the beginning as well as the different gradients of the decline at the end. However, the experimental pressure is overestimated by about 1%. The right-hand side of Fig. 2 shows a comparison of the simulated and the experimental wall heat flux. As the latter is measured with a caloric method, only averaged values within a segment are available. The simulation shows good agreement with the experiment, too. However, in the middle of the combustion chamber the wall heat flux is overestimated, whereas at the end it is underestimated. Simulated temperature profiles at multiple planes are shown in Fig. 3. At the injector lips, thin flames develop, indicating a weak mixing between methane and oxygen as well as a slow heat exchange in the direction perpendicular to the main flow. Hence, even at the end of the combustion chamber, the temperature distribution is non-homogeneous. Shortly downstream of the faceplate, the flames in the outer row deviate from their initially round shape. In addition, their centers move outward. Downstream of approximately 0.1 m the outer flames merge, while the middle flame is not affected by any other flame.

Test Case Description and Computational Setup
The PennState preburner combustor is a frequently simulated combustion chamber, which has been investigated both with RANS/URANS [9,10,13] and LES [15,17]. It was experimentally examined at the Pennsylvania State University with the goal to provide data for the verification and validation of numerical codes [16].
The combustion chamber is axisymmetric and exhibits a single coaxial injector. The chamber length is 285.75 mm, its diameter 38.1 mm and the diameter of the throat 8.2 mm. Gaseous oxygen and gaseous hydrogen are preburned in an oxidizer and fuel preburner respectively and then are supplied to the combustion chamber with a O/F of ∼6.6. A more detailed description can be found in [16] or [13].
The experimental data set consists of wall temperature and wall heat fluxes. The measurements revealed an chamber pressure of 54.2 bar. However, Ivancic et al. [9,10] observed some inconsistencies in the data set, which can be explained by an incomplete preburner combustion.
The computational grid consists of 19.5 million volumes, divided into 39 blocks. The experimental wall temperature is used to prescribe the temperature at the combustion chamber wall. The injector walls as well as the faceplate are assumed to be adiabatic. At the inlet, a mass flow rate is specified. A supersonic outflow condition is imposed at the outlet. Additionally, the hydrogen oxidation kinetic reaction mechanism of Ó Conaire [18] is utilized which consists of 8 species and 19 reactions. In a previous 2D hybrid Lagrangian transported PDF simulation [6] it has been shown, see Fig. 4, that chemistry close to the injector is in a chemical non-equilibrium. These figure shows scatter plots of particle hydroxyl and oxygen mass fractions over temperature along a vertical line located closely behind the injector (x = 0.5 mm). The different symbols indicate different regions along that line. Both figures, and especially the one with the oxygen mass fraction, exhibit a strong scattering illustrating the occurrence of non-equilibrium effects in the near-faceplate region. Thus, in the following iDDES, finite-rate chemistry and the computationally more efficient assumed PDF approach are used to model combustion and turbulence-chemistryinteraction.
The iDDES is performed with a constant time step of t = 10 −8 s. The results are time-averaged over approximately 0.0112 s, which corresponds to around 1.35 flow-through times as defined by Tucker et al. [26]. This is a rather short averaging period. Nevertheless, the general trend in the results is already observable.

Results and Discussion of the iDDES
As the compresssible simulation includes the nozzle, the occurring pressure is a result of the simulation and depends on the wall heat transfer and the combustion process. A pressure of around 49 bar is reached. This value is significantly smaller than the Simulation Experiment nominal experimental pressure, but fits well to other simulations [10]. Figure 5 shows the wall heat flux profiles. As in the experiment, the simulation predicts a steep ascent of the heat flux shortly downstream of the faceplate. Besides, the location of the maximal wall heat flux fits. However, the maximal value is underpredicted by ∼25%. Furthermore, the decay of the heat flux is nearly constant in the simulation, which is in contrast to the experiment. Only further downstream, towards the end of the chamber, simulation and experiment agree. This deviation is difficult to explain as the wall heat flux depends on the temperature field as well as the thermal conductivity. These in turn are functions of the species composition close to the wall.  In the shear layer downstream of the injector lip, typical Kelvin-Helmholtz instabilities occur, leading to small-scale vortical structures. The time-averaged temperature field is similar to other high-fidelity simulations, such as the ones conducted by Oefelein [26] or Ma et al. [15]. For example, the flame shapes look alike. However, some differences appear, too. The temperature in the recirculation area is predicted differently, which could explain the deviation to the measured wall heat flux. Figure 7 shows instantaneous and time-averaged values for the calculated temperature variance σ T . This transported variable is required for the APDF approach. Close to the flame boundaries, and especially shortly behind the recessed post, σ T approaches high values. This is caused by the high temperature gradients in axial direction between the hot flame and the comparably cold injection jets.

Single Injector Combustion Chamber (BKC)
In order to analyze pressure oscillations, the combustion chamber BKC was experimentally investigated at the German Aerospace Center [24]. Hydrogen and oxygen are injected at cryogenic conditions in the axisymmetric combustion chamber. In addition, a hydrogen cooling film is injected near the wall. Thus, the global O/F is equal to 1. The length of the chamber is 430 mm, its diameter 50 mm and the diameter of the nozzle throat 16.8 mm. At the simulated operating point, the combustion chamber exhibits a weak oscillatory behavior with amplitudes of the pressure oscillations clearly below 1%. More details on the BKC can be found in [12,24]. The computational grid for the iDDES consists of around 10.1 million volumes. The constant time step size is set to t = 5 × 10 −9 s. Again, the hydrogen oxidation kinetic reaction mechanism of Ó Conaire [18] is used. The simulated mean pressure value reaches a value of 60.9 bar, which is slightly above the experimental one of 59.4 bar. To investigate the oscillatory behavior, Fig. 8 compares the discrete Fourier transform of the pressure signal at two different monitor points with the experiment. The first monitor point (MP1) is located near the faceplate in direct vicinity to the wall, the second one (MP2) is located further downstream. MP1 correctly predicts the frequencies of the first two longitudinal modes. However, the amplitudes are underestimated. In contrast, MP2 only reproduces the second mode. This is due to its location in the middle of the chamber where the fundamental mode is not existent. Higher modes are hardly resolved in the simulation. It should be noted that the extremely high amplitude for higher frequencies seems to be an artifact of the measurement technique [24]. Figure 9 shows experimental and simulated shadowgraphs of the near injector region. The latter can be determined by utilizing the second derivative of the density and performing a line of sight integration. Experiment and simulation exhibit a similar behavior. This includes the spreading of the hydrogen jet and the length of the undisturbed, cold oxygen jet.

Conclusions
An up to fifth order multi-dimensional limiting process with low diffusion has been described. Often, high-order schemes yield results that a low-order scheme can only deliver with a clearly refined grid. This MLP ld approach has been applied to three different rocket combustion chamber configurations with different turbulence model complexity. First, a steady-state RANS simulation of a seven element methane/oxygen combustion chamber has been carried out. The results show good agreement with experimental data. Second, the PennState preburner combustor has been simulated using iDDES. Although the time-averaged temperature distribution fits to other simulations, the comparison with the experiment reveals an underprediction of the wall heat flux. The time-accurate simulation of the BKC is in accordance with experimental data regarding pressure fluctuations as well as the flow field. These test cases highlight the necessity of high-order methods, especially in time-accurate simulations