Swirled injector modeling for cavitating multiphase flow

The goal of this study is to present a comparison between different approaches to multiphase injection modeling of self-pressurized rocket engine propellant. Swirled, tangential orifice injector of nitrous oxide, for an “N” class hybrid rocket motor is the object of the study. A brief descriptionof the injector purpose and geometry is provided, followed by a description of different approaches for flow modeling. Examined techniques range from 0D, Homogeneous Equilibrium Model (HEM) to 3D multiphase with mass and heat exchange between phases. Results of analyses are provided and compared with experimental data. The discrepancies between results are of significant magnitude but expected nature. Co clusions about most feasible approaches for engineering calculations are drawn.


Introduction
The object of this study, is an open type swirled injector, being designed for an "N" class hybrid rocket engine (HRE) using nitrous oxide as an oxidizer. This kind of injector and it's more complex versions are commonly used in rocket engines, especially bipropellant engines requiring e↵ective mixing and atomization [1]. It's flowpath comprises of six tangentially oriented, o↵set from injector main axis circular holes referred to as injector orifices of 2.5 mm diameter, placed at 5.5 mm radius (o↵set between orifice axis and injector axis). Cavitating fluid is injected into the swirl chamber tangentially. Introduced angular momentum causes the fluid to follow cylindrical swirl chamber walls and build up axial velocity. Fluid leaving the swirl chamber has therefore both angular and axial momentum and forms a spray cone.
The rationale for using nitrous oxide as an oxidizer for small hybrid rockets is its capability of being stored in a liquid state at ambient temperatures (critical temperature of 304.9 K) and self-pressurizing capability (it's saturation pressure is 50.2 bar at 20 • C). Material properties of nitrous oxide pose a significant difficulty in numerical modeling of flow in question. In self-pressurizing oxidizer tank, the fluid is stored in the saturated state and starts to cavitate, once the fluid static pressure decreases along the ⇤ e-mail: bartosz.ziegler@put.poznan.pl flowpath. Additionally, the fluid parameters of state are very close to the critical point where liquid phase has a significant compressibility and the gas phase state cannot be modelled accurately with ideal gas law (or even most of the real gas laws). Therefore, in the design process of nitrous oxide injectors, only experiments provide reliable results. Nitrous oxide is easily substituted with carbon dioxide having nearly identical thermophysical properties as N 2 O, while being more cost efficient and safer in handling. In this paper three approaches for injector flow modeling, based on real fluid saturation line are compared and discussed in relation to measured data. All modeling approaches, although intended for nitrous oxide are compared on carbon dioxide flow case as also experiments were performed with use of CO 2 as a substitute for N 2 O. All computations, CFD analyses and the experiment were performed with the use of carbon-dioxide. Thermophysical data of CO 2 along saturation line were taken from NIST webbook [2] and all values of interest were approximated with high accuracy as functions of reduced temperature (ratio of temperature to critical temperature). The functions were fitted to the same form as highly accurate analytical forms found for nitrous oxide [3]. Such approach allows to use the same frameworks for both media, changing only the coefficient tables.

HEM based model of swirled injector
Basic assumptions of HEM is that the interphase mass and heat transfer rates are high enough to keep two phases in thermodynamic equilibrium, ensuring that both phases are always in saturated state and that they are mixed with each other to an extent that the slip velocity is negligible. Those assumptions make both phases share static pressure and temperature as well as velocity. HEM based models for 1D and multidimensional analyses [4,5] has been developed in the past and used for modeling cavitating flows, liquid depressurization etc. and as a comparative model for nonequilibrium methods. Most of those models are inherently time-transient and were developed for water hence nonequilibrium models, for engineering calculations nitrousoxide flashing flows are poorly available. In this work, greatly simplified approach of zero-dimensional injector model, with HEM based expansion work integration along saturation line between characteristic points in the injector flow field is applied as its low computational cost makes it feasible for engineering calculations. For an adiabatic flow, according to the first law of thermodynamics: , for an expansion process can be also rewritten as: With η being the expansion efficiency. For any given initial state (T, p, X, u), the expansion work (defining the change in kinetic energy) can be integrated numerically along the saturation line, where pressure is a known function of temperature, as are also enthalpies and densities of both phases. The bulk mixture enthalpy and density can be calculated for a given gas mass fraction: Using the bulk mixture density (equation 5), integration of mixture enthalpy change (equation 2) for assumed efficiency can be performed along the saturation line. For known saturation h gas and h liq equilibrium vapour mass fraction can be determined from equation 3. The HEM based model of swirled injector, is based on described expansion work integration between characteristic points (Fig. 1) of the injector. Between stations 1 and 2, an isentropic expansion is assumed, which for a given hole area and discharge coefficient Cd defines the mass flow rate and angular momentum of the fluid entering the swirl chamber. After station 2, all the walls are axisymmetric, therefore the angular momentum is assumed constant. Another expansion happens between stations 2 and 3, where axial momentum is gained by the fluid and angular momentum is assumed constant. Between points 3 and 4, radial equilibrium is assumed, allowing integration of the liquid film thickness for given temperature range as: Generalized reduced gradient (GRG) optimizer [6], is applied to the model to find temperatures (for all characteristic points) for which mass flows are matched and the pressure in point 4 (see Fig. 1) of the model is matched to the given backpressure.

CFD approach
The flow in question is of highly transient nature. Because the pressure drop in the injector orifice happens very quickly (on a distance of around 1 mm), and the flow has in that place velocities of the order of 35 m/s, roughly estimated timestep of 7e-7 seconds would be required for CFL = 1 at injector orifice inlet. Transient simulations of other researchers suggest use of even smaller time steps of 1e-7 s [7]. In mentioned research, it was required to perform 300 000 time steps for two orifice flow-throughs, what required 1500 core-hours on a 2D axisymmetric 20 000 element mesh. In this study meshes two orders of magnitude bigger, with longer required flowthrough (because of longer fluid path in the injector) times are used what makes transient approach impractically costly. For design purposes the main result of interest is the quasi steady mass flow rate through the injector so transient results are not required.
Therefore the following approach is taken. The flow is computed with a steady-state analysis with inlet and outlet pressure boundary condition representing a given timestamp of the experimental flow. Pressure values for inlet and outlet are obtained from the experiment. For described study, time of 4, 5.5 and 7 seconds are chosen for analysis. Two numerical setups are considered for the CFD analyses: a multiphase Eulerian analysis with heat and mass transfer between phases and an incompressible analysis.
As the orifice Reynolds number can be as high as 700 000, in all CFD cases Menter k-ω SST [8]. turbulence model was used. This turbulence model was chosen because of its known versatility and accuracy in predicting separated regions and shear layers [9]. Calculations are performed in ANSYS Fluent.

Mesh
The computational grid used for the study is a fully hexagonal grid of a 60 degree slice of the injector pre-chamber, swirl chamber and pressure chamber with one tangential orifice channel. The grid is created as a multiblockstructured, and converted to an unstructured mesh before calculations. The geometrical dimensions of the domain resemble the injector and pressure chamber of the stationary hybrid rocket motor used by PUTRocketlab Students Research group for evaluation of fuel mixtures. The mesh was created in Ansys ICEM meshing software. Fully structural mesh was obtained by using a yblock topology for triangular geometries combined with an o-grid topology for swirl chamber and pressure chamber walls. The grid is made with conformal periodic interfaces. Although the flow is mostly axisymmetric and swirled, during preliminary calculations, it was found that utilization of 2D grids with axisymmetric swirl momentum formulation results in highly unstable analyses. Results obtained from a 2D cases were also regarded as nonconclusive as the most important feature of the flow field determining injector mass flow rate, is the tangential injection port, which cannot be properly modeled in 2D. Wall spacings in injector orifice and prechamber were adjusted for a non-dimensional first element spacing (y+) parameter of 30, calculated for liquid CO 2 and for an average velocity expected in the injector orifice as this region is expected to be inertially dominated and orifice flow separation (being a feature mostly determining the orifice discharge coefficient) should be easily reproduced with lower viscous sublayer mesh resolution. In the swirl chamber and downstream of the injector, y+ was set to 1, to ensure possibly accurate solution of the boundary layer parameters in the swirl chamber. This was considered important because the flow travels much longer distance along that wall. Total element count in the final grid is 1 929 390 elements with maximum non-orthogonality of 71.6 degrees and worst skewness criterion of 0.2. The mesh quality metrics are directly imposed by sharp angles between orifice and swirl chamber walls resulting in sharp internal element angles. Despite that, no problems with computation stability were encountered.

Multiphase model
A two phase, Eulerian multiphase model with Lee [10] evaporation-condensation model, and fixed to saturation temperature heat transfer model (which assumes that the vapor phase is at saturation temperature, and that the heat transport is driven solely by mass transport between phases) was used in all multiphase calculations. Real fluid thermophysical properties for both phases are strongly dependent on pressure and temperature, for this calculations however, it is of highest importance, not to predict material properties in the flow field, but to mimic as closely as possible the mass flux and phase composition along the expansion process. Therefore, both phases thermophysical parameters (liquid density, standard temperature difference in specific enthalpies, individual gas constant and specific heat capacities) were determined with the generalized reduced gradient method to minimize the relative square deviation between mass fluxes for an isentropic expansion process with phase change for constant properties and for real CO 2 values.
Comparison of mass flux and vapour mass fraction with thermodynamic equilibrium assumed between phases (HEM model) for the real fluid data and devised model are depicted in Fig. 5 for one of the fluid initial states from the experiment (timestamp t = 5.5 s). Mass flux for devised model is very close to the one calculated for real fluid data (above 30 bar, relative RMS error is below 10-3). The model loses some of mass flux accuracy below 30 bar pressure, but analyses performed in this study are not expected to reach that pressure region.

Incompressible model
For comparison with other results also incompressible liquid flows were calculated on the same mesh and boundary conditions. For each analysis inlet boundary condition, fluid density and viscosity were taken from real fluid data. This incompressible analysis presents the opposite side of the density changes spectrum compared to HEM analyses, as it assumes zero rate of vapour production while HEM with its equilibrium condition assumes infinitely e↵ective heat and mass transfer bringing both phases to thermodynamic equilibrium with no delay.

Experimental setup
A set of experiments were conducted with CO 2 as the injected fluid, with similar to HRE physical setup (Fig. 6). The chamber had no fuel grain, and the engine nozzle was replaced with much lower area orifice to pressurize the chamber in a similar manner to its intended operation (35 bar chamber pressure). Carbon dioxide tank is suspended on a load cell and connected via bending connection to the injector screwed into a chamber which is closed with a 5.5 mm sharp edged orifice nozzle to pressurize it. p1 and p2 pressure transducers measure the tank and chamber pressures respectively. The flow is initiated by opening an electrically operated ball valve. Flow starting causes substantial inertial forces to be exerted on the tank, what causes mechanical vibrations in the system (weight fluctuations on Fig 7). Additionally cavitating flow even in a quasi-steady state introduces more transient forcing. Resulting vibrations introduce significant noise into the CO 2 weighing system which has to be filtered out, therefore to estimate the mass flow rate, raw load cell data are fitted to a smooth curve which analytic derivative is used as a mass flow rate indicator. From figure 7, it can also be seen that the first 3 seconds of the test (roughly two first seconds of the flow) experience also some peculiar behaviour of the tank pressure. After valve opening pressure drop resembling pressure drops for sole gas discharges occurs, followed by a sudden spike in pressure around 2.5 s of the test. It is suspected that this phenomenon might be a result of big vapor cavities creation and collapse in the fluid before the feed line is filled, or even with local solid phase formation. This phenomenon was observed in most of the tests, and remains unexplained. Because of that only time range from 4 to 8 seconds is used for curve fit and mass flow rate estimation and all timestamps used for comparisons are from this period.

Results and discussion
Main result of interest for each analysis technique is the injector mass flow rate, as this is the parameter most directly influencing the HRE operational parameters. Mass flow rates for performed analyses are depicted on Figure 8 and compared with the mass flow rate calculated as an analytical derivative of the carbon dioxide tank weight fitting curve.
Incompressible CFD results show highest deviation from the experimental mass flow rate (15.2% on average). This overprediction is due the fact that neglecting cavitation, results in high overestimation of the bulk fluid density in the injector orifice. The 3D incompressible analysis, was used to determine the injector orifice discharge coefficient equal in average 0.6715 (di↵erences between Cd for all three cases where of the order of 0.001). Such small discharge coefficient, results from extensive flow separation on one side of the injector orifice, which can be seen on mid orifice velocity fields (Figures 9 through 11) for both incompressible and multiphase cases as it is an inherent feature of given geometry. Eulerian multiphase results, even though obtained with Lee model evaporation coefficient an order of magnitude higher than mentioned in Fluent User's Guide [11] coefficient range upper bound, overpredict mass flow rate significantly (6.1% on average). HEM model on the other hand, if discharge coefficient from the incompressible CFD is applied to injector orifice, gives results with significant mass flow rate underestimation (-14.4% on average) even though the orifice flow in the HEM model was assumed isentropic. This is mostly, due to the equilibrium assumption. In real flow, the mass transfer between phases happens with a finite rate, limited by bubble dynamics and heat transfer within the fluid. The orifice of presented injector introduces very high pressure drop on a very short distance (at the centerline of the orifice inlet, pressure gradient reaches values above 360 MPa/m, with flow velocities in that region above 25 m/s) what results in huge pressure drop rate above 9 GPa/s for the fluid in that region. This value even left without context, shows why on one side, the HEM assumption of equilibrium between phases is inaccurate, and also why for a multiphase calculation impractically high evaporation coefficients are required for which probably the only way to determine them is with experimental results for a physically similar case. The same can be done with the simplest (HEM) approach by simply increasing the orifice discharge coefficient above its typical (geometric) value, based on experimental results, what is shown on Fig. 8 with HEM results for injector orifice discharge coefficient of 0.8. A way for predicting mass flow rates based on HEM procedure, was also proposed by J. Dyer et. al. in [12]. This method, developed for high aspect ratio orifices (high length/diameter) basically calculates the mass flow rate as a weighted average of HEM and incompressible predictions for given geometric discharge coefficient. Results calculated with Dyer et. al. equation, are also shown on Fig. 8 for discharge coefficient Cd = 0.6715 (deter-  Table 1. Flow patterns in the injector orifice midplane are shown in figures 9 through 11, they are of the same nature (separation region position and size, velocity distribution in the orifice) but in general incompressible flow analysis underestimates the flow velocity in comparison to the multiphase analysis (which is an obvious result of higher bulk fluid density). Moreover, a discrepancy between velocity of the phases for multiphase analysis (non-zero slip velocity). Figure 12 presents liquid phase volume fraction at the injector orifice. Liquid fraction increases almost evenly along the orifice while almost whole orifice pressure drop happens near its inlet (Fig. 13). This is an evidence of deviation of mixture properties from homogeneous equilibrium assumptions. Therefore it can be generalized, that in the orifice: This relation makes the non-equilibrium multiphase models over predicting mass flux of the mixture compared to models with homogeneous equilibrium assumption while at the same time underpredicting the flux compared to incompressible model. The value of flux lies between that for HEM and incompressible models, and where between those, is solely dependent for given case, on phase change model coefficients. Therefore the multiphase model, even though qualitatively represents flow physics in highest detail, is as accurate in predicting flashing orifice mass flow rate, as is its choice of coefficient value, which should be fine-tuned with experimental data. In all cases similar spray cone angles where achieved (Fig.  14) with typical for swirled injector swirl chamber with liquid phase rich mixture following swirl chamber walls, and almost completely gas filled core of the swirl chamber.

Conclusions
Performed analyses suggest, that mass flow rate prediction for a cavitating, tangential orifice swirled injector can be a challenging task even with commercial CFD tools which usually allow for very accurate mass flow prediction for single phase gas or liquid applications. This is due to the fact that each approach taking phase change into account, requires a case specific coefficient to address the finite rate of mass and heat transfer between phases and resulting deviation from thermodynamic equilibrium between phases in the injector orifice. On one side of the spectrum is the incompressible analysis which significantly overestimates the mass flow rate by neglecting the phase change, on the other side is the HEM approach which assumes thermodynamic equilibrium at saturation parameters and no slip velocity between phases. The relative di↵erence between those two limiting methods exceeds 30% for presented cases. All presented approaches expect the incompressible single phase analysis seem feasible. But they require a case specific coefficient to predict the mass flow rate (Dyer et. al. method presents a way to estimate its coefficient analytically). 3D CFD analysis gives more information about the flow spatial distribution, but for sole mass flow rate vs. boundary conditions characteristics 0D HEM with experimentally established injector orifice e↵ective discharge coefficient or the Dyer et. al. [12] method, seem more feasible because of negligible computational e↵ort compared to 3D CFD analysis. Dyer et. al. [12] formula, appears to be most feasible as this approach requires only the "geometric" discharge coefficient which for many practical geometries can be found in literature.