Numerical and Experimental Investigations of the Cavitating Behavior of an Inducer

A robust CFD model is described, suitable for general three-dimensional flows with extensive cavitation at large density ratios. The model utilizes a multiphase approach, based on volume-scalar-equations, a truncated RayleighPlesset equation for bubble dynamics, and specific numerical modifications (in a finite-volume solution approach) to promote robust solutions when cavitation is present. The model is implemented in the CFD software CFX TASCflow 2.12. The validation of the model was done on an inducer designed and tested at LEMFI. First, The physical model and the numerical aspects are described. Then, the experimental and numerical methodologies, at cavitating regime, are presented. Finally, for several flow rates, the comparisons between experimental and simulated results on the overall performances, head drop and cavitation figures, are discussed. For a range of flow rates, good agreement between experiment and prediction was found.

2002; Frobenius et al., 2002).These methods have been introduced only recently as a design tool and therefore there is yet very little experience in their use, especially in turbomachines (Dupont and Casartelli, 2002;Combes and Archer, 2000;Visser, 2001;Coutier-Delgosa et al., 2002).
For these approaches, the solution of general multidimensional cavitation problems via CFD has, in a very broad sense, fallen into two distinct categories: 1. Solutions based on single-phase conservation equations along with an equation-of-state (applied with the assumption of local thermodynamic equilibrium between the phases) to establish the two-phase mixture conditions (Reboud et al., 1998) 2. Solutions based on solving conservation equations for each of the phases (Athavale et al., 2002), with a variety of modelling assumptions applied in the development of non-equilibrium heat/mass/momentum exchange terms between the phases In either category, the emphasis is on solving the mass, momentum and sometimes energy transport, in a continuous fashion, over the non-cavitating and cavitating regions of the flow (Tremante et al., 2002).It should be mentioned here that interface tracking methods (Hirschi et al., 1998a) have not been included in the above two categories since, in general, they do not attempt to solve the full transport equations in both the cavitating and non-cavitating regions.As an example of the first approach, the solution of the mass and momentum equations, along with some assumption regarding the thermodynamic path of the flow (for example constant entropy or enthalpy) and a corresponding barotropic relationship (generally of the form P = f (ρ) where ρ = α l ρ l + α v ρ v ) for the mixture, allows the vapor content to be evaluated at any point in the flow.This approach provides a 16 F. BAKIR ET AL.
simple, and yet effective means to implement a cavitation model since no empirical constants are required other than equilibrium phase data for the substance (Avva et al., 1995;Saurel et al., 1999;Bunnell and Heister, 2000).
The second approach becomes useful as the equilibrium assumption breaks down for flows where the accelerations are large.In such flows the heat and mass transfer associated with any phase change will be slow relative to dynamic changes in the local pressure.From a modelling perspective, an equilibrium model will predict spatial gradients in the mixture density that are excessively large.Such situations make numerical solutions extremely difficult, particularly if the liquid/vapor density ratio is high (as is almost always the case for water).This situation is alleviated with the non-equilibrium models as rate limitations are placed on the vaporization (and condensation) process.For flows with low liquid/vapor density ratios (for example cryogenic fluids such as LOX), use of non-equilibrium models are also preferred since larger quantities of the liquid phase must be vaporized to support the volume growth of the vapor during cavitation.This increases the amount of energy to be removed from the surrounding liquid and results in larger thermal boundary layers.This process, termed the thermodynamic effect, is particularly important when operating at higher temperatures and can result in cavitation suppression (Billet et al., 1981).
The present model uses a non-equilibrium approach with simplifying assumptions to reduce the number of equations to solve in a multiphase system.The numerical adjustments, to solve the final multiphase system in a robust and accurate manner, will be presented around this framework.It is common practice in cavitation model development, to assume mechanical equilibrium between the phases and thereby reduce the number of momentum equations to a single set (Senocak and Shyy, 2001;Yuan and Schnerr, 2001;Medvitz et al., 2002;Ahuja et al., 2001).For non-equilibrium descriptions of the dynamics of phase change, volume-fraction equations are retained in order to solve for the distribution of the phases in the flow (Von Dirke et al., 1999).From this point, many approaches can be taken to arrive at the final form of the governing equations depending how the micro-scale phenomena is treated.Assumptions regarding the heat/mass transfer between the phases and the distribution and size of the vapour bubbles are required and, depending on the nature of the cavitation process, some approximations will be more appropriate than others.For this paper emphasising performance prediction, a cavitation model based on an assumed thermal equilibrium between the phases and a truncation of the Rayleigh-Plesset equation, governing bubble dynamics, has been developed.In this model, the vapour volume fraction in the flow is determined by a rate equation, introduced into the source terms of the VOF equations.In order to validate the VOF model, the numerical results were compared to the experimental results obtained on an inducer tested at LEMFI (Bakir et al., 1998).The comparison was done for both over-all performances and cavitating regime for a range of flow rates.

Volume Fraction Equations
The governing equations describing the cavitation process assumes a two-phase three-component system, with no interphase slip and thermal equilibrium between any of the components and phases.The three components in the flow are: 1. non-condensable gas (dispersed phase) 2. vapour (primary phase) 3. liquid (primary phase).
The relative quantity of each of the components is described by a volume fraction scalar α.Particular volume-fractions have the subscripts: d-dispersed, v-vapour in primary phase, and l-liquid in primary phase.The sum of all volume fractions must sum to one: where the volume fractions are related to the mass fractions, y, for each component through the relations: with the additional constraint that y d + y l + y v = 1.In the above relations ρ k is the density of the individual components (subscripts k = d, v, and l), or for the mixture when designated as simply ρ.For the general cavitation case, compressibility in the vapor and liquid regions are important and therefore an equation of state is required to determine the component density.Depending on the component considered, the equation of state required is: • The non-condensable component is assumed to be a gas (air) and its density can be determined from an ideal gas equation of state using local pressure and temperature information from the solution.
• The primary phase vapor component density can also be calculated using an ideal gas relation for the vapor, or determined on the basis of the saturated vapour curve for the primary fluid.For low speed problems, the vapor density can be assumed a constant evaluated at the vapor pressure and nominal fluid temperature.
• The liquid component of the primary phase, can be evaluated from a liquid EOS such as the Tait equation (Avva et al., 1995) if compressibility is important.In the present paper liquid compressibility is not significant and therefore a constant liquid density is applied.
Based on the above description of a multiphase system only, two volume-fraction equations need be solved since the distribution of the third phase can be determined from Equation (1).In many cavitation problems, the mass fraction associated with the noncondensable phase can be assumed to be well mixed in the liquid phase with a constant y d .On this basis, the mass fractions y l and y d can be combined and treated as one.This is the approach taken here and the volume scalar α m is introduced for this where α m = α l + α d and the density associated with α m becomes: Choosing the scalar α m to solve the governing volume-fraction equation for the primary liquid phase with non-condensable gas becomes: where The source terms have units of (kg/s) where the source terms Ṡv and Ṡl account for mass exchange between the vapour and liquid in the primary phases during cavitation.The form of the source term Ṡv for this paper, has been derived by considering the Rayleigh-Plesset equation for bubble dynamics which will be described in a subsequent section.

Mixture Mass, Momentum, and Turbulence
Since a basic assumption to the model is that all phases share the same velocity, the governing mixture equations for mass and momentum are: where the density is defined through volume fractions as ρ = α m ρ m + α v ρ v and the velocity field, u j , calculated from the momentum equations.Momentum conservation is described by: where τ i j is the deformation tensor and is further simplified by the assumption of Newtonian fluid behaviour.The flows in the devices presented in this paper are turbulent and require turbulence models to estimate the enhanced mixing effect on the flow.Most turbulence models have been developed for singlephase flows and for very simple geometries.Little work has been done to estimate the influence of phase change along a surface, in particular that associated with cavitating flows.Due to the limited experimental work to provide guidance, the standard two-equation high-Reynolds number k-ε model was employed for the liquid and cavitating regions without modification.The governing equations will not be presented here while the reader can refer to the paper of Launder and Spalding (1974) for further details.

Cavitation Source Term: Rayleigh-Plesset Model
A cavitation model has been implemented based on the use of the Rayleigh-Plesset Equation (RPE) to estimate the rate of vapor production.For a vapor bubble nucleated in a surrounding liquid, the dynamic growth of the bubble can be described by the RPE as follows: Where R is the radius of the bubble, p v the vapor pressure in the bubble, p the pressure in the surrounding liquid, and the liquid density ρ l .For the present model we assume that p v is at the vapor pressure corresponding to the temperature of the liquid.The above nonlinear ordinary differential equation is difficult to implement within an Eulerian-Eulerian framework for multi-phase flows, therefore a first order approximation is used where: In what follows we develop a source term for the generation and destruction (vaporization and condensation) of vapor bubbles.
The basic assumptions for the model are: • the growth or collapse of a bubble follows the RPE neglecting higher order terms • bubbles are assumed to grow from an initial average radius of R b and return (when condensing) to bubbles of the same size • there are no thermal barriers to the droplet growth.
and the number of bubbles per unit volume of the mixture, N , available as nucleation sites (during vaporization) is given by: and during condensation, Since the volume fraction α m is to be solved representing well mixed liquid and non-condensable gas, Equation ( 9) is applied as a liquid source term: Equation ( 13) is used for the vaporization and condensation of liquid where the definition of N (Equation [10] or [11]) changes depending on the direction of the phase change.In practice the vaporization and condensation processes have different time scales, the condensation process typically being the slower one.An empirical constant, F c , is introduced to allow for these constraints: with a typical default of F c = 50 when vaporization occurs ( p v − p > 0) and F c = .01with condensation ( p v − p < 0) based on numerical testing using the experimental data of Shen and Dimotakis (1989) and Gerber (2002) for leading edge and midchord hydrofoil cavitation.Note that in this model for cavitation, two physical parameters are required, namely: 1. the mass fraction of non-condensable which, when assumed as spherical bubbles, provides an estimate of the number of nucleation sites available 2. a typical radius for the nuclei which has been assumed to be R b = 10 −6 m.

Numerical Adjustments for Cavitating Flows
The numerical prediction of cavitating flows for improved performance prediction is a difficult problem on several fronts.As already, it is difficult to mathematically model cavitation in a manner including all the significant scales involved without a heavy reliance on empirical data, which in turn are often unknown or lack generality.Secondly, the solution of the flow equations, in this case the full Navier-Stokes equations for viscous flow, is a highly non-linear system within which the cavitation model, which can be described as a violent vaporization process, must be implemented.The rapid vaporization leads to large spatial gradients in the flow mixture density field, up to four orders of magnitude for some flows, which pose considerable difficulties for convergence behavior.Thirdly, in flows where operation of the system is under extensive cavitation the two-phase system may operate close to or at the speed sound.Under such conditions it is important to include compressibility in the model that reflects the speed of sound characteristics of the fluid, as well as handling potential supersonic two-phase flow.For the present model specific enhancements required for the solution of cavitation problems were as follows: • Linearization of density in the mass equation (transient term only), using a speed of sound definition for a twophase mixture.• Due to the extremely large variations in mixture density, the momentum equations use a blending function to switch from second-order to first-order discretization in the cavitation zones.This adds robustness to the solution, particularly in high-density ratio cases and in areas of high fluid accelerations.• Since the majority of the flow regime is in the liquid phase (i.e., α m = 1) the discretization of the advecting velocities in the volume scalar equations uses upwind differencing (UDS).In the cavitation zone, where α m changes greatly, the many assumptions inherent in the cavitation model do not justify a higher order advection scheme other than the UDS scheme.• For the methodology to be fully conservative in mass and volume the evaluation of densities at controlvolume faces was based on a UDS scheme, this also has the added feature that in the case where the flow becomes supersonic, density behaves has a one-way variable.Applying the density in this manner is consistent with the UDS discretization for the volume scalar equations.It is important for mass/volume conservation that the face mass flows used in all mixture equations use a density evaluated consistent with the advection scheme employed for the volume scalar solution.
The above enhancements enabled stable and converged CFD solutions to be achieved at much higher density ratios.For example it is possible to obtain reliable solutions with density ratios up to 40,000 (corresponds to the density ratio present in 300 K water).
It also allows for a robust solution in cases where large vapor fractions are present.The cavitation model described here has been implemented in the commercial CFD code CFX-TASCflow (AEA, 2002), which follows a finite-volume/finite-element discretization procedure (Raw, 1985) in the development of the linearized equations.The solution of the pressure-velocity system is based on a coupled approach (mass and momentum solved simultaneously) and implemented within an Algebraic Multigrid strategy (Hutchinson et al., 1988;Ruge and Stüben, 1987).All discrete equations are second-order accurate in space.

Inducer Geometry
LEMFI Inducer (Figure 1) is designed for nominal operating conditions corresponding to a flow coefficient of ϕ = 0.38, a rotational speed of 1450 rpm and a pressure coefficient of ψ = 0.15.It consists of three blades without camber and with constant thickness and rounded leading edge.The main constructive parameters of the inducer are shown in Table 1.

General Overview of LEMFI Test Rig
The LEMFI-Paris pump test rig (Figure 2) is composed of the following main elements: • Two storage tanks with a capacity of 4 m 3 each, connected by a pipe of 350 mm in diameter.They can be loaded and emptied by means of two electrical control valves.
• A liquid ring vacuum pump is used to control the pressure at the free surface inside the storage tanks.• A 22 kW alternate current motor powered by a variable frequency controller was used to drive the tested inducer.The electric efficiency of the motor is given by the manufacturer.The rotational speed is measured using a magnetic tachometer (accuracy 0.1%).• A motorized control valve serves to adjust precisely the flow rate.
• The inducer equipped with a transparent acrylic cover.
• A circulation centrifugal pump installed in series with the impeller in order to overcome the circuit losses.• Various measurement instruments and devices: Ultrasonic flowmeter (accuracy 1%), placed at the inducer inlet.Two piezoresistif manometers (accuracy 1%).They are positioned at the inlet and outlet sections and measure the average tip pressure on four points.A temperature probe (accuracy 1%): the average temperature during the tests presented below is 18 • C.
The experimental cavitation procedure was done as fellow: at given rotational speed, the flow rate, measured by the flowmeter, is set to the operating value using the motorized control valve.The inlet pressure drop is realized by the liquid vacuum pump.For certain operating values of flow rates and inlet pressure, some specified types of cavitation were observed.For all tested flow rates, a sudden head drop happens with very low inlet pressure values.

COMPUTATIONAL METHODOLOGY
After a grid dependence solution study, a 250,000 structured mesh for a single blade passage was used for all the computations (Figure 3).This mesh was created using the mesh generator CFX-TurboGrid.The boundary conditions used are: total pressure at the inlet and mass flow at the outlet.The connection between the periodic faces is made by a periodic connections.The standard two-equation k-ε model was used for modeling the turbulence.The convergence criteria were set to 1e-4 on the maximal residual.For all the numerical simulations, the tip clearance was not taken into account.For each flow rate, the head drop curve was created as follow: First, an incompressible solution without cavitation is computed.From this noncavitating solution, the VOF model is turn on while the total pressure at the inlet is decreased by a constant step of 10,000 Pa.Close to the drop zone, this step is reduced by a factor of 10 and more to overcome the high instability of the solution due to the strong non-linear behavior of the cavitation.Thus, the head drop curve is created gradually.We note finally that the time consuming for the creation of a whole one head drop curve is about 8 h on a DEC machine having 2 CPU.

FIGURE 3
Typical mesh used for the simulation.

Overall Performances
For noncavitating regime, Figure 4 shows the evolution of the experimental and predicted values of the pressure coefficient versus the flow rate.As shown, a good agreement between the two results is obtained at nominal and high flow rates.At low

FIGURE 5
Experimental performances (uncertainty in N ± 2 rpm, in efficiency ±1%, in σ ± 0.002, in Q/Qn ± 0.02, in ψ ± 0.002).flow rates, the difference between the two results is due to the tip clearance effects.For several flow rates, Figures 5 and 6 show, respectively, the experimental and predicted values of the pressure coefficient versus the cavitation number σ .Figure 7 shows the experimental and predicted values of the critical cavitation number σ i versus the flow rate.This critical cavitation number corresponds to 3% of the head drop.

FIGURE 7
Evolution of the Incipient cavitation number vs flow rate.At low flow rates (qv/qn = 0.55 to qv/qn = 0.79), the drop curve occurs smoothly and slightly before the experimental one.Close to design flow rates (qv/qn = 0.91 to qv/qn = 1.09), the drop curve occurs suddenly and simultaneously with the experimental one.The agreement between the two results is very satisfactory.At high flow rates (qv/qn = 1.15 to qv/qn = 1.27), the drop curve occurs suddenly and slightly after the experimental one.For the very high flow rate (qv/qn = 1.27), the blockage phenomenon has occurred experimentally and not yet numerically.The VOF model under predicted the blockage phenomenon and could be attribute to the turbulence modeling and also to the presence of the coefficient Fc in Equation ( 13).

Cavitation Visualizations
For the computed results, the vapor rate is specified.For example on Figure 8 a pocket attached to blade is presented at design flow rate and for a vapor rate of 95%.The hub is colored in green, the blade in red, and the cavitation zone in blue.

Development of Cavitation with the Cavitation Number
The size of the predicted cavitation pocket corresponds to 10% of vapor in the mixture.This is leads to an incertainty of the experimental vapor percentage in the mixture and 10% was a realistic value to be used.Others authors (Coutier-Delgosa et al., 2002;Medvitz et al., 2002) propose, with succes, 5% regarding their works done on centrifugal pumps.
The shape of the cavitation pockets was compared to the experimental visualizations at flow rate qv/qn = 1.09.Figure 9 illustrates the general development of the cavitation in the inducer corresponding to the three operating points.As shown in Figure 9a, the vapor first appears near the shroud at the leading edge of the suction side of the blade.While the NPSH decreases, the cavity remains attached to the suction side of the blade but growing towards the blade to blade channel and down towards the hub (Figure 9b).At Figure 9c, the vapor pocket passes to the pressure side of the blade inducing a performance breakdown, a sudden drop of the performances appears.This is due to the blockage phenomenon generated by the presence of the pocket in the blade to blade channel from the tip to the hub and also in both side of the blade.The comparison between predicted and experimental results is well satisfactory notably for the localization of the cavitation pocket and its evolution.For the size of the pocket, it clearly depends on the percentage of the vapor in the mixture.

Evolution of the Cavitation vs the Flow Rate
The predicted cavitation behaviors are typical and in general in conformity with those visualized experimentally.presents the most representative images.One can thus identify: • In partial flow rate, backflow vortex cavitation returning upstream of the inducer (Figure 10a).• At the nominal flow rate, the cavitation attached to the blade (Figure 10b).• For high flow rate, stable cavities developed on both sides of the blade, characterizing the blockage phenomenon.The development of cavitation is almost identical in the three inducer channels (Figure 10c).

CONCLUSION
A cavitation model has been developed for flows with extensive two-phase vaporization and condensation.The model has been implemented in the commercial CFD software CFX-TASCflow 2.12.The validation of the model was performed, at various flow rates, on the cavitating behavior of an inducer tested at LEMFI.The validation concerned the head drop curves and the cavitation evolution.In general, a good agreement was obtained: head drops predictions are comparable to the experimental measurements and the size and location of the cavitation pockets observed during experiments are similar to the ones predicted by the numerical model.For high flow rate where the blockage occurs, the model under-predict the head drop location.Thus, the cavitation model itself requires careful testing for the determination of empirical constants relevant to the flow conditions.It is expected that the model will be useful in the design of hydraulic machines operating under conditions where large vapor cavities are present.Future studies will focus on the effect of the empirical constants on the performance predictions.The turbulence models effect will be also studied.For hydraulic devices when operating with large vapor cavities at high velocities, and two-phase speed of sound behavior must be considered.Such studies will necessarily focus on extended cavitation models that include non-equilibrium effects such as slip and thermal non-equilibrium not considered in this study.

FIGURE 2
FIGURE 2Hydrodynamic test bench of the LEMFI-Paris.
FIGURE 9Cavitation evolution vs cavitation number.(Predicted figures correspond to 10% of vapor in the mixture.) FIGURE 10Cavitation evolution vs flow rate.(Predicted figures correspond to 10% of vapor in the mixture.)

TABLE 1
Main Constructive Parameters of the Inducers g acceleration due to gravity [m 2 /s] NPSP = p t(inlet)p v net positive suction pressure [Pa] NPSH = NPSP