# Unsteady Cavitating Flows in Turbomachinery: Comparison of Two Numerical Models and Applications

Benoît Pouffary, Regiane . Fortes Patella, Jean-Luc Reboud

► To cite this version:

Benoît Pouffary, Regiane . Fortes Patella, Jean-Luc Reboud. Unsteady Cavitating Flows in Turbomachinery: Comparison of Two Numerical Models and Applications. ECCOMAS2004 Symposium, Jul 2004, Jyväskylä, Finland. hal-00211427

HAL Id: hal-00211427

https://hal.science/hal-00211427v1

Submitted on 2 Mar 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# UNSTEADY CAVITATING FLOWS IN TURBOMACHINERY: COMPARISON OF TWO NUMERICAL MODELS AND APPLICATIONS

Benoit Pouffary*, Regiane Fortes -Patella* and Jean-Luc Reboud†

*Institut National Polytechnique de Grenoble - LEGI -
BP 58 - 38041 Grenoble cedex 9 - France

e-mail: benoit.pouffary@hmg.inpg.fr, fortes@hmg.inpg.fr

† University Joseph Fourier Grenoble - CNRS-LEMD -
BP 166 - 38042 Grenoble Cedex 9 - France

e-mails: jean-luc.reboud@grenoble.cnrs.fr

**Key words:** cavitation, unsteady flows, pump, inducer, pressure correction method, artificial compressibility method

**Abstract.** Two numerical models for unsteady cavitating flows in turbomachinery are developed at The Turbomachinery and Cavitation team of LEGI (Grenoble). The first one was entirely developed in the laboratory for 2D unsteady flows, with the support of the French Space Agency CNES. It applies a pressure-correction method derived from the SIMPLE algorithm and a finite volume discretization on structured meshes. The second model for 3D flows has been developed in the FINE/TURBO™ CFD code in collaboration with the Numeca International society. The solver is based on an artificial compressibility method (preconditioning technique) with dual time stepping, adapted to the very large variation of the Mach number.

The models are based on a single fluid approach to describe the liquid-vapour mixture. Applications are mainly performed using a barotropic state law to manage the relation between the local static pressure and the mixture density. The alternative consisting in solving a supplementary equation for the void ratio including empirical terms for vaporization and condensation has been tested in the 2D model and some comparative results between both physical approaches are presented. The lecture presents also a comparison between results obtained by the two numerical models on 2D unsteady cavitating flows, making use of the barotropic approach. A discussion about the influence of the numerical and physical parameters is proposed. Finally, examples of applications to 3D cavitating flows in centrifugal pump and turbopump inducer geometry are presented.
# 1 INTRODUCTION

## 1.1 Background

As the first attempt to simulate unsteady cavitation phenomena Furness and Hutton [1] considered the deformation of the liquid-vapour interfaces to describe the early stage of re-entrant jet formation. For the simulation of complete cyclic behaviour, including vapour cloud shedding, the tracking of the liquid-vapour interface becomes a quite impossible challenge, because of multiple splitting of the main vapour structures and very quick vaporisation and condensation phenomena. An alternative approach has been proposed for about 15 years (Delannoy and Kueny, 1990 [2], Kubota et al., 1992 [3]). It consists in modelling the cavitating liquid as an homogeneous two-phase mixture of liquid and vapour. One main assumption in this case is to neglect the possible slip between the two phases, which leads to a single-phase fluid whose density may vary over a large range from pure liquid to pure vapour. This approach has been investigated in different ways:

- Delannoy and Kueny [2] proposed a formulation that strongly links the mixture density to the static pressure: they use a barotropic law $\rho(P)$, which describes the mixture density in both the incompressible parts of the flow field and the transition zone. This kind of model has been applied during the last decade by other researchers with different state laws (Song and He [4], Shin and Ikohagi [5], Ventikos and Tzabiras [6]...).

- Kubota et al. [3] proposed to relate the density evolution to the motion of bubbles in the flow. A given number of bubbles are settled at the inlet, and their evolution is governed by the Rayleigh-Plesset equation according to the pressure field. It is assumed that the bubbles are spherical, and remain not too close from each other. The void fraction is thus theoretically limited to a small value. However, that approach has been widely developed since that first work, for example by Chen and Heister [7], Grogger and Alajbegovic [8], Bunnel and Heister [9] or Sauer and Schnerr [10].

- Different authors proposed more recently to consider a transport equation model for the void ratio, with vaporization/ condensation source terms to control the mass transfer between the two phases (Shingal et al. [11], Merkle et al. [12], Kunz et al. [13], Senocak and Shyy [14]...). This method has the advantage that it can take into account the time influence on the mass transfer phenomena, through empirical laws for the source term.

The main numerical problem in multidimensional simulations is the simultaneous treatment of two very different flow conditions: two almost incompressible ones (pure liquid and pure vapour), and a highly compressible one in the transition between vapour and liquid. Most of the methods have serious difficulties when the ratio $\rho_v/\rho_l$ is lowered. Mainly two numerical methods have been developed to simulate unsteady cavitation:

- The first method adapted to incompressible fluids computations, is the pressure correction method, based on the SIMPLE scheme. Delannoy and Kueny [2], adapted this algorithm to cavitation by considering the mass equation as a transport equation for density, which depends on the pressure through a barotropic state law. This method, first developed for
inviscid fluids, considers a physical speed of sound in the vapour/liquid transition and thus
captures the very strong density gradient in the mixture [2, 6, 14].

- The second method is based on the adaptation of compressible time-marching algorithms to low Mach numbers: artificial compressibility method. This kind of resolution was originally devoted to highly compressible flows. In the case of low-compressible or incompressible simulations the efficiency of the numerical scheme decreases dramatically. That problem was solved (Turkel [15]) by multiplying the pseudo-time derivatives by a preconditioning matrix that modifies the equations and accelerates the convergence, without altering the result accuracy provided each time step is correctly converged. Choi and Merkle [16] implemented it in an implicit algorithm, and they obtained satisfactory results, either with a barotropic law or a two mass equations model. Kunz et al. [17] presented results obtained with this algorithm, using a three-fluid method based on two mass transfer equations and including the presence of a non-condensable gas.

**1.2 Development and applications of two numerical models**

The Turbomachinery and Cavitation team of LEGI (Grenoble-France) has been developing
models of unsteady cavitating flows for many years, with the constant support of the French
Space Agency CNES – Centre National d'Etudes Spatiales.

Firstly, a numerical model "IZ" was entirely developed in the laboratory for 2D unsteady
flows. Initially based on the Euler equations, it applies a pressure-correction method derived
from the SIMPLE algorithm and a finite volume discretization on structured meshes [2].
Successive improvements of the 2D model "IZ" and applications to different flow
configurations were performed during the last decade [18-27]. The development of a second
model for 3D flows is in progress, based on the 3D RANS FINE/TURBOTM CFD code in
collaboration with the Numeca International society. The solver is based on an artificial
compressibility method with dual time stepping, adapted to the very large variation of the
Mach number trough a preconditioning technique. It has already been applied to some
cavitating flows in test sections and turbomachinery [28-34]. That work, still in progress, has
the final objective to provide assistance to the design and prevision of operating range of
rocket engine turbopumps, taking into account the steady state and unsteady effects of
cavitation.

The two numerical models are based on a single fluid approach to describe the liquid-vapour mixture, with a mixture density $\rho$ varying in the flow field between the vapour density and the liquid density. Applications are mainly performed using a barotropic state law to manage the relation between the local static pressure and the mixture density. The alternative consisting in solving a supplementary equation for the void ratio including empirical terms for vaporization and condensation has been already tested in the 2D code, and some comparative results between both physical approaches are presented in the paper.

The lecture presents also a comparison between results obtained by the two numerical
models on 2D unsteady cavitating flows, making use of the barotropic approach. A discussion
about the influence of the numerical and physical parameters is proposed.
Finally, an overview of the applications to 3D cavitating flows in centrifugal pumps and turbopump inducers geometry is presented.

## NOMENCLATURE

<table><tr><td>c<sub>min</sub></td><td>minimum speed of sound in the medium</td><td>(m/s)</td></tr><tr><td>Cp</td><td>(P-P<sub>ref</sub>)/(ρ<sub>ref</sub>V<sub>ref</sub><sup>2</sup>/2) pressure coefficient</td><td>-</td></tr><tr><td>k</td><td>turbulent kinetic energy per mass unit</td><td>(m<sup>2</sup>s<sup>-2</sup>)</td></tr><tr><td>L<sub>ref</sub></td><td>geometry reference length</td><td>(m)</td></tr><tr><td>P</td><td>local static pressure</td><td>(Pa)</td></tr><tr><td>P<sub>ref</sub></td><td>reference pressure</td><td>(Pa)</td></tr><tr><td>P<sub>v</sub></td><td>vapour pressure</td><td>(Pa)</td></tr><tr><td>V<sub>ref</sub></td><td>reference velocity (=V<sub>upstream</sub>)</td><td>(ms<sup>-1</sup>)</td></tr><tr><td>α</td><td>void ratio</td><td>-</td></tr><tr><td>ε</td><td>turbulence dissipation per mass unit</td><td>(m<sup>2</sup>s<sup>-3</sup>)</td></tr><tr><td>μ<sub>l</sub>, μ<sub>t</sub></td><td>laminar and turbulent dynamic viscosity</td><td>(Pa.s)</td></tr><tr><td>ρ</td><td>mixture density</td><td>(kg/m<sup>3</sup>)</td></tr><tr><td>ρ<sub>l</sub> (=ρ<sub>ref</sub>), ρ<sub>v</sub></td><td>liquid (=ref), vapor density</td><td>(kg/m<sup>3</sup>)</td></tr><tr><td>σ</td><td>(P<sub>ref</sub>-P<sub>v</sub>)/(ρV<sub>ref</sub><sup>2</sup>/2) cavitation number</td><td>-</td></tr><tr><td>ω</td><td>specific dissipation rate</td><td>(s<sup>-1</sup>)</td></tr></table>

## 2 MODEL BASES

We use a single fluid model to describe the liquid-vapour mixture, with a mixture density ρ varying in the flow field between the vapour density and the liquid density. The equivalent void ratio α of the liquid-vapour mixture relates to the varying specific mass by ρ = α ρ<sub>v</sub> + (1 - α) ρ<sub>l</sub>. In the mixture, the velocities of liquid and vapour phases are the same and we obtain only one set of equations for the mixture mass, momentum, temperature or turbulence (k-ε or k-ω), written in their conservative form.

The void ratio α of the mixture depends on the local static pressure. The main work was performed in our laboratory, as initially proposed by Delannoy and Kueny, [2], using 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 ρ<sub>l</sub> and vapour density ρ<sub>v</sub> when the local static pressure P is around the vapour pressure P<sub>v</sub>. The law is characterised by its maximum slope at P = P<sub>v</sub>, which is related to the minimum speed of sound c<sub>min</sub> in the two-phase homogeneous medium. No significant effect of the ratio ρ<sub>v</sub>/ρ<sub>l</sub> is observed, provided its value is imposed smaller than 0.01. c<sub>min</sub> is then the only adjustable parameter of the model, found from comparisons with the experimental results of Stutz and Reboud [35-37] close to 1.5 to 2m/s. That order of magnitude is kept constant for the different applications with water.
<figure><img src="image_2.png" /><figcaption>Figure 1: Barotropic state law ρ(P). Water 20°C.</figcaption></figure>

The alternative consisting in solving a transport equation for the void ratio including empirical terms for vaporization and condensation was also implemented in the 2D model, and compared with the barotropic model. We use the following form derived from Kunz et al. [13]:

$$ \frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha U) = - \frac{1}{\rho v} (\dot{m}_{VL} + \dot{m}_{LV}) $$

$$ \dot{m}_{LV} = C_{LV} \cdot \rho_{L} (1 - \alpha) \text{Min}(0, C_{p} + \sigma) / \tau $$

$$ \dot{m}_{VL} = C_{VL} \cdot \rho_{L} \alpha (1 - \alpha)^2 / \tau $$

(τ: characteristic time, C<sub>LV</sub>, C<sub>VL</sub> empirical constants)

## 2.1 - 2D model “IZ”

Reynolds-averaged Navier-Stokes equations are solved in the case of this single fluid with variable density. A finite volume spatial discretization is applied in curvilinear orthogonal co-ordinates on a staggered mesh. An iterative resolution based on the SIMPLE algorithm was developed to deal with quasi-incompressible flow (α=0 and α = 1) and highly compressible flow (0<α<1) [2]. The liquid-vapour interfaces are described by high gradients of the mixture density ρ, which was made possible by using a conservative approach and the HLPA non-oscillatory second order MUSCL scheme proposed by Zhu [38]. An estimation of the term dp/dP is needed to apply the pressure correction method to the mixture with the strong variations of compressibility. The barotropic state law ρ(P) directly gives that term, in the other case an estimation of dp/dP is deduced from the vaporization term of the transport equation for the void ratio.
To solve the time-dependant elliptic problem first or second order fully implicit methods are available. Turbulent flows are calculated by solving the Reynolds equations using a $k-\epsilon$ RNG turbulence model with standard laws of the wall or the $k-\omega$ model of Wilcox [39].

Classically, the mass flow rate is imposed constant as upstream boundary condition. The difference between the static pressure imposed as downstream boundary condition and the vapour pressure, is slowly decreased from a high value leading to steady non-cavitating conditions, down to the constant value required to assure the expected cavitation number $\sigma$. Phase change then occurs spontaneously in the regions where the pressure decreases close to the vapour pressure $P_v$.

## 2.2 - 3D Model

To calculate 3D cavitating flows the barotropic cavitation model was implemented in the commercial code FINE/TURBO™ developed by NUMECA International. FINE/TURBO™ is a three-dimensional structured mesh code that solves the time dependant Reynolds-averaged Navier-Stokes equations, with artificial compressibility method. A detailed description of the code is given in Hakimi [40]. 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 code resorts to a multi-grid strategy to accelerate the convergence in non-cavitating conditions, associated with a local time stepping and an implicit residual smoothing. In the case of low-compressible or incompressible flows, the pseudo-time derivatives are multiplied by a preconditioning matrix (Hakimi [40]) based on the studies of Turkel, [15], Choi and Merkle, [16].

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. 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.

# 3. SELF OSCILLATION BEHAVIOUR IN 2D VENTURI TEST SECTION

## 3.1 Geometry and conditions

Many numerical simulations were performed on a Venturi type section, the convergent and divergent angles of which are respectively about 18° and 8°. The shape of the Venturi bottom downstream from the throat simulates an inducer blade suction side with a bevelled leading edge geometry and a chord length $L_{ref} = 224$ mm. According to experimental observations in this geometry cavitation sheets develop from the throat and show a typical self-oscillation behaviour [35, 37], with quasi-periodic fluctuations. For a cavitation number $\sigma$
of about 2.4 (based on the time-averaged upstream pressure) and an inlet velocity $V_{ref} = 7.2$ m/s, the vapour shedding frequency observed experimentally is about 50 Hz ($\approx 1.5 V_{ref}/L_{ref}$) for a maximum attached cavity length of 45mm ($\approx 0.2 L_{ref}$). +/- 5mm. That gives the classical Strouhal number, based on the cavity length: $St \approx 0.3$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Curvilinear- orthogonal mesh of the Venturi type section (160x50 cells). $L_{ref}=224$ mm</figcaption></figure>

## 3.2 Self oscillation behaviour, barotropic model

Computations are performed with a k-ε RNG model, with laws of the wall and an empirical modification consisting in reducing the turbulent viscosity in the liquid-vapour mixture regions [19, 24]. Numerical simulation of several cloud shedding cycles (Fig. 3) leads to a reliable simulation of the Venturi unsteady cavitation behaviour. The shedding frequency and mean cavity development are found in fairly good agreement with experimental ones [25]. The influence of the speed of sound is illustrated on Fig. 4, and the difference observed between the results obtained with $c_{min}$ 1.5 and 2.25 m/s remains very small.

Results are consistent with those obtained with a k-ω model, including compressibility effects as proposed by Wilcox [39]. According to the numerical calculations, the fluid compressibility in the mixture regions seems to play a significant role in the turbulence structure, and must be taken into account to simulate the local unsteadiness generated by cavitating flows [24].

The satisfactory results obtained with the modified k-ε RNG model have been confirmed by simulations performed in other geometries, such as a hydrofoil [20], a foil cascade [23], and another Venturi type section leading to a more stable cavitation sheet [25]. In all cases, the general experimental behaviour is correctly obtained, and oscillation frequencies are well predicted in unsteady configurations. Moreover, velocity and void ratio distribution within the
cavity compare favourably with optical probes measurements performed in the same geometry
[19, 25, 35-37] :

<figure><img src="image_4.png" /><figcaption>Figure 3: Time evolution (in abscissa) of the cavity length (graduated in ordinate). Instantaneous density distribution of attached and cloud cavities are drawn on the left at T = 26 T<sub>ref</sub>. Calculation conditions: σ≈2.4; V<sub>ref</sub> = 7.2 m/s; L<sub>ref</sub> = 0.224 m, c<sub>min</sub> = 1.5 m/s, mesh = 173x61 - time step Δt = 0.003T<sub>ref</sub> (T<sub>ref</sub> = L<sub>ref</sub> / V<sub>ref</sub>), modified k-ε RNG model.</figcaption></figure>

However, in the calculations, the static pressure is imposed as the downstream boundary
condition, while the reference pressure is measured in a section upstream to the throat of the
Venturi. The losses generated in the cavitation tunnel by viscous and unsteady effects have to
be taken into account in the calculation to compare and adjust a *posteriori* the experimental
and numerical cavitation numbers. Moreover, because of the boundary condition fixed, the
pressure in the inlet section (where the flow rate is imposed) and the flow rate in the outlet
section (where the static pressure is imposed) are fluctuating with respect to the cavitation
cloud shedding cycles and the corresponding variations of the total vapour volume (Fig. 4). A
time-averaged value of the fluctuating upstream pressure is then used to compute the effective
cavitation number σ of the simulated case. One can observe that the order of magnitude of the
total vapour volume relative fluctuations is +/- 10% and that the associated downstream flow
rate fluctuations remain smaller than +/- 1%. Contrarily, the amplitude of the upstream pressure fluctuation is found very high, with a range $\Delta Cp$ close to 1.5.

<figure><img src="image_5.png" /><figcaption>Figure 4: Barotropic model :
Fluctuation of : (a) pressure coefficient (inlet in blue, outlet in rose)
(b) Flow rate (inlet in blue, outlet in rose)
(c) Total vapour volume</figcaption></figure>

Comparison of two values of the celerity of sound $c_{min}$ : 1.5 m/s (left) and 2.25 m/s (right)

To explain the high pressure fluctuations, we can note on the one hand, that the barotropic state law does not take into account the influence of the characteristic time of vaporization and condensation processes, and on the other hand, that the circuit impedance is ignored by the steady boundary conditions imposed at the tunnel upstream and downstream ends. Consequences are that the collapses of the vapour clouds generated during the cavitation cycles are very strong and that the associated pressure fluctuations, in the quasi-incompressible liquid, are not properly distributed on both ends of the tunnel. To solve that problem, we are testing two modifications of the model:

- the barotropic approach can be replaced by the transport equation for the void ratio presented previously. Because it takes into account a characteristic time limiting the phase
change, a sudden increase of pressure does not lead so quickly to the total disappearance of all the vapour structures.

- the effect of the hydro-elastic behaviour of upstream and downstream pipes can be taken into account, and the procedure used is described thereafter.

## 3.2 Influence of the cavitation model

The barotropic state law is replaced by the resolution of a supplementary equation for the void ratio with empirical source terms:

$$ \frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha U) = \frac{1}{\rho v} (\dot{m}_{VL} + \dot{m}_{LV}) $$

$$ \dot{m}_{LV} = C_{LV} \cdot \rho_{L} (1 - \alpha) \text{Min}(0, C_{P} + \sigma) / \tau \quad (\tau: \text{characteristic time}) $$

$$ \dot{m}_{VL} = C_{VL} \cdot \rho_{L} \alpha (1 - \alpha)^2 / \tau \quad (C_{LV}, C_{VL} \text{ empirical constants}) $$

<figure><img src="image_6.png" /><figcaption>Figure 5: Time evolution (in abscissa) of the cavity length (graduated in ordinate). Instantaneous density distribution of attached and cloud cavities are drawn on the left at T = 26.2 T<sub>ref</sub>. Calculation conditions: σ≈2.4; V<sub>ref</sub> = 7.2 m/s; C<sub>VL</sub>/τ = 3.5, C<sub>LV</sub>/τ = 1.5, mesh = 173x61 - time step Δt = 0.003T<sub>ref</sub> (T<sub>ref</sub> = L<sub>ref</sub> / V<sub>ref</sub>), modified k-ε RNG model.</figcaption></figure>
An example of result of such a modification is illustrated in figure 5. With the parameters chosen (transport equation of $\alpha : C_{LV}/\tau \approx 3.5$, $C_{VL}/\tau \approx 1.5$, $\rho v/\rho l = 10^{-3}$) it can be noticed that the cycle is less regular than with the barotropic model, the cavitation structures obtained are about two times larger (fig. 3) and the cycle frequency about 50% smaller. The instantaneous void ratio is smaller, more particularly in the upstream part of the cavity because of the characteristic time limiting the vaporisation, but the averaged vapour volume (fig. 6) remains larger because of the larger cavity obtained.

<figure><img src="image_7.png" /><figcaption>Figure 6: Transport equation for the void ratio (legend as fig. 4)
Comparison of two sets of evaporation-condensation parameters $C_{LV}/C_{VL}$
(left : $C_{VL}/\tau = 3.5$, $C_{LV}/\tau = 1.5$ ; right : $C_{VL}/\tau = 32$, $C_{LV}/\tau = 16$)</figcaption></figure>

Influence of the parameters is complex because vaporization and condensation terms can be adjusted separately but are in interaction in the right hand of the transport equation. An example is drawn on the right part of fig 6 with $C_{LV}/\tau \approx 32$, $C_{VL}/\tau \approx 16$. The average vapour volume appears smaller than in the previous case, because of the larger vaporisation and condensation constants, but the cloud shedding cycles remain more irregular than those obtained with the barotropic model.
### 3.3 Influence of the upstream and downstream conditions

To take into account the circuit impedance we chose to modify the fixed inlet and outlet boundary conditions with respect to the upstream and downstream pipes hydro-elastic behaviour. Because the computation is purely transient, we used a simplified 1D hydro-elastic model that solves the Allievi's equations taking into account the speed of propagation of pressure and flow rate perturbations in the pipes, function of their metal wall thickness and diameter (Longatte [41]). Non-reflecting boundary conditions are applied at one end of each pipe, the other end being coupled with the 2D computational domain of the cavitation tunnel. So at each time step an iterative procedure allows coupling the inlet and outlet flow rate and pressure with respect to the cavitation behaviour (pressure and vapour volume fluctuations due to the cloud shedding process) and the propagation at the given speed of sound of the pressure and flow rate fluctuations in the pipes (Figure 7).

So pressure and flow rate fluctuations can leave the cavitation tunnel on either side, with respect to the hydro-acoustic characteristic of the connected pipes. Slow pressure or flow rate variations (i.e. with a characteristic time much larger than the cavitation cycle period) are estimated during the computation by local time averaging, and are corrected to impose a given time-averaged value of flow rate and cavitation number.

<figure><img src="image_8.png" /><figcaption>Figure 7: Coupling between the 3 domains: cavitation tunnel, upstream and downstream pipes</figcaption></figure>

Influence of the pipes is easily noticeable on figure 8: Pressure and flow rate are fluctuating at the two ends of the Venturi. Approximately one third of the pressure fluctuation leaves the cavitation tunnel by its downstream pipe, while the vapour volume variations create approximately opposite flow rate fluctuation amplitude at the inlet and outlet of the tunnel. Amplitude of pressure, flow rate and vapour volume fluctuations are found about two times smaller than in the previous cases (fig 4 and 6, left).

Influence of the pipes characteristics can be studied and the hydro-elastic model itself might be improved to take into account more complex hydraulic characteristics of the circuit instead of non-reflecting boundary conditions.
<figure><img src="image_9.png" /><figcaption>Figure 8: Coupling with upstream and downstream pipes (legend as fig. 4)
Pipes diameter 0.24 m, speed of waves propagation in the pipes 1000 m/s.</figcaption></figure>

- left : barotropic model, $C_{min} = 1.5$ m/s (to be compared with Fig. 4, left)

- right : transport equation for $\alpha : C_{VL}/\tau = 3.5$, $C_{LV}/\tau = 1.5$ (to be compared with Fig. 6, left)

## 3.4 Comparison with the 3D model

Computation of unsteady 2D flows is one of the test cases of the 3D model. A first application to unsteady cavitating flows was performed on the 2D hydrofoil proposed as a test case for the workshop on physical models and CFD tools for computation of cavitating flows held in 2003 within the 5th Int. Symposium on Cavitation [34]. On this blind-test, a general good agreement with the results of the other numerical codes was found. Now looking more in details to the local behaviour, we note that the good convergence of instantaneous flow rate upstream and downstream from the strongly and rapidly fluctuating cavitation area remains a difficult challenge with the dual time-stepping approach, when trying to use a limited number of pseudo-time steps. That work is still in progress, and comparison between our two numerical models remains up to now only qualitative.
<figure><img src="image_10.png" /><figcaption>Figure 8: Comparison of instantaneous density field
a) pressure correction method ; b) pseudo compressibility approach</figcaption></figure>

4. 2D inducer blade cascades.

2D blade cascade are obtained from inducer geometry, by a cut at constant radius (Fig. 9). Because of the high rotation speed of such axial pumps, the design angle and the corresponding blade angle in the cascade are very low (typically 10 to 15degrees).

<figure><img src="image_11.png" /><figcaption>Figure 9 : From the 3D inducer to the blade cascade.</figcaption></figure>

Computations are performed with the two numerical models on a single channel with periodicity conditions. In the present case the two codes apply a standard k-ε turbulence model with the result that the strong turbulence viscosity in the cavity wake prevents the appearance of the self-oscillation behaviour [24]. Meshes are different: we use for the pressure correction method (2D model) an H-mesh (465x60) around the blade and for the 3D model a mesh (369x61) of one channel between two blades (figure 10). Because of the low blade angle, an important influence of the mesh on the non cavitating performance of the cascade was found. That is why we present only relative head curve to compare the cavitation
behaviour obtained with the two models. The cavitation number is calculated *a posteriori* in the two cases with the upstream static pressure as reference pressure.

As can be noticed on the figure 10 and 11, the development of the cavitation sheets on the blade and associated head drop are found very similar with the two models. Moreover, very large cavities corresponding to more than 20% head drop can be predicted with a good convergence level.

<figure><img src="image_12.png" /><figcaption>Figure 10: Development of cavitation sheet on the blade with decreasing cavitation number
a) pressure correction method ; b) pseudo compressibility approach</figcaption></figure>

One source of unsteadiness in turbo-pump inducers consists in rotating cavitation behaviour, characterized by different cavity shapes on the different blades, and leading to super or sub synchronous disturbances. We simulated this phenomenon in the case of 2D blade cascades corresponding to typical four-blade inducers [22, 26]. Four-channel computations were performed with the 2D model, with the standard $k-\epsilon$ turbulence model to prevent the self-oscillation of the cavitation sheets, and the results were compared with those of the single-channel computations. Non-symmetrical unstable flow patterns were obtained for certain operating conditions, and limits of stability according to the mass flow rate and the cavitation number are obtained [22, 26]. Such unsteady behaviour is basically different from the self-oscillation behaviour presented before, because it involves strong interactions between the cavity volume and the global structure of the liquid flow around the blades. Work is still in progress to compare the two models on such configurations.
<figure><img src="image_13.png" /><figcaption>Figure 11 : Comparison of relative head drop vs. cavitation number
a) pressure correction method ; b) pseudo compressibility approach</figcaption></figure>

## 4. 3D STEADY CAVITATING FLOWS IN TURBOMACHINERY

3D computation of unsteady cavitating flows in turbomachinery is the final objective of the development of the cavitation model implemented in FINE/TURBO™. In parallel to the previous studies of unsteady flows in 2D configuration, applications have already been preformed to steady 3D cavitating flows in turbomachinery. As in the blade cascade case, steadiness of the cavitation behaviour is obtained with the code by using a standard $k-\epsilon$ turbulence model on single channel meshes with periodicity conditions.

Illustrations presented hereafter concern results obtained in the SHF centrifugal pump tested at Electricité de France [30], and in inducer geometry.

### 4.1 Centrifugal pump

The SHF pump is a shrouded seven blades centrifugal runner. The scale and the calculation configuration are 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 3000 rpm, and the nominal flow rate is 0.157 $m^3$/s. One blade-to-blade channel is mesh of 250000 cells is used. Cavitating behaviour obtained for the centrifugal pump at nominal flow rate shows thin attached sheet cavities on the suction side of the impeller, their length increasing with the decreasing cavitation number (Fig. 12). Head drop begins approximately when the cavitation sheet reaches the position of the previous blade leading edge, and increases rapidly to more than 10% (Fig. 13).
<figure><img src="image_14.png" /><figcaption>Figure 12: Comparison experiment/calculation of the suction side cavity extension for several NPSH values, SHF pump. Experimental visualizations by EDF [30].</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Head drop charts at nominal flow rate. SHF pump: Experimental head drop by EDF [30].</figcaption></figure>

## 4.2 Turbopump inducer

To achieve high rotational speed and low inlet pressure, rocket engine feed pumps are equipped with an inducer stage, in which developed cavitation occurs at nominal point of operation. Cavitation observed in turbopump inducers is characterized by complex, three-dimensional, two-phase structures. These may appear as cavitation sheets attached to blades suction side or in rotating flows at the periphery of the runner (blade tip vortices or inlet back flow). First computations in inducers geometry show that the 3D flow and cavitation structures around the leading edge are very complex [28]. As a future test case for unsteady cavitation and the related instabilities, we are calculating the cavitating flow on a inducer with pure radial leading edge shape (Fig. 14) that as been studied experimentally at the University of Osaka [42]. First result presented in figure 14 illustrates the
development of attached cavitation sheet on the suction side of the impeller. That result is obtained with a mesh of 250000 cells for one blade-to-blade channel, without taking into account the tip leakage. Despite that simplification of the real geometry, the order of magnitude of the critical cavitation number is found in good agreement with the experimental range with a final head breakdown observed for $\sigma$ close to 0.03.

<figure><img src="image_16.png" /><figcaption>Figure 14: Inducer geometry and suction side cavity for $\sigma = 0.05$</figcaption></figure>

## ACKNOWLEDGMENTS

The authors wish to express their gratitude to the French space agency CNES and SNECMA-Moteurs for their support.

## REFERENCES

[1] Furness R.A., Hutton S.P. (1975), "Experimental and Theoretical Studies of two-dimensional Fixed-Type Cavities", *J. Fluid Engineering*, vol. 97, pp. 515-522.

[2] Delannoy, Y., Kueny, J.L. (1990), "Two phase flow approach in unsteady cavitation modelling", Cavitation and Multiphase Flow Forum, ASME-FED vol.98, pp. 153-158.

[3] Kubota A, Kato H., Yamaguchi H. (1992), "A new modelling of cavitating flows: a numerical study of unsteady cavitation on a hydrofoil section", *J. Fluid Mech.*, vol. 240, pp. 59-96.

[4] Song C. and He J., "Numerical simulation of cavitating flows by single-phase flow approach"; *3rd Int. Symp. on Cavitation*, Grenoble, France, April 1998.

[5] Shin B.R. and Ikohagi T, 1999, "Numerical analysis of unsteady cavity flows around a hydrofoil", ASME-FEDSM 99-7215, San Francisco.

[6] Ventikos, Y., and Tzabiras, G., 2000, "A numerical Method for the Simulation of Steady and Unsteady Cavitating Flows," *Comput Fluids*, 29, 63-88.

[7] Chen Y. and Heister S.D. (1995), "Modeling hydrodynamic Non-Equilibrium in Bubbly and Cavitating Flows", *J. Fluids Eng.*, vol. 118, n°1, pp172-178.

[8] Grogger H.A. and Alajbegovic A. (1998), "Calculation of the Cavitating Flow in Venturi Geometries Using Two Fluid Mode", ASME FEDSM 99-7364.

[9] Bunnell R.A. and Heister S.D. (2000) "Three-dimensional unsteady simulation of cavitating flows
in injector passages", *J. Fluid Eng.* Vol 122, pp 791-797

[10] Sauer J., and Schnerr G.H. (2000), Unsteady cavitating flow - A new cavitation model based on modified front capturing method and bubble dynamics, Proceedings of ASME FEDSM'00 – Boston.

[11] Singhal, A.K., Vaidya, N., and Leonard, A.D., 1997, "Multi-dimensional Simulation of Cavitating Flows Using a PDF Model for Phase Change," *ASME FEDSM* 97-3272

[12] Merkle, C.L., Feng, J., Buelow, P.E.O. (1998), "Computational modeling of the dynamics of sheet cavitation", *3rd Int. Symp. on Cavitation*, Grenoble, France.

[13] Kunz, R., Boger, D., Chyczewski, T., Stinebring, D., Gibeling, H. (1999), "Multi-phase CFD analysis of natural and ventilated cavitation about submerged bodies", *ASME FEDSM* 99-7364, San Francisco.

[14] Senocak, I., and Shyy, W., 2002, "A Pressure-Based Method for Turbulent Cavitating Flows," *J. Comput. Phys.*, 176, pp. 363-383.

[15] Turkel, E. (1987), "Preconditioning methods for solving the incompressible and low speed compressible equations", *J. of Comp. Phys.*, vol 72, pp. 277-298.

[16] Choi, D., Merkle, C. L. (1993), "The application of preconditioning in viscous flows", *J of Comp. Phys.*, vol. 105, pp. 207-223.

[17] Kunz, R.F., Boger, D. A., Stinebring, D. R., Chyczewski, T.S., Lindau, J.W., and Gibeling, H. J., 2000 "A Preconditioned Navier-Stokes Method for Two-phase Flows with Application to Cavitation," *Comput Fluids*, 29, pp.849-875.

[18] Reboud, JL., Delannoy, Y (1994), "Two-phase flow modelling of unsteady cavitation", *2nd Int. Symp. on Cavitation*, Tokyo.

[19] Reboud, JL., Stutz, B., Coutier, O. (1998), "Two phase flow structure of cavitation : experiment and modeling of unsteady effects" *3rd Int. Symp. on Cavitation*, Grenoble, France.

[20] Reboud, JL., Fortes-Patella, R., Hofmann, M., Lohrberg, H., Ludwig, G., Stoffel, B. (1999), "Numerical and experimental investigations on the self-oscillating behavior of cloud cavitation", *ASME-FEDSM* 99-6755/7259, San Francisco.

[21] Coutier-Delgosha O., Reboud J-L., Albano G. (2000) "Numerical Simulation of the Unsteady Cavitating Behaviour of an Inducer Blade Cascade", *ASME FEDSM00*, Boston

[22] Joussellin F., Courtot Y., Coutier-Delgosha O., Reboud JL. (2001) "Cavitating inducer instabilities : experimental analysis and 2D numerical simulation of unsteady flow in blade cascade", *4th Int. Symp. on Cavitation*, Pasadena.

[23] Lohrberg H., Stoffel B., Fortes-Patella R., Coutier-Delgosha O. Reboud JL. (2002) : "Numerical and Experimental Investigations on the Cavitating Flow in a Cascade of Hydrofoils", *Experiments in Fluids*, 33/4 : pp 578-586

[24] Coutier-Delgosha O., Fortes-Patella, R., Reboud J.L. (2003), "Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation", *J of Fluids Eng.*, Vol. 125 pp. 38-45.

[25] Coutier-Delgosha O., Reboud J.L., Delannoy Y. (2003): "Numerical simulation of the unsteady behaviour of cavitating flows, *Int. J. Num. Methods in Fluids*, Vol. 42, Issue 5, pp. 527-548.

[26] Coutier-Delgosha O., Perrin J., Fortes-Patella R., Reboud J.L (2003) " A numerical model to predict unsteady cavitating flow behaviour in inducer blade cascades ", *Fifth int. Symp. on Cavitation*, Osaka, Japan.

[27] Coutier-Delgosha O., Reboud J.L., Fortes-Patella R. (2002) "Numerical study of the effect of the leading edge shape on cavitation around inducer blade sections", *JSME International Journal – serie B Vol. 45, No. 3, 2002*

[28] Coutier-Delgosha O., Morel P., Fortes-Patella R., Reboud JL. (2002) : "Numerical simulation of
turbopump inducer cavitating behavior", *ISROMAC-9*, Honolulu.

[29] Coutier-Delgosha O., Fortes-Patella R., Reboud JL., Pouffary B. (2002) : "3D numerical simulation of pump cavitating behavior", *ASME-FEDSM* 02-31188, Montreal, Canada.

[30] Coutier-Delgosha O., Pouffary B., Fortes-Patella R., Reboud J.L., Archer A., Combes J-F. (2002) "Cavitation performance of a centrifugal pump: numerical and experimental investigations" 21st IAHR Symp. on Hydraulic Machinery and Systems, Lausanne, Sept. 2002.

[31] Pouffary B., Coutier-Delgosha O., Fortes-Patella R., Reboud J.L., Laffite S., Nguyen Duc J-M. (2002) "Evaluation of the effects of cavitation on the flow in a centrifugal pump: analysis of the breakdown mechanisms", 21st IAHR Symp. on Hydraulic Machinery and Systems, Lausanne, September 2002.

[32] Coutier-Delgosha O. - Fortes-Patella R. - Reboud J.L. - Hofmann M. - Stoffel B. (2003) : "Experimental and numerical studies in a centrifugal pump with 2D-curved blades in cavitating conditions", *J of Fluids Eng*. Vol. 125, Issue 6, pp. 970-978

[33] Pouffary B., Lorrain E., Fortes-Patella R., Reboud JL. (2003) "Numerical simulation of cavitating flows in space pumps", *MACSI-NET workshop on industrial challenges in the simulation of evolving interfaces*, Brussels, Belgium, September 2003.

[34] Pouffary B., Fortes-Patella R., Reboud JL. (2003) "Numerical simulation of cavitating flow around a 2d hydrofoil: a barotropic approach", *Fifth Int. Symp. on Cavitation*, Osaka, Japan.

[35] Stutz B., Reboud JL., (1997), "Experiments on unsteady cavitation", *Experiments in Fluids*, 22, pp. 191-198.

[36] Stutz, B., Reboud, JL (1997), "Two-phase flow structure of sheet cavitation", *Phys. Fluids* 9(12).

[37] Stutz, B., Reboud, JL (2000), "Measurements within unsteady cavitation", *Experiments in Fluids*, n°29, pp545-552

[38] Zhu J. (1991), "A low diffusive and oscillation-free convection scheme", *Comm. in Applied Num. Methods*, vol. 7.

[39] Wilcox D. (1998), *Turbulence modeling for CFD*, DCW Industries, Inc., La Canada, California, USA.

[40] Hakimi, N. (1997): "Preconditioning methods for time dependent Navier-Stokes equations", Ph.D.Thesis, Vrije Univ. Brussels.

[41] Longatte, F. (1998), "Contribution à l'Analyse Phénoménologique des Ecoulements Instationnaires dans les Turbomachines: Etude du Couplage Pompe-Circuit et Rotor-Stator", Ph.D. Thesis, INPGrenoble, France.

[42] Yoshida Y., Tsujimoto Y., Kataoka D, Horiguchi H., Wahl F. (2001), "Effects of Alternate Leading Edge Cutback on Unsteady Cavitation in 4-Bladed Inducers", *J of Fluids Eng.*, vol.123, pp.762-770