Dynamic characteristics of a running away pump-turbine with large head variation: 1D + 3D coupled simulation

Abstract

The dynamic characteristics of running away pump turbines (PTs) with a large head variable amplitude have not been understood thus far, primarily because of two difficulties in simulation and analysis. The first is how to provide accurate time-varying boundary conditions for transient simulation of the turbine runaway process (TRP). The other is how to determine the specific appearance time of each frequency component of the complex pressure fluctuations. This study presented a one- and three-dimensional (1D-3D) coupled approach considering waterway dynamics to provide accurate unsteady boundary conditions for the transient flow simulation of a PT with a large head variable amplitude. The short-time Fourier transformation (STFT) approach was adopted to analyse the time-frequency characteristics of the transient pressures and impeller forces. The study found that the fluctuations of pressures and impeller forces during the TRP of the PT with a large head variable amplitude contained two exclusive fluctuation frequency components. The former was approximately three times the rated rotational frequency of the impeller. The later was a series of integer fold transient rotational frequencies of the impeller, which was irrelevant to the rotor-stator interactions. The findings have important value for controlling the pressure fluctuations during the TRP of PTs.

KEYWORDS:

• Large head variable amplitude pump-turbine
• turbine runaway process
• one- and three-dimensional coupled simulation
• pressure fluctuation
• flow pattern evolution

1. Introduction

Pumped-storage hydropower (PSH) technology is recognized as the most mature and the only widely applied flexible regulation power supply at the grid scale to suppress the power fluctuations of intermittent renewable energy (Liu et al., 2021). Currently, several countries are still devoted to vigorously developing PSH technology (Pérez-Díaz et al., 2015). The hydraulic instability for pump-turbines (PTs) caused by high amplitude pressure fluctuations is still a key scientific problem. While a PT operates in off-design conditions, unstable flows such as flow separation, vorticities, backflow, and cavitation induce various pressure fluctuations and cause further hydraulic instabilities in a PT (Lu et al., 2022; Zhao et al., 2021; Zuo et al., 2015). During the dynamic transformations of operation conditions, the flow pattern transitions are more complex than those under static conditions. Complex flow-pattern evolutions induce intensive pressure fluctuations. Therefore, investigating the pressure fluctuation and flow pattern evolution mechanisms during the transient processes of PT is necessary.

The turbine runaway process (TRP) is one of the most representative transitions of PTs. The pressure fluctuations and flow patterns are very complex at runaway transitions. To explore turbine runaway instability mechanisms, researchers have performed numerous investigations. According to the differences in the adopted simulation approaches, these investigations are divided into three categories as listed in Table 1.

Table 1. Three simulation methods of turbine runaway process (TRP) investigations.

Method	Computational domain	Pipeline simulation method	PT simulation method
First	PT	–	3D
Second	PT + Pipeline	3D	3D
Third	PT + Pipeline	1D	3D

In the first type of simulation investigations, the TRP of PT was dynamically simulated adopting three-dimensional (3D) calculation approach without considering the pipeline system. Liu et al.(Liu et al., 2016) quantitatively investigated the backflow near the impeller inlet of PT while adopting the three-dimensional flow simulation approach. The investigation found that the backflow is more apparent at a larger wicket gate opening. Lu et al. (Lu et al., 2021) studied pressure fluctuations and flow patterns under turbine runaway conditions and low-frequency high-amplitude pressure fluctuations were found to be related to the unsteady evolution of vorticities in the guide vane. Xia et al. (2017) captured the local backflow and rotating stall in the impeller inducing the pressure fluctuations and impeller radial thrusts during the PT runaway process. Yang et al. (2020) also captured the apparent backflow in TRPs of PTs with varying specific speeds by adopting the three-dimensional flow simulation approach. The low-frequency high-amplitude pressure fluctuations were primarily attributed to the local backflow near the impeller inlet. Guo et al. (2012) found that a small guide vane can induce an apparent water ring in the vaneless PT. Zhang et al. (2019) suggested that the simulation approach can accurately predict the dynamic evolution process of the impeller backflow and the induced pressure fluctuations during the TRP of PTs. The 3D flow simulation approach relies on experimental data and is suitable for investigating the transient flow mechanism but not for engineering applications.

In the second type of simulation investigation, the TRP of the entire hydraulic system was dynamically simulated by adopting a full 3D flow simulation approach to provide accurate time-dependent boundary conditions at the entrances and exits of PTs. Zhou et al. (2018) simulated the TRPs of PTs with different impeller blades by adopting a simulation approach. The influence mechanisms of the impeller blade shape on the transient performances of the PT were revealed. Owing to the large calculation expense, the full 3D flow simulation approach was found to be unsuitable for engineering applications.

In the third type of simulation investigation, the TRP of the entire hydraulic system of the pump-storage hydropower plant was dynamically simulated by adopting the one- and three-dimensional (1D-3D) coupled flow calculation approach. The unsteady flow in the pipelines was simulated by adopting a 1D flow simulation approach. The turbulence flow in the PT was simulated by adopting a 3D computational fluid dynamics (CFD) approach. Fu et al. (2022) investigated the influences of guide vane opening and rotor inertia on the pressure fluctuations of middle- and high-head PTs by adopting a simulation approach. Zhang et al. (2016) revealed the generation mechanisms of loop dynamic characteristic curves during TRP by adopting the 1D-3D coupled calculation approach. Additionally, the effects of rotor inertia on impeller hydraulic forces were investigated by adopting this approach (Zhang et al., 2020).

By reviewing the references, apparent progress has been made on the turbine runaway of PTs with normal head variable amplitudes. However, the dynamic characteristics during TRP of the PTs with a large head variable amplitude remain unclear primarily due to two challenges in simulation and analysis, namely accurately simulating the TRP of a PT with a large head variable amplitude and analysing the complex frequency characteristics of pressure fluctuations during the TRP. To consider the waterway dynamics and provide accurate time-varying boundary conditions, this study adopted the 1D-3D coupled simulation approach to obtain the real dynamic characteristics during the TRP of a PT with a large head variable amplitude. To accurately determine the specific appearance time of each frequency component of the pressure fluctuations in the TRP, the short-time Fourier transformation (STFT) approach was adopted to analyse the time–frequency characteristics of the transient pressures and impeller forces.

The remainder of this paper is organized as follows. Section 2 presents the investigated pump-storage hydropower plant and transient flow simulation approach and its validation, section 3 presents an analysis and discussion of the simulation results and section 4 draws the research conclusion of the work.

2. Materials and methods

2.1. Flow governing equations and geometry model

In this study, a PSH plant with a large-head variable-amplitude ratio of 1.24 and capacity of 255 MW was investigated. The rated head of the investigated PT was 259 m. This study primarily focused on the transient flow in the hydraulic system of the plant. Figure 1 shows a simplified schematic of the hydraulic system of the investigated hydropower plant. The hydraulic system primarily includes two parts, namely hydraulic pipelines and PT. Considering that the hydraulic pipelines are excessively long (approximately 2200 m), the 3D characteristics of the flow in pipelines are not apparent. Therefore, the flow in the pipelines was described by adopting the 1D governing equation (Eqs. (1)(1) $\begin{array}{rl}& Ag\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{Q}{A}\frac{\mathrm{\partial }Q}{\mathrm{\partial }x}+\frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\frac{f_w}{2DA}Q|Q|=0\end{array}$(1) –(2)(2) $\begin{array}{rl}& \frac{Q}{A}\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{\mathrm{\partial }H}{\mathrm{\partial }t}-\frac{Q}{A}\sin \alpha +\frac{a^2}{gA}\frac{\mathrm{\partial }Q}{\mathrm{\partial }x}=0\end{array}$(2) ) of the transient flow. (1) $\begin{array}{rl}& Ag\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{Q}{A}\frac{\mathrm{\partial }Q}{\mathrm{\partial }x}+\frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\frac{f_w}{2DA}Q|Q|=0\end{array}$(1) (2) $\begin{array}{rl}& \frac{Q}{A}\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{\mathrm{\partial }H}{\mathrm{\partial }t}-\frac{Q}{A}\sin \alpha +\frac{a^2}{gA}\frac{\mathrm{\partial }Q}{\mathrm{\partial }x}=0\end{array}$(2) where H denotes mean water head in the pipe across the section (m), Q denotes the flowrate of the pipe (m³/s), A denotes the area of the pipe across the section (m²), D denotes the inner diameter of the pipe (m), g denotes gravitational acceleration (m/s²), f_w denotes the friction coefficient of the pipe wall (-), α denotes the angle of the pipe axis relative to the horizontal plane (°), x denotes the spatial coordinate along the pipe axis (m), t denotes the time (s), and a denotes acoustic speed (m/s), which was calculated using Eq. (3)(3) $\begin{array}{r}a=\frac{\sqrt{K/\rho }}{\sqrt{1+\phi \frac{KD}{Ee}}}\end{array}$(3) . (3) $\begin{array}{r}a=\frac{\sqrt{K/\rho }}{\sqrt{1+\phi \frac{KD}{Ee}}}\end{array}$(3) where K denotes the elastic modulus of water (Pa), ρ denotes water density (kg/m³), E denotes the elastic modulus of the pipe wall (Pa), e denotes thickness of the pipe wall (m), and φ denotes a coefficient influenced by the elasticity and connection characteristics of the pipe (-). For an entirely fixed pipe, φ = 1-σ², for a fixed pipe only at the upstream end, φ = 1-σ²/2, and for a pipe connecting adopting expansion joint, φ = 1, where σ denotes Poisson's ratio of pipe wall material (-). The time-step size Δt was 0.0017s, and the spatial-step size Δx of 1D grids can be determined according to the Courant condition, Eq. (4)(4) $\begin{array}{r}\frac{\left(V\pm a\right)\Delta t}{\Delta x}\le 1\end{array}$(4) . Consequently, the total discretization node amount of the 1D hydraulic pipelines was 1142. (4) $\begin{array}{r}\frac{\left(V\pm a\right)\Delta t}{\Delta x}\le 1\end{array}$(4) where V denotes the mean flow velocity in the pipe across the section (m/s). The 1D flow partial differential equations (Eqs. (1)(1) $\begin{array}{rl}& Ag\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{Q}{A}\frac{\mathrm{\partial }Q}{\mathrm{\partial }x}+\frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\frac{f_w}{2DA}Q|Q|=0\end{array}$(1) –(2)(2) $\begin{array}{rl}& \frac{Q}{A}\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\frac{\mathrm{\partial }H}{\mathrm{\partial }t}-\frac{Q}{A}\sin \alpha +\frac{a^2}{gA}\frac{\mathrm{\partial }Q}{\mathrm{\partial }x}=0\end{array}$(2) ) can be transformed into ordinary differential equations and solved by adopting the method of characteristics (MOC) (Ghidaoui et al., 2005). The 1D MOC is an explicit numerical simulation approach that can apparently reduce the simulation expense.

Figure 1. Hydraulic system schematic of the pumped-storage hydropower (PSH) plant.

Figure 2 shows the 3D geometric model of the investigated PT. The flow domain primarily includes four parts: the spiral casing, distributor (stay/guide vane), impeller, and draft tube. The rated specific speed n_s (Eq. 5(5) $\begin{array}{r}{n}_s=\frac{n_r\sqrt{P}}{{H_r}^{1.25}}\end{array}$(5) ) of the PT was 145.8. (5) $\begin{array}{r}{n}_s=\frac{n_r\sqrt{P}}{{H_r}^{1.25}}\end{array}$(5) where n_r denotes the rated rotation speed of the investigated PT (r/min), P denotes the rated power of the PT (kW), H_r denotes the rated head of the PT (m).

Figure 2. 3D flow domain model of the investigated PT.

The unsteady flow inside the PT was the research focus of this study and was described by adopting the 3D compressible flow governing equation, Eqs. (6)(6) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\frac{\mathrm{\partial }\left(\rho u_i\right)}{\mathrm{\partial }{x}_i}=0\end{array}$(6) –(7)(7) $\begin{array}{rl}& \frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}=f_i-\frac{1}{\rho }\frac{\mathrm{\partial }p}{\mathrm{\partial }{x}_i}+\nu {\mathrm{\nabla }}^2{u}_i+\frac{\nu }{3}\frac{{\mathrm{\partial }}^2{u}_j}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}\end{array}$(7) (Shadab et al., 2022). (6) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\frac{\mathrm{\partial }\left(\rho u_i\right)}{\mathrm{\partial }{x}_i}=0\end{array}$(6) (7) $\begin{array}{rl}& \frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}=f_i-\frac{1}{\rho }\frac{\mathrm{\partial }p}{\mathrm{\partial }{x}_i}+\nu {\mathrm{\nabla }}^2{u}_i+\frac{\nu }{3}\frac{{\mathrm{\partial }}^2{u}_j}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}\end{array}$(7) where u denotes the flow velocity component (m/s), p denotes the pressure (Pa), f_i denotes the mass force (m/s²), ν denotes the kinematic viscosity (m²/s), x denotes the spatial coordinates (m), and the subscripts i and j denote the Cartesian coordinate components (i, j = 1, 2, 3). The water density ρ was modified according to the pressure distribution using the weak compressibility model, Eqs. (8) (Zhang & Cheng, 2012). (8) $\begin{array}{r}\rho ={\rho }_0{e}^{\left(p-p_0\right)/{\rho }_0{a_0}^2}\end{array}$(8) where subscript 0 denotes initial state, ρ₀ = 998.2 kg/m³, p₀ = 0 Pa, a₀ = 1400 m/s. The 3D flow-governing equations, Eqs. (6)(6) $\begin{array}{rl}& \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\frac{\mathrm{\partial }\left(\rho u_i\right)}{\mathrm{\partial }{x}_i}=0\end{array}$(6) –(7)(7) $\begin{array}{rl}& \frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }t}+u_j\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}=f_i-\frac{1}{\rho }\frac{\mathrm{\partial }p}{\mathrm{\partial }{x}_i}+\nu {\mathrm{\nabla }}^2{u}_i+\frac{\nu }{3}\frac{{\mathrm{\partial }}^2{u}_j}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}\end{array}$(7) are complicated nonlinear high-order partial differential equations that can only be solved by adopting the iteration approach. In this study, it was solved by adopting the finite volume method (FVM). Although the large eddy simulation (LES) approach can accurately simulate small-scale turbulence flow, its simulation expenses are excessively large for engineering applications. In this study, the turbulence flow inside the PT was modelled and solved by adopting the delayed detached eddy simulation (DDES) approach. The DDES approach takes some advantages of LES and Reynolds averaged Navier–Stokes simulation (RANS) (Gritskevich et al., 2012). The LES and RANS approaches can be selected according to the turbulence flow scale and grid scale to improve the simulation accuracy (Spalart 2009).

2.2. Transient flow simulation scenario

1D and 3D simulation grids were arranged in a partly overlapped manner, as shown in Figure 3, to achieve data delivery between these two flow domains (Zhang & Cheng, 2012). The 1D-3D coupled flow simulation considered the interactions between the transient flow in the pipeline and turbulent flow pulsations in the PT. Hence, the model provided accurate time-dependent boundary conditions for flow simulations inside pipelines and PT during TRP. Given the effects of the boundary condition styles on the flow simulation convergence, the velocity entrance and pressure exit were set up at the entrance and exit of the 3D PT. The 1D transient flow parameters in the upper and lower pipelines were solved by adopting the user-defined functions of ANSYS Fluent at each time step. The 1D transient flow parameters at the entrance and exit of the PT were incorporated into medium data files. Subsequently, the Fluent solver read the medium data files to obtain the flow velocity and pressure at the entrance and exit of the PT for each time step. The rotation of the impeller flow domain was simulated by adopting the slide mesh approach to consider the rotor-stator interaction (RSI). The time-varying impeller rotation speed was calculated according to the following explicit difference equation: (9), which is derived from the rotation equation of a single-degree rigid body. (9) $\begin{array}{r}{n}_{N+1}=n_N+\frac{30}{\pi }\frac{M_N}{J}\Delta t\end{array}$(9) where n denotes the rotation speed of the investigated PT (r/min), M denotes the impeller hydraulic torque (N·m), J denotes the rotor inertia (kg·m²), and π denotes the circular constant (-). The subscript N denotes the sequence index of the time step. Each of the two neighbouring components of the PT were connected by adopting the interface pairs. To decrease the influence of the interface on the simulation results, the mesh density distribution at each side of an interface pair was generated as consistently as possible.

Figure 3. 1D and 3D mesh arrangement strategies.

The 1D and 3D flow simulations were both performed in the simulation platform ANSYS. The 3D flow was directly simulated using the ANSYS Fluent model. The 1D flow was simulated by adopting the second development interface of Fluent. For the three-dimensional flow simulation setup in Fluent, the bounded central differencing scheme was adopted to discretize momentum. A second-order implicit scheme was employed as a transient formulation. The pressure interpolation scheme was the second-order. The pressure and velocity coupling equations were numerically solved by adopting the SIMPLEC algorithm. Under-relaxation factors for each flow variable were set as one-third of the default values in Fluent to obtain satisfactory convergence of flow simulation for all flow conditions during the TRP. The unsteady simulation was initialized by adopting an interpolation approach based on the steady simulation results at the initial conditions of the TRP. Obtaining satisfactory convergence is challenging in the initial stage of transient flow simulation compared with that in the later stage of transient flow simulation. In the first 100 time steps, the maximal iteration number per time step was 150. Subsequently, the maximal iteration number per time step was 30 to obtain the same convergence as the first 100 time-steps.

2.3. Calculation grids

All components of the investigated PT were discretized by adopting a hexahedron cell. The grids near the wall were refined, as shown in Figure 4, to resolve the small-scale flow structures under boundary conditions. To study the effect of the grid node number on the simulation results, five set structure grids were generated, and the node number in each component is summarized in Table 2. The initial condition of the TRP was simulated by adopting five grids. Figure 5 shows that after the grid node number increased to 5.3 × 10⁶, the relative deviations of the flowrate were less than 3% compared with the simulated flowrate of the fifth grid with 11.22 × 10⁶ nodes. The simulated flowrates for the third and fourth grids were proximate to that of the fifth grid. The simulation results have unclear variations with an increase in the grid node number. Finally, the fourth grid with 9.3 × 10⁶ nodes was adopted to simulate the TRP for balancing the simulation expense and accuracy. The time-step size of the 3D flow simulation is the same as that of the 1D simulation.

Figure 4. Mesh distribution in each component of the PT.

Figure 5. Relationship between simulated flowrate and grid node number.

Table 2. Grid node number and quality in each component of the PT.

Component	First (10^4)	Second (10^6)	Third (10^6)	Fourth (10^6)	Fifth (10^6)	Quality
Guide vane	0.72	1.00	1.60	2.14	2.52	0.7
Stay vane	0.78	1.06	1.49	2.08	2.21	0.4
Impeller	1.18	1.69	2.15	2.60	3.02	0.4
Draft-tube	0.43	0.76	1.02	1.14	1.71	0.62
Spiral casing	0.39	0.79	1.09	1.34	1.76	0.25
Total number	3.50	5.30	7.35	9.30	11.22

2.4. Simulation scenario verification

The 1D-3D coupled flow simulation approach has been verified in many publications (Fu et al., 2021; Zhang & Cheng, 2012). This study also adopted available experimental data to validate the simulation scenario. The simulated dynamic trajectory of the TRP was compared with the static model-scale experimental result at the investigated guide vane opening provided by a manufacturer. Figure 6 shows the four-quadrant characteristic curves (n₁₁-M₁₁) of the PT for the simulation and experiment. The unit parameters n₁₁ and M₁₁ were defined as follows: (10) $\begin{array}{rl}{n}_{11}& =\frac{n{D}_t}{\sqrt{H_t}}\end{array}$(10) (11) $\begin{array}{rl}{M}_{11}& =\frac{M}{D_t^3{H}_t}\end{array}$(11) where D_t is the impeller inlet diameter (m) and H_t is the head of the PT (m). The simulated dynamic trajectory of PT always fluctuated near the static experimental curve. Finally, the dynamic trajectory fluctuated near the runaway condition (M₁₁ = 0). The dynamic trajectory and static experimental curve were not entirely consistent, which is possibly related to the accumulation error of the transient simulation and the complexity of the flow inside the PT near the runaway condition. Generally, the simulation results are in line with the experimental data. Therefore, the 3D flow simulation scenario adopted in this study is acceptable.

Figure 6. Simulated and experimental results of four-quadrant characteristic curves (n₁₁-M₁₁) for the PT.

3. Results and discussion

3.1. Transient performances

Simulated transient performances during the TRP of the PT with a large head variable-amplitude ratio of 1.24 are shown in Figure 7. The TRP of the PT was divided into turbine and turbine brake modes by no-load conditions (M = 0). The PT alternately operated in the turbine and turbine brake modes during the TRP. The transient performance fluctuated with the variation in rotational speed back and forth and finally tended to converge. The variation trends of the impeller torque and PT flowrate were similar, and those of the rotation speed and PT head were also consistent. The fluctuations of all the performances presented a low-frequency characteristic, which was the fluctuation frequency of the rotational speed. Additionally, the fluctuations in the impeller torque and PT head exhibited high-frequency characteristics, which were the pressure fluctuation frequencies. Before t = 5 s, the transient flowrate and impeller torque decreased gradually. After t = 5 s, the transient flowrate and impeller torque fluctuated. While the PT operated close to no-load conditions (M = 0), the PT head abruptly decreased, but the rotational speed continued to increase to the maximum runaway rotational speed. The abrupt decrease in the hydraulic head of the PT possibly caused a local backflow near the impeller high-pressure inlet under the centrifugal resultant force. The local backflow blocked the flow passages and resulted in fluctuations of the PT flowrate. This analysis is demonstrated in the subsequent internal flow analysis. The performance analysis suggested that the direct reason for the fluctuations in the PT flowrate and impeller torque was the abrupt decrease in the PT head. During the optimization of the PSH system, the abrupt variation in the PT head should be suppressed to avoid intensive fluctuations in the transient performance Figure 8.

Figure 7. Transient performances of the investigated PT during the TRP: (a) impeller torque, (b) rotational speed, (c) PT flowrate, and (d) PT head.

Figure 8. Transient impeller forces of the PT during the TRP: (a) radial force in X direction, (b) radial force in Y direction, and (c) axial force.

3.2. Impeller hydraulic thrusts

Near the abrupt transformation condition (t = 5 s), close to the maximum runaway rotational speed condition, the radial forces (F_X and F_Y) on the impeller resulted in instantaneous high-amplitude pulses. The impeller axial force (F_Z) increased abruptly. Subsequently, the radial and axial forces fluctuated severely at high frequencies, and the radial forces tended to converge. To investigate their time–frequency characteristics, STFT was performed for the impeller forces, as shown in Figure 9. The fluctuations of the impeller forces had four frequency components (f_RSI-S, f_RFI, f_BFV, and f_DV). The first two fluctuation frequency components (f_RSI-S and f_RFI) were the common frequencies for the radial and axial hydraulic thrusts. f_RSI-S existed in the entire TRP; however, f_RFI appeared after PT passed through the abrupt transformation condition (t = 5 s). The third fluctuation frequency component f_BFV only occurred near the abrupt transformation condition (t = 5 s) for the impeller radial forces (F_X and F_Y). The fourth fluctuation frequency component f_DV is prominent in the entire TRP of the PT, but only for the impeller axial force. The third fluctuation frequency component f_BFV was the high-amplitude fluctuation of radial forces and possibly induced the hydraulic exciting vibrations of the rotor. The fourth fluctuation frequency component f_DV was the high-amplitude fluctuation of the impeller axial force and possibly induced the axial elevation of the rotor. The two fluctuation frequency components of the impeller forces should be emphatically considered to improve the operating safety of the PSH plant with a large head variable amplitude.

Figure 9. Short-time Fourier transformation (STFT) results of impeller forces during the TRP of the PT: (a) F_X, (b) F_Y, and (c) F_Z.

3.3. Pressure fluctuations

To further analyse the formation mechanisms of the four fluctuation frequency components of the impeller forces, the fluctuating pressures at several representative monitoring points were analysed, as shown in Figure 10. Four monitoring points were located at the spiral casing end (SE), in the vaneless space (VL), near the impeller high-pressure inlet (IM), and draft-tube inlet (DT). The PT head was primarily influenced by the transient pressure at the SE; hence, the variation trends were similar. The axial hydraulic thrust on the impeller was closely related to the transient pressure at DT; thus, the variation trends were similar. Near the interfaces between the rotation impeller and stationary components of the PT, the flow patterns are usually complex and unstable under strong shear stress. Additionally, the monitoring VL and IM points are closer to the interfaces between the rotation impeller and stationary components of the PT than to the SE and DT. Therefore, the fluctuation frequencies and amplitudes of the transient pressures at VL and IM were much higher than those of SE and DT.

Figure 10. Transient pressure at each monitoring point during the TRP of the PT: (a) spiral casing end (SE), (b) vaneless space (VL), (c) impeller high-pressure inlet (IM), and (d) draft-tube inlet (DT).

The STFT approach was also employed to analyse the time–frequency distributions of pressure fluctuations, as shown in Figure 11. The pressure fluctuations consisted of five frequency components (f_RSI-S, f_RSI-R, f_RFI, f_BFV, and f_DV), which were consistent with those (f_RSI-S, f_RFI, f_BFV, and f_DV) of the impeller force fluctuations for the TRP of the PT with a large head variable amplitude. The first two fluctuation frequency components (f_RSI-S and f_RSI-R) existed in the entire TRP for each monitoring point. f_RSI-S is the 7-fold (impeller blade number) transient rotational frequency of the impeller and its harmonics and existed in the transient pressures of the monitoring points in the stationary components of the PT. f_RSI-R is the 20-fold (guide vane number) transient rotational frequency of the impeller and its harmonics, and existed in the transient pressures of the monitoring points in the rotational components of the PT. According to previous studies (Fu et al., 2020; Xia et al., 2017), the fluctuation frequency components (f_RSI-S and f_RSI-R) were induced by RSI phenomena. Therefore, the intensities of the pressure fluctuations caused by the RSI were higher and closer to the interfaces between the rotational impeller and stationary components of the PT.

Figure 11. STFT results of fluctuating pressures at each monitoring point during the TRP of the PT: (a) SE, (b) VL, (c) IM, and (d) DT.

3.4. Unsteady flow inducing pressure fluctuations

The third frequency component (f_BFV) of the pressure fluctuations primarily appeared near the abrupt transformation condition (t = 5 s), close to the maximum runaway rotational speed condition of the TRP. According to relevant studies (Fu et al., 2020; Liu et al., 2016; Zhang et al., 2019), the frequency component (f_BFV) of pressure fluctuations was induced by the dynamic evolution of local backflow vortex (BFV) proximate to the impeller high-pressure entrance. In this study, centrifugal backflow was also captured close to the impeller entrance, as shown in Figure 12. Near the abrupt transformation condition (t = 5 s) close to the maximum runaway rotational speed condition, the PT head abruptly decreased. However, the rotational speed continued to increase to the maximum runaway rotational speed, as shown in Figure 7. Therefore, the resultant force of the pressure difference related to the PT head and centrifugal force related to the rotational speed was in the centrifugal direction per unit mass of fluid close to the impeller entrance. The centrifugal backflow was driven by the resultant centrifugal force close to the impeller entrance. The backflow vortex was generated under the shear of a large velocity gradient between the centrifugal backflow and the centripetal inflow close to the impeller entrance. The backflow and backflow vortexes blocked the centripetal flow close to the impeller entrance and caused flowrate fluctuations, which induced a water hammer. The dynamic interactions between the backflow vortex and the solid wall of the PT resulted in a pressure fluctuation component (f_BFV). Near the abrupt transformation condition (t = 5 s), the flow passage sizes, blocked by the backflow and backflow vortexe, were much larger than those at other instants (e.g. t = 1 s). Consequently, the intensities of the pressure fluctuation components (f_BFV) near the abrupt transformation condition (t = 5 s) were higher than those at other conditions during the TRP of the PT with a large head variable amplitude. The pressure fluctuation component (f_BFV) primarily influenced the impeller radial forces and was closely related to the PT head and rotational speed during the TRP. To suppress the backflow vortex and pressure fluctuation components (f_BFV), the flow capacity of the impeller inlet can be improved by increasing the minimum PT head and decreasing the runaway rotational speed. Suitably increasing the guide vane height was a viable method for improving the flow capacity of the impeller inlet.

Figure 12. Velocity vector distributions in impeller meridian plane during the TRP of the PT.

The fourth frequency component (f_DV) of the pressure fluctuations was approximately three times the rated rotational frequency of the impeller and existed in the entire TRP of the PT with a large head variable amplitude. Moreover, the intensities of the pressure fluctuation component (f_DV) were higher and closer to the monitoring point SE. The excitation sources of the pressure fluctuation component (f_DV) were located in the spiral casing. Figure 13 shows the streamline distribution in a cross-section of the spiral casing. Apparent Dean vortexes (DV) always existed in the spiral casing at each instant during the TRP of the PT. These DV in the spiral casing periodically interacted with the rotational impeller and induced a pressure fluctuation component (f_DV). The evolution of the flow pattern in the spiral casing was influenced by the geometric shape of the spiral casing. Therefore, the pressure fluctuation component (f_DV) related to the DV in the spiral casing can be suppressed by optimizing the geometry shape of the spiral casing.

Figure 13. Streamline distributions at several typical instants in the cross-section of spiral casing during the TRP of the PT.

The fifth frequency component (f_RFI = nf_n) of the pressure fluctuations was a series of integral fold transient rotational frequencies of the impeller that occurred after the PT passed through the abrupt transformation condition (t = 5 s). Moreover, the intensities of the pressure fluctuation component (f_RFI) were higher and closer to the monitoring VL point. The exciting sources of the pressure fluctuation component (f_RFI) were possibly located in the distributor regions. Figure 14 shows the streamline and pressure distribution on the turbo-surface of the distributor and impeller with a span height of 0.5. Before the PT passed through the abrupt transformation condition (t = 5 s), the streamlines were smooth in the distributor and impeller regions, as shown in Figure 14 (a). The pressure varied smoothly. The pressure fluctuation component (f_RFI) did not appear. However, after the PT passed through the abrupt transformation condition (t = 5 s), the backflow occurred near the impeller inlet and induced the blockage of flow passage and complex vortexes in VL. The complex flow pattern in the VL changed the pressure distribution and induced some complex but regular high- and low-pressure regions in the distributor. The numbers (10 and 12) on the pressure cloud in Figure 14 are the numbers of low-pressure regions, and the numbers, 3, 4, and 6 are the numbers of high-pressure regions. These complex but regular high- and low-pressure regions were alternately distributed in the distributor. The periodical interactions between the rotational impeller and the complex but regular high- and low-pressure regions induced by the complex flow patterns in the distributor (RFI) induced the pressure fluctuation component (f_RFI) during the TRP of the PT with a large head variable amplitude. The periodic interaction mechanism of RFI is similar to the RSI mechanism in turbomachinery. The evolution of the vaneless flow pattern is constrained by the geometric shape of the guide/stay vanes and impeller. Therefore, the pressure fluctuation component (f_RFI) related to the RFI can be controlled by optimizing the geometric shape of the guide/stay vanes and the impeller. Additionally, the matching between the guide/stay vane and impeller blades was also a probable factor influencing the pressure fluctuation component (f_RFI).

Figure 14. Streamline and pressure distributions in the turbo-surface of distributor and impeller at several typical instances during the TRP of the PT.

4. Conclusions

This study dynamically simulated the TRP of a PT with a large head variable-amplitude ratio of 1.24 and capacity of 255 MW, adopting a 1D-3D coupled approach to provide accurate time-varying boundary conditions. This approach is considerably different compared with the previous simulation investigations adopting constant boundary conditions at the entrances and exits PTs. The time–frequency characteristics of the pressure fluctuations were analysed in detail. Two new pressure fluctuation frequencies (f_RFI and f_DV) were found during the TRP of the PT with a large head variable amplitude, which cannot be determined in PTs with a normal head variable amplitude. The former was approximately three times the rated rotational frequency of the impeller. The later was a series of integer fold transient rotational frequencies of the impeller, and their formation mechanisms were revealed.

In conclusion, the pressure fluctuation issues during the TRP of PTs with a large head variable amplitude contained two exclusive high-amplitude pressure fluctuation components in addition to common high-amplitude pressure fluctuations components and was more complex than that of PTs with a normal head variable amplitude. The exclusive high-amplitude pressure fluctuation components can aggravate the fluctuations of impeller hydraulic thrusts of PTs during the TRP. Therefore, the hydraulic instability of a runaway PT with a large head variable amplitude is highlighted in detail compared with that of a PT with a normal head variable amplitude.

The exclusive high-amplitude pressure fluctuation components of PTs with large head variable amplitudes are primarily caused by the unstable flow pattern evolution in the spiral casing and distributor. To control the flow pattern in the spiral casing and distributor and further suppress the exclusive high-amplitude pressure fluctuations, the hydraulic optimization of the spiral casing and distributor should be emphatically considered for the PT with a large head variable amplitude compared with the PT with a normal head variable amplitude. However, hydraulic optimization of the spiral casing and distributor was not conducted in the present study and should be the focus of future studies of this nature.

Notations

A	=	area of pipe across section (m2)
a	=	acoustic speed (m/s)
CFD	=	computational fluid dynamics
D	=	inner diameter of pipe (m)
Dt	=	impeller inlet diameter
DDES	=	delayed detached eddy simulation
E	=	elastic modulus of pipe (Pa)
e	=	thickness of pipe wall (m)
F	=	impeller hydraulic thrusts (N)
f	=	frequency (Hz)
fi	=	mass force (m/s2)
fBFV	=	frequency of pressure fluctuation related to backflow vortex (Hz)
fDV	=	frequency of pressure fluctuation related to DV in spiral casing (Hz)
fn	=	rotation frequency of PT impeller (Hz)
fRSI-S and fRSI-R	=	frequency of pressure fluctuation related to RSI (Hz)
fRFI	=	frequency of pressure fluctuation related to RFI (Hz)
FVM	=	finite volume method
fw	=	friction coefficient of pipe wall (-)
FX	=	impeller force in X direction (N)
FY	=	impeller force in Y direction (N)
FZ	=	impeller axial force (N)
g	=	gravitational acceleration (m/s2)
H	=	hydraulic head (m)
Ht	=	PT head (m)
Hr	=	rated head of investigated PT (m)
J	=	Rotor inertia (kg·m2)
K	=	elastic modulus of water (Pa)
LES	=	large eddy simulation
M	=	impeller torque (N·m)
M11	=	unit torque (N·m)
MOC	=	method of characteristic
n	=	impeller rotational speed (r min−1)
n11	=	unit speed (r min−1)
nr	=	rated rotational speed of impeller (r min−1)
ns	=	rated specific speed (r/min)
P	=	rated power of investigated PT (kW)
p	=	pressure (Pa)
Q	=	flowrate of the pipe (m3 s−1)
RANS	=	Reynolds average Navier-stokes simulation
RSI	=	rotor-stator interaction
STFT	=	short time Fourier transformation
t	=	time (s)
u	=	flow velocity component (m/s)
V	=	mean flow velocity in pipe across section (m/s)
x	=	spatial coordinate (m)
α	=	angle of the pipe axis relative to the horizontal plane (°)
ν	=	kinematic viscosity (m2/s)
ρ	=	water density (kg m−3)
φ	=	coefficient related to the elasticity and connection characteristics of the pipe (-)
σ	=	Poisson's ratio of pipe wall material (-)
π	=	circular constant (-)
Δt	=	time-step size (s)
Δx	=	spatial step size (m)

Subscripts

BFV	=	backflow vortex
DT	=	draft-tube inlet
DV	=	Dean vortexes
X, Y, Z	=	three directions in Cartesian coordinates
IM	=	high-pressure inlet of impeller
i, j	=	Cartesian coordinate components (i, j = 1, 2, 3)
N	=	sequence index of time step
r	=	at rated condition
RFI	=	periodical interactions between impeller and complex flow patterns in distributor
SE	=	spiral casing end
t	=	turbine
VL	=	vaneless space
0	=	initial state

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by China Postdoctoral Science Foundation: [Grant Number No. 2022M720948]; National Natural Science Foundation of China: [Grant Number No. 52079034, 52209108]; Sichuan Science and Technology Support Program: [Grant Number No. 2023YFQ0021].