# Experiments and modelling of cavitating flows in Venturi. Part I: stable cavitation

Julien Rolland, Guillaume Boitel, Stéphane Barre, Eric Goncalvès da Silva,
Regiane Fortes Patella

► To cite this version:

Julien Rolland, Guillaume Boitel, Stéphane Barre, Eric Goncalvès da Silva, Regiane Fortes Patella. Experiments and modelling of cavitating flows in Venturi. Part I: stable cavitation. CAV 2006 Symposium, Sep 2006, Wageningen, New Zealand. hal-00212021

HAL Id: hal-00212021

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

Submitted on 31 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
# EXPERIMENTS AND MODELLING OF CAVITATING FLOWS IN VENTURI
## PART I: STABLE CAVITATION

Julien ROLLAND - CNES, France : julien.rolland@hmg.inpg.fr

Guillaume BOITEL - LEGI - INPG, France

Stéphane BARRE - LEGI - INPG, France

Eric GONCALVES - LEGI - INPG, France

Regiane FORTES PATELLA - LEGI - INPG, France

## ABSTRACT

Correlated experimental and numerical studies were carried out to analyse cavitating flows in Venturi geometries. New double optical probe measurements were performed and special data processing methods were developed to estimate void ratio and velocity fields for cold water flows.

By applying a computational method previously developed in Laboratory LEGI based on the code Fine™/Turbo and on a barotropic approach (Delannoy and Kueny, 1990), several steady calculations were performed in cold water cavitating flows. Local and global analyses were proposed based on comparisons between experimental and numerical results.

The physical and numerical models were improved to take into account thermodynamic effects in cavitation phenomenon. A first application of the model is carried out for a Venturi geometry in the case of R114 flows. Numerical results were compared to experimental ones obtained previously by (Fruman et al., 1999).

## INTRODUCTION

The present work follows previous experimental and numerical studies carried out by the Turbomachinery and Cavitation team of LEGI (Grenoble, France). The purpose is a better understanding and modelling of cavitation phenomena (Reboud et al., 1998, 2003, Pouffary et al., 2003).

The studies are led in collaboration with the French space agency (CNES), the rocket engine division of Snecma, and Numeca International. The main final aims are the analyses of cavitating flows in the rocket engine turbopump inducers, where the run fluids are cryogenic fluids, as LH2 and LOx.

In this paper, we present experimental and numerical studies performed in the case of cold water and thermosensible fluid R114 ($C_2Cl_2F_4$). In order to facilitate and to carry out local measurements and analyses, cavitating flows inside of space turbopump inducers are simulated by Venturi geometries.

Section 1 presents the experimental device used for recent tests in cold water. Section 2 describes the new data processing methods developed to evaluate, from

measurements by double optical probes, the void ratio and velocity fields in cavitating flows.

The experimental approach has been applied for two Venturi geometries. In the first one, analysed in the present paper (named “Venturi 4°”), the cavitating flow is characterized by a quasi-steady behaviour. In the second one, named “Venturi 8°” and studied in the Part II of the work (Fortes Patella et al., 2006), we observe a quasi-periodic unsteady behaviour.

To gain further knowledge concerning cavitating flows, in parallel with experimental studies, two-dimensional and three-dimensional models of steady and unsteady cavitation have been developed in LEGI. One of the numerical works results from the integration of a cavitation model in the 3D Navier-Stokes code Fine™/Turbo developed by Numeca Int.. The cavitation model is represented by a barotropic state law, proposed by (Delannoy and Kueny, 1990) that strongly links the fluid density to the pressure variations.

The numerical model, presented in Section 3, has been applied to simulate water cavitating flow in “Venturi 4°” geometry. Global and local analyses based on the comparison between experimental and numerical results are presented in Section 4.

Moreover, the numerical model has been improved to take into account thermodynamic effects and to simulate thermosensible fluid cavitation. As a matter of fact, because of the latent heat of vaporization, heat transfer may strongly affect the cavitation behaviour of refrigerants and cryogenic fluids. Vapour production through cavitation extracts heat from surrounding liquid and produces a local cooling effect, which delays the cavitation phenomenon.

First studies about thermal effects were carried out by (Stahl 1956, Moore 1968), which defined and used the B-factor theory to characterize the sensitivity of fluids to thermodynamic effects. This theory was employed to provide semi-empirical correlations based on experimental data (Hord 1974). A similar theory, named the “entrainment theory” (Holl 1975) defined and correlated the temperature depression to flow parameters.

Numerous researchers have studied thermodynamic effects of cavitation: (Cooper 1975) presented a
barotropic law to simulate liquid hydrogen pump; (Rapposelli et al. 2003) also used a barotropic law to simulate warm water. The equilibrium model based on thermodynamic properties of each phases was used by (Ventikos 1995). More recently models (Sauer 2000, Saito et al. 2003, Hosangadi 2005, Utturkar et al. 2005), by using a liquid or vapour continuity equation with special attention to source terms, give results in good agreement with the experimental data obtained by (Hord 1974).

In this context, a new model including thermal effects in cavitation phenomenon has been implemented in Fine™/Turbo code. Section 5 presents first numerical results obtained for the “Venturi 4°” geometry in the case of R114 flows. Comparisons with experimental data given by (Fruman et al. 1999) are presented and analysed.

## NOMENCLATURE

c Venturi chord (m)

c<sub>min</sub> minimum speed of sound in the mixture (m.s<sup>-1</sup>)

c<sub>p</sub> specific heat at constant pressure (J.kg<sup>-1</sup>.K<sup>-1</sup>)

L latent heat of vaporization (J.kg<sup>-1</sup>)

T<sub>inlet</sub> inlet temperature (K)

U<sub>inlet</sub> inlet velocity (m.s<sup>-1</sup>)

$$ \alpha \quad \text{void ratio} \quad \alpha = \frac{V_v}{V} = \frac{(\rho - \rho_L)}{(\rho_v - \rho_L)} \quad (-) $$

$$ \rho \quad \text{mixture density} \quad \rho = \alpha \rho_v + (1-\alpha) \rho_L \quad (\text{kg.m}^{-3}) $$

$$ B = \frac{V_v}{V_l} = \frac{\rho_l c_{pl} T}{(\rho_v L)^2} \quad \Delta p_v = \frac{T_\infty - T_v(x)}{\rho_v L} \quad \text{B-factor} \quad (-) $$

$$ \sigma_{\text{inlet}} = \frac{p_{\text{inlet}} - p_{\text{vap}}}{\frac{1}{2} \rho_{\text{inlet}} U_{\text{inlet}}^2} \quad \text{Cavitation number} \quad (-) $$

$$ C_p = \frac{p - p_{\text{inlet}}}{\frac{1}{2} \rho_{\text{inlet}} U_{\text{inlet}}^2} \quad \text{Pressure coefficient} \quad (-) $$

# 1. EXPERIMENTAL DEVICE

The experimental device includes mainly a test loop and a sensor (a double optical probe) associated with an acquisition module.

## 1.1 The test loop

The Venturi type test section of the CREMHYG (INPG Grenoble) cavitation small tunnel was dimensioned and designed to simulate cavitating flows developing on the blades of space turbopump inducers. The hydraulic system (Figure 1) is composed of a circulating pump and of a free surface tank, used to impose the reference pressure in the circuit and to resorb dissolved gases. The flow rate, controlled by a computer, is measured by means of an electromagnetic flowmeter. The pressure within the flow is measured by two sensors located in the free surface tank and in the entry section of the Venturi (S<sub>E</sub>). The adjustment of the pressure in the circuit is obtained by managing the air pressure over the free surface of the tank by means of a vacuum pump. The measurement accuracies are evaluated to be:

$$ \Delta Q/Q = \pm 0,25\% \text{ for flow rate} $$

$$ \Delta P = \pm 0,05 \text{ bar for the pressure.} $$

The Venturi test section consists in parallel sidewalls generating rectangular cross sections. The bottom wall can be equipped with several interchangeable Venturi profiles allowing the study of cavitation on various geometries. These profiles are equipped to receive pressure sensors, temperature sensors, or double optical probe.

<figure><img src="image_2.png" /><figcaption>Figure N°1: Schematic view of the hydraulic set-up</figcaption></figure>
In the present study, we use a profile with a convergence angle of 4.3° and a divergence angle of 4°, illustrated in Figure 2. The edge forming the throat of the Venturi is used to fix the separation point of the cavitation cavity. This profile is characterized by the following geometrical data:

Inlet section: $S_E = 50 \times 44 \text{ mm}^2$ (where the reference pressure is measured);

Throat section: $S_{throat} = 43,7 \times 44 \text{ mm}^2$;

Length of the test section (chord): $L_{ref} = 252 \text{ mm}$.

This profile is equipped with five probing holes to take various measurements such as the local void ratio, instantaneous local speed and pressure. Their horizontal positions $X_1$ from the throat of the Venturi are:

$$X_1 = 5.1 \text{ mm}; X_2 = 20.9 \text{ mm}; X_3 = 38.4 \text{ mm}; X_4 = 55.8 \text{ mm}; X_5 = 73.9 \text{ mm}.$$

<figure><img src="image_3.png" /><figcaption>Figure N°2: Schematic view of the Venturi profile</figcaption></figure>

A photograph of this test section is shown in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure N°3: Venturi test section</figcaption></figure>

The selected operation point of the hydraulic system is the same used by (Stutz 1996, 1997 and 2003). It is characterized by the following physical parameters:

$U_{inlet} = 10.8 \text{ m/s}$: inlet speed

$Q = 0,02375 \text{ m}^3\text{/s}$: flow imposed in the circuit by the circulating pump

$P_{tank} = 0,713 \text{ bar}$: pressure measured in the tank

$\sigma_{inlet} = 0.547 \pm 0.005$: cavitation number in inlet section

$N_{pump} = 506 \text{ rpm}$: rotating speed of the circulating pump.

With these parameters, we obtain a cavity length $L$ ranging between $70 \text{ mm} \le L \le 85 \text{ mm}$, and having a relatively stable aspect, as shown on the Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure N°4: Photograph of the cavity</figcaption></figure>

### 1.b. The sensor

To evaluate void ratio and velocity fields inside of the cavity, we have used a double optical probe. The double optical probe is an intrusive sensor composed by two mono-probes functioning in a independent way but mechanically connected one to other. Each probe consists in an optical fiber ending in a frayed Sapphire point. This point is plunged in the flow. An infra-red radiation is injected at an end of fiber. The part of this signal that is reflected at the tip of the probe is converted into electric signal by an optoelectronic module.

The luminous ray is absorbed by the fluid when the tip of the probe is surrounded by liquid; in this case, the measured tension is minimal. When the tip is in a pure vapour zone, the measured tension is maximal. The signal analysis allows determining the phase (liquid or gas) of the flow structures around the probe. This information gives access to the local time fraction of the vapour phase (void rate) and the convection speed of the vapour structures. We use an optical probe made of two optical fibers of 80 µm diameter and with the frayed ends. The inter-tips distance is $1,05 \text{ mm} \pm 0,02$ (see Figure 5).

<figure><img src="image_6.png" /><figcaption>Figure N°5: Double optical probe</figcaption></figure>

The double probe is introduced into the cavitation sheet using a threaded interface assembled in the boreholes of the profile (also called "stations"). A graduated sliding gauge allows to measure the vertical displacement of the probe with an accuracy of $\Delta Y = \pm 0,05 \text{ mm}$. The optical signals obtained by the double probe are transformed into electric signals by an
optoelectronic module. The values of the gains and the zeros of the apparatus are tuned in order to obtain an exit signals of ± 1 V amplitude. Acquisitions are carried out on a fast A/D data acquisition board NI-DAQ PCI-6110E and a microcomputer. For each acquisition we record 6 million points per way (one way by optical probe), with a sampling rate of 100 Khz (either a time of 60 seconds acquisition per measurement). For each Y position of the optical probe in the Venturi (between ten and fifteen per station), five successive measurements are taken. It generates 30 million points of measurement per position in the cavity corresponding to an observation time of 300s.

## 2. DATA PROCESSING METHODS

### 2.a. Estimation of the local void ratio

The local void ratio $\alpha$ of a vapour/liquid mixture is defined here as the ratio between the cumulated attendance time of the vapour phase $\sum_{i=1}^{n} T_{vap\_i}$ and a given time of observation $T_{tot\_obs}$ in an operation point of the flow. It is written in the form:

$$ \alpha = \frac{\sum_{i=1}^{n} T_{vap\_i}}{T_{tot\_obs}} $$

The post processing algorithm enables us to process the signal from the optical probe, by distinguishing the vapour phase from the liquid phase according to the value of the measured tension. The maximum tension value $V_{max}$ corresponds to the vapour phase and the minimal value $V_{min}$ to the liquid phase (Figure 6). In order to be able to estimate the local void ratio, we must fix a threshold in tension $V_{threshold}$. It makes it possible to determine the phase of the fluid around the probe: liquid state if $V < V_{threshold}$ and in a vapour state if $V > V_{threshold}$. This threshold is fixed, in all this study, according to the parameter:

$$ \beta = \frac{V_{threshold} - V_{min}}{V_{max} - V_{min}} $$

The value of this parameter for this kind of flow was gauged by (Stutz 1996, 1997a, 1997b, 2000, 2003). They retained a value of $\beta= 0.1$ to give correct void ratios for this typical flow pattern. In the present work, we have used this value for the determination of the void ratio.

### 2.b. Estimation of instantaneous local velocity

The used optical probe, aligned in the direction of the flow, produces two distinct signals (S1 upstream and S2 downstream). It is possible to deduce the value of the tangential component (parallel with the profile) of the instantaneous local speed by considering the temporal shift between the trace of a bubble passing on the first probe and that of this same bubble passing on the second probe. For this, it is necessary to detect the vapour bubbles. The electric signals of the two optical probes are

thus analysed and broken up by introducing two tension thresholds (Figure 7). The low threshold $V_{low}$, taken on the rising part of the electric signal, leads us to detect the beginning of the crossing of a vapour bubble. The high threshold $V_{high}$ is considered to detect the end of the crossing of the bubble on the downward part of the signal. A bubble is taken into account in the analysis only when these two thresholds are crossed consecutively.

<figure><img src="image_7.png" /><figcaption>Figure N°6: Relative positions of the thresholds for void ratio analysis</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure N°7: Bubble localization with threshold method</figcaption></figure>

In the present study, the values of the low and high thresholds are given by considering the parameters:

$$ \beta_{low} = \frac{V_{low} - V_{min}}{V_{max} - V_{min}} \quad \text{and} \quad \beta_{high} = \frac{V_{high} - V_{min}}{V_{max} - V_{min}} $$

where $\beta_{low} = 0.2$ and $\beta_{high} = 0.8$ (Stutz 1996, 1997a, 1997b, 2003). An analysis of the $\beta$ value influence on the velocity evaluation has been performed. Some results are presented in Part II of this work (Fortes Patella et al., 2006).

After bubble detection by the both probes, we proceed to the evaluation of the instantaneous velocities. The algorithm used for obtaining these data, as well as the several statistical procedures of treatment used to evaluate the corresponding Probability Density Functions (PDF) are described hereafter.
First of all, the calculation algorithm reviews all the bubbles detected on the probe S1 (first probe in the direction of the flow). For each bubble seen on S1, a temporal analysis interval will be explored on the signal of the S2 probe. If *Ideb1* is the index corresponding to the beginning of a bubble considered on the S1 signal, then the analysis interval of the S2 signal will correspond to *Ideb1-Itau* < *t* < *Ideb1+Itau*, where *Itau* is the amplitude of the analysis interval. (Stutz 1996, 1997a, 1997b) proposed to determine *Itau* by using the value of the peak of the function of temporal correlation between the signals resulting from S1 and S2. In the present study, we rather choose an empirically constant evaluated from the various examinations obtained. This value is fixed here at *Itau* = 1000, which corresponds to a temporal analysis interval of $\Delta T = \pm Itau*Tech = \pm 10ms$, where *Tech* is the sampling period of the optical probe signals (*Tech*=1/*Fech*=10⁻⁵s) and *Fech* is the sampling frequency.

Thus, we seek on the S2 probe signal all the existing bubbles in the analysis interval defined. For each one of these bubbles, a function test is calculated as follows:

$$ \epsilon = \frac{t_1 - t_2}{t_1} \quad \text{where } t_1 \text{ is the duration of the} $$

bubble considered on S1 and $t_2$ the duration of the bubble found on S2 in the analysis interval. From the evaluation of all parameters $\epsilon$ obtained for the bubbles found on S2, the minimum value $\epsilon_{min}$ is determined. The bubble corresponding to the value of $\epsilon_{min}$ must now be validated or not according to the value of this test. In this study, the validation criterion corresponds to $\epsilon_{min}$ lower or equal to $\epsilon_{max} = 0.1$, which was empirically determined by a convergence test and a parametric study of the influence of $\epsilon_{max}$ on the velocity fields obtained.

By applying this criterion, we rejected approximately 60 % of the detected bubbles. If the bubble is validated on S2, we compute its temporal shift with the reference bubble of S1. We then obtain *Itauopt* = *Ideb2-Ideb1* where *Ideb1* and *Ideb2* are, respectively, the indices of beginning of bubbles on the signals of the two probes S1 and S2. From *Itauopt*, we obtain an instantaneous speed

$$ u = \frac{\Delta x}{Itauopt.Tech} \quad \text{where } \Delta x \text{ is the space} $$

difference between the two optical probes ($\Delta x$=1.05mm) and *Tech* is the sampling period of the optical probe signals.

By repeating this operation on all the detected bubbles on the S1 probe, we obtain an instantaneous velocity values population, which can be treated in a statistical way to study the dynamics of the flow. For the present study, we obtain typically between 5000 to 70000 values of instantaneous speeds per spatial position of the probes, and this, for an observation time of 300s. It is worth noting that higher is the void ratio, higher is the number of studied bubbles and better is the statistic quality. Following this process we can obtain, for each measurement point in the flow, the *Itauopt* and velocity PDF. An example of a statistical distribution obtained is given on the Figures 8 (*Itauopt*) and 9 (*velocity*). This

example corresponds to the point Y =1mm of the station N°2.

<figure><img src="image_9.png" /><figcaption>Figure N°8: PDF of *Itauopt* values initially obtained (station 2, Y=1mm)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure N°9: PDF of velocity values initially evaluated (station 2, Y=1mm)</figcaption></figure>

We notice an important amount of the *Itauopt* population around the peak representing the most probable value for the considered measurement point (Figure 8). Moreover, we note, for high absolute values of *Itauopt*, that the populations do not decrease and tend to reach an asymptotic value. From Figure 9, we can see, near velocity zero value, a regrouping of the occurrences corresponding to the two asymptotic branches of the *Itauopt* distribution. This effect of overpopulation is accentuated by the fact that speed being inversely proportional to *Itauopt*, so all the high values (in absolute value) for *Itauopt* correspond to a very small speed range around zero value. This population corresponds in fact to a statistical skew of the algorithm used to evaluate velocity values.

Indeed, when we seek on the S2 signal a bubble whose duration is almost equivalent to the basic bubble located on S1, the probability of finding a solution for the problem increases with the *Itau* amplitude of the analysis interval. As a matter of fact, if the signals of S1 and S2 were perfectly random, the velocity Probability Density Function would be a Poisson law centered across the zero
value speed. This law, which is superimposed on the PDF, must be eliminated. Its contribution on the PDF of the real flow can be removed by using a “neutral” *Itauopt* PDF, which considers only "background noise". It is determined thanks to the PDF asymptotic parts. Figure 10 illustrates an example of the PDF obtained after Poisson law correction.

<figure><img src="image_11.png" /><figcaption>Figure N° 10: *Itauopt* PDF obtained by subtraction of random distribution (Poisson law correction)</figcaption></figure>

In spite of this correction, we observe (Figure 11) numerous populations for relatively high speeds absolute values.

<figure><img src="image_12.png" /><figcaption>Figure N° 11: Velocity PDF obtained by subtraction of random distribution (Poisson law correction)</figcaption></figure>

This effect is related to the inaccuracy of the method for high speed. Indeed, because *Itauopt* is an integer number, the measurement method does not give an infinity of possible values for velocity. The accuracy on the velocity evaluation much depends on the *Itauopt* value. For better visualizing and quantifying this effect, the relative error function

$$ \phi = \frac{U_{Itauopt+1} - U_{Itauopt}}{U_{Itauopt}} \quad \text{was calculated for} $$

each value of *Itauopt* in the range corresponding to the considered flow. These results are presented in Figure 12

where the function φ is plotted for values of U between 1 and 20 m/s.

<figure><img src="image_13.png" /><figcaption>Figure N°12: Relative error function φ in the studied velocity range</figcaption></figure>

The accuracy is about 1% in the neighbourhood of 1 m/s, but it increases to reach 16% towards 20 m/s. Moreover, we note for U≥10m/s that the possible speed distribution is represented by few discrete solutions, which are separated each other by several m/s.

In fact, the velocity solutions are gathered according to these few possible values, what artificially increases their population. To correct this effect, we distribute the samples at high speed in a more continuous way by proceeding to an oversampling of the PDF. Several assumptions are formulated for this purpose:

* From |U| ≥17m/s we make a speed grid having a 1m/s step, which corresponds roughly to the obtained accuracy near to 10m/s.

* Between two successive points, we make the assumption that the PDF evolves linearly.

By considering these assumptions, we apply another correction algorithm to the PDF. It consists in generating, between two existing points, fictitious velocity values by respecting a 1m/s step. After that, we distribute all the events corresponding to these two basic points in a linear way between all the fictitious points. For this, we keep the same population ratio between the two basic points before and after correction. Thus, thanks to this correction, we obtain a new oversampled PDF. It should be noted that the form of this new function depends on the selected accuracy (here, 1m/s).

After these corrections, we obtain typical velocity PDF such as one presented in Figure 13.

By applying this treatment to measurements corresponding to the five stations of the experiment, we have evaluated the velocity PDF for all the probed points. From these data, we are able to describe the velocity field inside the cavity. For each point, we keep the most probable value of the PDF to define the local velocity.

Experimental results obtained for flow velocity and void ratio fields have been confronted with ones evaluated by numerical simulation. Main results and analyses are presented here above.
<figure><img src="image_14.png" /><figcaption>Figure N°13: Final Velocity PDF</figcaption></figure>

# 3. NUMERICAL CODE

## 3.a Multiphase Equation System

Calculations have been done with the Fine/Turbo™ code. It is a three-dimensional structured mesh code that solves the time dependant Reynolds-averaged Navier-Stokes equations. A detailed description of the initial code is given in (Hakimi 1998).

The governing equations are written for a homogeneous fluid. This fluid is characterized by a density that varies in the computational domain according to a state law and that is related to the void ratio defined

$$ \alpha = \frac{\rho - \rho_L}{\rho_V - \rho_L} $$

by the void ratio characterizes the volume of vapour in each cell: $\alpha=1$ means that the cell is completely occupied by vapour; inversely, a complete liquid cell is represented by $\alpha=0$. Liquid and vapour phases are characterized by their thermodynamic characteristics. On each cell, the unknowns are calculated by averaging them by the volume occupied.

In this model, the fluxes between the two phases are implicitly treated, and no supplementary assumptions are required. The two phases are considered to be locally (in each cell) in dynamic (no drift velocity) and in thermal (T<sub>V</sub>=T<sub>L</sub>) equilibrium.

The compressible Reynolds-Averaged Navier-Stokes equations are expressed as:

$$ \frac{\partial \rho_m}{\partial t} + \nabla \cdot (\rho_m \mathbf{u}_m) = 0 \\ \frac{\partial \rho_m \mathbf{u}_m}{\partial t} + \nabla \cdot (\rho_m \mathbf{u}_m \otimes \mathbf{u}_m) = -\nabla \cdot (p_m) + \nabla \cdot (\bar{\tau}_m) + \rho_m \mathbf{F}_m \\ \frac{\partial \rho_m E_m}{\partial t} + \nabla \cdot (\rho_m E_m \mathbf{u}_m) = -\nabla \cdot [\mathbf{q}_m] - \nabla \cdot [p_m \mathbf{u}_m] + \nabla \cdot [\bar{\tau}_m \mathbf{u}_m] + \rho_m \mathbf{F}_m \cdot \mathbf{u}_m $$

where $\bar{\tau}_m$ is the shear stress tensor, $\mathbf{q}$ is the heat flux vector, and $\rho_m$ is the mixture density, defined as

$$ \rho_m = \alpha \rho_v + (1-\alpha)\rho_L $$

The total energy is defined as:

$$ E_m = \frac{u^2}{2} + e_m \text{ where } e_m = c_m(T-T_{ref}) + e_m ref $$

$$ \text{with} \quad c_m = \frac{\alpha \rho_v c_v + (1-\alpha) \rho_L c_L}{\alpha \rho_v + (1-\alpha) \rho_L} $$

$$ \text{and} \quad e_m \text{ ref} = \frac{\alpha \rho_v e_v \text{ ref} + (1-\alpha) \rho_L e_L \text{ ref}}{\alpha \rho_v + (1-\alpha) \rho_L} $$

Each pure phase is considered incompressible.

The space discretization is based on a cell-centered finite-volume approach. The numerical fluxes are computed with the central scheme stabilized by the Jameson dissipation (Jameson et al., 1981).

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 explicit four stage Runge-Kutta time stepping procedure is used to advance the solution to steady state.

A complete description of the numerical scheme is presented by (Coutier-Delgosha et al, 2005a).

## 3.b. Preconditioned Navier-Stokes Equations

In the case of low-compressible or incompressible flows, the time-marching algorithm converges very slowly and the used of a low Mach number preconditioner in the Navier-Stokes equations is required (Turkel 1987, Hakimi 1998). It is based on the modification of the pseudo-time derivative terms in the governing equations. Such modifications have no influence on the converged results, 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.

Using this preconditioner the set of equations becomes:

$$ \iiint_V \bar{\Gamma}^{-1} \frac{\partial \mathbf{P}}{\partial \tau} dV + \frac{\partial}{\partial t} \iiint_V \mathbf{w} dV + \iiint_{\Sigma} \mathbf{F} \cdot \mathbf{n} d\Sigma = \iiint_V \mathbf{S} dV $$

introducing the preconditioning matrix $\bar{\Gamma}^{-1}$ and associated variables vector $\mathbf{P}$:

$$ \mathbf{P} = \begin{pmatrix} p_g \\ \mathbf{u} \\ E_g \end{pmatrix} \quad \text{and} \quad \bar{\Gamma} = \begin{pmatrix} \frac{1}{\beta^2} & 0 & 0 \\ \frac{(1+\alpha)\mathbf{u}}{\beta^2} & \rho & 0 \\ \frac{\alpha \mathbf{u}^2 + E_g}{\beta^2} & 0 & \rho \end{pmatrix} $$
$p_g$ and $E_g$ are respectively the gauge (relative) pressure and the gauge total energy, $\mathbf{u}$ is the velocity vector, $\alpha$ and $\beta$ are the preconditioning parameters.

The eigenvalues of the preconditioned system become:

$$ \lambda_{1, 2, 3} = \mathbf{u} \cdot \mathbf{n} \quad \text{et} $$

$$ \lambda_{4, 5} = \frac{1}{2} \left[ \mathbf{u} \cdot \mathbf{n}(1-\alpha) \pm \sqrt{(\mathbf{u} \cdot \mathbf{n}(1-\alpha))^2 + 4\mathbf{n}^2 \beta^2} \right] $$

where $\mathbf{n}$ is the normal vector to the elementary surface dS.

For more details concerning the used preconditionner, see (Coutier-Delgosha et al., 2005b).

### 3.c. Turbulence Models

In the present work, two different turbulence models have been applied : the Spalart-Allmaras model (Spalart & Allmaras 1992) and the Yang-Shih k-epsilon model (Yang & Shih 1993) with extended wall functions (Hakimi et al. 2000).

### 3.d. Barotropic Model

To model cavitation phenomenon and to enclosure the governing equations system, a barotropic state law introduced by (Delannoy 1990) has been implemented in Fine/TurboTM code (Pouffary 2003, Coutier-Delgosha et al. 2005a, 2005b). The fluid density (and so the void fraction) is controlled by a law $\rho(p)$ that links explicitly the mixture fluid density to the local static pressure as represented by Figure 14.

This law is mainly controlled by its maximum slope, which is related to the minimum speed of sound $c_{min}$ in the mixture. The parameter Amin controls this slope:

$$ AMIN^2 = \frac{\rho_L - \rho_V}{2} c_{min}^2 $$

<figure><img src="image_15.png" /><figcaption>Figure N°14: Influence of Amin parameter on barotropic law.</figcaption></figure>

Table 1 presents Amin values used for cold water and R114. For cases with R114, the same physical model is used. The thermodynamic effect is taken into account by a barotropic law given as a function of the local temperature.

<table><caption>Table 1: Physical parameters and tested Amin values</caption><thead><tr><th></th><th>T<br>(K)</th><th>ρ<sub>v</sub><br>kg/m<sup>3</sup></th><th>ρ<sub>ℓ</sub><br>kg/m<sup>3</sup></th><th>Amin<br>m/s.(kg/m<sup>3</sup>)<sup>0.5</sup></th><th>c<sub>min</sub><br>m/s</th></tr></thead><tbody><tr><td rowspan="3">water</td><td rowspan="3">293</td><td rowspan="3">0,0173</td><td rowspan="3">998</td><td>70</td><td>3.13</td></tr><tr><td>50</td><td>2.24</td></tr><tr><td>20</td><td>0.89</td></tr><tr><td rowspan="2">R114</td><td rowspan="2">293</td><td rowspan="2">13,53</td><td rowspan="2">1471</td><td>70</td><td>2.59</td></tr><tr><td>50</td><td>1.85</td></tr></tbody></table>

## 4. EXPERIMENTS-SIMULATION CONFRONTATION IN COLD WATER

### 4.a. Calculation Conditions

The numerical model is firstly compared to experimental data in cold water. The studied case, presented in Section 1, is a steady sheet with a length of about 80 mm, at $U_{inlet}=10.8$ m/s. The reference density is the liquid density equal to 1000 kg/m$^3$.

Steady calculations are led with the k-ε turbulence model, with extended wall functions. The 2D mesh, illustrated in Figure 15, contains 9861 nodes, 173 in the flow direction and 57 in the normal direction.

The y+ values vary from 18 to 50 in cold water. The boundary conditions are mass flow for inlet condition and static pressure imposed at outlet.

A schematic view of the Venturi geometry and of the reference and measurements stations is shown in Figure 16. Table 2 reminds the probing stations coordinates.

The defined sheet length used for representing the computational results has been determined by the length of the iso-line corresponding to a void ratio $\alpha=0.3$. This criterion is argued by the fact that it corresponds to the area of maximal value of the density gradient (see Figure 17). As a matter of fact, because of the strong density gradient at the enclosure of the sheet, the sheet length varies weakly as a function of the interface criterion. A test showed that, for a void ratio between 10% and 30%, the sheet length variation was about 2%.

<figure><img src="image_16.png" /><figcaption>Figure N°15: View of mesh for k-ε computational</figcaption></figure>

<table><caption>Table N°2: Station coordinates</caption><thead><tr><th>Station<br>1</th><th>Station 2</th><th>Station 3</th><th>Station 4</th><th>Station 5</th></tr></thead><tbody><tr><td>5.1mm</td><td>20.9mm</td><td>38.4mm</td><td>55.8mm</td><td>73.9mm</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>Figure N° 16: Venturi throat section and the five probing stations</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure N° 17: Magnitude of density gradient (0 to 10⁶ kg/m⁴) and iso-line α=0.3</figcaption></figure>

### 4.b. Global Analyses

Initially, we have analysed the cavity global behaviour, mainly the sheet length as a function of the cavitation number.

<figure><img src="image_19.png" /><figcaption>Figure N°18: Density field inside the cavity for Amin=50</figcaption></figure>

Three computations with different Amin values have been compared to experimental data. They have been conducted by trial and error to obtain almost the same sheet length. A smooth variation of the inlet cavitation number is observed. The closest result to experimental data is given using the smallest value Amin=20 (see Table N°3).

The influence of Amin on the sheet vapour volume seems weak, except for high values of Amin, where the volume becomes larger. It may be explained by the fact that for large Amin, the upper part of the sheet interface is more diffuse, as shown in Figure 18. Local analysis, presented here above, has confirmed these tendencies.

<table><caption>Table N°3: Comparison of experimental and simulation values of the sheet length (L) versus the cavitation parameter σ</caption><thead><tr><th>Amin</th><th>70</th><th>50</th><th>20</th><th>exp</th></tr></thead><tbody><tr><td>Sigma (inlet)</td><td>0.683</td><td>0.642</td><td>0.600</td><td>0.547</td></tr><tr><td>L<sub>sheet</sub> (mm)</td><td>79.8</td><td>80.2</td><td>81.5</td><td>80</td></tr><tr><td>Vapour volume (*10-4 m³)</td><td>4.005</td><td>3.741</td><td>3.699</td><td></td></tr></tbody></table>

### 4.b. Local Analyses

Complementary analyses concern local void ratio and velocity profiles comparisons inside the cavity. The void ratio and velocity profiles are obtained for five stations.

Figure 19 shows the evolutions of the longitudinal velocity for the present experiments and computations concerning three different Amin values. The overall agreement seems good between the experimental data and the Amin=20 simulation. This is especially true for stations N°1 and 2 where no re-entrant jet phenomena occurs (see Figures 19a and 19b). At these two positions, the velocity profiles are close to one observed in a classical turbulent boundary layer without any detachment or recirculation. Further downstream, for stations N°3, 4 and 5 (see Figures 19c, 19d and 19e) experimental observation indicates a recirculating behaviour with a re-entrant jet extending roughly half the sheet thickness. According to experiments, this flow configuration seems unstable and is smoothly time fluctuating. The steady computation technique used in this paper is not able to reproduce this situation. This result is important because it shows that we must use an unsteady simulation technique to describe correctly this kind of cavitating flow.

Figure 20 illustrates experimental and numerical results concerning void ratio. For “station 1”, closed to the throat, the vaporization phenomenon is mainly represented. In cold water, this is a relatively strong effect and void ratio value is almost equal to “1” near to the wall. As observed previously, the stiffness of the barotropic law (Amin=20) leads to a better interface representation, as shown in Figures 20a.

Downstream, at “station 2”, the void ratio is still high. The void ratio distribution is similar to one obtained for “station 1”, with an increase of the sheet thickness (see Figure 20b).

From "station 3", the re-entrant jet becomes noticeable, as observed before by the velocity field analyses. Then, the void ratio values are overestimated by computations. It is worth noting that, for high values of Amin, the simulated cavitation sheet appears thicker than experimental one.

At stations 4 and 5, the void ratio is highly overestimated by the numerical computation mainly due to the unsteady character of the sheet. Nevertheless, although this overestimation of void ratio inside the sheet for all values of Amin, the thickness of the sheet is quite well represented by the calculation with Amin= 20 (see Figures 20d and 20e).
These computations in cold water allow us to better evaluate the physical and numerical model. In steady situations (stations 1 and 2), results are in quite good agreements with experimental data. A stiffer barotropic law shows better results, firstly about the cavitation number, and secondly about the void ratio and velocity profiles inside the sheet. A better representation of re-entrant jet is obtained by unsteady calculations, as presented in Part II of this work.

## 5. EXPERIMENTS-SIMULATION CONFRONTATION IN THERMOSENSIBLE FLUID R114

To better understand the impact of thermal effects, a first study about cavitation in thermosensible fluid R114 has been also carried out.

In the case of thermosensible fluids, the cavitation phenomenon lead to a local cooling, which can be determined by a characteristic temperature depression. It

represents the decrease of temperature when the heat required to vaporize a volume of vapour is directly furnished by an identical volume of liquid :

$$ \Delta T_{car} = \frac{\rho_v L}{\rho_l C p_l} $$

As cryogenic fluids work close to the critical point, the ratio of liquid to vapour density is lower than for cold water and consequently more liquid mass has to vaporize to sustain an identical cavity.

The fluid thermal sensitivity is characterized also by the reaction to this cooling. For thermosensible fluids, near to the critical point, the slope of the saturation curve p=f(T) is more significant than for cold water (see Figure N°21 and Table N°4).

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

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

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

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

<figure><img src="image_24.png" /><figcaption>Figure N°19: Velocity profiles: Comparison between experimental and simulation results
a) station 1; b) station 2; c) station 3; d) station 4; e) station 5</figcaption></figure>
<figure><img src="image_25.png" /></figure>

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

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

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

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

<figure><img src="image_30.png" /><figcaption>Figure N°20: Void ratio profiles: Comparison between experimental and simulation results
a) station 1; b) station 2; c) station 3;
d) station 4; e) station 5</figcaption></figure>

<table><caption>Figure N°21: Thermodynamic phase diagram for water</caption><thead><tr><th></th><th>T<br>(K)</th><th>p<sub>v</sub><br>bar</th><th>ρ<sub>l</sub><br>kg/m<sup>3</sup></th><th>ρ<sub>v</sub><br>kg/m<sup>3</sup></th><th>ρ<sub>l</sub>/ρ<sub>v</sub></th></tr></thead><tbody><tr><td>water</td><td>293</td><td>0,023</td><td>998</td><td>0,0173</td><td>58000</td></tr><tr><td>R114</td><td>293</td><td>1,811</td><td>1471</td><td>13,53</td><td>108,7</td></tr><tr><td>LH2</td><td>22</td><td>1,73</td><td>68,4</td><td>2,18</td><td>31,4</td></tr></tbody></table>

<table><caption>Table N°4: Thermodynamic characteristics for three typical fluids</caption><thead><tr><th rowspan="2"></th><th rowspan="2">c<sub>p,l</sub><br>kJ/kg/K</th><th rowspan="2">L<br>kJ/kg</th><th colspan="2">ρ<sub>v</sub>L/ρ<sub>l</sub>c<sub>p,l</sub></th></tr><tr><th>K</th><th>Pa/K</th></tr></thead><tbody><tr><td>water</td><td>4,182</td><td>2454,3</td><td>0,01</td><td>146</td></tr><tr><td>R114</td><td>1,00</td><td>130,3</td><td>1,20</td><td>6070</td></tr><tr><td>LH2</td><td>11,0</td><td>433</td><td>1,25</td><td>44130</td></tr></tbody></table>
### 5.a. Calculation Conditions

The studied case is a steady sheet with a length of about 80 mm, at $U_{inlet}=14.4$ m/s and for a Freon reference temperature of $T_{inlet}=20^\circ C$. The geometry is similar to one studied with cold water. Measurement stations are localized at the same place, nevertheless only the first three have been used. The experimental database presented in this confrontation was obtained by (Fruman et al. 1999). The reference density is the liquid density, 1475 kg/m³. Calculations are led with the k-ε turbulence model, with extended wall functions and with the Spalart-Allmaras (SA) model. The 2D mesh is composed of 201x85 nodes for k-ε computations and of 173x121 nodes for SA computations. The y+ values vary from 25 to 95 with the k-ε model. For the Spalart-Allmaras model, y+ values vary from 0.25 to 6. Further work is needed to analyse the mesh convergence.

In a similar way of the water analyses, the sheet length has been evaluated from computations by using the void ratio criteria $\alpha=0.19$. It corresponds to the area of maximum value of the density gradient.

### 5.b. Global Analyses

Table 5 presents numerical results obtained by using different turbulence models, Amin values and the barotropic approach presented before. Results are compared with one obtained by numerical simulations without considering thermal effects.

The influence of turbulence modelling seems to be smooth; results are similar with SA and extended wall function k-ε model at the same Amin value. As in cold water, the cavitation number decreases as the Amin parameter for a given sheet length. Unfortunately, due to numerical instabilities, we are not able to test lower values of Amin yet.

Previous works (Merle 1994, Fruman et al. 1999) concluded to the universality of the law “sheet length versus cavitation number $\sigma_{cav}$” when $\sigma_{cav}$ is defined as

$$ \sigma_{cav} = \frac{p_{ref} - p_{vap}(T_{cav})}{\frac{1}{2} \rho_{ref} U_{ref}^2} $$

In cold water, where thermal effect is negligible, $\sigma_{cav} = \sigma_{ref}$. However, for a thermosensible fluid, we can write:

$$ \sigma_{cav} - \sigma_{ref} = \frac{p_{ref} - p_{vap}(T)}{\frac{1}{2} \rho_{ref} U_{ref}^2} - \frac{p_{ref} - p_{vap}(T_{ref})}{\frac{1}{2} \rho_{ref} U_{ref}^2} = \frac{p_{vap}(T_{ref}) - p_{vap}(T)}{\frac{1}{2} \rho_{ref} U_{ref}^2} $$

For the present study, with a temperature deficit of $\Delta T \sim 1.75$ K, we obtain

$$ \sigma_{cav} - \sigma_{ref} \approx 0.06 $$

This result illustrates that thermodynamic effect seems important for this cavitating flow, even for $\Delta T$ as weak as a few degrees.

<table><caption>Table N°5: Main characteristics and results of the computed test cases</caption><thead><tr><th>turbulence model</th><th>SA</th><th>SA</th><th>k-ε</th><th>k-ε without thermal model</th></tr></thead><tbody><tr><td>Amin</td><td>70</td><td>50</td><td>70</td><td>70</td></tr><tr><td>sigma ref</td><td>0.605</td><td>0.592</td><td>0.612</td><td>0.625</td></tr><tr><td>ΔT=T<sub>ref</sub>-T<sub>cav</sub></td><td>1.76</td><td>1.63</td><td>1.79</td><td></td></tr><tr><td>sigma cav</td><td>0.670</td><td>0.653</td><td>0.678</td><td></td></tr><tr><td>L sheet (mm)</td><td>82.9</td><td>80</td><td>81.4</td><td>83</td></tr></tbody></table>

### 5.c. Local Analyses

Due to thermal effects, the vaporization at the inlet of the cavitation sheet is much smooth. Figure 22 shows the void ratio evolution in the present sheet for the three probed stations. A comparison between experimental and computational results is made. The experimental void ratio at "station 1" is about 40% instead of 90% for cold water. The maximum value of numerical void ratio reaches about 50%. Comparisons between computations with and without thermodynamic modelling show a good description of the thermal effect by our model. Especially for "station 1", the void ratio goes from about 90% (without model) to 50% (with thermodynamic modelling). The overestimation of void ratio on the top of the sheet can be explained as in cold water by the slope of the barotropic law. It can be also shown that thermal effects thickened the sheet, increasing the diffusion at external boundaries.

Simulations overestimate the peak value of the void ratio for "station 1". At the contrary, the results are in good agreement at "station 2" for the peak value (concerning sheet thickness, it will be necessary to decrease Amin in order to sharpen the external edge of the sheet and reduce its thickness). At "station 3", the situation is different and simulations underestimate the peak void ratio value. These effects (for the three probed stations) may be linked to the fact that the applied barotropic model leads to instantaneous vaporization or condensation when the pressure tends towards the vapour pressure $P_v$. Indeed, a time delay, imposed by the physical characteristics of the studied fluid, exists for both vaporization and condensation. Thus, at "station 1", the experimental void ratio is half the computed one, probably because of the instantaneous barotropic model. At station 2, experimental and numerical peak values are comparable. The analysis of results concerning station 3 is more complex because of possible unsteadiness due to the re-entrant jet. In this case, the present stationary computation is unable to correctly describe the flow behaviour.
<figure><img src="image_31.png" /><figcaption>Figure N°22: Void ratio profiles: Comparison between experimental and simulation results</figcaption></figure>

a) station 1; b) station 2; c) station 3

The wall pressure profiles are plotted in Figure 23. For the model without taking into account thermal effects, the pressure remains at an almost constant value, which depends on the cavitation number. When thermodynamic effects are considered, the pressure is reduced at the leading edge because of the decrease of the local vaporization pressure due to temperature deficit. Downstream, local pressure and temperature values are similar to ones obtained without thermal model.

<figure><img src="image_32.png" /><figcaption>Fig N°23: Variation of wall pressure along the cavitation sheet axis</figcaption></figure>

The deficit temperature profile $\Delta T=T_{ref}-T$ calculated for "station 1" is shown in Figure N°24. The largest temperature deficit is approximately 1.75K in the sheet, which is about the magnitude found in previous work by (Fruman et al. 1999) from wall temperature measurements.

It is worth noting that, at the solid wall, the numerical results present a singular behaviour that is attributed to standard wall functions applied. Further work is required to better analyse and to improve calculations at the solid boundaries.

<figure><img src="image_33.png" /><figcaption>Figure N°24: Temperature profile at station 1</figcaption></figure>

## CONCLUSION

From double probe measurements and numerical calculations, we have analysed the global and the local behaviour of a cavitating flow through a Venturi geometry. For a quasi-steady cavitation sheet with a length of about 80 mm, we have evaluated void ratio and velocity fields for cold water and R114 cavitation.

Experimental measurements and data treatment have been improved for cold water tests. The new method proposed for data analysis leads to a better evaluation of local behaviour of steady and unsteady cavitation.
The applied barotropic model, if associated with a stiff slope (corresponding to $c_{min} \sim 1$m/s for cold water), seems to well predict local behaviour in the steady areas of cavitation. Although the global behaviour of the sheet seems steady, experimental data show a re-entrant jet which creates small cloud shedding. In the zones influenced by the re-entrant jet, unsteady calculations are required to better simulate cavitation behaviour. This topic is discussed in the second part of the work.

Concerning cavitation thermodynamic effects, we have implemented in Fine™/Turbo code a new model based on a barotropic approach. First results obtained for the considered Venturi geometry appears to be in good qualitative agreement with experimental and theoretical analyses. Nevertheless, the physical model has to be improved to take into account time delay during vaporization/condensation phenomena and to better predict cavitation behaviour. Moreover, further experimental works concerning local measurements in thermosensible fluid cavitation are capital to allow model calibration and validation.

## ACKNOWLEDGMENTS

The authors wish to express their gratitude to the French space agency CNES and to the SNECMA company to support this research. The authors wish also to express their gratitude to NUMECA International for its cooperation to the development of the numerical code.

## REFERENCES

Cooper P., "Analysis of single and two-phase flow in turbopump inducers", J. Eng. Power, 89 pp 577-588, 1967.

Coutier-Delgosha O., Fortes-Patella R., Reboud J.L., Hakimi N., Hirsch C., "Numerical simulation of cavitating flow in 2D and 3D inducer geometries", Int. Journal for Numerical Methods in Fluids, Vol. 48, Issue 2, pp. 135-167, 2005a.

Coutier-Delgosha O., Fortes-Patella R., Reboud J.L., Hakimi N., Hirsch C., "Stability of preconditioned Navier-Stokes equations associated with a cavitation model", *Computers & Fluids*, Vol. 34, Issue 3, pp. 319-349, March 2005b.

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

Fortes Patella R., Barre S., Reboud J-L, "Experiments and modelling of cavitating flows in Venturi, Part II: unsteady cavitation", *Sixth International Symposium on Cavitation, CAV2006*, Wageningen, The Netherlands, September 2006.

Fruman D. H., Reboud J-L., Stutz B., "Estimation of thermal effects in cavitation of thermosensible liquids", Int J of Heat and Mass Transfer, 42, 3195-3204, 1999.

Hakimi N., "Preconditioning methods for time dependent Navier-Stokes equations", Ph.D.Thesis, Vrije Univ., Brussels, 1998.

Hakimi N., Hirsch, C., Pierret S., "Presentation and application of a new extended $k-\epsilon$ model with wall functions", ECCOMAS, Barcelona, September 2000.

Holl J. W., Billet M.L., Weir D. S., "Thermodynamic effects on developed Cavitation", Journal of Fluids Engineering, 507-514, 1975.

Hord J., "Cavitation in liquid Cryogens", Vol 4, combined correlations for Venturi, Hydrofoil, Ogives and Pumps. NASA CR-2448, 1974.

Hosangadi A., Ahuja V., "Numerical study of cavitation in cryogenic fluids", ASME Fluids Engineering Divison Summer Meeting and Exhibition, Houston, June 2005.

Jameson A., Schmidt W., Turkel E., "Numerical solution of the Euler equations by finite-volume methods using Runge-Kutta time stepping schemes", AIAA Paper 81-1259, June 1981.

Merle L., "Etude expérimentale et modèle physique d'un écoulement cavitant avec effet thermodynamique", Ph. D. Thesis of Institut National Polytechnique de Grenoble, 1994.

Moore R.D., Ruggeri R.S., "Prediction of thermodynamic effects on developed cavitation based on liquid hydrogen and freon 114 data in scaled venturis", NASA TN D-4899, November 1968.

Pouffary B., Fortes-Patella R., Reboud J.L., "Numerical simulation of cavitating flow around a 2D hydrofoil: a barotropic approach", 5th Int. Symp. on Cavitation, Osaka, Japan, 2003.

Rapposelli E., d'Agostino L., "A barotropic cavitation model with thermodynamic effects", Fifth International Symposium on Cavitation, Osaka, Japan, November 1-4, 2003.

Reboud J-L., Stutz B. and Coutier O., "Two-phase flow structure of cavitation : experiment and modelling of unsteady effects", 3rd Int. Symp. on Cavitation, Grenoble, France, April 1998.

Reboud J-L., Coutier-Delgosha O., Pouffary B., Fortes-Patella R., "Numerical simulations of unsteady cavitating flows: some applications and open problems", Proceedings of CAV 2003 Symposium, Conférence Invitée, Osaka. November 2003.

Saito Y., Nakamori I., Ikohagi T., "Numerical analysis of unsteady vaporous cavitating flow around a hydrofoil", Fifth International Symposium on Cavitation, Osaka, Japan, November 1-4, 2003.

Sauer J., Schnerr G.H., "Unsteady cavitating flow - A new cavitation model based on modified front capturing method and bubble dynamics", Proceedings of FEDSM'00, Fluids Engineering Summer Conference, Boston, June 2000.

Stahl H.A, Stepanoff A.J., "Thermodynamic aspects of cavitation in centrifugal pumps", Transaction ASME Journal of Basic Engineering, Vol 78, p1691, 1956.

Spalart P.R., Allmaras S.R., "A one equation turbulence model for aerodynamic flows". AIAA paper 92-0439, 1992.

Stutz B., "Analyse de la structure diphasique et instationnaire de poches de cavitation", Ph. D. Thesis, INPG, CREMHYG Laboratory, 1996.
Stutz B., Reboud J.L., "Two-phase flow structure of sheet cavitation", Phys. Fluids vol 9 (12) pp 3678-3686, 1997a.

Stutz B., Reboud J.L., "Experiments one unsteady cavitation", Exp. In Fluids, Flight. 22, pp 191-198, 1997b.

Stutz B., Reboud J.L., "Measurements within unsteady cavitation", Experiments in Fluids, n°29, pp 545-552, 2000

Stutz B., "Influence of Roughness one the Two-Phase Flow Structure of Sheet Cavitation", J of Fluids Eng., vol. 125, pp 652-659, 2003.

Turkel E, "Preconditioning methods for solving the incompressible and low speed compressible equations", Journal of Comp. Phys., vol 72, pp.277-298, 1987.

Utturkar Y., Thakur S., Shyy W., "Computational modeling of thermodynamic effects in cryogenic cavitation", 43rd AIAA Aerospace sciences meeting an exhibit, Reno, Nevada, 10-13 Jan 2005.

Utturkar Y., Wu J., Wang G., Shyy W., "Recent progress in modelling of cryogenic cavitation for liquid rocket propulsion", Progress in aerospace sciences 41 pp 558-608, 2005.

Ventikos Y., Tzabiras G., "A numerical study of the steady and unsteady cavitation phenomenon around hydrofoils", International Symposium on Cavitation, May 1995, Deauville, France, 1995.

Yang Z., Shih T.H., "A k-e model for turbulence and transitional boundary layer", Near-wall Turbulent Flows, R.M.C. So., C.G. Speziale and B.E. Launder (Editors), Elsevier-Science Publishers B. V., pp165-175, 1993.