Modeling of gas–surface interface for paraffin-based hybrid rocket fuels in computational fluid dynamics simulations

Numerical simulations of the flowfield in a hybrid rocket engine are carried out with a multispecies chemically reacting Reynolds-averaged Navier–Stokes (RANS) solver which includes detailed gas–surface interaction (GSI) modeling based on surface mass and energy balances. The oxidizer is gaseous oxygen which is homogeneously fed into single-port cylindrical grains. The modeling of GSI already developed and validated for pyrolyzing fuels such as hydroxyl-terminated polybutadiene (HTPB), is extended to the case of liquefying fuels, such as paraffin wax. A simplified two-step global reaction mechanism is considered for the gas-phase chemistry to model the combustion process inside the chamber. Numerical simulations performed at different gas/melt-layer interface temperatures and oxygen mass fluxes show a considerable increase of fuel regression rate, in the range of 3 up to 5 times, for the liquefying fuel with respect to the pyrolyzing one. Results show that the regression rate enhancement is significant only when the gas/melt-layer interface of the liquefying fuel is close to the melting temperature. At increasing gas/melt-layer interface temperatures, the regression rate decreases following an inverse power law and gets close to that of a pyrolyzing fuel for the same operating conditions. Finally, regression rate behavior at varying oxygen mass flux of liquefying fuels is not substantially altered from that of pyrolyzing fuels as the oxidizer flux exponent remains rather constant.


INTRODUCTION
The intrinsic properties of hybrid propellant rockets in terms of performance, simplicity, safety, reliability, low development cost, reduced environmental pollution, and §exibility make them one of the envisaged future generation propul-3 Article available at https://www.eucass-proceedings.eu or https://doi.org/10.1051/eucass/201911003 sion systems. On the other hand, there are some technical challenges to be overcome before achieving the same level of maturity as solid and liquid traditional systems, such as low regression rate, reduced combustion e©ciency, and combustion instability. Therefore, the hybrid rocket engine development requires deeper understanding of the physico-chemical phenomena that control the combustion process and of the §uid dynamics inside the motor. The knowledge of the complex interactions among §uid dynamics, solid fuel regression, oxidizer atomization and vaporization, mixing and combustion in the gas phase, nozzle thermochemical erosion, particulate formation, and radiative characteristics of the gas and the §ame can only be improved by combined experimental and numerical research activities.
Especially, low fuel regression rate is one of the most signi¦cant shortcoming of conventional hybrid rockets. So far, many techniques (such as swirling oxidizer injection, insertion of diaphragms in the combustion port, use of complex multiport grain con¦gurations, etc.) have been tried to improve the regression rates of hybrid rockets but all of these methods su¨er from important drawbacks. One promising technique to solve the slow regression problem relies on the use of para©n-based fuels. The regression rate of para©n-based fuels is, in fact, about 35 times higher than that of conventional pyrolyzing fuels, such as HTPB. The low surface tension and low melt layer viscosity of para©n-based fuels generate a signi¦cant entrainment of melted fuel droplets into the oxidizer §ow and this additional mass transfer mechanism can potentially increase the regression rate by several times with respect to conventional pyrolyzing fuels. This is, mainly, due to the combined e¨ect of reduced e¨ective heat of gasi¦cation, decrease of blocking e¨ect in the boundary layer, and possibly higher heat transfer due to increased surface roughness. Di¨erently from conventional pyrolyzing fuels, in fact, the regression of para©n-based fuels, thanks to their low melting point, proceeds almost without vaporization and liquid fuel is mainly supplied to the §ame zone by the entrainment phenomena that plays a dominant role in determining the regression rate characteristics. For oxidizer mass §ux levels encountered in practical hybrid rocket applications, droplet entrainment can dominate direct fuel gasi¦cation. Hence, the regression rate model for these liquefying propellants must be based on a completely di¨erent mass transfer mechanism involving the entrainment of liquid droplets from the melt layer.
The theoretical [13], numerical [46], and experimental modeling [712] of the §ow in the combustion chamber and nozzle of a para©n-based hybrid propellant rocket has been the subject of considerable interest recently but adequate models are still lacking. The modeling e¨ort is, in fact, challenging and requires the ability to adequately describe the interaction between the reacting §ow and the solid surface, which includes the melt-layer, through suitable GSI modeling. The objective of this study is the numerical simulation of the §ow inside the combustion chamber and nozzle of a para©n-based hybrid rocket by solv-ing the RANS equations for multicomponent turbulent reacting §ows, including the required submodels in order to describe the homogeneous combustion in the gas phase and the §uidsurface interaction among the exhaust gases, the meltlayer, and the para©n-based solid propellant. The basis of the present analysis is the comprehensive theoretical model for HTPB pyrolysis and combustion in hybrid rocket engines developed and validated in previous studies [13,14], which is extended here to the case of para©n-based liquefying fuels.

THEORETICAL AND NUMERICAL MODEL
The study of hybrid rocket engine §ow¦elds requires a suitable modeling of both GSI and gas phase reactions. In the present study, the modeling of GSI already developed and validated for pyrolyzing fuels such as HTPB [1315] is extended to the case of liquefying fuels in hybrid rocket engines. In particular, a GSI model based on surface mass and energy balances is coupled with a three-dimensional (3D) chemically-reacting computational §uid dynamics (CFD) code. The CFD tool is a ¦nite-volume solver for 3D turbulent compressible multicomponent reacting §ows [16,17] with variable thermodynamic and transport properties. The thermodynamic properties of individual species are approximated by seventhorder polynomials of temperature and the transport properties are approximated by fourth-order polynomials [18]. Mixture properties for conductivity and viscosity are derived from the Wilke£s rule. The di¨usion model is based on an e¨ective di¨usion coe©cient obtained assuming a constant Schmidt number.
The SpalartAllmaras one-equation model [19] is used to compute the turbulent viscosity. Turbulent thermal conductivity and turbulent mass di¨usivity are computed from the turbulent viscosity, speci¦c heat at constant pressure, turbulent Prandtl number, and turbulent Schmidt number. The SpalartAllmaras turbulence model has been selected because of its simplicity and capability to model internal §ows in rocket motors, at least at a qualitative level. The chemistry is modeled assuming a zero-dimensional perfectly stirred reactor or laminar §ame submodel, i. e., turbulencechemistry interaction is neglected. The code solves the time-dependent conservation equations of mass, momentum, and energy for the chemical nonequilibrium §ow¦eld and adopts a standard ¦nitevolume Godunov-type formulation. It is second-order accurate in space and uses multiblock structured meshes. The system of equations is approximated by a cellcentered ¦nite-volume scheme. The viscous §uxes are approximated by centered di¨erencing, whereas the convective §uxes are computed by means of the solution of a Riemann problem whose left and right states are reconstructed by an interpolation procedure which uses the minmod limiter. The system of ordinary di¨erential equations is advanced in time by means of an explicit RungeKutta integration.

GasSurface Interface
Because of the chemically active surface, further physical modeling is necessary for the GSI. To complete the formulation of the theoretical model, boundary conditions must be speci¦ed at the gassolid interface, which describe the physics of the surface phenomena. The boundary conditions are detailed in the following for both liquefying fuels (such as para©n wax) and nonliquefying pyrolyzing fuels (such as HTPB).

Pyrolyzing fuels (hydroxyl-terminated polybutadiene)
If it is assumed that no material is being removed in a condensed phase (solid or liquid), then the general conservation laws at the gassolid interface over a pyrolyzing fuel grain can be written as [20]: where ' m w is the overall mass §ux injected from the pyrolysis process per unit surface area of the wall; ρ and v are the wall density and the normal-to-wall velocity in the gas phase due to pyrolysis products injection because of grain regression; and, ¦nally, ρ s and ' r are the solid material density and regression rate, respectively; where D im is the species-to-mixture di¨usion coe©cient; y i is the gas phase mass fraction of the ith species at the wall; N is the number of species; ' ω i is the rate of production of gas-phase species i at the surface due to fuel pyrolysis; and η is the inward (from solid to gas) coordinate normal to surface; and where h w is the enthalpy of the gas mixture at wall; h i is the enthalpy of the single gas species; h s is the enthalpy of the solid grain at the wall temperature; k is the gas conductivity; T is the gas temperature; and ' q ss cond is the conduction heat §ux inside the grain. Note that the term ' m w h s is the energy §ux entering the surface due to grain regression. The radiation contributions do not appear in the energy balance as they have not been modeled in the present study.
Concerning the energy transfer into the fuel grain, it is assumed that heat conduction is dominant in the direction normal to the grain surface. Although axial temperature gradients exist along the grain wall, they are generally small if compared with the normal heat conduction and represent a second-order e¨ect. In a moving local coordinate system tied to the receding fuel surface, the general one-dimensional (1D) form of the in-depth energy balance (assuming constant material properties) is: where the di¨erent terms (from left to right) represent the temporal variation of the grain sensible energy, the net conduction inside the material, and the convected energy due to the coordinate motion. The term α s = k s /(ρ s c s ) is the fuel grain thermal di¨usivity and T s is the grain local temperature. Assuming steady-state condition and integrating the stationary form of Eq. (3) between gasfuel interface w and a point su©ciently far from the wall so that an adiabatic condition can be assumed for the fuel at initial state i, the steady-state conduction term ' q ss cond can be expressed as: where c s indicates the fuel heat capacity; T w is the grain wall temperature; and T i is the grain initial temperature. The steady-state assumption is a reasonable approximation when the thermal lag in the solid is su©ciently small, which actually occurs in the operating conditions to be simulated, i. e., moderately high fuel regression rate and low fuel thermal di¨usivity. The energy balance expression can be recasted using Eq. (2) emphasizing the contributions due to convection from the gas, conduction into the solid propellant grain, and chemical reactions, as: The chemical heat §ux can be expressed as: where the term in parentheses represents the heat absorbed by the surface reactions and, hence, it represents the so-called heat of pyrolysis, -h pyro . The rate of production of gas-phase species i at the surface, ' ω i appearing in Eqs. (2) and (5), has to be estimated on the basis of the fuel pyrolysis model. In this work, according to [21], a single Arrhenius-type equation is used to model ' m w Table 1 Heterogeneous rate constants for HTPB [22] Surface reaction A, kg/(m 2 s) E, kJ/mol HTPBs → C4H6 2208 56.5 Table 2 The HTPB properties [23,24] Density ρs, Speci¦c heat cs, Thermal conductivity ks, Heat of pyrolysis, kg/m 3 J/(kg·K) W/(m·K) MJ/kg 960 1632 [24] 0.217 [24] 1.1 [23] assuming that the main product of HTPB pyrolysis is 1,3-butadiene (C 4 H 6 ). In particular, the rate of pyrolysis is obtained as where T w is the wall temperature and A and E a are the preexponential factor and the activation energy for the pyrolysis reaction, respectively. Coe©cients are taken as reported in Table 1, according to Arisawa and Brill [22]. As a single pyrolysis product is assumed, the species production rate at the fuel surface, ' ω i , is equal to ' m w for species C 4 H 6 and to zero for all the remaining gas-phase species. However, it is important to note that, although the only pyrolysis gas considered here is 1,3-butadiene, the gas mixture at the fuel surface is not entirely composed of C 4 H 6 as the other gaseous species (both oxygen and combustion gases) can actually reach/leave the surface due to di¨usion and convection induced by blowing. The mixture composition at the HTPB wall, in fact, is solved enforcing the species surface mass balance, Eq. (2) coupled with the surface energy balance, Eq. (5). This is very important as solving both the mass and species balances Eqs. (1) and (2) guarantees that the correct amount of fuel is injected in the §ow¦eld from the grain surface. In Table 2, the HTPB properties used for the simulations are listed. The heat conducted in the solid grain is modeled according to Eq. (4) where a constant speci¦c heat has been assumed for the solid HTPB, according to [24]. The solid fuel thermal conductivity k s is also assumed constant with temperature. Finally, a value of 1.1 MJ/kg has been assumed for the heat of pyrolysis of HTPB as performed in [25]. The heat of pyrolysis could be evaluated from Eq. (6) and the assumption that the pyrolysis gas is only composed of 1,3-butadiene, however, since this assumption is an approximation, in this work the heat of pyrolysis has been derived from the available experimental data [23].

Liquefying fuels (para©n wax)
If it is assumed that the material is dominantly removed by the entrainment phenomena (negligible evaporation, boiling, or pyrolysis fuel removal process), then the general mass conservation at the gas/melt-layer interface over a liquefying fuel grain can be written as where ' m ent is the overall mass §ux injected due to the entrainment phenomena per unit surface area of the wall; ρ l is the liquid material density; v l is the liquid velocity perpendicular to the surface in the melt-layer; and ' r l is the gas/meltlayer interface regression rate. From the melt-layer/solid interface mass balance, the liquid velocity in the melt-layer v l is expressed as (ρ s /ρ l − 1) ' r s where ' r s is the melt-layer/solid interface regression rate. If the steady-state condition is assumed (constant thickness of the melt-layer), then the gas/melt-layer interface regression rate and the melt-layer/solid interface regression rate are equal so that ' r l = ' r s = ' r and the entrainment mass §ux becomes Para©n wax (C 32 H 66 ) is assumed to decompose/gasify into C 2 H 4 and molecular hydrogen through the cracking reaction C 32 H 66 ⇒ H 2 + 16C 2 H 4 . As a singlephase approach is used in this study, the cracking of the entrained liquid droplets and the associated heat of cracking are not modeled, and gaseous C 2 H 4 and H 2 are directly injected from the fuel surface into the combustion chamber. It is worth noting that, despite this relevant simpli¦cation, the following GSI approach retains its validity as the decomposition of the entrained droplets and the corresponding heat absorption take place away from the fuel interface and inside the combustion chamber where plenty of thermal energy is available. Therefore, the overall surface mass balance becomes: where ρ and v are the wall density and the normal-to-wall velocity in the gas phase due to gaseous fuel products injection because of grain regression. In the assumption that the entrained droplets are directly injected from the fuel wall in the gaseous phase, the species mass balance becomes analogous to that of a pyrolyzing fuel (see paragraph 2.1.1): where ' ω i is nonzero only for C 2 H 4 and H 2 . Finally, the energy balance at the surface is: Equation (9) di¨ers substantially from Eq. (5) since it does not include the chemical contribution. The entrainment process, in fact, is a mechanical phenomena related to the shear stress at the gas/melt-layer interface and, therefore, there is no chemical energy associated with this kind of fuel mass removal. Again, the radiation contributions do not appear in the energy balance as they have not been modeled in the present study. Concerning the energy transfer into the fuel grain (melt-layer and solid), as done for the pyrolyzing fuel, it is assumed that heat conduction is dominant in the direction normal to the grain surface. The heat absorbed by the melted fuel §owing from the melt-layer/solid to the gas/melt-layer interface is accounted for through the liquid velocity in the melt-layer, v l . As far as the energy transfer in the melt-layer is concerned, the melted fuel §ow in the direction parallel to the fuel grain is neglected here. In a moving local coordinate system tied to the receding gas/melt-layer interface ( ' r l ) and melt-layer/solid interface ( ' r s ), the general 1D form of the in-depth energy balance (assuming constant material properties) is: (Energy in the solid fuel layer).
If the steady-state condition is assumed ( ' r = ' r l = ' r s ), then the stationary form of the in-depth energy balance is: α s ∂ 2 T s ∂η 2 = ' r ∂T s ∂η (Stationary energy in the solid fuel layer).
The terms α l = k l /(ρ l c l ) and α s = k s /(ρ s c s ) are the liquid and solid fuel grain thermal di¨usivities, respectively, and T l and T s are the grain liquid and solid local temperatures, respectively. Integrating Eqs. (10) and (11) within the liquid layer and solid layer thicknesses, respectively, and considering the surface energy balance at the melt-layer/solid interface, the steady-state conduction term ' q ss cond can be expressed as: where c l and c s indicate the  [26,27]; and (c) normal boiling (1) and melting (2) temperatures [26,27] melt-layer and solid fuel heat capacities, respectively. The terms -h melt and T m represent the para©n wax melting enthalpy and melting temperature, respectively. The term T w is the melt-layer temperature at the melt-layer/gas interface and T i is the solid fuel initial temperature.
The rate of melt-layer entrainment at the surface, ' m ent appearing in Eqs. (7) and (12), has to be estimated on the basis of a fuel entrainment model, which is typically based on semiempirical correlations. In this work, however, the attention is rather focused on the relationship between the surface temperature and the entrainment mass §ux so that a fuel entrainment model is not used. The entrainment mass §ux, hence, is determined from the surface energy balance, Eq. (9), coupled with the steady-state heat conduction inside the fuel melt-layer and solid layer, Eq. (12), assuming an isothermal wall condition. From the knowledge of the gas/melt-layer interface temperature, which is, in turn, Figure 2 Liquid para©n wax properties as functions of carbon number: (a) density [26,27]; (b) thermal conductivity [26,28]; (c) heat capacity [26,28]; 1 ¡ T = 300 K; 2 ¡ 400; 3 ¡ 500; 4 ¡ 600; and 5 ¡ T = 700 K related to the melt-layer thickness, the entrainment mass §ux can be obtained solving Eqs. (9) and (12) iteratively coupled with Eqs. (7) and (8), which are needed to solve for the gas surface chemical composition. In Table 3, the para©n wax properties (at solid state) used for the simulations are listed. The heat conducted in the solid grain is modeled according to Eq. (12) where a constant spe-ci¦c heat has been assumed for the solid para©n. The solid wax thermal conductivity k s is also assumed constant with temperature. A value of 0.17 MJ/kg has been assumed for the melting enthalpy of C 32 H 66 , shown in Fig. 1a, as performed in [3]. The critical pressure and temperature and the melting and normal boiling temperatures as a function of carbon number are shown in Figs. 1b and 1c and have been estimated using the asymptotic behavior correlation (ABC) method as discussed in [2628]. Finally, the liquid wax density, thermal conductivity, and heat capacity as functions of temperature and carbon number are shown in Fig. 2 and have been estimated using the ABC method [2628]. The melt-layer properties are assumed constant over the melt-layer thickness for simplicity and have been evaluated at the average temperature between the gas/melt-layer interface temperature (surface temperature) and the melt-layer/solid interface temperature (melting temperature). As shown in Fig. 1b, the critical pressure of C 32 H 66 is 6.5 bar so that, as the computed chamber pressure for the analyzed conditions is always above this value, liquid C 32 H 66 cannot evaporate or boil, being at a supercritical pressure. Also, as in the performed numerical simulations, the gas/melt-layer interface temperature is kept below the C 32 H 66 critical temperature (860 K), a supercritical liquid layer cannot be formed. Consequently, the only additional mass transfer mechanism other than that of entrainment is the pyrolysis of the melt-layer. In the current analysis, the pyrolysis of the melt-layer is neglected and the the gas/melt-layer interface temperature is con¦ned between the melting temperature (340 K), which corresponds to a null melt-layer thickness, and the normal boiling temperature (740 K) of C 32 H 66 , which corresponds to the maximum melt-layer thickness. Melt-layer pyrolysis analysis for C 32 H 66 has been performed in [3] showing that, for surface temperatures lower than ≈700 K, the pyrolysis fuel regression rate becomes negligibly small. Therefore, for the current analyzed conditions, the assumption that the dominant mass removal mechanism is due to the entrainment phenomena is justi-¦ed.

Gas Phase Reactions
Finite-rate gas phase reactions are modeled by global reaction mechanisms, because detailed chemical kinetics mechanisms would include many species and would be, on the one hand, computationally heavy and, on the other hand, beyond the scope of the present study, whose aim is to focus on GSI. Therefore, ¦nite-rate gas phase reactions are modeled by a simpli¦ed two-step global reaction mechanism which involves two global reaction steps. In case of pyrolyzing fuels, such as HTPB, the ¦rst irreversible global reaction step is between C 4 H 6 , which is the HTPB pyrolysis product, and molecular oxygen to form CO and H 2 O and is considered ¦rst order in both fuel and oxidizer. The second reversible global step accounts for the formation of CO 2 : The resulting rates of production and destruction of species are: ; Table 4 Reaction rate constants for reactions (13) and (14) Reaction rate A n Ea/R, K k f 1 [29] 4.9486 · 10 9 0.0 15200 k f 2 [30] 2.2400 · 10 6 0.0 5032.7 k b 2 [30] 1.1000 · 10 13 −0.97 39456.5 The forward and backward reaction rates k f and k b for the two reactions are expressed as Arrhenius functions in the form k = AT n exp(−E a /(RT )), and the values of the constants used in this study are tabulated in Table 4. The ¦rst global step reaction rate is taken from [29] and the second step is taken from [30], which is a modi¦ed version of the one used in [21,29] to account for oxy-fuel combustion conditions. When para©n-based liquefying fuels are considered, para©n wax (C 32 H 66 ) is assumed to decompose/gasify into C 2 H 4 and molecular hydrogen through the reaction C 32 H 66 ⇒ H 2 + 16C 2 H 4 . As a single-phase approach is used in this study, the pyrolysis/vaporization of the entrained liquid droplets is not modeled and gaseous C 2 H 4 and H 2 (assumed as inert) are directly injected from the fuel surface into the combustion chamber.
The simpli¦ed two-step global reaction mechanism for C 2 H 4 is: The resulting rates of production and destruction of species are: The forward and backward reaction rates k f and k b for the two reactions are expressed as Arrhenius functions in the form k = AT n exp(−E a /(RT )) and the values of the constants used in this study are tabulated in Table 4. In particular, the forward rate for the ¦rst irreversible global reaction step is assumed identical for both C 2 H 4 and C 4 H 6 .

TEST CASE DESCRIPTION
The results have been obtained on a typical cylindrical port geometry with a prechamber and an aft-mixing chamber whose dimensions and §ow conditions are shown in Fig. 3. The computational domain refers to a simpli¦ed geometrical Figure 3 Computational domain: numerical grid (not to scale) and boundary conditions representation of a real hybrid rocket chamber geometry where the prechamber, fuel grain, and postchamber have been schematized with constant cross section. The length of the domain portion corresponding to the fuel grain is 574 mm, the port diameter is equal to 42 mm, and the throat diameter is equal to 16 mm. Gaseous oxygen is injected throughout the entire fuel port at a total temperature of 300 K and at di¨erent mass §uxes. The fuel gaseous mixture which is injected from the wall in a gaseous oxygen stream is assumed to be C 4 H 6 in case of HTPB and a mixture of C 2 H 4 and H 2 in case of para©n wax (C 32 H 66 ).

RESULTS AND DISCUSSION
The results presented here are restricted to single-phase system; therefore, in case of liquefying fuels, the entrained fuel liquid droplets are not modeled and the fuel is directly injected in a gaseous state from the gas/melt-layer interface. This is an approximation of the real phenomena as the liquid droplets do not vaporize directly at the surface but are rather transported inside the §ow¦eld and, hence, they should not contribute to the blockage e¨ect, which is typical of pyrolyzing fuels. With a single-phase approach, however, the gaseous injected fuel is contributing to the blockage e¨ect. This e¨ect, which reduces the heat transfer to the fuel grain, is possibly counteracted by the absence of the endothermic cracking reaction and the associated heat of cracking. All the computations presented here are two-dimensional (axisymmetric) and at the steady-state condition. For each value of oxygen mass §ux enforced, the CFD simulation provides the fuel mass §ux (from the surface mass and energy balance) and the chamber pressure level as a consequence of the choked condition at the nozzle throat. Therefore, the chamber pressure level attained in the motor depends on the mixing and combustion process and also on the interaction with the pyrolyzing (HTPB) or liquefying (para©n wax) fuel. In this work, the e¨ect of the fuel port variation is not analyzed, so that all simulations have been performed at a ¦xed port diameter as shown in Fig. 3. The following results are obtained for para©n wax C 32 H 66 at di¨erent gas/melt-layer interface temperatures and for various oxygen mass §uxes and are compared with those obtained using HTPB at the same operating conditions.

Role of Surface Temperature
The role of the gas/melt-layer interface temperature on entrainment mass §ux is analyzed on a reference con¦guration with an oxygen mass §ux of 94 kg/(m 2 s) for C 32 H 66 para©n wax. A reference simulation using HTPB is also performed. The gas/melt-layer interface temperature is assigned in the computation and assumed constant throughout the entire 574-millimeter fuel grain length. It is varied from the melting temperature (340 K) to the normal boiling temperature (740 K) with constant intervals of 100 K. In case of HTPB, the surface temperature is not imposed at the fuel interface but it is calculated solving the surface energy balance, Eq. (5), and, hence, it is varying along the fuel grain length typically ranging from 750 to 850 K. Figures 4 and 5 shows the motor internal ballistics as functions of gas/meltlayer interface temperature at constant port area and oxygen mass §ux. As seen, the gas/melt-layer interface temperature is signi¦cantly a¨ecting the §ow¦eld with a much higher production of the gaseous fuel (which is predominantly C 2 H 4 ) at the lower interface temperatures. Figure 6a shows the typical regression behavior as predicted by classical boundary-layer theory with an initial regression rate decrease due to the boundary-layer growth and, then, an increase due to the mass and heat addition. The highest fuel regression rate shown in Fig. 6a is obtained at the lowest interface temperature (i. e., the wax melting temperature) due to the increase of the convective heating from the hot gases and the decrease of the conduction inside the fuel grain (see Eqs. (9) and (12)). Due to the higher regression rate of para©n wax with respect to HTPB, the §ame is shifted away (see Fig. 4) from the surface by the gaseous fuel injected into the fuel port, the chamber pressure is increased (see Figs. 5 and 6b), and more oxygen is consumed (see Fig. 5). The peak §ame temperature and the axial velocity inside the chamber, on the other hand, are not signi¦cantly a¨ected by the increasing wax regression rate. As shown in Fig. 6a, as the interface temperature approaches the maximum temperature (normal boiling temperature), the regression rate decreases and becomes similar to that of HTPB at the same port area and oxygen mass §ux. At the melting temperature, the wax regression rate is more than 4 times higher than HTPB while at the normal boiling temperature, the wax regression rate is only 1.3 times that of HTPB. Figure 7 shows the para©n wax spatially averaged regression rate as a function of gas/melt-layer interface temperature for di¨erent oxygen mass §uxes (ranging from 47 to 282 kg/(m 2 s)). As shown in Fig. 7b, the para©n wax regression rate variation with the gas/melt-layer interface temperature is very well matched by an inverse power law of the form ' r = a/T n w with the index n varying from 1.45 to 1.6. The para©n wax regression rate increase with decreasing gas/melt-layer interface temperature (see Fig. 7a) is following an inverse power law as the decrease of the interface temperature produces an increase of the convective heat §ux from the combustion gases and a concurrent decrease of the conductive heat §ux into the propellant grain and both contribute to increase the fuel regression rate. Figure 7b clearly shows that a signi¦cantly higher regression rate for para©n wax with respect to HTPB (from 3 to more than 4 times) can only be obtained when the interface temperature is not too far from the wax melting temperature. The results also show an increase of the para©n wax regression rate enhancement with respect to HTPB at increasing oxygen mass §ux.    At the lowest oxygen mass §ux of 47 kg/m 2 s, in fact, the para©n wax regression rate at the melting temperature is 4 times than that of HTPB, while it increase up to 5 times at the maximum considered oxygen mass §ux of 282 kg/(m 2 s).

Role of Oxygen Mass Flux
The role of the oxygen mass §ux on the entrainment mass §ux at di¨erent ¦xed gas/melt-layer interface temperatures is analyzed in Fig. 8. The results of Fig. 8b show the typical regression rate power law relationship of the form ' r = aG n ox . The  gas/melt-layer interface temperature is directly a¨ecting the a parameter while it is only slightly modifying the oxidizer §ux exponent n, which remains close to a classical value in the range of 0.8. At the boiling temperature, in fact, the oxidizer §ux exponent is equal to 0.79 and it increases up to 0.86 at the melting temperature, hence showing a limited di¨erence from the obtained HTPB value of 0.75. Most of the e¨ect, in fact, is on the a parameter which is increased by 2.4 times passing from the boiling temperature to the melting temperature. Therefore, the variation of the gas/melt-layer interface temperature, although signi¦cantly a¨ecting the fuel regression rate, does not alter in a substantial way the regression behavior at varying oxygen mass §ux as the oxidizer §ux exponent remains rather constant and close to the value of HTPB. In conclusion, the para©n wax regression rate enhancement with respect to HTPB is shown to be mostly a¨ected by the gas/melt-layer interface temperature (and, hence, by the melt-layer thickness) and only partially in §uenced by the oxidizer mass §ux, although a limited increase is also shown with increasing G ox due to a small increase of the oxidizer §ux exponent.

CONCLUDING REMARKS
A CFD approach to the internal ballistics of a hybrid motor with integrated GSI modeling capabilities for pyrolyzing fuels has been extended to the case of liquefying fuels, such as para©n wax. The §uid dynamic equations are the RANS equations with additional transport equations for chemical species and turbulence, which are coupled to simpli¦ed two-step global reaction mechanisms to model the combustion process inside the combustion chamber. The gas-phase equations are solved and coupled to the liquid/solid fuel phase using surface balances of mass and energy, which enable the determination of the fuel regression rate as a part of the §ow¦eld solution. A test case represented by the combustion in a cylindrical single-port hybrid rocket burning gaseous oxygen has been simulated. The considered fuel is either HTPB or para©n wax. In case of paraf-¦n wax, simulations have been performed varying the gas/melt-layer interface temperature to analyze its role on the entrainment mass §ux at di¨erent oxygen mass §uxes. Results have shown that the entrainment regression rate increases with decreasing interface temperature following an inverse power law. The obtained inverse power law is due to the increase of the convective heating from the hot gases and the concurrent decrease of the conductive heating inside the fuel grain at decreasing interface temperatures. The entrainment regression rate can reach levels up to 4-5 times higher than those of HTPB at the same operating conditions when the interface temperature is close to the melting temperature, i. e., the melt-layer thickness gets close to zero. Di¨erently, at the highest considered interface temperature which is assumed equal to the normal boiling point temperature, the entrainment regression rate is not signi¦cantly higher with respect to HTPB. Finally, the role of the oxygen mass §ux on the entrainment mass §ux at di¨erent ¦xed gas/melt-layer interface temperatures has been analyzed showing that the oxidizer §ux exponent of the typical regression rate power law relationship is not substantially altered when liquefying fuels are considered. The oxidizer §ux exponent, in fact, is only slightly increasing with decreasing interface temperatures but still it remains close to a classical value in the range of 0.8, as for HTPB. In conclusion, the para©n wax regression rate enhancement with respect to HTPB is shown to be mostly a¨ected by the gas/melt-layer interface temperature and, hence, by the melt-layer thickness, and only partially in §uenced by the oxidizer mass §ux due to the limited increase of the oxidizer §ux exponent. This work represents the ¦rst step towards the modeling of the complex interaction between the combustion gases and a liquefying fuel in hybrid rocket engines but still many areas require a modeling e¨ort to enhance the entrainment phenomena modeling. Future work will include an entrainment model in order to remove the assumption of a ¦xed interface gas/melt-layer interface temperature, a simpli¦ed two-phase treatment to describe the evolution of the entrainment liquid droplets inside the chamber, better combustion models, and the inclusion of radiation exchange between the gas and the fuel.