Test case n°30: Unsteady cavitation in a Venturi type section (PN)

The challenge consists in simulating the unsteady behavior of cavitation pockets formed at the throat of a converging-diverging channel. Comparisons with existing numerical predictions and experimental data are suggested.


Practical significance and interest of the test-case
A calculation of the cavitating flow field occurring in a Venturi type duct is presented in the present paper. Experimental results obtained in the same flow configuration by Stutz & Reboud (2000) are also reported. Convergent and divergent angles of the lower wall of the Venturi type section are respectively about 18 • and 8 • (see figure 2). According to experimental observations in this geometry the flow is characterized by unsteady cavitation behavior, 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 throat. 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.
The challenge consists in simulating correctly this unsteady behavior. Two tests are proposed to evaluate the consistency of the numerical solution with the experiments: • the evaluation of overall parameters (mean volume of vapor, standard deviation, frequency of the periodic fluctuations, phase average of the cavity shape evolution), • the description of the flow inside the sheet of cavitation (time-averaged values and standard deviations of velocities and void fraction) The numerical simulation of this problem requires a coupling between the Navier-Stokes equations, a model of turbulence (the Reynolds number in the proposed configuration equals 1.6 10 6 ), and a physical model of cavitation to predict the inception of cavitation and the behavior of the liquid/vapour mixture.
2 Definitions and physical model description

Physical model of cavitation
The present work considers a single fluid model: the fluid density ρ varies in the computational domain according to a barotropic equation of state, ρ(P ), that links the density to the local static pressure (see figure 1). This equation considers that phase change occurs within a small range of pressure, ∆P vap , centered on the saturation pressure P vap . When the pressure in a cell is larger than P vap + ∆P vap /2, the fluid is supposed to be pure liquid, the entire cell is occupied by liquid, and its density ρ l is calculated by the Tait equation (Knapp et al. , 1970).
where P T ref is the pressure at the domain outlet, ρ ref is the liquid density, and for water, P 0 = 3 10 8 Pa and the exponent is n = 7. If the pressure is lower than P vap − ∆P vap /2, the cell is full of vapor and its density ρ v is given by the perfect gas law (isotherm approach), where R = 462 J/K/kg for water vapor. In the other situations, the cell is occupied by a liquid/vapor mixture, which is considered as one single fluid with a variable density ρ. This one is directly related to the void fraction α = (ρ(P ) − ρ l )/(ρ v − ρ l ) corresponding to the local density of the fluid.
To model the mixture state, the barotropic equation of state presents a smooth transition in the vapor pressure value neighborhood, in the range P vap ± ∆P vap /2. In direct relation with the range ∆P vap , this equation is characterized mainly by its maximum slope 1/C min 2 , where C 2 min = ∂P/∂ρ. C min can thus be interpreted as the minimum speed of sound in the mixture. Its calibration was done in previous studies (Coutier-Delgosha et al. , 2003b). The optimal value was found to be independent of the hydrodynamic conditions, and is about 1.5 m/s for cold water (20 o C), with P vap = 0.023 bar, and corresponding to ∆P vap ≈ 0.06 bar. These values are used here throughout the presented results.
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 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 (Merkle et al. , 1998). The momentum transfer between the phases is thus strongly linked to phase change.

Numerical resolution
To solve the time-dependant Reynolds-averaged Navier-Stokes equations associated with the barotropic equation of state presented here above, the numerical code applies, on 2D structured curvilinear-orthogonal meshes, the SIMPLE algorithm (Patankar, 1981)), modified to take into account the cavitation process. It uses an implicit method for the time-discretization, and the HLPA non-oscillatory second order convection scheme proposed by Zhu (1991). The numerical model is detailed in Coutier-Delgosha et al.

Turbulence model
In our previous studies (Coutier-Delgosha et al. , 2002, 2003b, either the k-ω model proposed by Wilcox (1998) or the k-ǫ RNG model presented by Yakhot et al. (1992) were applied to model cavitating flows. Results obtained have demonstrated that for both models, corrections of the influence of the vapor/liquid mixture compressibility on the turbulence should be taken into account to obtain the unsteady effects due to cavitation. For the present test case, the modified k-ǫ RNG model presented in Coutier-Delgosha et al. (2003a) is applied.

Geometry and boundary conditions
The main features of the geometry can be seen in figure 2. The precise description of the Venturi section is given as a list of coordinates for the lower wall of the Venturi (see table 1) and the upper wall of the Venturi (see table 2).
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. Details of the prescribed values are given in section 3.3.

Grid
The computational grid is composed of 160×50 orthogonal cells (figure 2). 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 fifty grid points are used in this direction to model the 45 mm long mean cavity obtained by numerical calculations in the case σ = 2.4 where σ is defined in (3)(results presented hereafter). In the other direction, a contraction is also applied close to the walls, to obtain at the first grid point the nondimensional parameter y + of the boundary layer varying between 30 and 100, and to use standard laws of the wall. 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 contain about thirty cells across their thickness.

Initial conditions
To start unsteady calculations, the following numerical procedure is applied: first of all, a stationary step is carried out, with an outlet pressure value 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 σ defined by Vapor appears during the pressure decrease. The cavitation number is then kept constant throughout the computation.

Calculations
Calculations are performed with non-dimensional variables based on the following reference parameters :  are almost periodic, and their frequency as well as the mean cavity length depends on the cavitation parameter, σ.
As the final cavity obtained is fundamentally unstable, it cannot be characterized by its final shape or the final void fraction distribution. The comparisons are thus based on the transient evolution of the cavitating flow. This evolution can be defined at each time by the vapor quantity present in the domain or by the cavity shape (length, volume). We propose to focus on the vapor volume oscillations. For each computation, the timeaveraged vapor volume and its standard deviation are indicated in table 3. The cavitation cycle frequency can be calculated by using a FFT analysis of the inlet pressure signal and is also given in  The transient evolution observed for σ=2.4 during the unsteady calculation is presented in figure 3. Figure 3(a) illustrates at a given time and for each cross section of the Venturi type duct the value of the minimal density present in the section. 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 easily calculated.
Moreover, it also shows the minimum density, i.e. the maximum void ratio in each section. The two other curves shown in figures 3(b) and 3(c) represent respectively the total vapor volume and the inlet pressure evolutions.

Comparison with experiments 4.1 Overall behavior
First, the evolution of the cavity shape at a given cavitation number is compared with pictures obtained from experiments. Two phase-averaged cavitation cycles are presented in figure 4. The right one results from experimental visualizations: Video frames acquired during a 100 ns exposure time under Laser sheet light are identified and digitized in 256 grey levels. A sampling technique is applied to classify them in 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 drawing a sequence of phase-averaged images, from an initial data set of 300 frames. The left part of the figure 4 corresponds to the same sequence obtained by numerical simulation (calculation duration equal to 20 T ref , i.e. about 30 cycles). The same sampling technique is applied: 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.
From the results of table 3 concerning the effect of the cavitation number, a comparison can be proposed with results reported by Stutz & Reboud (2000): the frequency of the self-oscillation behavior is drawn with respect to the ratio V ref /L cav . In both the numerical simulation and the experiments L cav is chosen as the maximum length of the attached cavity.

Flow field inside the sheet of cavitation
Local comparisons are proposed in the case σ = 2.4 with experimental data obtained by double optical probes measurements. This technique and the results are presented in detail in Stutz & Reboud (1997, 2000. This is an intrusive technique, which allows measurements of the local void ratio and the velocities of the two-phase structures inside the cavitation sheet. Four data profiles located respectively at x = 1.410 −2 m, 3.110 −2 m, 4.910 −2 m, and 6.510 −2 m, are available. The time-averaged and standard deviation values of the velocity u and the void ratio α are compared along the four profiles in figure 6.
The three main features that should be obtained are: • presence of the re-entrant jet, characterized by negative or zero mean values of the velocity u close to the wall.
• the rather low mean void ratio observed in the main part of the cavitation sheet: it does not exceed 25%, excepted in the upstream end of the cavity.
• the general high level of velocity and void fraction fluctuations of the same order of magnitude than the mean values.
These three characteristics of the flow are strongly related to the overall unsteady oscillation of the cavitation sheet. So an accurate numerical simulation of these features is linked to the capability of the model to predict the unsteadiness of the cavitating flow.