Numerical simulation of cavitating flow in 2D and 3D inducer geometries

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


INTRODUCTION
In rocket engine turbopumps, disturbances due to cavitation phenomena can lead to substantial performance losses and=or strong unsteady forces acting on the pump components. The understanding of the unsteady behaviour of cavitation is therefore of primary importance for design purposes and has motivated the present development.
The cavitating ows observed in spatial turbopump inducers are characterized by complex, three-dimensional, two-phase structures, which are often localized close to the blades, as for example in the case of inducer tip cavitation. They are usually composed of a liquid=vapour mixture, whose structure can only be described macroscopically by dening a void ratio. Moreover, the cavitation process is often unsteady. Experimental results point out two main types of cavitation instability: • a self-oscillation behaviour of cavitation sheets, whose mechanism was studied in cavitation tunnels and analysed in detail by many authors [1][2][3][4]. • a rotating cavitation pattern, showing dierent sizes of the cavitation structures in the dierent blade-to-blade passages of the machine and leading to super-or sub-synchronous instabilities [5].
The diculty of predicting inducer cavitating behaviour is aggravated by these complex properties of the ow eld: three-dimensional, rotating, turbulent, unsteady, weakly compressible in the single-phase areas, highly compressible in two-phase ow conditions, and characterized by continuous phase changes. Previous numerical and theoretical models of cavitating ows in pump inducers proposed for example by Kueny et al. and Von Kaenel et al. [6,7] are based on the numerical simulation of steady cavitation sheets, attached to the blades. The cavities are usually considered as pure vapour, and the ow eld inside of them is not modelled. These models predict only the mean shape of sheet cavities on the blades. Physical cavitation models proposed by Delannoy [8] and Kubota et al. [9] are more eective for inducer applications, since they allow a complete 3D unsteady description of the cavitating conditions in the machine: vapour is created in all low pressure areas. In the rst case local uid density is controlled by a barotropic state law, while in the second case it is governed by the evolution of bubbles in the pressure gradient, on the basis of the Rayleigh-Plesset equation [10]. Three-dimensional Navier-Stokes codes based on these physical descriptions have been developed recently by Takasugi et al. [11], Alajbegovic et al. [12], Kunz et al. [13], Bunnel and Heister [14], and rst applications to pump geometries have been presented by Combes and Archer [15] and Medvitz et al. [16].
To gain further knowledge concerning cavitating ows in inducers, a three-dimensional model of steady and unsteady cavitation for inviscid or viscous uids is developed in the LEGI laboratory. That work is performed in co-operation with the Rocket Engine Division of SNECMA Moteurs and the French Space Agency CNES, with the nal objective to provide accurate simulations of unsteady cavitating ows in the inducers of rocket engine turbopumps. The numerical model results from the integration in the 3D code FineTurbo developed by Numeca Int. [17][18][19] of a physical cavitation model. This one is based on previous work performed in LEGI, concerning 2D numerical simulations of unsteady cavitating ows [3,8,20,21]. The two-phase aspects of cavitation are treated by introducing a barotropic state law that strongly links the uid density to the pressure variations. The liquid-vapour mixture is considered as a single uid in which the density varies from the liquid one to the vapour one, with respect to the local static pressure. The main numerical challenge of this approach results from the diculty to manage both an almost incompressible state in the pure vapour or pure liquid phases, and a highly compressible state in the liquid/vapour transition zone [22][23][24][25].
The present paper focuses on the computational method applied to treat eciently this twophase ow simulation. The use of the barotropic state law in association with a time-marching algorithm is not straightforward: rapid space and time density uctuations deteriorate the numerical stability, and the standard low-speed preconditioner is only devoted to incompressible or low-compressible ow. Although the preconditioner is here necessary to treat pure liquid and pure vapour congurations, its eciency is thus very poor in the two-phase compressible regions. A complete stability analysis was performed in Reference [26] to improve the understanding of these problems. The main results of this study are reported here, and solutions ensuring as well stability as satisfactory convergence rate and speed are presented in details. They imply local modications of the preconditioner in the vapour=liquid mixture, associated with a control of the density evolution at each step of the numerical resolution.
The presented computational method is validated for the prediction of the quasi-static behaviour of turbopump inducers. Two test cases are successively considered, respectively a two-dimensional Venturi type section, and a complete three-dimensional rotating inducer. The rst one was designed so that the ow is subjected to the same pressure eld as in an inducer blade-to-blade channel. It is thus representative of the cavitating conditions occurring in a real inducer, although its geometry is considerably simpler and phenomena associated with the machine rotation are not present. The geometry considered is characterized by a 4 • divergence angle, leading to a rather stable behaviour that has been already studied experimentally and numerically in LEGI [3,[27][28][29]. The second test case is a spatial turbopump inducer, whose experimental behaviour was investigated previously in the CREMHyG laboratory.
Time accurate simulations are performed on these two geometries. The aim in the rst case is to test the numerical model in a simple conguration, to check the accuracy of the results and the inuence of the numerical parameters. The second test case is a rst attempt to predict the cavitating ow eld in a complete 3D geometry, and compare it to experimental results. This validation of the code concerns presently only quasi-static two-phase ow conditions: pronounced unsteady behaviour associated with cavitation is out of the scope of this study. Current work is in progress to simulate eciently such congurations.
The main interest of the presented numerical model, compared to the previous ones, is its capability to predict the whole cavitating structure (including for example tip cavitation) that appear in an inducer in operating conditions, instead of only giving the extension of the attached cavities on the blades.
The paper presents the physical approach in Section 1, and the main features of the numerical method in Section 2, with emphasis on the special treatment required by the two-phase ow model. Section 3 is devoted to the modellization of the Venturi type section behaviour. Details of its geometry are presented in Section 3.1. A large range of cavitation number was investigated, in order to evaluate the model eciency. The inuence of the numerical and physical parameters was studied, and the conclusions are proposed in Section 3.2. Quantitative comparisons with experimental measurements are also analysed in Section 3.3.
The three-dimensional rotating inducer calculation is presented in Section 4. The quasi-static behaviour of the cavitating ow eld is investigated in a large range of NPSH (net pressure suction head). Starting from initial non-cavitating operating conditions, the NPSH is decreased slowly, down to highly cavitating conditions leading to the machine performance breakdown. Head drop charts obtained by the calculation are compared to the experimental ones, and the extension of the two-phase structures in the inducer are compared to visualizations performed previously at the CREMHyG laboratory.

THE PHYSICAL MODEL
Generally, two-phase ow models are based on the assumption that the uid is present in the computational domain both as liquid and vapour. The vapour is characterized by a density v , and the liquid by a density l . On each cell of the mesh, the unknowns are calculated for each phase, by averaging them on the volume occupied, respectively, by liquid and gas [30].
Neglecting the thermal eects, the number of balance equations in 3D is eight because of the two phases. These equations govern the behaviour of the two-phase structures larger than the cells, whereas the smaller structures are modelled by closure laws, which calculate empirically the uxes of mass and momentum between the two phases. The diculty of this kind of approach is to evaluate the transfer terms of mass and momentum in the balance equations.
In the present work, we apply a single uid model based on previous numerical and physical work developed in LEGI [8,20]. It assumes that only one uid is considered. This uid is characterized by a density that varies in the computational domain according to a state law ( Figure 1). When the density in a cell equals the liquid one l , the whole cell is occupied by liquid, and if it equals the vapour one v , the cell is full of vapour. Between these two extreme values, the cell is occupied by a water=vapour mixture that we still consider as one single uid. The void fraction =( − l )=( v − l ) can be dened as the local ratio of vapour contained in this mixture. If the cell is full of vapour, then = 1. If a cell is totally occupied by liquid, = 0. In this simple model, the void ratio is related to the state law, the uxes between the phases are treated implicitly, and no supplementary assumptions are required. Indeed, concerning the momentum uxes, our model assumes that locally (in each cell), velocities are the same for liquid and for vapour: in the mixture regions gas structures are supposed to be perfectly carried along by the main ow (the friction forces are high, compared to the buoyancy forces). This hypothesis is often assumed for this problem of sheet-cavity ows, in which the interface is considered to be in dynamic equilibrium [13,22]. The momentum transfer between the phases is thus strongly linked to the mass transfer.
Mass uxes resulting from vapourization and condensation processes are governed by the state law. In the present work, an empirical barotropic law (P) that links the density to the local static pressure is used. When the pressure is, respectively, higher or lower than the vapour pressure, the uid is supposed to be purely liquid or purely vapour, according to respectively the Tait equation [9] and to the perfect gas law. In this last case, temperature T 0 is imposed, since thermal eects associated with vapourization and condensation are not considered in the present study. = P r pg T 0 (perfect gas law for the pure vapour) (1) where P ref and ref are reference pressure and density, and r pg is the constant for the perfect gas. For water, P 0 =3× 10 8 and n =7. The two uid states are joined smoothly in the vapour-pressure neighbourhood by the following law: with A b and B b constants depending from the connection between the three parts of the state law, at the points (P l ; l ), and (P v ; v ).
It results in the evolution law presented in Figure 1, characterized mainly by its maximum slope 1=C 2 min , where C 2 min =(@P=@) min . C min can thus be interpreted as the minimum speed of sound in the two-phase mixture.

THE COMPUTATIONAL METHOD
The numerical model used to simulate cavitating ow elds is based on the 3D commercial code FineTurbo developed by Numeca Int. The cavitation model presented in the previous section was integrated into the algorithm, and some specic numerical treatments were applied to make this integration fully ecient. FineTurbo is a three-dimensional structured mesh code that solves the time dependent Reynolds-averaged Navier-Stokes equations, and whose a detailed description can be found in References [17][18][19]. We present here the main features of the numerical resolution, with emphasis on the modications introduced in the code to treat cavitating ow elds.

Steady-state calculations
The solution strategy of steady state problems is based on a time marching algorithm [31], with explicit Runge-Kutta time stepping. In the proposed cavitation model, thermal eects are neglected in the vapourization and condensation phenomena (see Section 1), so temperature does not appear in the state law of the uid. As a consequence, the energy equation is disconnected from the others: the temperature eld has no inuence on the resolution of the mass and momentum equations. Its resolution is thus of no use, and it is omitted hereafter. So the resulting governing equations, written in an integral conservative form in a control volume V whose surface is S, are: where U =(; u; v; w); F inv and F v are the inviscid and viscous uxes across the frontier S of V (normal n) Sce is the source term, and is the numerical time-step. Convergence is reached when the time derivative terms, which have no physical meaning, vanish.

Unsteady calculations
The time-accurate resolution procedure uses the dual time stepping approach (see References [31,32]). Compared to the steady-state formulation, the physical time derivative terms are added and the governing equations become At each physical time-step, a steady state problem is solved in pseudo-time and convergence is enhanced by a Multigrid strategy [19] coupled with local time stepping and implicit residual smoothing [19]. The above formulation is devoted to highly compressible ows. In the case of low-compressible or incompressible ows, its eciency decreases dramatically. As a matter of fact, at low Mach number the convective speed v of the ow becomes small, compared with the acoustic speed c. In other words, the eigenvalues of the system become much dierent, and the local time-step, which is based on the highest of these eigenvalues, is not adapted to the lowest waves, whose speed is v. This can be seen in the expression for the eigenvalues, expressed hereafter for a propagation direction n: As a result, the time-marching algorithms converge very slowly in such conguration.

Low-speed preconditioner
This well-known problem has been addressed by many authors and solved by introducing a preconditioner [33][34][35][36][37] in the Navier-Stokes equations. It is based on the modication of the pseudo-time derivative terms in the governing equations. Such modications have no inuence on the converged result, since these terms are of no physical meaning, and converge to zero. The resulting preconditioned system is controlled by pseudo-acoustic eigenvalues much closer to the advective speed, reducing the eigenvalue stiness and enhancing the convergence.
The simplest development consists in adding in the continuity equation a term (1= 2 p )@P=@. This solution was initially proposed by Chorin [38] to turn the system into an hyperbolic one, so that compressible time-marching algorithms can be applied. It was more recently improved by Choi and Merkle [33], who tuned the optimal value for p , and by Turkel [34][35][36] who proposed for a two-dimensional inviscid problem a supplementary term p v= 2 p @P=@ in the momentum equation: The term (1=c 2 ) u @P @x + v @P @y is derived from u @ @x + v @ @y , considering the particular case of perfect gases. Turkel analysed the inuence of the two preconditioning parameters p and p in this case, in order to remove the stiness of the system. Further investigations were presented by Choi and Merkle [33] to extend the capability of his preconditioner to treat viscous applications.
In the present work, the preconditioner of Hakimi [17] is applied. It is based on these two previous studies, and it consists in multiplying the pseudo-time derivative terms by a matrix −1 called the preconditioning matrix. The resulting governing equations are solved in the following conservative form: where Q(P g ;u;v;w) are the new dependent variables, P g is the gauge pressure (denoted simply P hereafter), F inv and F v are the same inviscid and viscous uxes as previously Accuracy is highly improved by choosing the gauge pressure as dependent variable. Preconditioning coecients are, respectively, equal to p = 1 and p = 0 · V ref where 0 is a constant whose default value is 3 and V ref is the reference velocity. As a result of this modication, the acoustic speed is replaced by a new one, which is of the same order of magnitude as u. So the stiness of the system is removed, and the convergence rate is considerably improved in the low speed ow areas [17,35].
However, the present application of cavitating ows implies the simultaneous treatment of very dierent ow congurations: according to the barotropic state law presented in Section 1, the uid is highly compressible in the two-phase ow areas (the Mach number can be as high as 4 or 5) and it is almost incompressible in the pure vapour or pure liquid regions. These two states of the uid have to be computed eciently by the same algorithm, without creating any spurious discontinuity in the ow eld. So a supplementary attention was paid to the preconditioning, which is unfavourable in compressible ow conguration, while it is necessary in weakly compressible or incompressible conditions. A modication of the algorithm is proposed hereafter to overcome this limitation.

Eigenvalues
The eigenvalues of the preconditioned system dened by Equation (9) are determined from the Euler part and can be obtained easily if the equations are expressed in terms of the primitive variables Q =(P; u; v; w). Suppressing source terms and viscous terms for sake of simplicity, a non-conservative form of Equation (9) is given by Since the density only depends on pressure, the sound celerity is dened as c = dP=d. As a consequence, @=@x, @=@y can be replaced respectively by 1=c 2 @P=@x,1 =c 2 @P=@y, and the mass equation becomes The system can be turned into the following form, involving the non-conservative Jacobians A; B; C: The nal expression for A; B; C depends on the expression for c, issued from the dierent parts of the barotropic state law (P).
As mentioned in the previous section, the preconditioner may be unfavourable in the compressible ow areas, i.e. in the two-phase mixture. So we focus hereafter on this particular part of the barotropic state law. In this conguration, the expression for the eigenvalues associated with propagation direction n(n x ;n y ;n z ) can be obtained from |An x + Bn y + Cn z − l| =0: and It was clearly shown in Reference [26] that the stiness of these four eigenvalues, which mainly depends on the term 2 p M 2 C2=B b u 2 , is much increased in such a two-phase ow conguration, compared with weakly compressible ow conditions. Thus, a modication of the preconditioner is suitable in compressible areas, to remove locally its action. Preconditioning parameters must be adapted progressively, to nally reach the following values in the most compressible ow parts (corresponding to =0:5): (i.e. p is the physical sound celerity in highly compressible regions of the ow: the action of the preconditioner is in this case locally completely removed).
To ensure a smooth evolution of parameters p and p , the following expressions are nally applied in the whole ow eld: The expression for p ensures for all ow conditions the same order of magnitude for the three eigenvalues 1 , 2 , and 3 . In weakly compressible or incompressible congurations, p is maintained at the constant value 0 · V ref , which limits the weight of the term 4 2 p · n 2 in 2 and 3 , while in highly compressible regions, p is reduced down to the physical sound celerity, to control the terms 2 p M 2 C2=B b u 2 . This improvement is applied in all the two-dimensional calculations reported in the present paper. Concerning the 3D calculations of rotating machines, work is in progress to enhance the preconditioner modication and avoid spurious computation divergences.

Discretization
Equation (9) is discretized in its conservative form with a nite volume approach. The discrete form of Equation (9) over a computational cell volume becomes where (F · n) * is the numerical ux at the cell interfaces. A centred approximation is applied to the viscous uxes, while the inviscid ones are calculated with a central convection scheme associated with articial dissipation. The resulting expression for the numerical ux along the i-direction on the right side of the cell is The rst right hand term is the centred approximation of the inviscid ux, and it is treated as The articial dissipation D art is composed of two terms, respectively of second and fourth order, as initially proposed by Jameson et al. [32]. The following form of D art leads to a central second order accurate convection scheme: where i is the spectral radius along the i-direction multiplied with the cell face area. d 2 and d 4 are coecients associated, respectively, to second order and fourth order articial dissipations. Their expression is given by The variables i and i are sensors to activate the second-dierence dissipation in regions of strong gradients, and to de-activate it elsewhere. i (proposed initially by Jameson et al. [32]) applies in pressure gradients areas, such as shocks, and i is a supplementary sensor added for cavitating ows: The standard values for 2 p and 2 are, respectively, 0.5 and 1. The second term has been added to ensure the numerical stability in cavitating areas, without modifying the numerical scheme for other ow congurations. The standard value for 4 p is 0.1. The inuence of these numerical parameters on the result of calculations in cavitating conditions is discussed in Section 3.
Note that all expressions above correspond to the numerical ux in the x direction. Similar treatments are applied for the two other directions.
The pseudo-time integration is made by a four steps Runge-Kutta procedure [31]. The physical time-derivative terms are discretized with a second order backward dierence scheme that ensures a second order accuracy in time:

Stability analysis of the numerical scheme
The numerical stability of the resolution has been studied in the three ow congurations dened above. It is based on the Von Neumann stability analysis described by Hirsch [31], and on the Fourier footprint representations used previously by Hakimi [17]. Details of this investigation are presented in Reference [26], and only the main results are given here. For sake of simplicity, the study is performed in two dimensions, and pure vapour as well as pure liquid are considered fully incompressible. The inuence of several physical and numerical parameters was tested, and the eciency of the preconditioner was analysed in several ow congurations. Nice convergence and stability properties are demonstrated in the incompressible case (which conrms the eciency of the preconditioner in such congurations), while poor characteristics are obtained concerning the central part of the barotropic law. It is also found that the maximum slope of the barotropic law (controlled by the value of C min ) strongly inuences the convergence rate and the stability: the reduction of C min induces a severe decrease of both convergence speed and numerical scheme stability. The other main physical parameter, namely the orientation of ow velocity inside the cells, has the same kind of impact: the conguration v=u = 1 is unfavourable for convergence and stability, compared with conguration v=u = 0 (ow in the direction of the mesh).
Concerning the numerical parameters, a diminution of the fourth order articial dissipation 4 p induces an important convergence slow down. So does also an increase of the mesh aspect ratio.

Under-relaxation of the density
The cavitation phenomenon is a very sharp and rapid physical process, which results in two main numerical diculties: • Vapourization and condensation phenomena induce strong density gradients. The interface between pure liquid and pure vapour can thus be very thin. It may lead to successive cells with respectively pure vapour and pure liquid inside, even in the case of very ne meshes. To ensure the local conservation in this conguration is complicated, and such discontinuity is highly unfavourable for convergence rate and scheme stability. • The two-phase structure of cavitation is fundamentally unsteady, particularly in the closure area of the cavitation sheet. The collapse of vapour structures happens so suddenly that particles of uid can turn from pure vapour to pure liquid during only one time step. A smooth description of this phenomenon would impose a physical time step about hundred times smaller than the one used currently. Considering the frequency of typical unsteady processes associated with cavitation. such as cavity self-oscillation, it appears that the use of a very small time step is unreasonable, because of the total computational time it would induce. This particularity leads to severe numerical instabilities.
For these two reasons, an ecient control of the density time uctuations and space distribution must be applied. Density variations are thus under-relaxed at each iteration i to prevent too sudden changes in a single pseudo time step: where the standard value of under-relaxation coecient is 0.2. This value has to be lowered down to 0.05 in the case of 3D calculations in turbopump inducers (Section 4). After such treatment on the density, pressure P i+1 is updated according to the barotropic relation (P), so that the uid state law is always respected. Underrelaxing the density is the same as multiplying the local time step by in the mass equation. As a matter of fact, the convergence signicantly slows down when underrelaxation is activated. Nevertheless, the computational time increases much more if a global reduction of the pseudo time step is applied also in the momentum equations. So underrelaxing the density was found to be a satisfactory choice to improve stability without being confronted to prohibitive computation times.

Boundary conditions
The boundary condition setting is based on a system of dummy cells. Classical incompressible type of boundary is applied: imposed velocities at the inlet, and an imposed static pressure at the outlet. Numerical studies have been performed in LEGI to improve these conditions, mainly by taking into account the test rig inuence [39]. They are not yet applied in the present application.

Initial transient treatment
In experiments, the classical machinery cavitation test consists in setting initially a relatively high pressure P ref in the ow eld, for which no vapour appears. Then, the pressure is decreased until the cavitation number =(P ref = P vap )=V 2 ref corresponding to operating conditions is reached.
A similar numerical procedure was developed, but instead of decreasing the reference pressure P ref in the ow eld, this one is kept constant, and a pseudo vapour pressure P vap is increased. First of all, a steady step is carried out, with a pseudo vapour pressure low enough to ensure no vapour presence in the whole computational domain. Then, this vapour pressure is increased smoothly at each new time step, up to the physical value. Vapour appears progressively during this rising of the vapour pressure. The cavitation number is then kept constant throughout the computation.
Another method would consist in decreasing the outlet pressure, with a constant vapour pressure. The drawback is however that the physical parameters in the whole ow eld change because of this continuous boundary condition variation. With the rst method, only the twophase areas are aected by the cavitation number evolution.

Turbulence model
Two turbulence models are applied for the simulations presented in this paper: • a standard Baldwin-Lomax model is applied in most of the computations. Details about this model can be found in Reference [40] • the numerical tests presented in Section 3 in the case of the Venturi type section include the use of a k-model with extended wall functions.
The Baldwin-Lomax model requires at the walls a rst cell located at a non-dimensional distance y+ close to 1, while in the case of the k-model, extended wall functions enable the use of various values of y + at the walls. So the same grids are used for both models of turbulence. More calculations considering k-turbulence models are in progress to improve physical analyses. An investigation of the respective eciency of dierent turbulence models associated with the single uid modellization of cavitation is also proposed by Coutier-Delgosha et al. in Reference [41].

mesh
Numerical simulations are rst performed on a two-dimensional Venturi type section, characterized by small convergent=divergent angles (4:3 • =4 • ) and a small contraction ratio at the throat. The shape of the Venturi bottom downstream of the throat simulates an inducer blade suction side with beveled leading edge geometry (see Reference [42]). According to experimental observations performed previously in the LEGI laboratory [28,29], a stable sheet of cavitation is obtained, with only small-scale uctuations in its downstream part. Although the considered geometry is much simpler than a real inducer one (3D eects and rotating phenomena are not present), its cavitating behaviour is globally representative of sheet cavitation phenomena observed in such a machine.
The Venturi prole was tested previously in LEGI [28,29], and the experimental data are used in the present work to validate the numerical simulations.
Three grids of dierent size (respectively, 108×36, 160×60, and 248×88 cells) are tested hereafter (Figure 2). A special grid contraction is applied in the areas of cavitation, i.e. after the Venturi throat and especially near the leading edge. A supplementary contraction is also imposed close to the foil and along the tunnel walls, to be consistent with the y + values required for the use of the turbulence models. It results in practice in y + values comprised between 0.5 and 5 at the centre of the rst cell at the walls.

Validation of the global cavitating behaviour
Steady-state calculations have rst been performed, and a satisfactory convergence level was obtained in the case of small cavities. For a lower cavitation number, some important numerical diculties appear, associated with large pressure uctuations in the Venturi section. As a matter of fact, cavitation is always unsteady at some scale, and the uctuations observed in the rear part of the cavity in experiments are very prejudicial to the convergence of steady-state computations. Time-accurate computations have thus been carried out, and a much better convergence rate was obtained for a large range of cavitation numbers. Comparisons with previous steady-state results showed a close agreement, so far the convergence levels were identical. Consequently, all the results presented hereafter were obtained with time-accurate computations.
The stable cavitating behaviour is correctly simulated by the model: in all the analysed cases, we obtain, after a transient uctuation of the cavity length, a quasi-steady behaviour of the cavitation sheet, which globally stabilizes. Figure 3 presents the shape of the converged cavity in the reference simulation. It corresponds to a cavitation number =0:4. The inlet velocity equals 10:8 m=s. We obtain in this case a stabilized cavity whose length is about 40 mm (that is to say L cav =L ref =0:16). 18 s of simulation were performed, that is to say 600 T ref , where T ref is a reference time, corresponding to the time necessary for the ow eld to cover the length L ref of the prole, with a speed V ref . The outlet pressure is set arbitrarily to 85 000 Pa, which leads to the pressure distribution drawn in Figure 3.
The transient evolution observed during this unsteady calculation is shown in Figure 4. In Figure 4(a), the colours represent the density values: white for the highest ones (pure liquid) and darker (red and then black) as the void ratio increases. The time evolution is represented in abscissa, and the position in the tunnel of cavitation is graduated in ordinate. So an horizontal line corresponds to the time evolution of the density in the corresponding section of the cavitation tunnel. Figure 4(a) illustrates at a given time and for each section the value of the minimal density present in the section. In this stable conguration the main information provided by this gure is the time evolution of the cavity length, associated with the maximum void ratio in each section. We also analysed the evolution of the vapour volume ( Figure 4(b)) and the inlet pressure uctuation (Figure 4(c)).
The length of the cavitation sheet (40 mm) is a little under estimated, since the experimental value for this cavitation number is 55 mm. Nevertheless, it can be observed in Figure 5, that a general good agreement is obtained with experiments concerning the cavity extension. Two experimental results are drawn: they are related to two values of the reference velocity V ref (9 m=s and 14:4 m=s) close to the value V ref =10:8 m=s applied for the calculation.
The length of the cavitation sheet is correctly predicted, down to =0:4. As the cavitation number continues to decrease, the numerical model leads to a cavity length smaller than the values reported by the experiments. But we must keep in mind that the experimental estimation of the cavity extension cannot be precise in this case, because of the uctuating rear part of the liquid=vapour structures.

Numerical tests
The inuence of numerical and physical parameters of the model is tested in the previous ow conguration. The objective is to determine the eect of their variation on the convergence rate and on the result. Table I presents the investigated parameters, with the corresponding range, and the standard value applied for the simulations. Inuence of numerical parameters indicates the accuracy level of the numerical scheme, while inuence of the physical ones gives an evaluation of the code robustness in a cavitating conguration.
Results of the tests are reported in Table II. As the nal cavity obtained is almost stable, it can be characterized by its nal volume, and the vapour volume contained inside. So for each case, these two results are indicated. The nal volume of the cavity is calculated by considering the volume of the cells where ¿4%.
3.3.1. Inuence of numerical parameters. The inuence of the numerical parameters, such as t; 2 ; 4 P ;, CFL, and the mesh size is globally satisfactory: • A too high time-step t generates some little oscillations of the cavity. For lower values of t, its inuence is acceptable: for t =10 −4 s, the cavity volume is by 7% smaller

Inuence of the physical parameters.
An increase of C min or v = l leads to an augmentation of the cavity thickness. As a matter of fact, these two parameters control the shape of the barotropic law. The more they increase, the more the maximum slope is reduced, resulting in a smoothing of the density variations. The general trend associated with this evolution is an increase of the diusive eects in the two-phase ow areas. In the case of v = l this inuence is small. On the contrary, the value of C min directly inuences the cavity shape. The calibration of this parameter is performed on the basis of 2D previous studies performed in LEGI [10]. Its reference value is chosen so that the thickness of the cavity is consistent with experimental results, for several dierent cavity lengths. The optimal value was found to be independent of both L ref and V ref , and is about 2 m=s.
• The turbulence model has only a little inuence on the result, in this rather stable cavitation conguration. • Numerical stability is ensured down to very low values of the cavitation number . For example, in the conguration =0:33, the sheet of cavitation is so large, that the blockage of the Venturi type section is almost reached (cf. Figure 6). Nevertheless, the computation goes on satisfactorily, in spite of important pressure uctuations in the channel. This behaviour, although it is associated with a noticeable drop of the convergence rate, indicates that the numerical model is properly robust in highly cavitating conditions. It suggests that the computation of inducer performance breakdown, as blade-to-blade channels are almost blocked by the volume of vapour, is a realistic objective.
These tests validate the choice of the reference values for all the numerical parameters used here after.

Comparison with measurements inside the cavity
Comparisons are investigated with experimental data obtained by double optical probes. This technique and the results are presented in detail in References [27,28]. This is an intrusive captor, which allows measurements of the local void ratio and the velocities of the two-phase structures inside the sheet of cavitation. Five data proles are available: their position is indicated in Figure 7. The mean values of the void ratio are presented for these ve proles in Figure 8. The dotted line corresponds to the experimental external shape of the cavity, obtained from image processing.
A good agreement is obtained between experiments and the calculation concerning the void ratio distribution along the x-direction: the low density measured in the upstream part of the cavity, and its smooth increase in the rear part, are correctly predicted. Nevertheless, the cavity thickness calculated by the numerical model is clearly over-estimated, mainly along Sections 2-4. That discrepancy is related to the fact that the standard value for the parameter C min is slightly higher than the one obtained in previous calibration of the physical model (2:25 m=s in the present case, 1:5-2 m=s reported in References [3,20,41]).

CALCULATION ON A FOUR-BLADE INDUCER
Simultaneously with the previous 2D studies on xed geometries, we performed some rst computations on a rotating four-blade inducer designed at the SNECMA Moteurs Rocket Table II. Results of the tests: inuence of numerical parameters.   Engine Division (Figure 9). This machine was tested experimentally, and the visualizations show a pronounced unsteady behaviour: depending on the value of the cavitation number, we can observe alternating blade cavitation or super-synchronous regimes. Moreover, the two-phase structures are strongly uctuating, particularly in the cavitating vortex of the tip area.
In the present study, only the quasi-static behaviour of the runner is simulated. As for the previous 2D Venturi type section, a non-cavitating steady computation is rst performed. Then, the unsteady calculation is started, and the NPSH is slowly decreased at each physical time step, as in the experiments. Vapour structures spontaneously appear and grow during this process, in the regions of low static pressure. In the presented computations, the NPSH is decreased until the inducer head drop, due to highly cavitating conditions, is reached. The physical time-step is high enough to only simulate the quasic-static eects of the NPSH decrease. So the unsteady eects are not modelled in the computations presented hereafter. For a given pseudo-vapour pressure value, the inlet pressure does not uctuate so far the convergence is reached, and it is directly used to compute the NPSH values reported on the charts.

Mesh
Only one blade-to-blade channel is considered, and periodicity conditions are applied to its frontiers in azimutal directions to simulate the presence of the contiguous blades. It implies that the four channels are supposed to have the same behaviour, which is not penalizing so far only quasi-static eects of cavitation are investigated. The mesh is composed of ve 3D blocks ( Figure 10): the upstream part of the bulb, the blade to blade channel, the tip gap, and two blocks, respectively, upstream from the blade leading edge and downstream from its trailing edge. The whole mesh is composed of 420 000 cells.
The main diculty to mesh an inducer, compared with a typical centrifugal impeller, is related to the low b angles in the inlet and in the outlet. This kind of geometry can lead to a much distorted mesh, with triangular cells in the vicinity of the leading edge. Indeed, the computational domain inlet and outlet are at constant coordinate z (cf. Figure 10), so the mesh is highly rounded in the blade-to-blade representation: a rst time around the leading edge and, a second time in the other direction around the trailing edge. This distortion implies non-matching cells along the frontiers connected by periodicity conditions, which complicates the boundary conditions and induces some supplementary numerical errors due to the variable interpolations.
Blocks 4 and 5 are used presently to treat correctly the sharp leading and trailing edges. On the basis of the automatic mesh generation, special attention was then paid to the grid in the cavitating areas: an axial stretching is applied on the blade suction side, near the leading edge, and the aspect of the mesh just upstream from the blade is improved to reduce the aspect ratio of the cells. Despite these eorts, the poor description of the blade tip could not be completely removed. The grid in the tip gap is composed of eight cells in the radial and azimuthal directions, so that the ow from the pressure side to the suction side is correctly modelized. A contraction of the mesh is systematically applied close to the walls, to be consistent with the requirements of the Baldwin-Lomax turbulence model. Figure 11 presents some details of the mesh: two zooms around the leading edge, respectively near the hub and close to the shroud, and a 3D view of the grid on the hub and on the blade tip.

Boundary conditions
The calculations are performed with water ow conditions. The rotating velocity is 6000 rpm, and the Reynolds number based on the shroud diameter, the inlet tangential tip velocity, and the water viscosity, equals 10 7 . A standard Baldwin-Lomax turbulence model is applied. The boundary conditions are the following ( Figure 12): • An axial velocity is imposed at the inlet, on the basis of the nominal inducer ow rate. • A radial equilibrium is xed at the outlet for the static pressure.
• The hub and the blades are rotating, the shroud is immobile. The gap between the blade tip and the shroud is taken into account.

Numerical parameters
l = 1000 kg=m 3 and v = l =0:01 (standard value), CFL number: 1.5, =0:05, t =10 −4 s (about 1=100 e of the inducer rotation period, to only simulate the quasi-static eects), The p preconditioning coecient is based on the tangential inlet tip velocity, All other parameters are set to their default value indicated in Section 3.3.

Results: performance charts
Cavitation conditions in the inducer and its performance are usually characterized by the following non-dimensional numbers: V ref is in this conguration the tangential tip velocity, and P inlet lot , P inlet , and P outlet are, respectively, the total pressure upstream from the inducer, and the static pressure upstream and downstream from the inducer.
A decrease of the cavitation parameter ′ i was performed from a non-cavitating conguration to the inducer performance drop. The resulting evolution of the static pressure elevation of the inducer is compared in Figure 13, to experimental measurements performed in the CREMHyG laboratory. The static pressure elevation is obtained in experiments with pressure sensors located on the shroud, from inlet to outlet. Pressure values issued from the calculations are taken at the same points, so that the comparison presented in Figure 13, is relevant. The static pressure elevation obtained by the computation is in excellent agreement with the experimental results. Concerning the head drop location, the cavitation parameter ′ i corresponding to the performance breakdown is over-predicted by the model: the gap between experimental and numerical results equals ′ i =1. This discrepancy must be considered circumspectly. In the case of centrifugal pumps [43], the order of magnitude of ′ i is notably higher, so that the same absolute error results in a much lower relative error (less than 5% in most of the simulations [44]). Since an inducer has much better suction capacities ( ′ i at the performance breakdown is about 1), the order of magnitude of ′ i is in the present case the same as the one of ′ i . So, the over-estimation of ′ i 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.
It can be also noticed that pressure uctuations observed in experiments before the performance drop are not obtained in the calculation. This discrepancy is mainly due to the fact that a high time step is applied in the present computation, so that only the quasi-static effects associated with cavitation are modelized. As a matter of fact, these oscillations are often correlated with unsteady cavitation behaviour, such as rotating cavitation eects. This kind of phenomenon is not treated here, since only one blade-to-blade channel is computed.
The general good agreement obtained between the experimental and the numerical performance charts (also it is not ecient enough to predict the inducer head drop location) is promising. It allows to study more closely the development of vapour inside the blade-to-blade channels, and to compare it to experimental visualizations. Figure 14. Evolution of the cavitating areas in the inducer as ′ i , decreases (the hub is coloured in grey, the blades in red, and the iso-surfaces = 950 kg=m 3 are represented in yellow).

Visualizations of the cavitating areas
The evolution of the two-phase structures during the ′ i decrease is presented in Figure 14. The hub is represented in grey, the blades in dark red, and the surfaces corresponding to a 5% void ratio are coloured in yellow. The operating points related to the nine visualizations are indicated in Figure 13, by the green squares. Two cavitating structures appear: a cavitation sheet attached to the blade suction side, and a tip cavitating vortex. This one is rst conned on the blade suction side, close to the shroud (view 2). Then, it rapidly lls the gap, and it progressively joins the pressure side of the adjacent blade (view 5). This tip cavitation continues after that to grow downstream, while the suction side cavity is almost stable (view [6][7][8]. Two comparisons between experimental visualizations and numerical results are presented in Figure 15. They correspond to two operating points located at the beginning of the performance breakdown. Since the experimental and numerical head drop are slightly different, the cavitating structures in the inducer are not compared at the same value of ′ i , but at the same performance drop (respectively, 0.5 and 1.3%).
A reliable qualitative agreement between experiments and calculation can be observed on these two examples: both the attached cavity and the tip cavitating structure are correctly obtained in the two cases. The shape of the cavitation sheet predicted by the numerical model in the vicinity of the hub is particularly close to the one observed in experiments. The tip cavitation, which seems to be characterized experimentaly by a rather low void ratio, is also qualitatively well simulated by the model. It is not signicant that its precise shape is not the same in experiment and in the calculation, since this structure is highly uctuating in time and in space. Nevertheless, some discrepancies can be seen at mid-height of the channel, where the cavity predicted by the computation is notably shorter than in experiments.
So, the experiments and the calculation appear to be consistent for these two operating conditions. It suggests that the 3D numerical model applied is ecient to predict the main features of the quasi-static development of cavitation in a spatial turbopump inducer. The remaining discrepancies may explain the overestimation of ′ at the performance drop. Nevertheless, this rst attempt to predict the cavitating behaviour (including attached cavity and detached structures) of such a complex geometry is very promising. A continuing work is pursued to improve the simulations. An improvement of the mesh aspect ratio in the cavitating areas (blade tip, leading edge) is an important way of investigation.

CONCLUSION
A 3D model for unsteady cavitation was presented in this paper. It is composed of a time-marching algorithm associated with a barotropic state law that governs the uid density evolution in the cavittaing regions. Particularities of this computational method, such as a modication of the low speed preconditioner and a necessary density uctuations control, were detailed. The accuracy and robustness of the numerical model have been tested rst with a 2D Venturi type section representing an inducer blade-to-blade channel. The global shape of the cavitation sheet is correctly predicted for rather small cavities, and the model seems to slightly under-estimate the vapour/liquid structures extension for smaller cavitation numbers. Further investigations are performed to improve the modellization. A complete investigation of the inuence of physical and numerical parameters was also presented, and reference values were dened for further calculations. The structure of the cavity was nally compared to experimental results, and a satisfactory agreement was found concerning the external shape of the cavity and the void ratio distribution inside. A rst attempt of a 3D simulation concerning a four-blade turbopump inducer was also presented. The quasi-static results show a promising agreement with expeimental measurements and visualizations. Work is continuing to improve the prediction of inducer performance breakdown, and simulate pronounced unsteady phenomena associated with cavitation, such as cavity self-oscillation and rotating cavitation.