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

Olivier Coutier-Delghosa, Regiane Fortes Patella, Jean-Luc Reboud, Nouredine Hakimi, Charles Hirsch

► To cite this version:

Olivier Coutier-Delghosa, Regiane Fortes Patella, Jean-Luc Reboud, Nouredine Hakimi, Charles Hirsch. Numerical simulation of cavitating flow in 2D and 3D inducer geometries. International Journal for Numerical Methods in Fluids, 2005, 48, pp.135-167. 10.1002/fld.820. hal-00021402

HAL Id: hal-00021402

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

Submitted on 17 Jan 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
Numerical simulation of cavitating flow in 2D and 3D
inducer geometries

O. Coutier-Delgosha¹,*,†, R. Fortes-Patella¹,‡, J. L. Reboud²,§,
N. Hakimi³,¶ and C. Hirsch⁴,||

¹LEGI/INPG, BP 53, 38041 Grenoble cedex 9, France

²LEMD/UJF, BP 166, 38042 Grenoble cedex 9, France

³NUMECA, 5 av. F. Roosevelt, 1050 Brussels, Belgium

⁴VRIJE Universiteit Brussel Peinlaan 2, 1050 Brussels, Belgium

SUMMARY

A computational method is proposed to simulate 3D unsteady cavitating flows in spatial turbopump
inducers. It is based on the code FineTurbo, adapted to take into account two-phase flow phenomena.
The initial model is a time-marching algorithm devoted to compressible flow, associated with a low-
speed preconditioner to treat low Mach number flows. The presented work covers the 3D implementation
of a physical model developed in LEGI for several years to simulate 2D unsteady cavitating flows. It is
based on a barotropic state law that relates the fluid density to the pressure variations. A modification
of the preconditioner is proposed to treat efficiently as well highly compressible two-phase flow areas
as weakly compressible single-phase flow conditions.

The numerical model is applied to time-accurate simulations of cavitating flow in spatial turbopump
inducers. The first geometry is a 2D Venturi type section designed to simulate an inducer blade suction
side. Results obtained with this simple test case, including the study of its general cavitating behaviour,
numerical tests, and precise comparisons with previous experimental measurements inside the cavity,
lead to a satisfactory validation of the model. A complete three-dimensional rotating inducer geometry
is then considered, and its quasi-static behaviour in cavitating conditions is investigated. Numerical
results are compared to experimental measurements and visualizations, and a promising agreement is
obtained.

KEY WORDS: cavitation; preconditioner; barotropic model; inducer; two-phase flow

<footnote>*Correspondence to: O. Coutier-Delgosha, ENSAM Lille/LML, 8 Bld. Louis XIV, 59046 Lille cedex, France.</footnote>

<footnote>†E-mail: olivier.coutier@lille.ensam.fr</footnote>

<footnote>‡E-mail: fortes@hmg.inpg.fr</footnote>

<footnote>§E-mail: jean-luc.reboud@grenoble.cnrs.fr</footnote>

<footnote>¶E-mail: nouredine.hakimi@numeca.be</footnote>

<footnote>||E-mail: hirsch@stro.rub.ac.be</footnote>

<footnote>Contract/grant sponsor: Education French Ministry MERT</footnote>

<footnote>Contract/grant sponsor: SNECMA Moteurs, Rocket engine division (DMF)</footnote>

<footnote>Contract/grant sponsor: CNES (Centre National d'Etudes Spatiales)</footnote>
## 0. 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 flows 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 defining 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–4].

* a rotating cavitation pattern, showing different sizes of the cavitation structures in the different blade-to-blade passages of the machine and leading to super- or sub-synchronous instabilities [5].

The difficulty of predicting inducer cavitating behaviour is aggravated by these complex properties of the flow field: three-dimensional, rotating, turbulent, unsteady, weakly compressible in the single-phase areas, highly compressible in two-phase flow conditions, and characterized by continuous phase changes. Previous numerical and theoretical models of cavitating flows 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 flow field 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 effective 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 first case local fluid 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 first applications to pump geometries have been presented by Combes and Archer [15] and Medvitz et al. [16].

To gain further knowledge concerning cavitating flows in inducers, a three-dimensional model of steady and unsteady cavitation for inviscid or viscous fluids 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 final objective to provide accurate simulations of unsteady cavitating flows in the inducers of rocket engine turbopumps. The numerical model results from the integration in the 3D code FineTurbo developed by Numeca Int. [17–19] of a physical cavitation model. This one is based on previous work performed in LEGI, concerning 2D numerical simulations of unsteady cavitating flows [3, 8, 20, 21]. The two-phase aspects of cavitation are treated by introducing a barotropic state law that strongly links the fluid density to the pressure variations. The liquid-vapour mixture is considered as a single fluid 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 difficulty 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–25].

The present paper focuses on the computational method applied to treat efficiently this two-phase flow simulation. The use of the barotropic state law in association with a time-marching algorithm is not straightforward: rapid space and time density fluctuations deteriorate the numerical stability, and the standard low-speed preconditioner is only devoted to incompressible or low-compressible flow. Although the preconditioner is here necessary to treat pure liquid and pure vapour configurations, its efficiency 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 modifications 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 first one was designed so that the flow is subjected to the same pressure field 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–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 first case is to test the numerical model in a simple configuration, to check the accuracy of the results and the influence of the numerical parameters. The second test case is a first attempt to predict the cavitating flow field in a complete 3D geometry, and compare it to experimental results. This validation of the code concerns presently only quasi-static two-phase flow conditions: pronounced unsteady behaviour associated with cavitation is out of the scope of this study. Current work is in progress to simulate efficiently such configurations.

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 flow 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 efficiency. The influence 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 flow field 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.

# 1. THE PHYSICAL MODEL

Generally, two-phase flow models are based on the assumption that the fluid 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 effects, 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 fluxes of mass and momentum between the two phases. The difficulty 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 fluid model based on previous numerical and physical work developed in LEGI [8, 20]. It assumes that only one fluid is considered. This fluid 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 fluid. The void fraction $α=(ρ - ρ_l)/(ρ_v - ρ_l)$ can be defined 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 fluxes between the phases are treated implicitly, and no supplementary assumptions are required.

<figure><img src="image_2.png" /><figcaption>Figure 1. Barotropic state law $ρ(P)$.</figcaption></figure>
Indeed, concerning the momentum fluxes, 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 flow (the friction forces are high, compared to the buoyancy forces). This hypothesis is often assumed for this problem of sheet-cavity flows, 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 fluxes resulting from vapourization and condensation processes are governed by the state law. In the present work, an empirical barotropic law $\rho(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 fluid 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 effects associated with vapourization and condensation are not considered in the present study.

$$ \rho = \frac{P}{r_{pg}T_0} \quad (\text{perfect gas law for the pure vapour}) \quad (1) $$

$$ \rho = \rho_{\text{ref}} \cdot \sqrt[n]{\frac{P + P_0}{P_{\text{ref}} + P_0}} \quad (\text{Tait law for the pure liquid}) \quad (2) $$

where $P_{\text{ref}}$ and $\rho_{\text{ref}}$ are reference pressure and density, and $r_{pg}$ is the constant for the perfect gas. For water, $P_0 = 3 \times 10^8$ and $n = 7$.

The two fluid states are joined smoothly in the vapour-pressure neighbourhood by the following law:

$$ \rho = A_b + B_b \times \sin\left(\frac{1}{B_b} \frac{P - P_{\text{vap}}}{C_{\min}^2}\right) \quad (3) $$

with $A_b$ and $B_b$ constants depending from the connection between the three parts of the state law, at the points $(P_1, \rho_1)$, and $(P_v, \rho_v)$.

$$ A_b = \frac{\rho_1 + \rho_v}{2} \quad \text{and} \quad B_b = \frac{\rho_1 - \rho_v}{2} \quad (4) $$

It results in the evolution law presented in Figure 1, characterized mainly by its maximum slope $1/C_{\min}^2$, where $C_{\min}^2 = (\partial P/\partial \rho)_{\min}$. $C_{\min}$ can thus be interpreted as the minimum speed of sound in the two-phase mixture.

## 2. THE COMPUTATIONAL METHOD

The numerical model used to simulate cavitating flow fields 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 specific numerical treatments were applied to make this integration fully efficient. 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–19]. We present here the main features
of the numerical resolution, with emphasis on the modifications introduced in the code to treat cavitating flow fields.

## 2.1. 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 effects are neglected in the vapourization and condensation phenomena (see Section 1), so temperature does not appear in the state law of the fluid. As a consequence, the energy equation is disconnected from the others: the temperature field has no influence 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:

$$ \iiint_V \frac{\partial \mathbf{U}}{\partial \tau} dV + \iint_S (\mathbf{F}_{\text{inv}} + \mathbf{F}_v) \cdot \mathbf{n} dS = \iiint_V \text{Sce} dV \quad (5) $$

where

$$ \mathbf{U} = (\rho, \rho u, \rho v, \rho w), $$

$\mathbf{F}_{\text{inv}}$ and $\mathbf{F}_v$ are the inviscid and viscous fluxes across the frontier $S$ of $V$ (normal $\mathbf{n}$) $\text{Sce}$ is the source term, and $\tau$ is the numerical time-step.

Convergence is reached when the time derivative terms, which have no physical meaning, vanish.

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

$$ \frac{\partial}{\partial t} \iiint_V \mathbf{U} dV + \iiint_V \frac{\partial \mathbf{U}}{\partial \tau} dV + \iint_S (\mathbf{F}_{\text{inv}} + \mathbf{F}_v) \cdot \mathbf{n} dS = \iiint_V \text{Sce} dV \quad (6) $$

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 flows. In the case of low-compressible or incompressible flows, its efficiency decreases dramatically. As a matter of fact, at low Mach number the convective speed $\mathbf{v}$ of the flow becomes small, compared with the acoustic speed $c$. In other words, the eigenvalues of the system become much different, and the local time-step, which is based on the highest of these eigenvalues, is not adapted to the lowest waves, whose speed is $\mathbf{v}$. This can be seen in the expression for the eigenvalues, expressed hereafter for a propagation direction $\mathbf{n}$:

$$ \begin{aligned} \lambda_{1,2} &= \mathbf{v} \cdot \mathbf{n} \\ \lambda_{3,4} &= \mathbf{v} \cdot \mathbf{n} \pm c |\mathbf{n}| \end{aligned} \quad (7) $$

As a result, the time-marching algorithms converge very slowly in such configuration.
## 2.3. Low-speed preconditioner

This well-known problem has been addressed by many authors and solved by introducing a preconditioner [33–37] in the Navier–Stokes equations. It is based on the modification of the pseudo-time derivative terms in the governing equations. Such modifications have no influence 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 stiffness and enhancing the convergence.

The simplest development consists in adding in the continuity equation a term $(1/\beta_p^2)\partial P/\partial\tau$. 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 $\beta_p$, and by Turkel [34–36] who proposed for a two-dimensional inviscid problem a supplementary term $\alpha_p v/\beta_p^2 \partial P/\partial\tau$ in the momentum equation:

$$ 
\begin{aligned}
& \frac{1}{\rho \beta_p^2} \frac{\partial P}{\partial \tau} + \frac{1}{\rho c^2} \left( u \frac{\partial P}{\partial x} + v \frac{\partial P}{\partial y} \right) + \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \\
& \frac{\alpha_p u}{\rho \beta_p^2} \frac{\partial P}{\partial \tau} + \frac{\partial u}{\partial \tau} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + \frac{1}{\rho} \frac{\partial P}{\partial x} = 0 \\
& \frac{\alpha_p v}{\rho \beta_p^2} \frac{\partial P}{\partial \tau} + \frac{\partial v}{\partial \tau} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + \frac{1}{\rho} \frac{\partial P}{\partial y} = 0
\end{aligned}
\quad (8) $$

The term $(1/c^2)(u \frac{\partial P}{\partial x} + v \frac{\partial P}{\partial y})$ is derived from $u \frac{\partial \rho}{\partial x} + v \frac{\partial \rho}{\partial y}$, considering the particular case of perfect gases. Turkel analysed the influence of the two preconditioning parameters $\alpha_p$ and $\beta_p$ in this case, in order to remove the stiffness 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 $\Gamma^{-1}$ called the preconditioning matrix. The resulting governing equations are solved in the following conservative form:

$$ \iiint_V \Gamma^{-1} \frac{\partial \mathbf{Q}}{\partial \tau} dV + \iint_S (\mathbf{F}_{\text{inv}} + \mathbf{F}_v) \cdot \mathbf{n} \, dS = \iiint_V \mathbf{S} \mathbf{c} e \, dV - \frac{\partial}{\partial t} \iiint_V \mathbf{U} \, dV \quad (9) $$

where $\mathbf{Q}(P_g, u, v, w)$ are the new dependent variables, $P_g$ is the gauge pressure (denoted simply $P$ hereafter), $\mathbf{F}_{\text{inv}}$ and $\mathbf{F}_v$ are the same inviscid and viscous fluxes as previously

$$ \Gamma^{-1} = \begin{bmatrix} \frac{1}{\beta_p^2} & 0 & 0 & 0 \\ \alpha_p \frac{u}{\beta_p^2} & \rho & 0 & 0 \\ \alpha_p \frac{v}{\beta_p^2} & 0 & \rho & 0 \\ \alpha_p \frac{w}{\beta_p^2} & 0 & 0 & \rho \end{bmatrix} $$
Accuracy is highly improved by choosing the gauge pressure as dependent variable. Preconditioning coefficients are, respectively, equal to $\alpha_p = 1$ and $\beta_p = \beta_0 \cdot V_{\text{ref}}$ where $\beta_0$ is a constant whose default value is 3 and $V_{\text{ref}}$ is the reference velocity. As a result of this modification, the acoustic speed is replaced by a new one, which is of the same order of magnitude as $u$. So the stiffness of the system is removed, and the convergence rate is considerably improved in the low speed flow areas [17, 35].

However, the present application of cavitating flows implies the simultaneous treatment of very different flow configurations: according to the barotropic state law presented in Section 1, the fluid is highly compressible in the two-phase flow 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 fluid have to be computed efficiently by the same algorithm, without creating any spurious discontinuity in the flow field. So a supplementary attention was paid to the preconditioning, which is unfavourable in compressible flow configuration, while it is necessary in weakly compressible or incompressible conditions. A modification of the algorithm is proposed hereafter to overcome this limitation.

## 2.4. Eigenvalues

The eigenvalues of the preconditioned system defined 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 $\bar{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

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \frac{1}{\beta_p^2} \frac{\partial P}{\partial \tau} + u \frac{\partial \rho}{\partial x} + \rho \frac{\partial u}{\partial x} + v \frac{\partial \rho}{\partial y} + \rho \frac{\partial v}{\partial y} + w \frac{\partial \rho}{\partial z} + \rho \frac{\partial w}{\partial z} &= 0 \\
\frac{\partial u}{\partial t} + \frac{\alpha_p u}{\rho \cdot \beta_p^2} \frac{\partial P}{\partial \tau} + \frac{\partial u}{\partial \tau} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} + \frac{1}{\rho} \frac{\partial P}{\partial x} &= 0 \\
\frac{\partial v}{\partial t} + \frac{\alpha_p v}{\rho \cdot \beta_p^2} \frac{\partial P}{\partial \tau} + \frac{\partial v}{\partial \tau} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} + \frac{1}{\rho} \frac{\partial P}{\partial y} &= 0 \\
\frac{\partial w}{\partial t} + \frac{\alpha_p w}{\rho \cdot \beta_p^2} \frac{\partial P}{\partial \tau} + \frac{\partial w}{\partial \tau} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z} + \frac{1}{\rho} \frac{\partial P}{\partial z} &= 0
\end{aligned}
\quad (10) $$

Since the density only depends on pressure, the sound speed is defined as $c = \sqrt{\frac{dP}{d\rho}}$.

As a consequence, $\partial\rho/\partial x$, $\partial\rho/\partial y$ can be replaced respectively by $1/c^2 \partial P/\partial x$, $1/c^2 \partial P/\partial y$, and the mass equation becomes

$$ \frac{\partial P}{\partial \tau} + u \frac{\beta^2}{c^2} \frac{\partial P}{\partial x} + \rho \beta^2 \frac{\partial u}{\partial x} + v \frac{\beta^2}{c^2} \frac{\partial P}{\partial y} + \rho \beta^2 \frac{\partial v}{\partial y} + w \frac{\beta^2}{c^2} \frac{\partial P}{\partial z} + \rho \beta^2 \frac{\partial w}{\partial z} = 0 \quad (11) $$
The system can be turned into the following form, involving the non-conservative Jacobians
A, B, C:

$$
\frac{\partial \bar{Q}}{\partial \tau} + A \frac{\partial \bar{Q}}{\partial x} + B \frac{\partial \bar{Q}}{\partial y} + C \frac{\partial \bar{Q}}{\partial z} = 0 \quad (12)
$$

$$
A = \begin{bmatrix}
u \frac{\beta_p^2}{c^2} & \rho \cdot \beta_p^2 & 0 & 0 \\
\frac{1}{\rho} - \frac{\alpha_p u^2}{\rho \cdot c^2} & u(1-\alpha_p) & 0 & 0 \\
-\alpha_p \frac{u \cdot v}{\rho \cdot c^2} & -\alpha_p v & u & 0 \\
-\alpha_p \frac{u \cdot w}{\rho \cdot c^2} & -\alpha_p w & 0 & u
\end{bmatrix},
\quad
B = \begin{bmatrix}
v \frac{\beta_p^2}{c^2} & 0 & \rho \cdot \beta_p^2 & 0 \\
-\alpha_p \frac{u \cdot v}{\rho \cdot c^2} & v & -\alpha_p u & 0 \\
\frac{1}{\rho} - \frac{\alpha_p v^2}{\rho \cdot c^2} & 0 & v(1-\alpha_p) & 0 \\
-\alpha_p \frac{v \cdot w}{\rho \cdot c^2} & 0 & -\alpha_p w & v
\end{bmatrix},
$$

$$
C = \begin{bmatrix}
w \frac{\beta_p^2}{c^2} & 0 & 0 & \rho \cdot \beta_p^2 \\
-\alpha_p \frac{u \cdot w}{\rho \cdot c^2} & w & 0 & -\alpha_p u \\
-\alpha_p \frac{v \cdot w}{\rho \cdot c^2} & 0 & w & -\alpha_p v \\
\frac{1}{\rho} - \frac{\alpha_p w^2}{\rho \cdot c^2} & 0 & 0 & w(1 - \alpha_p)
\end{bmatrix}
$$

The final expression for A, B, C depends on the expression for c, issued from the different
parts of the barotropic state law ρ(P).

$$
\text{Pure liquid: } c = \sqrt{\frac{dP}{d\rho}} = \sqrt{\frac{n}{\rho_{\text{ref}}} (P_{\text{ref}} + P_0)} \left(\frac{P+P_0}{P_{\text{ref}}+P_0}\right)^{(n-1)/2} \quad (13)
$$

$$
\text{Pure vapour: } c = \sqrt{\frac{dP}{d\rho}} = \sqrt{rT_0} \qquad (14)
$$

$$
\text{Two-phase mixture: } c = \sqrt{\frac{dP}{d\rho}} = C_{\min} \left[ \cos \left( \frac{P - P_{\text{vap}}}{B_{\text{b}} \cdot C_{\min}^2} \right) \right]^{-1/2} \quad (15)
$$

As mentioned in the previous section, the preconditioner may be unfavourable in the com-
pressible flow areas, i.e. in the two-phase mixture. So we focus hereafter on this particular part
of the barotropic state law. In this configuration, the expression for the eigenvalues associated
with propagation direction **n**(n<sub>x</sub>, n<sub>y</sub>, n<sub>z</sub>) can be obtained from |An<sub>x</sub> + Bn<sub>y</sub> + Cn<sub>z</sub> - λ| = 0:

$$
\lambda_{1,2} = \mathbf{v} \cdot \mathbf{n} \tag{16}
$$

$$
\lambda_3 = \frac{1}{2} \left[ \mathbf{v} \cdot \mathbf{n} \left( 1 - \alpha_p + \beta_p^2 \frac{M^2 C 2}{B_b u^2} \right) \right. \\
\left. \qquad + \sqrt{ \left( \mathbf{v} \cdot \mathbf{n} \left( 1 - \alpha_p + \beta_p^2 \frac{M^2 C 2}{B_b u^2} \right) \right)^2 + 4 \beta_p^2 \left( \mathbf{n}^2 - (\mathbf{v} \cdot \mathbf{n})^2 \frac{M^2 C 2}{B_b u^2} \right) } \right] \tag{17}
$$

$$
\lambda_4 = \frac{1}{2} \left[ \mathbf{v} \cdot \mathbf{n} \left( 1 - \alpha_p + \beta_p^2 \frac{M^2 C 2}{B_b u^2} \right) \right. \\
\left. - \sqrt{ \left( \mathbf{v} \cdot \mathbf{n} \left( 1 - \alpha_p + \beta_p^2 \frac{M^2 C 2}{B_b u^2} \right) \right)^2 + 4 \beta_p^2 \left( \mathbf{n}^2 - (\mathbf{v} \cdot \mathbf{n})^2 \frac{M^2 C 2}{B_b u^2} \right) } \right]
$$

where

$$
C2 = B_b \times \cos\left(\frac{P - P_{\text{vap}}}{B_b \cdot C_{\text{min}}^2}\right) \quad (18)
$$

and

$$
M = \frac{u}{C_{\min}} \tag{19}
$$

It was clearly shown in Reference [26] that the stiffness of these four eigenvalues, which mainly depends on the term βₚ² M²C2/B_b u², is much increased in such a two-phase flow configuration, compared with weakly compressible flow conditions. Thus, a modification of the preconditioner is suitable in compressible areas, to remove locally its action. Preconditioning parameters must be adapted progressively, to finally reach the following values in the most compressible flow parts (corresponding to α = 0.5):

$$
\begin{align}
\alpha_p &= 0 \notag \\
\beta_p &= \frac{u\sqrt{B_b}}{M\sqrt{C2}} \tag{20}
\end{align}
$$

(i.e. β<sub>p</sub> is the physical sound celerity in highly compressible regions of the flow: the action
of the preconditioner is in this case locally completely removed).

To ensure a smooth evolution of parameters α<sub>p</sub> and β<sub>p</sub>, the following expressions are finally applied in the whole flow field:

$$
\begin{align}
\alpha_p &= \begin{cases} 1 & \text{if } \alpha = 0 \text{ or } \alpha = 1 \\ |2\alpha - 1| & \text{if } 0 < \alpha < 1 \end{cases} \tag{21} \\
\beta_p &= \min \left( \frac{u\sqrt{B_b}}{M\sqrt{C2}}, \beta_0 \cdot V_{\text{ref}} \right) \notag
\end{align}
$$
The expression for $\beta_p$ ensures for all flow conditions the same order of magnitude for the three eigenvalues $\lambda_1$, $\lambda_2$, and $\lambda_3$. In weakly compressible or incompressible configurations, $\beta_p$ is maintained at the constant value $\beta_0 \cdot V_{ref}$, which limits the weight of the term $4\beta_p^2 \cdot \mathbf{n}^2$ in $\lambda_2$ and $\lambda_3$, while in highly compressible regions, $\beta_p$ is reduced down to the physical sound celerity, to control the terms $\beta_p^2 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 modification and avoid spurious computation divergences.

## 2.5. Discretization

Equation (9) is discretized in its conservative form with a finite volume approach. The discrete form of Equation (9) over a computational cell volume becomes

$$V \Gamma^{-1} \frac{\partial \mathbf{Q}}{\partial \tau} + \sum_{\text{Faces}} (\mathbf{F} \cdot \mathbf{n})^* \Delta S = V \cdot \mathbf{S} \mathbf{c} \mathbf{e} - V \frac{\partial \mathbf{U}}{\partial t} \quad (22)$$

where $(\mathbf{F} \cdot \mathbf{n})^*$ is the numerical flux at the cell interfaces. A centred approximation is applied to the viscous fluxes, while the inviscid ones are calculated with a central convection scheme associated with artificial dissipation. The resulting expression for the numerical flux along the i-direction on the right side of the cell is

$$(\mathbf{F} \cdot \mathbf{n})_{i+1/2,j,k}^* = [\mathbf{F}_{\text{inv}} \cdot \mathbf{n}]_{i+1/2,j,k} + [\mathbf{F}_v \cdot \mathbf{n}]_{i+1/2,j,k} - \mathbf{D}_{i+1/2,j,k}^{\text{art}} \quad (23)$$

The first right hand term is the centred approximation of the inviscid flux, and it is treated as

$$[\mathbf{F}_{\text{inv}} \cdot \mathbf{n}]_{i+1/2,j,k} = \frac{1}{2} ([\mathbf{F}_{\text{inv}} \cdot \mathbf{n}]_{i,j,k} + [\mathbf{F}_{\text{inv}} \cdot \mathbf{n}]_{i+1,j,k}) \quad (24)$$

The artificial dissipation $\mathbf{D}^{\text{art}}$ is composed of two terms, respectively of second and fourth order, as initially proposed by Jameson et al. [32]. The following form of $\mathbf{D}^{\text{art}}$ leads to a central second order accurate convection scheme:

$$\mathbf{D}_{i+1/2,j,k}^{\text{art}} = \frac{1}{2} d_{i+1/2,j,k}^2 \Lambda_{i+1/2,j,k}^i \delta \mathbf{Q}_{i+1/2,j,k} + d_{i,j,k}^4 \Lambda_{i+1/2,j,k}^i \delta^2 \mathbf{Q}_{i,j,k} - d_{i+1,j,k}^4 \delta^2 \mathbf{Q}_{i+1,j,k} \quad (25)$$

$$\delta \mathbf{Q}_{i+1/2,j,k} = \mathbf{Q}_{i+1,j,k} - \mathbf{Q}_{i,j,k} \quad (26)$$

$$\delta^2 \mathbf{Q}_{i,j,k} = \mathbf{Q}_{i+1,j,k} - 2\mathbf{Q}_{i,j,k} + \mathbf{Q}_{i-1,j,k} \quad (27)$$

where $\Lambda^i$ is the spectral radius along the i-direction multiplied with the cell face area. $d^2$ and $d^4$ are coefficients associated, respectively, to second order and fourth order artificial dissipations. Their expression is given by

$$d_{i+1/2,j,k}^2 = \varepsilon_p^2 \max(\delta_{i-1,j,k}, \delta_{i,j,k}, \delta_{i+1,j,k}, \delta_{i+2,j,k}) + \varepsilon_p^2 \max(\gamma_{i-1,j,k}, \gamma_{i,j,k}, \gamma_{i+1,j,k}, \gamma_{i+2,j,k}) \quad (28)$$

$$d_{i+1/2,j,k}^4 = \max(0, \varepsilon_p^4 - d_{i+1/2,j,k}^2)$$

The variables $\delta_i$ and $\gamma_i$ are sensors to activate the second-difference dissipation in regions of strong gradients, and to de-activate it elsewhere. $\delta_i$ (proposed initially by Jameson et al. [32])
applies in pressure gradients areas, such as shocks, and $\gamma_i$ is a supplementary sensor added for cavitating flows:

$$ \delta_{i,j,k} = \left| \frac{P_{i+1,j,k} - 2P_{i,j,k} + P_{i-1,j,k}}{P_{i+1,j,k} + 2P_{i,j,k} + P_{i-1,j,k}} \right| \quad \text{and} \quad \gamma_{i,j,k} = \left| \frac{\rho_{i+1,j,k} - 2\rho_{i,j,k} + \rho_{i-1,j,k}}{\rho_{i+1,j,k} + 2\rho_{i,j,k} + \rho_{i-1,j,k}} \right| \quad (29) $$

The standard values for $\epsilon_p^2$ and $\epsilon_p^4$ 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 flow configurations. The standard value for $\epsilon_p^4$ is 0.1. The influence 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 flux 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 difference scheme that ensures a second order accuracy in time:

$$ \left[ \frac{\partial}{\partial t} \iiint_V \mathbf{U} \, d\mathbf{V} \right]^{n+1} = \frac{1.5\mathbf{U}^{n+1} - 2\mathbf{U}^n + 0.5\mathbf{U}^{n-1}}{\Delta t} \mathbf{V} \quad (30) $$

## 2.6. Stability analysis of the numerical scheme

The numerical stability of the resolution has been studied in the three flow configurations defined 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 influence of several physical and numerical parameters was tested, and the efficiency of the preconditioner was analysed in several flow configurations. Nice convergence and stability properties are demonstrated in the incompressible case (which confirms the efficiency of the preconditioner in such configurations), 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 influences 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 flow velocity inside the cells, has the same kind of impact: the configuration $v/u=1$ is unfavourable for convergence and stability, compared with configuration $v/u=0$ (flow in the direction of the mesh).

Concerning the numerical parameters, a diminution of the fourth order artificial dissipation $\epsilon_p^4$ induces an important convergence slow down. So does also an increase of the mesh aspect ratio.

## 2.7. Under-relaxation of the density

The cavitation phenomenon is a very sharp and rapid physical process, which results in two main numerical difficulties:

* 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 fine meshes. To ensure the local conservation in this configuration 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 fluid 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 efficient control of the density time fluctuations 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:

$$ \rho^{i+1} = \rho^i + \eta(\rho(P^{i+1}) - \rho^i) \quad (31) $$

where the standard value of under-relaxation coefficient $\eta$ 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 $\rho(P)$, so that the fluid state law is always respected.

Underrelaxing the density is the same as multiplying the local time step $\Delta\tau$ by $\eta$ in the mass equation. As a matter of fact, the convergence significantly 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.

## 2.8. 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 influence [39]. They are not yet applied in the present application.

## 2.9. Initial transient treatment

In experiments, the classical machinery cavitation test consists in setting initially a relatively high pressure $P_{ref}$ in the flow field, for which no vapour appears. Then, the pressure is decreased until the cavitation number $\sigma = (P_{ref} - P_{vap})/\rho V_{ref}^2$ corresponding to operating conditions is reached.

A similar numerical procedure was developed, but instead of decreasing the reference pressure $P_{ref}$ in the flow field, 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 flow field change because of this continuous boundary condition variation. With the first method, only the two-phase areas are affected by the cavitation number evolution.

## 2.10. 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-\epsilon$ model with extended wall functions.

The Baldwin–Lomax model requires at the walls a first cell located at a non-dimensional distance $y^+$ close to 1, while in the case of the $k-\epsilon$ 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-\epsilon$ turbulence models are in progress to improve physical analyses. An investigation of the respective efficiency of different turbulence models associated with the single fluid modellization of cavitation is also proposed by Coutier-Delgosha et al. in Reference [41].

# 3. CALCULATION IN A VENTURI TYPE SECTION

## 3.1. mesh

Numerical simulations are first 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 fluctuations in its downstream part. Although the considered geometry is much simpler than a real inducer one (3D effects and rotating phenomena are not present), its cavitating behaviour is globally representative of sheet cavitation phenomena observed in such a machine.

The Venturi profile 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 different 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 first cell at the walls.
<figure><img src="image_3.png" /><figcaption>Figure 2. Mesh of the Venturi type section (248×88).</figcaption></figure>

## 3.2. Validation of the global cavitating behaviour

Steady-state calculations have first been performed, and a satisfactory convergence level was obtained in the case of small cavities. For a lower cavitation number, some important numerical difficulties appear, associated with large pressure fluctuations in the Venturi section. As a matter of fact, cavitation is always unsteady at some scale, and the fluctuations 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 fluctuation 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 $\sigma = 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 flow field to cover the length $L_{ref}$ of the profile, 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 configuration the main information provided by this figure is the time evolution of the cavity length, associated with
<figure><img src="image_4.png" /><figcaption>Figure 3. Pressure and density fields°-stable cavity for $\sigma = 0.4 - V_{\text{ref}} = 10.8 \text{ m/s}$.</figcaption></figure>

the maximum void ratio in each section. We also analysed the evolution of the vapour volume (Figure 4(b)) and the inlet pressure fluctuation (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_{\text{ref}}$ (9 m/s and 14.4 m/s) close to the value $V_{\text{ref}} = 10.8 \text{ m/s}$ applied for the calculation.

The length of the cavitation sheet is correctly predicted, down to $\sigma = 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 fluctuating rear part of the liquid/vapour structures.

### 3.3. Numerical tests

The influence of numerical and physical parameters of the model is tested in the previous flow configuration. The objective is to determine the effect 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. Influence of numerical parameters indicates the accuracy level of the numerical scheme, while influence of the physical ones gives an evaluation of the code robustness in a cavitating configuration.

Results of the tests are reported in Table II. As the final cavity obtained is almost stable, it can be characterized by its final volume, and the vapour volume contained inside. So for each case, these two results are indicated. The final volume of the cavity is calculated by considering the volume of the cells where $\alpha > 4\%$.

#### 3.3.1. Influence of numerical parameters.
The influence of the numerical parameters, such as $\Delta t$, $\epsilon_p^2$, $\epsilon_p^4$, $\eta$, CFL, and the mesh size is globally satisfactory:

* A too high time-step $\Delta t$ generates some little oscillations of the cavity. For lower values of $\Delta t$, its influence is acceptable: for $\Delta t = 10^{-4} \text{ s}$, the cavity volume is by 7% smaller
<figure><img src="image_5.png" /><figcaption>Figure 4. Time evolution of the cavity given by numerical calculation: (a) Colors represent the density value: white for the pure liquid, and from yellow to dark red for the two-phase mixture. It illustrates the temporal evolution (in abscissa) of the cavity length (graduated in ordinate); (b) Time evolution of the volume of vapour in the flow field; and (c) Time evolution of the inlet pressure.</figcaption></figure>

than for $\Delta t = 3 \times 10^{-4}$ s. The standard value of the time-step has been chosen so that $\Delta t = T_{\text{ref}}/100$.

* $\epsilon_p^2$ and $\eta$ have no noticeable influence on the result: nevertheless, a too small value of $\eta$($<0.1$) leads to an underestimation of the cavity development. The simulation should be still continued to finally lead to a standard converged cavity size. On the contrary, a
<figure><img src="image_6.png" /><figcaption>Figure 5. Numerical and experimental [3] results $L(\sigma)$.</figcaption></figure>

<table><caption>Table I. Tested parameters.</caption><thead><tr><th>Parameter</th><th>Tested range</th><th>Standard value</th></tr></thead><tbody><tr><td colspan="3"><strong>Numerical parameters</strong></td></tr><tr><td>Mesh size (three dimensions were tested)</td><td>108×36 – 160×60 – 248×88</td><td>160×60</td></tr><tr><td>Physical time step Δt (s)</td><td>10<sup>-4</sup> → 6 × 10<sup>-4</sup></td><td>3 × 10<sup>-4</sup></td></tr><tr><td>Under-relaxation coefficient η</td><td>0.05 → 0.8</td><td>0.2</td></tr><tr><td>Fourth order artificial viscosity coeff. ε<sub>p</sub><sup>4</sup></td><td>0.02 → 1</td><td>0.1</td></tr><tr><td>Second order artificial viscosity coeff. ε<sub>p</sub><sup>2</sup></td><td>0.5 → 3</td><td>1</td></tr><tr><td>CFL</td><td>0.1 → 5</td><td>3</td></tr><tr><td colspan="3"><strong>Physical parameters</strong></td></tr><tr><td>Cavitation number σ</td><td>0.22 (blockage) → 0.46 (little cavity)</td><td>0.38</td></tr><tr><td>Minimum sound celerity C<sub>min</sub></td><td>1.5 m/s → 4 m/s</td><td>2.25</td></tr><tr><td>Ratio ρ<sub>v</sub>/ρ<sub>l</sub></td><td>0.001 → 0.1</td><td>0.001</td></tr><tr><td>Turbulence model</td><td>Baldwin–Lomax/k-ε</td><td>k-ε</td></tr></tbody></table>

too high value (>0.4) results in the computation divergence. If $ε_p^2$ is smaller than 0.5, numerical instabilities also increase up to the calculation divergence.

* $ε_p^4$ has only a little influence on the final shape of the cavity: the mean vapour volume remains almost constant for all the values of this parameter.

* The CFL number has no influence on the result, so far it is small enough to ensure the stability of the numerical scheme.

* The mesh size has almost no influence on the result, so far it is fine enough: the difference between the vapour volumes obtained with the two coarsest meshes (108×38 and 160×60) is significant, but it becomes very weak if we compare the intermediate mesh with the finest one (248×88). The same behaviour is observed concerning the final cavity volume. Thus, the intermediate mesh is found to be fine enough to avoid any mesh effect in the computed flow field.
### 3.3.2. Influence of the physical parameters.
An increase of $C_{\min}$ or $\rho_v/\rho_1$ 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 diffusive effects in the two-phase flow areas. In the case of $\rho_v/\rho_1$ this influence is small. On the contrary, the value of $C_{\min}$ directly influences 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 different cavity lengths. The optimal value was found to be independent of both $L_{\text{ref}}$ and $V_{\text{ref}}$, and is about 2 m/s.

* The turbulence model has only a little influence on the result, in this rather stable cavitation configuration.

* Numerical stability is ensured down to very low values of the cavitation number $\sigma$. For example, in the configuration $\sigma = 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 fluctuations 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.

## 3.4. 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 profiles are available: their position is indicated in Figure 7. The mean values of the void ratio are presented for these five profiles 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]).

# 4. CALCULATION ON A FOUR-BLADE INDUCER

Simultaneously with the previous 2D studies on fixed geometries, we performed some first computations on a rotating four-blade inducer designed at the SNECMA Moteurs Rocket
<table><caption>Table II. Results of the tests: influence of numerical parameters.</caption><thead><tr><th rowspan="2">σ</th><th rowspan="2">C<sub>min</sub> (m/s)</th><th rowspan="2">ρ<sub>v</sub>/ρ<sub>l</sub></th><th rowspan="2">Δt×10<sup>-4</sup> s</th><th colspan="5">Parameters</th><th colspan="2">Results</th></tr><tr><th>Mesh</th><th>η</th><th>CFL</th><th>ε<sub>p</sub><sup>2</sup></th><th>ε<sub>p</sub><sup>4</sup></th><th>Turbulence model</th><th>Volume<sub>vapour</sub> (×10<sup>-4</sup> m<sup>3</sup>)</th><th>Volume<sub>cavity</sub> (×10<sup>-3</sup> m<sup>3</sup>)</th></tr></thead><tbody><tr><td colspan="12" style="text-align: left; font-weight: bold;">Standard values</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.5</td></tr><tr><td colspan="12" style="text-align: left; font-weight: bold;">Numerical parameters</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>248×88</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.5</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.5</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>108×36</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>2.1</td><td>0.9</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.3</td><td>1.4</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>3.9</td><td>1.3</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.5</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>6</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>3.0/10.1</td><td>0.8/3.6</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.0</td><td>1.5</td></tr><tr><td>0.46</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.05</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.6</td><td>1.5</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.15</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.5</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.6</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.4</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.6</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.5</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>DIV</td><td>DIV</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.5</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.7</td><td>1.5</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.2</td><td>1.4</td></tr><tr><td>0.48</td><td>2.25</td><td>0.01</td><td>3</td><td>160×60</td><td>0.2</td><td>3</td><td>1</td><td>0.1</td><td>Baldwin-Lomax</td><td>4.3</td><td>1.4</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6. Mixture density field: blockage due to the volume of vapour for $\sigma = 0.33$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Experimental data sections.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Time-averaged void ratio: num. (-) and exp. (○). Comparison between numerical results (lines) and optical probes measurements (points): $V_{\text{ref}} = 10.8 \, \text{m/s}$-Baldwin-Lomax model-$C_{\min} = 2.25 \, \text{m/s}$. Stable cavity: $L_{\text{cav}} = 80 \, \text{mm}$ (cavity external shape in dotted line from image processing) Ratio 3 between vertical and horizontal scales.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9. Inducer geometry.</figcaption></figure>

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 fluctuating, 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 first 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 effects of the NPSH decrease. So the unsteady effects are not modelled in the computations presented hereafter. For a given pseudo-vapour pressure value, the inlet pressure does not fluctuate so far the convergence is reached, and it is directly used to compute the NPSH values reported on the charts.

## 4.1. Mesh

Only one blade-to-blade channel is considered, and periodicity conditions are applied to its frontiers in azimuthal 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 effects of cavitation are investigated. The mesh is composed of five 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 difficulty to mesh an inducer, compared with a typical centrifugal impeller, is related to the low $\beta_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,
<figure><img src="image_11.png" /><figcaption>Figure 10. Mesh structure: meridional and blade-to-blade representations.</figcaption></figure>

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 first 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
<figure><img src="image_12.png" /><figcaption>Figure 11. Details of the grid: (a) Blade-to-blade representation, near the hub; (b) Blade-to-blade representation, near the shroud; (c) 3D view of the mesh on the hub and on the blade tip.</figcaption></figure>

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 efforts, 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 flow 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.
<figure><img src="image_13.png" /><figcaption>Figure 12. Boundary conditions.</figcaption></figure>

## 4.2. Boundary conditions

The calculations are performed with water flow 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⁷. 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 flow rate.

* A radial equilibrium is fixed 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.

## 4.3. Numerical parameters

$ρ_1 = 1000 \text{ kg/m}^3$ and $ρ_\infty/ρ_1 = 0.01$ (standard value),

CFL number: 1.5,

$η = 0.05$,

$Δt = 10^{-4}$ s (about 1/100e of the inducer rotation period, to only simulate the quasi-static effects),

The $β_p$ preconditioning coefficient is based on the tangential inlet tip velocity,

All other parameters are set to their default value indicated in Section 3.3.

## 4.4. Results: performance charts

Cavitation conditions in the inducer and its performance are usually characterized by the following non-dimensional numbers:

$$ \tau'_{i} = \frac{\tau_{i}}{\tau_{0}} \quad \text{with} \quad \tau_{i} = \frac{P_{\text{lot}}^{\text{inlet}} - P_{\text{vap}}}{\rho_{1} \cdot V_{\text{ref}}^{2}} \quad \text{and} \quad \tau_{0} \text{ a reference value of } \tau_{i} \text{ (cavitation parameter)} $$

$$ \Psi = \frac{P_{\text{outlet}} - P_{\text{inlet}}}{\rho_{1} \cdot V_{\text{ref}}^{2}} \quad (\text{static pressure elevation coefficient}) $$

$V_{\text{ref}}$ is in this configuration the tangential tip velocity, and $P_{\text{lot}}^{\text{inlet}}$, $P_{\text{inlet}}^{\text{inlet}}$, and $P_{\text{outlet}}^{\text{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 configuration 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
<figure><img src="image_14.png" /><figcaption>Figure 13. Experimental and numerical performance charts. Green squares indicate the operating conditions illustrated in Figure 14.</figcaption></figure>

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 $\tau'_i$ corresponding to the performance breakdown is over-predicted by the model: the gap between experimental and numerical results equals $\Delta\tau'_i = 1$.

This discrepancy must be considered circumspectly. In the case of centrifugal pumps [43], the order of magnitude of $\tau'_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 ($\tau'_i$ at the performance breakdown is about 1), the order of magnitude of $\tau'_i$ is in the present case the same as the one of $\Delta\tau'_i$.

So, the over-estimation of $\tau'_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 fluctuations 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 effects. 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 efficient 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><img src="image_15.png" /><figcaption>Figure 14. Evolution of the cavitating areas in the inducer as $\tau'_i$ decreases (the hub is coloured in grey, the blades in red, and the iso-surfaces $\rho = 950 \text{ kg/m}^3$ are represented in yellow).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15. Cavitation in the inducer at two operating points: Comparison between the experimental and numerical results.</figcaption></figure>

### 4.5. Visualizations of the cavitating areas

The evolution of the two-phase structures during the $\tau'_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 first confined on the blade suction side, close to the shroud (view 2). Then, it rapidly fills 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–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 dif-
ferent, the cavitating structures in the inducer are not compared at the same value of $\tau'_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 significant that its precise shape is not the
same in experiment and in the calculation, since this structure is highly fluctuating 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 operat-
ing conditions. It suggests that the 3D numerical model applied is efficient 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 $\tau'$ at the performance drop.
Nevertheless, this first 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.

5. 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 fluid
density evolution in the cavittaing regions. Particularities of this computational method, such
as a modification of the low speed preconditioner and a necessary density fluctuations control,
were detailed. The accuracy and robustness of the numerical model have been tested first 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 influence of physical and numerical parameters was also presented, and reference val-
ues were defined for further calculations. The structure of the cavity was finally compared to
experimental results, and a satisfactory agreement was found concerning the external shape
of the cavity and the void ratio distribution inside. A first attempt of a 3D simulation con-
cerning 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.

NOMENCLATURE

$C_{\min}$ minimum speed of sound in the mixture (m/s)

c local sound celerity (m/s)
<table>
  <tr>
    <td>L<sub>ref</sub></td>
    <td>length of the Venturi type section (m)</td>
  </tr>
  <tr>
    <td>NPSH</td>
    <td>net pressure suction head (m)</td>
  </tr>
  <tr>
    <td>P</td>
    <td>local static pressure (Pa)</td>
  </tr>
  <tr>
    <td>P<sub>l</sub>, P<sub>v</sub></td>
    <td>limit pressures of pure liquid and pure vapour in the barotropic law (Pa)</td>
  </tr>
  <tr>
    <td>P<sub>ref</sub></td>
    <td>reference pressure = inlet pressure (Pa)</td>
  </tr>
  <tr>
    <td>P<sub>tot</sub></td>
    <td>total pressure (Pa)</td>
  </tr>
  <tr>
    <td>P<sub>g</sub></td>
    <td>gauge pressure (Pa)</td>
  </tr>
  <tr>
    <td>P<sub>vap</sub></td>
    <td>vapour pressure (Pa)</td>
  </tr>
  <tr>
    <td>T<sub>ref</sub></td>
    <td>L<sub>ref</sub>/V<sub>ref</sub> = time reference (s)</td>
  </tr>
  <tr>
    <td>v(u, v, w)</td>
    <td>local flow velocity (m/s)</td>
  </tr>
  <tr>
    <td>V</td>
    <td>control volume (m<sup>3</sup>)</td>
  </tr>
  <tr>
    <td>V<sub>ref</sub></td>
    <td>reference velocity (m/s)</td>
  </tr>
  <tr>
    <td>x, y, z</td>
    <td>cartesian co-ordinates (m)</td>
  </tr>
  <tr>
    <td>z, r</td>
    <td>axial and radial co-ordinates in a rotaing machinery (m)</td>
  </tr>
  <tr>
    <td>α</td>
    <td>local void fraction (dimensionless)</td>
  </tr>
  <tr>
    <td>α<sub>p</sub></td>
    <td>preconditioning parameter (dimensionless)</td>
  </tr>
  <tr>
    <td>β<sub>p</sub></td>
    <td>preconditioning parameter = pseudo-sound celerity (m/s)</td>
  </tr>
  <tr>
    <td>β<sub>0</sub></td>
    <td>constant parameter used to calculate β<sub>p</sub> (dimensionless)</td>
  </tr>
  <tr>
    <td>β<sub>b</sub></td>
    <td>geometrical blade angle in a (z/r, θ) frame (rad)</td>
  </tr>
  <tr>
    <td>ε<sub>p</sub><sup>4</sup></td>
    <td>fourth order artificial dissipation coefficient (dimensionless)</td>
  </tr>
  <tr>
    <td>ε<sub>p</sub><sup>2</sup>, ε<sub>p</sub><sup>2</sup></td>
    <td>second order artificial dissipation coefficients (dimensionless)</td>
  </tr>
  <tr>
    <td>η</td>
    <td>under-relaxation coefficient for density (dimensionless)</td>
  </tr>
  <tr>
    <td>μ, ν</td>
    <td>dynamic, cinematic viscosity (Pa s, m<sup>2</sup>/s)</td>
  </tr>
  <tr>
    <td>ρ<sub>v</sub></td>
    <td>vapour density corresponding to pressure P<sub>v</sub> (kg/m<sup>3</sup>)</td>
  </tr>
  <tr>
    <td>ρ<sub>l</sub></td>
    <td>liquid density corresponding to pressure P<sub>l</sub> (kg/m<sup>3</sup>)</td>
  </tr>
  <tr>
    <td>ρ</td>
    <td>αρ<sub>v</sub> + (1 - α)ρ<sub>l</sub> local density of the mixture (kg/m<sup>3</sup>)</td>
  </tr>
  <tr>
    <td>θ</td>
    <td>revolution angle in a rotaing machinery (rad)</td>
  </tr>
  <tr>
    <td>σ</td>
    <td>(P<sub>ref</sub> - P<sub>vap</sub>)/ρV<sub>ref</sub><sup>2</sup> cavitation number (dimensionless)</td>
  </tr>
  <tr>
    <td>τ</td>
    <td>pseudo time-step (s)</td>
  </tr>
  <tr>
    <td>τ<sub>i</sub></td>
    <td>(P<sub>inlet</sub><sup>tot</sup> - P<sub>vap</sub>)/ρ<sub>l</sub> · V<sub>ref</sub><sup>2</sup> cavitation parameter used for an inducer (dimensionless)</td>
  </tr>
  <tr>
    <td>τ<sub>i0</sub></td>
    <td>reference cavitation parameter (dimensionless)</td>
  </tr>
  <tr>
    <td>τ'<sub>i</sub></td>
    <td>τ<sub>i</sub>/τ<sub>i0</sub> (dimensionless)</td>
  </tr>
  <tr>
    <td></td>
    <td>(P<sub>outlet</sub><sup>tot</sup> - P<sub>inlet</sub>) / ρ · V<sub>ref</sub><sup>2</sup> static pressure elevation coefficient (dimensionless)</td>
  </tr>
</table>

ACKNOWLEDGEMENTS

This research was supported by a doctoral grant from the Education French Ministry MERT and by SNECMA Moteurs DMF. The authors wish also to express their gratitude to CNES (Centre National d'Etudes Spatiales, France) for the continuous support.

REFERENCES

1. Kubota A, Kato H, Yamaguchi H, Maeda H. Unsteady structure measurement of cloud cavitation on a foil section using conditional sampling technique. *Journal of Fluids Engineering* 1989; **111**:204–210.

2. Le Q, Franc JP, Michel JM. Partial cavities: global behaviour and mean pressure distribution. *Journal of Fluids Engineering* 1993; **115**:243–248.

3. Reboud JL, Stutz B, Coutier O. Two phase flow structure of cavitation: experiment and modelling of unsteady effects. *Proceedings of the 3rd International Symposium on Cavitation*, Grenoble, France, 1998.
4. Reboud JL, Fortes-Patella R, Hofmann M, Lohrberg H, Ludwig G, Stoffel B. Numerical and experimental investigations on the self-oscillating behaviour of cloud cavitation. *Proceedings of the 3rd ASME/JSME Joint Fluids Engineering Conference*, San Francisco, 1999.

5. de Bernardi J. Experimental analysis of instabilities related to cavitation in a turbopump inducer. *Proceedings of the 1st International Symposium on Pump Noise and Vibrations*, Paris, France, 1993.

6. Kueny JL, Schultz JL, Desclaux J. Numerical prediction of partial cavitation in pumps and inducers. *Proceedings of IAHR Symposium*, Trondheim, 1988.

7. Von Kaenel, Maitre T, Rebattet C, Kueny JL, Morel P. Three dimensional partial cavitating flow in a rocket turbopump inducer: numerical predictions compared with laser velocimetry measurements. *Proceedings of Cav'95 International Symposium*, Deauville, France, 1995.

8. Delannoy Y. Modélisation d'écoulements instationnaires et cavitants. Ph.D. Thesis, INPG, Grenoble, 1989.

9. Knapp RT, Daily JT, Hammit FG. *Cavitation*. McGraw Hill: New York, 1970.

10. Plesset MS, Chapman RS. Collapse of initially spherical vapour cavity in the neighbourhood of a solid boundary. *Journal of Fluid Mechanics* 1971; **47**:283-290.

11. Takasugi N, Kato H, Yamagushi H. Study on cavitating flow around a finite span hydrofoil. *Proceedings of the Cavitation and Multiphase Flow Forum*, ASME-FED, vol. 153. 1993; 177-182.

12. Alajbevoic A, Grogger H, Philipp H. Calculation of transient cavitation in nozzle using the two-fluid model. *Proceedings of the 12th Annual Conference on Liquid Atomization and Spray Systems*, Indianapolis, 1999.

13. Kunz R, Boger D, Chyczewski T, Stinebring D, Gibeling H. Multi-phase CFD analysis of natural and ventilated cavitation about submerged bodies. *Proceedings of the 3rd ASME/JSME Joint Fluids Engineering Conference*, San Francisco, 1999.

14. Bunnell RA, Heister SD. Three-dimensional unsteady simulation of cavitating flows in injector passages. *Journal of Fluids Engineering* 2000; **122**:791-797.

15. Combes JF, Archer A. Etude de la cavitation dans la pompe SHF. *Proceedings of the Colloque Machines Hydrauliques: instationnarités et effets associés*, Société Hydrotechnique de France, Chatou, France, 2000.

16. Medvitz RB, Kunz RF, Boger DA, Lindau JW, Yocum AM, Pauley LL. Performance analysis of cavitating flow in centrifugal pumps using multiphase CFD. *Proceedings of ASME FEDSM 2001*, New Orleans, U.S.A., 2001.

17. Hakimi N. Preconditioning methods for time dependent Navier-Stokes equations. Ph.D. Thesis, Vrije Universiteit Brussels, Belgium, 1997.

18. Numeca International. User guide of FineTurbo, 2002.

19. Zhu ZW. Multigrid operation and analysis for complex aerodynamics. Ph.D. Thesis, Vrije Universiteit Brussel, Belgium, 1996.

20. Coutier-Delgosha O, Reboud JL, Delannoy Y. Numerical simulation of unsteady cavitating flow. *International Journal for Numerical Methods in Fluids* 2003; **42**(5):527-548.

21. Delannoy Y, Kueny J-L. Two phase flow approach in unsteady cavitation modelling. *Proceedings of the 1990 ASME Cavitation and Multiphase Flow Forum* 1990; **98**:153-158.

22. Merkle CL, Feng J, Buelow PEO. Computational modelling of the dynamics of sheet cavitation. *Proceedings of the 3rd International Symposium on Cavitation*, Grenoble, France, 1998.

23. Song C, He J. Numerical simulation of cavitating flows by single-phase flow approach. In *Proceedings of the 3rd International Symposium on Cavitation*, Michel J-M, Kato H (eds), vol. 2, Grenoble, France, 1998; 295-300.

24. van der Heul DR, Vuik C, Wesseling P. Efficient computation of flow with cavitation by compressible pressure correction. In *Proceedings of the European Congress on Computational Methods in Applied Science and Engineering, ECCOMAS 2000*, Onate E, Bugeda G, Suarez B (eds), Barcelona, Spain, 2000.

25. Ventikos Y, Tzabiras G. A numerical method for the simulation of steady and unsteady cavitating flows. *Computers and Fluids* 2000; **29**:63-88.

26. Coutier-Delgosha O, Fortes-Patella R, Reboud JL, Hakimi N, Hirsch C. Stability of preconditioned Navier-Stokes equations associated with a cavitation model. *Computers and Fluids*, 2004, to be published.

27. Stutz B. Analyse de la structure diphasique et instationnaire des poches de cavitation. Ph.D. Thesis, Institut National Polytechnique de Grenoble, France, 1996.

28. Stutz B, Reboud JL. Experiments on unsteady cavitation. *Experiments in Fluids* 1997; **23**:191-198.

29. Stutz B, Reboud JL. Two-phase flow structure of sheet cavitation. *Physics of Fluids* 1997; **9**(12):3678-3686.

30. Bestion D. The physical closure laws in the CATHARE code. *Nuclear Engineering and Design*, vol. 124. North Holland: Amsterdam, 1990.

31. Hirsch C. *Numerical Computation of Internal and External Flows*. Wiley: New York, 1990.

32. Jameson A, Schmidt W, Turkel E. Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes. *AIAA Paper*, vol. 81. 1991; 1259.

33. Choi D, Merkle CL. The application of preconditioning in viscous flows. *Journal of Computational Physics* 1993; **105**:207-223.

34. Turkel E. Preconditioning methods for solving the incompressible and low speed compressible equations. *Journal of Computational Physics* 1987; **72**:277-298.
35. Turkel E, Fiterman A, van Leer B. Preconditioning and the limit of the incompressible flow equations. *Journal of Applied Numerical Mathematics*, 1994.

36. Turkel E, Radespiel R, Kroll N. Assessment of preconditioning methods for multidimensional aerodynamics. *Computers and Fluids* 1997; **26**:613–634.

37. Weiss J, Maruszewski J, Smith W. Implicit solution of preconditioned Navier–Stokes equations using algebraic multigrid. *AIAA Journal* 1999; **37**:29–36.

38. Chorin AJ. A numerical method for solving incompressible viscous flow problems. *Journal of Computational Physics* 1967; **2**:12–26.

39. Longatte F. Contribution à l'Analyse Phénoménologique des Ecoulements Stationnaires dans les Turbomachines: Etude du Couplage Pompe-Circuit et Rotor-Stator. *Ph.D. Thesis*, Institut National Polytechnique de Grenoble, France, 1998.

40. Hakimi N, Pierret S, Hirsch C. Presentation and application of a new k-eps model with wall functions. *Proceedings of Eccomas*, 2000.

41. Coutier-Delgosha O, Fortes-Patella R, Reboud JL. Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation. *Journal of Fluids Engineering* 2003; **125**:38–45.

42. Kueny JL, Reboud JL, Desclaux J. Analysis of partial cavitation: image processing and numerical prediction. *Proceedings of the Cavitation'91*, ASME-FED, vol. 116. 1991; 55–60.

43. Hoffman M, Stoffel B, Coutier-Delgosha O, Fortes-Patella R, Reboud JL. Experimental and numerical studies on a centrifugal pump with 2D curved blades in cavitating conditions. *Journal of Fluids Engineering* 2003; **125**:970–978.

44. Pouffary B. Modélisation numérique des écoulements dans une roue centrifuge: analyse des mécanismes associés à la chute de performances en cavitation. *Master Report*, E.C. Lyon, 2001.

45. Coutier-Delgosha O, Fortes-Patella R, Reboud JL. Numerical simulation of turbopump inducer cavitating behaviour. *International Journal of Rotating Machinery*, 2004, to be published.

46. Hirsch C, Hakimi N. Preconditioning methods for time-marching Navier-Stokes solvers. *Solution Techniques for Large-Scale CFD Problems* (CMAS). Wiley: New York, 1995; 333.

47. Kubota A, Kato H, Yamaguchi H. A new modelling of cavitating flows: a numerical study of unsteady cavitation on a hydrofoil section. *Journal of Fluid Mechanics* 1992; **240**:59–96.