Numerical Simulation of Turbopump Inducer Cavitating Behavior

In the present study a numerical model of 3D cavitating flows is proposed. It is applied to investigate the behavior of a spatial turbopump inducer in noncavitating and cavitating conditions. Experimental and numerical results concerning inducer characteristics and performance breakdown are compared at different flow conditions. The cavitation development and the spatial distribution of vapor structures within the inducer are also analyzed. The results show the ability of the code to simulate the quasisteady cavitating behavior of such a complex geometry. Discrepancies concerning the breakdown prediction are also discussed.


INTRODUCTION
To achieve high rotational speed and low inlet pressure, rocket engine feed pumps are equipped with an inducer stage, in which developed cavitation occurs at nominal point of operation.Cavitation observed in turbopump inducers is characterized by complex, three-dimensional, two-phase structures.These may appear as cavitation sheets attached to blades suction side or in rotating flows at the periphery of the runner (blade tip vortices or inlet back flow).They are usually composed of a liquid/vapor mixture, whose structure can only be described macroscopically by defining a void ratio.Moreover, the cavitation process is often unsteady.Experimental results point out two main types of cavitation instability: This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
(i) a self-oscillation behavior of cavitation sheets, whose mechanism was studied in cavitation tunnels and analyzed in detail by many authors (e.g., Kubota et al. [1], Le et al. [2], Stutz and Reboud [3]); (ii) a rotating cavitation pattern, mainly observed in inducers, showing different sizes of the cavitation structures in the different blade-to-blade passages of the machine and leading to super-or sub-synchronous instabilities (Kamijo et al. [4], Tsujimoto et al. [5], de Bernardi et al. [6]).
In this context, a numerical model, aiming to take into account the cavitation phenomenon in inducers geometry, Reference density (kg/m 3 ) Reference cavitation parameter τ τ/τ 0 has been developed in collaboration with SNECMA Moteurs and the French space agency CNES.The final objective is to provide assistance to the design and prevision of operating range of rocket engine turbopumps, taking into account the steady-state and unsteady effects of cavitation.The model is based on the commercial code FINE/TURBO, developed by NUMECA International (Hakimi [13]).To simulate the cavitation process, a single-phase flow model is implemented (Coutier-Delgosha et al. [14]) based on previous studies in 2D (Delannoy and Kueny [7], Reboud and Delannoy [15], Reboud et al. [16]).The present paper focuses on quasi-steady effects associated with cavitation in a four-blade turbopump inducer.The numerical experience achieved by calculations on centrifugal pumps (Hofmann et al. [17], Pouffary [18]) allows us now to treat that more complex geometry, and to obtain a simulation of its cavitating behavior.Experimental measurements performed in the CREMHyG laboratory, from Institut National Polytechnique de Grenoble, provide data to perform a first validation of the calculation.
Section 2 presents the characteristics of the inducer geometry and the experimental setup.The main features of the physical and numerical models are proposed in Section 3, as well as the calculation conditions.Noncavitating and cavitating computations are then reported in Sections 4 and 5, and comparisons with experiments are performed, mainly concerning head drop charts at three flow rates and cavity extension.Notations used are defined in Table 1.

INDUCER GEOMETRY AND EXPERIMENTAL SET UP
The considered turbopump inducer is a four-blade axial runner especially designed to operate in cavitating conditions (Figure 1).
The experimental study of the inducer behavior was investigated in the CREMHyG laboratory in water flow field conditions.Time-averaged pressure measurements were performed along the shroud with seventeen pressure taps, whose position is indicated in Figure 2, connected to pressure transducers.Relative accuracy of the pressure measurements is 0.2%.The mass flow rate is measured by an electromagnetic flow meter, with a relative accuracy of 0.3%, the rotation speed is controlled with a relative accuracy of 0.2%.The reference velocity is the tangential velocity at the blades tip.
The combination of inaccuracies on pressure and velocity measurements ensures on the elevation coefficient ψ and on the cavitation parameter τ a global uncertainty of about ±1%.
A large range of flow rates was investigated around the nominal point of operation in noncavitating and cavitating conditions.Noncavitating experiments consist in setting a high static pressure level, to avoid any apparition of vapor, and then to vary the flow rate from 0.6 Q n up to 1.4 Q n .In cavitating conditions, the flow rate is kept constant, and the static pressure is decreased slowly to enhance vapor development in the inducer and reach the performance breakdown.

PHYSICAL AND NUMERICAL MODEL
The main features of the physical and numerical models are summarized in the present paper.More details are given in by Coutier-Delgosha et al. [14].

Cavitation model
Cavitating flows are described by a single-fluid model, based on previous numerical and physical work developed in LEGI (Delannoy and Kueny [7], Reboud et al. [16]).This fluid is characterized by a density ρ that varies in the computational domain: when the density in a cell equals the liquid density (ρ l ), the whole cell is occupied by liquid, and if it equals the vapor one (ρ v ), the cell is full of vapor.Between these two extreme values, a liquid/vapor mixture, still considered as one single fluid, occupies the cell.The void fraction α = (ρ − ρ l )/(ρ v − ρ l ) can thus be defined as the local ratio of vapor contained in this homogeneous mixture.Velocities are assumed to be locally the same for liquid and for vapor.An empirical state law is used to manage the mass fluxes resulting from vaporization and condensation processes.That barotropic law ρ(P) links the density to the local static pressure.When the pressure is higher or lower than vapor pressure, the fluid is supposed to be purely liquid or purely vapor, according to the incompressible state law ρ = ρ l or ρ = ρ v .The two fluid states are joined smoothly in the vapor-pressure neighborhood.It results in the evolution law presented in Figure 3, characterized mainly by its maximum slope 1/c min 2 , where c min 2 = ∂P/∂ρ • c min can thus be interpreted as the minimum speed of sound in the mixture.

Numerical resolution
The numerical model of cavitating flows based on that physical description is developed from the commercial code FINE/TURBO developed by NUMECA International.FINE/TURBO is a three-dimensional structured mesh code that solves the time-dependant Reynolds-averaged Navier-Stokes equations.Time-accurate resolutions of the equations use the dual time stepping approach.Pseudo-time derivative terms are added to march the solution towards convergence at each physical time step.The range of application is extended to low-compressible or incompressible flows by introducing a preconditioning matrix (Hakimi [13]).
The discretization is based on a finite volume approach.Convection terms are treated by a second-order central scheme associated with artificial dissipation terms.The pseudo-time integration is made by a four-step Runge-Kutta procedure.The physical time-derivative terms are discretized with a second-order backward difference scheme.
The code resorts to a multigrid strategy to accelerate the convergence, associated with a local time stepping and an implicit residual smoothing.
The numerical model was adapted to treat the cavitation process (Coutier-Delgosha et al. [14]).The key point of this adaptation is the modification of the state law of the fluid.The applied barotropic law implies the simultaneous treatment of two different cases: the fluid is highly compressible in the liquid/vapor mixture (the Mach number can be as high as 4 or 5) and it is almost incompressible in the pure vapor or pure liquid areas.So the main difficulty consisted in managing these two different states of the fluid, without creating any spurious discontinuity in the flow field.Besides, cavitation phenomenon is characterized by very sharp and rapid processes.The density variations in time and space are smoothed to avoid numerical instabilities.

Calculation conditions
The model was applied to study the flow in the turbopump inducer geometry presented above.A 500 000-cell multiblock mesh of one single blade-to-blade channel was used (Figure 4).Conditions applied for the simulations are the following.(i) Turbulence model.We use for the simulations presented in this paper a Baldwin-Lomax turbulence model.
(ii) Boundary conditions.Constant velocity is imposed at the inlet of the suction pipe, about 2 diameters upstream from the runner, and a constant static pressure is imposed at the outlet.Laws of the wall are imposed along solid boundaries.The relative motion between the inlet pipe walls and the inducer is taken into account.The tip gap is also modelized.Periodic conditions are applied at the lateral boundaries of the blade-to-blade channels.
(iii) Initial transient treatment.First of all, a steady step is carried out, with a pseudo vapor pressure low enough to ensure non cavitating conditions in the whole computational domain.Then, the cavitation parameter τ is slowly lowered by increasing smoothly the pseudo vapor pressure at each new time step up to the physical value.Vapor structures spontaneously appear and grow during that process, in the regions of low static pressure.The final value of τ, depending on the outlet static pressure imposed, is then kept constant throughout the computation.
(iv) Cavitation parameters.The calculation was performed in hydrogen.The minimum speed of sound c min in the two-phase mixture, which is necessary to determinate the maximum slope of the barotropic law, is calibrated to ensure the similarity between the calculations (in H 2 ) and the experiments (in water).The parameter c min /V ref is thus the same in both cases, which leads for the calculations to a value c min = 34 m/s.That result is consistent with the celerity range calculated with Jakobsen's law (Jakobsen [19]).

NONCAVITATING RESULTS
Experimental tests and numerical calculations were performed considering a large range of flowrates, in noncavitating conditions.Figure 5 illustrates the computational result at nominal flow rate: the total pressure increase is represented on the hub. Figure 6 presents a comparison between the numerical and experimental performance charts.Directly based on the measurements performed, the head of the pump is defined as the difference between downstream and upstream static pressures.The quantities P , P tot , and ψ   have the following expressions: with P 0 an arbitrary reference value, where the inlet and outlet static pressure are both considered at the sensors location, in the vicinity of the shroud (sensors 1 and 16 on Figure 2), and the reference velocity is the tangential tip velocity.A reliable agreement with the experimental performance chart is obtained in the whole flow range considered.That result confirms that the numerical model well predicts the general noncavitating behavior of the inducer.

CAVITATION BEHAVIOR
Numerical simulations of the inducer are performed in cavitating conditions at three flow rates.We present first some qualitative results, consisting in a visualization of the vapor/liquid structures in the pump at nominal flow rate and for several τ values.Then, a quantitative analysis is performed by comparing the experimental and numerical results: head drop charts are studied at two flow rates, and the extension of the vaporized area near the shroud is compared to the experimental tip pressure measurements during the τ decrease.

Qualitative results
The head drop chart ψ s (τ ) obtained by the calculation at nominal flow rate is drawn on Figure 7.The green line corresponds to the apparition of vapor in the flow field.The eight red squares indicated on the chart are related to the eight visualizations of the cavitating flow field presented on Figure 8 hereafter.The quantity τ has the following expression: τ 0 is an arbitrary reference value. ( The head drop is correctly obtained by the computation: numerical instabilities associated with the simulation of the breakdown are controlled, and about 20% reduction of the inducer head is predicted. Figure 8 shows the development of cavitation corresponding to the eight operating points indicated on the chart.It illustrates the apparition and the growing of vapor/liquid areas on both faces of the blades, and finally the progressive filling of blade-to-blade channels by the vapor.Blades are colored in grey, the hub in blue, and the external shape of the two-phase areas (corresponding to a 5% void ratio) is colored in yellow.
It can be observed that vapor first appears in the vicinity of the shroud, near the leading edge.That phenomenon slightly affects the inducer performance, which starts to decrease, even if only a little volume of vapor is present in the flow field.Vapor remains localized near the shroud and in the tip gap (Figures 8a-8d) until τ reaches 0.3.That behavior is consistent with the smooth decrease of the inducer head observed at the same time.
The development of cavitation on the blade suction side (Figure 8e) directly induces a severe performance breakdown.This behavior can be explained by the important blockage generated immediately by the vapor in the upstream part of the blade-to-blade channels.It is a very sudden phenomenon (Figures 8f-8g), which leads to the presence of vapor structures, characterized by a low void ratio, from hub to shroud and from pressure side to suction side at the channel inlet.

Comparisons with experimental performance charts
Numerical head drop charts are compared to the experimental ones at three flow rates, namely, 0.8 Q n , Q n , and 1.2 Q n (Figure 9).The cavitation parameter τ corresponding to the performance breakdown is in all cases over-predicted by the model: the gap between experimental and numerical results equals ∆τ = 0.1.
That discrepancy must be considered circumspectly.In the case of centrifugal pumps (Hofmann et al. [17]), the order of magnitude of τ is notably higher, so that the same absolute error results in a much lower relative error (less than 5% in most of the simulations (Pouffary [18])).Since an inducer has much better suction capacities (τ at the performance breakdown is about 0.1), the order of magnitude of τ is in the present case the same as the one of ∆τ .
So, the over-estimation of τ at the breakdown is consistent with the absolute discrepancies observed in the calculation of centrifugal pumps.As a matter of fact, it is much more penalizing in the case of an inducer, and the simulation does not allow at the present time a correct prediction of the performance for such an axial runner.
The influence of the barotropic law maximum slope on the result obtained has been studied, by reducing the value of c min (minimum speed of sound in the mixture).Decreasing c min increases the maximum slope of the law, and thus theoretically reduces the cavitating areas.A wrong tuning of that parameter could be responsible for an over-prediction of the vapor volume in the blade-to-blade channels, and thus could explain the overestimation of τ breakdown .
c min was initially fixed to 34 m/s (value obtained by similarity rules between calculations in hydrogen and experiments in water).Another test was performed with c min = 17 m/s, and the comparison between the two results is proposed in Figure 10.
That result confirms that reducing c min has in the present case almost no influence on the global chart, and particularly on the head drop location.The reason of the discrepancy still remains to be found and a continuing work is pursued to improve the simulation.

Comparison of the vaporized areas extension
The size of the two-phase flow areas predicted by the calculation in the vicinity of the shroud is compared to the experimental results.The evolutions of the time-averaged pressure P measured at nominal flow rate are reported as functions of τ on Figure 11.Only 12 sensors are considered, since the first five ones are located upstream from the inducer leading edge (see Figure 2).Sensors 7 to 17 are all successively affected by the presence of vapor, as can be seen in the Figure 11.For each sensor, the value of τ corresponding to the presence of vapor can be obtained with an uncertainty of ±0.05.It leads to a curve τ appear (z ) where z = z/z 0 with z the distance between a plane z = constant and the leading edge (at the hub), and z 0 is an arbitrary reference scale, τ appear is the value of τ corresponding to the apparition of vapor in that plane.
The same type of curve is obtained from the numerical result, with a maximum uncertainty on each τ value of 0.01, and the comparison between the two curves is proposed in Figure 12.
An excellent agreement is observed between the two results: the discrepancies mostly lie beyond the uncertainty level, excepted for high values of τ .In that particular case, the calculation predicts very soon a presence of vapor at the leading edge, near the shroud, which is not reported by the experimental measurements (Figures 8a-8c).
That error must be associated with the very low pressure level that is predicted by the code in noncavitating conditions in the gap: this nonphysical prediction leads during the cavitating calculation to very high local density gradients, which generate strong numerical instabilities.Efforts are pursued to overcome this difficulty.
The good general agreement with the experimental results suggests that the axial vapor extension in the inducer is correctly simulated.The remaining discrepancy concerning  the breakdown location could thus result from physical aspects that are not modelized in the present calculation, such as unsteady effects for example.Numerical errors associated with the high aspect ratio of the mesh in some areas (leading edge, tip of the blades) could also be responsible for some discrepancies in the prediction of the static pressure field.They induce for example the systematic apparition of vapor near the leading edge, even for a high cavitation number, which is not observed experimentally.

CONCLUSION
A numerical model of 3D cavitating flows, based on the 3D code FINE/TURBO, has been developed to predict the cavitation behavior in turbomachinery (Coutier-Delgosha et al. [14,20]).Numerical results were presented in this study, concerning a four-blade inducer.Noncavitating and cavitating conditions were investigated.Calculations were found to be    in good agreement with experimental measurements in noncavitating conditions.The qualitative development of cavitation in the inducer was simulated and the extension of vapor near the shroud was compared to experimental results.A good general agreement was obtained.Experimental and numerical performance charts including in both cases the breakdown due to cavitation were also compared.The numerical uncertainty concerning the breakdown occurrence, which was found small in the case of centrifugal pumps (Hofmann et al. [17], Pouffary [18]), leads in the present application to an important relative error.The vapor extension as well as the breakdown mechanisms seem however to be correctly simulated by the numerical model.A continuing work is pursued to improve the prediction of the cavitation head drop in the particular case of inducers.

Figure 2 :
Figure 2: Position of the pressure taps along the shroud.

Figure 4 :
Figure 4: View of the mesh on hub side of the inducer (the entire inducer geometry is reconstructed by rotation of a single blade-toblade channel).

Figure 5 :
Figure 5: Total pressure increase on the hub at nominal flow rate.

Figure 6 :
Figure 6: Characteristics ψ s (Q/Q n ) of the inducer in noncavitating conditions.

Figure 7 :
Figure 7: Head drop chart at nominal flow rate.The red squares indicate the cavitating conditions visualized on Figure 8.The green line corresponds to the apparition of vapor.

Figure 8 :
Figure 8: Development of the two-phase areas as the cavitation number decreases (corresponding to red squares reported on Figure7).

Figure 11 :
Figure 11: Evolution of the pressure measured by the taps located near the shroud.

Figure 12 :
Figure 12:  Comparison between numerical simulation and experimental measurements.The nondimensional inducer axial length z is given in abscissa.(The reference z = 0 is the leading edge position near the hub.) τ , corresponding to the apparition of vapor in the plane z = constant, is reported in ordinate.