Abstract
The existence of runner interblade vortices can cause instability problems in Francis turbines, for example, pressure fluctuations, vibrations of runners, noise, and so forth. It is favorable in engineering practice to have the knowledge of the appearance of the incipient and developed inter blade vortex lines on the Hill diagram in unit parameters. Most numerical research on the inter blade vortices has been focused on the study of the characteristics of pressure fluctuations by single-phase flow calculations. However, since the two vortex lines are distinguished by observations of visible vortices, which contain cavitating flows, it is clear that cavitation calculations are needed for their predictions. A method by solving RANS equations with RNG k – ∊ turbulence model and ZGB cavitation model was proposed for the predictions in a Francis model turbine. Modifications of the turbulence viscosity were made for better simulation results. Vorticity criterion was chosen to identify the vortices. The fact that the results of cavitation calculations have a better agreement with experimental observations than single-phase calculations proves the validity of this prediction method.
1. Introduction
The presence of vortices is one of the main causes of hydraulic instability, that is, pressure fluctuations in hydraulic reaction turbines. The phenomenon of runner inter blade vortices, which is induced by high incidence at the runner inlet and driven through the blade channel toward the runner outlet at excessive off-design operating conditions, has been recognized to be responsible for several engineering accidents in Francis turbines. For instance, severe vertical vibrations (over 200 μm) in the guide bearings, as well as major horizontal cracks (∼1 m) in the draft tube cones, have been observed on unit numbers 13 and 14 of Tarbela power station in Pakistan at high head after about 6 months of their commissioning, which were pointed out by the manufacturer, DBS-Escher Wyss in Canada, to be associated with inter blade vortices at high head in the units [1]. Alstom has specified the cracks on the runners of two hydroelectric power stations in Brazil as induced by the inter blade vortices at low head, through the experimental analysis of the model and prototype runners [1]. Extensive noise in the rehabilitated turbines in Gongzui power station has been reported, and the source has been verified to be the inter blade vortices at low load operating conditions [2].
Due to the strong influence of the existence of the inter blade vortices on the stable operation of Francis turbines, it becomes prevalent to monitor the incipience and development of the vortices in the model acceptance tests of the runner, especially for large scale units. A rule of thumb has been applied to determine the incipience of the inter blade vortices in the model acceptance tests, by identifying visible vortices in two or three inter blade flow channels simultaneously. And the conditions when visible vortices occur in all the inter blade flow channels are defined as with developed inter blade vortices. Lines connecting the operating points on the hill diagram in unit parameters (Q11 ∼ n11, Q11 is the unit discharge, and n11 is the unit rotating speed) with incipient and initial occurrence of developed inter blade vortices are called the incipient inter blade vortex line, and developed inter blade vortex line, respectively, as shown in Figure 1. Normally it is suggested to avoid extensive inter blade vortices, that is, the developed inter blade vortex line in the normal operating range of the runner. Hence it is favorable to acquire knowledge of the characteristics of the vortices, as well as the appearance of the incipient and developed vortex lines beforehand.

Hill diagram in unit parameters.
Research has been carried out on the behavior of the inter blade vortical flows. Grindoz has made observations of the visible inter blade vortices of a model Francis turbine [3]. A series of experimental studies on the flows and pressure fluctuations in the model runner for the Three Gorges power station has been undertaken, which conclude that the inter blade vortices occur from the runner inlet at high head working conditions [4, 5]. No distinct frequency of the vortices was identified through the spectrum analysis of the results of pressure fluctuations. However, it was indicated that the existence of the vortices could induce alternating loads on the runner blades, resulting in a fatigue failure of them. Numerical studies have been focused on the predictions of the appearance of the vortices, and the analysis of the characteristics of the induced pressure fluctuations. Most of these research were carried out through single-phase flow simulations, and no conclusive frequency characteristics of the pressure fluctuations induced by the inter blade vortices have been obtained [6–8].
It has been pointed out that the visible vortices in the inter blade flow channels in the model tests are cavitating flows, that is, inter blade vortex cavitation [9]. Thus it is clear that single-phase flow simulations are not able to model the complex process of this phenomenon. Kurosawa et al. have provided a numerical model to predict the major characteristics, including the occurrence of the inter blade vortices, by solving Reynolds Averaged Navier-Stokes (RANS) equations combined with Reynolds Stress model and Volume of Fluid (VOF) method [10]. Reasonable agreements have been achieved between numerical and experimental results.
This paper is focused on the prediction of the incipient and developed runner inter blade vortex lines of a model Francis turbine. A numerical method based on RNG k – ∊ turbulence model and ZGB cavitation model has been applied for this purpose. Modifications were made on the turbulence viscosity for better simulation results in cavitation regions. Comparisons with single-phase flow simulations have been made, and the results are presented below. As far as the authors are aware, this is the first publication addressing the prediction method of the two vortex lines in hydraulic machinery.
2. Numerical Methods
2.1. Cavitation Simulation
Simulations of cavitating flows can be achieved by assuming the vapor/liquid mixture as a multiphase single fluid, with variable densities. No slip exists between the two phases in this mixture model. Thus only one set of momentum equations are needed to describe the pressure and the velocities of the mixture [11, 12]. Different forms of cavitation model have been proposed to account for the mass transfers between the two phases in the additional continuity equation of vapor phase, for example, cavitation model by Merkle et al. [13], Kunz et al. [14, 15], Singhal et al. [16], Senocak and Shyy [17], and so forth. The so-called ZGB (Zwart-Gerber-Belamri) cavitation model was applied in this study, which considers the variation of the volume fraction of the cavitation nuclei and evaluates the rate of mass transfers through a simplified Rayleigh-Plesset equation [18, 19]. The continuity equation for vapor phase and the expressions of the two source terms, R e and R c (specifying the mass transfer of the evaporation and condensation processes, resp.), of the model are as follows:
where sgn in sgn(p v – p) is the sign function; α v , ρ v represent the volume fraction and density of the vapor, respectively; p v denotes the vapor pressure; αnuc is the nucleation site volume fraction; R B is interpreted as the radius of a nucleation site; F e , F c are two empirical calibration coefficients.
As stated above, RANS equations were solved by adapting RNG k – ∊ turbulence model and ZGB cavitation model in this study. Modifications of the turbulence viscosity have been applied for predicting different flows [20]. It has been indicated that turbulence viscosity is usually overestimated in cavitation simulations by using mixture model and modifications of temporal and spatial discretization schemes alone are not able to neutralize this effect [21]. Referring to the research of Coutier-Delgosha [22, 23], a modification of the turbulence viscosity μ t from Boussinesq eddy viscosity assumption was adopted in this study,
where k, ∊ are turbulent kinetic energy and its dissipation rate, respectively; ρ v , ρ l , and ρ m represent the density of the vapor, the liquid, and the mixture, respectively; function f(ρ m ) is introduced to consider the influence of the variations in density on the turbulence viscosity. As shown in Figure 2, value of exponent n controls the shape of the f(ρ m ) curve against ρ m . n = 10 was chosen in this study as suggested by reference mentioned above.

f(ρ m ) ∼ ρ m .
2.2. Model Francis Turbine
The numerical studies were conducted with a model Francis turbine, whose parameters are listed in Table 1. The computational domain is shown in Figure 3.
Parameters of the model Francis turbine.

Computational domain.
The model gives a 94.54% level of peak efficiency. The performance characteristics determined by model tests under experimental head of 30 m can be seen in Figure 1. Systematic uncertainty on the efficiency measurement of the test rig is less than ±0.25%.
2.3. Computational Details
In order to investigate the validity of the methodology of predicting the incipient and developed inter blade vortex lines of the turbine by cavitation simulations, single-phase flow simulations by RANS RNG k – ∊ turbulence model were also carried out for comparisons. Due to the heavy computational loads, calculations of both cavitating and single-phase flows were performed in steady state.
The model's mesh system was composed of unstructured grids in the spiral case, the stay vanes and the runner, and structured grids in the guide vanes, the draft tube, and the inlet pipe. Commercial software package ICEM was used for mesh discretization. Local refinements to the boundary layer in the runner were applied, to ensure that the values of y+ are compatible with the chosen turbulence model (y+ < 200 on runner blades). Verification of grid independence was performed with respect to runner efficiency, with total number of grid cells that varied from 2 million to 8 million, as shown in Figure 4. It can be seen that the degree of computational precision is satisfactory when the number of cells is more than 6 million. A mesh with about 7.6 million cells in total was chosen for the simulations. The number of nodes and cells in each part is shown in Table 2.
Number of cells and nodes of each subdomain.

Verification of grid independence.
The three dimensional steady calculations were realized with the MRF (Multiple Reference Frame) model in Ansys CFX. A second-order upwind scheme was used for the discretization of the convective terms, and a second-order centered scheme was used for the source terms. Total pressure, initial values of the turbulent kinetic energy, and turbulent dissipation rate were set at the inlet boundary, while static pressure was set as the outlet boundary condition. For cavitation calculations, cavitation number σ = 0.086 was specified according to the experimental conditions.
3. Results and Discussions
In order to predict the incipient and developed inter blade vortex lines, it is essential to evaluate the compatibility of various vortex identification methods with our simulation cases. A number of vortex identification methods/criteria have been proposed [24], including vorticity criterion, Q criterion, swirling strength criterion, and helicity criterion. By comparing the isosurfaces of the four quantities at corresponding operating points with n11 = 59 r/min and n11 = 72.5 r/min, along the incipient and developed inter blade vortex lines, respectively, it can be concluded that the vorticity criterion is most suitable for distinguishing the vortices. As shown in Figures 5 and 6, distinct differences in the distribution of the isosurfaces of the vorticity can be observed, with vorticity = 250 s−1 with single-phase calculation and vorticity = 398 s−1 with cavitation calculation, between operating conditions along the incipient and developed inter blade vortex lines. Vortices in 2–3 interblade runner flow channels can be seen along the incipient line, and in most flow channels along the developed line.

Inter blade vortices on the incipient line.

Inter blade vortices on the developed line.
In order to predict the two vortex lines, seven and nine n11 values were chosen, each with 3∼6 operating conditions for single-phase and cavitation calculations, as shown in Tables 3 and 4, respectively. The incipient and developed inter blade vortex lines were asymptotically derived with vorticity criterion. The results of the predictions with both single-phase and cavitation calculations were compared with experimental observations, as shown in Figure 7. It can be concluded that better agreements were achieved by using the prescribed cavitation calculations. It is also seen that the inflection points on both lines were successfully predicted with cavitation calculations.
Operating conditions for single-phase calculations.
Operating conditions for cavitation calculations.

Predictions of the incipient and developed inter blade vortex lines.
Figure 8 shows the velocity triangles at runner inlet, with different operating heads. Higher head corresponds to larger value of absolute velocity c1, hence smaller n11; lower head corresponds to smaller value of c1, hence larger n11. The flow angle at optimum n11 is denoted as β1, which is approximately equal to the blade angle βb1. The angle of attack i = β1 – βb1 ≈ 0 at optimum n11 (no. 7 in Table 5); i > 0 at small n11, where separations on suction sides of the blades occur; i < 0 at large n11, where separations on pressure sides of the blades occur. Bigger angles of attack exist more along the developed inter blade vortex line than incipient line, as shown in Table 5. The existence of the inflection points on the two lines can be roughly explained as below. The occurrence of the inter blade vortices is associated with the variations of the flow angles at runner inlet. This angle tends to be positive with operating conditions at high head, causing flow separations at suction sides of the blades near crown (Figure 9(a)), and it tends to be negative at low head, causing flow separations at pressure sides of the blades near crown (Figure 9(c)). In the case of inflection point (Figure 9(b)), a larger area of separation exists at the crown compared to the case of high head, which induces inter blade vortices at larger Q11.
Angles of attack.

Velocity triangles at runner inlet at different n11 (head).

Distributions of vorticity on runner.
4. Conclusions
In order to predict the incipient and developed inter blade vortex lines, a method of cavitation calculations was proposed, by solving RANS equations with RNG k – ∊ turbulence model and ZGB cavitation model. Modification of the turbulence viscosity was made to compensate its overestimation in cavitation area. The predictions of a Francis model turbine were derived by an asymptotic approach. By comparing the results with single-phase calculations and experimental observations, it was concluded that the prescribed cavitation calculations could achieve better predictions for the two vortex lines.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Footnotes
Acknowledgments
The authors would like to thank National Natural Science Foundation of China (no. 51076077) and National Science and Technology Ministry of China (ID: 2008BAC48B02) for their financial supports.
