Numerical Simulations of the Internal Ballistics of Parafﬁn–Oxygen Hybrid Rockets at Different Scales

: Hybrid rockets are considered a promising future propulsion alternative for speciﬁc appli-cations to solid or liquid rockets. In order to raise their technology readiness level, it is important to perform predictive numerical simulations of their internal ballistics. The objective of this work is to describe and validate a numerical approach based on Reynolds-averaged Navier–Stokes simulations with sub-models for ﬂuid–surface interaction, radiation, chemistry, and turbulence. Particular atten-tion is given to scale effects by considering two different parafﬁn–oxygen hybrid rocket engines and a simpliﬁed grain evolution approach from the initial to the ﬁnal port diameter. Moreover, a mild sensitivity of the computed regression rate to parafﬁn’s melting temperature, surface radiation emissivity, and Schmidt numbers is observed. Results highlight the increasing importance of radiation effects at larger scales and pressures. A numerical rebuilding of regression rate and pressure is obtained with simulations at the time-space-averaged port diameter, producing a reasonable agreement with the available experimental data, but a noticeable improvement is obtained by considering the grain evolution in time.


Introduction
Research in hybrid rocket engines is gaining momentum in recent years owing to experimental and numerical advances, and to an increasing number of test flights [1][2][3][4][5][6]. Hybrid rockets, which are propulsion devices burning a solid fuel and a gaseous or liquid oxidizer, are considered a promising alternative to liquid or solid propulsion for applications, including sounding rockets, space engines, auxiliary power units, and boosters [7,8]. They can be preferred to solid rockets owing to their higher specific impulse, increased safety, and throttling capabilities, and to liquid rockets because they are cheaper and simpler, and with higher average propellant density.
Paraffin-based hybrids are currently under investigation by the community due to their higher regression rate with respect to conventional fuels. However, there is still a relatively incomplete understanding of the relevant physical phenomena, including fluidsurface interaction, radiation, combustion, and turbulence, which occur simultaneously inside paraffin-based hybrid rockets. In fact, numerical simulations of the internal ballistics of paraffin-based hybrid rocket engines are becoming increasingly important to raise their technology readiness level. The most used numerical approach in the literature is based on Reynolds-averaged Navier-Stokes (RANS) simulations employing different sub-models [6], which either rely on a prescribed fuel mass flow rate [9][10][11][12][13] or employ a parametric gassurface interaction model [14,15]. Radiation effects on the regression rate, which can be important in HRE operating conditions [16][17][18][19], are also usually neglected when dealing with paraffin-based engines.
The objective of this work is to provide an additional step toward obtaining a predictive tool for the internal ballistics of paraffin-based HREs. To this end, a sensitivity study is performed on various aspects of the modeling strategy recently adopted and validated against experimental data of a lab-scale engine [4]. In particular, the effects of a change in paraffin's melting temperature, surface radiation emissivity, and laminar and turbulent Schmidt numbers are evaluated. Moreover, results on a different test case considering a larger scale are provided, highlighting scale effects on the internal ballistics. Finally, simulations at different diameters are analyzed in order to characterize convection and radiation heat flux changes obtained from a change in scale.
The manuscript starts with the description of the theoretical and numerical model in Section 2, followed by the illustration of the engine configuration and firing tests chosen for validation in Section 3. Then, results and discussion are shown in Section 4, and conclusions are summarized in Section 5.

Theoretical and Numerical Model
Numerical simulations were performed by solving the RANS equations for axisymmetric, turbulent, compressible, and reacting flows [20]. The specific isobaric heat and enthalpy dependencies on temperature were expressed, for all species apart from paraffin, with the seventh-order polynomials taken from the Chemical Equilibrium with Applications (CEA) database [21], which also provides the species standard heat of formation. The fourth-order polynomials of temperature for transport properties reported in [21] were employed for all species apart from paraffin, and Wilke's rule [20] was used to obtain mixture molecular transport properties. A constant Schmidt number equal to 0.7 imposed the same molecular diffusivity for all species. The standard Spalart-Allmaras model [22] was used to evaluate the turbulent viscosity, with turbulent Schmidt and Prandtl numbers equal to 0.7 and 0.9, respectively. Properties of paraffin will be discussed later.
In this work, a global reaction mechanism was employed to model the combustion of ethylene, which is the main product obtained from paraffin thermal cracking [23,24]. Due to the lack of detailed data in the literature, a reaction mechanism developed for butadiene combustion [25] was adapted here for use in ethylene combustion. Seven reactions and ten species were considered ( Table 1). The chemical source terms were obtained through Arrhenius-type forward reaction rates and backward rates calculated as the ratio between forwarding rates and the equilibrium constant evaluated from thermodynamic data taken from [26]. The reaction rate constant for the thermal cracking of C 32 H 66 was taken as the one for the liquid C 16 H 34 [24] due to the lack of data in the literature. More details on the reaction mechanism can be found in [4]. Table 1. Chemical reactions involved in the global reaction mechanism used for paraffinoxygen combustion.
The RANS simulations were performed with an in-house computational fluid dynamics (CFD) solver that has been validated in different operating conditions against experimental data [4,27]. The finite-volume computational tool is second-order accurate in space, and employs a Roe Riemann solver [28]. The Strang operator-splitting technique [29] for time integration was adopted: a second-order Runge-Kutta scheme integrates the convective and diffusive terms, whereas an implicit integrator for stiff ordinary differential equations is used for the chemical source terms [30]. Note that a local time-step approach was adopted, as steady-state solutions were sought. The code has been recently validated against experimental data of paraffin wax-oxygen hybrid rockets in [4,31].
The paraffin taken into account in the present manuscript is C 32 H 66 . Differently from conventional fuels, when paraffin wax is heated, it does not pyrolyze but rather melts, producing a liquid layer over the fuel grain. However, since the typical chamber pressure of hybrid rockets is higher than paraffin's critical pressure, equal to 6.5 bar [32], the melted paraffin is assumed to be at supercritical conditions, where no surface tension or boundary for droplets can be defined [33,34]. It is, therefore, reasonable to assume that in these conditions, the entrainment of the supercritical species is part of the turbulent mixing process. The melted species was modeled with the simplified dense fluid approach described in [4,31], with thermodynamic and transport properties taken from [35]. As described in [4], the fluid-surface interaction sub-model is based on mass and energy balances, which reduce to where q w,conv and q w,rad are the convective and radiative wall heat fluxes, respectively, and the paraffin density, melting enthalpy, specific heat, melting temperature, and initial temperature are, respectively, ρ s = 920 kg/m 3 , ∆h melt = 169.83 kJ/kg, c s = 1946.03 J/(kg · K), T melt = 343 K, and T s,in = 298.15 K. The radiative heat flux was computed with in-house software, which has already been described and used in [4,19,[36][37][38]. The main details of the numerical approach are reported below for the sake of completeness. The software solves the radiative transfer equation with the discrete transfer method for generic axisymmetric geometries, gray/diffuse boundaries, and inhomogeneous gray/nonscattering media. The evaluation of the radiative heat was carried out only at the boundaries given its small relevance, compared to the whole thermal power generated within the thrust chamber. Radiation from hydroxyls and from soot was not considered in the employed model. Absorption of radiative energy was assumed proportional to the pressure and to the absorption coefficients of H 2 O, CO 2 , and CO, which are considered as the major and only participating species to radiation, weighted with their molar fraction. A discretization consisting of 256 rays for each calculation point and a step of 1 mm along each ray were used after performing convergence analyses for both parameters. A wall emissivity equal to 0.91 was assumed for the paraffin wax grain by using the emissivity model proposed in [39] with a refractive index of 1.43 according to [40]. The CFD and radiation codes were coupled through the repeated evaluation of the radiative wall heat flux, then of the regression rate, and finally of the resulting flow field, until convergence was reached.

Engine Configuration and Firing Tests
The engine configuration for the simulations of this work is the axisymmetric thrust chamber depicted in Figure 1, which has already been used in [4]. It is composed by a pre-chamber (0 < x < x 0 ), a chamber with the paraffin grain (x 0 < x < x 1 ), a post-chamber (x 1 < x < x 2 ), and finally a converging-diverging nozzle (x > x 2 ). The nozzle conical converging and diverging sections are connected by circular arcs with each other and with the cylindrical post-chamber. The cylindrical pre-chamber, post-chamber, and fuel grain have all the same diameter. Note that, in general, the radius of the pre-chamber and post-chamber is different from that of the cylindrical grain port, which changes during the burn. However, this difference influences only slightly the numerical regression rate [41]; hence, it was omitted in the present study, and a simplified geometry with a single radius (equal to that of the cylindrical grain port) was considered. The computational setup is completed by the boundary conditions outlined below. On the left-hand side of the setup, an inlet boundary condition imposing gaseous oxygen mass flow rate, uniform static temperature equal to 300 K, and uniform turbulent viscosity equal to 10 −6 Pa · s is set up for 0 ≤ r ≤ r inj , while an adiabatic wall is set up for r inj < r ≤ R. The influence of the inlet turbulence viscosity on the flowfield is assessed as minimal. On the top side, Equation (1) is computed on the paraffin grain section, and adiabatic walls are imposed in the pre-chamber, post-chamber, and nozzle. Finally, a supersonic outlet boundary condition is considered on the right-hand side. Depending on the test conditions, a first simulation is carried out from an initial uniform flowfield with oxygen at 12 bar and 3000 K until steady state, which serves as the initial condition for the following simulations.
Aiming to validate the predictive capabilities of the present approach, two sets of literature experimental firing tests were considered (see Table 2). Seven test cases were chosen from the lab-scale experimental tests carried out at the University of Naples "Federico II" [13] (named set 1), while five were from experiments performed at the hybrid combustion facility of Nasa Ames [42] (named set 2). Oxidizer mass flow rates ranged from 29 to 60 g/s for set 1, and from 2.05 to 4.4 kg/s for set 2. The injection pattern of the experiments in both sets was correctly replicated by the abovementioned inlet boundary condition, as it was composed of a single circular injector. The experimental setups were of different sizes, the latter being quite larger than the former. For all tests of set 1, A single simulation at the time-space-averaged port diameter was performed for each test case. In addition, simulations at the initial and final diameters were also performed for tests 3 and 4 of set 1. Table 2. Setup conditions for the two hybrid engines considered in this study. Set 1 is taken from [13], while set 2, from [42].

Set
Testṁ ox (g/s) R (mm) r t (mm) Computational meshes were created considering a proper grid clustering in the injection, near-wall, and throat regions to sufficiently resolve the mixing layer, the boundary layer, and the transition through sonic conditions, respectively, (see Figure 2). The mesh employed for the numerical simulations of tests in set 1 was composed of 160 grid cells in the axial direction and 100 grid cells in the radial direction. The resolution employed yielded a maximum cell height on the grain of 12.34 µm, and a corresponding dimensionless wall distance for wall-bounded flows y + of the order of 1. Additional details and adequate mesh convergence studies have been reported previously in [4]. For simulations of tests in set 2, a topologically similar mesh of 600 × 320 elements ensured approximately the same wall resolution as for set 1.  Table 2).

Results
The results from the Navier-Stokes simulations show several peculiarities of hybrid rockets burning paraffin wax and gaseous oxygen. We start with a discussion of the results from the small-scale engine (set 1), followed by the analysis of the firing tests of the larger engine (set 2), concluding with the numerical rebuilding of experimental regression rate and chamber pressure.

Results on Set 1 Tests
The tests in set 1 are characterized by experimental oxidizer-to-fuel ratios from 0.77 to 1.26; hence, they are all in fuel-rich conditions (the stoichiometric O/F is equal to 3.44). Figures 3 and 4 show numerical results for test 4 of set 1 obtained by increasing the port diameter from the initial to the final one with the average diameter as an intermediate step. This was performed to simulate the time evolution during a firing test. A common feature for all diameters is the presence of a wide diffusive flame developing from the pre-chamber throughout the whole combustion chamber. The inhomogeneous presence of reactants and combustion products is observed at the nozzle entrance, where combustion still takes place.  At the initial diameter, the recirculation zone is mostly confined upstream of the fuel grain, which entails complete combustion of the cold liquefied paraffin before reaching the injector plate; hence, hot combustion gases are filling the pre-chamber. On the other hand, as the port diameter grows, the vortex penetrates towards the grain and brings upstream more paraffin, which cools down the pre-chamber and allows the flame to be anchored at the interface between the oxygen injector and the injector plate wall. The amount of paraffin present in the pre-chamber is, in fact, strongly dependent on the extension of the vortex, and post-firing inspections of the injector plate confirmed the presence of unburnt paraffin [13]. Downstream of the reattachment point, appearing at varying axial locations on the fuel grain depending on the port diameter, the liquid paraffin accumulates with increasing axial distance, covering as a thick film the post-chamber and partially also the nozzle walls. At the final diameter, a lower mass fraction of both oxygen and paraffin is observed at the nozzle entrance. Moreover, the flame reaches the symmetry axis earlier with respect to the other two smaller diameters, yielding a more complete mixing and combustion in the core region. The mass fraction of molecular oxygen reduces along the length of the engine because of combustion; however, part of the injected oxygen directly flows away from the nozzle at the average diameter, indicating a lower combustion efficiency since this test is very fuel-rich (O/F = 1.3). Combustion products such as water vapor or carbon dioxide cannot reach the walls due to a thick layer made mainly of paraffin and carbon monoxide. This actually also prevents nozzle erosion, which was not observed in the experiments [13].
The regression rate is greatly influenced by the vortex extension over the grain ( Figure 5): As the vortex reattachment point moves downstream, the regression rate shows a corresponding local minimum or change in slope. Furthermore, due to the reduction of the total mass flux, the overall magnitude of the regression rate tends to decrease with increasing port diameter. Moreover, pressure axial nonuniformity vanishes with increasing port diameter, due to the reduced Mach number in the port caused by the larger contraction ratio (see Figure 6).

Model Sensitivity Analysis
In order to establish the effects of model parameters/constants on the computed regression rate and pressure, in this section, we consider only the simulation of test 4 of set 1 at its average diameter.
An increase in the melting temperature of paraffin T melt from the reference value of 343 K is expected to reduce the convective wall heat transfer, yielding a lower regression rate [14]. This is indeed observed by comparing numerical simulations at three different values of T melt (Figure 7). The sensitivity of the regression rate is shown to be −5% every 15 K. The chamber pressure also decreases for higher values of T melt , at a rate of approximately −0.65% every 15 K. The pressure rate of change is lower than the one of the regression rate because a higher characteristic velocity balances the lower total mass flow rate. Regression rate and pressure are unchanged for variations of the paraffin grain emissivity ε w (which equals absorptivity due to the gray wall assumption) in the range 0.91-1.0 ( Figure 8). This leads to the possibility of treating the paraffin grain as a black body for this fuel-rich case and shows that the results are not sensitive to increased grain wall emissivities in the range 0.91-1.0. As the turbulent Schmidt number is reduced, mixing is enhanced due to the increased turbulent diffusivity; hence, higher regression rate and pressure are obtained (Figure 9). Different turbulent Schmidt numbers entail 1-3% changes in average regression rates and chamber pressures; however, it is worth noting that typical values of turbulent Schmidt numbers used in the literature are in the range 0.7-0.9. On the other hand, regression rate and pressure are insensitive to changes in the laminar Schmidt number (Figure 10). In fact, molecular diffusion phenomena are not expected to play a significant role in this fuel-rich test case, as next to the wall, there is mostly paraffin. The small sensitivity of regression rate and pressure to the laminar and turbulent Schmidt numbers justifies the adoption of Fick's law and the assumption of constant Schmidt numbers in this work and allows us to avoid using other multicomponent diffusion models.

Results on Set 2 Tests and Effect of Radiation
The tests in set 2 are characterized by a larger scale with respect to set 1, by experimental oxidizer-to-fuel ratios from 1.70 to 2.69 (whereas in set 1, O/F ranges from 0.77 to 1.26), and by a reduced post-chamber length (for test L01 (x 2 − x 1 )/R = 1.1, while for test 4 of set 1 (x 2 − x 1 )/R = 4.6). Nevertheless, the same qualitative internal combustion phenomena of set 1 are observed: the axial injection mode yields flame anchoring at the interface between the oxygen injector and the injector plate wall with a recirculation zone in the pre-chamber, which is filled by unburnt paraffin (Figure 11). Due to the short post-chamber, mixing is not favored, and the flame does not reach the centerline. In addition, a strong CO presence is observed above the flame but not close to the wall, which is protected by a thick layer of cold paraffin up to the nozzle walls. Moreover, a large fraction of the oxygen does not burn but is just accelerated through the nozzle, decreasing significantly the efficiency of this fuel-rich regime. It is interesting to note that ethylene mass fraction is mostly null, indicating that when the cracking reaction of paraffin is triggered, the ethylene produced is burnt right away and does not accumulate (the same holds for set 1 results).
During the simulations of all tests in set 2, mild numerical instabilities involving the whole flow field are observed. This is due to either numerical or physical reasons. In fact, numerical results showed oscillations around a constant average flow field, which is taken as reference in the present work, and the deviations from it are considered as numerical error. The physical nature of these oscillations should be confirmed by time-accurate simulations, which are, however, outside of the scope of this work because of their too high computational effort. Note that in the firing tests of set 2, experimental pressure oscillations between 4 and 12% of the mean pressure are observed [42]. Their peak in the pressure spectrum is at 30, 100, and 350 Hz, corresponding to boundary-layer interaction with solid thermal lag (low-frequency hybrid instability), bulk mode, and acoustic half-wave in the combustion chamber. Such oscillations are not observed either physically or numerically in set 1. A significant effect obtained by increasing the size and the operating pressure of the engine from set 1 to set 2 is the higher relevance of radiation in the wall heat flux balance, which directly controls the paraffin's regression rate (Table 3). For set 2, this contribution reaches values as high as 88%, whereas for set 1, it is at a maximum of 62%. In fact, the radiation contribution is found to be proportional to the product between chamber pressure and port radius, except for the L04 case, which is the one characterized by the highest port mass flux, yielding an increased share of convection heat transfer. These results depend on the employed models for turbulence, combustion, and radiation, and need to be verified with experimental data on radiation heat transfer in paraffin-based hybrid rockets, which are lacking in the literature.  Table 2). Table 3. Wall heat flux contributions on the grain surface for all tests at the average port diameter.

Set
Test p c R (bar·m) q w,tot (MW/m 2 ) q w,rad /q w,tot

Numerical Rebuilding
The integral average of the axial profiles of the regression rates and the post-chamber pressures (x = 0.28 m) are extracted from the numerical simulations at the average diameter and are compared to the respective experimental values in Figure 12. Very close agreement is found for chamber pressures lower than 20 bar, whereas an error of about 15% is observed for the higher chamber pressures of set 2. Average regression rates are generally in good agreement with experimental results, with maximum errors of 15%. Table 4 shows a comparison between numerical and experimental O/F, pressures, and combustion efficiencies for all tests. The combustion efficiencies are underestimated by the numerical simulations, only slightly (by 1-5%) for tests in set 1 and by larger amounts (6-22 %) for tests in set 2, also due to the larger errors in the pressure rebuilding.
Overall, the model, which does not require information from any firing test, shows an acceptable prediction of the experimental data of the two engines at different scales without any dedicated tuning. This is achieved by making simplifying assumptions that do not make the model computationally expensive and do not alter significantly its predictive capabilities. Table 4. Experimental and numerical oxidizer-to-fuel ratio, chamber pressure, characteristic velocity, and combustion efficiency for all tests in Table 2.   In order to further validate the computational model and improve numerical predictions, it is important to verify that results at evolving diameters are coherent with the experimental data. To this purpose, simulations of test 3 of set 1, for which pressure probe data in time are available from [13], were performed at the initial, average, and final port diameters. Flow fields of test 3 are not shown because they are very similar to the ones already shown for test 4 of set 1 (Figures 3, 4, and 6). As previously observed for test 4, increasing diameters result in decreasing regression rates, due to lower mass fluxes in the port, and more axially uniform pressure profiles due to lower Mach numbers ( Figure 13). The O/F resulting from the simulations at the different diameters are 0.90, 1.22, and 1.24, showing that test 3 is in fuel-rich conditions for the whole burn duration, which is consistent with the absence of throat erosion observed in the experiments. In addition, the characteristic velocities obtained are 1218 m/s, 1388 m/s, and 1501 m/s, from the lowest to the highest diameter, respectively, corresponding to c * -efficiencies of 86%, 85%, and 91%. It is interesting to observe how the regression rate is influenced by the different contributions of convection and radiation. At the initial diameter, the dominant share of the wall heat flux is provided by the convective heat flux, with radiation actually helping in making the total wall heat flux more uniform in the axial direction ( Figure 14). However, the radiation wall heat flux becomes more significant by increasing the engine diameter (Table 5). This is quite an interesting result as, despite the significant change of the share of radiative/convective heat flux, the numerical pressure is rather constant, in line with the experimental measurement ( Figure 15). In addition, it is noted that the radiation contribution starts to become non-negligible for oxidizer mass fluxes lower than G ox of the order of 100 kg/(m 2 s), which is consistent with the value of 140 kg/(m 2 s) obtained by [7]. This gives further confidence to the numerical models used in this work. q w,conv D 0 q w,conv D ave q w,conv D final q w,rad D 0 q w,rad D ave q w,rad D final Figure 14. Effect of port diameter on wall heat flux for test 3 of set 1.  Numerical times in Figure 15 are obtained by considering the average between two successive regression rates and diameters and were axially shifted by 1 s to match the experimental ignition transient. By considering just one simulation at the average port diameter, the average experimental regression rate of test 3 is underestimated by 11.8%. On the other hand, the integral average of the numerical regression rate obtained with the three successive simulations yields an underestimation of just 2.2%, due to the higher regression rate obtained at the initial diameter. This underlines the need to numerically simulate the evolution of the grain in time for more accurate numerical predictions.

Conclusions
A numerical approach based on Reynolds-averaged Navier-Stokes simulations including fluid-surface interaction, radiation, combustion, and turbulence, was used to analyze hybrid rocket engines at different scales. Due to the uncertainty of some of the model parameters, a sensitivity analysis was carried out for paraffin's melting temperature, surface radiation emissivity, and Schmidt numbers, showing a negligible or mild effect on the computed regression rate and pressure. Although no tuning or calibration was carried out, the results compare well with experimental data, and, owing to appropriate assumptions, the model retains its predictive capabilities without compromising the simplicity and inclusiveness of the model strategy.
A quite interesting result is that numerical predictions are in line with experimental measurements despite the wide range of the share of convective and radiative heat flux. This provides confidence that both convection and radiation models are well suited to numerically simulate the internal ballistics of paraffin-based hybrid rocket engines. In particular, the share of the radiative over total wall heat flux ranges from 10% to 88%, depending on the chamber pressure, port radius, and mass flux. As a further demonstration of the capability of the approach, it is found that by evolving the port diameter, the nonlinear time evolution of the experimental regression rate is correctly captured.