Numerical Simulations of Fuel Shape Change and Swirling Flows in Paraffin/Oxygen Hybrid Rocket Engines

The objective of this work is to describe and validate a numerical axisymmetric approach for the simulation of hybrid rocket engines (HREs), based on Reynolds-averaged Navier–Stokes simulations, with sub-models for fluid–surface interaction, radiation, chemistry, and turbulence. Fuel grain consumption is considered on both radial and axial directions and both axial and swirl injection of the oxidizer are simulated. Firing tests of two different paraffin–oxygen hybrid rockets are considered. First, a numerical rebuilding of fuel grain profile, regression rate and pressure for axial-injected HREs is performed, yielding a reasonable agreement with the available experimental data. Then, the same numerical model is applied to swirl-injected HREs and employed to analyze both the flowfield and the regression rate variation with swirl intensity. A validation of the model through the rebuilding of small-scale firing tests is also performed.


Introduction
Hybrid rocket engines (HREs) are propulsion devices usually employing a solid fuel and a gaseous or liquid oxidizer, stored physically separated from each other. If compared to solid rocket engines, HREs are safer during fabrication, storage, and operations. They allow throttling, shutdown, and restart capabilities, and they present less ambient temperature sensitivity, higher crack robustness of fuel grain, and higher specific impulse. On the other hand, with respect to liquid rocket engines, they are much simpler and cheaper to build, more reliable, and have higher average propellant density. HREs are therefore considered one of the potentially preferred options for specific future-generation propulsion systems [1] as they have already shown some promising results with the successful flight of SpaceShipOne [2] and SpaceShipTwo, and with the ongoing work on the SL1 launcher by Hympulse [3]. However, one of the most important drawbacks of conventional HREs using classical pyrolyzing fuels such as hydroxyl-terminated poly-butadiene is the low regression rate of the grain, which entails low thrust levels especially in combination with high-performance oxidizers, such as oxygen, which require lower O/F for maximum efficiency. This shortcoming can be mitigated with the application of different techniques, such as multiport grains or the introduction of fuel additives. However, the most promising one to date is the use of paraffin-based fuels. In fact, contrary to conventional polymeric fuels, which pyrolyze before burning, paraffin-based grains in HREs exhibit a liquid or supercritical fluid layer, depending on the operating pressure, leading to the entrainment of droplets into the gaseous mixture stream [4,5]. This mechanism allows for a continuous spray of fuel along the port, leading to an additional mass transfer toward the melt layer and the flame front, with most of the fuel vaporization occurring around the droplets. Regression rates up to three to four times higher than the conventional values were first observed in laboratory-scale motors and then confirmed in scale-up tests with different oxidizers [6,7]. High regression rates allow one to design high-volumetric-loading single-port combustion chambers, avoiding complex and inefficient multiport grains. In addition, paraffin-based fuels are nonhazardous, nontoxic, and easy to handle.
The performance of paraffin-based HREs may be further improved by the use of swirl injection. In fact, the resulting vortex motion increases the convective heat flux to the grain [8] and the regression rate, and improves the mixing between the propellants, which may lead to an increase in combustion efficiency. Several firing tests have shown how such injection results in a regression rate increase up to 5 times [9][10][11] with respect to axial injection.
CFD analysis of HREs with swirl injection has been the subject of many works in the past years, mostly employing 3D [12,13] or LES [14] simulations, wich have high computational cost, or making use of simplified injection approaches [15]. Most works also neglect the effects of radiation, which is instead known to give an important contribution to the total heat flux.
In this work, RANS simulations of HREs burning gaseous oxygen and paraffin-wax, employing both axial and swirl injection, are performed. All simulations use an axisymmetric approach, to reduce the computational time and to be used as a quick design tool.
The manuscript is organized as follows: first, the numerical model is briefly described (Sect. 2), then it is validated through the rebuilding of firing tests with axial injection (Sect. 3). In Sect. 4 the main features of hybrid rockets with swirl injection are investigated and the regression rate prediction capabilities of the numerical approach are validated though the rebuilding of small-scale firing tests.

Numerical Model
The numerical simulations have been performed solving the Reynolds-averaged Navier-Stokes (RANS) equations [16], with submodels accounting for the effects of turbulence, chemistry, gas-surface interaction and thermal radiation. In all simulations, an axisymmetric approach with a periodic boundary condition ensuring the tangential symmetry of the flow has been employed. The numerical model is extensively described in [17] and some of the main details are reported below for the sake of completeness. The adopted turbulence model is Spalart-Allmaras, which allows for a reduced computational cost with respect to other RANS turbulence models. A correct modeling of turbulence is needed to correctly assess the wall heat flux and the mixing of the propellants. In the case of swirl injection, conventional isotropic RANS models are known to fail to correctly predict the anisotropic characteristics of turbulence [18]. However, works such as [19] have shown that isotropic models yield satisfactory results near the wall, thus ensuring a good prediction of the fuel regression rate. The chemical reactions are solved with a finite-rate approach, employing the global chemical model shown in [17]. The interaction between turbulence and chemistry is neglected, which can lead to underestimation of the combustion efficiency and of chamber pressure. This is however deemed to yield a model uncertainty comparable to other model assumptions employed.
As described in [17], the fluid-surface interaction submodel regarding pure paraffin 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, ̇r the fuel regression rate and the paraffin melting enthalpy, specific heat, melting temperature, and initial temperature are, respectively, Δ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 paraffin density s is 929 kg/m 3 for the setup analyzed in Sect. 3 and 870 kg/m 3 for the one of Sect. 4.
The radiative heat flux is computed with in-house software, already used in [17,[20][21][22], which solves the radiative transfer equation with the discrete transfer method for generic axisymmetric geometries, gray/diffuse boundaries, and inhomogeneous gray/nonscattering media. The radiative heat is evaluated only at the boundaries, given its small relevance if 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. A discretization consisting of 256 rays for each calculation point and a step of 1 mm along each ray were used, after a convergence analysis for both parameters. A wall emissivity equal to 0.91 was assumed for the paraffin wax grain, with a refractive index of 1.43. 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 reaching convergence. More details on the radiation model are present in [17].

Shape Change
A typical approach used to perform rebuilding of experimental data consists in a set of three simulations performed at the initial, average and final radius. This approach is referred to as simplified shape change. This allows to simulate three different steady-state conditions of the firing, considering a variation of the engine setup in the radial direction.
In particular, starting from the averaged values of oxidizer mass fluxes and oxidizer mass flow rate experimentally measured, it is possible to compute the average radius that corresponds to that condition and, as consequence, the final radius respectively as: To perform a better rebuilding of experimental data or to provide a more precise design tool, it is important to perform a series of steady-state simulations in which the shape of the propellant grain is adjourned each time, taking into account the computed regression rate. Since the characteristic time of grain regression is much larger than the one of the gas dynamic and chemical phenomena in the engine, approximating an engine firing with a series of steady-state simulations is acceptable. The displacements of the computational grid nodes are not uniform all over the grain length but, rather, they vary according to the different calculated values of the regression rate. As the regression rate is defined in the direction normal to the fuel surface, the displacement of a generic point occurs along both the radial and axial directions [23], due to the local surface inclination. This approach is referred to as real shape change.
The fuel grain surface displacement in the axial and radial direction are computed with where Δt = t n+1 − t n is the time interval between times n + 1 and n, n = (n x , n r ) is the normal to the segment joining nodes i − 1 and i + 1 and the indexes i and i + 1 are used to indicate adjacent grain surface nodes, as described in [21].

Axial Injection
In this section, an experimental firing test from [24], named L5, characterized by the highest oxidizer mass flow rate, is rebuilt, simulating the grain shape change evolution in time. In particular, both the simplified shape change and real shape change approaches are carried out, to underline the advantage of the latter in the prediction of space-time averaged regression rate and final grain profile. A description of the flowfield, a grid sensitivity analysis and the rebuilding of the time-averaged experimental regression rate, chamber pressure, and fuel grain shape evolution in time are shown in this section.

Computational Setup
A simple geometrical axisymmetric configuration is considered, where the combustion chamber has a constant cross x n+1 i = x n i +̇r n i Δt n x r n+1 i = r n i +̇r n i Δt n r section ( r = r P ). A portion of the motor is occupied by the paraffin C 32 H 66 fuel grain ( x 0 < x < x 1 ), while the injected oxidizer is gaseous oxygen O 2 . The pre-chamber and postchamber are located respectively before and after the fuel grain ( x < x 0 and The converging-diverging nozzle is drawn with conical sections connected by circular arcs with each other and with the post-chamber. The last portion of the motor, consisting in the divergent part of the nozzle, is not useful to evaluate performances in this work, because thrust and expansion mechanism are not investigated, thus the expansion ratio is truncated to A e ∕A t = 1.5 . An approximation is made on prechamber and post-chamber diameters, since they are considered to be equal to the port diameter, without considering the two cavities for flow mixing. This choice is supported by previous studies [21,25], which demonstrate that, for an axial injection, the presence of a pre-chamber cavity only alters pressure and regression rate values up to a maximum of 2-3%.
In this section, the injection arrangement is made by a central hole and series of circular injectors set around a circumference, which lead to vortex generation in the first part of the motor.
On the left-hand side of the setup, a subsonic inflow boundary condition imposing mass flow rate and static temperature simulates an oxygen injector, with the injector plates between the axial injector and the ring one and above the ring one modeled as adiabatic walls. On the top side, adiabatic walls are imposed outside of the paraffin grain section. Finally, a symmetry condition is applied at the center line because of the axisymmetric flow assumption, and a supersonic outflow is assumed at the outlet section, as shown in Fig. 1. For this test case, r P ,in = 20 mm, [24].
The reference mesh employed is divided into 400 grid points in the axial direction and 200 grid points in the radial direction.
Along the x-axis, there is a refinement at the beginning of the fuel grain and around the throat section, to capture fluid dynamics phenomena even in the presence of strong gradients. In addition, mesh clustering is also used near the injection region, in order to resolve the small recirculating vortex due to the presence of the wall between the two injectors and, above all, the main vortex originating in the prechamber region.
On the r-axis, cells gradually get smaller moving from the symmetry axis to the wall, to capture boundary layer phenomena and guarantee a suitable maximum value of y + always lower than 3 in the fuel grain region.
In the pre-chamber region, a cubic transition is used to follow vortex contour and catching flame anchoring point on fuel grain wall, as shown in Fig. 2.

Grid Sensitivity Analysis
Two different grids have been considered for the grid convergence analysis, which is performed at the initial radius. In fact, the smaller the radius, the higher the values of oxidizer mass flux (and thus Reynolds number) and  nondimensional wall distance y + , which are both critical parameters for the grid resolution. In fact, the maximum value of y + is 2.3 for the initial radius, and 1.65 and 1.03 for the average and final radii, respectively. The fine grid, which is made by 600 × 300 cells, is obtained by multiplying by 1.5 the number of cells in both axial and radial directions with respect to the reference grid. The characteristic spatial discretizations Δx and Δr are therefore always divided by a factor 1.5 by refining the mesh.
A fairly accurate representation of the flowfield can be obtained with the reference grid, which has 80,000 control volumes. In fact, it was observed that the fuel regression rate is weakly sensitive to wall resolution if y + < 3.
The profiles of regression rate and pressure obtained with the fine mesh at the same average port diameter of the reference one are reported in Fig. 3. The profile of ̇r varies with respect to the reference discretization level only slightly after x = 0.16 m, and its integral average is about the same.
The flowfield does not show significant changes with increasing mesh resolution from the medium to the fine level. Therefore, the 400 × 200 mesh is considered acceptable.

Flowfield Description
Considering the simplified shape change approach, the three simulations, at initial, average and final radius ( r P ,in = 0.02 m, r P ,ave = 0.0332 m, r P ,f = 0.0465 m), show several peculiarities of hybrid rockets burning paraffin wax and gaseous oxygen.
In particular, as shown in Fig. 4, a wide diffusion flame is observed throughout the engine. The flame is anchored to the injector edge, due to the recirculation region entailed by the axial injection of oxygen. The vortex develops on the paraffin grain, allowing a large fraction of prevalently cold and unburnt paraffin fuel to circulate in the pre-chamber region. This causes the melted paraffin wax injected from the solid grain, and generally present mainly in a narrow layer close to the grain surface, to move upstream toward the injector. The flame widens and gets closer to the grain and then moves away from the wall further downstream because of the progressive consumption of the oxygen injected. The flame reaches the center line in proximity of the mid grain. The vortex reattachment point on the paraffin grain moves towards the end of the grain during the burn, according to the increase of port diameter due to fuel consumption. Indeed, 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 prechamber and allows the flame to be anchored at the interface between the oxygen injector and the injector plate wall [22].

Experimental Rebuilding
First, post-chamber pressures (at x = 0.44 m) are extracted from the numerical simulations, both for simplified shape change and real shape change, and are compared to the respective experimental values in Fig. 5a. Very good agreement is found for post-chamber pressures p posCC , for which an error of about 11% is observed. It is very interesting to notice that both the two approaches give a good  The integral average of the axial profiles of the regression rates, shown in Fig. 5c, is generally in good agreement with experimental results. Averaged values decrease during the firing, starting from the maximum value at the beginning of the firing, when the mass flux is maximum and port diameter minimum. The maximum point of the curves shifts towards right, according to vortex enlargement in the combustion chamber. After time t = 4.5 s, the curves maintain the same shape, only scaling down the average value because of the port diameter increase. Thus, the regression rate indicates a slight coning effect.
Considering then the space-time-averaged regression rate, the real shape change approach improves the prediction capability of the simplified shape change. In particular, as evident in Fig. 5b, considering the time evolution of the fuel grain, the percentage error in ̄̇r space−time estimation decreases up to 10% with respect to 13% of the former.
The grain profile evolution can be seen in Fig. 5d, in which the fuel grain shape at different times of the firing is reported. The experimental average final radius is matched with good agreement.

Swirl Injection
In this section, simulations performed on a lab-scale hybrid rocket with swirl injection are presented. The test cases, taken from [9], are chosen because of their high rotational intensity and the multiple injectors employed, which allow to investigate the main features of such engines and to provide a first validation of the numerical approach.
A swirl injector (see Fig. 6) feeds the oxidizer inside the engine with a non-zero tangential velocity. Each swirling flow is characterized by its rotational intensity, quantified by the swirl number SN [26], defined as the dimensionless ratio between the axial flux of the angular momentum and the axial flux of axial momentum where u and w are the axial and tangential velocities, r w the wall radius and S the cross section.
Assuming the conservation of angular momentum, the swirl number is usually rewritten as a function of the sole injector geometry, obtaining the geometrical swirl number [27] where r inj = r w − r h is the distance of the injector holes from the axis, r h the injection channel radius and A inj is the injection area. However, the simplification of SN in SN g does come with a loss of generality: SN g does not assume the same value of the integral swirl number, which is typically lower.
All simulations have been performed employing an axisymmetric approach, replacing the discrete injectors with an equivalent annular one, with the same injection area and setting the flow injection angles to obtain the same swirl number of the tests analyzed. An appropriate periodic boundary condition has been implemented and set on the lateral faces of the cells to ensure the symmetry of the flow.

Computational Setup
The computational setup employed for the simulations is shown in Fig. 7. The engine dimensions are r t = 5.25 mm, r preCC = 32.5 mm, r P ,in = 15 mm, L preCC = x 1 = 88.5 mm, For all simulations an oxidizer mass flow rate ṁ ox = 33.2 g/s at 300 K is employed. Three different injectors are simulated (named S06, S18 and S31), with SN g = 45, 135 and 236 respectively. The setup employed is quite similar to the one of Fig. 1, with a swirl injector ( 0 < x < x 0 ) instead of an axial one. In the prechamber, an isothermal boundary condition (with wall temperature equal to oxidizer injection temperature) is imposed, to ensure that, when the flame attaches at the wall, non-realistic values of the radiative heat flux are avoided.
The constant-radius setup employed for the computations requires to impose at the injection a modified swirl number SN inj , instead of SN g . In fact, the swirl intensity of the flow is affected by cross-section variations, and decreases if the cross section reduces, due to flow acceleration [26]. Consequently, when using a constant-radius setup, which neglects to consider the correct prechamber radius, a reduced swirl number has to be used, to model the swirl intensity reduction between the prechamber and the fuel grain. This is of course a simplification, which, however, allows one to use simpler computation grids, with lower computational cost. Assuming that the flow is steady and inviscid, the swirl number varies linearly with the wall radius and thus one gets The computational grid employed for all simulations is made by 310 × 80 cells, with clustering near the injection region (Fig. 8), at the prechamber-grain interface (where the flame attachment point is expected to be) and at the nozzle throat. The grid allows to have y + < 1 for all simulations, guaranteeing a correct resolution of the boundary layer over the fuel grain.

Discussion and Results
Simulations at three port radii ( r P = 16.3 mm, r P = 21.4 mm and r P = 27.8 mm) have been performed for all three (7) SN inj = SN g r P r preCC Fig. 6 Schematics of a swirl injector: r h is the injection channel radius, r w the wall radius and r inj = r w − r h the distance of the injector channels from the axis injectors. The selected port radii are representative of the initial, average and final working conditions of the engine.
The temperature flowfields for r P = 16.3 mm are shown in Fig. 9. Increasing the swirl intensity the combustion zone shifts towards the grain leading edge and the mean flow temperature is reduced, because the increased fuel mass flow rate leads to a more fuel-rich O/F ratio. Figure 10 shows the projection of the streamlines on the symmetry plane. A first recirculation zone (top left of the image) is formed just next to the injector. This vortex is caused by the injection setup, which gives the oxidizer a non-zero radial velocity, being not possible, in an axisymmetric approach, to have only tangential velocities at the injection. Due to the centrifugal forces, however, the flow is quickly pushed back towards the wall, and a second, larger vortex is formed on the axis. Being the swirl intensity larger than the critical value of SN = 0.6 , the formation of a central toroidal recirculation zone (CTRZ) is indeed expected [26]. Similar structures have also been found by [12].
The strong flow recirculation leads to a significant decrease in the swirl number, which can be computed from the simulations using its integral definition. As shown in Table 1, the actual swirl intensity at the grain leading edge SN grain is significantly lower than the injection swirl number. The decay is stronger for small port radii and high swirl intensities, since both lead to higher flow velocities and  increased viscous dissipation. This shows how, depending on the setup, the geometric swirl number may not be representative of the real conditions inside the engine. Fig. 11 shows the swirl intensities at the prechamber-grain transition for all simulations. The simulations with same SN g yield similar swirl intensities, although the injection swirl numbers SN inj are quite different, due to the strong decay phenomena. Furthermore, the swirl intensity decays rapidly along the grain, due to the effect of combustion and mass adduction, which increase the axial velocity of the flow and cause a reduction of the rotational intensity.
The local fuel regression rates are shown in Fig. 12a for all injectors at r P = 16.3 mm. The regression rate increases with swirl intensity, due to the higher convective heat flux (Fig. 12b). The radiative heat flux, instead, becomes on average lower with higher swirl numbers, due to the more fuel rich O/F ratio and the lower average flow temperature, but maintains the same qualitative profile, with the peak shifting towards the grain leading edge. This is due to the increased regression rate, which leads to the maximum temperature being reached at a smaller abscissa, and to the stronger prechamber vortex, which pulls the hot gases towards the prechamber and generates a thicker flame in the first part of the grain. In Fig. 13b, the radiative and convective heat fluxes for the S06 injector are shown. The former increases with the port radius, due to the increased volume of emitting gas, and the latter decreases with r P , due to the lower mass flux. These opposing effects explain why, at least in the first half of the grain, the regression rate remains almost constant during the burn (Fig. 13a).
In Fig. 14 the time-averaged regression rates obtained from the simulations are compared with the experimental ones, obtained from the aG n ox laws reported in [9]. In spite of the simplifying hypotheses made on the engine geometry and of the complexity of the test case, the regression rate is predicted with good accuracy for the S06 and S18 injectors. For what concerns the S31 injector, the regression rate is instead overpredicted by ≈ 30% . It must be also noted that the S31 experimental regression rate is almost identical to the one of the S18 injector, likely due to the high swirl decay in the prechamber, which leads to almost identical conditions on the fuel grain in both cases. This is indeed suggested also by the simulations: Table 1 shows how at r P = 16.3 mm SN grain is the same between the two configurations, however this is not the case for r P = 21.4 mm and r P = 27.8 mm. The inability of the numerical model to correctly predict this phenomenon may be caused by the simplifications made in the engine geometry, which do not allow to correctly model the interaction between the swirling flow and the prechamber-grain transition, or by the adopted turbulence model. In fact, RANS simulations are known to lack the precision to model the details of complex, high-swirl phenomena, mainly because of the presence of anisotropic characteristics in the turbulent stress tensor [26]. However, considering the large computational cost of more complex RANS turbulence models (such as the Reynolds Stress Model) or of Large Eddy Simulations, the adopted approach is deemed a good compromise between computational time and predictive capabilities.

Conclusion
In this work, a numerical approach based on Reynoldsaveraged Navier-Stokes simulations, with sub-models for fluid-surface interaction, radiation, chemistry, and turbulence has been validated. The effects of fuel grain consumption, on both radial and axial directions, and swirl injection have been taken into account. The numerical rebuilding of regression rate and pressure yields a reasonable agreement with the available experimental data in axial-injected HREs, highlighting noticeable improvement in space-time-averaged regression rate estimation by considering the grain shape evolution in time with respect to the three diameters approach. Despite the geometrical simplifications, the absence of a turbulence-chemistry interaction model, and the axisymmetric hypothesis neglecting the 3D features of the flow, which can be important near the injection zone, a good estimation of the regression rate and of chamber pressure is provided, with an error of only 10%. The same numerical approach has been applied also to HREs with swirl injection. Significant features of such engines have been highlighted and the increase of the fuel regression rate with swirl injection has been analyzed. The rebuilding of experimental regression rate data has been performed, obtaining good results at moderate swirl intensities, where an error of about 10% is considered acceptable, given the simplifications made in the engine geometry and due to the simplified turbulence model adopted. On the contrary, at higher swirl intensities, the results are not satisfactory, likely because the simplifying hypotheses do not allow to correctly predict the strong swirl decay in the prechamber and thus the actual rotational intensity of the flow.
Funding Open access funding provided by Università degli Studi di Roma La Sapienza within the CRUI-CARE Agreement.

Conflicts of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.