Numerical Simulation of 3D Cavitating Flows: Analysis of Cavitation Head Drop in Turbomachinery

The numerical simulation of cavitating flows in turbomachinery is studied at the Turbomachinery and Cavitation team of LEGI (Grenoble - France) in collaboration with the French space agency (CNES) and the rocket engine division of SNECMA Moteurs. A barotropic state law is proposed to model the cavitation phenomenon and this model has been integrated in the commercial CFD code Fine/TurboTM, developed and commercialized by Numeca International. The numerical aspects of the work are mainly focused on numerical stability and reliability of the algorithm, when introducing large density variations through the strongly non linear barotropic state law. This research conducted first to changes in the way preconditioning parameters are calculated. Internal flows in turbomachinery have been deeply investigated. A methodology allowing the numerical simulation of the head drop induced by the development of cavitation has been proposed on the basis of computations in inducers and centrifugal pumps. These simulations have allowed the characterization of the mechanisms leading to the head drop and the visualization of the effects of the development of cavitation on internal flows.


INTRODUCTION
Cavitation is an important phenomenon for the design of rocket engine turbopumps and remains a major problematic for the space industry. That is why the French space agency CNES and the rocket engine division of the SNECMA-Moteurs society have supported research for many years, in order to progress in the understanding and prediction of the mechanisms associated with cavitation. In this frame, one main purpose of the research performed at the Turbomachinery and Cavitation team of LEGI (Grenoble -France) consists of developing and validating a 3D RANS code allowing to simulate the cavitating operation of turbopumps. Therefore, a cavitation model was developed and implemented in the commercial CFD code Fine/Turbo TM developed by Numeca International. The applied numerical model is very similar to the one proposed by Kunz et al. [1][2] to calculate cavitating flows in centrifugal pumps. The main difference between the present work and those studies is the physical model applied to describe cavitation phenomenon. In this study, the barotropic cavitation model initially proposed by Delannoy and Kueny [3] was considered and will be presented in the first section of the paper.
Then we discuss about the modifications of the numerical algorithm that are needed when introducing large density variations through the strongly non linear barotropic state law. The main part of the paper will then focus on the applications of the model to one centrifugal pump and two different inducer geometries, for different operating points in non-cavitating and cavitating conditions. The observations of the pumps performance and associated flow fields led us to carry out a detailed analysis of the influence of cavitation in the three different geometries. Finally, from these results, the origin of the cavitation head drop is discussed and an analysis of its mechanisms is proposed.

BAROTROPIC STATE LAW
We use a single fluid model to describe the liquid-vapor mixture, with a mixture density ρ varying in the flow field between the vapor density and the liquid density. The equivalent void ratio α of the liquid-vapor mixture relates to the varying specific mass by ρ = α ρ v + (1 -α) ρ l . In the mixture, the velocities of liquid and vapor phases are the same and we obtain only one set of equations for the mixture mass, momentum, or turbulence (k-ε), written in their conservative form.
The void ratio α of the mixture depends on the local static pressure. As initially proposed by Delannoy and Kueny, [3], we use a barotropic state law ρ(P) to manage the relation between pressure and mixture density. A smooth arbitrary law was chosen, ρ rapidly varying between liquid density ρ l and vapor density ρ v when the local static pressure P is around the vapor pressure Pv. The law is characterized by its maximum slope at P = Pv, which is related to the minimum speed of sound c min in the two-phase homogeneous medium [4]. In the present work, the ratio ρ v /ρ l is 0.1. Based on experimental previous work concerning Venturi type sections (Stutz and Reboud [5,6]), the parameter c min is taken equal to 2.35 m/s. That order of magnitude is kept constant for different applications with water.

NUMERICAL MODEL
The barotropic cavitation model was implemented in the commercial code Fine/Turbo TM developed by NUMECA International. Fine/Turbo TM is a three-dimensional structured mesh code that solves the time dependant Reynolds-averaged Navier-Stokes equations, with artificial compressibility method. A detailed description of the initial code is given in Hakimi  The space discretization is based on a finite volume approach. A second order central scheme is associated with two artificial dissipation terms, respectively of second and fourth order. The first one is activated in the strong pressure and density gradient areas. The other one is used in the whole domain, and it results in a second order space accuracy.
In the case of low-compressible or incompressible flows, the pseudo-time derivatives are multiplied by a preconditioning matrix, based on the studies of Turkel, [8], Choi and Merkle, [9]. The modifications of the preconditioning technique developed for cavitating flows will be discussed hereafter. 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 that ensures a second order accuracy in time. In non-cavitating conditions, the code resorts to a multi-grid strategy to accelerate the convergence associated with a local time stepping and an implicit residual smoothing.

Modification of the preconditioning technique
Under the preconditioned dual time stepping form, the set of equation becomes, as detailed by Pouffary [10]: introducing the preconditioning matrix Γ -1 and associated variables vector Q: The eigenvalues of the preconditioned system in noncompressible flow are: V is the velocity vector, n the vector normal to the surface dS. α is generally chosen equal to -1, and the spectral radius (maximum eigenvalue) is : In the compressible case, with the speed of sound c, the last eigenvalues are modified in: The second order artificial dissipation term and the local time step, respectively proportional and inversely proportional to the spectral radius, are then modified in the cavitating regions, where ρ d dp c = is calculated from the barotropic state law.
The preconditioning parameter β is expressed using a reference velocity U ref : To ensure the numerical stability in cavitating regions β must be reduced, but to save computational time this reduction is applied only in those regions. Two different ways to reduce β parameters in cavitating zones have been investigated (Pouffary et al. [11]): -the first method is based on a "local velocity scaling" strategy. In this case, U loc is the local velocity, U min is the minimum velocity allowed to ensure sufficient numerical dissipation, and 3 2 0 = β is a constant. As a matter of fact, the relative velocity in the cavitating zones, close to the pump blades, is often smaller than U ref , and the coefficient β is reduced as a function of the local velocity values.
-the second method consists in modifying the reference velocity U ref with respect to the local static pressure and the barotropic state law. The relation β(P) drawn on figure 2 was proposed and the influence of the different curve parameters was widely investigated. Both methods have been tested on different cases in 2D and 3D steady or unsteady flows (Pouffary et al. [12], [13]). The different computations presented hereafter have been performed in the steady case, without taking into account the real time derivative term. Depending on the case, one of the two methods generally led to a better stability and was then kept for all the computations performed on the given geometry. It is worth noting that these numerical modifications led to major improvements in accuracy and numerical stability in comparison with our previous computations of cavitating flows in turbomachinery (for example in Coutier-Delgosha et al., [14,15,16] or Pouffary et al. [17]).

Centrifugal pump
The first considered geometry is the SHF centrifugal pump tested by Combes and Archer [18]. The SHF pump has a shrouded seven-blade centrifugal impeller. The scale and the calculation configuration are chosen as close as possible to the experimental conditions to simplify comparisons in cavitating conditions. The outlet radius is 330 mm, the inlet pipe length is 200 mm, the rotation speed equals to 3000 rpm, and the nominal flow rate is 0.157 m 3 /s. After comparison of different meshes in non cavitating conditions, a mesh of 250000 cells for one blade-to-blade channel was chosen (Fig. 3).

Turbopump inducers
To achieve high rotational speed and low inlet pressure, rocket engine pumps are equipped with an axial stage. Two different four-blade inducers have been tested. Inducer-1 was designed by SNECMA Moteurs (Fig. 4a). The influence of the mesh and of the turbulence model (i.e., Baldwin-Lomax and three different k-ε models) was firstly studied in non cavitating conditions, mainly by comparing the static pressure profile at the shroud with experimental measurements performed at different flow rates at the CREMHyG Laboratory (Grenoble -France). An example at nominal flow rate is given on figure 5. From these comparisons, a mesh of about 500,000 cells for one blade-to-blade channel was chosen, associated with a k-ε model with wall function. Tip leakage was taken into account in noncavitating conditions, but to increase the convergence speed it was suppressed for the cavitating flows computations, by locally modifying the mesh. Analysis of the influence of the tip leakage on cavitating flows is in progress. The geometry of Inducer-2 shows a pure radial leading edge shape (Fig. 4b) and has been studied experimentally at the University of Osaka [19]. The mesh has 250000 cells for one blade-to-blade channel, and do not take into account the tip leakage.

ANALYSIS OF CAVITATING FLOWS
Steady cavitating flows were calculated in these geometries for different cavitation numbers. Two sets of upstream and downstream boundary conditions have been compared: for the first one, flow rate is imposed upstream and static pressure downstream (with or without the radial equilibrium condition), for the second one, total pressure is imposed upstream and flow rate downstream. Influences of the rotation of hub and shroud boundaries upstream and downstream from the impeller have also been analyzed. With a good convergence level of the computations, no major influences of the boundary conditions chosen were observed on the results.
In cavitating conditions, the boundary conditions are not modified. The different cavitation numbers (or corresponding NPSH) are obtained by varying the pseudo-vapor pressure P v involved in the definition of σ, while the reference pressure is kept constant. It is worth noting that each calculated steady state solution was carefully converged. Through the modifications of the numerical scheme, a very good convergence level can be obtained, even in the cases of very large cavitation structures and associated strong modifications of the pumps head. Indeed, for the three pumps, operating points in cavitating regime leading to 10% head drop or more could be accurately simulated.

Centrifugal pump
Numerical simulation of the SHF centrifugal pump was performed using the "Local velocity scaling" method. The cavitating behavior obtained at nominal flow rate Qn shows thin attached sheet cavities on the suction side of the impeller, their length increasing with the decreasing cavitation number (Fig. 6). Head begins to drop approximately when the cavitation sheet reaches the throat, and decreases rapidly from more than 10% of its non cavitating value (Fig. 7). The numerical prediction of head breakdown is found for the critical NPSH close to 6.3 m, which corresponds to the NPSH 3% found experimentally at Electricité de France [18]. Flow rates over the nominal one lead to larger critical NPSH, head breakdown being observed at 1.22 Qn and 1.3 Qn for NPSH values respectively about 18 m and 26 m.

Experimental result (Qn)
Numerical results obtained at nominal flow rate were analyzed in detail. It can be seen first on figure 8, that the head drop can be mainly associated to a strong decrease of the torque (for the critical NPSH value, the calculated torque decreases from more than 10% of the non cavitating value). Otherwise, the cavitation head drop corresponds to a much smaller decrease of the hydraulic efficiency, which is found only 2% smaller during the head breakdown than in the non cavitating regime.
The torque decrease can be illustrated by Figure 9, where the static pressure distributions around a blade section at midspan are drawn for different NPSH values. When the net pressure Pref-P v becomes small, a long sheet cavity develops along the blade suction side, reducing the suction effect near the blade leading edge, without modifying significantly the pressure distribution along the wetted surfaces, and thus clearly reducing the blade load. At head breakdown, it can be observed that cavitation also appears on the upstream part of the blade pressure side (red curve). Such observations can also be made along other sections of the pump blades. An analysis of the flow within the pump have been performed in eight different flow sections, from cut I near the blade leading edge, to cut VIII near the trailing edge ( fig. 10). First of all, the total pressure rise between successive sections is drawn. One can clearly see on Figure 11 that the total pressure breakdown occurs principally in the upstream sections I to III, while the downstream ones (IV to VIII) remain approximately unaffected by the development of cavitation.   (1 to 6).
A local analysis of the influence of cavitation on the flow structure in the different sections has been also performed. For example, the secondary flows can be characterized by the dimensionless relative helicity, drawn for the four upstream sections in figures 12a in non cavitating condition and 12b near to the cavitation breakdown. Blue colors represent anticlockwise vorticity and yellow to red ones indicate clockwise vorticity. In non cavitating conditions, one can observe that the secondary flows in the rotating frame develop mainly along the blade suction side (blue color), progressively evolving to a hub vortex (section IV). Clockwise vorticity (red color) can also be observed along the shroud. In cavitating conditions, secondary flows are strongly modified along the cavitation sheet surface (suction side, section I). Then the cavity wake generates a large anti-clockwise structure that can be observed in sections II and III, and is progressively dissipated.  Comparison of the numerical cavitation performance chart with experiments performed at the CREMHyG Laboratory is drawn on figure 14. For a high cavitation number the numerical value is found about 5% higher than the experimental inducer head, which can be attributed mainly to the fact that the tip leakage flow is not taken into account by numerical calculations. Successive cavitating flows are computed and a good convergence can be obtained even with more than 10% head loss. The critical cavitation parameter is very well predicted. Calc.

Fig. 14: Cavitation head drop curves at nominal flow rates. Comparison between computation and experiments performed at the CREMHyG Laboratory.
Static pressure distributions around the blade are drawn on figure 15 in non cavitating and strongly cavitating conditions: the blade load decrease is mainly due to the modification of the pressure side distribution near the leading edge, observed when the cavitation sheet developed on the previous blade suction side reaches the entry of the blade-to-blade channel. Perturbations of the trailing edge pressure distribution can also be observed. Those perturbations can be associated to the important cavity wake entering into the blade-to-blade channel. The total pressure rise in the inducer, drawn on figure 16, shows a very progressive behavior, which agrees with the fact that the blade loads remain quite constant from the leading edge to trailing edge of the blades ( fig. 15). During head drop the pressure rise keeps mainly that regular behavior.  (1 to 8).
The secondary flows in cut sections IV and V (presented on figure 13) are analyzed in figure 17. In non cavitating condition, the main effect that can be noticed is the centrifugation of the boundary layer along the blade pressure side, leading to clockwise vorticity (red color). In cavitating conditions, near to the head breakdown, the important cavity wake structure can be clearly observed approximately in the middle of section V.

Inducer -2
The second inducer geometry exhibits a quite different behavior in cavitating condition. First of all, as illustrated by figure 18, the development of the cavitation structures appears much more pronounced at the inducer periphery. We can observe also a large back-flow structure in front of the blades. On the blade suction side, the cavitation sheets show a triangular shape, increasing from hub to periphery, and merging with the cavitating back-flow structure.
On the performance chart, figure 19, the head drop appears much more progressively than for Inducer-1. The last point calculated shows a 20% head drop, mainly due to the decreasing of the torque, as for the centrifugal pump (the efficiency decrease is about 3%). The fact that the cavitation head drop is much more progressive can be explained by looking to the blade load evolution drawn on figure 20, and comparing with figure 15. In non cavitating condition, the blade loads of the two inducers appear very different: Inducer-2 blades are loaded mainly near to the leading edge, while the load remains approximately constant along the Inducer-1 blades. In cavitating regime, the development of the sheet cavity from the leading edge of the Inducer-2 blade suction side modifies progressively the blade load, while in Inducer-1 the blade load decreases only when the cavity closure reaches the blade to blade passage. It can clearly be seen on figure 21 that the main part of the total pressure rise in Inducer-2 is provided by the upstream part of the inducer. As a consequence, cavitation head drop occurs mainly with a very progressive decrease of the pressure rise between section I and II.

CONCLUSION
A numerical model of 3D cavitating flows, based on the CFD code Fine/Turbo TM , has been developed to predict the cavitation behavior in turbomachinery. Recent improvements of numerical stability and reliability of the algorithm were performed through modification of the preconditioning technique. As a consequence, very good convergence and accuracy of the numerical simulation of cavitating flows in turbomachinery can now be obtained, even for operating conditions very close to the machine head drop.
Numerical results were presented and analyzed for a centrifugal pump and two different four blades inducers. Noncavitating and cavitating conditions were investigated and compared. The vapor extension and associated pumps head breakdown seem to be correctly simulated by the numerical model and compare favorably with experimental results. Local analyses of the results provide interesting information on the breakdown mechanisms: -in the three geometries, head breakdown is clearly associated with a decrease of the torque provided by the impeller, through a strong variation of the blade load. Otherwise, only a small decrease of the machine efficiency can be associated with the development of the cavitation structures.
-an important difference between the cavitation behavior of the two inducers can be observed, that can be explained by an analysis of their respective blade load.
-local visualizations of the flow fields and their evolution with cavitation can also be obtained. For example, the analysis of dimensionless relative helicity maps gives useful information about secondary flows in the machines and allows identifying the location and influence of the cavity wakes in the blade-toblade channels.
It is now demonstrated that numerical simulation can provide useful information for the design of turbomachinery, and more particularly of inducers. Applications performed on unsteady flows are in progress and first results presented in Pouffary et al. [18] are very promising for the analysis and prediction of unstable cavitation behavior in turbomachinery. Improvement of the cavitation model will be proposed in the future, in particular to take into account the thermodynamic behavior of the fluids as hydrocarbons or cryogenic propellants.