Prediction of the Cavitation over a Twisted Hydrofoil Considering the Nuclei Fraction Sensitivity at 4000 m Altitude Level

: Cavitation phenomenon is important in hydraulic turbomachineries. With the construction of pumping stations and hydro power stations on plateau, the inﬂuence of nuclei fraction on cavitation becomes important. As a simpliﬁed model, a twisted hydrofoil was used in this study to understand the cavitation behaviors on pump impeller blade and turbine runner blade at different altitude levels. The altitudes of 0 m, 1000 m, 2000 m, 3000 m and 4000 m were comparatively studied for simulating the plateau situation. Results show that the cavitation volume proportion f cav increases with the decreasing of cavitation coefﬁcient C σ . At a speciﬁc C σ , high altitude and few nuclei will cause smaller size of cavitation. The smaller C σ is, the higher the sensitivity ∆ f cav is. The larger C σ is, the higher the relative sensitivity ∆ f cav* is. On the twisted foil, ﬂow incidence angle increases from the sidewall to mid-span with the decreasing of the local minimum pressure. When C σ is continually decreasing, the size of cavitation extends in spanwise, streamwise and thickness directions. The cavity is broken by the backward-jet ﬂow when C σ becomes small. A tail generates and the cavity becomes relatively unstable. This study will provide reference for evaluating the cavitation status of the water pumps and hydroturbines installed on a plateau with high altitude level.


Introduction
Cavitation is a commonly seen phenomenon in hydrodynamic turbomachinery cases such as pumps and hydroturbines [1][2][3]. It usually has bad influences including noise [4], vibration [5] and material damage [6]. For this reason, cavitating flow in hydraulic turbomachinery has been widely studied by researchers to understand its mechanism, behavior and influence. As a simplification of the blade of hydraulic turbomachinery, hydrofoil is usually used in the study of cavitation and cavitating flow. Arabnejad et al. [7] and Tao et al. [8] investigated the leading-edge cavitation over symmetric hydrofoils. Different types of leading-edge cavitation are discussed in detail. Dreyer et al. [9] and Guo et al. [10] used hydrofoil to have a better study of the tip-leakage cavitation of blade. Cavitation happening in the tip-leakage vortex was studied to help the anti-cavitation design of axial-flow hydraulic runners. Escaler et al. [11] and Melissaris et al. [12] used hydrofoil models to study the cavitation erosion behavior of hydraulic runner blades. The erosion risk can be well evaluated.
In these previous studies, cavitation on hydrofoil surfaces is mainly found related to the local flow separation and pressure drop [13]. Under different flow conditions, cavitation performs in different types. If the size of cavitation is small and flow is stable, the cavitation bubble will also stably attach on the surface [14]. If local flow is unstable and the size of cavitation becomes large, the cavitation bubble will move with flow as a cloudy cavitation [15]. Generally, the status of cavitation on impeller blade or hydrofoil is strongly related to the blade (foil) profile, local flow stability and the size of cavitation itself. As another important factor, the nuclei fraction will also affect the status of cavitation [16]. However, this factor is usually ignored, because the nuclei fraction in water medium has only a slight influence on cavitation.
Currently, more and more pumping stations and water power stations are constructed on plateau [17][18][19]. There are many stations built above 3000 m or even above 4000 m. Thus, the nuclei fraction in water medium is completely different from the situation in low-altitude area. It may have a strong influence on cavitating flow. As a robust and effective way in studying cavitation, the computational fluid dynamics (CFD) tool is widely applied especially for water pumps and hydroturbines [20][21][22][23][24]. Based on CFD, different cavitation models can be supplied to predict the cavitating two-phase flow. Currently, most of the popular cavitation models are based on the Rayleigh-Plesset equation [25], which describes the bubble dynamics under the influence of pressure field. Kubota et al. [26] builta cavitation model based on the Rayleigh-Plesset equation to calculate the cavitation bubble size and vapor fraction. Considering the influence of compressibility and pressure pulsation caused by turbulence kinetic energy, Singhal et al. [27] supposed the full cavitation model as an improvement in the simulation of cavitating flow. Schnerr and Sauer [28] focused on the modeling of the mass transfer rate between gas phase and liquid phase. The Schnerr-Sauer model is built to treat the cavitating flow as a mixture of liquid and a large number of gas bubbles. Zwart, Gerber and Belamri [29] improved the cavitation prediction method based on the Kubota model. In the Zwart model, the volume fraction terms are corrected because gas volume will vary with the medium density. Kunz et al. [30] considered the difference between the form and collapse of cavitation bubble. They introduced the Kunz model, which describes the liquid-gas phase change and gas-liquid phase change in different ways. In general, all these models consider the influence of nuclei fraction on cavitation. Hence, it will be feasible to discuss the influence of nuclei fraction in simulating the cavitation at different altitude levels.
However, there is no article that analyzes the influence of nuclei fraction in simulating the cavitation in hydraulic turbomachineries. Cavitation in hydraulic turbomachinery is strongly relative to the internal complex flow regimes like vortex, back flow, jet-wake, stallcells and other pressure drop cases [31]. In this case, a twisted hydrofoil [32] was used as a simplified case to have a sensitivity analysis of nuclei fraction in simulating the cavitation at different altitude levels. Considering the altitude of the Qinghai-Tibet Plateau, which is over 4000 m, the air pressure is strongly different from that in the low-altitude area [33]. According to Henry's law, the nuclei fraction strongly varies from low-altitude area to a 4000 m plateau. This study will comparatively analyze the cavitation at high altitude level. With the nuclei fraction at 4000 m, the turbulent flow and development of cavitation over the twisted hydrofoil was analyzed in detail. It will provide a good reference for water pumps and hydroturbines for high-altitude levels considering the anti-cavitation issue.

Numerical Method for Turbulent Flow
Considering the strong adverse pressure gradient in the main flow region and the shear flow near the, the Shear Stress Transport (SST) model-based Detached Eddy Simulation (DES) method [34,35] is used as the numerical method for turbulent flow in CFD simulation. The k and ω equations of SST model can be written as follows:

∂(ρω) ∂t
where ρ is density, P is the turbulence production term, µ is dynamic viscosity, µ t is eddy viscosity, σ is model constant, C ω is the turbulence dissipation term, F 1 is the zonal blending function, l k−ω is called the turbulence scale that l k−ω = k 1/2 β k ω, where β k is the model constant. The DES method uses a zonal treatment term that min(l k−ω , C DES ·L mesh ), where L mesh is the maximum mesh element dimension and C DES is a constant. When l k−ω is larger than C DES ·L mesh , Large Eddy Simulation will be activated. Otherwise, SST model is used in the Reynolds-averaged mode.

Cavitation Model
The Zwart model [29] is a widely used cavitation model in current CFD simulations of cavitating two-phase flow. It is based on the Rayleigh-Plesset bubble dynamic equation to describe the mass transfer. It has the advantage that the sensitivity of cavitation vapor volume due to density change is considered. The requirement of computational cost when reasonable is used for engineering simulations. Therefore, it is applied as the cavitation model in this case. The rate of mass transfer . m can be expressed as: where p v is the saturation pressure, ρ v is the vapor density, f vnuc is the nuclei volume fraction, R B is the nuclei average radius, F e and F c are coefficients of evaporization and condensation which are commonly F e = 50 and F c = 0.01.

Vapor Volume Fraction
According to Henry's law, the gas dissolved as nuclei in liquid is proportional to the gas pressure over the liquid surface. Based on the data statistics on the Qinghai-Tibet Plateau [33], the atmosphere pressure P atm and altitude H alt have the approximate relationship as: where C d are constants that C d1 = 1.013 × 10 5 , C d2 = 1.259 × 10 1 , C d3 = 6.476 × 10 −4 . Based on Henry's law, the nuclei volume fraction f vnuc will be different if the altitude H alt is different. The original f vnuc at H alt = 0 m is 5 × 10 −4 . In this case, the altitude conditions that H alt = 0 m, 1000 m, 2000 m, 3000 m and 4000 m are comparatively studied. Therefore, the values of f vnuc are listed in Table 1.

Important Dimensionless Parameters
In this study, the hydrofoil case was discussed. To investigate the cavitating flow, two dimensionless parameters are defined as follows. Firstly, the fundamental parameter in the cavitation description is the cavitation coefficient C σ : where p v is the saturation pressure, p ref and v ref are reference pressure and velocity, respectively, which are usually measured at the upstream of hydrofoil, ρ is density. Therefore, C σ can be adjusted by changing the value of p ref . Secondly, the pressure coefficient C p can be defined as: where p is pressure. Therefore, the same C σ values can be compared for different altitudes with considering that the altitude caused a difference on nuclei volume fraction f vnuc .

Flow Domain of Hydrofoil
In this study, the Delft Twist 11 hydrofoil [32,36] was used based on the NACA0009 profile. This symmetric NACA 4-digit foil profile can be expressed as follows: where c is the total camber length, X is the position along camber direction, Y is the thickness, Y m is the maximum thickness of foil. Parameter C a and C b are constants that C a = 0.2, C b1 = 0.2969, C b2 = 0.126, C b3 = 0.3516, C b4 = 0.2843 and C b5 = 0.1015. This twisted hydrofoil has a different installation angle α at a different span. It has the law of α as: where s is the non-dimensionalized spanwise position against the chord length c and here 0 ≤ s ≤ 2. α m is the maximum installation angle at mid-span, which is 11 degrees. α w is the installation angle on the wall side, which is −2 degrees in this case. The 3D flow domain around the twisted hydrofoil is shown in Figure 1, which location in the Cartesian x-y-z coordinate is used for CFD simulation. As shown in Figure 2

CFD Setup
The CFD was conducted following the sequence of meshing, setting and solving. The mesh was independence-checked, as shown in Figure 2. The criterion is that the residual of lift/drag ratio is less than 0.001. The final mesh is in the structured type with 2,449,280 nodes and 2,380,744 hexahedral elements. Mesh in the near-wall region was controlled for wall-function where y + was from 0.47 to 23.75. As introduced above, the DES method and Zwart cavitation model were used in this numerical study. The fluid is water at 20 • C. As indicated, boundary conditions are set on the domain including a velocity inlet, a pressure outlet, a symmetry boundary at mid-span, a no-slip wall on foil surface, slip wall boundaries on upper wall, lower wall and side wall. The inlet velocity v in is 6.97 m/s, which means that the Reynolds number Re is 1.05×10 6 . The reference location for C p and C σ is the inlet boundary. To have a better study of cavitation, the same C σ situations are compared for different altitude levels. For a specific value of v in , the inlet-outlet pressure difference is almost unchanged. Thus, p ref at inlet can be adjusted to an expectable value by setting a specific pressure value at outlet. Steady-state simulation is firstly conducted. It will converge after the RMS residual of momentum and continuity equation is less than 1 × 10 −4 or finish after 1000 iterations. Transient simulation is conducted based on steady-state simulation. The total time is 1s and the time step is 1 × 10 −5 s. The maximum iteration number for each time step is 10. The convergence criterion is also RMS residual less than 1 × 10 −4 .

Numerical-Experimental Verification
Before analyzing the CFD results, it is necessary to verify the simulation. Figure 3 shows the comparison of pressure coefficient C p between CFD prediction and experimental data [12,36]. Three different spanwise positions in which s = 0.6, 0.8 and 1.0 are compared, especially focusing on the low-pressure side where cavitation usually occurs. The CFD predicted C p curves are accurate on the three spanwise surfaces. The CFD simulation can be used for further analyses of the cavitating flow in the following sections.

Variation Law
To study the sensitivity of nuclei fraction at different altitudes H alt , the cavitation vapor proportion f cav in the fluid domain is defined as: where V cav is the cavitation vapor volume and V fluid is the fluid domain volume. Figure 4 shows the comparison of cavitation vapor proportion f cav among different H alt at different C σ . With the decreasing of C σ from 2.713 to 1.071, f cav continually increases to a high level. The smaller C σ is, the quicker f cav increases. It represents the increasing of cavitation vapor in the entire fluid domain. However, there are differences among different H alt situations. Figure 5 shows the cavitation vapor proportion among different H alt at specific C σ values. The tendency is similar in all the 9 situations that f cav decreases with the increasing of the H alt level.

Sensitivity Analysis
Because the size of cavitation is different at different H alt , it is necessary to analyze the sensitivity of f cav on H alt and f vnuc . The difference between maximum and minimum f cav among different H alt at a specific C σ is defined as the sensitivity ∆f cav . For a better comparison, the relative sensitivity ∆f cav * is defined as: where f cav a is the average vapor proportion of all the 5 H alt situations. Figure 6 shows the sensitivity analysis of f cav at different C σ . The H alt -average vapor proportion f cav a shows the same variation tendency as in Figure 4. f cav a increases with the decreasing of C σ . The sensitivity ∆f cav has almost the similar tendency of f cav a . The smaller C σ is, the greater the difference is among different altitudes. However, there is a special local peak region when ∆f cav drops to a low level around C σ = 1.8. It means that the difference between maximum and minimum f cav is locally higher.
When the size of cavitation is small at high C σ , the absolute difference of f cav among different H alt is also small. Therefore, it is necessary to compare the variation of relative sensitivity ∆f cav * . As shown, the tendency is completely different. It is in a W-shape with two slowly variating regions and two rapidly rising regions, as indicated. The first rapidly rising region is about C σ = 1.5~1.9. The secondly rapidly rising region is about C σ = 2.5~2.7. Generally, when the size of cavitation is small, the relative sensitivity ∆f cav * is higher. In the high altitude plateau area, it is necessary to consider the influence of altitude level on cavitation inception.

Pressure Distribution Law on Foil Surface
After analyzing the influence of altitude level on the simulation of cavitating flow, it is necessary to study the flow behaviors around the hydrofoil. As is commonly known, cavitation on hydrofoil is strongly related to the pressure distribution. Figure 7 shows the distribution of pressure coefficient C p on different spanwise positions (0 ≤ s ≤ 1) of foil without considering cavitation. The maximum pressure coefficient C pmax is similar (about 1.0) for different s. This high pressure is because of the local flow striking on the foil lower surface, as shown in Figure 8. The minimum pressure coefficient C pmin varies with s. The larger s is, the smaller C pmin is. This is because of the local flow separation on the foil upper surface, as is also shown in Figure 8.  To have a better comparison, Figure 8 shows the minimum pressure coefficient C pmin and installation angle α at different spanwise s positions. There is a significant inverse relationship between C pmin and α. The larger the installation angle is, the larger the flow incidence angle is. It indicates the stronger and stronger flow separation and pressure drop when incidence angle is increasing.

Turbulent Flow around Foil
To have a better understanding of the pressure drop, Figure 9 shows the velocity coefficient C v at different spanwise s positions with indication of vectors. The uniformed velocity C v is defined as:  (11) where v is velocity and v in is the velocity at inlet. At s = 0.6, installation angle α is about 5.13 degrees. Flow attaches well on the foil surface. There are two typical low C v regions. Firstly, it is the local flow striking region on the leading-edge on the foil lower surface. Secondly, it is the wake region downstream to the foil trailing-edge. At s = 0.8, installation angle α increases to about 7.86 degrees. An obvious flow separation region occurs on the foil upper surface with low C v . The leadingedge striking region and the trailing-edge wake region are wider. At s = 1.0 (mid-span), installation angle α increases to 9 degrees, which is relatively large. It is obvious that the flow separation region on the foil upper surface is much wider. The leading-edge striking region and the trailing-edge wake region are also wider.
Considering the flow separation and sudden pressure drop at s = 1.0 (mid-span), it is necessary to analyze the local vortex shedding pattern, especially on the foil upper surface. In this case, the velocity helicity H v is used [37]. It can be non-dimensionalized to the velocity helicity coefficient C vhe by: where g is the acceleration of gravity. Figure 10 shows the velocity helicity coefficient C vhe on the mid-span plane and foil upper surface. The vortex-shedding phenomenon can be seen with the indication of the vortex-shedding route (VSR). On the mid-span plane, obvious VSR can be found from leading-edge (LE), along the upper surface, to trailing-edge (TE) and towards downstream. On the foil upper surface, VSR is complex and mainly along the diagonal line from the LE-mid-span corner to the TE-wall corner. Generally speaking, local vortical flow moves along the slope of the twisted foil surface.

Development of Cavitation at H alt -= 4000 m
Considering the plateau environment that H alt = 4000 m and f vnuc = 3.01×10 −4 , the cavitating flow is simulated and analyzed in Figure 11. The development of cavitation is comparatively studied from C σ = 2.713 to C σ = 1.071. The scale of cavitation continually increases in different directions. Figure 11 mainly shows the region of cavitation covering on the foil upper surface. Leading-edge (LE) is on the left side. Two parameters are defined to have a quantitative comparison. One is the length of cavity-covered area l cav and another is the width of cavity-covered area w cav . Figure 11 also includes the mid-span view of cavitation. Two parameters are defined in this view. One is the maximum thickness of cavity on mid-span t cav . Another is the total length of the attached part of the cavity on mid-span, which is denoted as l cav * . In general, l cav , w cav , t cav and l cav * increase with the decreasing of cavitation coefficient C σ . To have a better analysis, the variation of l cav , w cav , t cav and l cav * are compared in Figure 12. These four parameters are normalized against the foil chord length c. The growth rate dϕ/dC σ is also analyzed between each two conditions.
For the length of cavity-covered area on foil upper surface l cav , the growth rate is relatively low, within C σ = 2.173~1.784. The value of d(l cav /c)/dC σ is lower than 0.05. The value of l cav /c increases from about 0.026 to about 0.051. When C σ is smaller than 1.784, the growth rate of l cav becomes much higher. The value of d(l cav /c)/dC σ is about 0.09~0.40. From C σ = 1.784 to C σ = 1.071, the value of l cav /c strongly increases from about 0.051 to 0.202. Cavity covers about 1/5 of the foil surface along x direction at C σ = 1.071.
For the width of the cavity-covered area on foil upper surface w cav , the growth rate d(w cav /c)/dC σ is stable around 0.3. The increasing of w cav /c against C σ is almost linear. From C σ = 2.173 to C σ = 1.071, w cav /c increases from about 0.158 to about 0.636. At C σ = 1.071, cavity covers more than 3/5 of the half foil surface along the y direction. For the maximum thickness of cavity on mid-span t cav , the growth rate is relatively stable. The value of d(t cav /c)/dC σ is lower than 0.03. The increasing of t cav /c against C σ is almost linear, except in a small range between 1.784 and 1.480. In this range, a special phenomenon occurs in which a tail can be seen on the profile of cavity. This is because of the backward-jet flow (indicated in Figure 13 as an example), and the cavity becomes much thicker. From C σ = 2.173 to C σ = 1.071, t cav /c increases from about 0.001 to about 0.023. This thickness (t cav /c = 0.023 at C σ = 1.071) is about 1/4 of the maximum foil thickness. For the total length of the attached part of cavity on mid-span l cav * , the growth rate is similar to l cav . From C σ = 2.713 to C σ = 1.480, the value of d(l cav * /c)/dC σ is lower than 0.2. The value of l cav * increases from about 0.019 to 0.162. From C σ = 1.480 to C σ = 1.071, the value of d(l cav * /c)/dC σ is 0.36~0.44. The value of l cav * /c obviously increases from about 0.162 to 0.326. Comparing with l cav , both the value and the growth rate of l cav * is higher. This is also because of the tail on cavity. The total length of cavity is bigger than the cavity length covering on foil surface.
Generally, cavitation on a twisted hydrofoil extends along the streamwise direction, spanwise direction and thickness direction with the decreasing of cavitation coefficient C σ . When C σ is at a higher level, the small-scale cavity is attached on foil surface. When C σ becomes lower, cavity on large-incidence-angle spans is broken by the backward-jet from small-incidence-angle spans. A tail is generated on the cavity and the cavity becomes relatively unstable. The growth of cavity becomes quicker, especially in streamwise (length) direction.

Conclusions
Based on the above studies and discussions, the conclusions can be drawn as the following two main points: (1) With the decreasing of cavitation coefficient C σ , the scale of cavitation continually increases and the increasing is quicker and quicker. The nuclei volume fraction f vnuc has obvious influence on cavitation. The size of cavitation is different at different altitude levels. If the altitude is higher within 0~4000 m, the f vnuc is lower and the size of cavitation is smaller. The difference of the size of cavitation among altitude levels is bigger when C σ is small. That is, the sensitivity ∆f cav is high. On the contrary, the relative sensitivity ∆f cav * , which is the ratio between ∆f cav and the absolute cavitation fraction f cav , is high when C σ is large. When C σ is 1.071, the ∆f cav * between 0 m and 4000 m altitudes is about 4.6%. When C σ increases to 2.713, the ∆f cav * can be up to about 22.8%. It means that the cavitation volume fraction sensitivity should be considered in judging the inception cavitation of water pumps and hydro-turbines in the plateau environment.
(2) For this twisted hydrofoil, the installation angle and flow incidence angle are different at different spans. The incoming flow will cause local high pressure on the lower surface of hydrofoil. There will be a local low pressure site on the foil upper surface due to flow separation. This low pressure will cause cavitation. From sidewall to mid-span, the installation angle increases and the minimum pressure decreases. With the decreasing of C σ , the size of cavitation extends along the spanwise direction, streamwise direction and thickness direction. The growth rate is high in the spanwise (cavity width) and streamwise (cavity length) directions and low in thickness direction. When the size of cavitation is large enough, it will be broken by backflow-jet flow. A tail generates and the cavity becomes relatively unstable.
In general, this study focused on the sensitivity and influence of nuclei fraction on cavitation. The cavitating flow on a twisted hydrofoil was studied in detail. It is helpful for the anti-cavitation design of water pump units and hydro-turbine units installed on plateau.

Data Availability Statement:
The data used in this study are available upon request, from the corresponding author.