A Finite Element Model of Unsteady Cavitating Fluid Flow around a Hydrofoil

The present study deals with the unsteady dynamics of cavitation around the NACA 0015 hydrofoil in a channel. A finite element model is proposed to solve the governing equations of momentum and mass conservation. Turbulent flows around the hydrofoil are described by the Prandtl-Kolmogorov model. The cavitation phenomenon is modeled through a mixture model in-volving liquid and vapor flows and the Zwart-Gerber-Belamri (ZGB) model is considered to evaluate the transport of the water vapor fraction. The variational finite element model formulation includes the mixing of the characteristic method and the finite element. Also, at the open sides of the channel flow, an open boundary condition is imposed. Numerical experiments are performed for cavitation numbers 0.8 and 0.4. The presented model predicts the essential features of unsteady cavitating flows, the generation of vapor cavities, the time-dependent oscillations of the variables and the presence of vortical flow structures associated to vapor volume concentrations during the shedding process.


Introduction
The fluid flow around submersed bodies includes many fluid mechanics phenomena. Cavitation is one of these phenomena, which is a sequence of vaporization and condensation processes during the pressure oscillation of the flow forced by high velocities around a body. The cavitation is not well understood and many questions are open.
The fluid circulation around bodies produces complex processes. The result-ing fields of pressure and velocity around the body are modified due to the geometry. When the Reynolds number increases, the fluid flow begins to separate with the generation of unsteady vortex motions mainly behind the body [1]. The turbulent behavior of the fluid is an open question and exists different options to model the turbulence. The more used turbulent formulations are the family of k-ε and k-w models. A problem of these kinds of models is that they are dependents on the geometry of the case considered and also suffer from the deficiencies of the gradient ansatz [2]. An additional problem of these models is the presence of many constants, which are not universal constants. Also, it is well known the overprediction of the turbulence viscosity [3] that numerically dampens the unsteadiness of the cavity bubbles. In spite of this, it is frequently used in many software packages. Another option is the use of zero or one equation turbulent models (e.g. Smagorinsky model, Prandtl-Kolmogorov model). The use of different turbulent models leads to discrepancies in the pattern results.
Also, the deficiency of the standard models was also reported by different authors [4]. The eddy viscosity depends on the non-uniform characteristics of the flow velocity field and the Prandtl-Kolmogorov turbulence model describes this concept in a consistent way. Here, this model is adopted.
The mixture model of water and vapor uses a transport equation to describe the rate of change of the water vapor fraction. The use of transport equation models has an advantage because can predict the dynamic influence of momentum on vapor cavities deformation and drift of bubbles. This kind of model applies different condensation and evaporation empirical coefficients to regulate the mass exchange of water and vapor, which is the case among others of the Singhal model [5] and the Zwart-Gerber-Belamri model [6].
The numerical modeling of Cavitation is a challenge because numerical uncertainties of the models produce important changes in the solutions [7]. The unsteady generation and collapse of the vapor cavities induce an oscillatory behavior, which is approximately periodic in time as reported in the literature [8].
Here is necessary to remark, that exist discrepancies between the numerically calculated oscillations frequencies reported by different authors during the cavitation phenomena [8]- [13]. Also, laboratory experiments, have reported that oscillatory frequencies change in function of the cavity length [14]. Different frequency oscillations during cavitation were also observed [10] [13] [15].
The difficulties of modeling cavitation flows are mostly associated with the non-permanent dynamics of the problem, the spatial precision of vapor bubbles and the resulting velocity field (vortices, reentrant flows). Some articles, comparing the performance of cavitation models [16] [17] and capturing cavity morphology [18], provide a recent overview of the numerical modeling of cavitation flows.
In the literature, many numerical solutions were reported about flow motions over hydrofoils. Most of them use finite differences and finite volume techniques [11] [19] [20]. The applications of finite element solutions are mostly unex- This paper is a step in the study of such problems using a developed model based on a Characteristic Galerkin finite element formulation [25]. This FEM option is very efficient and easy to implement while the accuracy of the results obtained by the algorithm is still ensured [26] [27].
Here a finite element model is presented. It is capable of describing the hydrodynamic behavior of a flow around a hydrofoil NACA0015 and the resulting cavitation process from its initial state. The Prandtl-Kolmogorov turbulence model is utilized and the numerical performance in the two-dimensional cavitation flow of the Zwart-Gerber-Belamri model is evaluated. The model predicts the generation of vapor cavities and vortex motion structures during the cavitation.

The Hydrodynamic Model
In the study domain Ω, the two-dimensional hydrodynamic incompressible flow around a hydrofoil is described by the momentum and continuity equations of a turbulent fluid of density ρ in a vertical cartesian coordinate system ( ) In the present paper, the eddy viscosity is modeled following the Prandtl-Kolmogorov turbulent model written as  where c is a constant usually equal to 0.54, k is the turbulent kinetic energy and l represents the characteristic mixing length. The kinetic energy k is calculated according to: and ε is the dissipation of kinetic energy approached as 3 2 . k c l ε ε ρ = (9) and c ε is a constant. The flow dynamic is forced by a specified input boundary flow at the entrance of the channel, upstream from the leading edge of the hydrofoil. Along the entrance (inflow side), a Dirichlet boundary condition u u ∞ = is imposed.
Mixture flows could be studied considering a simple single-fluid approach and a transport equation for the vapor mass fraction f could be written as where ρ is the mixture density. The relation between the density mixture ρ and the vapor mass fraction f is described by where l ρ dis the density of liquid and v ρ is the density of the vapor. And the volume fraction of vapor phase v α is described according to: .
In the transport equation model, source terms are based on the Rayleigh-Plesset equation for bubble dynamics and the Zwart-Gerber-Belamri model [6] is applied. The source terms of the transport model are ( ) where ,

The Numerical Model
For the numerical solution, the two-dimensional spacial domain Ω is partitioned in 1 e N triangular subdomains e Ω , with nod N total nodes. For the time domain, an ordered partition of time levels is defined Δt. For a generic variable U(t), a linear approach between the two time levels n and n+1, is expressed as ( ) ( ) Also, the total time derivative is approached by including a characteristic estimation for ˆn U [25]. Here, the parameter θ was fixed to equal 1. In this way, the governing equations read where ˆ, ,

Experiments
The experiments conducted in the present study are based on the flow behavior over a NACA0015 symmetric hydrofoil in a water tunnel. The hydrofoil has a chord length c = 0.1 m and in the present section is studied numerically. A tunnel of length 10c and height 4c is considered. Figure 1 shows The cavitation number σ, which describes the state conditions related to the saturation pressure p v , is defined as   Figure 5 shows the instantaneous field, at two different time instants for the pressure coefficient C p . The dynamical flow behavior around the hydrofoil is shown in Figure 6. A main clockwise vortex motion is produced capturing the vapor which is separated from the sheet cavity tail. Therefore, the cloud cavity, which grows in time is

Solutions When σ = 0.8
The numerical solution shows a cavity of vapor volume fraction generated near the leading edge with a limited length. The time-dependent response of the cavitation phenomena is presented in Figure 7 for the α v and C p values, at two control points located on the upper surface of the hydrofoil. These results show      Figure 8 and Figure 9 show the instantaneous field, at two different time instants, of vapor volume fraction α v and pressure coefficient C p . In the present case, only half of the hydrofoil upper surface is covered with vapor. Additionally, the corresponding velocity fields are presented in Figure 10.     . Calculated relative cavity length compared with experimental data (2D, 3D) extracted from Arndt [15].

Summary and Conclusion
The

Conflicts of Interest
The author declares no conflicts of interest regarding the publication of this paper.