Abstract

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

Keywords

Cavitation, Finite Elements, Hydrofoil, Unsteady Flows, Zwart-Gerber-Belamri Model

Share and Cite:

1. 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 resulting 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 unexplored. The finite element method is a powerful numerical technique to solve engineering problems [21] [22] [23]. Complex geometries and irregular shapes are easier to describe. Comparisons with another method [24] have verified the better behavior of finite element method (FEM) options compared to the finite volume method (FVM).

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.

2. 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 $\left(x_1,x_2\right)$ with turbulent velocities $u=\left(u_1,u_2\right)$, pressure p, $\bar{g}=\left(0,g\right)$ and eddy viscosity ${\nu }_T$ :

$\rho \left(\frac{\partial u}{\partial t}+u\nabla u\right)+\nabla p-\rho {\nu }_T\nabla \nabla u+\rho \bar{g}=0,$ (1)

$\nabla u=0,$ (2)

where $\nabla =\left(\partial /\partial x_1\partial /\partial x_2\right)$ is the Nabla operator. The study domain is composed of solid and open boundaries. Non-slip boundary conditions are prescribed on the surface of hydrofoil,

$u=0,$ (3)

free-slip boundary conditions for the velocity and natural boundary conditions for the pressure are defined on the channel boundary walls

$u_n=0,$ (4)

$\frac{\partial p}{\partial x_n}=0$ (5)

where subindice n indicates the normal direction to the wall. Also, weakly reflective dynamic conditions must be satisfied on the open side boundaries of the domain, this condition could be written for the present case as

$\rho \frac{\partial u_n}{\partial t}+{\nabla }_{n}p=0$ (6)

In the present paper, the eddy viscosity is modeled following the Prandtl-Kolmogorov turbulent model written as

${\nu }_t=c\rho l\sqrt{k},$ (7)

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:

$\begin{array}{c}\frac{\partial k}{\partial t}+u\nabla k-\nabla {\nu }_T\nabla k-R+\epsilon =0\end{array},$ (8)

where $R=\frac{{\nu }_T}{2}{\left|\nabla u+\nabla u^{\text{T}}\right|}^2$ and ε is the dissipation of kinetic energy approached as

$\epsilon =\rho c_{\epsilon }\frac{k^{3/2}}{l}.$ (9)

and $c_{\epsilon }$ 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_{\infty }$ 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

$\frac{\partial \rho f}{\partial t}+\nabla f\rho \mu -\left(S_e-S_c\right)=0$ (10)

where ρ is the mixture density. The relation between the density mixture ρ and the vapor mass fraction f is described by

$\frac{1}{\rho }=\frac{f}{{\rho }_v}+\frac{1-f}{{\rho }_l},$ (11)

where ${\rho }_l$ dis the density of liquid and ${\rho }_v$ is the density of the vapor. And the volume fraction of vapor phase ${\alpha }_v$ is described according to:

${\alpha }_v=f\frac{\rho }{{\rho }_v}.$ (12)

The presented Equations (1), (2), (8), (10) and (11) are considered to calculate the unsteady response at each time instant the variables $u=\left(u_1,u_2\right)$, pressure p, kinetic energy k, vapor mass fraction f and the density mixture ρ.

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

$S_e=C_e\frac{3{\alpha }_{nuc}\left(1-{\alpha }_v\right){\rho }_v}{R_B}\sqrt{\frac{2}{3}\frac{P_v-P}{{\rho }_l}},P<P_v$ (13)

$S_c=C_c\frac{3{\alpha }_v{\rho }_v}{R_B}\sqrt{\frac{2}{3}\frac{P-P_v}{{\rho }_l}},P>P_v$ (14)

where $S_e,S_c$ represent source terms for vapor generation and vapor condensation respectively. The $R_B$ is the Bubble radius and ${\alpha }_{nuc}$ is the nucleation site volume fraction. The constants $C_e,C_c$ are evaporation and condensation coefficients respectively. Assuming that all the bubbles in a system have the same size, the Zwart-Gerber-Belamri cavitation model proposed that the total interface mass transfer rate per unit volume is calculated using the bubble density numbers.

3. The Numerical Model

For the numerical solution, the two-dimensional spacial domain Ω is partitioned in $N_{e1}$ triangular subdomains ${\Omega }_e$, with $N_{nod}$ total nodes. For the time domain, an ordered partition of time levels is defined

$0=t^0<t^1<\cdots <t^n<t^{n+1}<\cdots <t^M=T,$ (15)

where T is the end time instant. A time interval is denoted by $\left[t^n,t^{n+1}\right]$ of length Δt. For a generic variable U(t), a linear approach between the two time levels n and n+1, is expressed as $U\left(t\right)=\theta U^{n+1}+\left(1-\theta \right)U^n$, where $\theta =\left(t-t^n\right)/\left(t^{n+1}-t^n\right)$. Also, the total time derivative is approached by $\text{d}U\left(t\right)/\text{d}t=\left(U^{n+1}-{\hat{U}}^n\right)/\Delta t$ including a characteristic estimation for ${\hat{U}}^n$ [25]. Here, the parameterθ was fixed to equal 1. In this way, the governing equations read

$L_M=\frac{\rho u^{n+1}-\rho {\hat{u}}^n}{\Delta t}+\nabla p^{n+1}-\rho v_T\nabla \nabla u^{n+1}+\rho \bar{g}=0$ (16)

$L_C=\nabla u^{n+1}=0$ (17)

$L_K=\frac{k^{n+1}-{\hat{k}}^n}{\Delta t}-\nabla v_T\nabla k^{n+1}-R+\epsilon =0$ (18)

$L_F=\frac{\rho f^{n+1}-\rho {\hat{f}}^n}{\Delta t}-S_e+S_c=0$ (19)

where ${\hat{u}}^n,{\hat{k}}^n,{\hat{f}}^n$ are defined using a characteristic approach. For a generic variable U(t), the term ${\hat{U}}_j^n$, is obtained in the function of the vector field u and the particle path $s\left(x_1,x_2\right)$, such that ${\hat{U}}_j^n=U\left(s-u\left(s\right)\Delta t\right)$. The mixing method of characteristics and finite element method [25] is a good way which gives satisfactory solutions.

To develop a variational formulation, it is necessary to define discrete functional spaces $V^h=\left\{v\in V^h:{v|}_{{\Omega }_e}={\left(P_b{}^1\left({\Omega }_e\right)\right)}^2\right\}$ and $Q^h=\left\{q\in Q^h:{q|}_{{\Omega }_e}=P^1\left({\Omega }_e\right)\right\}$, where $P_b^1\left({\Omega }_e\right)$ is a linear piecewise continuous finite element plus bubble and $P^1\left({\Omega }_e\right)$ a linear piecewise continuous finite element.

The variational weak formulation of the unsteady hydrodynamic boundary value problem reads: Find u in a functional space V^h and $\left[p,k,f\right]$ in space Q^h, such that for $\forall w^u\in V^h$ and $\forall \left(w^p,w^k,w^f\right)\in Q^h$, satisfy:

${\displaystyle {\int }_{\Omega }{w}^u{L}_M\text{d}\Omega }=0,$ (20)

${\displaystyle {\int }_{\Omega }{w}^p{L}_C\text{d}\Omega }=0,$ (21)

${\displaystyle {\int }_{\Omega }{w}^k{L}_K\text{d}\Omega }=0,$ (22)

${\displaystyle {\int }_{\Omega }{w}^f{L}_F\text{d}\Omega }=0,$ (23)

and also the boundary constrains (3), (4), (5) and (6). The resulting system of linear equations at each time level is solved by the direct method of decomposition L-U. The FreeFem++ Language [23] was used to implement the presented finite element model. In the present paper, the solutions are obtained strictly from the numerical approach of the governing equations and boundary conditions.

4. 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 description of the domain. The NACA 0015 hydrofoil is located in the middle of the water tunnel between 0.45 m ≤ x ≤ 0.55 m and two control points located on the upper surface of the hydrofoil are considered. One of them, point A is at x = 0.46 m, near the leading edge point. The other one, point B, is at x = 0.53 m, near the trailing edge point. In the present paper, the solutions are obtained strictly from the numerical approach of the governing equations and boundary conditions.

The two-dimensional study domain is partitioned into linear triangular elements. The obtained finite element mesh is composed of 14,120 triangular elements and 7220 nodes (Figure 2). Slip boundary conditions are imposed in the upper and lower tunnel walls. Non-slip conditions are imposed on the surface of the hydrofoil. Additionally, the considered reference pressure p_∞ increases hydrostatically with depth.

At the open boundaries, an adequate condition is imposed. In this way, the

contamination due to wave reflections generated during the experiment is controlled. Also, at the entrance boundary (left open side of the tunnel) an inflow uniform velocity u_∞ = 6 m/s is considered, whereas at the output boundary (right open side) the hydrostatic pressure is imposed.

At instant t = 0, a hydrostatic pressure increasing from the upper to the lower tunnel wall is imposed in the tunnel. An initial velocity is defined at all interior points. Therefore, during the initial calculation stage, the dynamics show strong changes.

The following parameters are used in the experiments: the density of water at 25˚C is fixed as ρ_l = 997 kg/m³, the water dynamic viscosity is equal to μ_l = 8.91 × 10⁻⁴ Pa s, the vapor water density is ρ_v = 0.02308 kg/m³, the vapor dynamic viscosity is μ_v = 9.8626 × 10⁻⁶ Pa s and the vaporization pressure is equal to P_v = 3169 Pa. The mixing length is fixed equal to l = c/200 with a chord length of c = 0.1 m. Computations are integrated forward in time step by step up to 0.5 s. The time step is Δt = 0.0005 s.

In the literature, the following model parameters [6] work well for a variety of fluids and devices, the bubble radius R_B = 10⁻⁶ m, the nucleation site volume fraction α_nuc = 5 × 10⁻⁴, the evaporation coefficient C_e = 50 and condensation coefficient C_c = 0.01. But, these parameters are nonuniversal constants and need calibration. That is the case of a reported experiment [6] where values of C_e = 0.4 and C_c = 0.001 were used, indicating changes of two orders of magnitude. In the experiments developed here, the following values were used: R_B = 10⁻⁶ m, α_nuc = 5 × 10⁻⁴, C_e = 41 and C_c = 0.0000081.

For the evaluation of the C_p field, it is considered that the pressure coefficient at the stagnation point is maximal and equal to 1. The pressure coefficient is defined as

$C_p=\frac{p-p_{\infty }}{\frac{1}{2}\rho u_{\infty }^2}$ (24)

The cavitation number σ, which describes the state conditions related to the saturation pressure p_v, is defined as

$\sigma =\frac{p-p_v}{\frac{1}{2}\rho u_{\infty }^2}$ (25)

The numerical solutions are performed for two cavitation numbers, σ = 0.4 and σ = 0.8. The evaluated solutions presented in the present section show the results after t = 0.1s of calculation. During the initial phase (from t = 0 s to 0.1 s), the response is dominated by strong changes of the variables adjusting the initial state. The model solutions show time-dependent oscillations of the pressure and vapor volume.

4.1. Solutions When σ = 0.4

For α_v and C_p at two control points, the unsteady response of the cavitation phenomena using the ZGB model can be appreciated in Figure 3. The solutions show the cavity of vapor becomes highly unsteady, oscillating cyclically. Continuous cycles of an increase and decrease of the pressure coefficient and volume fraction formation during the cavitation phenomena were observed. The intensity of the vapor fraction is high at the leading edge reaching values near saturation.

The oscillations of the cavity conditions are the most particular features of the simulation experiments. It is observed oscillations for α_v and C_p with a period of 0.021335 s (46.87 Hz). Also, a period oscillation of 0.1280 s (7.81 Hz) is observed

In this experiment, the cavitating flow develops a sheet cavity near the leading edge. The sheet cavity has a long tail that is lifted from the foil. Separated cavities are produced at the end of the long sheet cavity and the other cavity is generated near the trailing edge.

Figure 4 shows the evolution of the vapor cavities. 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 formed and is shed downstream. The clockwise vortical flow structure induces a returning flow (re-entrant flow) and forces a secondary counterclockwise vortex (Figure 6) very close to the trailing edge of the hydrofoil. This vortex also captures vapor. This vapor cloud grows (Figure 4(d) to Figure 4(h)) and then it is shed downstream.

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

oscillations of pressure coefficient and vapor fraction of 0.016 s (62.5 Hz) near the leading edge point. At the control point near the trailing edge of the hydrofoil, there is not a significant presence of α_v. The C_p time evolution at the control points shows values around −0.53 at point A, near the leading edge point, whereas at point B the C_p values are around −0.15.

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.

The frequencies calculated in the present work are a particular feature of the unsteady solutions. The literature had reported numerical results predicting many frequencies, for example, 24 Hz and 40.9 Hz [8], 9 Hz [9], 32.1 Hz and 65 Hz [10], 7.75 Hz [11], 17 Hz [12], 120 Hz and 250 Hz [13]. Otherwise, different experimental frequency oscillations were observed, for example, 18 Hz [15] and 120 Hz, 285 Hz [13]. Experiments, indicating that oscillatory frequencies change

according to the cavity length [14], were also reported. A comparison between the calculations presented and existing publications of experimental and numerical results were performed. The comparison (Figure 11) between the relative cavity length (l/c) of numerically simulated flow and experimental results is

plotted versus σ/2α (for σ equal to 0.4 and 0.8) and compared to experimental data extracted from Arndt [15]. The parameter α is the angle of attack in radians. For the ZGB model, when σ = 0.4, the relative length l/c = 1.35 and σ/2α = 1.90985977 and when σ = 0.8, l/c = 0.75 and σ/2α = 3.81971955.

5. Summary and Conclusion

The flow dynamics around a hydrofoil are investigated using a proposed finite element model in a domain with open boundaries. The model includes the Prandtl-Kolmogorov turbulence formulation and the Zwart-Gerber-Belamri (ZGB) model to evaluate the transport of vapor water. Also, the model includes a characteristic scheme to approach the advection terms and at the open sides of the channel flow, open boundary conditions are imposed. The experiments show the generation of vapor bubbles in the upper side of the hydrofoil by the decrease in pressure. The bubble formation has an unsteady behavior. Solutions obtained for the cavitation number σ = 0.4, predict a sheet vapor cavity with a very long tail, the generation of a cloud cavity associated to a main clockwise vortex which induces a counterclockwise secondary vortex on the trailing edge. The intensity of the vapor volume generation is stronger when σ = 0.4 compared to the solutions with σ = 0.8. The solutions are non-permanent with low pressures located on the upper side. Unsteady behavior of the cavitating flows shows the main frequencies of 54.7 Hz for σ = 0.4 and 62 Hz for σ = 0.8.

Acknowledgements

This work is supported by founds of the Universidad Nacional Mayor de San Marcos (UNMSM) of Perú within the activities of the Research Group of Numerical Modeling in Fluid Mechanics.

Conflicts of Interest

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

References