Influence of the Cavitation Model on the Simulation of Cloud Cavitation on 2D Foil Section

For numerical simulations of cavitating flows, many physical models are currently used. One approach is the void fraction transport equation-based model including source terms for vaporization and condensation processes. Various source terms have been proposed by different researchers. However, they have been tested only in different flow configurations, which make direct comparisons between the results difficult. A comparative study, based on the expression of the source terms as a function of the pressure, is presented in the present paper. This analytical approach demonstrates a large resemblance between the models, and it also clarifies the influence of the model parameters on the vaporization and condensation terms and, therefore, on the cavity shape and behavior. Some of the models were also tested using a 2D CFD code in configurations of cavitation on two-dimensional foil sections. Void fraction distributions and frequency of the cavity oscillations were compared to existing experimental measurements. These numerical results confirm the analytical study.


INTRODUCTION
In liquid flows, cavitation generally occurs if the pressure drops below the vapor pressure; it can be observed in a wide variety of propulsion and power systems like pumps, nozzles, injectors, marine propellers, and hydrofoils. The cavitation development may be the origin of several negative effects, such as noise, vibrations, performance alterations, erosion, and structural damages.
Numerical studies and simulations of cavitation have been pursued for years, but it is still a very difficult and challenging task to predict such complex unsteady and twophase flows with an acceptable accuracy.
In the first category, the cavity region is generally assumed to have a constant pressure equal to the vapor pressure of the corresponding liquid and the computations are performed only for the liquid phase.
In the second category, the single-fluid modeling approach is employed for both phases. Mass and momentum transfers between the two phases are managed either by a barotropic state law or by a void fraction transport equation.
In the present study, the attention is focused on this second category, which has demonstrated in the last 20 years its capability to reproduce the main features of steady and unsteady cavitating flows [3]. Various source terms for the void fraction transport equation have been proposed in the literature, but most of the time, they have been applied to different flow configurations. Direct comparisons between the models are thus difficult to handle. However, the test case proposed for the Fifth International Symposium on Cavitation (CAV2003) in Osaka has revealed some similarities between results obtained with different models [8][9][10][11].
The present work consists of a comparative study between the different vaporization and condensation terms proposed for the void fraction transport equation. The ability of each one to reproduce the behavior of cavitating unsteady turbulent flow is discussed on the basis of the  analytical expression of the source terms and 2D numerical simulations.
In Section 1 of this paper, the expressions of the vaporization and condensation terms proposed in the literature are given. Then, both source terms are expressed as a function of the pressure and compared in Section 2. In Section 3, some of the models are evaluated in two configurations of two-dimensional foil sections. 2D simulations are performed and the results are compared to existing experimental data. Comparisons are based on the maximum length of the attached cavity, the frequency of its oscillations, and the void fraction distribution within the sheet cavity.
In Section 4, the influence of some parameters of the vaporization and condensation terms, such as the empirical constants and the exponents of P − P v and α l , is discussed. The objective is to assess their respective effect on the shape and behavior of the sheet cavity.

CAVITATION MODELING
The attention is focused in the following analysis on the single-fluid modeling approach associated with the Reynolds-averaged Navier-Stokes (RANS) equations.
The general balance equations for the single-fluid modeling can be found in [12]. In this previous paper, the general assumptions common to nearly all the current cavitation models have been clarified: (i) no slip between the vapor and the liquid phases, which leads to identical velocities in both phases, (ii) equal local pressure in the vapor and in the liquid. This leads to the following simplified four equations:  (equations derived from the mass balance, where A is the mass transfer between the two phases) (equations derived from the momentum balance, where B is the momentum transfer between the two phases). It has been shown in [12] that a third assumption is inherent to all current cavitation models: (4) is never resolved, which implicitly leads to the following expression for momentum transfer term B: The first term in the expression of B corresponds to a momentum transfer directly linked to the mass transfer during vaporization and condensation processes: as a matter of fact, let us consider inside a volume V an elementary mass δm of liquid moving at velocity u. If these liquid particles vaporize, δm is transferred from the liquid phase to the vapor one, and the vapor phase also gains the momentum δm · u.
The second term in the expression of B results from the momentum fluxes due to the inertial effects, that is, the stresses applied by one phase on the other when this last one accelerates or decelerates. The present form of this term means that if a volume of liquid is accelerated, the resistance applied by the adjacent volumes of vapor is transmitted exclusively by the local stress tensor. The forces between vapor and liquid are thus only due to the pressure and the classic viscous stresses, as in a single-phase flow.
So, in practice, (1) and (3), which are the classical RANS equations for a single-phase flow, are coupled either via  supplementary void fraction transport (2) or simply with a barotropic state law that governs the density evolution according to the local pressure variation. In the second case, the density is directly linked to the pressure, which prevents from taking into account any vaporization delay or separate treatments of the vaporization and condensation sequences. Moreover, coupling the barotropic state law with the presence of noncondensable gas and/or two-phase inlet flow is quite complicated. Conversely, these limitations can be avoided by using the void fraction transport equation.
In the present study, several models based on the void fraction transport equation are compared. They are characterized by different expressions of term A in (2). In the literature, A is usually divided into two terms: A =ṁ + +ṁ − whereṁ + andṁ − stand for the vaporization (vapor production) and condensation (vapor destruction) processes, respectively. Table 1 presents an overview of several forms oḟ m + andṁ − proposed by different researchers. All the models use empirical constants C prod and C dest to calibrate the mass transfers.

COMPARISON BETWEEN THE MODELS
Most of the source terms depend mainly on the difference between the local pressure and the vapor pressure P − P v . Thus, the following comparison between the models is based on the expression of the source terms as a function of P − P v .
However, the void fraction α usually also appears in the expression of the source terms. To express them as a function of P − P v only, the barotropic state law of Delannoy is used (Figure 4). Such a state law enables to link α and P − P v , and thus to suppress α in expressions ofṁ + andṁ − . It is clear that the choice of this particular state law is arbitrary. However, it enables to obtain here expressions depending only on P − P v , which is necessary for comparison. It must be also noticed that this transformation has only minor influence on the differences between the models: using another state law would lead to slightly different results, but the agreement or disagreement between the models would remain nearly the same.
In literature, empirical factors are determined through numerical/experimental results and are adjusted for different geometries and different flow conditions; C prod = 9 * 10 5 and C dest = 3 * 10 4 for the Kunz model [13], C prod = 0.02 and C dest = 0.01 for the Singhal model [7], and C prod = 0.1 and C dest = 0.1 for the Saito model [8].  The empirical factors are adjusted here to obtain the same maximum value for the source terms, in order to make the comparison of the models easier. The empirical factors adopted in this study include constants parameters of the source terms (ρ l , ρ v , U ∞ , . . .); so, it is not possible to compare our empirical constants with those proposed by the authors.
The evolution of the source terms according to P − P v is presented in Figures 1 and 2. Large similarities between both vaporization and condensation terms are clear. Some of the models (Merkle and Reboud or Sauer and Singhal) give identical expressions ofṁ + andṁ − and should thus result in exactly the same prediction of cavitating flows in numerical simulations; so far they would be associated to the same numerical model.
It can be observed that the transfer termsṁ + andṁ − in Table 1, although they seem quite different at first sight, often involve similar expression of P −P v and α, which explains the resemblance of the charts in Figures 1 and 2. For example, m + is usually a combination of P − P v and α l , with various exponents. Several constants are systematically included iṅ m + andṁ − , which leads to various analytical expressions, but does not change fundamentally the model. Adjusting them to obtain here the same maximal value ofṁ + andṁ − makes this point clearer. We also notice that the slopes of the source terms are different; the more the slope of the source term of vaporization (condensation) is increased, the more the vaporization (condensation) is.
In the next section, a numerical study in two configurations of cavitation on two-dimensional foil sections is performed. Some of the previous models were tested and the results are compared to existing experimental data.

NUMERICAL STUDY
Two configurations of cavitating flow on 2D foil sections are considered in this section: a convex foil with an upper flat surface and an NACA 66 foil. Experimental investigations have been performed previously in both cases, and data are used in the present study for comparison with numerical simulations.
Calculations of the cavitating flow around the two foil sections are performed with the code IZ, which was developed previously in the LEGI laboratory (Grenoble, France), with the support of the CNES (French space agency). Details concerning the numerical model can be found in [17], only the main features are given hereafter.

Reboud and Stutz [5]
Kunz et al. [13,14]ṁ Singhal et al. [7]ṁ Schnerr and Sauer [16]ṁ  The code is based on a finite volume discretization on curvilinear two-dimensional orthogonal mesh. The numerical resolution consists of a pressure correction method derived from the SIMPLE algorithm. Each physical time step consists of successive iterations which march the solution toward convergence. A finite volume discretization and a second order implicit time integration scheme are applied as follows:   where Φ is 1 in the case of the mass equation, and the velocity components u or v in the case of the momentum equations. A modified k-ε RNG turbulence model is applied: the turbulent viscosity is ; n is set to the value of 10 proposed by Coutier-Delgosha et al. in 2003 [18].
The classical boundary conditions for incompressible flows are applied: imposed inlet velocity and fixed outlet pressure. A C-type orthogonal mesh is used in the case of the NACA 66 foil, and an H-type orthogonal mesh is applied in the case of the convex foil ( Figure 3).
A 630 × 50 C-type orthogonal mesh for the NACA66 foil and 140 × 70 orthogonal cells for the convex foil.

Study of the sheet cavity unsteady behavior
The geometry considered in this section is the NACA 66 foil. Its chord is 150 mm and the angle of attack is 8 • . The reference velocity at the inlet of the test section is V ref = 5.33 m/s and the cavitation number is σ = 1.3.
A calculation based on the barotropic state law of Delannoy and Kueny [6] is first performed to be compared to the results that will be obtained with the various models based on the void fraction transport equation. In this case, (2) is not solved, but the state law presented in Figure 4 is coupled to (1) and (3) to govern the flow density variations according to the local pressure variations.
For a pressure much higher than the vapor pressure P v , the flow is composed of pure liquid and the Tait state law is applied: ρ/ρ ref = n P + P 0 /P ref + P 0 with P 0 = 3108 Pa and n = 7 for water.
For a pressure much lower than P v , the fluid is locally completely vaporized and the density is governed by the perfect gas law: P/ρ = cste.
These two low compressible configurations are joined in the vapor pressure neighborhood by the central part of the chart, whose high slope models the high compressibility of the liquid/vapor. C 2 min = ∂P/∂ρ. C min is the minimum speed of sound in the mixture.
The minimum density (indicated by the color) in each section (denoted by the position X/L ref in ordinate) as a function of the time (in abscissa).
Results of this calculation are presented in Figures 5, 6, and 7. An unsteady periodical behavior of the sheet cavity is obtained, including vapor cloud shedding. The maximum length of the attached cavity can be estimated to l max /c = 0.8. The cavity oscillation frequency is about 11 Hz (Figure 7), which gives a Strouhal number based on the maximum length of the attached cavity equal to 0.25.
During this cycle, a sheet cavity is formed, which successively grows and breaks off. The detachment of the rear part of the cavity is due to a reentrant jet that flows upstream close to the foil surface. The vapor cloud is convected downstream and finally collapses, while the residual cavity starts to grow again. Three transport equation-based models with source terms of Reboud and Stutz [5], Kunz et al. [13,19], and Singhal et al. [7] are now tested, and the results are presented in Figure 19. The maximum length of the attached cavity, the oscillation frequency, the Strouhal number, and the cavity behavior during one cycle are indicated for each model. The empirical factors were adjusted to obtain the same maximum value for the source terms, as in the previous section. Figure 19 shows a large resemblance between the four tested models; the unsteady behavior of the sheet cavity is similar in all cases: development of a cavity up to the maximum length (about 80% of the chord), cavity breakoff, cloud detachment, and growth of the residual cavity.
Sheet cavity oscillation frequencies are also all similar, leading to Strouhal numbers Str based on the maximum cavity length comprised between 0.25 and 0.28. These results are in close agreement with the experimental measurement, which gave Str = 0.28 in the configuration of a slightly smaller cavity (L/L ref = 0.6). These results confirm the similarities between the three models based on the void fraction transport equation. They also suggest that using the barotropic state law instead of (2) leads to very similar predictions of the sheet cavity behavior.
The colors indicate the fluid density, white for pure liquid and from red to blue when the void ratio increases.

Study of the void fraction distribution
The second foil geometry is a convex foil characterized by a flat upper surface. Its chord is 150 mm and the present angle of attack is 4 • 7. The reference velocity at the inlet of the test section is V ref = 6 m/s. Calculations of the cavitating flow on the foil are performed for two sheet cavity lengths: L/c = 0.9 and L/c = 0.5, respectively. In both simulations, an unsteady behavior of the cavity is obtained, as previously in the configuration of the NACA66 foil section. Discussion is focussed hereafter on the study of the time-averaged two-phase flow structure of the liquid/vapor mixture inside the cavity.
An X-ray absorption device was applied previously to investigate the local volume fraction of the vapor phase inside the sheet cavity. The experiments were carried out by Coutier-Delgosha et al. [17] in the scope of collaboration between the ENSTA laboratory and the French Atomic Energy commission (CEA).
From these experiments, the instantaneous and timeaveraged distribution of void fraction inside the cavity was obtained, with a space resolution of 6 mm in the flow direction and 3 mm in the vertical direction. In the present study, the time-averaged experimental void fraction measurements at three positions (X = 15 mm, 55 mm, and 95 mm downstream from the leading edge) are considered (Figure 8).
The first studied configuration is V ref = 6 m/s and L/c = 0.9. Void fractions profiles obtained with the barotropic state law of Delannoy and the void fraction transport equationbased models of Kunz and Reboud are compared to the experimental data in Figure 9. Note that the empirical constants in the models have been adjusted as in the previous sections.
The void fraction evolutions obtained with the barotropic state law and the Reboud model are similar. They also exhibit a nice agreement with the X-ray measurements. Conversely, The use of the Kunz model leads to smaller void fractions in all the cavity, although its maximum length is close to L/c = 0.9, as for the other calculations and the experimental data. It suggests that this particular model underestimates the flow vaporization and/or overestimates the flow condensation. Based on Figure 1, the vaporization  Figure 10.
The void fraction profiles obtained with the different cavitation models are very similar, and they are generally in close agreement with the X-ray measurements. The only discrepancy still concerns the Kunz   predictions of this model at stations 2 and 3 are very close to the experimental data, which confirms that the mean length of the cavity is correctly estimated. The nice agreement obtained between all the present data confirms the similarities observed previously in the expression of the vaporization and condensation terms. However, some slight differences between the models suggest that some constants inṁ + andṁ − , such as the exponents of P − P v or α, may have some significant influence on the predicted flow structure. The objective in the next section is to study this influence in order to assess the particular effect of each variable.

INFLUENCE OF THE MODEL PARAMETERS
Four parameters are considered: the empirical constants C prod , C dest , and the exponents of P − P v and α l .

Influence of the empirical constants C prod and C dest
Two calculations around the NACA 66 foil section were performed with the Kunz model and several values of C prod and C dest were applied in the source terms. The results are presented in Figures 11 and 12. When C prod is increased, more vapor is created in the upstream part of the sheet cavity, so the maximum void fraction increases ( Figure 12). Consequently, the inertia of the cavity increases and the oscillation frequency f decreases. Conversely, it can be noticed that C prod has nearly no influence on the cavity thickness. The increase of C dest causes more intense condensation of vapor; so the maximum void fraction decreases, the inertia of the cavity decreases, and the oscillation frequency increases.
One calculation around the convex foil section was performed with the Kunz model and with 0.5 * C dest instead of C dest . It can be observed in Figure 13 that a bigger sheet cavity is obtained. The maximum void fraction has increased at the leading edge (X = 15 mm). This confirms that the Kunz model overestimates the flow condensation.
We also obtain a correct void fraction in the upstream part of the cavity without deteriorating the downstream one because of the form of the term of condensation; this term has much impact on the upstream part of the cavity, which is the zone of lower pressure.

Influence of the exponent of P − P v
One calculation around the NACA 66 foil section was performed with the Kunz model and a different exponent of P − P v in the expression of the vaporization term: the exponent is set to 1.1, instead of 1. The empirical factors were adjusted as in the previous sections. This modification results in a diminution of the slope of the vaporization source term (Figure 14), which suggests that less vapor may be obtained in the sheet cavity. However, it can be observed in Figure 15 that increasing this parameter seems to result in a bigger sheet cavity. Indeed, drawing the void fraction profiles at stations x/c = 0.1 and 0.5 ( Figure 16) enables to confirm that the cavity is thicker, while the maximum void fraction has decreased.

Influence of the exponent of α l
Calculations around the NACA 66 foil section are performed with the Kunz model and two different values of the exponent of α l in the expression of the vaporization term. The empirical factors were adjusted as in the previous sections.
A clear influence of this parameter can be observed in Figure 17: the more the exponent is increased the higher is the slope of the vaporization term, which suggests that more vapor may be obtained in the sheet cavity.
Drawing the void fraction profiles at stations x/c = 0.1 and 0.5 ( Figure 18) enables to confirm that the void fraction has increased in the whole height of the sheet cavity for x/c = 0.5, while the maximum void fraction has decreased for x/c = 0.1, and the cavity thickness has not changed much. It suggests that the vaporization inception due to the fluid inertia, which occurs at the foil leading edge, is slightly reduced (see station x/c = 0.1), whereas the vapor bubble expansion, which mainly occurs in the upstream part of the sheet cavity, is more intense when the exponent of α l is increased.

CONCLUSIONS
The present work is a contribution to the physical modelling and numerical simulation of cavitating flows. An analytical approach was first performed in order to compare transport equation models proposed by various authors. The resemblance between most of the models was observed.
A numerical study using a 2D CFD code was also performed, and several models were tested and compared. Sheet cavity lengths, oscillation frequencies, and cavity behaviors obtained with the different models are found in close agreement.
The comparison between the numerical and the experimental void fraction distributions in a second configuration of 2D foil section demonstrates also the resemblance between the models, and thus it confirms the results of the analytical study.
The influence on the results of some arbitrary constants used in the models was also studied. It was shown that modifying these parameter enables to change the cavity shape and structure. The objective of the work is now to link these different constants with physical considerations in order to construct vaporization and condensation terms which would take into account more physics of the vaporization and condensation processes.

c:
Chord of the foil C min : Minimum speed of sound in the medium C dest : Empirical constant in the condensation term Liquid volume fraction α l = 1 − α γ: Surface tension ρ: Mixture density ρ v , ρ l : Vapour and liquid densities ρ ref : Reference density (outlet liquid density) C μ : Turbulence model constant σ: Cavitation parameter, Stress tensor in the liquid/vapor mixture.