Benoît Pouffary

Regiane Fortes

Jean-Luc Reboud

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/Turbo TM , 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. [START_REF] Kunz | A Preconditioned Navier-Stokes Method for Two-phase Flows with Application to Cavitation[END_REF][START_REF] Medvitz | Performance analysis of cavitating flow incentrifugal pumps using multiphase CFD[END_REF] 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 [START_REF] Delannoy | Two phase flow approach in unsteady cavitation modelling[END_REF] 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.

NOMENCLATURE

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, [START_REF] Delannoy | Two phase flow approach in unsteady cavitation modelling[END_REF], 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 [START_REF] Pouffary | Numerical analysis of cavitation instabilities in inducer blade cascade[END_REF]. In the present work, the ratio ρ v /ρ l is 0.1. Based on experimental previous work concerning Venturi type sections (Stutz and Reboud [START_REF] Stutz | Two-phase flow structure of sheet cavitation[END_REF][START_REF] Stutz | Measurements within unsteady cavitation[END_REF]), 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 [START_REF] Hakimi | Preconditioning methods for time dependent Navier-Stokes equations[END_REF]. Time accurate resolutions use the dual time stepping approach. Pseudo-time derivative terms are added to the equations. They march the solution towards convergence at each physical time step. The initial set of mass and momentum equation can then be written under the following form:

∫∫∫ ∫∫ ∫∫∫ ∫∫∫ = ⋅ + ∂ ∂ + ∂ ∂ V S V V STdV dS F UdV t dV U τ (1) where ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = w v u U ρ ρ ρ ρ is the vector of the conservative variables, v c F F F - =
notes the convective and diffusive fluxes across the surface dS = n dS, ST is the source term, and τ the pseudo-time.

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, [START_REF] Turkel | Preconditioning methods for solving the incompressible and low speed compressible equations[END_REF], Choi and Merkle, [START_REF] Choi | The application of preconditioning in viscous flows[END_REF]. 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 [START_REF] Pouffary | Simulation numérique d'écoulements 2d/3d cavitants, stationnaires et instationnaires : analyse spécifique pour les turbomachines[END_REF]:

∫∫∫ ∫∫ ∫∫∫ ∫∫∫ = ⋅ + ∂ ∂ + ∂ ∂ Γ - V S V V STdV dS F UdV t dV Q τ 1 (2)
introducing the preconditioning matrix Γ -1 and associated variables vector Q:

⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = w v u p Q g and ( ) ( ) ( ) ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + + + = Γ - ρ β α ρ β α ρ β α β 0 0 1 0 0 1 0 0 1 0 0 0 1
) ( ) ( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + - ⋅ ± - ⋅ 2 2 2 4 1 1 2 1 β α α n n V n V (3)
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 :

2 2 max β λ + + = V V . ( 4
)
In the compressible case, with the speed of sound c, the last eigenvalues are modified in:

( ) ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⋅ - + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - ⋅ ± ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - ⋅ 2 2 2 2 2 2 2 2 2 4 1 1 2 1 β β α β α c n V n c n V c n V (5)
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 : -the first method is based on a "local velocity scaling" strategy. In this case, Both methods have been tested on different cases in 2D and 3D steady or unsteady flows (Pouffary et al. [START_REF] Pouffary | Numerical simulation of cavitating flow around a 2D hydrofoil: a barotropic approach[END_REF], [START_REF] Pouffary | Unsteady cavitating flows in turbomachinery: comparison of two numerical models and applications[END_REF]). 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., [START_REF] Coutier-Delgosha | Numerical simulation of turbopump inducer cavitating behavior[END_REF][START_REF] Coutier-Delgosha | Experimental and numerical studies in a centrifugal pump with 2D-curved blades in cavitating conditions[END_REF][START_REF] Coutier-Delgosha | Cavitation performance of a centrifugal pump: numerical and experimental investigations[END_REF] or Pouffary et al. [START_REF] Pouffary | Evaluation of the effects of cavitation on the flow in a centrifugal pump: analysis of the breakdown mechanisms[END_REF]).

STUDIED GEOMETRIES Centrifugal pump

The first considered geometry is the SHF centrifugal pump tested by Combes and Archer [START_REF] Combes | Etude de la cavitation dans la pompe SHF[END_REF]. 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).

Fig.3: Mesh of the SHF centrifugal pump.

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. Non-cavitating conditions.

The geometry of Inducer-2 shows a pure radial leading edge shape (Fig. 4b) and has been studied experimentally at the University of Osaka [START_REF] Yoshida | Effects of Alternate Leading Edge Cutback on Unsteady Cavitation in 4-Bladed Inducers[END_REF]. 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 [START_REF] Combes | Etude de la cavitation dans la pompe SHF[END_REF]. 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.

Meridian coordinate

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.

Fig. 11: Repartition of the total pressure rise between the sections I to VIII of the SHF pump, in non cavitating condition (nc: left) and with decreasing NPSH (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.

3m). Plane I to IV

Inducer-1 Development of cavitation structures in the Inducer-I geometry starts from the blade tip, and extend rapidly along the leading edge to the hub side. Then the cavitation sheet length increases progressively until it reaches the throat (fig. 13). Head drop then occurs abruptly when the cavity enters into the bladeto-blade channel.

Fig. 13: Inducer-I, nominal flow rate, for decreasing cavitation coefficient. Blue color corresponds to a 5% equivalent void ratio. (Cut sections are numerated from I to VIII, VII corresponding to the hidden trailing edge of blade A).

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. 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. [START_REF] Combes | Etude de la cavitation dans la pompe SHF[END_REF] 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.

Fig. 1 :

 1 Fig. 1: Barotropic state law ρ(P). Water 20°C.

  -p ref is the gauge (relative) pressure, α and β are the preconditioning parameters.The eigenvalues of the preconditioned system in noncompressible flow are: n V ⋅ and

(

.

  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.[START_REF] Pouffary | Numerical simulation of cavitating flows in space pumps[END_REF]):

  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 figure2was proposed and the influence of the different curve parameters was widely investigated.

Fig. 2 :

 2 Fig.2: Evolution of the parameter β as a function of the static pressure and comparison with the barotropic state law ρ(P)

Fig. 4 :Fig. 5 :

 45 Fig.4: Inducers geometries: a -Inducer-1 b-Inducer-2

Fig. 6 :Fig. 7 :Fig. 8 :

 678 Fig. 6: SHF centrifugal pump, nominal flow rate, for decreasing NPSH = (P upstream -P v )/(ρg) ≈ 9.5, 7 and 6.3 m.Yellow color corresponds to a 5% equivalent void ratio.

Fig. 9 :

 9 Fig. 9: Nominal flow rate. Blade load at mid span, for decreasing NPSH ≈ 15m, 9.5m, 7m and 6.3 m (with about 8% head drop).

Fig. 10 :

 10 Fig. 10: Location of the analyzed flow sections in the SHF pump.

Fig. 12a :Fig. 12b :

 12a12b Fig. 12a : Secondary flows : dimensionless relative helicity, varying from -1 (blue, maximum anti-clockwise vorticity) to 1, red, maximum clockwise vorticity) Non cavitating flow. Plane I to IV

Fig. 15 :

 15 Fig. 15: Nominal flow rate. Inducer-1. Blade load at mid span, for decreasing cavitation parameters.

Fig. 16 :

 16 Fig. 16: Repartition of the total pressure rise between the sections I to VIII of Inducer-1, in non cavitating condition (nc: left) and with decreasing NPSH (1 to 8).

Fig. 17 :

 17 Fig. 17 : Secondary flows (dimensionless relative helicity), Non cavitating (a) and cavitating (b) flows. Plane IV and V.

Fig. 18 :

 18 Fig. 18: Inducer-2, nominal flow rate, for decreasing cavitation coefficient. Blue color corresponds to a 5% equivalent void ratio. (Cut sections are numerated from I to VIII, VII corresponding to the hidden trailing edge of blade A).

Fig. 19 :

 19 Fig. 19: Cavitation head drop curve.Nominal flow rate, Inducer-2.

Fig. 20 :

 20 Fig. 20: Nominal flow rate. Inducer-2. Blade load at mid span, for decreasing cavitation parameters.

ACKNOWLEDGMENTS

This research was supported by the French space agency CNES (Centre National d'Etudes Spatiales) and SNECMA Moteurs. The authors wish also to express their gratitude to Electricité de France and CREMHyG for providing experimental results, to prof. Tsujimoto for providing the Inducer-2 geometry, and to Numeca International for its cooperation to the development of the numerical code.


Benoît Pouffary

Regiane Fortes

Jean-Luc Reboud

Pierre-Alain Lambert

Numerical Analysis of Cavitation Instabilities in Inducer Blade Cascade

The cavitation behavior of a four-blade rocket engine turbopump inducer was simulated by the computational fluid dynamics (CFD) code FINE/TURBO™. The code was modified to take into account a cavitation model based on a homogeneous approach of cavitation, coupled with a barotropic state law for the liquid/vapor mixture. In the present study, the numerical model of unsteady cavitation was applied to a four-blade cascade drawn from the inducer geometry. Unsteady behavior of cavitation sheets attached to the inducer blade suction side depends on the flow rate and cavitation number �. Numerical simu-lations of the transient evolution of cavitation on the blade cascade were performed for the nominal flow rate and different cavitation numbers, taking into account simultaneously the four blade-to-blade channels. Depending on the flow parameters, steady or unsteady behaviors spontaneously take place. In unsteady cases, subsynchronous or supersynchronous regimes were observed. Some mechanisms responsible for the development of these instabilities are proposed and discussed.

Introduction

One of the most prejudicial consequences of the cavitation in the rocket engine turbopump inducers is the generation of system and machinery instabilities. Machinery instabilities may appear for a particular range of the cavitation number, near to the performance drop. These unsteady phenomena are associated with different rotating nonsymmetrical cavitation patterns, which are characterized by sub-or supersynchronous frequencies.

Cavitation instabilities lead to downstream and upstream pressure fluctuations, which may be at the origin of system instabilities. They induce also axial and radial forces on the blades and shaft, which can disturb the bearings working and endanger the turbopump running.

For these reasons, inducer cavitation instabilities have been studied for many years in the aim of observing, better understanding, quantifying, and predicting such phenomena. Extensive studies have been carried out mainly from the viewpoint of system instability ͓1-3͔. Several experimental works have been performed to visualize global cavitation patterns and to measure the radial loads on the shaft due to cavitation instabilities ͓4-8͔. More recently, theoretical and numerical models have been developed to predict these instabilities and to analyze their origin ͓9-15͔.

In spite of those relevant works, more detailed cavitating flow analyses are needed for the understanding of local cavitation instabilities.

Through collaborations between the laboratory LEGI-Laboratoire des Ecoulements Géophysiques et Industriels, the Rocket Engine Division of Snecma, the French Space Agency CNES, and NUMECA International, a cavitation physical model was implemented in the commercial code FINE/TURBO™. The code allows steady and unsteady cavitating flows calculations in twodimensional ͑2D͒ or three-dimensional ͑3D͒ geometries, as presented in Refs. ͓16-20͔, and it can be used as a tool to analyze inducer instabilities.

In this context, and in complement to previous work ͓13,14͔, the main target of the present study is to propose a qualitative analysis of cavitation instabilities. The presented analyses are based on some numerical results obtained with the unsteady numerical code in a blade cascade geometry, derived from a 3D inducer. Subsynchronous and supersynchronous cavitation regimes have been simulated and analyzed. This paper illustrates numerical results obtained and proposes qualitative analyses of mechanisms responsible for these instability phenomena.

Models

The applied numerical code solves the conservative unsteady Reynolds-averaged Navier-Stokes ͑RANS͒ equations of a homogeneous fluid ͓16-21͔. In the present study, the Spalart-Allmaras model was used to simulate turbulent behavior ͓22͔.

Associated with these basic equations, the barotropic state law proposed by Delannoy and Kueny ͓23͔ was implemented in the code to model mass transfer in cavitating flows. This simple model involves in the hypotheses of instantaneous vaporization and condensation phenomena, as well as the no-slip condition between liquid and vapor phases. It does not take into account thermal effects, i.e., the energy equation is not considered in the calculation. The liquid-vapor mixture is characterized by a mixture specific mass given as a function of void ratio ␣, by the relation

= ␣ v + ͑1-␣͒ l .
This barotropic law depends on the minimum speed of sound, c min , in the mixture and on the value of the density ratio, v / l , as illustrated in Figs. 1 and2.

In the present study, the Reynolds number is in the order of 10 6 . The value of the minimum speed of sound c min is taken equal to 4.7 m / s, which is higher than the values used in our numerical previous works concerning cavitating flows in a Venturi apparatus ͓24͔. The density ratio v / l is imposed to be 0.1. These choices were made to improve the robustness of these first feasibility computations with FINE/TURBO™ code.

For the same blade cascade geometry considered in the present work, some numerical tests were carried out in a previous work to evaluate the influence of the speed of sound c min , the Reynolds number, and the used time step and mesh on the cavitation instability pattern. Most important results obtained from these influence analyses are described in Ref. ͓13͔.

Regardless of considered hypotheses, previous studies ͓24-27͔ pointed out that the barotropic approach can succeed to quantita-tive and qualitative predictions of cavitating flow global parameters ͑i.e., characteristic frequencies, vapor structure size, and pump performance͒ as well as mean local parameters ͑void ratio and velocity distribution͒.

These physical models have been implemented in the code FINE/TURBO™ developed by NUMECA International. This threedimensional structured mesh code solves the time dependent Reynolds-averaged Navier-Stokes equations. Time accurate resolutions use the dual time stepping approach proposed by Jameson ͓28͔. Pseudotime derivative terms are added to the equations. They march the solution toward convergence at each physical time step.

This kind of resolution is devoted to highly compressible flows. In the case of low-compressible or incompressible flows, its efficiency decreases dramatically. This well-known problem has been addressed by many authors and solved by introducing a preconditioner ͓21͔. This one is based on the studies presented in Refs. ͓29,30͔, and has been modified to take into account the cavitation model.

The discretization is based on a finite volume approach. We use a second order central scheme that must be associated with two artificial dissipation terms, respectively, of second and fourth orders as initially proposed by Jameson ͓28͔. The first one is activated in the strong pressure and density gradient areas. The other one is used in the whole domain. The applied form of the artificial dissipation term leads to a central second order accurate convection scheme. The pseudotime 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.

A meaningful numerical work has been performed recently by Pouffary ͓17͔ to improve the preconditioner and the stability of numerical code for calculations of cavitating flows. A more detailed description of the code is given in Refs. ͓19-21͔. Some influence tests of numerical parameters, including the effects of second and fourth order dissipation terms, are presented in Ref.

͓19͔.

An example of studied inducer geometry is given in Fig. 3. Because 3D unsteady calculations in inducer geometries are very time consuming, a 2D approached geometry was adopted ͑even if it cannot represent completely the original three-dimensional case͒.

The studied inducer was designed to operate, even at a nominal flow rate, with nonzero incidence angle. Nevertheless, in the range of running conditions, the angle of attack remains very small for the considered inducer geometry, because of the variable thread and the associated curvature of the blades. In these conditions, the inlet backflow is not relevant and a 2D geometry approach can be used for qualitative analyses. It is important to note that tip clearance and backflow effects are neglected and, consequently, the applicability of the approach proposed in this paper might be limited to high inducer flow rates.

Hence, computations are performed in a four-blade cascade derived from an entire inducer. The transformation from the 3D geometry to the 2D blade cascade leads us to neglect the peripheral cavitation in the inducer: Only cavitation sheets attached on the blades will be considered. To obtain the 2D geometry, the 3D inducer ͑Fig. 3͒ is cut at a constant radius R c equal to 70% of the tip radius R. The resulting shape of the blade cut in the plane ͑R c , Z͒ is given in Fig. 4.

We use a 215ϫ 31 structured mesh per blade-to-blade channel, giving 26,660 internal nodes when calculating the four-blade cascade. The boundary conditions are imposed velocity at the mesh inlet, imposed static pressure at the outlet, and periodicity or connection conditions between the different channels of the blade cascade ͑Fig. 4͒. The distributions of the inflow velocity and the downstream pressure are uniform.

The time step, mesh, and turbulence model are chosen to put the attention on the low frequency fluctuations of the attached cavity, more than to the local unsteadiness in the cavitation sheet wake ͑cloud shedding͒. v Several four-channel computations are performed at the nominal flow rate, for different cavitation numbers ͑based on outlet reference pressure͒ varying from low cavitating conditions down to the final performance drop of the cascade. The corresponding head drop chart is drawn in Fig. 5 at a nominal flow rate.

Stable configurations correspond to cavitation numbers lower than 0.65 or higher than 0.8 ͑square dots in the figure͒. For Ͼ 0.8, cavitation sheets are small and identical on the four blades. For these flow conditions, the performance of the runner is only slightly affected by the presence of vapor. For Ͻ 0.65, the cavitation sheets are much more developed, and they are responsible for the important performance drop observed. The performance drop of the cascade is directly related to the evolution of the blade load illustrated in Fig. 6: The appearance of cavitation structures Circular dots in Fig. 5 correspond to unsteady regimes and define the range of cavitation instabilities of the cascade. Three kinds of instabilities have been observed from our calculations: a supersynchronous behavior, a subsynchronous configuration, and an unstable alternate cavitation.

Calculations are done in relative referential: Sub-or supersynchronous configurations are determined by the unsteady visualization of the cavitation sheets. For example, if the nonsymmetrical cavitation sheet pattern rotates in the same direction as the inducer, the characteristic frequency of the instability is added to the inducer rotation frequency and a supersynchronous configuration is established. In another way, subsynchronous regimes are characterized by cavitation patterns that rotate in the opposite direction as the inducer.

v

For downstream ϳ 0.75, two kinds of cavitation configurations have been observed. The first one is the alternate blade cavitation ͑Fig. 7͒. The system here is symmetrical, characterized by two radial opposite large vapor sheets and two opposite small ones.

For this cavitation number, the calculations pointed out also the appearance of another cavitation configuration, associated with a supersynchronous behavior. Calculations indicate that this cavitation pattern rotates in the blade cascade in the same direction as the inducer rotation, but more quickly.

v

For downstream ϳ 0.7, another type of alternate blade configuration has been observed. As illustrated in Fig. 9, for this cavitation number, the cavitation pattern is initially characterized by two opposite large sheets in Channels 1 and 3, separated by two channels without cavitation. During two periods of inducer rotation, this configuration shifts and the large vapor structures change to Channels 2 and 4.

This kind of configuration was studied in a previous work presented in Ref. ͓13͔.

For smaller cavitation numbers ͑ϳ0.65͒, a subsynchronous configuration appears ͑Fig. 10͒. In this case, the cavitation pattern is not symmetrical: We observe four sheets of different sizes. The largest sheet occurring in Channel 3 obstructs the flow in this channel and leads to the increase of the cavitation sheet in the upper channel. This configuration rotates in the opposite direction as the inducer rotation.

v

Based on the results presented here above, we propose in this paper a qualitative analysis of the mechanisms responsible for the instability phenomena observed.

According to our calculations, and by making an analogy with rotating stall phenomenon observed in compressor inlet ͓31͔, the subsynchronous configuration seems to be associated with the obstruction of a channel by the largest sheet, as schematized in Fig. 11. This obstruction deflects the flow and modifies the angle of attack of the neighboring upper channel, which leads to a pressure decrease and to the increase of the cavitation sheet size in this channel. Otherwise, as flow rate is imposed as a constant in the inducer, the flow angle of attack in the channel below the obstructed channel should decrease, which leads to the reduction of vapor structure size in this channel. Hence, this kind of rotating cavitation configuration would be strongly coupled with the modification of the angle of attack due to the flow rate fluctuations in the channels.

Further works are needed to improve the comparison between cavitating subsynchronous configuration and rotating stall phenomena.

To identify and analyze the mechanism responsible for the supersynchronous instabilities, we will adopt the nomenclature introduced in Fig. 11.

First, we will study the time evolution of the inlet and outlet flow rates in each channel during an inducer rotation period ͑Fig. 12͒, as well as the corresponding transient evolution of the total pressure variation ⌬P tot ͑Fig. 13͒. The difference between outlet and inlet flow rates in a channel ͑Fig. 12͒ represents the variation of volume of the vapor structure in the channel. These results concern the unsteady cavitation pattern illustrated in Fig. 8.

These results point out the following.

-The rise of the cavitation sheet size in Channel 2 ͑between Blades 1 and 2͒ leads to a partial obstruction of this channel, which induces the mass flow rate decrease. -As a consequence of the cavity length increase, we can observe the augmentation of total pressure variation in Channel 2. As a matter of fact, the cavitation sheet in Blade 1 does not reach the channel throat ͑i.e., the channel is not completely obstructed͒. The ⌬P tot ͑or blade load͒ increase is due to the rise of the angle of attack at Blade 1, which is associated with the flow rate decrease in Channel 2. -During the rotation of the small sheet from Channel 2 to Channel 1, the flow rates and ⌬P tot in Channels 3 and 4 are slightly modified.

Contrarily to the subsynchronous configuration, the propagation of the supersynchronous instability seems not to be related only to the inlet flow velocity fields. Indeed, the instability moves in the opposite direction as the velocity fields, as schematized in Fig. 14.

Hence, we think that the mechanism responsible for this kind of instability should be related mainly to the pressure field in the blade cascade. To study pressure fields, Fig. 15 illustrates the time evolution of the blade loads during the propagation of supersynchronous instability for a complete inducer rotation period ͑corresponding to Fig. 8͒.

In agreement with the flow rates and ⌬P tot time evolution analyses, we note that the loads of Blades 2 and 3 are slightly modified. Concerning Blade 2 ͑Fig. 15͑a͒͒ the pressure along the blade pressure side decreases during the development of the cavi- The loads concerning Blades 1 ͑Fig. 15͑c͒͒ and 4 ͑Fig. 15͑d͒͒ are strongly modified by the cavitation instability. The increase of the vapor sheet length at the suction side of Blade 1 is associated with a strong pressure rise, on the pressure side, near to the leading edge. This peak pressure could be the mechanism responsible for the decrease of the sheet in Channel 1.

To summarize, from the analyses of the whole of these unsteady results, we can propose the following scenario in order to explain the cavitation supersynchronous configuration ͑Figs. The Laboratory LEGI, in collaboration with the French Space Agency CNES and the rocket engine division of the Snecma, develops, for many years, numerical codes to simulate and analyze cavitation phenomena, mainly occurring in turbopump inducer geometries. In the scope of collaboration with NUMECA International, a cavitation model was implemented in the commercial code FINE/TURBO™. It has been successfully applied to perform calculations of 3D steady cavitating flows ͓18͔, as well as unsteady 2D simulations of cavitating flows ͓16͔.

In this paper, we presented first results obtained by 2D unsteady calculation in a four-blade cascade drawn from one inducer geometry. Computations were performed at the nominal flow rate for several cavitation numbers. The resulting performance chart was presented. A good agreement was found with previous work performed for the same geometry, by using another numerical code based on a different numerical scheme ͓13͔. characterized by four identical cavitation sheets -stable alternate blade cavitation ͑for ϳ 0.75͒ -supersynchronous configuration ͑also for ϳ 0.75 by in-termittency͒ -unstable alternate cavitation ͑for ϳ 0.7͒ -subsynchronous configuration ͑for ϳ 0.65͒ Based on the qualitative analyses of numerical results ͑i.e., density, velocity, pressure fields, and time evolution of the mass flow rates and of the blade loads͒, the mechanisms responsible for cavitation instabilities were analyzed and described. An original scenario concerning supersynchronous instability origin was more particularly proposed and discussed.

The duration of the calculations presented in this paper does not allow us to perform reliable quantitative evaluation of the characteristic frequencies. In previous work ͓13͔, by using the same physical model and a different computational method, many timeconsuming calculations have been carried out to evaluate the characteristic frequencies of instabilities in this same geometry. The main contribution of the present paper, in complement of the previous one, is to propose physical scenarios to describe cavitation instabilities. This study was also a test case used to evaluate the capability of the new 3D numerical code to perform unsteady calculations and to simulate cavitation instability. This step of the work is very important to prepare 3D unsteady calculation in the entire inducer geometries.

Further work is now needed to assess the prediction capability of the model. Presently, a detailed study is being performed to improve the applied physical models by taking into account thermodynamic effects ͓32͔. Numerical model aspects are also in progress to succeed 3D unsteady cavitating calculations in inducer real geometries ͓33͔.

This research was supported by the French Space Agency CNES ͑Centre National d'Etudes Spatiales͒ and the Snecma company. The authors wish also to express their gratitude to NU-MECA International for its cooperation to the development of the numerical code.

  l l o o o o o l o o o o . 3 I 3 o o l l

  side of the blade leads to the decrease of the applied mechanical torque, and consequently to the performance drop.

Figure 8

 8 Figure 8 illustrates this kind of instability. It is characterized by a nonsymmetrical rotating pattern of vapor sheets: We observe four different sizes of cavitation sheets in the four channels.Calculations indicate that this cavitation pattern rotates in the blade cascade in the same direction as the inducer rotation, but more quickly.

  l o , o o o o o , o o l o l o o o

  V l o , o o o o o , o o l o l o o o

  Channel 2. The load of Blade 3 ͑Fig. 15͑b͒͒ is slightly modified at the blade suction side ͑probably due to the decrease of the cavitation sheet in Channel 1͒.

  8 and 16͒. -The growth of the cavitation sheet on the Blade 1 suction side leads to the increase of the angle of attack on this blade. -A consequence of this phenomenon is the rise of the pressure at the Blade 1 pressure side, near the leading edge. -The pressure field in Channel 1, located below Blade 1, is modified by this pressure rise, which leads to the decrease of the cavitation sheet attached to Blade 4. -The presence of the cavitation sheet in Channel 4 will change the Blade 4 angle of attack, inducing the progressive rise of the attached sheet cavitation, and the continuation of the instability propagation.

  configurations ͑for Ͻ 0.65 or Ͼ0.8͒

c

  min ϭ minimum sound speed in the medium ͑ms -1 ͒ P ϭ local static pressure ͑Pa͒ P ref ϭ reference pressure ͑Pa͒ P ! ϭ vapor pressure ͑Pa͒ R c ϭ cut radius for the passage from three dimensions to two dimensions ͑m͒ T 0 ϭ inducer rotation period ͑s͒T ref ϭ reference time= blade passage time= T 0 / 4 ͑s͒ V ref ϭ reference velocity ͑=R c ⍀͒͑ ms -1 ͒ y+ ϭ nondimensional distance to the boundary ␣ ϭ void ratio ⍀ ϭ angular rotation speed of the inducer ͑rd s -1 ͒ ϭ mixture density ͑kg/ m 3 ͒ l ͑= ref ͒, ! ϭ liquid ͑=ref͒, vapor density ͑kg/ m 3 ͒ ϭ ͑P ref -P ! ͒ / ͑V ref 2 / 2͒ cavitation number ϭ static pressure coefficient= ͑P -P ref ͒ / ͑ L V ref 2 / 2͒ ⌬P tot ϭ totalpressure variation in a channel; ⌬P tot = P downstream -P upstream ͑Pa͒