Numerical investigations on cavitation intensity for 3D homogeneous unsteady viscous flows

The cavitation erosion remains an industrial issue. In this paper, we deal with the cavitation intensity which can be described as the aggressiveness - or erosive capacity - of a cavitating flow. The estimation of this intensity is a challenging problem both in terms of modelling the cavitating flow and predicting the erosion due to cavitation. For this purpose, a model was proposed to estimate cavitation intensity from 3D unsteady cavitating flow simulations. An intensity model based on pressure and void fraction derivatives was developped and applied to a NACA 65012 hydrofoil tested at LMH-EPFL (École Polytechnique Fédérale de Lausanne) [1]. 2D and 3D unsteady cavitating simulations were performed using a homogeneous model with void fraction transport equation included in Code_Saturne with cavitating module [2]. The article presents a description of the numerical code and the physical approach considered. Comparisons between 2D and 3D simulations, as well as between numerical and experimental results obtained by pitting tests, are analyzed in the paper.


Introduction
The prediction of cavitation and material erosion remains an issue for hydraulic machinery manufacturers and users. High flow velocities cause regions of low pressure where vapour structures are generated. These cavitating structures collapse rapidly after reaching a region of higher pressure and are able to cause performance loss, vibration and can damage the material.  The main problem of simulating cavitation erosion is the fact that it deals with several length and time scales phenomena (see Figure 1) and involves both fluid and mechanical behavior. The cavitation intensity -or cavitation aggressivness -represents the mechanical loading imposed by the cavitating flow to the material. Erosion, defined as mass loss, can then be deduced from this quantity using methods as in [3] and [4].
In this paper, a non-exhaustive state of the art will be done. Then, we will describe our cavitation intensity model and we will apply it to a flow around a NACA hydrofoil.

State of the art
Many numerical studies have been carried out to predict cavitation erosion. Table 1 summarizes different approaches proposed by some authors. Note that a homogenous model is most often used and that the energy equation is not taken into account except for [5]. 3. The cavitation intensity model Our cavitation intensity model is based on the idea of Fortes-Patella et al. [3], who proposed a cavitation erosion model where "potential" energy variations of the cavitation structures are considered as the main factor that generates erosion. This approach is applied in the present study as a sub-mesh model (i.e. post-processing model) using U-RANS calculation with Code Saturne with cavitating module.
3.1. Code Saturne with cavitating module main features Code Saturne is a free open source CFD software developped by EDF [12]. It carries out 2-or 3-D, steady or unsteady, incompressible, laminar or turbulent simulations on any sort of mesh.
It is based on a colocalised finite volume method. Some modules can be added on Code Saturne in order to describe additionnal phenomena, such as compressible, rotor-stator or cavitation phenomena. Code Saturne with cavitating module enables the mean resolution of a homogeneous mixture model with void fraction (α) transport. Pure phases have constant properties (density, ρ l/v and dynamic viscosity µ l/v ) following the relations : It is assumed that the mixture dynamic is ruled by the Navier-Stokes equations (mass (1) and momentum (2) conservation) with a void fraction transport equation (see equation (3)) : with Γ v the vaporisation source term.
A standard k-ε turbulent model with Reboud correction [14] is used. The resolution scheme is based on a co-located fractional step scheme, which is associated with the SIMPLEC-type algorithm (see [12] for more details).

Sub-mesh modelisation 3.2.1. Energy approach
Based on the idea of Pereira [15], we can calculate a "potential" energy (E pot ) of vapour structures (see equation (4)) : Then, we can deduce the potential power (P pot ) of those structures for each cell and separate it in two parts (see equation (5)). The first one takes into account the contribution of the void fraction derivative, and the second one deals with the pressure derivative influence.
V cell is the volume of a cell. We assume that the vapour structure is aggressive if P pot > 0. Since we have, by definition, . Moreover, we know, from mass conservation that dρ/dt + ρ div(u) = 0. Then : . Fortes-Patella et al. [3] decided to ignore the pressure derivative part of the potential power since, as we will see in our application, ||P pot | α=cst || << ||P pot | p=cst ||. A comparison of this model [3] with the Nohmi's one [9] is presented in [16].
We will analyse in this paper the influence of both terms on the amplitude and location of the cavitation intensity.

Solid angle
On the basis of Krumenacker's work [6] and by using the analytic exact expression of the solid angle (Ω) for a planar triangle [17] (see equation (6)) , we can deduce the potential power applied on the material surface (P mat pot ) (see Figure 2 and equation (7)), which defines the quantity we will name the "instantaneous cavitation intensity". The use of the solid angle quantifies the distance and angle dependancies of the potential energy from the cell source to the surface.
On each element j of the hydrofoil surface and for all i fluid cells : 4. Application to a cavitating flow around a hydrofoil 4.1. General description Our prediction model has been applied to a NACA 65012 hydrofoil (chord length is 100 mm and span 150 mm) tested in the cavitation tunnel of the LMH-EPFL (École Polytechnique Fédérale de Lausanne) [1] (see Figure 3).  Experimental conditions tested by Pereira [15] and simulated in this paper are summarized in the Table 2, where C ref describes the mean axial flow velocity at the inlet of the tunnel, i the attack angle of the hydrofoil, σ the inlet cavitation number (see equation (8)), l the cavitation sheet length and L the hydrofoil chord.
with p in the inlet pressure.  Influence tests concerning mesh size and time parameters will be dealt with in subsequent work. The y + dimensionless wall distance (y + = u * y/ν with u * the friction velocity, y the wall distance and ν the kinematic viscosity) varies between 50 and 70 under non-cavitating conditions. Inlet flow velocity (C ref ) and outlet pressure (p out ) are imposed (see Figure 3 for others boundary conditions).

Hydrodynamic results
In order to validate the cavitating flow behaviour, we will first calibrate the cavitation sheet length (by iteration on the outlet pressure) and then compare the cavitating structures shedding frequency of the experimental results with the simulated one. A void fraction isosurface at 10% is used to calibrate the cavitating sheet length (see Figure  5).
Then, we check the shedding frequency by doing a Discrete Fourier Transform for the inlet pressure signal (see Figure 6).  In this example, the first natural frequency (58 Hz) characterizes the periodic transverse oscillation (oscillation in the span direction). The second one (118 Hz) characterizes the shedding frequency (f c ). Harmonics are also present. Finally, we can compare our results with experimental ones in terms of shedding frequency (see Figure 7).  Shedding frequency fonction of reduced frequency for experimental (with linear regression at S t = 0.3) [15], 2D and 3D simulations results.
2D and 3D simulations based on the same boundary conditions have different shedding frequencies. Even if the cavitating sheet length is the same (l/L = 40%) the 3D dynamic is quite faster than the 2D one. One notes that we do not have exactly the same inlet σ for 2-D and 3-D cases.
In conclusion, 3D simulated cases are in good agreement with the experimental hydrodynamic behavior (for cavitation sheet length and shedding frequency) whereas 2D ones are a bit less accurate.

Cavitation intensity results
The case of the NACA 65012 hydrofoil at 6°and C ref = 15 m.s −1 will be developped here. The calculation duration is set up to 0.2 s and we keep the same transit time (0.05 s) before applying post-processing. This study can be extended to all the other simulated cases.
We first calculate the instantaneous potential power in the fluid (see Figure 8). One can note that P pot /V cell is higher at the cavitating sheet closure (where dα/dt reaches maximum values). One can note that the potential power coming from the void ratio derivative is about 5 times higher than the one coming from the pressure derivative. Then we use the solid angle value to evaluate the potential power on the hydrofoil surface. Figure 9 shows, for a given time, the instantaneous received surface power on the hydrofoil (P mat pot /∆S) and the iso-surface at 10% of void fraction. One can note that, in agreement with P pot /V cell , the maximum value of P mat pot /∆S is located at the cavitating sheet closure. We can finally add up all the received surface power (P mean pot /∆S) by each surface (see Figure  10) and divide the result by the number of time steps to have a mean loading (see equation 9), which can be used as a qualitative representation of the eroded region.
with N the number of time step considered (here N = 150.000).

Cavitation intensity analysis
Close to the tunnel wall, the flow seems to be twice less agressive. Even if we can't compare this statement with the experimental results (all samples are taken in the middle of the hydrofoil) we can physically explain this phenomenon : Near the wall, only the half domain can damage the surface compared to the median region where the fluid domain all around it can damage the surface. The reflection of the potential power due to the (plexiglass) side wall is not taken into account because of the lack of knowledge on the part of absorbed or transmitted potential power. Figure 11 shows the mean surface power function of the hydrofoil chord. A comparison is done with the experimental volume damage rate given by pitting tests (V d i.e. the deformed volume divided by the analyzed sample surface area and test duration) [3] for 2D and 3D simulations, using the same algorithm. Even if the 2D frequency behaviour is not well fitted with the experimental one, the maximum location of P mean pot /∆S between 2-D and 3-D simulations are nearly the same. Quantitavely, the 3-D simulation is a bit more erosive for the hydrofoil because of the transverse elements contribution (3-D effect). Figure 12(a) shows the mean cavitation intensity for different velocities. The maximum of P mean pot /∆S increases with the velocity but its location is nearly the same for each case. Figure 12(b) shows the relation between the maximum cavitation intensity value and the inlet velocity. We find P mean pot /∆S = 6.2 C ref 3 for the 3D case, which agrees with [18]. Indeed, we have  In comparison with experimental results, the maximum erosive contribution is shifted of 10% downstream. This difference could be explained by several reasons. The first one is the arbitrary 10% of void fraction taken to calibrate the cavitation sheet length and therefore impacts the choice of σ in the CFD. The second one is that we do not know how the experimental sheet length is found and what is the precision of this measurement. By taking cavitation sheet length of 50% (σ = 1.34), 3D simulation matches much better with experimental results (see Figure 11).

Conclusion
Based on the literature and on previous works carried out in the scope of scientific collaborations between the University of Grenoble and EDF R&D, a cavitation intensity model has been developed using Code Saturne with cavitating module. The model was applied to evaluate the aggressiveness of cavitating flows around a NACA hydrofoil. Comparisons between numerical and available experimental results allow the qualitative and quantitative validation of the proposed approach concerning the prediction of the flow unsteady behavior, of the location of erosion area and of the influence of flow velocity on the cavitation intensity. The comparative analyses of 2D and 3D numerical results indicated that 3D effects should be taken into account to obtain reliable quantitative evaluations of the potential power applied on the foil.
In further work, various tests will be done to evaluate the influence of numerical (mesh, time steps, convergence levels) and physical parameters (constants in the source terms of the cavitation model) on the results obtained. A local model (at the bubble scale) is in progress to better understand the microscopic phenomena and to improve our sub-mesh model. This approach will then be applied to evaluate the cavitating flow damage in a centrifugal pump [19].