Investigation of Pressure Fluctuation and Pulsating Hydraulic Axial Thrust in Francis Turbines

: Under the circumstances of rapid expansion of diverse forms of volatile and intermittent renewable energy sources, hydropower stations have become increasingly indispensable for improving the quality of energy conversion processes. As a consequence, Francis turbines, one of the most popular options, need to operate under o ﬀ -design conditions, particularly for partial load operation. In this paper, a prototype Francis turbine was used to investigate the pressure ﬂuctuations and hydraulic axial thrust pulsation under four partial load conditions. The analyses of pressure ﬂuctuations in the vaneless space, runner, and draft tube are discussed in detail. The observed precession frequency of the vortex rope is 0.24 times that of the runner rotational frequency, which is able to travel upstream (from the draft tube to the vaneless space). Frequencies of both 24.0 and 15.0 times that of the runner rotational frequency are detected in the recording points of the runner surface, while the main dominant frequency recorded in the vaneless zone is 15.0 times that of the runner rotational frequency. Apart from unsteady pressure ﬂuctuations, the pulsating property of hydraulic axial thrust is discussed in depth. In conclusion, the pulsation of hydraulic axial thrust is derived from the pressure ﬂuctuations of the runner surface and is more complicated than the pressure ﬂuctuations.


Introduction
Presently, renewable energy sources have increasingly contributed to global electricity generation [1]. As a result, the dependence upon fossil fuel is gradually decreasing [2]. In this context, environmental pollution caused by the combustion of fossil fuels is being alleviated accordingly. In the course of developing various renewable and new resources with a long history of exploitation, hydropower holds the largest share and has become one of the most suitable and cost-effective solution for energy supply in both urban and rural areas [2,3]. On a global scale, the hydropower-installed capacity is still growing, and related conservancy projects are expected to prosper, especially in East Asia and the Pacific, South America, and South and Central Asia [1,2].
Hydraulic turbines, one of the critical components of a hydropower plant, are responsible for converting hydro energy into mechanical energy; subsequently, electricity can be generated via a generator. To date, a miscellaneous assortment of hydraulic turbines has been designed to operate reliably and efficiently under different conditions of head and discharge at the actual site [4,5]. Due to desirable advantages, such as a wide operating range and high efficiency, Francis turbines are widely used compared to other types of hydraulic turbines [6,7]. With the advantage being able to regulate operating conditions for a rather low cost and short time, Francis turbines are able to meet real-time electricity demands and, thus, are frequently employed to stabilize power grid operation [6,8,9]. However, Francis turbines are generally designed to operate at the best efficiency points (BEPs) and, theoretically, the flow leaving the runner is almost axial, with little swirl entering the draft tube. Francis turbines are inevitably operated from a partial to full load for the sake of accommodating the real-time demand of end-users, which is likely to bring about unfavorable issues that undermine hydraulic stability, including Karman vortex at the trailing edge of guide vanes, blade channel vortices, and helical vortex rope in the draft tube [10,11]. In addition, being one of the crucial factors in hydraulic stability, hydraulic axial thrust represents severe fluctuations and poses a great threat to the bearing and powerhouse structure under off-design conditions.
Hydraulic axial thrust is mainly comprised of the axial force on the surface of runner blades and the inner and outer surface of the crown and band. Numerous studies can be found regarding theoretical calculation of the hydraulic axial thrust for the design of the thrust bearing and powerhouse structure [12]. In these studies, a constant was determined for specific operating points, and this has satisfied the design criteria. However, the theoretical equations fail to take fluctuation of hydraulic axial thrust into account, which can no longer be ignored given the increasing attention toward the vibration problem of Francis turbines. To achieve this, both laboratory measurements and computational fluid dynamics (CFD) can be used to obtain the evolution of hydraulic axial thrust over time. The latter is favored by researchers for its benefits of robustness, convenience, and capturing of detailed fluid information. To date, CFD methodology has been adopted by researchers and engineers not only to elaborate such complex dynamic fluid flow effects, but also to carry out the optimization of turbine components during design or rehabilitation [13,14].
Over the last several decades, both experiments and numerical simulations have been performed to investigate flow instabilities as well as the induced pressure fluctuations inside Francis turbines [15,16]. Arpe et al. [17] presented a detailed experimental study about the unstable pressure field that developed under low-discharge operating conditions within the elbow-type draft tube of a high specific speed Francis turbine. They found that pressure fluctuations resulted from the vortex precession within the range of 0.2-0.4 f n (runner rotational frequency), which were generated only in the conical part of the draft tube and then disappeared in the straight diffuser. Müller et al. [18] showed the high-speed visualization of the runner blade channel flow pattern during one pressure oscillating period and concluded that the critical cavitation evolution on the runner blades modified their hydrodynamic performance to some degree, so as to result in a reduction in the momentum transferred from flow to the turbine shaft. Frunzverde et al. [19] illustrated a case of runner failure in a low head Francis turbine, where a crack was discovered at the junction between the blade leading edge and the crown. Their discussion revealed that both the rotor stator interaction (RSI) and the vortex rope were capable of giving rise to dynamic stresses in the runner blades. Trivedi and Cervantes [20] reviewed the studies conducted on RSI within hydraulic turbines and pointed out that the frequencies in the stationary and rotating domains were computed using f r = n·Z gv and f s = n·Z b , respectively. In the two equations, n referred to the runner speed in revolutions per second, and Z gv and Z b were the number of guide vanes and rotating blades, respectively. Xia et al. [21] investigated the pressure fluctuations and runner loads on a Francis pump-turbine runner during the runaway process, and found that the variation of discharge was the main reason for the intense fluctuation in hydraulic thrust during this process. However, limited studies are associated with the pulsating property of hydraulic axial thrust under the impact of multiple hydraulic excitation forces. Being one of the crucial factors in hydraulic stability, hydraulic axial thrust is the resultant force in the axial direction of the turbine runner exerted by the flow in the turbine passage, which is then transmitted to the load-bearing frame base of the powerhouse via thrust bearing, and consequently, poses a threat to the bearing and the powerhouse. As Francis turbines are more liable to operating under off-design conditions than before, it is imperative to have a deeper understanding of hydraulic thrust for the safe and stable operation of hydropower stations. To bridge the gap that exists in the available state-of-the-art studies on hydraulic axial thrust, the aim of this paper was to numerically analyze the pulsating property of hydraulic axial thrust under the impact of multiple hydraulic excitation forces. Given that the size of the labyrinth seal gap is much smaller than that of the main turbine passage, it is hard to take the labyrinth seal into consideration in the whole turbine model simulation. In addition, a small variation in seal clearance would have a significant influence over the thrust value and, therefore, only axial force on the surface of the runner blades and inner surface of the crown and band were investigated in this paper. After conducting the performance evaluation, the Francis turbine operated at four different guide vane openings (GVOs) was explored in order to fully study the pressure fluctuations in the draft tube, runner, and vaneless space, as well as to clarify the fluctuation mechanism of hydraulic axial thrust.

Computational Domain
The Francis turbine investigated in this study is a prototype that operates in a hydropower plant in China. Its design head is 116.0 m, the rated volume flow rate is 87.72 m 3 /s, the rated power output is 90 MW, and the designed rotational speed is 250 r/min. Figure 1a shows the whole turbine geometry in the following simulations. Figure 1b-e presents the discretization of a spiral case, a tandem cascade, a runner blade, and a draft tube, respectively. The computational domain includes five subdomains, i.e., a spiral case, a region with 24 stay vanes, a region with 24 guide vanes, a region with 15 runner blades, and a draft tube. Two kinds of mesh were created for these five domains via TurboGrid 14.5 and ANSYS ICEM 14.5 to produce high-quality grids. The spiral case, stay and guide vanes, and draft tube were discretized using hexahedral grids, while tetrahedral grids were generated in the runner. As for the number of grids, eight different sizes of grids with element numbers varying from 0.98 × 10 6 to 8.32 × 10 6 were created at 18 deg guide vane opening to conduct the mesh sensitivity analysis (see Figure 2). After confirmation of the variation of average volume flow rate and runner torque with the grid number and taking into account the trade-off between the required computational resources and the calculation accuracy, a total number of 4.12 × 10 6 grids was used for the subsequent investigation (see Table 1). To bridge the gap that exists in the available state-of-the-art studies on hydraulic axial thrust, the aim of this paper was to numerically analyze the pulsating property of hydraulic axial thrust under the impact of multiple hydraulic excitation forces. Given that the size of the labyrinth seal gap is much smaller than that of the main turbine passage, it is hard to take the labyrinth seal into consideration in the whole turbine model simulation. In addition, a small variation in seal clearance would have a significant influence over the thrust value and, therefore, only axial force on the surface of the runner blades and inner surface of the crown and band were investigated in this paper. After conducting the performance evaluation, the Francis turbine operated at four different guide vane openings (GVOs) was explored in order to fully study the pressure fluctuations in the draft tube, runner, and vaneless space, as well as to clarify the fluctuation mechanism of hydraulic axial thrust.

Computational Domain
The Francis turbine investigated in this study is a prototype that operates in a hydropower plant in China. Its design head is 116.0 m, the rated volume flow rate is 87.72 m 3 /s, the rated power output is 90 MW, and the designed rotational speed is 250 r/min. Figure 1a shows the whole turbine geometry in the following simulations. Figure 1b-e presents the discretization of a spiral case, a tandem cascade, a runner blade, and a draft tube, respectively. The computational domain includes five subdomains, i.e., a spiral case, a region with 24 stay vanes, a region with 24 guide vanes, a region with 15 runner blades, and a draft tube. Two kinds of mesh were created for these five domains via TurboGrid 14.5 and ANSYS ICEM 14.5 to produce high-quality grids. The spiral case, stay and guide vanes, and draft tube were discretized using hexahedral grids, while tetrahedral grids were generated in the runner. As for the number of grids, eight different sizes of grids with element numbers varying from 0.98 × 10 6 to 8.32 × 10 6 were created at 18 deg guide vane opening to conduct the mesh sensitivity analysis (see Figure 2). After confirmation of the variation of average volume flow rate and runner torque with the grid number and taking into account the trade-off between the required computational resources and the calculation accuracy, a total number of 4.12 × 10 6 grids was used for the subsequent investigation (see Table 1).    Table 2. Most of the simulation processes and setup were consistent with the author's previous paper [8]. However, a smaller time step of 0.0005 s corresponding to 0.75 degree of runner rotation was chosen in this paper to better capture the high-frequency dynamic effects inside the whole flow passage [22,23]; the total time of calculation was set to 10 s (about 40 runner revolutions). The flow rates of four different partial load conditions were achieved by changing the guide vane opening while ensuring the same inlet and outlet pressure to maintain the water head design. Due to the inconsistency of grids between adjacent subdomains, four interfaces were created for transferring computational information, namely between the spiral case and the stay vane, between the stay vane and the guide vane, between the guide vane and the runner, and between the runner and the draft tube.

Recording Points and Sections
Fourteen points located in the vaneless space, the runner, and the draft tube (see red crosses in Figure 3) were selected to record the time-resolved pressure. Three monitoring points of R1p, R2p, and R3p were located on the pressure side of one blade from near the runner inlet to near the runner outlet, while three others, R1n, R2n, and R3n, were positioned in the corresponding locations on the suction side of the same blade. Four monitoring points, defined as H1 and H2, and S1 and S2, were located on the runner hub and shroud, respectively. For the recording points belonging to stationary walls, they were created by specifying each individual Cartesian coordinate. Meanwhile, for the monitoring points attached to the rotating walls (such as runner blades), there were two steps: firstly, an intersecting line between the blade vanes and a cylindrical surface was drawn; then, the corresponding points can be created by seeking the intersection of this intersecting line with a z-axis   Table 2. Most of the simulation processes and setup were consistent with the author's previous paper [8]. However, a smaller time step of 0.0005 s corresponding to 0.75 degree of runner rotation was chosen in this paper to better capture the high-frequency dynamic effects inside the whole flow passage [22,23]; the total time of calculation was set to 10 s (about 40 runner revolutions). The flow rates of four different partial load conditions were achieved by changing the guide vane opening while ensuring the same inlet and outlet pressure to maintain the water head design. Due to the inconsistency of grids between adjacent subdomains, four interfaces were created for transferring computational information, namely between the spiral case and the stay vane, between the stay vane and the guide vane, between the guide vane and the runner, and between the runner and the draft tube.

Recording Points and Sections
Fourteen points located in the vaneless space, the runner, and the draft tube (see red crosses in Figure 3) were selected to record the time-resolved pressure. Three monitoring points of R 1p , R 2p , and R 3p were located on the pressure side of one blade from near the runner inlet to near the runner outlet, while three others, R 1n , R 2n , and R 3n , were positioned in the corresponding locations on the suction side of the same blade. Four monitoring points, defined as H 1 and H 2 , and S 1 and S 2 , were located on the runner hub and shroud, respectively. For the recording points belonging to stationary walls, they were created by specifying each individual Cartesian coordinate. Meanwhile, for the monitoring points attached to the rotating walls (such as runner blades), there were two steps: firstly, an intersecting line between the blade vanes and a cylindrical surface was drawn; then, the corresponding points can be created by seeking the intersection of this intersecting line with a z-axis plane. Five sections, ascending from P 1 to P 5 (indicated in Figure 3), were chosen to analyze the pressure distribution in the draft tube passage. To investigate hydraulic axial thrust, the monitoring was conducted by plane. Five sections, ascending from P1 to P5 (indicated in Figure 3), were chosen to analyze the pressure distribution in the draft tube passage. To investigate hydraulic axial thrust, the monitoring was conducted by recording the time-resolved net lift force that acts on runner surfaces, including the runner blade, hub, and shroud.

Performance Evaluation and Assessment of Unsteady 3D Flow Simulations
To confirm the quality of the numerical methodology of the unsteady three-dimensional flow simulations performed in this paper and to assess the performance of this Francis turbine, four partial load operating points (including 13 deg (O1), 16 deg (O2), 18 deg (O3), and 19 deg (O4) GVOs) were selected for analyzing the oscillating properties in terms of runner torque, hydraulic efficiency, vortex rope, and other parameters. It should be noted that the corresponding average volume flow rates of these four working conditions are 44.58 (O1), 56.38 (O2), 65.71 (O3) and 70.13 (O4) m 3 /s.
Considering that there exists many dynamical effects (including RSI, vortex rope, cavitation, etc.) inside a large-sized Francis turbine, it is more appropriate to adopt unsteady calculations of flow, which is a more accurate way to capture the cyclic interactions [24]. Figure 4a depicts the variation of runner torque with time, and this plot gives a general idea about the characteristic of power output due to periodic flow instabilities inside the flow passage. The runner torque presents two types of periodic behavior: one is local high-frequency oscillation; the other is long-distance low-frequency oscillation. The dominant relative frequencies ⁄ were calculated via a fast Fourier transform (FFT), which was also used in the following analysis of pressure fluctuation and hydraulic axial thrust. As revealed by the FFT spectrum (Figure 4b), five distinct peaks could be observed that are 0.24, 4.0, 15.0, 24.0, and 30.0 times that of the runner rotational frequency. As for the reason for the appearance of 24.0 times that of the runner revolution, this corresponds to the number of guide vanes during one revolution of the runner, which was reported, in detail, in papers [24,25]. However, the other four dominant frequencies are rarely observed. Given that dynamic asymmetrical pressure distributions on runner blades may result in unstable runner torque, the sources of these four frequencies are explored from the perspective of pressure fluctuation in the subsequent sections. Above results help prove the simulations of unsteady 3D Flow are reliable to capture the expected high frequencies (such as blade/guide vane passing frequencies) and their harmonics so as to obtain fine time resolution for the dynamic analysis. It should be noted that only the results for 16 deg GVO case are presented in some figures (including following sections) as it is a typical partial load condition wherein vortex rope, pressure fluctuations and other instability pattern are very prominent and distinct.

Performance Evaluation and Assessment of Unsteady 3D Flow Simulations
To confirm the quality of the numerical methodology of the unsteady three-dimensional flow simulations performed in this paper and to assess the performance of this Francis turbine, four partial load operating points Considering that there exists many dynamical effects (including RSI, vortex rope, cavitation, etc.) inside a large-sized Francis turbine, it is more appropriate to adopt unsteady calculations of flow, which is a more accurate way to capture the cyclic interactions [24]. Figure 4a depicts the variation of runner torque with time, and this plot gives a general idea about the characteristic of power output due to periodic flow instabilities inside the flow passage. The runner torque presents two types of periodic behavior: one is local high-frequency oscillation; the other is long-distance low-frequency oscillation. The dominant relative frequencies ( f / f n ) were calculated via a fast Fourier transform (FFT), which was also used in the following analysis of pressure fluctuation and hydraulic axial thrust. As revealed by the FFT spectrum (Figure 4b), five distinct peaks could be observed that are 0.24, 4.0, 15.0, 24.0, and 30.0 times that of the runner rotational frequency. As for the reason for the appearance of 24.0 times that of the runner revolution, this corresponds to the number of guide vanes during one revolution of the runner, which was reported, in detail, in papers [24,25]. However, the other four dominant frequencies are rarely observed. Given that dynamic asymmetrical pressure distributions on runner blades may result in unstable runner torque, the sources of these four frequencies are explored from the perspective of pressure fluctuation in the subsequent sections. Above results help prove the simulations of unsteady 3D Flow are reliable to capture the expected high frequencies (such as blade/guide vane passing frequencies) and their harmonics so as to obtain fine time resolution for the dynamic analysis. It should be noted that only the results for 16 deg GVO case are presented in some figures (including following sections) as it is a typical partial load condition wherein vortex rope, pressure fluctuations and other instability pattern are very prominent and distinct.  After obtaining runner torque, the hydraulic efficiency, h  , of the turbine can be acquired with where N is the runner output,  is the water density, g is the gravity constant, Q is the volume flow rate, and H is the net pressure head.
The runner output can be obtained with Equation (2).
where M is the runner torque and w is the angular velocity of the runner.
The total pressure head, H, is the pressure difference between the spiral case inlet and the draft tube outlet. According to the Bernoulli equation, H can be acquired with Equation (3).
where in P and out P are the pressure at the inlet of the spiral case and the outlet of the draft tube, respectively; in V and out V are the velocity at the inlet of the spiral case and the outlet of the draft tube, respectively; ∆z denotes the vertical distance between centerline of the inlet and the outlet of the turbine. In this study, ∆z = 6.86 m for the Francis turbine was investigated.
The hydraulic efficiency and the runner output of the turbine can be calculated in accordance with Equations (1)-(3), as shown in Figure 5. The aforementioned analysis indicates that hydraulic efficiency and runner output decrease as the operating condition deviates from the BEP. Please note that for the O1 operating point, the hydraulic efficiency drops dramatically, therefore Francis turbine should be operated within a reasonable operating range, so as to avoid operation which is too far away BEP and corresponding hydraulic efficiency is able to maintain at a desirable level. On the other hand, the output of Francis turbine should satisfy the demand of end-users, some measures should be taken to promote the consumption of excess power produced by Francis turbine. For example, support the construction of transmission line to large factories that can function with varying power input (e.g., aluminum smelter); enhance the planning and construction of power transmission channels for the mutual benefit of hydropower and coal-fired power; adjust the pricing strategy for the consumption of excess power and increase local absorption and electricity delivering. After obtaining runner torque, the hydraulic efficiency, η h , of the turbine can be acquired with Equation (1).
where N is the runner output, ρ is the water density, g is the gravity constant, Q is the volume flow rate, and H is the net pressure head. The runner output can be obtained with Equation (2).
where M is the runner torque and w is the angular velocity of the runner. The total pressure head, H, is the pressure difference between the spiral case inlet and the draft tube outlet. According to the Bernoulli equation, H can be acquired with Equation (3).
where P in and P out are the pressure at the inlet of the spiral case and the outlet of the draft tube, respectively; V in and V out are the velocity at the inlet of the spiral case and the outlet of the draft tube, respectively; ∆z denotes the vertical distance between centerline of the inlet and the outlet of the turbine. In this study, ∆z = 6.86 m for the Francis turbine was investigated. The hydraulic efficiency and the runner output of the turbine can be calculated in accordance with Equations (1)-(3), as shown in Figure 5. The aforementioned analysis indicates that hydraulic efficiency and runner output decrease as the operating condition deviates from the BEP. Please note that for the O1 operating point, the hydraulic efficiency drops dramatically, therefore Francis turbine should be operated within a reasonable operating range, so as to avoid operation which is too far away BEP and corresponding hydraulic efficiency is able to maintain at a desirable level. On the other hand, the output of Francis turbine should satisfy the demand of end-users, some measures should be taken to promote the consumption of excess power produced by Francis turbine. For example, support the construction of transmission line to large factories that can function with varying power input (e.g., aluminum smelter); enhance the planning and construction of power transmission channels for the mutual benefit of hydropower and coal-fired power; adjust the pricing strategy for the consumption of excess power and increase local absorption and electricity delivering.

In Draft Tube
The flow instabilities in a Francis turbine draft tube under partial load conditions are closely related to the formation of a helical processing vortex, also known as the vortex rope. The iso-surfaces of swirling strength, λ3D = 0.01 s −1 , were adopted to illustrate the evolution of the vortex structure in this paper (see help of CFD-Post for derivation of λ3D). The spiral vortex ropes generated under the four different GVOs are presented in Figure 6. A distinct vortex rope is observed for operating conditions of 16, 18, and 19 deg GVOs; it could be deduced that the volume of the vortex rope becomes smaller when the operating point is closer to the BEP. For the operating condition of 13 deg GVO, the vortex shedding is so intense that no major vortex rope is formed, instead resulting in a fragmented vortex core. The rotating vortex rope induces an uneven pressure distribution in the same cross section, as displayed in Figure 7. There exists a negative pressure zone in the central region as a result of low rotational velocity, while the peripheral rotational velocity is comparatively large. The pressure increases gradually along the radial direction, and the radial distribution of pressure is relatively uniform within the draft tube cone. The above analysis indicates that vortex rope and pressure distribution are caused by the large tangential component of the velocity at the runner outlet under the partial load condition. Therefore, countermeasures for alleviating the pressure fluctuations within draft tube could be proposed by reducing the tangential velocity of the flow in the draft tube inlet.

In Draft Tube
The flow instabilities in a Francis turbine draft tube under partial load conditions are closely related to the formation of a helical processing vortex, also known as the vortex rope. The iso-surfaces of swirling strength, λ 3D = 0.01 s −1 , were adopted to illustrate the evolution of the vortex structure in this paper (see help of CFD-Post for derivation of λ 3D ). The spiral vortex ropes generated under the four different GVOs are presented in Figure 6. A distinct vortex rope is observed for operating conditions of 16, 18, and 19 deg GVOs; it could be deduced that the volume of the vortex rope becomes smaller when the operating point is closer to the BEP. For the operating condition of 13 deg GVO, the vortex shedding is so intense that no major vortex rope is formed, instead resulting in a fragmented vortex core. The rotating vortex rope induces an uneven pressure distribution in the same cross section, as displayed in Figure 7. There exists a negative pressure zone in the central region as a result of low rotational velocity, while the peripheral rotational velocity is comparatively large. The pressure increases gradually along the radial direction, and the radial distribution of pressure is relatively uniform within the draft tube cone. The above analysis indicates that vortex rope and pressure distribution are caused by the large tangential component of the velocity at the runner outlet under the partial load condition. Therefore, countermeasures for alleviating the pressure fluctuations within draft tube could be proposed by reducing the tangential velocity of the flow in the draft tube inlet.
During the partial load operation of a Francis turbine, the violent pressure fluctuations are associated with the formation of the vortex rope in the draft tube. Figure 8a presents the variation in wall pressure (normalized by subtracting the average pressure and same practice was applied to following time domain figures) with time at point D 1 located 2.35 m downstream of the draft tube inlet. As can be seen from Figure 8a, there exists two kinds of periodic behaviors which can be separated by FFT analysis. As demonstrated in Figure 8b, the main pressure fluctuation and its harmonics brought about by the vortex rope exhibit low-frequency and high-amplitude characteristics. The precession frequency is 0.24 times that of the runner rotational frequency, which is consistent with the value of Energies 2020, 13, 1734 8 of 16 0.2-0.4 times that of the runner rotational frequency presented in earlier studies [17,22,26,27]. It should be noted that a frequency of 15.0 times that of the runner rotational frequency is also found, which corresponds to the number of runner blades and signals that the pressure wave is able to travel from the runner domain to the draft tube wall in the whole hydraulic circuit. This amplitude of RSI frequency (f runner ) was also confirmed in the experimental and numerical studies performed by Trivedi et al. [28,29].  During the partial load operation of a Francis turbine, the violent pressure fluctuations are associated with the formation of the vortex rope in the draft tube. Figure 8a presents the variation in wall pressure (normalized by subtracting the average pressure and same practice was applied to following time domain figures) with time at point D1 located 2.35 m downstream of the draft tube inlet. As can be seen from Figure 8a, there exists two kinds of periodic behaviors which can be separated by FFT analysis. As demonstrated in Figure 8b, the main pressure fluctuation and its harmonics brought about by the vortex rope exhibit low-frequency and high-amplitude characteristics. The precession frequency is 0.24 times that of the runner rotational frequency, which is consistent with the value of 0.2-0.4 times that of the runner rotational frequency presented in earlier studies [17,22,26,27]. It should be noted that a frequency of 15.0 times that of the runner rotational frequency is also found, which corresponds to the number of runner blades and signals that the  During the partial load operation of a Francis turbine, the violent pressure fluctuations are associated with the formation of the vortex rope in the draft tube. Figure 8a presents the variation in wall pressure (normalized by subtracting the average pressure and same practice was applied to following time domain figures) with time at point D1 located 2.35 m downstream of the draft tube inlet. As can be seen from Figure 8a, there exists two kinds of periodic behaviors which can be separated by FFT analysis. As demonstrated in Figure 8b, the main pressure fluctuation and its harmonics brought about by the vortex rope exhibit low-frequency and high-amplitude characteristics. The precession frequency is 0.24 times that of the runner rotational frequency, which is consistent with the value of 0.2-0.4 times that of the runner rotational frequency presented in earlier studies [17,22,26,27]. It should be noted that a frequency of 15.0 times that of the runner rotational frequency is also found, which corresponds to the number of runner blades and signals that the pressure wave is able to travel from the runner domain to the draft tube wall in the whole hydraulic  Figure 9 shows the normalized frequency spectrum obtained from D 1 to D 6 for 16 deg GVO. The pressure amplitude in each individual part of the draft tube (consisting of cone, elbow part, and diffuser) indicates that the amplitude in the draft tube cone is the largest; then, the pressure oscillation gradually weakens downstream. This phenomenon can be explained by the tangential velocity being largest in the cone area and decreasing downstream on account of viscous effects. Figure 10 presents the normalized frequency spectrum of D 1 under the four different operating conditions. It can be concluded that the pressure fluctuation of the draft tube reaches its maximum Energies 2020, 13, 1734 9 of 16 in the operating condition of 16 deg GVO; the corresponding amplitude decreases with decreasing GVO. For 13 deg GVO, its small amplitude of pressure fluctuation is a result of the non-existence of a major vortex rope. These results are closely related to the size of the vortex rope and can be clarified by Figure 6.
Energies 2020, 13, x FOR PEER REVIEW 9 of 16 Figure 9 shows the normalized frequency spectrum obtained from D1 to D6 for 16 deg GVO. The pressure amplitude in each individual part of the draft tube (consisting of cone, elbow part, and diffuser) indicates that the amplitude in the draft tube cone is the largest; then, the pressure oscillation gradually weakens downstream. This phenomenon can be explained by the tangential velocity being largest in the cone area and decreasing downstream on account of viscous effects. Figure 10 presents the normalized frequency spectrum of D1 under the four different operating conditions. It can be concluded that the pressure fluctuation of the draft tube reaches its maximum in the operating condition of 16 deg GVO; the corresponding amplitude decreases with decreasing GVO. For 13 deg GVO, its small amplitude of pressure fluctuation is a result of the non-existence of a major vortex rope. These results are closely related to the size of the vortex rope and can be clarified by Figure 6.    Energies 2020, 13, x FOR PEER REVIEW 9 of 16 Figure 9 shows the normalized frequency spectrum obtained from D1 to D6 for 16 deg GVO. The pressure amplitude in each individual part of the draft tube (consisting of cone, elbow part, and diffuser) indicates that the amplitude in the draft tube cone is the largest; then, the pressure oscillation gradually weakens downstream. This phenomenon can be explained by the tangential velocity being largest in the cone area and decreasing downstream on account of viscous effects. Figure 10 presents the normalized frequency spectrum of D1 under the four different operating conditions. It can be concluded that the pressure fluctuation of the draft tube reaches its maximum in the operating condition of 16 deg GVO; the corresponding amplitude decreases with decreasing GVO. For 13 deg GVO, its small amplitude of pressure fluctuation is a result of the non-existence of a major vortex rope. These results are closely related to the size of the vortex rope and can be clarified by Figure 6.     Figure 9 shows the normalized frequency spectrum obtained from D1 to D6 for 16 deg GVO. The pressure amplitude in each individual part of the draft tube (consisting of cone, elbow part, and diffuser) indicates that the amplitude in the draft tube cone is the largest; then, the pressure oscillation gradually weakens downstream. This phenomenon can be explained by the tangential velocity being largest in the cone area and decreasing downstream on account of viscous effects. Figure 10 presents the normalized frequency spectrum of D1 under the four different operating conditions. It can be concluded that the pressure fluctuation of the draft tube reaches its maximum in the operating condition of 16 deg GVO; the corresponding amplitude decreases with decreasing GVO. For 13 deg GVO, its small amplitude of pressure fluctuation is a result of the non-existence of a major vortex rope. These results are closely related to the size of the vortex rope and can be clarified by Figure 6.

In Runner
The recording points located in the runner surface rotate with the runner; in this way, the data are able to reflect the pressures on the blade during operation. The time domain plot of the pressure fluctuation of R 1p for the 16 deg GVO condition is presented in Figure 11a. It is evident that the periodic behavior of the pressure fluctuation is more complicated for the points belonging to the runner blade than those of the draft tube wall. After conducting FFT analysis, three kinds of dominant frequencies Energies 2020, 13, 1734 10 of 16 could be observed in Figure 11b. For the low frequencies of 0.24, 0.72, 0.96, and 1.92 times that of the runner rotational frequency, they correspond to the rotating vortex rope (RVR) frequency and its harmonics. It is worth noting that the highest amplitude corresponds to the fourth harmonic of the RVR frequency. The high frequency of 15.0 times that of the runner rotational frequency pertains to the blade passing frequency. This RSI frequency (f runner ) was also confirmed in the numerical simulations of Francis turbines performed by Zhu et al. [30]. However, it could not be found in the numerical results of Francis pump-turbines conducted by Trivedi [31]. The other high frequency of 24.0 times that of the runner rotational frequency belongs to the guide vane passing frequency. This particular RSI frequency (f gv ) was also detected by Trivedi et al. [28,32].
pertains to the blade passing frequency. This RSI frequency (frunner) was also confirmed in the numerical simulations of Francis turbines performed by Zhu et al. [30]. However, it could not be found in the numerical results of Francis pump-turbines conducted by Trivedi [31]. The other high frequency of 24.0 times that of the runner rotational frequency belongs to the guide vane passing frequency. This particular RSI frequency (fgv) was also detected by Trivedi et al. [28,32].
The FFT spectrums of the pressure-time signals for the six recording points on the runner blade are shown in Figure 12a. The plot indicates that the pressure amplitude on the suction side is higher than that on the pressure side, which agrees well with the findings in [33,34]; the pressure fluctuation amplitudes of the recording points near the runner outlet are much larger than those near the runner inlet. The frequency spectrum for the other four recording points of the runner hub and shroud is illustrated in Figure 12b. As shown in this picture, the highest amplitude of the two recording points near the runner inlet corresponds to the guide vane passing frequency, fgv; while the highest amplitude of the other two recording points near the runner outlet belong to the RVR frequency. Figure 13 displays the normalized frequency spectrum obtained from R1p under four different operating conditions. It can be deduced that the pressure fluctuation in the runner reaches its peak in the operating condition of 16 deg GVO. In addition, the other basic law of change related to the working condition is consistent with that of the monitoring point in the draft tube (see Figure 10).  The FFT spectrums of the pressure-time signals for the six recording points on the runner blade are shown in Figure 12a. The plot indicates that the pressure amplitude on the suction side is higher than that on the pressure side, which agrees well with the findings in [33,34]; the pressure fluctuation amplitudes of the recording points near the runner outlet are much larger than those near the runner inlet. The frequency spectrum for the other four recording points of the runner hub and shroud is illustrated in Figure 12b. As shown in this picture, the highest amplitude of the two recording points near the runner inlet corresponds to the guide vane passing frequency, f gv ; while the highest amplitude of the other two recording points near the runner outlet belong to the RVR frequency. Figure 13 displays the normalized frequency spectrum obtained from R 1p under four different operating conditions. It can be deduced that the pressure fluctuation in the runner reaches its peak in the operating condition of 16 deg GVO. In addition, the other basic law of change related to the working condition is consistent with that of the monitoring point in the draft tube (see Figure 10).

In Vaneless Space
The diagrams regarding the time and frequency domain of the pressure fluctuation for V1 are displayed in Figure 14, respectively. The highest amplitude corresponds to the runner blade passing frequency, and the second highest amplitude is the RVR frequency. These two frequencies were also observed in the calculation of Francis pump-turbines performed by Iliev et al. [35]; whereas, only 15.0 times that of the runner rotational frequency was reported in some of the papers by Trivedi et al. [36][37][38]. Based on the simulation results in this study, it can be concluded that the pressure wave of RVR is able to travel upstream (from the draft tube domain to the runner domain, and then to vaneless space), although the amplitude of fluctuation decreases, indicating that the oscillating intensity gradually recedes. Finally, the recording points in the vaneless space have the lowest peak pressure amplitude compared to those of the runner blade and the draft tube wall. Figure 15 shows the normalized frequency spectrum obtained from V1 under the four different operating conditions, which reveals that there is not much difference between the various operating points in terms of the highest amplitude of pressure fluctuation.

In Vaneless Space
The diagrams regarding the time and frequency domain of the pressure fluctuation for V 1 are displayed in Figure 14, respectively. The highest amplitude corresponds to the runner blade passing frequency, and the second highest amplitude is the RVR frequency. These two frequencies were also observed in the calculation of Francis pump-turbines performed by Iliev et al. [35]; whereas, only 15.0 times that of the runner rotational frequency was reported in some of the papers by Trivedi et al. [36][37][38]. Based on the simulation results in this study, it can be concluded that the pressure wave of RVR is able to travel upstream (from the draft tube domain to the runner domain, and then to vaneless space), although the amplitude of fluctuation decreases, indicating that the oscillating intensity gradually recedes. Finally, the recording points in the vaneless space have the lowest peak pressure amplitude compared to those of the runner blade and the draft tube wall. Figure 15 shows the normalized frequency spectrum obtained from V 1 under the four different operating conditions, which reveals that there is not much difference between the various operating points in terms of the highest amplitude of pressure fluctuation.

Hydraulic Axial Thrust
The pulsating hydraulic axial thrust is transmitted to the load-bearing frame base of the powerhouse via thrust bearing, and consequently poses a threat to the bearing and the powerhouse structure. During certain transition processes, accidental lifting of the rotating part will occur when the thrust bearing is subjected to upward force; on the other hand, excessive downward load on the thrust bearing may increase the temperature of the thrust bearing pads, thus causing damage to them [39]. In addition, the pulsation characteristics of axial hydraulic thrust make it one of the main sources of noise and strong vibration occurring in the powerhouse structure. There are many hydropower plants where prominent vibration problems frequently take place during off-design and transient operations, such as the Three Gorges Hydropower Station, the Ertan Hydropower Station, and the Guangzhou Pumped-Storage Power Station [40,41]; therefore, it is necessary to explore, in depth, and understand the pulsation properties of the hydraulic axial thrust of Francis turbines.
The time-resolved signals of hydraulic axial thrust were analyzed via FFT; Figure 16 corresponds to the results of four different operating conditions; the frequency spectrum plots (a) to (d) denote the total hydraulic axial thrust, hydraulic axial thrust in the runner blade, hydraulic axial thrust of the runner hub, and hydraulic axial thrust of the runner shroud, respectively. There are three kinds of dominant frequencies and their corresponding harmonics, including the RVR frequency, the RSI frequency fgv and frunner. For the operating point of 13 deg GVO, the highest pulsation amplitude belongs to the guide vane passing frequency, fgv; while for the other three cases, the highest pulsation amplitude corresponds to the RVR frequency. It can also be concluded that the highest pulsation amplitude of total hydraulic axial thrust increases with the deviation of the operating point from the BEP. With regard to the runner components, the fluctuation amplitude of hydraulic axial thrust in the runner hub makes up a large part of the total hydraulic axial thrust; the proportions for the four

Hydraulic Axial Thrust
The pulsating hydraulic axial thrust is transmitted to the load-bearing frame base of the powerhouse via thrust bearing, and consequently poses a threat to the bearing and the powerhouse structure. During certain transition processes, accidental lifting of the rotating part will occur when the thrust bearing is subjected to upward force; on the other hand, excessive downward load on the thrust bearing may increase the temperature of the thrust bearing pads, thus causing damage to them [39]. In addition, the pulsation characteristics of axial hydraulic thrust make it one of the main sources of noise and strong vibration occurring in the powerhouse structure. There are many hydropower plants where prominent vibration problems frequently take place during off-design and transient operations, such as the Three Gorges Hydropower Station, the Ertan Hydropower Station, and the Guangzhou Pumped-Storage Power Station [40,41]; therefore, it is necessary to explore, in depth, and understand the pulsation properties of the hydraulic axial thrust of Francis turbines.
The time-resolved signals of hydraulic axial thrust were analyzed via FFT; Figure 16 corresponds to the results of four different operating conditions; the frequency spectrum plots (a) to (d) denote the total hydraulic axial thrust, hydraulic axial thrust in the runner blade, hydraulic axial thrust of the runner hub, and hydraulic axial thrust of the runner shroud, respectively. There are three kinds of dominant frequencies and their corresponding harmonics, including the RVR frequency, the RSI frequency f gv and f runner . For the operating point of 13 deg GVO, the highest pulsation amplitude belongs to the guide vane passing frequency, f gv ; while for the other three cases, the highest pulsation amplitude corresponds to the RVR frequency. It can also be concluded that the highest pulsation amplitude of total hydraulic axial thrust increases with the deviation of the operating point from the BEP. With regard to the runner components, the fluctuation amplitude of hydraulic axial thrust in the runner hub makes up a large part of the total hydraulic axial thrust; the proportions for the four operating conditions are 88.26%, 83.44%, 90.60%, and 86.83%, respectively. This regularity makes sense because the runner hub results in the change of direction of the flowing water (from radial to axial direction). The fluctuation amplitude of hydraulic axial thrust in the runner shroud accounts for the smallest part of the fluctuation amplitude of total hydraulic axial thrust compared to the other runner components. This could probably be explained by the small area of vertical projection of the runner shroud. As for the runner blade, its insignificant fluctuation amplitude of hydraulic axial thrust is due to both the superposition of positive and negative pressures on the pressure and suction sides (see Figure 17) as well as the tiny area of vertical projection. Figure 18 shows detailed dimensions of the runner, which may help interested researchers figure out the effect of variation of hydraulic axial thrust on the bearings and loadings. runner components. This could probably be explained by the small area of vertical projection of the runner shroud. As for the runner blade, its insignificant fluctuation amplitude of hydraulic axial thrust is due to both the superposition of positive and negative pressures on the pressure and suction sides (see Figure 17) as well as the tiny area of vertical projection. Figure 18 shows detailed dimensions of the runner, which may help interested researchers figure out the effect of variation of hydraulic axial thrust on the bearings and loadings.  runner components. This could probably be explained by the small area of vertical projection of the runner shroud. As for the runner blade, its insignificant fluctuation amplitude of hydraulic axial thrust is due to both the superposition of positive and negative pressures on the pressure and suction sides (see Figure 17) as well as the tiny area of vertical projection. Figure 18 shows detailed dimensions of the runner, which may help interested researchers figure out the effect of variation of hydraulic axial thrust on the bearings and loadings.

Conclusions
The pressure fluctuation and hydraulic axial thrust pulsation of a prototype Francis turbine were numerically investigated under four different conditions for partial load operation. Its design head is 116.0 m, the rated volume flow rate is 87.72 m 3 /s. The input and output pressures for the numerical

Conclusions
The pressure fluctuation and hydraulic axial thrust pulsation of a prototype Francis turbine were numerically investigated under four different conditions for partial load operation. Its design head is 116.0 m, the rated volume flow rate is 87.72 m 3 /s. The input and output pressures for the numerical simulation are 1.18 MPa and 0.04 MPa, respectively. Unsteady three-dimensional (3D) flow field simulation of the whole flow passage was carried out and validated by means of comparison with experimental and numerical results presented in the previously published literature. A good agreement among all parameters was obtained, indicating that the configuration adopted in the current CFD simulation was practical and reliable for the subsequent investigation of pressure fluctuations. Pressure fluctuations in the draft tube, runner, and vaneless space were fully analyzed, which aids the investigation of the pulsating property of hydraulic axial thrust in the Francis turbine. The vortex rope in the draft tube was observed to appear in different forms at varying operating points; the volume of the vortex rope became smaller when the operating point was closer to the BEP. The frequency of the highest amplitude of rope-induced pressure fluctuation in the draft tube was 0.24 times that of the runner rotational frequency. The pressure fluctuation of the draft tube reached its maximum in the operating condition of 16 deg GVO. As for the pressure fluctuation in the runner blade, the highest amplitude corresponded to the fourth harmonic of the RVR frequency; as a result, the basic law of change related to the working condition is consistent with that of the monitoring point in the draft tube. It should be noted that the RSI frequencies of both f gv and f runner were found in the runner blade, hub, and shroud. Regarding the pressure fluctuation in the vaneless space, the highest amplitude belonged to the runner blade passing frequency and there was not much difference between the various operating points in terms of the highest amplitude of pressure fluctuation. In the end, the frequency spectrum plots of pulsating hydraulic axial thrust were analyzed from the perspectives of total, runner blade, runner hub, and runner shroud. Three kinds of dominant frequencies and their corresponding harmonics exist, including the RVR frequency, the RSI frequency f gv , and f runner . For the operating point of 13 deg GVO, the highest pulsation amplitude belonged to the guide vane passing frequency, f gv , while for the other three cases, the highest pulsation amplitude corresponded to the RVR frequency. With regard to the runner components, the fluctuation amplitude of hydraulic axial thrust in the runner hub made up a large part of the total hydraulic axial thrust; the proportions for the four operating conditions were 88.26%, 83.44%, 90.60%, and 86.83%, respectively. The results presented in this paper provide a comprehensive understanding of the pressure fluctuations caused by flow instabilities and the resulting pulsating hydraulic axial thrust, which is expected to provide reference for issues associated with hydraulic axial thrust in the process of relevant engineering design and safe operational planning. Future work is to consider conducting on-site test of the prototype to validate the findings in this paper.