A modified density based cavitation model for time dependent turbulent cavitating flow computations

The predictive capability of transport equation-based cavitation models including the Kubota cavitation model (Model-1) and interfacial dynamics cavitation model (Model-2), is evaluated for the attached turbulent cavitating flows. In this study, the test problem is the unsteady cloud cavitating flows around a Clark-Y hydrofoil. Based on the evaluations of existing models, we iden-tified the differences between these two vaporization and condensation processes in the affected region, and provided a modified density based cavitation model (Model-3). The numerical results of the cavity shapes, velocity distributions and dynamics of the cavity oscillations were compared to existing experimental data. Compared with the other cavitation models, a significant improvement for the numerical results of unsteady cavitating flows has been obtained with the new model. Our study provides the information for further modeling development.

Cavitation plays an important role in the design and operation of fluid machinery and devices because it causes performance degradation, noise, vibration, and erosion. Cavitation involves complex phase-change dynamics, large density ratios between phases, and multiple time scales [1]. Navier-Stokes computations of turbulent cavitating flows have been progressively adopted because of advances in computational capabilities and understanding of the physics of these problems. However, these issues pose challenges with respect to accuracy, stability, efficiency and robustness because of the complex unsteady interaction associated with cavitation dynamics and turbulence.
Wang and Su [2] simulated the cavitating flows inside diesel injection nozzle holes using a two-fluid model. Another common approach in cavitation modeling is the use of homogeneous flow theory. The mixture density concept is introduced and a single set of mass and momentum equations is solved in this theory. The differences between the *Corresponding author (email: wangguoyu@bit.edu.cn) various models in this category mostly come from the relationship that defines the variable density field: Delannoy and Kueny [3] proposed a baroscopic state law that strongly links the mixture density to the static pressure, which describes the mixture density in the incompressible parts. Chen and Heister [4] derived a time and pressure dependent differential equation for density. Ventikos and Tzabiras [5] introduced the water-vapor state laws for this purpose. They considered the whole enthalpy along with the N-S equation. In a similar fashion, Edwards et al. [6] used the Sanchez-Lacombe equation in addition to the N-S equations, because most cavitating flows are isothermal.
Another popular approach is the transport equation-based model (TEM) [7][8][9][10][11][12]. The transport equation for either mass or volume fraction is solved, with appropriate source terms to regulate the mass between phases in the TEM. Different source terms have been proposed by different researchers. Kubota et al. [7] coupled the Rayleigh-Plesset equation to the flow solver and computed the void fraction based on the bubble radius. Merkle et al. [8] and Kunz et al. [9] em-ployed the artificial compressibility method with special attention given to the preconditioning formulation. Singhal et al. [10] and Cao et al. [11] considered the effect of noncondensable gas in their "full cavitation model". Senocak and Shyy [12] developed an interfacial dynamics-basedcavitation model (IDM), which offers direct interpretation of the empirical parameters in the existing transport-equationbased models.
Turbulence cavitating flow computations need to address cavitation modeling issues because, phenomenologically, cavitation often involves complex phase-change dynamics. Until now, these physical mechanisms were not well understood because of the complex cavitating flow structures. Different transport equation-based cavitation models from the literature are summarized below and a modified density based cavitation model using the unsteady characteristics of cavitating flow is developed. To enable comparison with available experimental data, the new model is evaluated for unsteady cavitating flow around a Clark-Y hydrofoil. This study provides the useful information for further modeling modifications.

Numerical methods in cavitation computations
The set of governing equations employed in this study consists of the conservative form of the Reynolds averaged Navier-Stokes equations, plus a volume fraction transport equation to account for the cavitation dynamics, and the turbulence closure.

Continuity and momentum equations
The Navier-Stokes equations in their conservative form governing a Newtonian fluid without body forces and heat transfers are presented below in Cartesian coordinates: where ρ m is the mixture density, u is the velocity, p is the pressure, μ l and μ t are the laminar and turbulent viscosity, and subscripts i, j, k are the directions of the axes. In this study, a modified RNG k-ε turbulence model is applied, and the turbulent viscosity is where n is set to the value of 10 proposed by Coutier-Delgosha et al. [13].

Evaluation of transport equation-based models (TEM)
A cavitation process is governed by the kinetics of the vapor and water phase change in the system. These issues are modeled using a transport equation with the source terms regulating the evaporation and condensation of the phases. The source term m + and sink term m − in eq. (3) represent the condensation and evaporation rates. In this study, three different transport cavitation models are considered.
(i) Model-1: Kubota cavitation model [7]. A popular phenomenological model was originally proposed by Kubota et al. [7] and investigated by numerous researchers. In this model, the growth and collapse of a bubble cluster were given by a modified Rayleigh-Plesset equation, which provided the rate equation controlling vapor generation and condensation. The Rayleigh-Plesset equation is given by where ρ l is liquid density, R B is the bubble diameter, p v is pressure in the bubble (assumed to be the vapor pressure in the local temperature), p is pressure in the liquid surrounding the bubble, and σ is the surface tension coefficient between the liquid and vapor. Ignoring the second order terms and the surface tension, this equation reduces to The representative liquid-vapor evaporation and condensation rates for this category are shown as follows: where C prod and C dest are two empirical coefficients, designed to account for the fact that they may occur at different rates (condensation is usually much slower than vaporization). To obtain an interphase mass transfer rate, further assumptions regarding the bubble concentration and radius are required. The Kubota cavitation model uses the following defaults for the model parameters: (ii) Model-2: Interfacial dynamics cavitation model-IDM [12]. The interfacial dynamics-based cavitation model [12] is adopted for cavitation modeling. The development of this model is based on the analysis of the mass and normal momentum conservation at a liquid-vapor interface within the context of homogeneous equilibrium flow theory. The interfacial dynamics-based cavitation model reads as follows: In the next step, it is considered that the phase changes takes place between the mixture and vapor phases across clear interfaces by simply replacing the liquid phase with the mixture phase. The above equation defines the value of the liquid volume fraction. In the context of TEM, it is necessary to couple the above interfacial condition as a source term into the transport equation of the liquid volume fraction. For this purpose, normalized with a characteristic time scale, the evaporation and condensation are given by The subscripts I, l, v represent the interface, the liquid phase and the vapor phase. For time-dependent computations, the model requires that an interface should be constructed to obtain the interface velocity (U I,n ), as well as the normal velocity of the vapor phase. The vapor phase normal velocity is the product of the velocity and the normal vector: .
The interface velocity is only needed for time-dependent problems, and requires an additional method to track the movement of the interface. Development of such an interface tracking scheme for turbulent two-phase flows in body-fitted curvilinear coordinates deserves a separate study.
In this study, the interface velocity is estimated based on a simplified approach, using the mass conservation condition as follows: (iii) Evaluation of the cavitation models with experi-mental results. Wang et al. [14] performed a series of experiments to study time-dependent cavitating flows around a Clark-Y hydrofoil. Figure 1 shows the shapes of cavities during an oscillation cycle. The unsteady cavitation structure consists of two parts, which are the attached and detached cavity, respectively. The attached cavity is located in the leading edge of the hydrofoil, while the detached cavity is formed because of the re-entrant jet and overlaps with the recirculation zone near the trailing edge. Cloud cavitation has the following qualitative features: the initiation of the cavity, growth towards the trailing edge, and subsequent shedding, which are observed experimentally. The capability for simulating unsteady cavitating flows with the different cavitation models is further investigated. The computational domain and boundary conditions are given according to the experimental setup [14], which is shown in Figure 2. The Clark-Y hydrofoil is located at the center of the test section with the angle of attack equal to 8°. The two important scaleless parameters used are the Reynolds number Re and the cavitation number σ, which are defined as where U ∞ is free stream velocity, ν is kinematic viscosity, and c is the chord length of the hydrofoil.
Computations have been done under cloud cavitation conditions (σ=0.8). The velocity is imposed at the inlet (U ∞ =10 m/s in the present case), and the corresponding Reynolds number is 7×10 5 .
The comparisons of the time evolutions of the attached cavity length predicted by the above mentioned models and experimental data are presented in Figure 3. The cavity   length l is made non-dimensional using the chord length c of the hydrofoil. As shown in the experimental results, one cycle of the attached cavity oscillation can be divided into two stages: the first stage is the growth of the attached cavity (during 0 to 37.5% of the cycle), and the second stage is the development of re-entrant flow in the rear of the cavity (during 37.5% to 100% of the cycle). At the beginning of the second stage, the attached cavity length does not change within about 1.0c. When the adverse pressure gradient is strong enough to overcome the weaker momentum of the flow confined by the near-wall region, the re-entrant flow forms and pushes the flow towards the leading edge during the development of the attached cavity. Finally, the cavity is cut into two parts by the re-entrant flow, and the attached cavity length suddenly drops at 66.5% of the cycle.
In the first stage, a good agreement between the simulated and experimental results is obtained concerning the attached cavity length. However, some differences are observed in the second stage: Model-1 produces cavity breakup and detachment too early, while Model-2 fails to predict the process of the detached cavity, because the fact that the interaction between the turbulence and the vapor collapse in the closure region is not considered in this model.
Model-1 and Model-2 are both transport equation-based models. No single model choice performs better in all respects: Model-1 can explain the interaction between viscous effects including vortices and cavitation bubbles. In this model, the Navier-Stokes equations including cavitation bubble clusters are solved, where the growth and collapse of a bubble cluster are given by a modified Rayleigh's equation. As a result, this model effectively expresses the mechanism of cavitation cloud generation and larger-scale vortices, which are the main features of cloud-type cavitation shedding from the trailing edge of the attached cavities [7]. The development of Model-2 is based on the analysis of the mass and normal momentum conservation at a liquid-vapor interface. The model has been applied to different cavitating flow problems and has produced results which are in good agreement with the experimental data, especially predicting the unsteady interface between vapor and mixture phases [12].

Modified density based cavitation model-DMBM (Model-3)
In the following study, a modified density based cavitation model is developed to capture the features of unsteady cloud cavitation structures: the interface between the liquid and vapor at the front of an attached cavity and detailed process of the vortex shedding of the detached cavity. In this new model, the liquid-vapor evaporation and condensation rates are shown as Here, a blending function χ(ρ m /ρ l ) is used to combine the mass transfer rates in different cavitation areas, which is given as where C 1 and C 2 are chosen to be 4 and 0.2. The hybrid function χ is shown in Figure 4, which will blend Model-1 and Model-2 based on the local mixture density. In the following, we will utilize Model-1, Model-2 and Model-3 to investigate the impact of cavitation models on cavitating flow simulation results.

Analysis and modeling
The instantaneous contours of liquid volume fraction are compared with experimental data side by side in Figure 5. The cavity visualization is placed according to 10%, 37.5%, 66.5%, and 95.0% of each corresponding cycle.
For Model-1 shown in Figure 5(b), before the cavity is detached, the density inside the attached cavity region has already become very low, i.e. α l <0.1. At 95.5% of the cycle, the density inside the detached cavity still contains 80% of the vapor phase. However, the attached periodic cavity seems to be underestimated by this model in the second stage of cloud cavitation. For Model-2 in Figure 5(c), a good agreement between CFD and the experimental results regarding the attached part of the cavity is achieved. The attached cavity is divided into two sub-areas during 37.5% to 66.5% of the cycle: when the vapor phase is located in the front part of the cavity, a two phase mixture area is in the rear region, and a distinct interface between the two zones can be seen. The interface location is highly unsteady, which indicates a reverse motion from the twophase mixture area to the vapor area, which is similar with the experimental views in Figure 1. However, the IDM model cannot capture the detached cavity at 95.0% of the cycle.
Model-3 combines the advantages of both Model-1 and Model-2 in simulating unsteady cloud cavitating flows. As shown in Figure 5(d), the features of every stage in the experiment can be captured, including the global structure of both the attached part and the detached cavity in the trailing edge of the last stage. Figure 6 shows the time evolution of the attached cavity length with the alternative models. The numerical result for  Model-3 is better agreement with the experimental data, and an abrupt drop in cavity lengths due to the sudden shedding is clearly predicted. Figure 7 presents the time evolutions for the water vapor fractions of the cavity from the leading edge to the trailing edge in each case, with different cavitation models. Figure  7(a) shows an example of the numerically predicted flow. A reasonable agreement is obtained between the two results including the cavity growth, the vapor cloud detachment, and its convection downstream. The numerical simulations with the alternative cavitation models are implemented in Figures 7(b) For the alternative cavitation models, different periodic times for the cavity shapes are obtained. Model-2 fails to capture a large scale structure that sheds clouds of vapor, and results in a longer period as expected, while compared The time-averaged drag, lift coefficients and frequencies obtained by numerical simulations are also provided for comparison with experimental data in Table 1. The agreement between CFD and experimental data in terms of mean lift and drag coefficient is good except for Model-2. In particular the mean drag coefficient for Model-2 overpredicts C l and C d .
It can be seen that these are highly dependent on the features of the cavity shapes and the dynamic behaviors. The frequencies obtained from Model-1 and Model-3 are almost the same as the experimental data which are close to 24 Hz. The main factor that explains the close similarity between Model-1 and Model-3 involves the re-entrant jet, which triggers the shedding and unsteady motion, and basically consists of a high liquid volume fraction. Model-1 is more influential than Model-2 in this area. The time-averaged x-direction velocity of the flow field is illustrated in Figure 8. These time-averaged velocity profiles are tracked along the vertical direction at different locations, namely x/c = 40%, 60%, 100%, and 120%. The difference between prediction and measurement is substantial, especially at the cavity closure region. In fact, recirculation breaks up the cavity and tilts the rear portion of the cavity upward, which makes the cavity thicker in comparison to the experiment. Although computations do not match the experimental data quantitatively, the overall trends are in agreement.
As discussed above, Model-1 captures the detached cavity with less time in one cycle, and the re-entrant flow is not substantial because of the fluid forward velocity lasting a relatively long time. For Model-2, the computations fail to predict the detached cavity, and the local velocity which will sweep the detached cavity downstream is noticeably underestimated. All the above content can account for different velocities between Model-1 and Model-2 in the near wall region. Model-3 tends to be similar to Model-2 near the wall region, where cavitation can occur easily, so there is high similarity in the velocity distribution between Model-2 and Model-3 in this area. Reasonable agreement is observed in the velocity profiles for Model-3.

Discussion and conclusion
In this study, we illustrate the issues associated with flow computations for unsteady cloud cavitating flows. First, we presented widely-used cavitation models to highlight the merit and weakness of these models, and also provided a survey of cavitation studies in terms of unsteadiness and flow structure. We then presented our efforts in developing a modified density based cavitation model using the unsteady characteristics of cavitating flow.
The Kubota, IDM and DMBM models are considered as cavitation models for the computation of unsteady cavitating flows around a Clark-Y hydrofoil. The study suggests that the existing versions of the TEM have merits and limitations. For the Kubota model, the attached cavity tends to stay for a shorter time in a period cycle, and results in less substantial re-entrant flow. However, for the IDM model, the cavity cutting and cloud shedding is not obviously captured. For the DMBM model, the features of every stage in the experiment can be captured successfully, including the attached cavity development, the detached cavity in the trailing edge of the last stage and the form of vortex shedding flow which quickly moves downstream. Furthermore, the same dynamic characteristics as the experimental measurements are obtained. Overall, the new modified density based closure model shows considerable improvement in its predictive capability. In summary, this study has applied a concept to blend different cavitation models and uses this idea to simulate the cavitating flows. Compared with other existing models, better agreement with experimental data is obtained. Our current study provides useful information for further modeling modifications.