1. 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 two-phase flows with an acceptable accuracy.

These studies can be put mainly into two categories: interface tracking methods [1, 2] and homogeneous equilibrium flow models [3–7].

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–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 and , is discussed. The objective is to assess their respective effect on the shape and behavior of the sheet cavity.

2. 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 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 of liquid moving at velocity . If these liquid particles vaporize, is transferred from the liquid phase to the vapor one, and the vapor phase also gains the momentum .

The second term in the expression of 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: where and stand for the vaporization (vapor production) and condensation (vapor destruction) processes, respectively. Table 1 presents an overview of several forms of and proposed by different researchers. All the models use empirical constants and to calibrate the mass transfers.

Table 1: Overview of some cavitation models existing in the literature.

3. Comparison between the Models

Most of the source terms depend mainly on the difference between the local pressure and the vapor pressure . Thus, the following comparison between the models is based on the expression of the source terms as a function of .

However, the void fraction αusually also appears in the expression of the source terms. To express them as a function of only, the barotropic state law of Delannoy is used (Figure 4). Such a state law enables to link αand , 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 , 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; and for the Kunz model [13], and for the Singhal model [7], and and 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 ; so, it is not possible to compare our empirical constants with those proposed by the authors.

The empirical factors have the following values: , for the Kunz model, , for the Merkle model, , for the Singhal model, , for the Visonneau model, , for the Saito model, and for the Saito model.

The evolution of the source terms according to 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.

Figure 1: Source terms for the vaporization process.

Figure 2: Source terms for the condensation process.

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 and α, which explains the resemblance of the charts in Figures 1 and 2. For example, is usually a combination of and , with various exponents. Several constants are systematically included in 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.

4. 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.

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 RNG turbulence model is applied: the turbulent viscosity is , where ; 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).

Figure 3: General view of the grid, flow is from the left.

Figure 4: Barotropic state law ; water 20^°C.

A 630  ×  50 C-type orthogonal mesh for the NACA66 foil and 140  ×  70 orthogonal cells for the convex foil.

4.1. 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  m/s and the cavitation number is .

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: with  Pa and 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: .

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. . is the minimum speed of sound in the mixture.

The minimum density (indicated by the color) in each section (denoted by the position 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 . 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.

Figure 5: Time evolution of the sheet cavity.

Figure 6: (a) Time evolution of the cavity shape during a single cycle of vapor shedding. The colors indicate the fluid density, white for pure liquid and from red to blue when the void ratio increases. (b) Velocity field and fluid density during cavity breakoff show a reentrant jet (picture 6 in (a)).

Figure 7: (a) Time evolution of the pressure coefficient at and (b) spectral analysis.

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 (). 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.

4.2. 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  m/s.

Calculations of the cavitating flow on the foil are performed for two sheet cavity lengths: and , 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 time-averaged 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 ( mm, 55 mm, and 95 mm downstream from the leading edge) are considered (Figure 8).

Figure 8: Position of the experimental void fraction measurements ( m/s, %).

The first studied configuration is  m/s and . Void fractions profiles obtained with the barotropic state law of Delannoy and the void fraction transport equation-based 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.

Figure 9: Comparison between numerical and experimental void fraction distributions (,  m/s).

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 , 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 sources for Kunz and Reboud are nearly identical; this result suggests that the Kunz condensation model is too large in magnitude and it overestimates the flow condensation.

The second flow configuration is  m/s and . The barotropic state law of Delannoy and the void fraction transport equation-based models of Kunz, Reboud, and Singhal are applied to these flow conditions, and the computed void fraction profiles are compared to the experimental data in Figure 10.

Figure 10: Comparison between numerical and experimental void fraction distributions (,  m/s).

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 model, which predicts a too low void fraction close to the leading edge. However, 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 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.

5. Influence of the Model Parameters

Four parameters are considered: the empirical constants , , and the exponents of and .

5.1. Influence of the Empirical Constants and

Two calculations around the NACA 66 foil section were performed with the Kunz model and several values of and were applied in the source terms. The results are presented in Figures 11 and 12.

Figure 11: Influence of the empirical constants on the oscillation frequency and the shape of the time-averaged cavity ().

Figure 12: Void fraction distributions (,  m/s).

When 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 has nearly no influence on the cavity thickness. The increase of 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 instead of . It can be observed in Figure 13 that a bigger sheet cavity is obtained. The maximum void fraction has increased at the leading edge ( mm). This confirms that the Kunz model overestimates the flow condensation.

Figure 13: Void Fraction distributions (,  m/s).

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.

5.2. Influence of the Exponent of

One calculation around the NACA 66 foil section was performed with the Kunz model and a different exponent of 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 and 0.5 (Figure 16) enables to confirm that the cavity is thicker, while the maximum void fraction has decreased.

Figure 14: Vaporization term for two values of the exponent of .

Figure 15: Influence of the exponent of on the behavior of the time-averaged cavity, .

Figure 16: Void Fraction distributions (,  m/s).

5.3. Influence of the Exponent of

Calculations around the NACA 66 foil section are performed with the Kunz model and two different values of the exponent of 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.

Figure 17: Vaporization term for two values of the exponent of .

Drawing the void fraction profiles at stations and 0.5 (Figure 18) enables to confirm that the void fraction has increased in the whole height of the sheet cavity for , while the maximum void fraction has decreased for , 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 ), whereas the vapor bubble expansion, which mainly occurs in the upstream part of the sheet cavity, is more intense when the exponent of is increased.

Figure 18: Void fraction distributions (,  m/s).

Figure 19: Numerical results obtained with different models (, ,  s).

6. 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.

Acknowledgments

The authors would like to express their appreciation for the continuous support of the Ecole Navale and the Ministry of Defense (France). Concerning the numerical part of this work, the authors wish to express their gratitude to the CNES (French space agency) for its continuous support.