Thermodynamic Effect on a Cavitating Inducer in Liquid Hydrogen

This study was led in collaboration with the French space agency (CNES) and the rocket engine division of Snecma. The main aims were the simulations and the analyses of cavitating ﬂows in the rocket engine turbopump inducers, where the operating ﬂuids are LH 2 and LOx under cryogenic conditions. A ρ ( P, T ) state law modeling the cavitation phenomenon was integrated by the labora-tory LEGI in the commercial CFD code Fine/Turbo TM , developed by Numeca International. Various 3D numerical results are given for an inducer geometry and comparisons are made with experimental data (head drop curves) obtained by NASA.


Introduction
The simulation and the prediction of cavitation in cryogenic fluids is of critical importance for the efficient design and performance of turbopumps in rocket propulsion systems. In temperature sensitive fluids as LH2 and LOx, thermal effects and strong variations in fluid properties are observed, which alter the cavity characteristics.
For cryogenic fluids, the liquid-vapor density ratio is lower than in cold water and consequently more liquid mass has to vaporize to sustain a cavity. Therefore evaporative cooling effects are more pronounced and the temperature of the liquid in the immediate vicinity of the liquid-vapor interface is depressed below the free-stream temperature. Because of the strong variation of thermodynamic properties (mainly the vapor pressure), the temperature depression, negligible in water, is quite substantial. The local cooling effect delays the cavitation phenomenon and reduces the local vapor pressure of the fluid, which leads to a lower observed cavity pressure. Typically this results in improved mean performance of cryogenic pumps.
First studies about thermal effects were generally focussed on obtaining correlations for temperature depression as a function of flow conditions and liquid properties. Particular methods of importance include the B-factor theory [11,16,21] to characterize the sensitivity of fluids to thermodynamic effects, and the entrainment theory [1,10]. The B-factor method is based on the ratio of vapor volume to liquid volume affected by the vaporization process. A simple heat balance between the two phases can estimate the scale of temperature difference ∆T * caused by thermal effects. The B-factor is estimated as the ratio between the actual temperature drop ∆T and ∆T * : where L vap is the latent heat, ρ L and ρ V densities for the liquid and vapor phases, respectively. C p L represents the specific heat.
Some numerical models have been developed to investigate thermodynamic effects in cavitating flows with one-fluid method. The one-fluid method treats the cavitating flows as a mixture of two fluids behaving as one. The governing equations are composed by three conservation laws written for the mixture. These models are based on the assumption of local kinematic equilibrium between phases (the local velocity is the same for both phases) and local thermodynamic equilibrium between the two components (local temperature and pressure equality between phases). This Homogeneous Equilibrium Model (HEM) cannot reproduce strong thermodynamic or kinetic non-equilibrium effects but, due to its simplicity, it has often been used for the numerical simulation. Different equations of state (EOS) have been used to define the thermodynamic behavior of the mixture to simulate cavitation: in liquid hydrogen [2], in warm water [18], in freon R-114 [6,19] and in octane [4].
A four-equation model, very popular to simulate cavitating flows in water, has been used in cryogenic applications [12,13,22,25]. It is obtained by adding a mass equation for the 3 vapor or liquid density including a cavitation source term. The main difficulty is related to the formulation of the source term and the tunable parameters involved for the vaporization and condensation process (different sets of parameters are presented in [22]).
This paper proposes an equation of state based on a modified barotropic law to model the mass and heat transfer exchanges. In the following, we first summarize the essential elements of the governing equations, the modeling concepts and the numerical methods, then present the three-dimensional results on an inducer.

Models and numerics 2.1 Equations and models
The solver is based on the compressible Reynolds-averaged Navier-Stokes equations, expressed in a relative frame of reference. Cavitation is modeled through a single fluid approach with a ρ(P, T ) equation of state. In the cavitation region, the liquid-vapor mixture is described with its averaged density, velocity and pressure. The density is assumed to be constant in the pure phases, while in the mixture region, it is directly linked to the pressure and the temperature with the EOS. The temperature T is computed with the total energy equation and constant thermal capacities. The EOS is based on a barotropic law, proposed by Delannoy and Kueny for cold water [3], modified by introducing thermal dependency of vapour pressure [20]:

Numerical methods
The cavitation model has been implemented in the Fine/Turbo T M solver by Rolland [20].
This code is based on a cell-centered finite-volume discretization for multi-domain structured meshes [5].
For the mean flow, the convective flux density is computed with the space-centered Jameson scheme stabilized with scalar artificial dissipation [14]. The diffusive flux density is computed with a second-order space-centered scheme.
A well-known problem concerns the stiffness on the solution convergence in incompressible areas. In this situation, the dominance of convection terms renders the system stiff and compressible solvers converge slowly. To overcome this difficulty, a preconditioned method is necessary. The physical acoustic waves are replaced by pseudo-acoustic modes that are much closer to the advective velocity, reducing the stiffness and enhancing the convergence.
The preconditioned method developped by Hakimi [8] is used.
The code marches the solution in time to steady state, using an explicit 4 step Runge-Kutta algorithm. The convergence is enhanced through local time stepping, implicit residual smoothing and multigrid acceleration techniques (for non cavitating computations).

Turbulence model
The Yang-Shih k − ε model [24] with two-layer wall functions [9] is used.

Computation procedure to simulate head drop
The cavitating behavior of the inducer is investigated through the technique described by Pouffary et al. [17] for cold water. The pump flow is first computed with incompressible liquid hydrogen. Then cavitation is progressively introduced by increasing the vapor pressure.

Geometry and experimental conditions
The tested geometry is a NASA three-bladed flat plate helical inducer ( Figure 1  An enlargement of the blade leading edge is shown in Figure 2, where the bevel is drawn. For this geometry, cavitation can appear along the bevel, on the leading edge and also on the trailing edge. Cavitation appears also at the back end of the bevel (the point B), which corresponds to the attachment point of sheets.

Properties of liquid hydrogen
The LH2 is characterized by large variations of thermodynamics properties with the temperature (especially the vaporization pressure), as shown in Table 1. A change of temperature corresponds to a change in vapor pressure of a factor of three.
Moreover, it increases also the ∆T * value. As illustrated by the experimental results presented in Figure 4, the temperature has a major influence on the NPSH breakdown value.

Mesh and boundary conditions
The grid is a HOH-type topology. Only one blade-to-blade channel is considered. The mesh contains around 740000 nodes. The tip leakage is considered. The y + values at the walls vary between 0.5 and 100. A view of the complete mesh is given in Figure 3.
Calculations are performed with the following boundary conditions. The mass flow is imposed at the inlet, and the static pressure is imposed at the outlet. The hub and the blades are rotating, the shroud is immobile. Wall functions are imposed along solid boundaries, assumed to be adiabatic. 7

Results and analyses
For this geometry, according to non cavitating numerical results, the flow rate corresponding to the best hydraulic efficiency of the pump is about φ = 0.1. The two studied flow rates, φ = 0.098 and φ = 0.108, are respectively lower and higher. For the second flow rate, the angle of attack of fluid is negative, which induces also the development of cavitation in the pressure side of the blade.  Table 1).

Head drop curve
In this table, we introduced a characteristic difference of pressure ∆P * defined as: Both characteristic differences ∆T * and ∆P * are higher for LH 2 at T ref =23K. During the cavitation process, the local cooling down and P vap variations are more relevant at this temperature. These thermodynamic effects delay the cavitation phenomenon.
To deep the analyses, it is interesting to evaluate also the dimensionless pressure gradient Table 2 For the other operating point (φ=0.108, T ref =23K), the breakdown is more abrupt, with small variations in the instabilities zone. The numerical performance breakdown is not so close to the experimental data, the error is around 22%. This inaccuracy is the same order of one obtained in [7] for a cold water inducer.
The physical analyses proposed here below will be based on numerical results obtained for T ref =23K, because the thermodynamic effects is more important than the one at T ref =19K.

Visualization of the vapor/liquid structures
The computed mixture vapor/liquid structures are presented in Figure 5, during the NPSH decrease, for the operating point (φ=0.108, T ref =23K).
For NPSH=108m (top of the figure), a small cavity sheet appears on the leading edge bevel, near the shroud, and also an attached cavity from the point B (see Figure 2). When the NPSH value decreases, the small cavity on the bevel reduces whereas the principal cavity grows and detached from the suction side. Finally, for NPSH=30m, a large cavity reaches the throat of the channel.

Void ratio and temperature distributions
The void ratio and temperature distributions are studied in different cylindrical cutting planes, at constant radius (see Figure 6), presented in a blade-to-blade view, for the operating point (φ=0.108, T ref =23K). This point is chosen because it corresponds to the strongest thermal effect.

Mid-span cut
The void ratio contours on a mid-span cutting plane are plotted in Figure 7 when the NPSH value decreases (from left to right). We observe a short cavity on the suction side, attached after the end of the bevel, which grows substantially. The cavity sheet thickness is relatively important (in comparison with the other case T ref =19K). The void ratio is lower than 30%.
We can see also a small cavity on the trailing edge bevel. The development of cavitation on the pressure side is clearly observed. When the NPSH decreases, the cavity sheet on the suction side reaches the throat, generates large blockage to the oncoming flow and the performance breaks down.
The temperature distribution on the same cutting plane, is showed in Figure 8. The temperature is well correlated to the void ratio, plotted in Figure 7. The cooling effect due to the vaporization is evidenced. The cooling maximum value estimated with the reference temperature is about 1K.

Near-shroud cut
The void ratio contours on a cutting plane close to the shroud are plotted in Figure 9.
In comparison with the mid-span contours, for high values of NPSH, the cavity sheets are thicker, the void ratio reaches 50%, and we do not observe cavity sheets on the pressure side.
Rapidly, an extended cavity interacts with the leading edge of the neighboring blade. For low values of NPSH (NPSH=30m), cavitation appears on the pressure side and the channel is completely obstructed. We observe small cavities on the trailing edge bevel, in which the maximum void ratio value is about 50%. When the cavitation develops, we can see that the maximum void ratio value moves from the leading edge to the trailing edge.
Moreover, for this case with strong thermal effects, during the NPSH decrease, we can note that the void ratio first increases then decreases. When the cavity becomes larger, the void ratio and the temperature depression increase. Then, when the cavity reaches the throat section of the inducer, the void ratio and the temperature depression decrease slightly, which is probably due to the interaction between the cavity and the flow around the leading edge of the adjacent blade. This situation was also commented in [23].
To compare, the void ratio contours, in the same cutting plane, are presented in Figure 10 for the operating point (φ=0.098, T ref =19K). For this point, thermal effects is less significant, therefore the void ratio is higher. Indeed, the maximum value is around 80%.
Small cavities appear on the pressure side (for NPSH=34m). For NPSH=27m, the cavity do not obstruct the channel. The performance breakdown is certainly due to the torque reduction.
The temperature contours close to the shroud, for the operating point (φ=0.108, T ref =23K), are presented in Figure 11. The temperature is well correlated to the void ratio, as noticed previously. The location of the maximum temperature depression corresponds to the location to the maximum values of void ratio. We observe a heating effect along the blade, especially on the pressure side and at the end of the trailing edge. We interpret this phenomenon as a viscous warming effect along the wall. To study it, the wall temperature is analysed on the next section.

Wall temperature
The wall temperature on the suction side is plotted in Figure 12, for the operating point For the non cavitating flow, we observe a viscous warming effect on the wall, which increases with the radius. It is around 1.2K close to the shroud, and it is on the same order of the temperature depression ∆T * . The same heating effect can be evidenced on the pressure side.
For the cavitating regime, there is a "competition" between these two phenomena : the viscous heating effect versus the vaporization cooling effect. We can clearly see the thermal "footprint" of the cavity sheet, attached after the end of the leading edge bevel (dotted line).
The lowest temperature is around 22K; it corresponds to a cooling effect, estimated with the reference temperature, between 0.5K (hub) and 1K (shroud). The absolute cooling value can reach 2K.
Due to the high rotation velocity, the Eckert number is high (about a factor 100 compared to an usual inducer in cold water) and the viscous heating becomes no more negligible. Values of Eckert and Mach numbers are presented in Table 2, for hub and shroud radii. Mainly close to the shroud, the Mach number is about 0.2 in the pure liquid, compressibility effects should be not influential.  To check this interesting result, especially regarding the grid used, a finer mesh has been built. The new grid is a HOH-type topology, which contains around 1457500 nodes. The y + values at the walls are on the same order than the coarse grid. Views of the leading and trealing edges are given in Figure 17 for both meshes.  Figure 18 shows the wall temperature on the suction side in non cavitating regime (on the left) and in cavitating regime (on the right). As previously, the opposition of a viscous 14 heating effect and a vaporization cooling effect are clearly observed.

Conclusions
Only few numerical models allowing the 3D simulation of viscous cavitating flows under cryogenic conditions exist in the literature. In the present work, a new homogeneous onefluid cavitation model taking into account thermal effects was proposed and implemented in the Fine/Turbo T M code. This original approach, based on a ρ(P, T ) state law, is an improvement of our previous works concerning a barotropic model and cold water applications [7,17]. In those studies, we have obtained a good prediction of the head drop charts for different pump geometries at several flow rates.
In the present paper, the studied application concerns LH2 cavitating flows in a NASA inducer geometry, running at very high rotation velocity.