# Simulation of unsteady cavitation with a two-equations turbulence model including compressibility effects

Olivier Coutier-Delghosa, Regiane . Fortes Patella, Jean-Luc Reboud

► To cite this version:

Olivier Coutier-Delghosa, Regiane . Fortes Patella, Jean-Luc Reboud. Simulation of unsteady cavitation with a two-equations turbulence model including compressibility effects. Journal of Turbulence, 2002, 3, 10.1088/1468-5248/3/1/058 . hal-00211194

HAL Id: hal-00211194

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

Submitted on 3 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
# Simulation of unsteady cavitation with a two-equation turbulence model including compressibility effects†

O Coutier-Delgosha¹, R Fortes-Patella² and
J L Reboud³

¹ ENSTA UME/DFA, Chemin de la Hunière, 91761 Palaiseau cedex, France

² LEGI-INPG, BP 53, 38041 Grenoble Cedex 9, France

³ ENISE-LTDS, 58, rue Jean Parot, 42023 St Etienne, France

E-mail: coutier@enstay.ensta.fr, regiane.fortes@inpg.fr and
reboud@enise.fr

**Abstract.** Unsteady effects associated with cavitation were investigated by numerical simulations in three configurations. The simplest one was a Venturi-type section in which the cavitation sheet oscillates periodically with vapour cloud shedding. The second one was a hydrofoil whose unsteady cavitating behaviour depends on the angle of attack, and the most complex one was a cascade of three hydrofoils. In this last configuration, in addition to the unsteadiness associated with each cavity, a coupling between the three channels was also observed. These cavitating flows were simulated by 2D computations. Resolution of Reynolds-averaged Navier-Stokes equations was based on a finite-volume discretization associated with a pressure correction algorithm. Cavitation was simulated by using a barotropic vapour/liquid state law that links the fluid density evolution to the pressure variations. As standard $k-\epsilon$ RNG or $k-\omega$ turbulence models were found to be weakly efficient to simulate unsteady cavitation, influence of the compressibility of the two-phase medium on turbulence was considered. Both $k-\epsilon$ RNG and $k-\omega$ turbulence models, including corrections of these compressibility effects, were applied and results obtained were consistent with experiments: in the three configurations, the oscillation frequencies, the cavity length, the void ratio and the velocity fields obtained by numerical simulation were in reliable agreement with the available experimental data.
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Physical and numerical models</td><td>4</td></tr><tr><td>2.1</td><td>Physical approach</td><td>4</td></tr><tr><td>2.2</td><td>Numerical resolution</td><td>5</td></tr><tr><td>2.3</td><td>Turbulence model</td><td>6</td></tr><tr><td>3</td><td>Cavitation behaviour of the Venturi-type section</td><td>7</td></tr><tr><td>3.1</td><td>Experimental behaviour</td><td>7</td></tr><tr><td>3.2</td><td>Calculation with no effect of compressibility on the turbulence</td><td>8</td></tr><tr><td>3.3</td><td>Calculation including the compressibility effects</td><td>9</td></tr><tr><td>3.4</td><td>Discussion</td><td>11</td></tr><tr><td>4</td><td>Cavitation behaviour of a hydrofoil</td><td>12</td></tr><tr><td>4.1</td><td>Experimental behaviour</td><td>12</td></tr><tr><td>4.2</td><td>Calculation</td><td>13</td></tr><tr><td>5</td><td>Case of a cascade of three hydrofoils</td><td>14</td></tr><tr><td>5.1</td><td>Experimental behaviour</td><td>14</td></tr><tr><td>5.2</td><td>Calculation</td><td>14</td></tr><tr><td>6</td><td>Conclusion</td><td>18</td></tr></table>

## 1. Introduction

In the case of cavitation, the two-phase mixture is unsteady at the microscopic scale, since vapour and liquid structures are in continuous interaction through mass and momentum exchanges. At a larger scale, typically the size of the global cavitation area, unsteady behaviours are also observed in experiments. For example, the cavitation sheets observed on the suction side of turbomachinery blades often adopt a pulsating behaviour [1]. Rotating cavitation configurations are also reported, mainly in space turbopump inducers (see references [2, 3, 4, 5]). These consist of unbalanced attached cavities in the different channels, their distribution rotating faster than the inducer.

In such axial pumps, the tip vortex cavitation also appears highly unsteady, with no discernible arrangement. All these effects, directly due to cavitation, are responsible for spurious transient phenomena, such as flow rate fluctuations, noise, vibration, erosion, and radial or axial loads on the shaft or the casing [6].

The experimental study and the numerical modelling of these unsteady effects are thus essential to understand their mechanisms and to reduce their consequences in turbomachinery.

The pulsating behaviour of cavities attached to the blades is usually studied in simplified 2D non-rotating configurations, to focus on the physical effects associated with cavitation. To simulate a 3D pump geometry, we have considered three 2D simplified geometries as illustrated in figure 1. They are a cascade of hydrofoils, a single hydrofoil, and a Venturi-type section (whose lower wall represents the suction side of a inducer blade).

In these three configurations, at low-pressure conditions a cavitation sheet appears, respectively at the leading edge of the hydrofoils, and at the throat of the Venturi-type section. When hydrofoil angle of attack and Venturi divergence angle are high enough, this cavity becomes unstable. Three behaviours, indicated in figure 2, can be observed experimentally: they are
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Successive 2D simplifications of an inducer geometry.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Different types of unsteady cavitation observed on a foil.</figcaption></figure>

usually called respectively sheet cavitation (only the rear part of the cavity is unsteady, with a rather chaotic pulsating behaviour of low magnitude), cloud cavitation (the main part of the cavity is oscillating periodically with cloud shedding), and shear cavitation (in the case of high angle of attack, no attached cavitation sheet remains permanently at the leading edge, and only clouds of vapour appear, which are immediately carried out by the stream). An example of cloud cavitation observed experimentally on a hydrofoil suction side is given in figure 15.

The numerical simulation of cavitating flow is usually based on a single-fluid approach: the liquid/vapour mixture is considered as a homogeneous medium whose density varies from the
pure vapour one to the pure liquid one. So the local void fraction $\alpha$ is directly linked to this density evolution, and physical cavitation models mainly differ in the way these two related variables are calculated. They can be managed either by a state law (see references [7], [8], [9]), using a supplementary equation relating the void fraction to the dynamic evolution of bubble cluster (see references [10] and [11]), or by a multiple species approach with a mass transfer law between liquid and vapour [12]. The last model can be used in a full two-fluid approach, including relative motion between the phases [13].

This paper presents two-dimensional simulations of unsteady cavitating flows. The cavitation treatment is based on the physical model developed by Delannoy and Kueny [7] for inviscid fluids. Reynolds-averaged Navier-Stokes equations are resolved for the liquid/vapour mixture considered as a single fluid. The two-phase structure is defined by the local void ratio, which is directly governed by the fluid density, varying between the pure liquid one and the pure vapour one. Density evolution is related to local pressure variations through the fluid barotropic state law. This physical model for cavitation was implemented numerically into a pressure-correction scheme derived from the SIMPLE algorithm. Several physical modifications have been investigated successively by Reboud and Delannoy [14], Reboud et al [15] and Coutier-Delgosha et al [16] to increase the range of applications and improve the physical modelling.

In this physical approach, the high compressibility of cavitation areas, as well as the very low compressibility of pure liquid and pure vapour, are both taken into account. These two very different fluid properties are responsible for the main numerical difficulties that are encountered. Another crucial point is the modelization of turbulence: the effect of cavitation on the turbulence structure is not well known. Although it may have an important impact on the complex mechanisms that govern unsteady cavitation [15], no fully efficient model has been reported yet. However, a first attempt to take into account the influence of only one of a cavitating fluid properties, namely its high compressibility, has been proposed in a previous paper [17].

In the present cavitation model, the vapourization and condensation processes correspond mainly to highly compressible flow areas. Therefore, it is of primary importance to take into account the compressibility of the vapour/liquid mixture in the turbulence model. Several two-equation turbulence models were tested in the case of cavitation in a Venturi-type section, and important differences were noticed, whether empirical corrections of the compressibility effects were considered or not. It was finally found that correction terms proposed by Wilcox [18] in the case of compressible flows, as well as arbitrary modifications investigated by Reboud et al [15] were both adequate to obtain the cyclic cavitation behaviour reported in experiments.

This paper presents an interpretation of this conclusion, and results obtained by applying the compressibility corrections in three different configurations, from quite simple (the Venturi-type section) to more complicated (cavitation on a hydrofoil, and then in a cascade of three hydrofoils). In all cases, the unsteady behaviour reported experimentally is correctly predicted by the numerical model. Comparisons with experimental data focus on the length of the mean cavitation sheet for different flow conditions and oscillation frequencies. Several video animations are proposed to illustrate and to analyse the unsteady behaviour of cavitation.

# 2. Physical and numerical models

## 2.1. Physical approach

This work applies a single-fluid model based on previous numerical and physical studies (see references [7, 14, 15, 16, 17]). The fluid density $\rho$ varies in the computational domain according to a barotropic state law $\rho(P)$ that links the density to the local static pressure (figure 3). When the pressure in a cell is higher than the neighbourhood of the vapour pressure $P_v(P > P_v+(\Delta P_v/2))$, the fluid is supposed to be purely liquid. The entire cell is occupied by liquid, and its density
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Barotropic state law $\rho(P)$. Water $20^{\circ}$C.</figcaption></figure>

$\rho_l$ is calculated by the Tait equation [19]. If the pressure is lower than the neighbourhood of the vapour pressure ($P < P_v - (\Delta P_v/2)$), the cell is full of vapour and its density $\rho_v$ is given by the perfect gas law (isotherm approach). Between purely vapour and liquid states, the cell is occupied by a liquid/vapour mixture, which is considered as one single fluid with a variable density $\rho$. This is directly related to the void fraction $\alpha = (\rho - \rho_l)/(\rho_v - \rho_l)$ corresponding to the local ratio of vapour contained in the mixture.

To model the mixture state, the barotropic law presents a smooth link in the vapour pressure neighbourhood, in the interval $\pm(\Delta P_v/2)$. In direct relation with the range $\Delta P_v$, the law is characterized mainly by its maximum slope $1/A_{\min}^2$, where $A_{\min}^2 = \partial P/\partial \rho$. $A_{\min}$ can thus be interpreted as the minimum speed of sound in the mixture. Its calibration has been done previously [16]. The optimal value was found to be independent of the hydrodynamic conditions, and is about 2 m s$^{-1}$ for cold water, with $P_v = 0.023$ bar, and corresponding to $\Delta P_v \approx 0.06$ bar (chart in dotted line on figure 3). That value is applied for the computations presented hereafter.

Mass fluxes resulting from vapourization and condensation processes are treated implicitly by the barotropic state law, and no supplementary assumptions are required. Concerning the momentum fluxes, the model assumes that locally velocities are the same for liquid and for vapour: in the mixture regions vapour structures are supposed to be perfectly carried by the main flow. This hypothesis is often assessed to simulate sheet-cavity flows, in which the interface is considered to be in dynamic equilibrium [9]. The momentum transfers between the phases are thus strongly linked to the mass transfers.

## 2.2. Numerical resolution

To solve the time-dependent Reynolds-averaged Navier-Stokes equations associated with the barotropic state law presented above, the numerical code applies, on 2D structured curvilinear-orthogonal meshes, the SIMPLE algorithm, modified to take into account the cavitation process. This uses an implicit method for the time discretization, and the HLPA non-oscillatory second-order convection scheme proposed by Zhu [20]. The numerical model is detailed in Coutier-Delgosha et al [16]: a complete validation of the method was performed, and the influence of the numerical parameters was widely investigated.

In the code, the velocity field is imposed at the computational domain inlet, and the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Curvilinear-orthogonal mesh of the Venturi-type section (160 × 50 cells). (b) Zoom of the mesh in the throat region of the Venturi duct.</figcaption></figure>

static pressure is imposed at the outlet. Along the solid boundaries, the turbulence models are associated with laws of the wall (figure 4(a)).

To start unsteady calculations, the following numerical procedure is applied: first, a stationary step is carried out, with an outlet pressure high enough to avoid any vapour in the whole computational domain. Then, this pressure is lowered slowly at each new time step, down to the value corresponding to the desired cavitation number $\sigma$ defined as $(P_{\text{upstream}} - P_v)/(\rho V_{\text{upstream}}^2/2)$. Vapour appears during the pressure decrease. The cavitation number is then kept constant throughout the computation.

## 2.3. Turbulence model

In our previous studies [17], either the k-ω model proposed by Wilcox [18] or the k-ε RNG model presented by Yakhot *et al* [21] were applied to model cavitating flows. Results obtained have demonstrated that for both models, corrections of the influence of the vapour/liquid mixture compressibility on the turbulence should be taken into account to obtain the unsteady effects due to cavitation. The corrections proposed by Wilcox [18] in the k-ω model were implemented in the numerical code and the first results obtained were very promising. This section presents the main features of this turbulence model.

Equations are similar to the governing equations solved in the case of the k-ε model: all the modifications result from the use of the dissipation specific rate $\omega$ instead of the turbulence dissipation $\epsilon$. The turbulent viscosity expression and the k equation are modified, and the $\epsilon$
equation is replaced by a $\omega$ equation:

$$ 
\begin{aligned}
\nu_t &= \frac{k}{\omega} \\
\frac{\partial k}{\partial t} + u_j \frac{\partial k}{\partial x_j} &= \tau_{ij} \frac{\partial u_i}{\partial x_j} - \beta_w^* k \omega + \frac{\partial}{\partial x_j} \left[ (\nu + \sigma_w^* \nu_t) \frac{\partial k}{\partial x_j} \right] \\
\frac{\partial \omega}{\partial t} + u_j \frac{\partial \omega}{\partial x_j} &= \alpha_w \frac{\omega}{k} \tau_{ij} \frac{\partial u_i}{\partial x_j} - \beta_w \omega^2 + \frac{\partial}{\partial x_j} \left[ (\nu + \sigma_w \nu_t) \frac{\partial \omega}{\partial x_j} \right]
\end{aligned}
$$

where $\alpha_w = 5/9$, $\beta_w = 3/40$, $\beta_w^* = 9/100$, $\sigma_w = 0.5$, and $\sigma_w^* = 0.5$ are the coefficients of the model;

$$ \tau_{ij} = 2\nu_t S_{ij} \quad \text{with} \quad S_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}). $$

Influence of the compressibility is derived from an analysis in supersonic conditions: Wilcox [18] observed that the reduction of the mixing layer growth rate as a function of the Mach number could not be well simulated only by taking into account the mean density evolution through the mixing length, and that an explicit supplementary correction must be applied. He concluded that the Reynolds-stress transport equation is directly concerned with the compressibility effects, and introduced the influence of the Mach number in the turbulence equations.

The modifications proposed by Wilcox [18] are based on the previous studies of Sarkar et al [22] and Zeman [23], who aimed to take into account the compressibility effects in a $k-\epsilon$ turbulence model. He adapted this development to the $k-\omega$ model to propose the following equations giving the parameters $\beta_w$ and $\beta_w^*$ as functions of the turbulence Mach number defined as $M_t^2 = 2k/A^2$ (where A is the local speed of sound): $\beta_w^* = \beta_{wi}^*(1 + \xi^*F(M_t))$ and $\beta_w = \beta_{wi} - \beta_{wi}^*\xi^*F(M_t)$, where the index i indicates the values of the incompressible model given before,

$$ 
\begin{aligned}
\text{with } F(M_t) &= \begin{cases} (M_t^2 - M_{t0}^2) & \text{if } M_t > M_{t0} \\ 0 & \text{if } M_t \le M_{t0} \end{cases} \\
M_{t0} &= 0.25 \quad \text{and} \quad \xi^* = 1.5.
\end{aligned}
$$

We have applied this turbulence model to perform the calculations described above concerning the geometries of a Venturi, a hydrofoil and a hydrofoil cascade.

# 3. Cavitation behaviour of the Venturi-type section

## 3.1. Experimental behaviour

The first geometry considered is the Venturi illustrated in figure 4(a). Convergent and divergent angles of the lower wall of the Venturi-type section are respectively about 18° and 8° (figure 4(a)). The shape of the Venturi bottom downstream from the throat simulates an inducer blade suction side with a bevelled leading edge geometry [24] and a chord length $L_{ref} = 224$ mm.

According to experimental observations [25], in this geometry the flow is characterized by unsteady cavitation behaviour, with quasi-periodic fluctuations. Each cycle is composed of the following successive steps: the attached sheet cavity grows from the Venturi throat. A reentrant jet is generated at the cavity closure and flows along the Venturi bottom toward the cavity upstream end. Its interaction with the cavity surface results in the cavity break off. The generated vapour cloud is then convected by the main stream, until it collapses.

For a cavitation number $\sigma$ of about 2.4 (based on the time-averaged upstream pressure) and an inlet velocity $V_{ref} = 7.2$ m s$^{-1}$, vapour shedding frequency $f$ observed experimentally is about
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Time evolution of the cavity length. Abscissa, time; ordinate, the X position in the tunnel of cavitation. The colours represent density values: white for pure liquid and from red to dark blue for vapour. At a given point in time and position, the colour indicates the minimum density in the corresponding cross section of the cavitation tunnel. Calculation conditions: $\sigma \approx 2.4$; $V_{\text{ref}} = 7.2 \ \text{m s}^{-1}$; mesh = 160 × 50; time step $\Delta t = 0.005 T_{\text{ref}} (T_{\text{ref}} = L_{\text{ref}}/V_{\text{ref}})$ [17].</figcaption></figure>

50 Hz for a cavity length of 45 ± 5 mm [25]. A more detailed study concerning these calculations is presented in Coutier-Delgosha *et al* [17]. This section summarizes the main results obtained.

## 3.2. Calculation with no effect of compressibility on the turbulence

The computational grid is composed of 160 × 50 orthogonal cells (figure 4(a)). A special contraction of the mesh is applied in the main flow direction just after the throat, so that the two-phase flow area is efficiently simulated: about 50 grid points are used in this direction to model the 45 mm long mean cavity obtained hereafter (figure 4(b)). In the other direction, a contraction is also applied close to the walls, to obtain at the first grid point the non-dimensional parameter $y^+$ of the boundary layer varying between 30 and 100 and to use standard laws of the walls. The grid is finer in the bottom part of the Venturi section than in its upper part, to enhance the accuracy in the cavitation domain: cavities obtained in the following sections contain about 30 cells across their thickness.

The k-ω turbulence model presented in section 2.3 is applied with wall functions. No effect of two-phase flow compressibility on the turbulence structure is considered in this section. Thus, the fluid compressibility is only taken into account in the turbulence equations through the mean density $\rho$ change.

With this model, the unstable cavitating behaviour observed experimentally is not correctly simulated: after an initial transient fluctuation of the cavity length, the numerical calculation leads to a quasi-steady behaviour of the cavitation sheet, which globally stabilizes (figure 5).

The resulting cavity length is much too small, compared with the experimental observations (in this case, the error is larger than 50%). Moreover, comparisons with experimental data obtained by double optical probes by Stutz and Reboud [25, 26] show that the numerical mean void ratio is over-estimated in the main part of the cavity. Calculations give a high time-averaged void ratio in the upstream part of the cavitation sheet (>90%), abruptly falling to 0% in the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Transient evolution of unsteady cavitating flow in the Venturi-type duct. Temporal evolution (abscissa) of the cavity length (ordinate). Instantaneous density distribution of attached and cloud cavities is illustrated on the left at $T = 11T_{\text{ref}}$ (velocity vectors are drawn only 1 cell over 2 in the two directions) [17].</figcaption></figure>

wake, while the measured void ratio never exceeds 25% and decreases slowly from the cavity upstream end to its wake.

This poor agreement with the real configuration seems to be related to an over-prediction of the turbulent viscosity in the rear part of the cavity. The cyclic behaviour of the cloud cavitation process is strongly related to the re-entrant jet development from the cavity closure [19]. As a matter of fact, the main problem in the turbulent flow simulations consisted in the premature removal of the reverse flow along the solid wall: the re-entrant jet was stopped too early and it did not result in any cavity break off.

It is worth noting that numerical tests reported by Coutier-Delgosha et al [16] confirm that the mesh size, spatial scheme, and the time discretization applied in the model do not modify the results. Computations performed with a finer mesh (264 × 90) and first- or second-order accurate time discretization schemes still lead to the same complete stabilization of the cavitation sheet.

### 3.3. Calculation including the compressibility effects

As mentioned above, the proposed numerical simulations take into account a single-flow physical model to describe the cavitation phenomenon. According to the adopted barotropic state law, in the vapour/liquid mixture zones the speed of sound $A = \sqrt{\partial P/\partial \rho}$ is very low, and the fluid is locally highly compressible. Thus the flow often reaches Mach numbers larger than 5 in the vapourization or condensation areas. In fact, Birch and Eggers [27] showed that an increase of the Mach number could substantially modify the turbulence structure in some configurations, such as mixing layers. So the corrections proposed by Wilcox [18] and presented in section 2.3 are applied in the case of the Venturi-type section.

This leads to substantial changes in the simulation. The prediction of the unsteady re-
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Frequencies of the cloud shedding obtained by experimental and numerical studies.</figcaption></figure>

entrant jet is now obtained, and the vapour cloud shedding is well simulated. The transient evolution observed during this unsteady calculation is shown in figure 16 for the following flow conditions: $V_{\text{ref}} = 7.2 \text{ m s}^{-1}$, $\sigma = 2.4$.

For more quantitative comparisons with experiments, figure 6 illustrates at a given time and for each cross section of the Venturi-type duct the value of the minimal density in the section. By comparison with figure 5, it gives information concerning the vapour cloud shedding process: the part of the cavity that breaks off clearly appears, and the fluctuation frequency can be evaluated. Moreover, it also supplies the maximum void ratio in each section. We obtain in this case a fluctuating cavity whose mean length is about $L_{\text{cav}} = 45 \text{ mm}$ (i.e. $L_{\text{cav}}/L_{\text{ref}} = 0.2$, with $L_{\text{ref}}$ the geometry reference length). The mean fluctuation period is about $0.59T_{\text{ref}}$, and the corresponding shedding frequency equals 55 Hz. Both length and frequency are in good agreement with measurements (i.e. 45 mm and 50 Hz respectively). The transient evolution is almost periodic, with some random disturbances affecting the oscillation regularity. This behaviour is consistent with the experimental observations, which pointed to a quite regular cavitation cycle, whose frequency fluctuated around a central value [25].

The reliable agreement between experiments and calculation is confirmed by figure 7: the frequency of the self-oscillation behaviour is drawn with respect to the ratio $V_{\text{ref}}/L_{\text{cav}}$. The maximum length of the attached cavity given by the numerical simulation is chosen as an estimate of the visual length reported from the experiments. It can be seen that the cavitation oscillation frequency obtained by the calculation is very close to the measured one for a large range of cavity lengths. Dividing the shedding frequency by the ratio $V_{\text{ref}}/L_{\text{cav}}$ gives an almost constant Strouhal number $Str = f L_{\text{cav}}/V_{\text{ref}} = 0.27 \pm 0.02$, which is a classical value for vapour cloud shedding (see references [1, 25, 28]).

Figure 17 presents two phase-averaged cavitation cycles. The upper one results from experimental visualizations: video frames acquired during a 100 ns exposure time under laser sheet light are identified and digitized into 256 grey levels. A sampling technique is applied to classify them into nine sets corresponding to the different states of the recorded quasi-periodic pressure signal. Then, averaging the grey levels pixel per pixel for each set, allows one to draw a sequence of phase-averaged images, from an initial data set of 300 frames. Maximum standard deviation computed for each set is about 2/3 of the maximum mean grey level, and is mainly observed along the attached cavity and cloud boundaries.

The bottom part of figure 17 corresponds to the same sequence, obtained by numerical
<figure><img src="image_9.png" /><figcaption>Figure 8. Modification of the mixture turbulent viscosity (n = 10).</figcaption></figure>

simulation. This one lasts 20T<sub>ref</sub>, i.e. about 30 cycles. The same sampling technique is used: the computational result is decomposed into 30 × 9 short sequences corresponding to the nine steps of the cavitation cycle and the phase-averaging process is applied. We observe a good agreement between the two results concerning the external shape and global structure of both attached cavity and vapour cloud shedding.

3.4. Discussion

Results obtained with the k-ω turbulence model point to the relevant role of the compressibility effects on the turbulence: unsteadiness due to cavitation is obtained in the computation only if the corrections proposed by Wilcox to consider these effects are applied.

In fact, these corrections are derived from considering possible density fluctuations: a supplementary term appears from the averaged equations, increasing the turbulent dissipation. The final effect is a reduction of the turbulent viscosity in the compressible areas, i.e. in our particular case in the vapour/liquid mixture zones. The same result is obtained by reducing directly the turbulent viscosity in the two-phase areas, as proposed by Reboud et al [15] with a k-ε RNG turbulence model. This model is originally devoted to fully incompressible fluids, and no particular correction is applied in the case of the highly compressible two-phase mixture. To improve the turbulence modelling and to try to better simulate the re-entrant jet behaviour and the vapour cloud shedding, we modified arbitrarily the standard model simply by reducing the mixture turbulent viscosity, mainly in the low void ratio areas, as represented in figure 8:

$$
\mu_t = f(\rho)C_{\mu}k^2/\epsilon \qquad C_{\mu} = 0.0085
$$

where

$$
f(\rho) = \rho_v + \left( \frac{\rho_v - \rho}{\rho_v - \rho_l} \right)^n (\rho_l - \rho_v) \quad n \gg 1.
$$

The results obtained were almost identical to those presented previously for the k-ω model
including compressibility effects [17].

The reduction of turbulent viscosity in the two-phase area, which appears to be the direct
cause of the simulation improvement, is not confirmed yet by any measurements. However, this
trend is consistent with the one reported by Wang *et al* [29], who mentions for bubbly flows
an increase of turbulent dissipation with void ratio, and thus a possible diminution of velocity
fluctuations.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Geometry of the hydrofoil and computational grid ($L_{\text{ref}}$ = 15 cm).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Mean shape of the cavitation sheet observed on the hydrofoil tested in the ENSTA test rig (Angle of attack 5°, $V_{\text{ref}}$ = 6.8 m s⁻¹, $\sigma$ = 1.65, length scale indicated on the foil each centimetre.).</figcaption></figure>

# 4. Cavitation behaviour of a hydrofoil

## 4.1. Experimental behaviour

The second geometry is a two-dimensional foil of 150 mm chord and 80 mm span. Its cross section is composed of a flat upper surface and a convex lower surface of 195 mm radius, as illustrated by figure 9. Experiments were performed in the ENSTA cavitation tunnel whose test section is 150 mm height, 80 mm width, and 640 mm length. The foil was located at mid-height, with 5° angle of attack, so a cavitation sheet appears on the upper face, when the pressure is decreased in the tunnel. Compared with the previous configuration, the cavitation behaviour is modified by the possible interaction between foil pressure and suction sides.

For an inlet velocity equal to 6.8 m s⁻¹ and a cavitation parameter $\sigma$ = 1.65, a cavitation sheet is obtained on the upper surface. Its mean length is about 60 ± 5 mm, but its behaviour is unsteady with large periodical vapour cloud shedding. Figure 10 shows the mean cavitation sheet.

It is worth noting that the length of this mean cavity (about 90 mm) includes the vapour
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison between numerical and experimental oscillation frequencies concerning the hydrofoil geometry.</figcaption></figure>

cloud before it begins to collapse. So it is notably higher than the mean length of the attached cavity reported during the experiment.

The cavity oscillation frequency, measured from the cavitation noise with a hydrophone, equals 33 ± 3 Hz. It was reported in previous studies [30] that in this configuration the vapour shedding is directly driven by the re-entrant jet that periodically flows under the cavity towards the leading edge.

## 4.2. Calculation

The numerical model is applied with the *k*-ω model including compressibility effects. The computational domain is composed of 140 × 70 cells, with a special contraction in the expected cavitating area, as indicated in the case of the Venturi-type section (figure 9). The experimental flow conditions are imposed ($V_{\text{ref}} = 6.8 \text{ m s}^{-1}$, $\sigma = 1.65$).

The predicted flow behaviour is then in good general agreement with experiments: the cavitation sheet adopts a periodical unsteady behaviour, with a mean frequency of 30/35 Hz and a mean cavity length equal to 60 mm. These values are in reliable agreement with experimental data.

Figure 18 presents a comparison between phase-averaged cavity shapes, obtained respectively in the upper part from the calculation, and in the bottom part from pictures of the cavitation sheet. It can be noticed that the evolution of the cavity shape in both cases is very similar. In the calculation, a small-pressure side cavitation sheet is also obtained, whereas it is observed in experiments only for slightly lower pressure or higher angle of attack.

Several flow conditions were investigated by varying the cavitation number. The oscillation frequencies obtained for different mean cavity lengths are reported in figure 11, as well as the corresponding Strouhal numbers. These are all close to the standard value of 0.3 reported by many authors in such configurations (see [1, 25]). The experimental measurement for $V_{\text{ref}}/L_{\text{cav}} = 194 \text{ Hz}$ must be considered circumspectly, since the visual estimate of the mean cavity length is not precise for such small sheets of cavitation. The high value of the Strouhal number in this case ($Str = 0.36$), as well as the discrepancy with the numerical result, could be related to this uncertainty. For other flow conditions, calculation and experimental results are
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Geometry of the cascade of hydrofoils tested in the Laboratory for Turbomachinery and Fluid Power, Darmstadt.</figcaption></figure>

in good general agreement, which indicates that the physical cavitation model, associated with the $k-\omega$ turbulence model (including compressibility effects) is efficient to predict the behaviour of the cavitation sheet in this more complex configuration.

# 5. Case of a cascade of three hydrofoils

## 5.1. Experimental behaviour

The last test case—the most complex one—is a cascade composed of three hydrofoils, located in a test section whose size is 30 × 100 mm (figure 12). The experiments were performed at the Laboratory for Turbomachinery and Fluid Power of the Darmstadt University of Technology. The angle of attack of all foils is 5° in the results presented hereafter, and the inlet velocity equals 14 m s⁻¹. A piezoresistive pressure transducer (Kulite XTM-190M) with a frequency range up to 125 kHz is used to investigate the unsteadiness of cavitation. The sensor membrane is connected to the flow by a small opening. All details concerning the experimental setup and the pressure signal treatment are given by Lohrberg *et al* [31].

According to the experimental study, the cavitation phenomenon, for different values of $\sigma$, is characterized by two frequencies: the higher one is found to be related to the self-oscillation behaviour of the cavitation sheet on the upper foil, while the lower one corresponds to a coupling between the different channels.

## 5.2. Calculation

This section presents results obtained in this geometry by calculations for $\sigma = 1.95$. More details concerning numerical simulations in this configuration are presented in reference [31].

The turbulence model applied in this configuration is the modified $k-\epsilon$ RNG model, whose modifications were found in section 3 to be almost equivalent to the treatment of the compressibility effects in the $k-\omega$ model. Vapour structures, characterized in the model by small values of the density $\rho$, are generated along the suction side of the foils from the leading edge. After the initial transient, a special unsteady behaviour is observed (figure 13):

* On the upper foil suction side, typical cloud cavitation behaviour takes place. However, the size and shedding period of the cavitation structures show strong variations between the successive cycles.
<figure><img src="image_14.png" /></figure>

Density field at t = 160 and 172 ms (=45 and 48 tref)

Variation with time (Tref = Lref/Vref = 3.6 ms) of the size of the cavitation structures in the channels.

**Figure 13.** Time evolution of the cavitation structures on the three hydrofoils ($V_{ref}$ = 14 m s⁻¹, σ = 1.95) [31].

* On the two other foils the cavities appear as almost stable sheets attached to the blades, but with large amplitude and low frequency length fluctuations.

The difference between the cavitation behaviour in channel 1 and the others is associated with the larger flow angle of attack observed in front of the upper hydrofoil. This effect is due to the presence of the upper wall. The divergence angle between the wall and the foil suction side, associated with the trailing edge condition, increases the flow rate passing over the foil. An opposite effect occurs between the third foil and the lower wall. With the larger angle of attack, typical unsteady cloud cavitation takes place, while stable sheet cavities with low relative thickness develop on the two other foils (channels 2 and 3). It can be seen in figure 13 that the fluctuations of the cavitation structures (i.e. cloud cavitation in channel 1 and attached cavitation sheets in channels 2 and 3) appear with completely opposite phases. Disappearance of the cavitation structures in the upper channel corresponds to a maximum development of the cavitation sheets on the two other foils, and vice versa.
<figure><img src="image_15.png" /></figure>

**Figure 14.** (a) Angle of velocity vector, 5% upstream of the upper foil leading edge. (b) Longitudinal velocity near the upper foil trailing edge (suction side) [31].

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

**Figure 15.** Example of an unsteady cloud cavitation on a foil suction side (see animation).

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

**Figure 16.** Numerical simulation of unsteady oscillation behaviour in the Venturi type section (see animation).

To evaluate the characteristic frequencies of the cavitating flow, the pressure signals obtained by numerical calculations have been analysed: the lowest characteristic frequency obtained equals 31.42 ± 0.98 Hz, while the frequency of cloud shedding process in channel 1 is estimated as 262 ± 35 Hz. Based on the corresponding mean attached cavity length (~0.3 L ref ~1.5 cm), the inlet velocity (V ref = 14 m s⁻¹) and a frequency of 262 Hz, the Strouhal number is found to be
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Comparison between numerical and experimental phase-averaged cavity shapes for the Venturi type section (see [animation](#)).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Comparison between numerical and experimental phase-averaged cavity shapes for the hydrofoil (see [animation](#)).</figcaption></figure>

about 0.28. These values present an order of magnitude in good agreement with the experimental measurements at $\sigma = 1.95$ (respectively about 40–50, and 110–180 Hz). The whole unsteady behaviour including the two different frequencies can be observed on figure 19.

To investigate the reason for the coupling at low frequency between the cavitating flows in the different channels, the time evolution of the local flow angle in front of the upper foil was first drawn (figure 14(a)). As explained before, the occurrence of large cloud cavitation in channel 1 corresponds to the increase of the incident flow angle. In contrast, periods with stable cavitation sheets on the second and third foils correspond to a small flow angle. Longitudinal velocity near the trailing edge, in the boundary layer of the suction side of the upper foil, is shown in figure 14(b) with respect to time. This shows that periods of developed cloud cavitation are followed by the passage at the foil trailing edge of travelling vortices, with negative longitudinal velocity. The interactions between these travelling flow structures and the foil trailing edge seem at to be at the origin of the decrease of the flow angle in front of the foil.

Another illustration of this phenomenon is given in figure 20, where the flow vorticity evolution is observed during several low-frequency cycles. The period of the low-frequency
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Numerical simulation of the unsteady cavitating behaviour of the cascade of hydrofoils (see [animation](#)).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Evolution of the vorticity field in the hydrofoils cascade (see [animation](#)).</figcaption></figure>

behaviour seems to be linked to the convection time of the large vortices generated by unsteady
cavitation, to the foil trailing edge.

**6. Conclusion**

Unsteady cloud cavitation has been investigated in three configurations, of increasing complexity.
Experimental measurements and visualization have been compared to the results of two-
dimensional computations, performed with a single-fluid cavitation model based on a barotropic
state law.

The first case was a simple Venturi-type section, for which several turbulence models were
applied. It was found that it is necessary to take into account corrections of the effects of
the vapour/liquid mixture compressibility on the turbulence structure to obtain the unsteady
cavitating behaviour observed experimentally. This influence of compressibility was tested as
well with a k-ω model proposed by Wilcox [18] for compressible flows as with a k-ε RNG model
modified arbitrarily. Results obtained are in reliable agreement with experimental data. In both
cases, the main reason for the improvement of the prediction is a reduction of the turbulent
viscosity in the two-phase flow area.

The ability of the physical cavitation model, associated with any of these two turbulence
models, to predict unsteady behaviour of cavitation sheets was confirmed by two complementary
test cases, namely a single hydrofoil and a cascade of three hydrofoils. In these two configurations,
the numerical model correctly predicts both cavity lengths and oscillation frequencies, for given
flow conditions.

Concerning the hydrofoil cascade, the numerical approach allowed us to obtain a detailed
analysis of the cavitating flow, and a reliable prediction of the characteristic frequencies of the
phenomena of cloud shedding and of the coupling between the three channels. This indicates that
not only single oscillation behaviours, but also other types of unsteady effects due to cavitation,
may be efficiently simulated by the proposed model.

## Acknowledgments

The authors wish to express their gratitude to the French Space Agency, Centre National d'Etudes Spatiales (CNES), for its continuous support. The study of the hydrofoil cascade was performed in collaboration with Dr H Lohrberg and Professor B Stoffel, from the Laboratory for Turbomachinery and Fluid Power of Darmstadt. The experimental results concerning the Venturi-type section were investigated by Dr B Stutz, at the LEGI laboratory.

## References

[1] Joussellin F, Delannoy Y, Sauvage-Boutar E and Goirand B 1991 experimental investigations on unsteady attached cavities *ASME 1991 Cavitation Symp. (Portland, USA, 1991)*

[2] Kamijo K, Shimura T and Watanabe M 1977 An experimental investigation of cavitating inducer instability
*ASME Paper 77-WA/FW-14*

[3] de Bernardi J, Joussellin F and Von Kaenel A 1993 Experimental analysis of instabilities related to cavitation
in turbopump inducer *Proc. 1st Int. Symp. on Pump Noise and Vibrations (Paris, France, 1993)* pp 91–9

[4] Tsujimoto Y, Kamijo K and Yoshida Y 1993 A theoretical analysis of rotating cavitation inducers *J. Fluids Eng.* **115** 135–41

[5] Pagnier P, Morel P, Spettel F, Henry C and Champagne J Y 1995 Conception and experimental study of an
inducer *Proc. Cav’95 Int. Symp. (Deauville, France, 1995)*

[6] Goirand B, Mertz A, Joussellin F and Rebattet C 1992 Experimental investigation of radial loads induced
by partial cavitation with liquid hydrogen inducer *Proc. 3rd Int. Conf. on Cavitation, ImechE C453/056*
*(Cambridge, UK, 1992)* pp 263–9

[7] Delannoy Y and Kueny J L 1990 Two-phase flow approach in unsteady cavitation modeling *Cavitation and Multiphase Flow Forum ASME-FED vol 98* pp 153–8

[8] Song C and He J 1998 Numerical simulation of cavitating flows by single-phase flow approach *Proc. 3rd Int. Symp. on Cavitation (Grenoble, France, 1998)* pp 295–300

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

[10] Kubota A, Kato H and Yamagushi H 1992 A new modelling of cavitation flows: a numerical study of unsteady
cavitation on a hydrofoil section *J. Fluid Mech.* **240** 59–96

[11] Chen Y and Heister S D 1996 Modeling hydrodynamic non-equilibrium in cavitating flows *J. Fluid Eng.* **118** 172–8

[12] Kunz R, Boger D, Chyczewski T, Stinebring D and Gibeling H 1999 Multi-phase CFD analysis of natural
and ventilated cavitation about submerged bodies *3rd ASME/JSME Joint Fluids Engineering Conf.*
*(San Francisco, CA, USA, 1999)*

[13] Alajbegovic A, Grogger H and Philipp H 1999 Calculation of transient cavitation in nozzle using the two-fluid
model *12th Annual Conf. on Liquid Atomization and Spray Systems (Indianapolis, USA, 1999)*
[14] Reboud J L and Delannoy Y 1994 Two-phase flow modeling of unsteady cavitation 2nd Int. Symp. on Cavitation (Tokyo, Japan, 1994)

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

[16] Coutier-Delgosha O, Reboud J L and Delannoy Y 2002 Numerical simulations in unsteady cavitating flows *Int. J. Numer. Methods Fluids* submitted

[17] Coutier-Delgosha O, Fortes-Patella R and Reboud J-L 2002 Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation *J. Fluid Eng.* at press

[18] Wilcox D 1998 *Turbulence Modeling for CFD* (La Canada, CA: DCW Industries Inc)

[19] Knapp R T, Daily J T and Hammit F G 1970 *Cavitation* (New York: McGraw-Hill)

[20] Zhu J 1991 A low diffusive and oscillation-free convection scheme *Commun. Appl. Numer. Methods 7*

[21] Yakhot V, Orszag S A, Thangham S, Gatski T B and Speziale C G 1992 Development of turbulence models for shear flows by a double expansion technique *Phys. Fluids A* **4** 1510–20

[22] Sarkar S, Erlebacher G, Hussaini M Y and Kreiss H O 1989 The analysis and modeling of dilatational terms in compressible turbulence *NASA CR-181959*

[23] Zeman O 1990 Dilatational dissipation: the concept and application in modelling compressible mixing layers *Phys. Fluids A* **2** 178–88

[24] Kueny J-L, Reboud J-L and Desclaux J 1991 Analysis of partial cavitation: image processing and numerical prediction *Cavitation'91* ASME-FED vol 116 pp 55–60

[25] Stutz B and Reboud J-L 1997 Experiments on unsteady cavitation *Exp. Fluids* **22** 191–8

[26] Stutz B and Reboud J-L 2000 Measurements within unsteady cavitation *Exp. Fluids* **29** 545–52

[27] Birch S F and Eggers J M 1972 Free turbulent shear flow *NASA SP-321* vol 1 pp 11–40

[28] Lush P A and Peters P I 1982 Visualisation of the cavitating flow in a Venturi-type duct using high speed cine photography *Proc. IAHR Conf. on Operating Problems of Pump Stations an Power Plants (Amsterdam)* vol 5

[29] Wang S K, Lee S J, Jones O C and Lahey R T 1987 3D turbulence structure and phase distribution measurements in bubbly two-phase flows *Int. J. Multiph. Flow* **13** 327–43

[30] Pham T M, Larrarte F and Fruman D H 1999 Investigation of unsteady sheet cavitation and cloud cavitation mechanisms *J. Fluid Eng.* **121** 289–96

[31] Lorberg H, Stoffel B, Fortes-Patella R, Coutier-Delgosha O and Reboud J-L 2002 Numerical and experimental investigations on the cavitating flow in a cascade of hydrofoils *Exp. Fluids* at press