Evaluation of the turbulence model influence on the numerical simulations of unsteady cavitation

Unsteady cavitation in a Venturi-type section was simulated by two-dimensional computations of viscous, compressible, and turbulent cavitating flows. The numerical model used an implicit finite volume scheme (based on the SIMPLE algorithm) to solve Reynolds-averaged Navier-Stokes equations, associated with a barotropic vapor/liquid state law that strongly links the density variations to the pressure evolution. To simulate turbulence effects on cavitating flows, four different models were implemented (standard k−e RNG; modified k−e RNG; k−ω with and without compressibility effects), and numerical results obtained were compared to experimental ones. The standard models k−e RNG and k−ω without compressibility effects lead to a poor description of the self-oscillation behavior of the cavitating flow. To improve numerical simulations by taking into account the influence of the compressibility of the two-phase medium on turbulence, two other models were implemented in the numerical code: a modified k−e model and the k−ω model including compressibility effects. Results obtained concerning void ratio, velocity fields, and cavitation unsteady behavior were found in good agreement with experimental ones. The role of the compressibility effects on turbulent two-phase flow modeling was analyzed, and it seemed to be of primary importance in numerical simulations.


Introduction
Cavitating flows in turbomachinery lead to performance losses and modifications of the blades load. These consequences can be related to the quasi-steady cavitation effects, mainly depending on the time-averaged shapes of the vaporized structures. Moreover, turbomachinery under cavitating conditions is also submitted to transient phenomena, such as compressibility effects, flow rate fluctuations, noise, vibration, erosion, in which the unsteady behavior and the two-phase structure of the cavitating flow are influent. The study and modeling of the unsteady behavior of cavitation are thus essential to estimate the hydraulic unsteadiness in turbomachinery and the associated effects.
In this context, some numerical models taking into account the two-phase structure and the unsteady behavior of cavitation have been developed. Most of them are based on a single fluid approach: the relative motion between liquid and vapor phases is neglected and the liquid vapor mixture is treated as a homogeneous medium with variable density. The mixture density is related directly to the local void fraction and is managed either by a state law ͑Delannoy and Kueny ͓1͔,Song and He ͓2͔,and Merkle et al. ͓3͔͒, using a supplementary equation relating the void fraction to the dynamic evolution of bubble cluster ͑Kubota et al. ͓4͔, and Chen and Heister ͓5͔͒, or by a multiple species approach with a mass transfer law between liquid and vapor ͑Kunz et al. ͓6͔͒. The last model can be used in a full two-fluid approach, including relative motion between the phases ͑Alajbegovic et al. ͓7͔͒.
Main numerical difficulties are related to the strong coupling between the pressure field and the void ratio and to the coexistence of the strong compressibility of the two-phase medium with the quasi-incompressible behavior of the pure liquid flow. Moreover, the influence of the turbulence on unsteady two-phase compressible flows is not yet well known. This paper presents unsteady cavitating flow simulations performed with a two-dimensional code. It is based on the model developed by Delannoy and Kueny ͓1͔ for inviscid fluids. Several physical modifications have been investigated by Reboud and Delannoy ͓8͔,Reboud et al. ͓9͔, to increase the range of applications and improve the physical modeling.
In the numerical code, the two-phase aspects of cavitation are treated by introducing a barotropic law that strongly links the fluid density to the pressure variations. From the numerical point of view, this physical approach is associated with a pressurecorrection scheme derived from the SIMPLE algorithm, slightly modified to take into account the cavitation process. The model has been validated on numerous cases, such as Venturi, ͓9,10͔, foils, ͓8,11͔, or blade cascades, ͓12,13͔. The numerical results showed a good agreement with experiments. Mainly, the code leads to a reliable simulation of the cyclic self-oscillation behavior of unsteady cavitating flows.
The complex unsteady mechanism that governs this cyclic cavitation behavior is strongly affected by the applied turbulence model, ͓9,10͔. The aim of this paper is to study the influence of different models ͑standard k-RNG; modified k-RNG; kincluding or not compressibility effects͒ on the numerical simulation of cavitating flows. Indeed, because of the barotropic state law adopted in the model, the vaporization and condensation processes correspond mainly to highly compressible flow areas. Therefore, it is of primary importance to take into account the compressibility of the vapor/liquid mixture in the turbulence model.
The paper describes the turbulence models applied in the numerical simulations and presents comparisons between different numerical results obtained in a two-dimensional Venturi-type section, whose experimental behavior has been already studied by Stutz and Rebouch ͓14,15͔. In this geometry, the flow is characterized by an unstable cavitation behavior, with almost periodical vapor cloud shedding.
In addition to the standard k-model, three other turbulence models were applied, to investigate their influence on the twophase turbulent flow. Mainly, the turbulence model proposed by Wilcox ͓16͔, which includes compressibility considerations, was applied, and promising results are obtained. Comparisons with experimental measurements of void ratio and velocities inside the cavity, ͓9͔, are presented.
The ability of these turbulence models to predict complex twophase flow is discussed, and the role of the compressibility effects is studied.

Physical Model
The present work applies a single fluid model based on previous numerical and physical studies, ͓1,8 -10͔. The fluid density varies in the computational domain according to a barotropic state law ( P) that links the density to the local static pressure ͑Fig. 1͒. When the pressure in a cell is higher than the neighborhood of the vapor pressure ( PϾ P v ϩ(⌬ P v /2)), the fluid is supposed to be purely liquid. The entire cell is occupied by liquid, and the density l is calculated by the Tait equation, ͓17͔. If the pressure is lower than the neighborhood of the vapor pressure ( PϽ P v Ϫ(⌬ P v /2)), the cell is full of vapor and the density v is given by the perfect gas law ͑isotherm approach͒. Between purely vapor and purely liquid states, the cell is occupied by a liquid/vapor mixture, which is considered as one single fluid, with the variable density . The density is directly related to the void fraction ␣ ϭ(Ϫ l )/( v Ϫ l ) corresponding to the local ratio of vapor contained in this mixture.
To model the mixture state, the barotropic law presents a smooth link in the vapor pressure neighborhood, in the interval Ϯ(⌬ P v /2). In direct relation with the range ⌬ P v , the law is characterized mainly by its maximum slope 1/A min 2 , where A min 2 ϭ‫ץ‬P/‫.ץ‬ A min can thus be interpreted as the minimum speed of sound in the mixture. Its calibration was done in previous studies, ͓10͔. The optimal value was found to be independent of the hydrodynamic conditions, and is about 2 m/s for cold water, with P v ϭ0.023 bar, and corresponding to ⌬ P v Ϸ0.06 bar ͑green chart on Fig. 1͒. That value is applied for the computations presented hereafter. Mass fluxes resulting from vaporization 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 vapor: In the mixture regions vapor structures are supposed to be perfectly carried along 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, ͓3͔. The momentum transfers between the phases are thus strongly linked to the mass transfers.

Numerical Model
To solve the time-dependent Reynolds-averaged Navier-Stokes equations associated with the barotropic state law presented here above, the numerical code applies, on two-dimensional structured curvilinear-orthogonal meshes, the SIMPLE algorithm, modified to take into account the cavitation process. It uses an implicit method for the time-discretization, and the HLPA nonoscillatory second-order convection scheme proposed by Zhu ͓18͔. The numerical model is detailed in Coutier-Delgosha et al. ͓10͔: A complete validation of the method was performed, and the influence of the numerical parameters was widely investigated. The present paper describes mainly the different turbulence models applied.

Boundary Conditions.
In the code, the velocity field is imposed at the computational domain inlet, and the static pressure is imposed at the outlet. Along the solid boundaries, the turbulence models are associated with laws of the wall.

Initial Transients Conditions.
To start unsteady calculations, the following numerical procedure is applied: First of all, a stationary step is carried out, with an outlet pressure high enough to avoid any vapor 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 . Vapor appears during the pressure decrease. The cavitation number is then kept constant throughout the computation.

The Geometry
Numerical simulations have been performed on a Venturi-type section whose convergent and divergent angles are, respectively, about 18 deg and 8 deg ͑Fig. 2͑a͒͒. The shape of the Venturi bottom downstream from the throat simulates an inducer blade suction side with a beveled leading edge geometry ͑Kueny et al. According to experimental observations, in this geometry the flow is characterized by unsteady cavitation behavior, ͓14͔, with quasi-periodic fluctuations. Each cycle is composed of the following successive steps: The attached sheet cavity grows from the Venturi throat. A re-entrant 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 vapor cloud is then convected by the main stream, until it collapses.
For a cavitation number of about 2.4 ͑based on the timeaveraged upstream pressure͒ and an inlet velocity V ref ϭ7.2 m/s, vapor shedding frequency observed experimentally is about 50 Hz for a cavity length of 45Ϯ5 mm, ͓14͔.
The standard computational grid is composed of 160ϫ50 orthogonal cells ͑Fig. 2͑a͒͒. A special stretching 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 ͑Fig. 2͑b͒͒. In the other direction, a stretching 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.

Turbulence Models
"a… Standard k-RNG Model. The first model applied in the numerical code is a standard k-RNG model, associated with laws of the wall, ͓9͔.
In this model, the effective viscosity applied in the Reynolds equations is defined as ϭ t ϩ l where t ϭC k 2 / is the turbulent viscosity and C ϭ0.085 ͑Yakhot et al. ͓20͔͒.
This model is originally devoted to fully incompressible fluids, and no particular correction is applied here in the case of the highly compressible two-phase mixture. Thus, the fluid compressibility is only taken into account in the turbulence equations through the mean density changes.
With this model, the unstable cavitating behavior observed experimentally is not correctly simulated: After an initial transient fluctuation of the cavity length, the numerical calculation leads to a quasi-steady behavior of the cavitation sheet, which globally stabilizes ͑Fig. 3͒.
The resulting cavity length is much too small, compared with the experimental observations ͑in this case, the error is larger than 50 percent͒. Moreover, comparisons with experimental data ob-tained by double optical probes by Stutz and Reboud ͓14,15͔ show that the numerical mean void ratio is overestimated in the main part of the cavity. Calculations give a high time-averaged void ratio in the upstream part of the cavitation sheet ͑Ͼ90 per-cent͒, abruptly falling to 0 percent in the wake, while the measured void ratio never exceeds 25 percent and decreases slowly from the cavity upstream end to its wake.
This poor agreement with the real configuration seems to be related to an overprediction of the turbulent viscosity in the rear part of the cavity. The cyclic behavior of the cloud cavitation process is strongly related to the re-entrant jet development from the cavity closure, ͓17͔. 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 in ͓10͔ 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.
"b… Modified k-RNG Model. To improve the turbulence modeling and to try to better simulate the re-entrant jet behavior  and the vapor cloud shedding, we modified the standard k-RNG model simply by reducing the mixture turbulent viscosity, mainly in the low void ratio areas: Indeed, according to the experimental results ͓15͔, the reentrant jet seems to be mainly composed of liquid (␣ϳ0), and thus the reduction of the mixture turbulence viscosity leads to substantial changes in the simulation. The prediction of the unsteady re-entrant jet is now obtained, and the vapor cloud shedding is well simulated ͑See Fig. 4.͒ The transient evolution observed during this unsteady calculation is shown in Fig. 5. Fig. 5͑a͒ 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 Fig. 3, it gives information concerning the vapor 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. Curve 5b presents the transient upstream pressure evolution.
The experimental self-oscillatory behavior of the cavitation sheet is correctly simulated ͑Fig. 5͒. We obtain in this case a fluctuating cavity whose mean length is about 45 mm ͑i.e., L cav /L ref ϭ0.2). The mean fluctuation period is about 0.59 T 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 reduced frequency based on the reference velocity and the mean cavity length equals 0.33, which only slightly overestimates the classical Strouhal number Stϭ0.3 ͑0.28 from the experimental measurements͒. The transient evolution is almost periodic, with some random disturbances affecting the oscillation regularity. This behavior is consistent with the experimental observations, which pointed out a quite regular cavitation cycle, whose frequency fluctuated around a central value, ͓14͔.
Several tests were performed to investigate the influence of the numerical parameters on the unsteady cavitating result. The whole study, including the influence of the grid size, the time-step, the order of the time discretization scheme, the ratio v / l , and the maximum slope of the barotropic law A min , is detailed in ͓10͔. Results concerning the grid resolution and the time step influence are reported in the present paper in Table 1. The effect of the mesh size appears to be small, so far it is fine enough: the cavity oscillation frequency is almost constant with the two finest grids. In contrast, the influence of the time-step cannot be completely removed: using ⌬tϭ0.005 T ref leads to a systematic slight overestimation of the cycles frequency. Nevertheless, it must be considered that the uncertainty on experimental measurements is of the same order of magnitude.
"c… k-Model: Compressibility Effects. As presented before, the proposed numerical simulations take into account a single-flow physical model to describe cavitation phenomenon. According to the adopted barotropic state law, in the vapor/liquid mixture zones the sound celerity Aϭͱ‫ץ‬ P/‫ץ‬ is very low, the fluid is locally highly compressible. Thus the flow often reaches Mach numbers larger than 5 in the vaporization or condensation areas. As a matter of fact, Birch and Eggers ͓21͔ showed that an increase of the Mach number could modify substantially the turbulence structure in some configurations, such as mixing layers. To ana-lyze and to evaluate the compressibility effects on the turbulence modeling of cavitating flows, we implemented in the numerical code the k-turbulence model proposed by Wilcox ͓16͔. Its major advantage, compared with the previous models, consists in the corrections that are proposed by Wilcox in the turbulence equations to model the effects of density fluctuations. The main features of the models are presented hereafter. All details concerning the turbulence model can be found in ͓16͔.

Basic Model Equations.
Equations are similar to the governing equations resolved in the case of the k-model: All the modifications result from the use of the dissipation specific rate instead of the turbulence dissipation . The turbulent viscosity expression and the k-equation are modified, and the equation is replaced by a equation: where ␣ w ϭ5/9, ␤ w ϭ3/40, ␤ w *ϭ9/100, w ϭ0.5, and w *ϭ0.5 are the coefficients of the model. i j ϭ2 t S i j with S i j ϭ 1 2 ͑ u i, j ϩu j,i ͒ Compressibility Effects. By analyzing the reduction of the mixing layer growth rate as a function of the Mach number, Wilcox ͓16͔ observed that the phenomenon 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 concludes 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 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-turbulence model.
He adapted this development to the kmodel to propose the following equations giving the parameters ␤ w and ␤ w * as functions of the turbulence Mach number defined as M t 2 ϭ2k/A 2 ͑where A is the local sound celerity͒.
Results. Figure 6 presents calculations performed with the k-turbulence model.
In the first case ͑Fig. 6͑a͒͒, the influence of the compressibility is not taken into account in the model. The numerical results are not satisfactory: As in the case of standard k-RNG model, the unsteady cavitation behavior is not correctly simulated and the cavity length obtained by simulations is more than 50 percent smaller than the experimental one. On the other hand, numerical results obtained by taking into account the compressibility effects are in good agreement with experimental ones. They are very similar to those presented previously in the case of the modified k-RNG model: As can be seen in Fig. 6͑b͒, the vapor cloud shedding process is correctly simulated. The mean cavity length is still correct, and the shedding frequency remains almost the same.
Comparisons are investigated with experimental data obtained by using double optical probes, ͓15͔. This intrusive technique allows measurements of the local void ratio and the velocities of the two-phase structures inside the cavitation sheet. Their timeaveraged and standard deviation values are presented for four profiles in Fig. 7.
These comparisons point out a general good agreement between experimental and numerical results. The time-averaged void ratio repartition is well predicted, and the re-entrant jet structure under the cavity appears qualitatively correct. The fluctuation levels are also in reliable close agreement with measurements. The negative velocity region at station xϭ0.065 m corresponds to the rear part of the cavity, alternatively affected by the re-entrant jet progression and by the vapor cloud shedding. Both models fail to simulate these negative velocities, because of a slight underestimation Fig. 6 Unsteady cavitation behavior calculated with the kmodel-"a… without and "b… with compressibility effects of the mean cavity length in this configuration. A second discrepancy can be observed at station 0.03 m, concerning the void ratio distribution close to the wall; nevertheless, the experimental estimation of velocities and void ratio in the re-entrant jet is associated to a high uncertainty in this area, due to important fluctuations level, ͓14͔. Thus, the comparison, along the wall, between experiments and the numerical result cannot be relevant.

Conclusion
Four turbulence models have been applied to simulate unsteady cavitating flows in a Venturi-type section.
Results obtained with the standard k-RNG model and with the k-model without taking into account compressibility effects are in poor agreement with experimental observations: the models fail to reproduce the vapor cloud shedding behavior observed experimentally. Moreover, numerical tests reported in ͓10͔ and the results presented here above indicate that the global level of dissipation induced by the numerical model seems not to be responsible for its inability to simulate the unsteady cavity behavior.
On the other hand, the modified k-RNG model and the kmodel including compressibility effects lead to a very reliable simulation of the Venturi unsteady cavitation behavior. The satisfactory results obtained with the modified k-RNG model have been confirmed by simulations performed in other geometries, such as a hydrofoil, ͓11͔, a foil cascade, ͓13͔, and another Venturitype section leading to a more stable sheet of cavitation, ͓10͔. In all cases, the general experimental behavior is correctly obtained, and oscillation frequencies are accurately predicted in unsteady configurations.
Results obtained in the present paper with kmodels pointed out the relevant role of the compressibility effects on the turbulence model. Indeed, the corrections proposed to treat compressibility effects take into account the 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 the vapor/liquid mixture zones.
The modifications proposed to improve the k-RNG model had been also based on the reduction of the turbulent viscosity in the mixture zones, which are characterized by a very low sound celerity and large Mach number. Indeed, results obtained with the modified k-RNG model are similar to the ones calculated by the compressible kmodel.
According to the numerical calculations, the fluid compressibility has a strong effect on the turbulence structure, and must be taken into account to simulate unsteady cavitating flows.