Cavitation erosion prediction based on a multi-scale method

In hydraulic turbines, cavitation has unwanted consequences such as vibrations, noise, or material degradation. Thus, developping a method which can anticipate this phenomenon with good accuracy is vital from the designer’s perspective. This work aims at predicting two aspects of cavitation erosion: determining areas most endangered, and quantifying cavitation erosion in a relative way.


Introduction
Cavitation is characterized by the apparition of bubbles which have an erosive power when they implode. This is damaging for hydraulic machinery and threatens the lifespan of every installation. Thus, predicting the impact of cavitation is justified for both technological and economical reasons. A cavitation erosion prediction model is implemented as a post-processing phase of a flow simulation with a commercial solver.
Numerous authors have studied cavitation erosion and its prediction in a variety of ways: experimental, statistical or numerical approaches have been developed [1,2,3,4]. Numerical models have the advantage to be based on the resolution of fluid dynamics equations rather than on models built thanks to observations. Furthermore, they do not need the set up of a new physical test case for each analysis. Generally, numerical cavitation erosion prediction models can be classified in two categories: either centered on a bubble and its evolution, or on a population of bubbles. An example of the first category is proposed by L. Krumenacker [5], who is concerned by the damage caused by the implosion of a single bubble. The results obtained by an acoustic approach are then extented to a larger number of bubbles. This method requires the resolution of Keller's equation [6] coupled with calculations using numerous parameters, such as the bubble radius or its surface tension. The other category of cavitation erosion modeling approaches avoids consideration of individual bubbles. C. Leclercq's model [7], based on R. Fortes-Patella's work [8] concerning the potential energy of vapor structures, only depends on the pressure, density and velocity fields obtained with URANS simulations.
In this work, an adaptation of Leclercq's model for a post-processing based on a Python script is proposed in order to determine the most endangered areas, and to quantify cavitation erosion in a relative way.

Method
The approach is based on a multi-scale method: at the macroscopic level, the cavitating flow is simulated, and at the microscopic scale, the impact of bubble implosion on material surfaces is evaluated. First, URANS simulations for cavitating flows are performed with CFX to get the needed pressure, density and velocity fields. From the CFD results, the cavitation power ψ, a quantity representative of erosion, and based on the mixture density, the pressure and the water velocity divergence, is calculated at each fluid point. Next, a Python-based script, using the solid angle method [9], calculates χ, the fraction of ψ that reaches the solid surface. This is computed for each point of the surface and each timestep. Finally, the erosion intensity coefficient I cav is calculated at each point of the surface, based on the time averaging of χ, using another Python script. This allows a visualization through colormap highlighting of critical areas in CFD-Post.
3. Macroscopic scale 3.1. Two-phase flow in CFX In CFX, there are two approaches to model multiphase flows: the homogeneous and inhomogeneous models. With the eulerian-eulerian approach, the two phases of the flow are continuous and can interpenetrate each other. Thus, the notion of volumetric fraction is used to model the quantity of the vapor and liquid phases in the mixture. According to the software documentation [10], the homogeneous model is based on the hypothesis that the two phases share the same velocity field, while the inhomogeneous model accounts for different velocity fields for each phase. In the homogeneous model, the Navier Stokes equations are solved for the mixture, while their resolution is done independently for each phase in the inhomogeneous model. The equations of both models are detailed in Stenmark's work [11].

Multiphase comparison
In order to determine which model should be used in the URANS simulations for the cavitating flow, a study was carried out on a NACA66 foil, the same used by Leroux [12] in his experimental work. The geometry is shown in figure 1 and the boundary conditions are in table 1.  A time step of 0.01s is used, corresponding to a RMS Courant number of 234. Furthermore, a SAS-SST turbulent model is used, and a 10 −6 RMS residual level is fixed as the convergence criterion. Kubota's model is adopted to determine the cavitating flow [13]. The Rayleigh Plesset mean bubble diameter is set at 2 · 10 −6 m, the saturation pressure p sat at 3169P a and the cavitation rate under-relaxation factor at 0.1. Numerical predictions of the pressure coefficient C + l are compared to experimental results in figure 2. Figure 2 shows that the inhomogeneous model gives better results concerning the shedding frequency of the cavitation sheet and the values of the pressure coefficient C + l . The formation and evolution of the cavitation sheet is presented in figure 3 to validate the model. CFX calculations with the inhomogeneous model give a good prediction of the development of the cavitation sheet and its shedding mechanism, even if a slight offset in time could be observed in comparison with Leroux's results. Thus, this method will be used in the following work.   Figure 3: Cavitation sheet evolution -a) Leroux's experimental (blue background) and numerical (white background) void ratio results [12], b) CFX void ratio results In this work, a 1.4 million element mesh is used, refined on the suction side of the foil and at the trailing edge. URANS SAS-SST simulations are carried out with a 0.01s timestep over a 0.6s period, which corresponds to a RMS Courant number of 742.

Calculation of the cavitation power
In Leclercq's model [7], the cavitation power at a given fluid point is defined as: where ρ g is the vapor density, ρ l the liquid density, ρ m the mixture density, p is the local pressure, p sat the vapor pressure, u the fluid local velocity and V the cell volume associated to the considered fluid point. The fluid velocity is defined as: where u g is the local vapor phase velocity and u l the local liquid phase velocity. Using a CFD-Post expression, the cavitation power ψ is calculated at each fluid point. Then, only the values of ψ V above a given threshold, Θ, are exported for the post-processing. This filtering is based on the assumption that the material is not endangered by low cavitation power values. Thus, this filter is a way to focus only on the areas of the solid surface that are exposed to the highest level of cavitation power. More details on Θ are given in section 4.2.

Calculation of the cavitation power per material surface area
In Leclercq's model, the cavitation power per material surface area χ is the cavitation power applied on a given material surface region, indexed j. To determine which fraction of ψ will reach the surface, the solid angle method is used and the solid surface has to be meshed with triangles.
For each timestep, and each triangle j, χ j is calculated as: where A j is the area of the triangle j of interest and Ω ij the solid angle associated with the triangle j and a fluid point i. Nearby fluid points i are cavitation sources which could directly reach the triangle j.
For the post-processing step, working with matrices rather than with equation (3) is advantageous for an efficient formulation in Python. The linearized formulation used is: with ∀j ∈ 1; n t , λ j = 1/A j where the "×" operator is a term by term product.
[Ω] is a n t × n f matrix composed of the solid angle values associated with every fluid-solid couple, [ψ] is a n f × 1 vector with the cavitation power values associated with every fluid point. n t and n f are the number of triangles in the material surface mesh and the number of fluid points respectively.
The results are shown in figure 4, where two different thresholds are used. As expected, using a high threshold value leads to the highlighting of areas under stronger erosion power. This is due to the fact that less cavitation sources are taken into account. In the following developments, the threshold is set at 5 · 10 6 W/m 3 , or a thousand times lower than the maximum value of ψ V in the domain.

Calculation of the cavitation erosion intensity
The coefficient I cav describing the time-averaged cavitation power per material surface area for N timesteps is defined as: This leads to the results presented in figure 5. One of the highlighted areas from the results of our approach is located at the end of the cavitation sheet, at about half the chord length along the suction side of the foil. This is in accordance with experimental observations and Leclercq's results. However, the CFX results predict erosion near the trailing edge and leading edge as well which does not correspond to the experiments. The deviation is due to high values of ψ near de leading and the trailing edge. This could be attributed to numerical issues due to the high value of the Courant number encountered in the simulations of the macroscopic scale.

Conclusion
A method for cavitation erosion prediction based on URANS simulations and Python-based calculations was developed. This method shows reasonable agreement with experimental observations concerning the most severe regions of cavitation damage on an airfoil. However, this implementation has also led to the detection of regions of low risk of cavitation damage, such as the ones located at the leading and trailing edges. The work presented here has been done for a steady cavitation sheet. A further step would be to study the cavitation erosion with the same model and an unsteady cavitation sheet.