A Reduced Order Model based on Large Eddy Simulation of Turbulent Combustion in the Hybrid Rocket Engine

A combined method of large eddy simulations for non-premixed combustion in a turbulent boundary layer coupled with proper orthogonal decomposition of instantaneous velocity, pressure and temperature fields is developed in order to obtain a reduced order model. First, we investigate a channel turbulent reacting flow using Large Eddy Simulations (LES) technique. Polypropylene/O2 has been considered as fuel/oxidant pair. The turbulence-combustion interaction is based on a combination of finite rate/eddy dissipation model applied to a reduced chemical mechanism with four reactions. The LES numerical results are analyzed with respect to RANS simulations and with other reference data. The second part of the paper refers to the derivation of a Reduced Order Model (ROM) based on proper orthogonal decomposition (POD) technique for the unsteady flow field. In order to achieve that, the eigenmodes of the flow are computed from several snapshots of the instantaneous fields uniformly spaced and the most energetic ones are used to set up the Reduced Order Model. Constant regression rate of the fuel grain is considered. The flow and thermal fields obtained with ROMs are compared with the ones obtained from the full simulation and an analysis on the number of modes required to achieve the desired accuracy is presented.


Introduction
The concept of hybrid rocket propulsion is not new and it was first introduced around 1937 in Russia by Andrussow according to Green [1]. In the early stages of rocketry research, the solid rocket fuel systems were extremely used because of their major advantages: simple to design, high thrust and low cost of the motor. On the other hand, the solid rocket booster cannot be throttled, cannot be shut down, can operate for a relatively short time and present a potential risk of explosion. This causes an increased cost as it requires special manufacturing, handling and using conditions and many restrictions for the launch campaign. Starting with 1960's, more and more attention is paid to hybrid rocket propulsion which fill the gap between solid and liquid systems, sharing the advantages from both configurations. They are relatively cheaper in design and manufacture; they produce a higher specific impulse than solid rocket motors and higher density-specific impulse than liquid bipropellant engines and smoothly change motor thrust over a wide range on demand. What is extremely important, hybrid engines present a negligible risk of explosion or detonation. There are three major disadvantages of hybrid rocket engines: the low regression rate, the lower densityspecific impulse than solid propellant systems and the casing around the fuel grain must be built to withstand to full combustion pressure and often to extreme temperatures as well. In the hybrid combustion case, however, nonuniformity of flow conditions is the principal characteristic. Not only does mass addition accelerate the flow, as in a solid-propellant rocket, but the streamwise variation in chemical composition and temperature also is considerable [2].
Typical fuels for a hybrid rocket are the polymeric synthetic rubbers based on the polybutadiene monomer, (PB with the formula C 4 H 6 ). The most popular of this group, based on cost and commercial availability, is HTPB (hydroxyl-terminated polybutadiene) with possess an energy density comparable to kerosene. Other hydrocarbons that have been used, mostly in smaller motors, are the paraffin waxes, polyethylene, plexiglass (Lucite), polypropylene or meta toluene diamine/nylon. The formulated fuels, PB polymers or paraffin waxes, have the virtue of allowing performance additives, like Al, AlH 3 , Li, LiH, Li 3 AlH 6 , B, B 10 H 14 . Commonly used liquid oxidizers are liquid oxygen (LO 2 or LOX), liquid fluorine (LF 2 ), nitrogen tetroxide (N 2 O 4 ), and nitric acid (HNO 3 ). In terms of propulsion performance, the theoretical vacuum specific impulse performance of LOX/HTPB hybrids exceeds 360 s. This is considerably higher than the best of all solid-propellant rockets (about 320 s) [3].
The main drawback of the hybrid rocket is that the combustion process relies on a relatively slow mechanism of fuel melting, evaporation and diffusive mixing. As a rough comparison, the regression rate in a solid rocket at a typical rocket combustion chamber pressure may be on the order of 1.0 cm/s whereas for a typical hybrid using a classical polymeric fuel may have a regression rate on the order of 0.1 cm/sec. To compensate for the low regression rate, a simple idea is to increase the surface area for burning using a multiport fuel grain. Most attempts to increase the regression rate involve methods for increasing the heat transfer rate to the fuel surface such as increasing turbulence levels in the port or adding roughness to the fuel grain. The regression rate determines the overall sizing, mass fluxes, and geometric configuration of the hybrid fuel ports. In classical hybrid rocket motors, liquid or gaseous oxygen flowing over the solid fuel reacts with the pyrolyzed gases close to the fuel surface and forms a turbulent diffusion flame. Convective and radiative heat transfer from the flame, provide the heat of pyrolysis for the thermal decomposition of the solid fuel. It is evident that the regression rate is a function of both the axial position in the grain and the stage in the burn.

Mathematical model
In the course of its operation, a hybrid motor undergoes several transient phenomena. Some of these transients, such as ignition and thrust termination, are inevitable, some are imposed by mission requirements, such as throttling, and the remainder are undesirable events, such as sudden pressure spikes or instabilities. Understanding the transient phenomena that take place in a motor is critical for the designing a propulsion system that can deliver the desired performance within the tolerance limits.
The fully Navier-Stokes equations describing the conservation of mass, momentum, total energy and conservation of N chemical species are [4]: In the above equations, m Y is the species mass fraction of the i-th species, , i m V is the diffusion velocity of the j-th species in the i-th direction. The mass reaction rate per unit volume of the j-th species is denoted by m ω  . The heat flux vector contains the thermal conduction, enthalpy diffusion (i.e. diffusion of heat due to species diffusion), the Dufour heat flux and the radiation heat flux. Dufour heat flux and radiation heat flux are neglected.

Large Eddy Simulation (LES) formulation
In turbulent reacting flows, the instantaneous range of velocity and thermal scales increases rapidly with turbulence intensity. It is obvious that for the DNS approach, the grid should have a smaller resolution, both, spatial and temporal, smaller than the physical scales. Generally, in engineering flows, approaches based on the Reynolds-averaged Navier-Stokes (RANS) equations are the most prevalent and involve computing one-point moments such as mean velocity and turbulent kinetic energy. The greatest downside of RANS solutions is that they drastically limit the amount of information provided by numerical simulation.
Since DNS is currently out of reach even for the most advanced computers in use, the only approach that is both reliable in terms of accuracy and feasible in terms of computational cost remains the LES. LES resolves both the large, geometry dependent turbulent scales (as RANS does as well) and a fraction of the smaller energy containing scales within the inertial range, up to a level dictated by the resolution of the numerical grid, and only the remaining scales are modelled. If the grid resolution is appropriately chosen, the unresolved scales, by Kolmogorov's hypothesis are isotropic and, therefore, more amenable to modelling [5]. The closure of subgrid terms is a major area of research and many approaches have been proposed in the past. In the present paper, the Smagorinsky model has been chosen for assessing the subgrid closure. The effects of the turbulence are generally advantageous for the efficiency of the combustion, since turbulence enhances the mixing of component chemical species and heat. Generic combustion models can be divided into three categories: topological model, reactor models and models with finite reaction rate. In our work the last class was used. This class of models calculates the rate of chemical reaction based on a reaction where, however, occur filtered values of temperature, pressure and molar fractions of mass and species. The finite rate models are versions of classic models, in the sense that it is considered that the process that controls the release of heat is the lowest response rate of turbulent mixing and chemical reaction. Thus, it is calculated, in addition to turbulent mixing rate and the rate of chemical reaction and the filtered reaction rate is determined as the lower of two rates.

Flow configuration
The problem studied is similar to that presented by G. Gariani, F. Maggi and L. Galfetti, considered as a simplified model of a hybrid rocket engine (Fig. 1) [6].

Fig. 1 Flow configuration
The case is related to the 2D hybrid combustion of polypropylene (C3H6) with nitrous oxide (N2O). The adopted mesh is a structured one having a constant mesh size in both space directions and containing 1562000 cells. The constant time step 5 10 t s − ∆ = s is imposed for the unsteady calculations. The imposed boundary conditions were: at the inlet section, the oxidant (considered pure oxygen) has an bulk velocity (V inlet = 10 m/s) and a constant temperature T= 300 K. On lower wall of the combustion chamber, a constant fuel injection ( 0.01 inj V = m/s) is considered. Finally, at the exit section, the second order extrapolation for all variables conditions has been imposed. The kinetic mechanism of propylene combustion has been adapted from Farbar et al. [7], Marzouk and Huckaby [8] and Westbrook and Dryer [9]. The turbulence-chemistry interaction has been provided through finite-rate/ED (eddy dissipation) mechanism. The inlet boundary condition for LES calculations were obtained from the RANS steady state solution (using standard − turbulence model with nonequilibrium wall functions) [10]. During this preliminary run, uniform inlet velocity for oxidant (Vinlet=10 m/s), uniform injection velocity for C3H6 and turbulence intensity, T u =10%, have been chosen. The operating pressure was 1 atm. The inlet mass fraction of O2 was 0.36 (as resulting from nitrous oxide (N2O) thermal decomposition), while along the fuel slab a constant mass fraction of 1 has been imposed for C3H6 and constant temperature of 1000K. The temperature of the oxidant stream was 300K. Lower solid walls were considered adiabatically isolated and for the upper walls a constant temperature of 900 K was imposed.

Results
To initiate the combustion, a spark ignition at 5 2 10 t − = ⋅ s was assumed. Due the very small chemistry time scales involved by the combustion, the flame develops very rapidly and, together with the blowing dominates the flow development in the channel (Fig. 2a). Anchored at the leading edge of the blowing wall the flame is lifted by the injection. This affects the boundary development inducing a sudden growth in the streamwise velocity (Fig. 2b). In this zone, in the middle of the channel the streamwise velocity decreases. This inclusion with strong velocity gradient along the channel width is then transported downstream and as Fig.  2a shows the streamwise velocity defect in the cannel core disappears. Fig. 2 a) Instantaneous streamwise velocity, temperature, vorticity magnitude and sound sources at t=0.001s; b) Instantaneous streamwise velocity, temperature, vorticity magnitude and sound sources at t=0.2s

a) b)
As previously mentioned, the sound sources can be related with the velocity and vorticity fields interaction (Fig. 2a-b). For large time values, far downstream, this interaction affects almost the entire channel width. Finally, it is clearly that the constant wall blowing induces a non-uniform pressure gradient in the normal to wall direction. Moreover, the Proper Orthogonal Decomposition analysis will demonstrate the persistence of the time pressure fluctuation, even the temperature and velocity fields reach a quasi-steady state.

Proper Orthogonal Decomposition Analysis
Historically the method of Continuous POD (or the classical method) of Lumley [11] proceeded by the Snapshot POD of Sirovich [12]. More information regarding the application of the proper orthogonal decomposition in the analysis of turbulent flows together with a detailed bibliography is given in Berkooz et al. paper [13]. In this paper, we used the Snapshot POD because it is much more efficient from the numerical point of view.
The POD is a method that reconstructs a data set from its projection onto an optimal base. Besides using an optimal base for reconstructing the data, the POD does not use any prior knowledge of the data set. It is because of this that the basis is only data dependent, and this is reason that the POD is used also in analysing the natural patterns of the flow field.
For the reconstruction of the dynamic behaviour of a system the POD decomposes the data set in two parts: a time dependent part, ak(t), that forms the orthonormal amplitude coefficients and a space dependent part, ψ k (x), that forms the orthonormal basis. For the variable ( , ), the reconstructed data set is: where M is the number of time instant observations in the data set. We denote the error of the reconstructed data set as: The base from which the data set is reconstructed is said to be optimal in the sense that the average least squares truncation error is minimized for any given number ( m M ≤ ) of basis functions over all possible sets of orthogonal functions: = ⟨( , )⟩ (7) where the ⟨⋅⟩is the ensemble average and (⋅)is the standard Euclidian inner product. It was shown that the minimization condition for error ( , ) translates into maximum condition for: This maximization can be proven to take place if the time independent base functions ψ(x) are obtained from the Fredholm integral equation: where R ij is the correlation kernel. In this way, we transform this into an eigenvalue problem and λ k is the eigenvalue corresponding of the eigenvector ψ k . Because we can consider the inner product as being the equivalent of an "energy", the value of λ k is linked to the energy contained in mode ψ k and the optimization process involved can be summarized as: the data set is projected onto a basis that maximizes the energy content. While in the classical approach of Lumley [11], the correlation matrix is constructed as a space correlation matrix and solving the eigenvalue problem, we obtain directly the eigenvectors as the spatial modes and then use them in order to obtain the time-dependent coefficients ( ) = ( ( , ), ( )) (10) in the Snapshot POD of Sirovich [12], the correlation matrix is a time correlation matrix: which is of the size of the square of the number of snapshots. From the time correlation matrix, we get the eigenvalues λ k and time dependent eigenvectors φ k (t). The spatial eigenmodes that are time independent, are computed according to the formula: where = � (13) For the reconstruction of u(x,t), we take into account only a small number of modes that contain the most energy . (14) The processed data are the variations of streamwise velocity magnitude, temperature static pressure fields. Fig. 3 a) Fraction of total energy for the most energetic modes for pressure, temperature and streamwise velocity; b) The time-dependent coefficients a k (t) corresponding to the first four most energetic modes of pressure field a) b) Fig. 4 a) The time-dependent coefficients a k (t) corresponding to the first four most energetic modes of temperature field; b) The time-dependent coefficients a k (t) corresponding to the first four most energetic modes of streamwise velocity field These variations were obtained from numerical simulations presented in the previous paragraph. We took 25 snapshots and the time between adjacent snapshots is of ∆t = 0.001s; therefore, the Snapshot POD of Sirovich yields 25 eigenmodes for each considered field.

a) b)
The very high efficiency of the proper orthogonal decomposition is clearly underlined in Fig. 4a. The first mode contains more that 94% of the total "energy" for all considered variables. For the pressure and temperature fields the second mode contains only 0.38% and 3.37% of total energy. Consequently, the temperature field is more sensitive to high order perturbations.
The analysis of time-dependent coefficients a k (t) (also called modal amplitudes or Fourier coefficients) allows to see if the neighbouring spatial modes could interact or the interaction among them is excluded. Comparing Fig. 3b and Figs. 4a and 4b it is obviously that the velocity and temperature fields tend to a quasi-steady state. However, the pressure field remains unsteady (the first mode) and very sensitive to high order perturbations represented by next 3 modes. This behaviour can be underlined also analysing the modes itself [14]. It is clearly that, in order to reconstruct the pressure field, we need to combine at least two more pressure modes since their contributions reflected by the values of ak(t) are very close. In conclusion, the pressure fluctuations cannot be directly correlated with streamwise velocity field or with temperature field. Thus, we presume the wall injection velocity as responsible to strong pressure fluctuations. We recall that in practical hybrid combustion the "injection" velocity is the result of the regression rate, which is modelled as function of the local pressure value. Appears that the local regression rate can strongly vary in time and space due the strong pressure changes.
The much lower variation of the temperature field is probably due, on one hand to the very small-time scales of chemical reactions and, on the other hand to the assumed chemical kinetic mechanism. In first approximation, the temperature is determined by the flame development, depending mainly on the oxidizer mass flux.
Considering the above results, we suggest investigating the pressure fluctuations evolution by increasing the interaction between streamwise velocity field and the normal to wall velocity field induced primarily by the wall injection. Because the interactions are very complex and extremely hard to understand we propose to increase the global turbulence intensity of the flow. However, the supplementary turbulent structures injected in the main flow must keep the low spatial and temporal scales in order to avoid strong interaction with the coherent high energetic modes.

Conclusions
The paper investigates the turbulent combustion in flows configurations compatible with hybrid rocket engine. Polypropylene (C3H6) with nitrous oxide (N2O) has been considered as fuel/oxidant pair. The turbulence-combustion interaction is based on a combination of finite rate/eddy dissipation model applied to a reduced chemical mechanism with four reactions. Using POD strategy, the series of complex interdependencies between velocity, thermal and pressure fields are outlined. In order to reduce the pressure oscillations a static control technique, is analysed. Two grids inserted in pre-chamber are designed to generate vortices at different scales, increasing the global mixing. Finally, the POD is used to identify the turbulence effects introduced by the fine grid.