Numerical Investigation for Steady and Unsteady Cavitating Flows

Cavitation is a particular two-phase flow with phase transition (vaporization/condensation) driven by pressure change without any heating. It can be interpreted as the rupture of the liquid continuum due to excessive stresses. Modeling of cavitating flows as a multi-fluid is a complex problem especially when the 3D consideration is adopted. However most recently research works are based on the mixture consideration of homogeneous fluid composed by two phases of liquid and vapor, which is also used by our model....


Introduction
Cavitation is a particular two-phase flow with phase transition (vaporization/condensation) driven by pressure change without any heating. It can be interpreted as the rupture of the liquid continuum due to excessive stresses. Modeling of cavitating flows as a multi-fluid is a complex problem especially when the 3D consideration is adopted. However most recently research works are based on the mixture consideration of homogeneous fluid composed by two phases of liquid and vapor, which is also used by our model….
Minimizing the nuisance of cavitation is a great challenge in the design phase of a marine propeller. For efficiency reasons, the propeller usually needs to be operated in cavitating conditions but one still needs to avoid the effects of vibrations, noise and erosion. However, cavitation is a complex phenomenon not yet neither reliably assessable nor fully understood. Experimental observations can only give a part of the answer due to the obvious limitations in the measurement techniques; one example is measuring reentrant jets and internal flow, where flow features are hidden for optical measurement techniques by the cavity itself. Standard simulation tools used in design typically include potential flow solvers, lifting surface or boundary element approaches, with strict theoretical limits on cavitation modeling that only in the hands of an experienced designer may give satisfactory propeller designs. Adding to the challenge is a lack of theoretical knowledge of the physical mechanisms leading to harmful cavitation and thus how to modify a design if some form of nuisance is detected.
The transport equation models of cavitation suggested by Alajbegovic, Grogger et Philipp [ 1] and Yuan, et al. [2], use the simplistic Rayleigh model. This model describes the limiting case of inertia-controlled growth of a spherical bubble in a liquid under a step variation in pressure of the surrounding liquid. However, this model cannot accurately describe bubble collapse and neglects a number of effects, which determine the behavior of cavitation bubbles.
Recently, progress has been made in the development of numerical models for calculation of cavitation flows. Though the models may differ in terms of realization (using the singlefluid or multi-fluid frame-work, the Eulerian-Eulerian or Eulerian-Lagrangian approaches), all of them are empirical to a certain level.
Modeling of cavitation flow as a multi-fluid is a complex problem which does not lead to satisfying results especially when the cavitation is modeled as 3D. Therefore in the engineering practice cavitation flow is often modeled as a single-fluid, where the cavitation area is handled as an area with the pressure lower then the vapour pressure. This approach always leads to the result, and the requirement of computer time is many times lower in comparison with multi-phase flow models. Moreover the steady solution of multiphase flow model may not be found at all due to the unsteady nature of cavitation flow.
Significant progress has been achieved recently in the development of homogeneousmixture (single-fluid ) models for the simulation of three-dimensional transient cavitating flows ( [8]). These models allow single-fluid solvers to be applied to the conservation equations for the mixture, without increase in computational cost due to the increase in the number of conservation equations when applying the multi-fluid flow concept.
The present work is an investigation to develop a relevant physical model to simulate the cavitating flow. The main goal is the development of computational methodologies which can provide detailed description of the numerical set up for modeling and simulation the cavitation with the CFD code.

Mathematical formulation
To simulate cavitating flows, the two phases, liquid and vapour, need to be represented in the problem, as well as the phase transition mechanism between the two. Here, we consider a one fluid, single-fluid (mixture), introduced through the local vapour volume fraction and having the spatial and temporal variation of the vapour fraction described by a transport equation including source terms for the mass transfer rate between the phases. The numerical model solves the Reynolds averaged Navier-Stockes equations, coupled with a localized vapour transport model for predicting cavitation.
The fundamental equations governing the flow are taken with incompressible fluid case as given by the Navier-Stockes equations, which control the transport of momentum within the fluid (2) in addition to the mass conservation constraint (1): The effective density and viscosity of the mixture are given (3) and (4) respectively : Where α is the vapor fraction (α 1: vapor and for α 0: liquid). The density and the viscosity of liquid and vapor are assumed to be constant. To compute the volume fraction we need a closure model, the distribution of values α was obtained on each cell of the computational domain will guide the attendance rate of the vapour.

Vapour fraction transport equation
The transport equation for the vapour scalar fraction α is given by: This transport equation of volume fraction of vapour, with appropriate source terms to regulate the mass transfer between phases (liquid/vapor), is solved.

The proposed model formulation
a. This work deals with a numerical simulation of cavitation process around a hydrofoil. The numerical approach is based on predicting of the collapse process of pocket of vapor due to liquid compressibility. The model based is the Rayleigh Plesset (dynamics of spherical bubbles) integrated in the term source. The interface velocity investigated during the collapse process is given by: finally, the proposed model can be express : with : The value of 0 n is a calibrate parameter with experimental results; it is define the bubble density.

Numerical study
To validate the proposed model, a confrontation to experimental measures and to the numerical models (Yuan et al., Schnerr and Sauer and EOS). The application is a NACA0009 hydrofoil, truncated at 90% of the original chord length. It has the final dimensions of 100 mm of chord length. The 3D test section is modeled by a quasi 2D domain, with three rows of cells in spanwise direction for the numerical domain. The same mesh and numerical setup is used for the computations with the two models.

Geometry
The domain is 9 blocks (See Fig. 1). The boundary conditions are set using a velocity inlet and a pressure at the outlet (the parameter which fixes the cavitation number). The equations of the non-truncated hydrofoil are:

Mesh quality
Numerical solutions of fluid flow and heat transfer problems are only approximate solutions. In addition to the errors that might be introduced in the course of the development of the solution algorithm, in programming or setting up the boundary conditions… The discretization approximations introduce errors which decrease as the grid is refined, and that the order of the approximation is a measure of accuracy. However, on a given grid, methods of the same order may produce solution errors which differ by as much as an order of magnitude. Theoretically, the errors in the solution related to the grid must disappear for an increasingly fine mesh [ HYPERLINK \l "Fer96" 9 ]. The pressure coefficient at =0.4 was taken as the parameter to evaluate six grids (See Fig. 2) and determine the influence of the mesh size on the solution. The selected convergence criteria were a maximum residual of 10 −4 . According to this figure, the grid with 53419 cells is considered to be sufficiently reliable to ensure mesh independence. In these calculations turbulence effects were considered using turbulence models, as the k-ε RNG models, with the modification of the turbulent viscosity. To model the flow close to the wall, standard wall-function approach was used. For this model, the used numerical scheme of the flow equations was the segregated implicit solver.

Turbulence model discussion
For steady state flows, there is a priori no difference between the two equations formulations. The use of SST in the case of 2D steady hydrofoil computations instead of k-ω or k- is justified in this way.
SST model works by solving a turbulence/frequency-based model (k-ω) at the wall and k-ε in the bulk flow. A blending function ensures a smooth transition between the two models. The SST model performance has been studied in a large number of cases. In a NASA Technical Memorandum, SST was rated the most accurate model for aerodynamic applications [11].

Pressure coefficient validation
A validation of the proposed model is described in this section. Firstly, we proceed with a distribution of the pressure coefficient around the surface of the hydrofoil and then we adjust with the profile velocity.   These figures show satisfactory results of the proposed models in predicting the pressure distribution on the hydrofoil. As expected, the cavity becomes larger with decreasing cavitation number. However the models exhibit very different flow behavior at the cavitation detachment and closure regions. This is the result of dynamic bubble which is implementing in the proposed model. Clearly the model can perfectly predict the cavity of the vapor.  Figure 6 show the predicted velocity profile of the main flow far from the wall is in good agreement with measurements. In the near-wall and in the wake region, the proposed model tends to reproduce the flow re-entrant jet.

2D configuration(NACA0009)
The vapor cavity is characterized by a thick main cavity and an important shed cavity volume in a disorganized way, such as the cavity can be divided into many small cavities. For one period of cavity creation and collapse, one can remark two distinct life cycles highlighted by the lift and drag signals. Firstly, starting from the maximum cavity length, the closure region is in the small adverse pressure gradient and the reentrant jet is too thin, such as the cavity closure region is continuously broken into small vapor volumes. Secondly, as the reentrant jet reaches the region of the high pressure gradient, the reentrant jet is more important and the whole cavity is extracted and shed downstream. We consider the non-truncated variant of the 2D NACA0009 hydrofoil at high incidence angles. The flow parameters are i=5° and =1.2, and reflected by a high pressure gradient over the hydrofoil leading to an unsteady state flow behavior and shedding of large transient cavities. The domain is taken larger than the experimental test section to avoid numerical problem mainly due to reflections on the boundaries.
Cavitating flow are highly sensitive to turbulent fluctuations present in the flow [14] [15]. The effects of modeling turbulence quantities have an enormous impact on the cavitation dynamics and the overall flow structure. The k- RNG turbulence model and LES model are adopted in this study.

Result and discussion
In this section we present the behavior of the proposed model with LES and the modified turbulence model. Most of characteristics of the unsteady flow is defined as a function of lift and drag coefficients, this coefficient can be expressed as: Numerical Investigation for Steady and Unsteady Cavitating Flows 89 For flows about 2D hydrofoils the dimensionless total vapour volume Vvap,2D is defined as : (10) where N is the total number of control volumes and α is the volume fraction of vapour in each control volume, with the volume Vi of the fluid in control volume.
The total vapour volume Vvap,2D, defined in equation (10), is a convenient parameter for understanding the transient evolution of the cavitating flow. The total vapour volume is calculated at each time step. After the start-up phase the growth and shedding of the vapour sheet and the collapse of the shed vapour cloud induce a self-oscillatory behavior, which is approximately periodic in time. The graphics of the analyzed variables (Fig.8) show a periodic signal, even if the identified periods are very different from each others. The use of Vvap,2D is the easiest way to identify the periodicity of the vapour formation and collapse during a life cycle. Furthermore, the time-history of the lift and drag coefficients are compared to those of the total vapour volume to correlate the occurring flow phenomena.The lift and drag coefficient has a cyclic time signal. Even if the periods are pretty clear, the fluctuations in a given period are very different from one to another. This phenomenon is mainly due to the dynamic of the shed cavities which are driven by the main flow-field downstream of the cavity closure. The growth and collapse mechanism of the cavity is driven by a cyclic phenomenon, whereas the dynamic of the cavity, when it is swept away can have very different non reproducible behavior. The vapour can be attached to the wall or far from it. This different ways of cavities shedding have an important influence on the pressure field at the hydrofoil wall, and there by on the lift and drag values.
In this section the cycles illustrated in figures 15 and 16 are considered. The solution for the volume fraction α above the hydrofoil is presented for a number of equidistant timeintervals during the cycle.
The vapour cavity is characterized by a thick main cavity and an important shed cavity volume in a disorganized way, such as the cavity can be divided into many small cavities. For one period of cavity creation and collapse, one can remark two distinct life cycles highlighted by the lift and drag signals.
Firstly, starting from the maximum cavity length, the closure region is in the small adverse pressure gradient and the reentrant jet is too thin, such as the cavity closure region is continuously broken into small vapour volumes. Secondly, as the reentrant jet reaches the region of the high pressure gradient, the reentrant jet is more important and the whole cavity is extracted and shed downstream.

3D configuration(NACA66mod)
These excellent results obtained for 2-D cavitating flow confirm the correct assumptions of the proposed numerical model in terms of vaporization and condensation processes, and to verify it's performance with 3-D considerations, we present in the next some numerical results compared to experimental measurement [16] for cavitating flow around hydrofoil NACA66 (Fig 10). In this section, a numerical simulation was performed to further study the performance of the proposed model and unsteady three-dimensional, the choice is justified because of the availability of experimental measurements. The profile used is NACA66 (mod) -312 a = 0.8 with a string of 0.15m, the width of the profile is 0.19m and is placed at 6 ° incidence relative to the upstream flow, not in an infinite medium -viscous. It has a thickness of approximately 12% and a camber on 2% to 45% and 50% of the leading edge along the rope. Theoretical points of this section have been interpolated by B-spline technique using the software mesh.
The study area is in 3D, it is made larger than the experimental test section to avoid numerical problems mainly due to reflections from the boundary conditions. The estate consists of 156 000 cells, the mesh topology is structured with a C-type conditions with wall on its border, the hydrofoil with the no-slip condition. (see Fig. 8). The steady state solution was used as an initial condition, the turbulence model used is k- RNG.

No slip Wall
Outlet Inlet

No slip Wall
Hydrofoil : no slip Wall Figure 11. Domain grid and boundary condition at the left, at the right, the validation of the pressure measurement [16] respectively at 50%,70% and 90% of the chord length with the proposed model, Cref= 5.3 m/s. The proposed model was validated quantitatively and qualitatively compared in twodimensional distribution of pressure coefficient and velocity profiles calculated with experimental measurements. We propose in the following one-dimensional confrontation, the application as we have described above is the NACA66 (mod), the advantage of selecting this application is the availability of measuring pressure along the profile for two cycles of detachment from the pocket of steam. We have reported three points of pressure measurement located respectively 50%, 70% and 90% chord.
The obtained numerical result presented by figure9 show that the model reproduce correctly the typical behaviour of partial cavity with development of re-entrant jet and the periodic shedding of cavitation clouds. Firstly, starting from the maximum cavity length, the closure region is in the small adverse pressure gradient and the re-entrant jet is too thin, such as the cavity closure region is continuously broken in to small vapour volumes. Secondly, as the re-entrant jet reaches the region of the high pressure gradient, the re-entrant jet is more important and the whole cavity is extracted and shed downstream. This different ways of cavities shedding have an important influence on the pressure field at the hydrofoil wall, and there are characterized by the temporal evolution of the lift coefficient Cl compared to experimental result [16]

Conclusion
A comprehensive theoretical approach is done, and detailed formulations of the proposed model are presented. Influence of the numerical parameter has been widely studied. A comparative study is made for the steady flow. Good agreement with measurements was obtained for the proposed model. We have shown the importance of liquid compressibility on the pocket of vapour and how it is modeling in the proposed model. Computations with the proposed model is compared with experimental data and other numerical models, it shows the ability of the model to reproduce the steady-state developed cavitation flow fields. Finally, the unsteady behaviour of the cavitating flow depends strongly on the turbulence model, a k- RNG model is adopted in this steady.

Author details
Hatem Kanfoudi, Hedi Lamloumi and Ridha Zgolli Laboratory of Hydraulic and Environmental Modelling, National Engineering School of Tunis, Tunis, Tunisia