Leakage Vortex Progression through a Guide Vane’s Clearance Gap and the Resulting Pressure Fluctuation in a Francis Turbine

: A clearance gap (CG) between guide vanes (GVs) and facing plates exists at both ends of a Francis turbine and allows the opening angle to be adjusted for varying operating conditions. Leakage ﬂow is induced through this gap due to the pressure difference between the two sides of the guide vanes. While some research works have used qualitative approaches to visualize and predict the strength of a leakage vortex (LV), this paper presents a method for quantifying vortices along a trajectory. In this paper, a prototype high-head Francis runner with speciﬁc speed of 85.4 is considered as a reference case. A systematic investigation across both space and time is carried out, i.e., analysis of the spatial temporal progression of LV for three operating conditions. While travelling from the CG to runner leading edge, LV evolution and trajectory data are observed and the values of vorticity and turbulent kinetic energy are calculated for the LV trajectory. Frequency spectrum analyses of pressure oscillations in the vaneless space, runner blade, and draft tube are also performed to observe the peak pressure pulsation and its harmonics. Unsteady ﬂuctuations of the runner output torque are ﬁnally studied to identify the patterns and magnitudes of torque oscillations.


Introduction
Francis turbines are the most prevalent reaction turbines, which are normally adopted for medium-head and medium-flow conditions. In these turbines, the flow rate is regulated with guide vanes (GVs) that exist upstream of a runner. Water leaving from the guide vanes enters the runner blades with the desired angle at a high velocity which is consequently converted into rotary motion. In addition, GVs are incorporated with a small clearance gap (CG) at both ends so that they can pitch around their axis, which enables the adjustment of the opening angle based on various operating conditions. Several studies [1][2][3][4][5] have suggested that the sizes of CGs increase with head cover deflection due to water pressure, which results in the disturbance of the runner inlet flow conditions and a loss in overall efficiency. Brekke [1] studied the influence of varying CG and GV sizes on flow behavior and turbine efficiency, observing that the size of the CG has a high influence on the overall turbine efficiency. Thapa [2] developed an experimental setup with a GV cascade representing the flow inside the distributor of a Francis turbine. Pressure and velocity data obtained from his experiments showed that all CG sizes greater than 1 mm induced a turbulent cross-flow jet which mixes with the main flow and disrupts the runner inlet flow conditions. A critical size for maximum leakage flow effects was identified to be 2 mm from his experiments. Chitrakar [3] investigated flow phenomena around GV parameters computationally and compared and validated their results with the experimental results from Thapa [2]. The results showed that flow on the suction side is significantly affected by pressure-to-suction side flow through the gap and that the reduction of the pressure gradient reduces leakage through GVs. Koirala [4] performed a study on GVs at the operating lifetime. Vortex evolution and propagation is a complicated phenomenon and there are several examples of vortex-induced instabilities in hydropower turbines [12]. The unsteady phenomenon in turbomachines can be classified into periodic and non-periodic types [13]. RSIs are associated with the periodic class, whereas vortex rope formation and cavitation, observed at the exit of the runner and in the cone of draft tube, which rely on the operating points, are associated with the non-periodic class. Vortex identification can lead to an increased understanding of complex flow phenomena. In complex flow scenarios, as in Francis turbines, there are several interacting vortices which may be bent or twisted, thus resulting in rotation in several different planes for each single vortex. Several systematic procedures for the identification of vortices have been developed and proposed [12][13][14][15][16]. Zhang [12] has performed a review of methods for vortex identification in hydraulic turbines. He has also discussed and summarized experimental techniques for vortex observation. Similarly, Qian [13] has experimentally investigated unsteady flows at the stator level in a Francis turbine model in his doctoral thesis. Chakraborty [14] proposed a local vortex identification criterion and requirements with the introduction of spiraling compactness for material orbits in vortices. Elsas [15] introduced an alternative vortex identification method which he termed as the "vorticity curvature criterion", which is fundamentally based on the local properties of the vortex field. He has carried out a realistic comparative study for vortex identification for the direct numerical simulation of a turbulent channel flow. Likewise, ANSYS ® has devised several methods to detect vortices as spatial regions with a specified set of equations [16]. The most popularly used criteria include the Q criterion [17], λ2 criterion, [18], Δ criterion [19], and swirling strength criterion (λci) [20]. The Q criterion is relevant for identification of a vortex in incompressible fluid, especially in large scale vortices in turbulent flow. The details of the Q criterion, which is the method of vortex identification in this paper, are discussed in the results section. The λ2 criterion corresponds to the pressure minimum in a plane when the contributions of unsteady irrotational straining and viscous terms in Navier-Stokes equations are discarded. The Δ criterion defines vortices in regions with a complex velocity gradient. A higher value of Δ denotes stronger spiraling within a vortex region. The λci criterion identifies a vortex structure based on the imaginary eigenvalues of the velocity gradient with the quantification of the swirling strength in the vortex. The spatiotemporal evolution of a leakage vortex (LV) in a Francis turbine is still unpredictable and the relationship between flow instability and LV progression has not been revealed completely. Hence, this work is intended to understand the spatiotemporal evolution of LVs and related steady and unsteady flow phenomena in a reference case of a Francis turbine. Moreover, the frequency spectra of pressure signals at various locations inside a turbine are analyzed to understand the RSI and vortex rope phenomena. Vortex identification can lead to an increased understanding of complex flow phenomena. In complex flow scenarios, as in Francis turbines, there are several interacting vortices which may be bent or twisted, thus resulting in rotation in several different planes for each single vortex. Several systematic procedures for the identification of vortices have been developed and proposed [12][13][14][15][16]. Zhang [12] has performed a review of methods for vortex identification in hydraulic turbines. He has also discussed and summarized experimental techniques for vortex observation. Similarly, Qian [13] has experimentally investigated unsteady flows at the stator level in a Francis turbine model in his doctoral thesis. Chakraborty [14] proposed a local vortex identification criterion and requirements with the introduction of spiraling compactness for material orbits in vortices. Elsas [15] introduced an alternative vortex identification method which he termed as the "vorticity curvature criterion", which is fundamentally based on the local properties of the vortex field. He has carried out a realistic comparative study for vortex identification for the direct numerical simulation of a turbulent channel flow. Likewise, ANSYS ® has devised several methods to detect vortices as spatial regions with a specified set of equations [16]. The most popularly used criteria include the Q criterion [17], λ2 criterion, [18], ∆ criterion [19], and swirling strength criterion (λci) [20]. The Q criterion is relevant for identification of a vortex in incompressible fluid, especially in large scale vortices in turbulent flow. The details of the Q criterion, which is the method of vortex identification in this paper, are discussed in the results section. The λ2 criterion corresponds to the pressure minimum in a plane when the contributions of unsteady irrotational straining and viscous terms in Navier-Stokes equations are discarded. The ∆ criterion defines vortices in regions with a complex velocity gradient. A higher value of ∆ denotes stronger spiraling within a vortex region. The λci criterion identifies a vortex structure based on the imaginary eigenvalues of the velocity gradient with the quantification of the swirling strength in the vortex. The spatiotemporal evolution of a leakage vortex (LV) in a Francis turbine is still unpredictable and the relationship between flow instability and LV progression has not been revealed completely. Hence, this work is intended to understand the spatiotemporal evolution of LVs and related steady and unsteady flow phenomena in a reference case of a Francis turbine. Moreover, the frequency spectra of pressure signals at various locations inside a turbine are analyzed to understand the RSI and vortex rope phenomena.

Methodology
This study uses a prototype turbine with nominal head size of 207 m, rated discharge of 4.33 m 3 /s, and rated speed of 750 rpm from the Bhilangana-III Hydroelectric Power Plant, Uttarakhand, India, which is mentioned as the reference turbine onwards. Figure 2 shows the methodology adopted for this research work.
Numerical studies were performed using ANSYS ® CFD TM by solving steady-state Reynolds-averaged Navier-Stokes equations. The calculations were conducted at 77%, 100%, and 130% loads and the respective flow rate values were 3.36, 4.34, and 5.64 m 3 /s [21]. A 100% load corresponds to the best efficiency point (BEP), whereas the 77% and 130% loads correspond to a partial load (PL) and full load (FL), respectively. Figure 3 shows a numerical model of the reference turbine, which has a runner rotating at 750 rpm and four  [22]. Shaft allowance was also included in the GVs for this study to make the numerical case as close as the prototype [22,23]. The full model of the turbine consists of 13 runner blades and 16 GVs, as shown in Figure 2. The draft tube mesh is not shown in the figure. Spiral casing and stay vanes were excluded in the domain during simulation.
Plant, Uttarakhand, India, which is mentioned as the reference turbine onwards. Figure 2 shows the methodology adopted for this research work.
Numerical studies were performed using ANSYS ® CFD TM by solving steady-state Reynolds-averaged Navier-Stokes equations. The calculations were conducted at 77%, 100%, and 130% loads and the respective flow rate values were 3.36, 4.34, and 5.64 m 3 /s [21]. A 100% load corresponds to the best efficiency point (BEP), whereas the 77% and 130% loads correspond to a partial load (PL) and full load (FL), respectively. Figure 3 shows a numerical model of the reference turbine, which has a runner rotating at 750 rpm and four stationary domains. The GVs of this powerplant were numerically studied by Acharya et al. earlier with a different runner [22]. Shaft allowance was also included in the GVs for this study to make the numerical case as close as the prototype [22,23]. The full model of the turbine consists of 13 runner blades and 16 GVs, as shown in Figure 2. The draft tube mesh is not shown in the figure. Spiral casing and stay vanes were excluded in the domain during simulation.

Mesh Generation and Boundary Conditions
Full passage modelling was carried out and the computational domain was discretized with a hexahedral structured mesh. ANSYS ® Turbogrid TM was applied to generate

Mesh Generation and Boundary Conditions
Full passage modelling was carried out and the computational domain was discretized with a hexahedral structured mesh. ANSYS ® Turbogrid TM was applied to generate the meshes in the runner and guide vanes, while the meshes in clearance gaps and draft tube were generated in ANSYS ® ICEM CFD TM . Meshes on the overall domain were well defined and the values of positive y-values were maintained below 30. A mass flow rate of 4330 kg/s at the circumferential inlet of guide vanes and average static pressure of 0 Pa at the outlet were chosen as boundary conditions for BEP case [24,25]. The designed mass flow rate of the turbine, representing 100% flow, was taken as the BEP. Interfaces between the stationary and rotational domain were considered with a frozen rotor method for steady-state analysis. This forms a local coupling which uses the rotating reference frame to save on computation by converting transient turbomachinery flow into a steady state. For transient analysis, a transient rotor-stator model was used as it accounts for all interaction effects between components that are in relative motion to each other and hence predicts the true transient interaction of the flow between a stator and rotor passage [26]. Transient simulations were carried out for total time corresponding to three runner revolutions, and at a time step of 1 • of angular rotation for runners. The turbulence model selected for this study was the shear stress turbulence model developed by Meter [27], which was determined based upon the literature [28][29][30]. The details of the simulation settings are listed in Table 1.

Mesh Sensitivity Analysis
The grid convergence index (GCI) method was used for the estimation of the discretization error and value extrapolation [31]. This particular technique has been found to be effective in predicting numerical uncertainties in the case of a Francis turbine [32,33]. For uncertainty analysis, three structured hexahedral meshes with different resolutions were created. Mesh refinement was performed by increasing the distribution in each direction. Uncertainties in the simulation were monitored by monitoring the torque and pressure inside the domain for three different points, i.e., VL1 (vaneless space), RV1 (runner blade) and DT1 (draft tube). The discretization error for the numerical model was determined as follows: (i) Average length of each element for a 3D mesh was determined as follows: (ii) Let h 1 < h 2 < h 3 and r 21 = h 2 /h 1 , r 32 = h 3 /h 2 . The apparent order was solved as in Equations (2)-(4) using a fixed point iteration method: (3) s = 1.sign(ε 32 /ε 21 ). (4) (iii) The extrapolated values were calculated as follows: (iv) The approximate and extrapolated relative errors were calculated as follows: Table 2 shows the uncertainties and extrapolated values. The numerical uncertainties in pressure at points 1 and 2, with a fine mesh and medium mesh, were 0.0151% and 0.2568%, respectively. Efficiency was calculated as the ratio of power output to input. Power output was calculated from the torque of the runner and its rotational speed, whereas power input was obtained from the net head. From the efficiency measurement, uncertainty for the fine mesh was 0.0307%.  Figure 4a shows the normalized velocity at the outlet of the GV along the circumferential position at the trailing edge of GVs. Velocity in this case is normalized in terms of net head of the turbine as expressed in Equation (9).
where v = local velocity at the outlet of the GV, which is given as For the numerical measurement of velocity, circumferential locations of the outlets of two GVs were considered such that the flow beneath the trailing edge of GV could be distinguished. A dimensionless term x/c was used in the x-axis, which denotes the position (x) from LE with respect to the chord length (c). In Figure 4, two different low velocity regions were observed due to the influence of the wake travelling from the trailing edge of the GV. Because of the influence of wakes and vortices travelling from the trailing edge, the uncertainties at these locations are higher, as shown in Figure 4b. As clearly observed in Figure 4a, the graph representing the coarse mesh has higher discrepancies as compared to the medium and fine meshes. Hence, the coarse mesh was discarded and the medium one was used for further numerical calculations. Figure 4a shows the normalized velocity at the outlet of the GV along the circumferential position at the trailing edge of GVs. Velocity in this case is normalized in terms of net head of the turbine as expressed in Equation (9). * [−] = (9) where v = local velocity at the outlet of the GV, which is given as = + .

Validation with Prototype Data
Steady-state simulations were conducted at designed and off-design conditions in the initial stage. Another set of simulations was carried out with the total pressure as the inlet boundary condition. The mass flow rate was calculated in a postprocessing step. Flow rate data available for the prototype turbine at different operating points [21] were compared with the results from the simulation as shown in Figure 5. It can be inferred from this figure that at the same guide vane opening angles, the resultant flow obtained from the CFD analysis closely matches the field data. Moreover, the efficiency comparison from the CFD calculation shows lower efficiency than the prototype measurement, which could be due to the overprediction of losses by the turbine with the use of the turbulence model [21]. It can be noticed at a deep partial load condition that the efficiency variation is higher than the designed point. In off-design conditions, the losses predicted by the numerical model are higher [7]. Figure 5 also shows the error bar which was constructed with the difference of numerical and experimental value for each GV opening angle. For the numerical measurement of velocity, circumferential locations of the outlets of two GVs were considered such that the flow beneath the trailing edge of GV could be distinguished. A dimensionless term x/c was used in the x-axis, which denotes the position (x) from LE with respect to the chord length (c). In Figure 4, two different low velocity regions were observed due to the influence of the wake travelling from the trailing edge of the GV. Because of the influence of wakes and vortices travelling from the trailing edge, the uncertainties at these locations are higher, as shown in Figure 4b. As clearly observed in Figure 4a, the graph representing the coarse mesh has higher discrepancies as compared to the medium and fine meshes. Hence, the coarse mesh was discarded and the medium one was used for further numerical calculations.

Validation with Prototype Data
Steady-state simulations were conducted at designed and off-design conditions in the initial stage. Another set of simulations was carried out with the total pressure as the inlet boundary condition. The mass flow rate was calculated in a postprocessing step. Flow rate data available for the prototype turbine at different operating points [21] were compared with the results from the simulation as shown in Figure 5. It can be inferred from this figure that at the same guide vane opening angles, the resultant flow obtained from the CFD analysis closely matches the field data. Moreover, the efficiency comparison from the CFD calculation shows lower efficiency than the prototype measurement, which could be due to the overprediction of losses by the turbine with the use of the turbulence model [21]. It can be noticed at a deep partial load condition that the efficiency variation is higher than the designed point. In off-design conditions, the losses predicted by the numerical model are higher [7]. Figure 5 also shows the error bar which was constructed with the difference of numerical and experimental value for each GV opening angle. Furthermore, validation of numerical model was also carried out by measuring the pressure around the middle GVs in three-GV cascade rigs developed by Chitrakar et al. [34]. The closed loop test setup was equipped with a pump delivering water from a lower reservoir. A flow meter was mounted at the outlet of the rig and two pressure taps were mounted at the inlet and outlet to acquire the correct operating points for measurement. The arrangement has a trapezoidal slot with holes with a diameter of 2 mm around the middle GV with a 2 mm offset from the GV surface. Piezoresistive pressure transducers Furthermore, validation of numerical model was also carried out by measuring the pressure around the middle GVs in three-GV cascade rigs developed by Chitrakar et al. [34]. The closed loop test setup was equipped with a pump delivering water from a lower reservoir. A flow meter was mounted at the outlet of the rig and two pressure taps were mounted at the inlet and outlet to acquire the correct operating points for measurement. The arrangement has a trapezoidal slot with holes with a diameter of 2 mm around the middle GV with a 2 mm offset from the GV surface. Piezoresistive pressure transducers were connected to these holes, which were used to determine GV loading. Figure 6 shows the experimental test section of the three-GV cascade rig that consists of a pressure measurement location at the middle GV (GV2). A pressure transducer was calibrated with a dead weight calibrator with an uncertainty of ±0.2%. Normalized pressure (Cp), which is the ratio of the pressure at a point to the pressure at an inlet as found by CFD analysis, was validated with the data obtained from the pressure transducer measurement around the circumferential positions from the leading edge to the trailing edge. Local pressure measurement at each measurement location was normalized by the pressure value at the leading edge of GV. Thus, the blade loading characterizing the flow around the GV was observed experimentally with a maximum normalized pressure where Cp = 1 at the leading edge [35]. In the CFD pressure measurement, at same experimental pressure measurement locations, pressure values were measured. It was seen that using the current numerical model, the maximum error in pressure measurement was 8.8% towards the leading edge in the pressure side of the GV. This might be due to the influence of GV1 and inappropriate stagnation as compared to the numerical solution. At all other locations of pressure measurement, the errors in the pressure measurement were less than 5%. The blade loading distribution of the current GV profile in the reference case is shown in Figure 7b.  Normalized pressure (Cp), which is the ratio of the pressure at a point to the pressure at an inlet as found by CFD analysis, was validated with the data obtained from the pressure transducer measurement around the circumferential positions from the leading edge to the trailing edge. Local pressure measurement at each measurement location was normalized by the pressure value at the leading edge of GV. Thus, the blade loading characterizing the flow around the GV was observed experimentally with a maximum normalized pressure where Cp = 1 at the leading edge [35]. In the CFD pressure measurement, at same experimental pressure measurement locations, pressure values were measured. It was seen that using the current numerical model, the maximum error in pressure measurement was 8.8% towards the leading edge in the pressure side of the GV. This might be due to the influence of GV1 and inappropriate stagnation as compared to the numerical solution. At all other locations of pressure measurement, the errors in the pressure measurement were less than 5%. The blade loading distribution of the current GV profile in the reference case is shown in Figure 7b. Normalized pressure (Cp), which is the ratio of the pressure at a point to the pressure at an inlet as found by CFD analysis, was validated with the data obtained from the pressure transducer measurement around the circumferential positions from the leading edge to the trailing edge. Local pressure measurement at each measurement location was normalized by the pressure value at the leading edge of GV. Thus, the blade loading characterizing the flow around the GV was observed experimentally with a maximum normalized pressure where Cp = 1 at the leading edge [35]. In the CFD pressure measurement, at same experimental pressure measurement locations, pressure values were measured. It was seen that using the current numerical model, the maximum error in pressure measurement was 8.8% towards the leading edge in the pressure side of the GV. This might be due to the influence of GV1 and inappropriate stagnation as compared to the numerical solution. At all other locations of pressure measurement, the errors in the pressure measurement were less than 5%. The blade loading distribution of the current GV profile in the reference case is shown in Figure 7b.   Figure 8a shows the overall pattern of the GV leakage vortex obtained from the transient simulation for the BEP condition whereas Figure 8b shows three distinct leakage flow (LF) areas. Part 1 shows a LF initiating from the leading edge of GV, whereas part 2originates from the GV shaft region. From Figure 8b, it is clearly observed that LF from part 1 mixes with part 2 which ultimately strikes the leading edge of runner. Similarly, another LF is set up from the trailing edge of GV, which is further divided in part 3a and part 3b in Figure 8b. The intensity of these vortices can be related with the guide vane loading curve shown in Figure 7b. The pressure difference between the two sides of GV is maximum in the stream-wise location (x/c) of 0.5-0.7. As a result, the leakage flow is driven with higher acceleration, increasing the intensity of the vortex filament. Due to the difference in number of GVs and runner blades in the turbine, the region of vortices hitting the runner is different for each blade and in every revolution.  Figure 8a shows the overall pattern of the GV leakage vortex obtained from the transient simulation for the BEP condition whereas Figure 8b shows three distinct leakage flow (LF) areas. Part 1 shows a LF initiating from the leading edge of GV, whereas part 2originates from the GV shaft region. From Figure 8b, it is clearly observed that LF from part 1 mixes with part 2 which ultimately strikes the leading edge of runner. Similarly, another LF is set up from the trailing edge of GV, which is further divided in part 3a and part 3b in Figure 8b. The intensity of these vortices can be related with the guide vane loading curve shown in Figure 7b. The pressure difference between the two sides of GV is maximum in the stream-wise location (x/c) of 0.5-0.7. As a result, the leakage flow is driven with higher acceleration, increasing the intensity of the vortex filament. Due to the difference in number of GVs and runner blades in the turbine, the region of vortices hitting the runner is different for each blade and in every revolution.  Figure 9 shows the LV trajectory obtained by transient simulation. Among several vortex identification methods as discussed in the introduction part above, the Q criterion is considered for this work. It is usually suitable for vortex identification in incompressible fluid, especially in large scale vortex in turbulent flow [12]. The Q criterion identifies vortices using second invariant of velocity gradient tensor which is defined as follows:

Leakage Vortex Progression
where is the velocity gradient tensor and tr is the trace of the matrix. A total of 13 points (P1 to P13) were located based on the vortex trajectory, as shown in Figure 9a, and vorticity was calculated on those points. Figure 9b shows the values of vorticity along LV trajectory for three operating conditions. Vorticity value was normalized in correspondence with BEP case, as shown in y axis of Figure 9b. Considering the propagation of LV from point P1, vorticity increases until point P2 for all cases and after point P3, it sharply drops until point P6. From point P6 to P9, it decreases slowly, and from P10 to P12, it is vorticity is almost lowest for the BEP and full load conditions as compared to the partial load condition from point P6 to P13 before striking the runner.  Figure 9 shows the LV trajectory obtained by transient simulation. Among several vortex identification methods as discussed in the introduction part above, the Q criterion is considered for this work. It is usually suitable for vortex identification in incompressible fluid, especially in large scale vortex in turbulent flow [12]. The Q criterion identifies vortices using second invariant of velocity gradient tensor which is defined as follows:

Leakage Vortex Progression
where ∇v is the velocity gradient tensor and tr is the trace of the matrix. A total of 13 points (P1 to P13) were located based on the vortex trajectory, as shown in Figure 9a, and vorticity was calculated on those points. Figure 9b shows the values of vorticity along LV trajectory for three operating conditions. Vorticity value was normalized in correspondence with BEP case, as shown in y axis of Figure 9b. Considering the propagation of LV from point P1, vorticity increases until point P2 for all cases and after point P3, it sharply drops until point P6. From point P6 to P9, it decreases slowly, and from P10 to P12, it is vorticity is almost lowest for the BEP and full load conditions as compared to the partial load condition from point P6 to P13 before striking the runner.  Transient simulations were performed for three operating conditions to inspect spatial-temporal progression of LV before it hits the leading edge of runner blades. Figure  11a-f show the evolution of the LV as one runner blade passes from one GV to the consecutive GV. It shows the vortex for a 10° revolution, capturing 2° at a time. For the investigation of the progression of the part II LV of the hub region, which strikes the leading edge of runner, it can be categorized into three stages, namely (a) the elongation stage, (b)  Transient simulations were performed for three operating conditions to inspect spatial-temporal progression of LV before it hits the leading edge of runner blades. Figure  11a-f show the evolution of the LV as one runner blade passes from one GV to the consecutive GV. It shows the vortex for a 10° revolution, capturing 2° at a time. For the investigation of the progression of the part II LV of the hub region, which strikes the leading edge of runner, it can be categorized into three stages, namely (a) the elongation stage, (b) Transient simulations were performed for three operating conditions to inspect spatialtemporal progression of LV before it hits the leading edge of runner blades. Figure 11a being slender and weaker as seen in Figure 11b,c. Figure 11d,e shows the disintegration stage, in which LV becomes more stretched with a decrease in volume and length. Only a few traces of the LV can be observed in Figure 11f, which can be described as the dissolving stage. The same process is expected to continue as the runner rotates. The LV at the shroud region is shown by the green ellipse in Figure 11a and is transported all the way until the mid-section of the blade before the majority of the vortex is dissolved.
Energies 2021, 14, x FOR PEER REVIEW 12 of 21 disintegration stage, and (c) dissolving stage. The elongation stage represents some shedding and the LV being slender and weaker as seen in Figure 11b,c. Figure 11d,e shows the disintegration stage, in which LV becomes more stretched with a decrease in volume and length. Only a few traces of the LV can be observed in Figure 11f, which can be described as the dissolving stage. The same process is expected to continue as the runner rotates. The LV at the shroud region is shown by the green ellipse in Figure 11a and is transported all the way until the mid-section of the blade before the majority of the vortex is dissolved. Figure 12a-f show the propagation of LV during 10° of rotation of the runner during the partial load condition. It also follows the same trend as in case of BEP, but the vortex disintegrates and dissolves earlier than BEP. For this condition, vortex almost disappear in Figure 12e but while compared to the same position for BEP, there is still some traces of vortex during that time. In partial load conditions, the pressure difference between the two sides of the GVs rises (because of the closing), so logically the intensity of the vortex should also rise. In this case, we saw that the intensity dropped instead, and this may be due to the reduced flow rates in the part load condition. Figure 12a-f show the propagation of LV during 10 • of rotation of the runner during the partial load condition. It also follows the same trend as in case of BEP, but the vortex disintegrates and dissolves earlier than BEP. For this condition, vortex almost disappear in Figure 12e but while compared to the same position for BEP, there is still some traces of vortex during that time. In partial load conditions, the pressure difference between the two sides of the GVs rises (because of the closing), so logically the intensity of the vortex should also rise. In this case, we saw that the intensity dropped instead, and this may be due to the reduced flow rates in the part load condition.
Similarly, the progression of the leakage vortex for full load conditions is depicted in Figure 13a-f. In this case, the strength and volume of LF at the hub and shroud is higher than in previous cases. Some amount of a vortex still exists before another set of LF from consecutive GV hits the runner blade. At full load conditions, due to opening, the pressure difference should be minimum; however, in this case, the flow rate is high, and this could be the reason why we see high intensity vortices compared to other cases.
The formation of complex vortices can represent significant influences on the performance and flow stability of hydraulic machinery. Insertion of points in the trajectory which corresponds to LV progression helped to compute the value of vorticity and TKE which are the vital parameters with major influence on the flow pattern. Intensity of these vortices can be related to GV loading curve. Similarly, the strength of turbulence in the flow can be related to the value of the TKE. Energies 2021, 14, x FOR PEER REVIEW 13 of 21 Similarly, the progression of the leakage vortex for full load conditions is depicted in Figure 13a-f. In this case, the strength and volume of LF at the hub and shroud is higher than in previous cases. Some amount of a vortex still exists before another set of LF from consecutive GV hits the runner blade. At full load conditions, due to opening, the pressure difference should be minimum; however, in this case, the flow rate is high, and this could be the reason why we see high intensity vortices compared to other cases. Similarly, the progression of the leakage vortex for full load conditions is depicted in Figure 13a-f. In this case, the strength and volume of LF at the hub and shroud is higher than in previous cases. Some amount of a vortex still exists before another set of LF from consecutive GV hits the runner blade. At full load conditions, due to opening, the pressure difference should be minimum; however, in this case, the flow rate is high, and this could be the reason why we see high intensity vortices compared to other cases. The formation of complex vortices can represent significant influences on the performance and flow stability of hydraulic machinery. Insertion of points in the trajectory which corresponds to LV progression helped to compute the value of vorticity and TKE which are the vital parameters with major influence on the flow pattern. Intensity of these

Pressure Pulsations Inside Runner
The stationary domain experiences pressure fluctuations at a frequency corresponding to the number of rotating blades (z b ) and the rotational speed of the runner (n), as shown in Equation (11). Likewise, rotating domain experience a frequency (f gv ) corresponding to number of guide vanes (z gv ) and rotational speed of runner (n), as shown in Equation (12) [33,36,37].
A total of sixteen points were located in the vaneless space, runner blade and draft tube, as shown in Figure 14a,b, to analyze frequency spectra. Among the sixteen points, six points from DT1 to DT6 were located in the draft tube passage from the runner outlet to the draft tube outlet. Similarly, four points, from VL1 to VL4, starting from hub area to shroud area, were defined in the vaneless space before the leading edge of runner. The remaining six points, namely, RV1 to RV6, starting from leading edge to trailing edge of the blade, were inserted along the runner blade. Pressure pulsation was calculated by subtracting the mean pressure (p) from the instantaneous pressure (p) and was normalized by the reference pressure (ρE) BEP [38]. In a high-head Francis turbine, as in the reference case of this study, the vaneless space between the guide vane and runner blade is small. Pressure fluctuations associated to RSI are a major concern in such an area, as several failures related to hydraulic turbines are related to this phenomenon [39]. Figure 15 shows the Fourier-transformed pressure pulsation for three operating conditions at vaneless space with a sampling frequency of 4500 Hz. The research case presented here is comprised of 16 guide vane blades. For the runner rotation of 750 rpm, the value of fgv for the runner is 200 Hz. The first peak observed in Figure 15 represents the guide vane passing frequency. Consequently, other successive harmonics occurs at 2fb = 400 Hz, 3fb = 600 Hz, and so on. The pulsating signal is seen in every point for each operating conditions. It shows that pressure pulsating amplitude is most significant while operating the runner in part load conditions. In a high-head Francis turbine, as in the reference case of this study, the vaneless space between the guide vane and runner blade is small. Pressure fluctuations associated to RSI are a major concern in such an area, as several failures related to hydraulic turbines are related to this phenomenon [39]. Figure 15 shows the Fourier-transformed pressure pulsation for three operating conditions at vaneless space with a sampling frequency of 4500 Hz. The research case presented here is comprised of 16 guide vane blades. For the runner rotation of 750 rpm, the value of f gv for the runner is 200 Hz. The first peak observed in Figure 15 represents the guide vane passing frequency. Consequently, other successive harmonics occurs at 2f b = 400 Hz, 3f b = 600 Hz, and so on. The pulsating signal is seen in every point for each operating conditions. It shows that pressure pulsating amplitude is most significant while operating the runner in part load conditions. space between the guide vane and runner blade is small. Pressure fluctuations associated to RSI are a major concern in such an area, as several failures related to hydraulic turbines are related to this phenomenon [39]. Figure 15 shows the Fourier-transformed pressure pulsation for three operating conditions at vaneless space with a sampling frequency of 4500 Hz. The research case presented here is comprised of 16 guide vane blades. For the runner rotation of 750 rpm, the value of fgv for the runner is 200 Hz. The first peak observed in Figure 15 represents the guide vane passing frequency. Consequently, other successive harmonics occurs at 2fb = 400 Hz, 3fb = 600 Hz, and so on. The pulsating signal is seen in every point for each operating conditions. It shows that pressure pulsating amplitude is most significant while operating the runner in part load conditions.  Figure 16 shows FFT analysis of pressure pulsation inside the draft tube. As demonstrated in the figure, peak pressure pulsation is at 4 Hz for all measured points which is   FFT plot of pressure fluctuation from points RV1 to RV6 along the runner blade is presented in Figure 17. It is assumed that the points located in the runner surface rotate with the runner, hence the value obtained is during runner operation. It is evident from the figure that the periodic behavior of pressure fluctuation is governed by both frequencies. The low frequencies of 0.32, 0.64, and 0.96 times the runner rotational frequency correspond to the rotating vortex rope (RVR) frequency and its harmonics [41]. The bigger frequencies of 200 Hz, 400 Hz, and 600 Hz correspond to the guide vane passing frequency and its harmonics.
Frequency spectrum of pressure time signals showed that operation outside the BEP can result in high amplitudes and vortex roping, which can cause serious damage to the runner and other components. A swirling component will occur in the draft tube outside BEP operation. The high velocities in vortex core can decrease the pressure, resulting in a FFT plot of pressure fluctuation from points RV1 to RV6 along the runner blade is presented in Figure 17. It is assumed that the points located in the runner surface rotate with the runner, hence the value obtained is during runner operation. It is evident from the figure that the periodic behavior of pressure fluctuation is governed by both frequencies.
The low frequencies of 0.32, 0.64, and 0.96 times the runner rotational frequency correspond to the rotating vortex rope (RVR) frequency and its harmonics [41]. The bigger frequencies of 200 Hz, 400 Hz, and 600 Hz correspond to the guide vane passing frequency and its harmonics.

Torque Oscillations
RSI between guide vanes and runner induces pressure fluctuations inside the runner, which causes torque fluctuation. Fluctuation in torque inside the runner is caused by the expanding wakes from the guide vanes, which progresses until runner blades [42]. Hydraulic torque (Th) from the water acts on the runner, thus accelerating it. Th is a function of hydraulic power (Ph) and angular velocity (ω), as expressed in Equation (14). Figure 18 shows the dynamic forces in runner blade showing an average torque for three operating conditions plotted against runner angular position.

=
[Nm]  Frequency spectrum of pressure time signals showed that operation outside the BEP can result in high amplitudes and vortex roping, which can cause serious damage to the runner and other components. A swirling component will occur in the draft tube outside BEP operation. The high velocities in vortex core can decrease the pressure, resulting in a vapor-filled cavity core. This can cause operational challenges when coinciding with dynamics oscillations in a system.

Torque Oscillations
RSI between guide vanes and runner induces pressure fluctuations inside the runner, which causes torque fluctuation. Fluctuation in torque inside the runner is caused by the expanding wakes from the guide vanes, which progresses until runner blades [42]. Hydraulic torque (T h ) from the water acts on the runner, thus accelerating it. T h is a function of hydraulic power (P h ) and angular velocity (ω), as expressed in Equation (14). Figure 18 shows the dynamic forces in runner blade showing an average torque for three operating conditions plotted against runner angular position.
Energies 2021, 14, x FOR PEER REVIEW 17 of 21 Figure 17. Frequency spectrum of pressure time signals along runner blade.

Torque Oscillations
RSI between guide vanes and runner induces pressure fluctuations inside the runner, which causes torque fluctuation. Fluctuation in torque inside the runner is caused by the expanding wakes from the guide vanes, which progresses until runner blades [42]. Hydraulic torque (Th) from the water acts on the runner, thus accelerating it. Th is a function of hydraulic power (Ph) and angular velocity (ω), as expressed in Equation (14). Figure 18 shows the dynamic forces in runner blade showing an average torque for three operating conditions plotted against runner angular position.

=
[Nm]  Torque is normalized with the value corresponding for the BEP conditions. The patterns of torque fluctuations across different operating conditions are quite similar but the amplitude varies for all. The highest fluctuation of torque is observed at part load condition of operating stage shown by red graph in Figure 18, which is followed by the full load condition shown by the blue graph. The lowest torque fluctuation of all can be seen in BEP condition as shown by the green graph. The larger fluctuations were obtained for the full runner as compared to single blade for all cases, which agrees with the work done by Nicolle et al. [43]. For the single blade, a total of 16 peaks can be seen which is due to RSI arising from pressure fields around the guide vane and it corresponds to guide vanes number; however, for the full runner, 13 peaks are seen for all cases. The reason behind this is that the number of runner blades is a prime number and has no common integer with the number of guide vanes. As such, the torque variation for the runner must be completed in one pitch, i.e., 360/13. This phenomenon can also be explained with the assumption that when we consider full runner, more than one blade is experiencing RSI at the same time.

Conclusions
This work features an exploration of the spatial temporal progression of a leakage vortex in a high-head Francis turbine. The LV evolution and trajectory, LV progression from the clearance gap to runner blades, vortex dynamics, and pressure fluctuations have been analyzed in this study based on numerical simulations.
LV pattern and trajectory data were acquired from transient simulations. Vortex core intensity in the leakage flow decreased along the progression. The progression of leakage flow at the hub area was classified into three forms, i.e., the elongation stage, disintegration stage, and dissolving stage. Leakage flow at the shroud area is transported until the midsection of the blade, i.e., before the majority of vortex is dissolved, and seems to pass from the outlet of the runner as well. The values of vorticity and TKE along the LV trajectory somehow gave a clear picture of the complex vortices and turbulence phenomena for varying operating cases.
Unsteady pressure fluctuations were investigated with the reference turbine at three operating conditions, i.e., BEP, partial loading, and full loading. A total of sixteen points were located in the vaneless space, runner blade, and draft tube to study the pressure pulsations in the stationary and rotating domains. The frequency spectra data of pressure signals in the vaneless space showed the highest amplitude at the guide vane passing frequency (fgv) and its successive harmonics, i.e., 2fgv, 3fgv, 4fgv and so on, due to RSI. Monitoring the points at the draft tube showed peak pressure pulsation that was 0.32 times that of the runner rotational frequency and subsequent peaks for its harmonics. The points located along the runner blade from the leading edge to the trailing edge showed that the frequency spectrum was affected by both available frequencies. Low frequencies of 0.32, 0.64, and 0.96 times the runner rotational frequency corresponded to the rotating vortex rope frequency and its harmonics, whereas higher frequencies of 200 Hz, 400 Hz, and 600 Hz corresponded to the guide vane passing frequency and its harmonics. Moreover, while taking torque fluctuations for different operating conditions into account, it was observed that the patterns were quite similar in terms of the variation of amplitude. The highest fluctuation of the torque was observed at partial load condition, followed by the full load condition, and least for the BEP condition.
In summary, important dynamic characteristics for leakage vortex evolution and progression have been revealed in this work. The frequency spectra of pressure signals at various locations in the turbine provide some information regarding pressure fields in the blade channel; however, the spatial temporal evolution of a vortex is complex phenomenon which is governed by several factors. Further investigation is needed for overall energy improvement and the elucidation of a novel technique for vortex suppression before hitting a runner. Grooves or shapes can be added and studied which could reduce the flow separation and pressure gradient between two sides of the CG and hence suppress the