**1. Introduction**

The hybrid rocket is a chemical propulsion engine in which the fuel and oxidizer are separated in different physical states [1]. In the classical system configuration, the fuel is stored in the combustion chamber in the solid state and a liquid or gaseous oxidizer is injected into single or multiple ports obtained in the solid fuel grain. When the two propellants are ignited, a diffusive flame is formed in the boundary layer developing in the grain port, relatively far from the fuel surface, and it is fed from the outer side by the oxidizer, which is transported from the free stream by turbulent diffusion mechanisms, and, from the inner side, by the products of fuel gasification and/or liquefaction (depending on the nature of the fuel itself) which is sustained by the flame. The combusted mixture expands through an exhaust nozzle, generating the required thrust.

In the last few years, a significantly growing interest has been addressed towards hybrid rocket propulsion, thanks to its numerous advantages [2] relative to traditional solid and liquid systems,

which include, among the main features, the re-ignition and throttling capabilities combined with the possibility of employing environmentally sustainable propellants and, most importantly, its intrinsic safety [3,4]. Although the potentialities of the hybrid rocket have been widely recognized to warrant the renewed research e fforts that are being devoted to its development, such propulsion technology still raises challenging tasks.

Engine performance is governed by the rate at which the fuel is consumed, i.e., by the fuel regression rate, as this determines the total mass flow rate and overall oxidizer-to-fuel mixture ratio, which, for a given chamber pressure, controls the motor thrust and the ideal specific impulse. For given operating conditions, the regression rate depends on the type of solid fuel employed. Fuels typically burned in hybrid rockets are classical polymers, such as high-density polyethylene (HDPE), hydroxyl-terminated polybutadiene (HTPB), and polymethylmethacrylate (PMMA), or polymers with metal additives to improve the density impulse [5]. The fuel regression of these polymers is determined by the ratio between the heat flux to the surface and the heat of the phase change, thus, it is limited by the heat and mass transfer mechanisms occurring from the flame to the fuel wall. The blowing of the fuel from the surface decreases the velocity gradient at the wall and the convective heat transfer for the so-called blocking e ffect [6]. Owing to this sort of "counter-balance" between heat flux and blowing, the regression rate marginally depends on the nature of classical polymeric fuels [6], and is usually relatively small; consequently, the volumetric fuel loading e fficiency may be too poor for space applications.

Compared to conventional polymers, the consumption mechanism of para ffin-based fuels is di fferent and allows for a significantly larger regression rate. Karabeyoglu et al. [7] have shown that these fuels display regression rates up to 3–4 times higher than those achieved with traditional hybrid fuels. For this reason, a lot of studies have been focused on the type of materials that belong to the class of the so-called liquefying fuels [8]. Their intrinsic characteristic is the onset of a thin liquid layer on the fuel grain surface, which may become unstable at low viscosity and surface tension [9,10], leading to the lift-o ff and entrainment of fuel liquid droplets or filaments into the main gas stream, increasing the fuel mass transfer rate. This mass-transfer mechanism essentially does not depend on the heat transfer to the surface and raises the fuel mass flow without entailing the consequent blocking on gaseous fuel blowing. As a result, the overall regression rate can be considered to be composed of two terms, one determined by classical fuel vaporization and the other by liquid entrainment. The entrainment phenomenon is strongly susceptible to fuel composition, its manufacturing process and the obtained thermo-mechanical properties, as well as to the engine operating conditions [11], which makes the prediction of the combustion chamber internal ballistics even harder than in the case of non-liquefying polymers. Hence, on the one hand, designers need to characterize the fuel with extended experimental campaigns [12,13] and, on the other, to carry out rocket static firings to measure the achieved engine performance.

Predicting the local, time-resolved fuel regression rate is one major challenge in the design of a hybrid rocket. A number of models have been developed over the years to describe the hybrid combustion mechanism, but often they lack some important aspects and/or fail in the prediction of the experimental results [14,15]. The first hybrid combustion theory accounting for the main underlying phenomena involved in the burning of classical polymeric fuels was developed in 1963 by Marxman and Gilbert [6,16] and it still remains the starting point of design calculations and experimental comparisons. Subsequent improvements [17] are, however, all based on the same assumption of a turbulent boundary layer with chemical reactions occurring over a fuel slab burning in an oxidant gas flow. Because of this fundamental hypothesis, those classical theories are unable to reproduce the effects of complex flows, such as the recirculation ensuing from the oxidizer injection, which may have a non-negligible impact even in standard-flow motors [18,19]. Analytical models later developed for liquefying fuels [7,20] are essentially modifications of the classical hybrid boundary-layer combustion theory for the entrainment mass transfer from the fuel grain and, consequently, yield the same limits as the original theory.

In this framework, the numerical modeling of the rocket internal thermo-fluid-dynamics, with predictive capabilities of the fuel regression rate and overall engine performance, is becoming a key tool both in the system design process and in the experimentally measured performance-analysis stage. In fact, a ffordable computational models provide a quick detailed representation of the phenomena governing the engine internal ballistics, on the one hand allowing for numerous motor optimization trials and, on the other, reducing the need of expensive experimental testing.

This paper presents the computational fluid dynamic (CFD) models defined by the authors for simulating the internal ballistics of a hybrid rocket burning either classical polymeric or para ffin-based fuels with gaseous oxygen, with an eye to their evolution. The presentation starts with an overview of the numerical techniques available in the competent literature and develops through the steps taken up to the current level. Numerical results are summarized and compared to the experimental data gathered by static firings of two laboratory-scale engines.

#### **2. State of the Art of CFD Techniques for Rocket Internal-Ballistics Simulation**

CFD approaches to the solution of the flowfield in the hybrid-propellant rocket chamber have been considerably developed recently [21–23], but comprehensive models to describe the complex interactions among fluid dynamics [24], solid fuel pyrolysis [25] or melting and entrainment [7], oxidizer atomization and vaporization, mixing and combustion in the gas phase [26], nozzle thermochemical erosion [27], particulate formation, and radiative characteristics of the flame [28] are still lacking and numerical simulations are often considered only as a qualitative tool to a fford the thermo-fluid-dynamics of the rocket. Most of the e ffort has been addressed, so far, to classical non-liquefying fuels. A common strategy is solving the Reynolds averaged Navier-Stokes (RANS) equations, with suitable turbulence closure and combustion models. The latter models are essentially split into two categories: finite-rate kinetics and mixing-limited combustion mechanisms. When the former philosophy is followed, usually gas phase reactions are modeled by global reaction mechanisms because detailed chemical kinetics would include many species requiring a huge computational cost and, based on the experience gained by the authors, would probably not e ffectively improve the gas-surface interaction simulation [29,30]. Whereas, following the approach that chemistry is infinitely fast, the overall rate of the reaction is controlled by turbulent mixing and the combustion can be described, for instance, with the eddy dissipation model [31,32] or with chemical equilibrium and an interaction model between turbulence and chemistry [33].

Justified by the fact that the chemical and fluid-dynamic characteristic times are much shorter than the fuel-regression time scale, the steady-state solution of RANS equations is generally sought [34]. An acceptable method to study the hybrid rocket internal ballistics can, therefore, be simulating the flowfield at di fferent times in the motor firing by updating the fuel port diameter [26]. Nevertheless, a single numerical simulation is often performed on the chamber geometry drawn at the time-space averaged port diameter [29,30]. The cases in which the port section longitudinal non-uniformity is considered are rare [35].

Depending on the specific target of the simulation, several techniques can be implemented; the fuel regression rate can be, on the one hand, an input to the numerical simulation, for instance, when one is interested in studying the propellant mixing and combustion processes in the motor [31] compared to experimental data; or, maybe in a pre-screening stage, to make trade-o ffs among several oxidizer injections and/or motor configurations as described by the authors in Reference [33]. In this latter case, the simple mass flow inlet boundary at the fuel wall is su fficient and the problem is simplified. On the other hand, the fuel regression rate is an outcome of the simulation and the proper fuel-surface boundary conditions are required to allow for its calculation.

In particular, with classical polymeric fuels, an energy balance at the solid fuel interface with gas, coupled with a semiempirical pyrolysis-rate relationship between the regression rate and surface temperature, is used. Typically, the regression rate is modeled with a one-step irreversible Arrhenius-type equation. Given a polymeric fuel, the laws of pyrolysis available in the open literature

are widely scattered because of significantly di fferent pre-exponential factors and activation energies. However, by means of a sensitivity analysis conducted with HDPE, in Reference [36], it has been demonstrated that there is no significant dependence of the calculated regression rate from the law of pyrolysis adopted, provided that each variable is consistently taken from a single set of data.

The computational domain can include the actual prechamber geometry or just a simplification of it [29,30] based on the assumption that, with a gaseous oxidizer, the real dump-plenum arrangemen<sup>t</sup> does not have a significant impact on the results. Additionally, the postchamber can be assimilated to a cylinder that is equal to the grain port if one is not concerned with the details of mixing the upstream of the nozzle. In other cases, one may exclude the motor aft-mixing chamber and exhaust nozzle if the chamber pressure and combustion e fficiency are not under investigation [35]. In that case, the pressure has to be an input to the problem and a pressure-outlet boundary condition is usually imposed.

When considering the energy balance at the wall for the regression rate calculation, other than the convective heat transfer, the heat exchanged by radiation with the combustion gases may play a non-negligible role for some polymers [36]; however, a clear conclusion has not been drawn yet.

The definition of a suitable and computational cost-e ffective strategy for liquefying fuels poses further complications because, in principle, to successfully simulate the para ffin-fuel consumption, two non-trivial tasks have to be accomplished. These are modeling, first, the formation, the breaking up and the entrainment of the liquid film on the melting fuel surface; and, second, the transformation of the melted fuel into a gaseous species participating in the combustion process. These demanding efforts have probably discouraged researchers, so drastic simplifications are usually introduced, such as giving the regression rate calculation away by assuming it from experiments [31,37–39], or limiting the analysis to one-dimensional integral-di fferential models [40]. In other cases, by observing that the melted para ffin wax is in the supercritical state under the hybrid rocket chamber characteristic conditions (thus, surface tension vanishes and the sharp distinction at droplets surface between gas and liquid phases disappear), the melted layer break up and, subsequently, the liquid para ffi n injection in the flowfield is disregarded, supposing that the entrainment is part of the turbulent mixing process [41,42]. However, those models are not successfully validated, still displaying significant deviations from the experimental data which, in some cases, are around 25% [42].

In the forthcoming sections, the physical and numerical models developed by the authors for the calculation of the regression rate of both classical and liquefying fuels are presented, and a summary of the main results obtained in several test cases is discussed with the complement of experimental data. All the numerical simulations are carried out with a commercial fluid dynamic solver with ad-hoc user-defined functions. The base numerical framework is firstly shown and, then, the details of the di fferent wall-boundary managemen<sup>t</sup> are reported.

#### **3. Physical and Numerical Models**

#### *3.1. Governing Equations*

The RANS equations for single-phase multicomponent turbulent reacting flows are solved; note that when dealing with liquefying fuels, the main underlying hypothesis here is that, being liquid para ffin in supercritical conditions [42], its viscosity and di ffusivity are close to those of a gas and the di ffusion processes are significantly faster than in the liquid phase, which, in a first approximation, may allow for neglecting the two-phase flow. Hence, in both the cases of polymers and para ffin wax burning, the fuel is supposed to enter the combustion chamber as 100% gaseous ethylene (which is expected to be the main product of both the HDPE thermal pyrolysis and para ffin-wax decomposition).

The Favre-averaged (i.e., density-weighted) equations of continuity and momentum can be expressed in Cartesian tensor form, with an understanding that repeated indices mean summation, as [43]

$$\frac{\partial \overline{\rho}}{\partial t} + \frac{\partial}{\partial x\_j} (\overline{\rho} \overline{u}\_j) = S\_m \tag{1}$$

*Aerospace* **2019**, *6*, 56

$$\frac{\partial}{\partial t}(\widetilde{\rho u\_i}) + \frac{\partial}{\partial \mathbf{x}\_j}(\widetilde{\rho u\_i u\_j}) = -\frac{\partial \overline{p}}{\partial \mathbf{x}\_i} + \frac{\partial \overline{\tau}\_{i\bar{j}}}{\partial \mathbf{x}\_j} + \frac{\partial}{\partial \mathbf{x}\_j}(-\overline{\rho u\_i' u\_j'}) \tag{2}$$

where *Sm* is the mass source term, which, as will be explained later, is used to model the liquid fuel mass entrainment occurring with para ffin wax. Here, the bar denotes conventional time averaging, while the tilde denotes density-weighted averaging; <sup>τ</sup>*ij* is the stress tensor that is defined as

$$\pi\_{i\bar{j}} = \mu \left[ \left( \frac{\partial u\_i}{\partial \mathbf{x}\_{\bar{j}}} + \frac{\partial u\_{\bar{j}}}{\partial \mathbf{x}\_{\bar{i}}} \right) - \frac{2}{3} \delta\_{i\bar{j}} \frac{\partial u\_k}{\partial \mathbf{x}\_k} \right] \tag{3}$$

where δ*ij* is the Kronecker delta. Symbols with a prime indicate the corresponding quantity fluctuation. The term <sup>R</sup>*ij* = −ρ*<sup>u</sup> i u j*, originating from the averaging operation, is known as the Reynolds stress tensor, and it needs to be modeled.

The Shear Stress Transport (SST) turbulence model [44] has been employed for its improved capability of predicting flows with separated regions. This is a combination of the robust and accurate *k*–<sup>ω</sup> model, developed by Wilcox [45] in the near-wall region, with the standard *k*–<sup>ε</sup> model implemented away from the wall using a blending function. With the SST model, the transport equations of the turbulence kinetic energy, *k*, and the specific dissipation rate, ω, are formulated as

$$\frac{\partial}{\partial t}(\overline{\rho}k) + \frac{\partial}{\partial \mathbf{x}\_i}(\overline{\rho}\overline{u}\_i k) = \frac{\partial}{\partial \mathbf{x}\_j} \Big[ (\mu + \mu\_l \sigma\_k) \frac{\partial k}{\partial \mathbf{x}\_j} \Big] + \mathcal{R}\_{i\bar{j}} \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_j} - \beta^\* \overline{\rho} \alpha k \tag{4}$$

$$\frac{\partial}{\partial t}(\overline{\rho}\omega) + \frac{\partial}{\partial \mathbf{x}\_l}(\overline{\rho}\overline{\mu}\_l\omega) = \frac{\partial}{\partial \mathbf{x}\_j} \Big[ (\mu + \mu\_l \sigma\_{\omega}) \frac{\partial \omega}{\partial \mathbf{x}\_j} \Big] + \overline{\rho}\frac{a}{\mu\_l} \Re\_{\overline{\eta}} \frac{\partial \overline{u\_l}}{\partial \mathbf{x}\_j} - \beta \overline{\rho}\omega^2 + 2(1 - F\_1)\overline{\rho}\sigma\_{\omega\_2} \frac{1}{\omega} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_j} \frac{\partial \omega}{\partial \mathbf{x}\_j} \tag{5}$$

in which the Reynolds stress is modeled using the Boussinesq approximation

$$\mathcal{R}\_{ij} = \mu\_l \left[ \left( \frac{\partial \widetilde{u}\_l}{\partial \mathbf{x}\_j} + \frac{\partial \overline{u}\_j}{\partial \mathbf{x}\_i} \right) - \frac{2}{3} \delta\_{ij} \frac{\partial \widetilde{u}\_l}{\partial \mathbf{x}\_l} \right] - \frac{2}{3} \overline{\rho} k \delta\_{ij} \tag{6}$$

The turbulent viscosity, μ*<sup>t</sup>*, is expressed as follows

$$
\mu\_t = \frac{\overline{\rho}k}{\omega} \frac{1}{\max\{1; \frac{\Omega F\_2}{0.31\omega}\}}\tag{7}
$$

where the function *F*2 is defined depending on the distance from the wall *y*, as

$$F\_2 = \tanh(\Phi\_2^2) \tag{8}$$

with

$$\Phi\_2 = \max\left(\frac{2\sqrt{k}}{0.09\omega y}; \frac{500\mu}{\overline{\rho}\omega y^2}\right) \tag{9}$$

The coe fficient α is given by

$$\alpha = \gamma \frac{1/9 + \text{Re}\_l/2.95}{1 + \text{Re}\_l/2.95} \tag{10}$$

where *Ret* = ρ*k*/μω is the turbulent Reynolds number. The blending function *F*1 takes the value of 1 on the wall and tends to zero at the boundary layer edge, being defined as

$$F\_1 = \tanh\left(\Phi\_1^4\right) \tag{11}$$

with

$$\Phi\_1 = \min\left[\max\left(\frac{\sqrt{k}}{0.09\alpha y}; \frac{500\mu}{\overline{\rho}\alpha y^2}\right); \frac{4\overline{\rho}\sigma\_{\alpha\alpha}k}{\overline{\sigma}D\_{k\omega}y^2}\right] \tag{12}$$

where *CDk*ω is the positive part of the last term in Equation (5) (cross-diffusion term):

$$\text{CD}\_{k\omega} = \max \left( 2\overline{\rho}\sigma\_{\omega\_2} \frac{1}{\omega} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_j} \frac{\partial \omega}{\partial \mathbf{x}\_j}; 10^{-20} \right) \tag{13}$$

The model coefficients σ*k*, σω, β, γ are defined by blending the corresponding coefficients of the original *k*–<sup>ω</sup> model, denoted with the subscript 1, with those of the transformed *k*–<sup>ε</sup> model that are denoted with the subscript 2:

$$
\begin{bmatrix} \sigma\_k \\ \sigma\_{a\nu} \\ \beta \\ \gamma \end{bmatrix} = F\_1 \begin{bmatrix} \sigma\_{k\_1} \\ \sigma\_{a\nu\_1} \\ \beta\_1 \\ \gamma\_1 \end{bmatrix} + (1 - F\_1) \begin{bmatrix} \sigma\_{k\_2} \\ \sigma\_{a\nu\_2} \\ \beta\_2 \\ \gamma\_2 \end{bmatrix} \tag{14}
$$

All the model constants are listed in Table 1.

**Table 1.** The values of the Shear Stress Transport (SST) model constants [43].


Assuming that the chemical kinetics are fast compared to the diffusion processes occurring in the motor for the typical mass fluxes and chamber pressures considered here [46], the non-premixed combustion of oxygen and gaseous ethylene is modeled by means of the Probability Density Function (PDF) approach coupled with chemical equilibrium [47]. Accordingly, the combustion is simplified to a mixing problem and the difficulties associated with closing non-linear reaction rates are avoided. In fact, under the hypothesis of equal diffusivities for all chemical species and assuming that the Lewis number is equal to 1, the species equations can be reduced to a single equation for the transport of the mixture fraction, which represents the sum of the element mass fractions contained in the fuel stream [48], *f* = 1/(1 + *OF*), where *OF* is the local oxidizer-to-fuel mass ratio for the equivalent non-burning field. The density-averaged mixture fraction equation is

$$\frac{\partial}{\partial t} \left( \overline{\rho} \,\widetilde{f} \right) + \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} \left( \overline{\rho} \overline{u}\_{\dot{j}} \,\widetilde{f} \right) = \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} \left( \frac{\mu\_{t}}{\mathrm{Sc}\_{t}} \frac{\partial \overline{f}}{\partial \mathbf{x}\_{\dot{j}}} \right) + S\_{\text{m}} \tag{15}$$

where *Sct* is the turbulent Schmidt number, which, for the hypothesis of the unity Lewis number, is equal to the turbulent Prandtl number, *Prt*; the latter is assumed to be equal to 0.85.

For the closure model describing the turbulence-chemistry interaction, the variance of the mean mixture fraction *f*<sup>2</sup> is introduced and an additional equation for this quantity is needed, which, according to Reference [49], and making use of the relation between ω, *k*, and ε, is written as

$$\frac{\partial}{\partial t} \left( \overline{\rho} \overline{f'^2} \right) + \frac{\partial}{\partial \mathbf{x}\_j} \left( \overline{\rho u\_j} \overline{f'^2} \right) = \frac{\partial}{\partial \mathbf{x}\_j} \left( \frac{\mu\_t}{Pr\_l} \frac{\partial \overline{f'^2}}{\partial \mathbf{x}\_j} \right) + 2 \frac{\mu\_t}{Pr\_l} \left( \frac{\partial \overline{f}}{\partial \mathbf{x}\_j} \right)^2 - 2 \beta^\* \overline{\rho} \omega \overline{f'^2} \tag{16}$$

The assumed shape of the PDF is based on the beta distribution [50]:

$$\mathbf{B}(f) = \frac{f^{a-1}(1-f)^{b-1}}{\int\_0^1 \mathbf{x}^{a-1} (1-\mathbf{x})^{b-1} d\mathbf{x}} \tag{17}$$

in which the two parameters *a* and *b* are functions of the mean mixture fraction and its variance

$$a = \overline{f} \left[ \frac{\overline{f} \left( 1 - \overline{f} \right)}{\overline{f'^2}} - 1 \right] \tag{18}$$

$$b = \left(1 - \overline{f}\right) \left[\frac{\overline{f}\left(1 - \overline{f}\right)}{\overline{f'^2}} - 1\right] \tag{19}$$

Finally, the energy equation needs to be solved because the system is non-adiabatic (the heat is exchanged at the fuel surface), which has an impact on the temperature and species of the reacting flow at the chemical equilibrium. The enthalpy form of the energy equation is selected for the obvious benefit connected with the equilibrium calculations; by neglecting the term of viscous dissipation (the Brinkman number, *Br* = *<sup>U</sup>*<sup>2</sup>*Prt*/Δ*Hw*, defined as the ratio between heat produced by viscous dissipation and heat transported by turbulent conduction at the wall, is on the order of 5 × <sup>10</sup>−3) as well as the variation of pressure (which is approximately constant all over the combustion chamber except through the discharge nozzle where, however, adiabatic walls are imposed), and combining the conduction and species diffusion terms, the latter can be written as

$$\frac{\partial}{\partial t} \left( \overline{\rho} \overline{H} \right) + \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} \left( \overline{\rho} \overline{\mu}\_{\dot{j}} \overline{H} \right) = \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}} \left( \frac{\mu\_{l}}{Pr\_{l}} \frac{\partial \overline{H}}{\partial \mathbf{x}\_{\dot{j}}} \right) + S\_{l\mathbf{r}} \tag{20}$$

where the source term *Sh* includes the volumetric heat of phase change. Here, the enthalpy of the fluid mixture is defined as

$$H = \sum\_{j} \mathbf{Y}\_{j} \Big[ H\_{j}^{0} \Big( T\_{ref,j} \Big) + \int\_{T\_{ref,j}}^{T} c\_{p,j}(T) dT \Big] \tag{21}$$

where *Yj* is the mass fraction of the *j*-th species and *<sup>H</sup>*0*jTre f*,*j*is its formation enthalpy at the reference temperature *Tre f*,*j*.

Density-weighted mean temperature and individual species mole fractions are then obtained as functions of *f* , *f*<sup>2</sup> and *H* as yielded by the following equation

$$\widetilde{\varphi} = \int\_0^1 \varphi(f, \overleftarrow{H}) \mathbf{B}(f) df \tag{22}$$

where <sup>ϕ</sup>*<sup>f</sup>*, *H*represents the generic functional dependency of one of the above-mentioned parameters on the mixture fraction and enthalpy, whereas the mean density is calculated as

$$\frac{1}{\overline{\mathcal{P}}} = \int\_0^1 \frac{\mathcal{B}(f)}{\rho(f, \widehat{H})} df \tag{23}$$

Note that the enthalpy turbulent fluctuations are assumed to be independent of the enthalpy level so that the selected PDF in Equation (17) does not change.

Chemical equilibrium conditions are calculated via the minimization of Gibbs free energy algorithm [51] in a number of *f* , *f*<sup>2</sup> and *H* points at a specified value of the chamber pressure (which here has been assumed equal to that measured), and a lookup table is constructed. Thus, once the mean mixture fraction, its variance and mean enthalpy are calculated at each point in the flowfield, the corresponding time-averaged species mole fractions, temperature and density are obtained by interpolating the values in the lookup table. The densities are, then, scaled with the values of the actual pressure field in the system.

Heat capacities, molecular weights, and enthalpies of formation for each species considered are extracted from the solver chemical database; the specific heat is determined via the mixing law. Molecular dynamic viscosities and thermal conductivities of each species are calculated as functions of local temperature, according to Reference [51].

#### *3.2. Computational Domain and Boundary Conditions*

Computational results presented in the following table refer to experimental firing tests performed with two laboratory-scale hybrid rocket engines, one in the 200-N thrust scale and the other in the 1-kN class. The two engines have a similar design. They have axisymmetric combustion chambers. A conical nozzle injector, whose exit-section diameter is 6 mm in the case of the smaller engine and 8 mm in the case of the larger one, has been employed to deliver gaseous oxygen into a single-port cylindrical fuel grain. Upstream and downstream of the solid grain, a dump plenum and an aft-mixing chamber are set up respectively. In both configurations, the combustion products are accelerated through converging-diverging discharge nozzles. The main dimensions are summarized in Table 2.


**Table 2.** The computational domains' characteristic dimensions.

For the flowfield axial-symmetry, numerical computations are conducted with two-dimensional grids built in the internal volume of the combustion chamber components; examples of the grids employed with two different scale engines are shown in Figure 1.

The prechamber of the smaller scale engine (Figure 1a) is subdivided into 40 × 90 grid cells in the axial and radial directions, respectively; the grain port into 240 × 40 grid cells, the postchamber into 80 × 90 cells and the nozzle into 60 × 40 cells; whereas in the larger scale engine (Figure 1b), the prechamber and the postchamber are both subdivided into 60 × 110 cells, the grain port into 310 × 50 grid cells and the nozzle into 140 × 50 cells. Cells are clustered towards the grain wall in such a way to ensure that the maximum value of *y*+ is around 1 ÷ 2 at the wall-adjacent cell all along the grain length. Additional axial clustering of cells is placed in the regions near the grain inlet and outlet edges, and near the prechamber, postchamber and nozzle inner surfaces. Grid sensitivity analyses assessing the grid-convergence of the numerical results are reported in Reference [52].

On the inner surface of both the prechamber and postchamber as well as on the nozzle wall, no-slip and adiabatic boundary conditions are imposed. At the injector exit section, a mass flow boundary condition is prescribed along with the temperature (equal to 300 K), the oxygen mass fraction and the turbulent quantities, and a supersonic outflow condition is set at the nozzle exit section for which all the variables are extrapolated from the domain interior.

(**b**) 1-kN class hybrid rocket engine. 

**Figure 1.** Examples of computational grid construction.

#### *3.3. Gas-Solid Fuel Interface Modeling*

The theoretical model formulation is completed by assigning the boundary conditions at the interface between the gaseous flow region and the solid fuel wall, which can properly describe the fuel consumption mechanism. The fuel surface is an inlet boundary along which the fuel mass flux, the temperature and the mixture fraction depend on the regression rate that is an unknown to be determined.

Now it is necessary to separate the case of classical polymers from the one of liquefying fuels. In fact, while regardless of the type of fuel, the gas/surface interface treatment relies on local mass, energy and mean mixture-fraction balances. In the former case, an additional pyrolysis rate equation is needed for modeling the regression rate, whereas, in the latter case, a different formulation is introduced to account for the entrainment of liquid paraffin from the unstable melt layer forming along the fuel surface, which is the dominant consumption mechanism.

#### 3.3.1. Classical Non-Liquefying Fuels

In the case of pyrolyzing fuels, since no material is removed from the surface in a condensed phase (neither solid, such as in the case of fuel loaded with metal particles; nor liquid, as in the case of paraffin wax analyzed later), the mass conservation at the gas-solid interface imposes that

$$(\rho v)\_w = \rho\_f \dot{r} \tag{24}$$

where ρ is the gas density at the wall and *v* is the normal-to-wall velocity component due to the pyrolysis products injection; ρ*f* is the solid fuel density and .*r* is the local regression rate.

The energy balance at the gas-solid interface, taking into account the convective heat transfer from the gas to the fuel surface, the steady heat conduction into the solid, and neglecting the radiation, leads to the following relationship between the convective heat flux to the wall, .*qw*, and the regression rate [29]

$$\dot{q}\_w = \left(k\_\% \frac{\partial T}{\partial n}\right)\_w = \rho\_f \left\dot{r} \left[\Delta l\_p + \mathbb{C}\_s (T\_{w\prime} - T\_a)\right] \tag{25}$$

where *n* is the coordinate normal to the surface oriented from solid to gas, *kg* is the gas thermal conductivity, *Cs* is the solid heat capacity per unit mass, Δ*hp* is the so-called heat of pyrolysis, *Tw* is the fuel surface temperature, and *Ta* is its initial temperature (which is assumed to be equal to the one of the external surface of the fuel). The term in brackets at the right-hand side represents the e ffective heat of gasification of the fuel, which accounts for the heat radially conducted into the solid grain, further than for the heat of pyrolysis. Note that Equation (25), by properly defining the heat of pyrolysis, implicitly takes into account the heat transfer to the wall by the species di ffusion [30].

The fuel pyrolysis is, finally, modeled with the following Arrhenius-type equation [25]:

$$\dot{r} = A \cdot \exp\left(-\frac{E\_d}{2RT\_w}\right) \tag{26}$$

where *A* is the pre-exponential factor, *Ea* is the activation energy and *R* is the universal gas constant.

The values of the constants appearing in Equations (25) and (26) considered for the HDPE fuel grains analyzed in this work are summarized in Table 3. The density, specific heat and heat of pyrolysis are taken from the work in Reference [53], while the values of the pre-exponential factor and the activation energy from Reference [25] are obtained by modifying the activation energy to match the surface temperature commonly observed in polymeric hybrid fuels (which is around 800 K) [54].

**Table 3.** The high-density polyethylene (HDPE) fuel properties and rate constants.


A dedicated treatment of the mean mixture-fraction boundary condition at the fuel wall is needed as well. The fuel regression rate is typically low in hybrid rockets; therefore, on the one hand, the normal convection of the fuel at the grain surface is relatively weak compared to the gas convection in the cells near the boundary. On the other hand, the di ffusive flux plays a dominant role in the mixture-fraction transport near the grain surface where the species concentrations rapidly change, so that a steep mixture-fraction gradient at the fuel wall is present. As a consequence, the simple Dirichlet boundary condition, *fw* = 1, on the gas-fuel interface (which is equivalent to imposing that the di ffusion coe fficient is equal to zero in the cells close to the fuel inlet boundary) is not adequate. In fact, first, it would imply a non-exact evaluation of the gradients in this zone, and, in particular, of the heat flux to the wall (which, for Equation (25), would produce a flawed regression rate), and, second, an extra mixture fraction would be di ffused into the flow a ffecting the global oxidizer to fuel ratio and the chemical equilibrium properties, which eventually would lead to an incorrect estimation of the characteristic exhaust velocity and chamber pressure. The correct solution to this problem is considering an additional equation for the mean mixture-fraction balance at the gas-solid interface [55], which, by rearranging Equation (15) at the wall, can be expressed as

$$(\rho v)\_w f\_w - \left(\frac{\mu\_t}{Sc\_t} \frac{\partial f}{\partial n}\right)\_w = \rho\_f \,\dot{r} \tag{27}$$

According to Equation (27), the total mass flux entering the gaseous domain due to the solid fuel regression, which appears on the right-hand side of the equation and represents the production term, is partially balanced by the convection and partially by the di ffusion of the fuel mass fraction.

The balance of the mixture-fraction variance, *f*2, at the wall is instead ignored and it is imposed to be zero.

Note that Equations (25)–(27), upon substitution of Equation (24) into Equation (27) for the wall mass flux, constitute a system of three algebraic equations for the three unknowns—regression rate, the surface temperature and the mixture fraction at the wall—which need the computation of the flowfield at each iterative step to be solved.

The required level of mesh refinement near the grain surface allows for the resolution of the viscous sub-layer, with the following boundary conditions for the turbulent kinetic energy and the specific dissipation rate, respectively:

$$\frac{\partial k}{\partial n} = 0, \; \omega\_{\!w \!v\!=} = \frac{6\rho}{\mu\beta} \left( \frac{\mu^\*}{\mathcal{Y}^+} \right)^2 \tag{28}$$

where *u*<sup>∗</sup> = \$<sup>τ</sup>*w*/<sup>ρ</sup> is the friction velocity.

> Both the source terms appearing in Equations (1), (15) and (20) are identically zero:

$$S\_{\mathfrak{m}} = S\_{\mathfrak{h}} = 0 \tag{29}$$

As we will see in the forthcoming section, they are instead calculated in the modeling of the paraffin-wax fuel combustion.

#### 3.3.2. Liquefying Fuels

The boundary conditions for the turbulent quantities applied in this case are equal to those described above for polymeric fuels.

The regression rate, .*r*, can be, now, assumed composed of two terms: the vaporization fraction, .*rv*, that is generated by the liquid thermal decomposition and later vaporization into the gas stream, and the entrainment fraction, .*rent*, that is related to the mechanical transfer of the liquid from the surface melt layer

$$
\dot{r} = \dot{r}\_{\upsilon} + \dot{r}\_{\text{ent}} \tag{30}
$$

A scheme of this peculiar fuel consumption mechanism is shown in Figure 2 in the case of the supercritical regime: part of the molten fuel on the solid surface undergoes pyrolysis and part leaves the surface in the form of a supercritical fluid that is entrained in the gas stream and burns farther from the wall. More precisely, the molten paraffin wax will behave like a dense fluid which, for simplicity, we call "liquid" and assume having the same properties as paraffin in the liquid state.

**Figure 2.** The liquefying fuel consumption mechanism schematic.

A set of equations is, thus, needed for the calculation of the regression rate and its two components, along with the resulting fuel mass flow rates, whose solution is to be incorporated in the fluid dynamic computation. As in the previous case, the mass, energy and mixture-fraction balances at the gas/fuel surface boundary are formulated, with a difference that an additional equation for the calculation of the entrainment component of the fuel mass flow rate is required.

Following the arguments in References [7,54], by coupling the energy balance equations at both the liquid-solid and gas-liquid interfaces (see Figure 2), at the steady state, the following relationship of the total surface heat flux with the total and the vaporization regression rate is obtained

$$\dot{q}\_{\rm av} = \left( k\_{\rm g} \frac{\partial T}{\partial n} \right)\_{\rm w} = \rho\_f \dot{r} [\mathbb{C}\_s (T\_m - T\_a) + L\_m + \mathbb{C}\_l (T\_w - T\_m)] + \rho\_f \dot{r}\_v L\_v' \tag{31}$$

where ρ*f* is the solid fuel density, *Cs* and *Cl* are the specific heats of the solid and liquid fuels (which are considered independent from the temperature here), respectively, *Tm* is the fuel melting temperature, and *Lm* and *L v* are the fuel heat of the fusion and the heat of pyrolysis, respectively.

Equation (31) represents the fact that the total heat flux transferred from the combusting gases to the fuel surface must be equal to the sum of heat conducted into the liquid layer and the heat required for the pyrolysis of the fuel vaporized fraction (the overall term in the square brackets and last term on the right-hand side of Equation (31), respectively). Note that, for the finite thickness of the liquid layer, the conductive heat flux at the liquid-gas interface is not equal to ρ*f* .*rCl*(*Tw* − *Tm*), which one might accidentally expect by analogy with the heat conducted in the solid appearing in Equation (25).

In the subcritical regime, the wall phenomena are governed by the evaporation process and *L v* = *Lv* (*Lv* being the heat of vaporization) and the wall temperature is equal to the vaporization temperature; whereas, in supercritical conditions, pyrolysis takes place on the surface. In the absence of clear para ffin pyrolysis data, for the sake of simplicity, and also in the supercritical case, the heat of vaporization reported in Reference [7] has been employed by neglecting the heat required for the thermal degradation of para ffin into gaseous ethylene monomers. The wall temperature has a significant influence on the fuel regression rate as it a ffects both the heat flux and the term *Cl*(*Tw* − *Tm*) appearing in the wall energy balance, Equation (31). Compared to a purely pyrolyzing polymer, due to the e ffect of the entrainment, the surface temperature is reduced [54]; however, this parameter can be hardly determined and. in Reference [56], a sensitivity analysis was performed for which the value of 675 K is the one allowing for the best fit of the experimental data. An isothermal boundary condition is set all over the fuel wall.

It is worth noting that regardless of the definition of the wall temperature and enthalpy (either vaporization or pyrolysis) entailing the vaporization fraction of the regression rate, Equation (31) produces a total regression rate that is marginally sensitive to the value of the enthalpy because. If the vaporization regression rate is decreased for a larger vaporization enthalpy, the total heat flux to the wall tends to rise and so does the entrainment fraction. This sort of balance is further explained later.

Material constants appearing in Equation (31) are listed in Table 4.

> **Table 4.** Para


ffin fuel properties.

According to the approach described in Reference [7], the following semiempirical relationship has been considered for modeling the entrainment component of the fuel regression rate as a function of the rocket chamber operating conditions and of an entrainment parameter that lumps the liquid fuel properties together:

$$
\dot{r}\_{cnt} = a\_{cnt} \frac{G^3}{\dot{r}^{1.5}} \tag{32}
$$

where *G* = 4 . *m*/π *D*<sup>2</sup> is the total mass flux in the local section of the grain port and *aent* is the entrainment factor depending on the physical properties of the selected fuel, primarily on the fuel liquid viscosity, and on the average gas density in the chamber as

$$a\_{\rm cntt} \propto \frac{1}{\mu\_l \rho\_\mathcal{S}^{1.5}} \tag{33}$$

Equation (33) is derived from a theoretical assessment of the surface liquid-layer fluid dynamic stability; the main result is that the susceptibility of a given fuel to the instability increases with decreasing viscosity and surface tension of the melt layer; the entrainment component of fuel regression rate is, therefore, roughly inversely proportional to viscosity at the characteristic temperature of the layer, while it depends directly on the dynamic pressure. A parametric analysis of the e ffect of the entrainment parameter on the regression rate components is reported in Reference [56]. From this analysis, the value of 2.1 × 10−<sup>13</sup> m8.5s0.5/kg<sup>3</sup> has been identified for the best fit of the experimental data in the reference Test P4. The latter value, considering the di fferent densities (average gas density in Test P4 is 1.62 kg/m3) is in good agreemen<sup>t</sup> with the one reported in Reference [7].

Once Equations (30)–(32) are combined, given the heat flux to the wall and the total mass flux, the three components of the fuel regression rate can be calculated. The fuel mass fluxes associated with the vaporization and entrainment components, respectively, are obtained as follows

$$G\_{f,v} = \rho\_f \dot{r}\_v \tag{34}$$

$$G\_{f, \text{ent}} = \rho\_f \dot{r}\_{\text{ent}} \tag{35}$$

Vaporization and entrainment components are handled di fferently.

The vaporization component is treated equally to the case of pyrolyzing fuels, considering the mass and mixture-fraction balance equations at the grain wall, given by

$$(\rho v)\_w = G\_{f,v} \tag{36}$$

$$(\rho v)\_w f\_w - \left(\frac{\mu\_t}{Sc\_t} \frac{\partial f}{\partial n}\right)\_w = G\_{f,v} \tag{37}$$

This allows for correctly taking the blocking e ffect on the heat transfer to the wall into account, while, as explained in the previous section, Equation (37) is needed to ensure the mixture fraction has a global balance.

The entrainment mass flux does not contribute to the blocking e ffect, thus, a specific treatment is adopted for the introduction of the entrainment component into the computational domain. For the sake of simplicity, we assume that despite the entrained para ffin initially in the liquid phase immediately gasifies because of the large combustion heat release. The local entrainment contribution is uniformly assigned as a mass production term in the local volume of the grain port corresponding to the surface cell of length Δ*x* through which the fuel mass enters the fluid domain, π *<sup>D</sup>*2Δ*x*/4:

$$S\_m = 4 \frac{G\_{f,out}}{D} \tag{38}$$

In order to satisfy the species balance, an equal production term is assigned also for the mean mixture fraction. Finally, the energy required by the pyrolysis of the liquid fuel mass flow rate is taken into account by assigning a corresponding negative energy source term in the same volume:

$$S\_{\rm li} = -4 \frac{G\_{f,cut}}{D} L'\_{\rm v} \tag{39}$$

Additionally, in this case, as the heat flux to the surface and the total mass flux needed for the calculation of the regression rate are outputs of the flowfield resolution, which, in turn, depends on the regression rate itself, an iterative procedure is needed for the problem solution.

#### *3.4. Solution Strategy*

In this context, steady-state solutions to the equations presented above are searched, even in the case of transient simulations, as explained in Section 3.5. The equations are solved with a control-volume-based technique and a pressure-based solver [57], i.e., Equations (1) and (2) are not solved directly, but are combined to derive a pressure equation in such a way that the velocity field, corrected by the pressure satisfies the continuity equation. The velocity and pressure fields are simultaneously solved with a coupled algorithm. The di ffusion terms are central-di fferenced and are second-order accurate; the convective terms are represented according to a quadratic upwind scheme. Convergence is obtained by iterating until the residuals drop by five orders of magnitude.

The solution procedure is approached through the following steps.

(1) Input the boundary conditions and a combustion pressure reference value; at the solid grain wall, trial values of the unknown required variables are used.

(2) Computation of the look-up table containing the time-averaged values of species mass fractions, density, and temperature as a function of mean mixture fraction, mixture fraction variance, and enthalpy by means of the integrations in Equations (22) and (23).

(3) The solution of the mass, momentum, turbulence and flow mixture fraction equations, Equations (1), (2), (4), (5), (16) and (20).

(4) Calculation of the spatial distribution of temperature, density and individual chemical species mass fractions by interpolating the values in the look-up table.

(5) Scaling of the density field with pressure with respect to the reference value. Dependence of the temperature and mixture composition on pressure is neglected.

(6) From the results of this simulation, the convective heat flux to the wall is evaluated, so that Equations (25)–(27) in the case of standard polymers, or Equations (30)–(32) and (37) in the other case of liquefying fuels, can be solved simultaneously to compute the new distributions of the variables along the grain surface and, accordingly, the mass flux distribution either from Equation (24), or Equations (34) and (35), respectively, along with the volume source terms in Equations (38) and (39) required in the latter case.

Steps from (3) to (6) are then iterated by adjusting the local values of the mentioned quantities until convergence is reached.

#### *3.5. Fuel Port Diameter Update with Time*

To capture the local regression-rate variations in the firing, transient simulations have been attempted. Here, for transient simulation, a series of steady-state computations is intended, which are carried out on di fferent port sizes calculated over several instants in the firing by updating the fluid-solid interface boundary. A quasi-steady approach is employed by neglecting the time-dependent terms in the balance equations and boundary conditions in the time interval between two successive flow-domain updates. This assumption is based on the consideration that the fluid dynamic characteristic time, which can be assumed to be proportional to the di ffusion time in the boundary layer, is shorter than the time required for the thermal profile adjustment in the solid grain due to a change in the port diameter and, thus, in the regression rate.

The displacements of the computational grid nodes are not uniform all over the grain length but, rather, they vary according to the di fferent calculated values of the regression rate. As the regression rate is defined in the direction normal to the fuel surface, due to the local surface inclination, the displacement of a generic point occurs along both the radial and axial directions.

Starting from a given port profile at the *n*-th time-step, defined by the axial and radial coordinates *x<sup>n</sup> i*, *yn i*of the grid nodes (where the subscript *i* indicates the *i*-th node), a CFD simulation is carried out with the models described above in order to compute the fuel regression rate distribution .*r<sup>n</sup>*(*xi*) at the same time step. A forward Euler integration of the local regression rate is then carried out to calculate the displacements after a fixed time-step, Δ*t* = *t<sup>n</sup>*+<sup>1</sup> − *t<sup>n</sup>*, which, for the *i*-th node can be expressed as

$$
\Delta\_i^n = \dot{r}^n(\mathbf{x}\_i)\Delta t\tag{40}
$$

By indicating with θ*ni* the local inclination of the fuel surface with respect to the axial direction in the *i*-th node, which is equal to the angle between the direction normal to the surface and the radial direction (see Figure 3), the coordinates of the node at the time *t<sup>n</sup>*+<sup>1</sup> are calculated as

$$\mathbf{x}\_{i}^{n+1} = \mathbf{x}\_{i}^{n} - \boldsymbol{\Delta}\_{i}^{n} \cos \theta\_{i}^{n} \tag{41}$$

$$y\_i^{n+1} = y\_i^n + \Lambda\_i^n \sin \theta\_i^n \tag{42}$$

which allow reconstructing the new grain port profile.

Once the new distribution is calculated, the fluid domain geometry is consequently modified, the computational mesh is adjusted to the new geometry and the steady-state numerical simulation at the new time-step is performed.

**Figure 3.** The schematic representation of the *i*-th node displacement calculation.

#### **4. Classical Polymeric Fuels: Numerical Results**

In this section, the numerical results obtained with the model described in Section 3.3.1 for the simulation of hybrid rockets burning classical polymeric fuels are presented.

#### *4.1. Experimental Test Cases*

Two motor firings, both performed by burning HDPE fuel grains with gaseous oxygen, are considered here as test cases.

Test HDPE-1 was performed with the subscale 200-N thrust class rocket. The fuel grain was 220 mm long and the initial port diameter was equal to *D0* = 15 mm. The firing duration was set to 12 s. A graphite nozzle with a 9.6-mm throat diameter and an area expansion ratio of 2.99 was employed.

Test HDPE-2 was performed with the 1-kN thrust engine. The grain length was equal to *L* = 570 mm and the initial port diameter was *D0* = 25 mm; a 44 s firing duration was set in this case. The large scale motor employed a copper water-cooled nozzle with a 16-mm throat diameter and an area ratio equal to 2.5. A detailed description of the rocket engine setup and of the experimental facilities can be found in Reference [58].

All the data averaged over the firings presented in this work have been obtained with the loss mass method; the measurement uncertainty quantification followed the procedure in Reference [5]. The average parameters measured in the two firings are summarized in Table 5. The overall mixture ratio, *O*/*F*, is the ratio of the average oxygen mass flow rate to the average fuel mass flow rate in the burn; the latter was calculated as the fuel mass loss divided by the burning time.


**Table 5.** The test operating parameters measured with HDPE fuel grains.

\* The different measurement uncertainties of the variables involving the oxidizer mass flow rate are due to the different system for measuring the mass flow rate.

#### *4.2. Internal Ballistics Steady Simulation: Comparison with Experimental Data*

A single numerical simulation of the Test HDPE-1 has been carried out by considering the time-spatially averaged grain port diameter measured in the burn (note that we call a steady simulation one in which the port diameter is not updated across the firing); the results are discussed in this section.

Figure 4 shows the plot of the calculated temperature contour with the streamlines overlapped on the top half (with respect to the axis of symmetry) and the fuel mass-fraction in the unburned mixture isolines drawn on the bottom half. From this picture, all the main features of the internal flowfield can be unveiled.

**Figure 4.** Test HDPE-1: The temperature contour plot with overlapped streamlines (top half) and mixture-fraction isolines (bottom half).

The combustor inlet flow is dominated by the development of the oxygen jet emerging from the axial nozzle injector (which is clearly distinguishable from the low-temperature region), which spreads almost linearly up to the impingement point on the grain surface. Upstream of the impingement point, in the entrance region of the grain, there is an extended recirculation region characterized by a main, broad counter-clockwise-rotating (if seen from the top half) vortex that is bounded, on the front side, by the zone of oxygen impingement. In the pre-combustion chamber, another large vortex, clockwise rotating, is formed by delimiting the main one on the backside. Finally, also in the aft-mixing chamber, a large trapped counter-clockwise-rotating vortex is formed, which further promotes the propellant mixing.

As a result of the flow recirculation generated at the motor head end, the propellant mixing is strongly promoted, and combustion takes place in the recirculation core; hot combustion gases are transported from the grain entrance region back to the prechamber, where the temperature is very high [5]. Downstream of the recirculation, the temperature distribution reflects the typical structure

of a diffusion flame, with a narrow region close to the fuel surface where the near-stoichiometric conditions are reached, and the temperature shows its maximum value. Anyway, as a consequence of the relatively high turbulent kinetic energy determined by the different vortices, relatively high temperatures also characterize the core flow.

Figure 5 shows a comparison between the computed fuel regression rate axial profile and the measured data. Experimental points are obtained by sectioning the grain transversally in a number of slices, and measuring the port diameter by means of a caliper; in each transversal section, the minimum, maximum and the average of eight diameter measurements have been recorded; the corresponding time-averaged regression rate is calculated by dividing the port radius variation in the firing by the burning time.

**Figure 5.** Test HDPE-1: The comparison between the computed and measured regression rate axial profile.

The regression rate axial trend yields a peak for the oxygen jet impingement, followed by a minimum point, after which it monotonically increases. This latter behavior is typical of the boundary layer heat transfer, for which the heat flux increase due to the mass addition down the port becomes dominant on the decrease due to the boundary layer growth from a certain axial distance.

Provided that a comparison is drawn between the regression rate calculated by simulating the flowfield at the average port diameter and experimental data retrieved from the port shape measured after the motor extinguishment, a good agreemen<sup>t</sup> is shown. The maximum deviation is yielded in the zone of maximum fuel consumption where, however, the experimental uncertainty is maximum because of the asymmetric consumption determined by the motor ignition device [35].

The numerical model is able to predict also the pressure measured in the aft-mixing chamber as well as an indication of the combustion efficiency. The latter is defined in terms of the characteristic exhaust velocity efficiency, η, as the ratio

$$
\eta = \frac{c^\*}{c\_{th}^\*} \tag{43}
$$

where *c\** is the characteristic exhaust velocity estimated with the pressure, total mass flow rate and nozzle throat area, *At*, either measured or numerically calculated

$$\mathcal{L}^\* = \frac{p\_c A\_t}{\left(\dot{m}\_{c\alpha} + \dot{m}\_f\right)}\tag{44}$$

and *c*∗*th* is the theoretical exhaust characteristic velocity obtained in adiabatic chemical equilibrium conditions at the given overall mixture ratio with the CEA code [51]. The fuel mass flow rate is calculated from the numerical simulations by integrating the fuel mass flux along with the fuel port; *O*/*F* is, hence, obtained as the ratio between the prescribed oxidizer mass flow rate entering in the computational domain and the total fuel mass flow rate. The calculated parameters with the relative

Overall mixture ratio,

*c*∗ efficiency

Aft-mixing chamber pressure, bar

*O*/*F*

deviations with respect to the measured values are reported in Table 6 and an excellent agreemen<sup>t</sup> is obtained.


5.54

 6.60

+1.6%

+1.7%

+1.1%


 0.966 \* Experimental regression rate to be intended as time-space averaged value.

#### *4.3. Internal Ballistics Transient Simulation: Comparison with Experimental Data*

It is clear that, although numerically simulating the flowfield in the hybrid rocket combustion chamber at the time-space average port diameter in the burn is an e fficient method to provide meaningful details, as well as the time-space, averaged regression rate, a transient simulation, in the sense specified above, is required to compute the *O*/*F* time shift and the fuel consumption distribution at the end of the burn, especially when the regression rate axial profile is uneven and strongly depends on the grain port geometry, as observed in Section 4.2.

The results of such a simulation carried out in the conditions of Test HDPE-2, by updating the local port diameter in 22-time-steps of 2 s each according to the procedure described in Section 3.4, are analysed.

Figure 6 shows the plots of the calculated temperature contours with the streamlines on the top half and the fuel mass-fraction in the unburned mixture isolines on the bottom half, which have been obtained at several time instants in the burn. First, note that the fluid domain enlarges with the time step because, of course, the port opens up and its shape changes as well for the axial regression-rate non-uniformity

In all the time frames, the oxygen jet spreading and the associated recirculation region in the grain inlet portion can also be clearly seen in this larger-scale engine configuration. The recirculation region at the motor start-up is confined to the prechamber, while, as the grain port enlarges, the oxygen jet impingement point on the fuel wall moves downstream and the recirculation region becomes larger.

The numerical results confirm that the forward axial shift of the impingement point at each time step, once the port diameter is known, can be easily estimated by referring to the characteristics of a free jet ensuing from the injector with a spreading angle of 8 degrees, as observed in Reference [59].

**Figure 6.** Test HDPE-2: The temperature contour plot with overlapped streamlines (top half) and mixture-fraction isolines (bottom half).

The lateral boundaries of the computational domain shown in Figure 6 sugges<sup>t</sup> that the maximum fuel consumption is achieved in the recirculation region, nearly upstream of the impingement point. Accordingly, the axial profiles of both the regression rate and local port diameter, displayed in Figure 7, yield maximum points, which move downstream in the burn, for the jet impingement shift described above.

**Figure 7.** Test HDPE-2: (**a**) the regression-rate and (**b**) the fuel-grain port diameter axial profiles at several time steps in the burn.

Moreover, because of the mass flux decrease, a reduction of the average regression rate can be observed.

The average parameters computed with the numerical simulation are summarized in Table 7 along with the deviation with respect to the corresponding experimental data. A good agreemen<sup>t</sup> between the numerical results and experimental data has been obtained, which validates the employed model. The time-space averaged regression rate is fairly predicted with a deviation with respect to the measured value of only 1.2%; note that, based on the definition of the average regression rate, the latter coincides with the error on the postburn space-averaged port diameter prediction.

The numerically predicted average values of the mixture ratio and of the chamber pressure are calculated as the ratio of the average oxidizer mass flow rate to the average fuel mass flow rate and as the arithmetic mean of the chamber pressure over all the time steps, respectively. The numerically predicted average combustion efficiency is calculated by means of Equation (44), considering the average values of the chamber pressure and mass flow rates.

The prediction of the aft-mixing chamber pressure is affected by the maximum deviation of 5%, which can be explained through the combined effect of the errors in the fuel mass flow rate (and, consequently, in the overall mixture ratio) and combustion efficiency that both affect the chamber pressure. Once the comparison is made in terms of the combustion efficiency, a better agreemen<sup>t</sup> is obtained, with a deviation of 3.7%.


**Table 7.** Test HDPE-2: The computed average parameters and deviation with respect to the experimental data.

Test HDPE-2 allows for an assessment of the calculated local time-resolved regression rate. The ultrasound pulse-echo technique was used to measure the local fuel web thickness and derive the regression rate. One ultrasonic transducer was located around the chamber mid-span [60], whose axial location along the grain is indicated in Figure 7.

Figure 8a shows the time trends of the local regression rate and port diameter. The numerically computed regression rate fairly reproduces the characteristic behavior shown by the measured data for which, around 25 s, there is a minimum consumption rate. This feature can be explained considering that the regression rate initially decreases for the mass flux reduction, then, owing to the shift of the

oxygen jet impingement point towards the probe location, the recirculation region enlarges and the consumption rate increases (see also Figure 7a).

Additionally, despite the bias relative to the measured trend, the increase of the aft-chamber pressure over the firing, which is due to the increase of the fuel mass flow rate, is well predicted (Figure 8b). The nearly constant difference between the numerical and experimental curves is essentially a consequence of the regression rate underestimation (see the deviation in *O*/*F* shown in Table 7), which remains constant at the different time steps.

**Figure 8.** Test HDPE-2: The comparison between the computed and measured (**a**) local regression rates and port diameters; (**b**) the pressure and mass flow rates

#### **5. Para**ffi**n-Based Fuels: Numerical Results**

The numerical model presented in Section 3.3.2 has been applied to the simulation of several test cases obtained from the static firing of the sub-scale hybrid rocket burning paraffin-based fuels; the results are discussed in this section.

#### *5.1. Experimental Test Cases*

Seven test cases are presented here, which were performed by varying the oxidizer mass flow rate and the firing time, with the aim of achieving a significant range of the average oxidizer mass flux, which is the main control variable in the hybrid rocket's operation (see Table 8). Paraffin-based fuel grains were made of a blending of a low-melting-point paraffin wax and a microcrystalline wax; gaseous oxygen was axially injected in the grain's single port. In all the firings, an exhaust nozzle having the throat diameter equal to 10.7 mm and an area expansion ratio equal to 2.4 was employed. The details of the test campaign from which the experimental data have been gathered can be retrieved in Reference [31].

Steady numerical simulations are carried out on the port geometry corresponding to the average diameter in the burn and the imposing oxygen mass flow rate (both reported in Table 8).


**Table 8.** The test operating parameters measured with paraffin-based fuel grains.

#### *5.2. Role of the Regression-Rate Entrainment Component*

The primary difference of the paraffin-wax consumption mechanism with respect to the classical polymer pyrolysis relies upon the liquid layer formation and entrainment in the gas stream. With the aim of highlighting the role played by the entrainment component of regression rate on the flowfield, two numerical test cases built with the input conditions of Test P4 (Table 8) have been considered.

First, an extreme case is analyzed in which the entrainment component is assumed to be zero, so that the overall regression rate, as with a standard polymer, is only due to vaporization, i.e., .*r* = .*rv*, (Equation (32) is, thus, not considered), and it is compared to a second case in which the entrainment term also is included in the calculations. In the second case, the entrainment parameter has been assumed to be equal to 2.1 × 10−<sup>13</sup> m8.5s0.5/kg3; the latter is the one allowing for the best fit of the experimental data obtained in test P4 as discussed in Reference [56]. In addition, the calculations are repeated with an HDPE fuel grain of equal dimensions by imposing an equal oxygen mass flow rate, i.e., with the model described in Section 3.3.1.

Figure 9a shows the obtained fuel regression rates and Figure 9b, the corresponding surface heat fluxes. Total regression rate and relevant heat flux have similar axial profiles. It is worth noting that, as the heat requested for HDPE pyrolysis is larger than that used in Equation (31) for modeling paraffin melting and vaporization (about 5500 kJ/kg against 1400 kJ/kg, respectively), in the first extreme case, without considering the entrainment contribution, the regression rate (see the dark grey continuous line in Figure 9a) is fairly higher than that obtained with HDPE (light grey continuous line in Figure 9a) at equal oxidizer mass fluxes, despite the fact that the enhanced blocking effect determines a significantly lower surface heat flux (see Figure 9b). However, the spatially-averaged regression rate obtained in this case is equal to 1.11 mm/s, which is still significantly lower than the corresponding measured value of 2.29 mm/s (see Table 8).

When the entrainment component is taken into account, the calculated regression rate is more than doubled (see the black continuous line in Figure 9a) because, with the set of parameters considered here, the most significant contribution is given by the entrainment itself (black dotted line), the vaporization component (black dashed line) being smaller than the entrainment fraction. With the larger mass flux due to the entrained fuel, as entrainment does not contribute to the heat-transfer blocking, the heat flux is raised (see the black line in Figure 9b). In particular, in the fore end of the grain (up to about 80 mm), where the effect of the mass addition is low, the vaporization regression rate yields values similar to the HDPE regression rate; accordingly, comparable equilibrium conditions between the heat transfer and the mass blowing at the grain wall is obtained, and the heat flux profiles are similar, whereas, downstream of that point, for the largely increased mass flux due to the entrainment component, the heat flux to the paraffin fuel surface significantly diverges from that achieved with HDPE.

**Figure 9.** The comparison between the results obtained with and without considering the entrainment.

Figure 10 shows the two contour maps of the temperatures calculated with paraffin (considering both vaporization and entrainment) and HDPE. It is worth noting that, in the case of paraffin fuel burning, the recirculation region is smaller for the modeling of entrained paraffin mass that, as mentioned above, is introduced in the port volume; the latter needs to be axially accelerated at the expense of the oxygen jet momentum, which decreases and causes a larger jet spreading. For the same reason, the hottest region in the flowfield of the HDPE grain port is attained close to the grain surface, whereas, in the paraffin-fuel port, it rapidly extends into the core flow because of the significant fuel mass addition largely due to the entrainment.

**Figure 10.** The temperature contour plots with overlapped streamlines (top half) and mixture-fraction isolines (bottom half), (**a**) paraffin fuel, (**b**) HDPE fuel.

#### *5.3. Comparison between Numerical and Experimental Results*

The firing tests presented in Section 5.1 have been simulated for the sake of model validation. The entrainment parameter used in each test case has been obtained by scaling the reference value of 2.1 × 10−<sup>13</sup> m8.5s0.5/kg3, identified for the best fit of the experimental data in Test P4, with the ratio ρ<sup>∗</sup>*g*/ρ*<sup>g</sup>*1.5, where ρ∗*g* is the average gas density in the grain port calculated in Test P4 and ρ*g* is the corresponding value calculated in the analyzed test case. This allows for considering the dependence of the entrainment parameter on the average gas density as prescribed by Equation (33).

Figure 11a shows the calculated fuel regression rates averaged along the grain compared with the measured time-space averaged ones; the maximum deviation of 11% is reached at the minimum mass flux (Test P1 in Table 8). The numerical prediction improves with higher mass fluxes, showing excellent agreemen<sup>t</sup> at the largest mass fluxes where the deviation is only 0.3% (still lower than that accepted at the reference, Test P4).

(a) Regression rate (b) Postchamber pressure 

**Figure 11.** The comparison between computed and measured (**a**) regression rate, (**b**) pressure.

In Figure 11b, the calculated chamber pressures in the whole set of test cases are compared with the measured values; data are retrieved in the rocket aft-mixing chamber. All the calculated pressures are overestimated by at most 10% except those of Test P1 and Test P2, whose deviations from the measured data fall around 20%. A detailed analysis of the factors of deviation between the computed and measured pressure is addressed in Reference [56]. However, the displayed deviation trend can be explained by observing that the critical pressure of paraffin wax is 6.5 bar [42] and that the chamber pressure attained in the test with the largest deviation is lower than the critical pressure (see Table 8). The agreemen<sup>t</sup> with experiments is improved as the pressure increases with the mass flux. In fact, below the critical pressure, neglecting the effects of the entrained liquid paraffin dynamics is a much less suitable assumption.

## **6. Conclusions**

A review of the computational thermo-fluid-dynamic models developed by the authors for the simulation of the internal ballistics of a hybrid rocket burning either classical polymeric or paraffin-based fuels with gaseous oxygen is presented along with some basic results compared to the experimental data.

The theoretical approach consists in solving the single-phase RANS equations with two additional transport equations for the average mixture fraction and its variance combined to the PDF combustion model and thermochemical equilibrium. In addition, for the prediction of the fuel regression rate, two integrated sub-models that are able to describe the interaction between the gaseous flow and the grain surface depending on the nature of the fuel itself—non-liquefying or liquefying—were implemented. The gas/surface interaction modeling is based on the local mass, energy and mixture-fraction balances, regardless of the type of fuel, but specialized treatments were needed to cope with the di fferent consumption mechanisms of classical polymers and liquefying fuels. In the former case, a pyrolysis rate equation correlating with the regression rate to the surface temperature was included for predicting the fuel burning rate; whereas, in the latter case, an equation relating to the additional regression rate fraction (which originates from the liquid entrainment phenomenon) with the total mass flux and the total regression rate through an entrainment parameter, was implemented.

A procedure for the grain port local section evolution in the burning was also defined, accounting for the regression rate axial unevenness.

A series of test cases gathered from the static firings of two di fferent laboratory-scales hybrid rockets that were burning gaseous oxygen and HDPE or para ffin-based fuels were numerically reproduced in order to assess the numerical model capability of predicting the regression rate and chamber pressure. With this purpose, on the assumption that combustion is limited by turbulent di ffusion (and, thus, the combustion e fficiency is governed by the mixing of propellants), the flow recirculation occurring upstream and downstream of the fuel grain needs to be simulated; accordingly, the computational domains were built on the exact thrust chamber geometry including the real preand post-chamber configurations and the exhaust nozzle.

Two tests with HDPE were analysed: one on the small-scale engine in which a steady simulation was carried out to compare the regression rate axial profile, and another with a larger-scale rocket in which an unsteady simulation allowed for making a comparison with the local time-resolved regression rate and chamber pressure that was experimentally measured. The model has been demonstrated to capture (not only from a qualitative standpoint) the fundamental features of the motor internal ballistics, which involve flow recirculation at the grain port inlet with the characteristic heat transfer mechanism, leading to a point of maximum regression rate that shifts downstream during combustion. The calculated regression-rate time trend, despite the relatively large time step of the port diameter integration, fairly reproduces the change of the wall heat transfer consequent on the recirculation region dynamics.

Seven firing tests with para ffin-based fuels were reproduced with steady simulations to evaluate the approach formulated to model the entrainment fraction of the regression rate, which is based on two main considerations: first, at a combustion pressure larger than para ffin critical pressure, one may disregard the fact that a significant fuel fraction is injected in the main flow in the liquid state by retaining the key aspect that the entrained fuel does not contribute to the heat transfer blocking; second, the entrainment parameter given the fuel formulation only depends on the mean gas density in the port.

In all the analyzed cases, a satisfactory agreemen<sup>t</sup> between the calculated regression rates and chamber pressure with the relevant measured data is obtained.

In particular, the chamber pressure is predicted with lower accuracy; a more accurate prediction likely requires modeling the two-phase flow dynamics ensuing from the melted para ffin injection into the chamber, which is the subject of future developments.

**Author Contributions:** Conceptualization and definition of the CFD models, G.D.D.M., C.C. and R.S.; implementation of numerical models, G.D.D.M.; running of numerical simulations, G.D.D.M. and S.M.; carrying out of experimental tests, G.D.D.M. and S.M.; writing—original draft preparation, C.C. and G.D.D.M.; supervision, C.C. and R.S.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
