For three-dimensional numerical simulation of unsteady flow in a pump turbine under wave disturbance, the governing equations involve the mass and momentum conservation equations.

The mass conservation equation is ∂ ρ ∂ t + ∂ ∂ x i ρ u i = 0 . (1)

The momentum conservation equation is ∂ ∂ t ρ u i + ∂ ∂ x j ρ u i u j = − ∂ p ∂ x i + ∂ 2 τ i j ∂ x j + S m i , (2) where u i and u j are velocity vectors, t is the flow time, ρ is the fluid density, p is the static pressure, τ i j is the stress tensor, and S m i is a generalized source term of the momentum equation.

Previous studies have demonstrated the high accuracy and stability of the SST k–ω turbulence model in the numerical simulation of turbomachinery (Menter, 2009). The SST k–ω turbulence model has the advantages of both the k–ε model and k–ω model. It adopts different functions according to the characteristics of fluid at different regions; for instance, it uses the k–ω model for simulating fluid in the regions with a low Reynolds number near the side wall and uses the k–ε model for simulating fluid in the regions with a high Reynolds number in the central part of the flume, which makes it possible to simulate the characteristics of fluid under complex conditions. The SST k–ω turbulence model involves two transport equations, namely, the k-transport equation and the ω-transport equation.

The k-transport equation can be written as ∂ ρ k ∂ t + ∂ ρ k u i ∂ x i = ∂ ∂ x i μ + σ k μ t ∂ k ∂ x i + P ∼ k − β ′ ρ k ω . (3)

The ω-transport equation can be written as ∂ ρ ω ∂ t + ∂ ρ ω u i ∂ x i = α ρ S 2 − β ρ ω 2 + ∂ ∂ x i μ + σ ω μ t ∂ ω ∂ x i + 2 1 − F 1 ρ σ ω 2 ω ∂ k ∂ x i ∂ ω ∂ x i , (4) where P k is the production term of turbulent kinetic energy, S is the invariant measure of the strain rate, μ is the dynamic viscosity, F 1 is the blending function, μ t is the turbulent eddy viscosity, and β , β ′ , σ k , and σ ω are the constants of the turbulence model.

Airy’s micro-amplitude wave theory is one of the most basic and widely used linearization theories to solve partial differential equations in wave problems (Janssen, 2008). It assumes that the wave amplitude is significantly smaller than the wavelength and water depth and considers the wave as a linearity wave. According to Airy’s micro-amplitude wave theory (Huang, 2011), the pressure p can be expressed by the linearized Bernoulli equation as follows: p = − ρ g z − ρ ∂ Φ ∂ t . (5)

Substituting the velocity of the micro-amplitude wave into Eq. 5, the pressure of a micro-amplitude wave can be expressed as p = − ρ g z + ρ g H 2 cosh ⁡ k z + d cosh ⁡ k d cos k x − ω t . (6)

Eq. 6 shows that the pressure consists of two parts: the hydrostatic pressure and hydrodynamic pressure. Considering K Z = cosh ⁡ k z + d cosh ⁡ k d ,   η k = ρ g H 2 cos k x − ω t , (7)

we can obtain p = ρ g η k K Z − Z , (8) where p is the pressure, ρ is the fluid density, η _k is the vertical distance between the instant free water surface and the average water surface, K Z is the pressure sensitivity coefficient, which decreases with the increase in the distance between the particle position and the still water surface, and Z is the elevation of the free water surface. For small-amplitude sinusoidal waves, the hydrostatic and hydrodynamic pressures follow sinusoidal or cosine function, and the amplitude depends on the water depth (Sullivan et al., 2007).

In this study, the Q-criterion (Lu et al., 2021) is used to identify the strength and structure of the vortex, which has high efficiency in vortex identification. The equation of the velocity gradient tensor is as follows: λ 3 − P λ 2 + Q λ − R = 0 , (9) where λ is the eigenvalue of the characteristic equation and P, Q, and R are three invariants of the velocity gradient tensor. Q can be expressed as Q = 1 2 Ω 2 − S 2 , (10) where Ω is the vorticity tensor and S is the strain rate tensor. When Q>0, the vorticity dominates the dynamics of the fluid, and when Q<0, it indicates that the strain rate or viscous stress dominates the dynamics of the fluid. The expression of Q can be written as Q = − 1 2 ∂ u ∂ x 2 + ∂ v ∂ y 2 + ∂ w ∂ z 2 − ∂ u ∂ y ∂ v ∂ x − ∂ u ∂ z ∂ w ∂ x − ∂ v ∂ z ∂ w ∂ y . (11)

Figure 1 shows the three-dimensional model of the pump turbine, which mainly includes five regions: spiral case, stay vanes, guide vanes, runner, and draft tube. The initial settings of the pump turbine are as follows: rated head H = 40 m, rated speed n = 760 r/min, rated flow Q_r = 0.28 m³/s, runner inlet diameter D₁ = 0.55 m, runner outlet diameter D₂ = 0.25 m, number of runner blades Z = 7, number of stay vanes Z_s = 20, and number of guide vanes Z₀ = 20.

Figure 2 shows the sketch of the grid division of the pump turbine. The grid division was conducted using Fluent Meshing software. To ensure the computing accuracy, unstructured grids are used to divide the model, and boundary-layer infilling is carried out on the walls of the runner blades, spiral case, and guide vanes. Considering the complexity of the flow pattern evolution in the runner, to accurately simulate the flow separation of the runner, the boundary-layer grid in the runner area is set with 18 layers with a minimum thickness of 0.02 mm, and the average y⁺ value of the whole runner wall is less than 5. To avoid the number of grids affecting the simulation results, the independence of the number of total grids to the Q_r hydraulic turbine is verified (See Figure 3). When the number of grids is larger than 6 million, the calculation results tend to be stable. In addition, reducing the number of grids could decrease the calculation time. Here, the total number of grids is set to 7 million, and the minimum orthogonal quality of the grids is 0.3. Sliding grid technology is used to cope with the dynamic static interface of the runner, and data exchange is carried out between the calculation domains and the interface. The equations for calculating the head and efficiency of the turbine are shown in Eq. 12 and Eq. 13, respectively. H = P in − P out ρ g , (12) η = M ω 9.81 Q H × 100 % , (13) where P in is the inlet pressure, P out is the outlet pressure, M is the runner torque, ω is the angular velocity, and Q is the water flow.

The operating point is determined based on the characteristic curve and similarity criterion of the pump turbine. As shown in Figure 4, this study provides insights into the working condition of the first quadrant of the pump turbine. Under the condition of rated head H = 40 m, the flow rate is adjusted by changing the guide vane opening so as to obtain different loads and outputs. As wind energy and solar energy are unstable and intermittent, it is necessary to operate the pumped-storage units under unbalanced load conditions to balance the power grid parameters. In order to compare the influence of wave disturbance on the operating status of the pump under different operating loads, this study selects three levels of guide vane openings, namely, small opening α₀ = 6°, medium opening α₀ = 15°, and rated opening α₀ = 24°. The corresponding flows are 0.3 Q_r, 0.7 Q_r, and rated flow Q_r, respectively.

For the wave disturbance, we select the characteristics of waves near an island in the East China Sea during the last 10 years as an example. To simplify the complexity of natural waves and quantify the effect of each influencing factor of natural waves, we assume the wave as a regular wave based on the characteristics of natural waves rather than using monitoring data of natural waves directly. According to the micro-amplitude wave theory, hydrostatic pressure and hydrodynamic pressure can be considered following sine or cosine functions. That is, the waveform basically follows a sine function, the maximum effective wave height is 10 m, and the wave period is 6 s (Yu et al., 2020). According to the relevant literature (Nanjing Hydraulic Research Institute, 2002), the pressure scale and frequency scale can be expressed as follows: λ p = λ k = l p l m , (14) λ f = λ k − 1 2 , (15) where λ k is the wavelength scale, l p is the wavelength of the wave in the real world, l m is the wavelength of the model, λ p is the pressure scale, and λ f is the frequency scale. In this study, the wavelength scale λ k is set to 10.

After a scale conversion between the modeled wave and real wave, the tail water level disturbance caused by wave disturbance in the pump turbine model is ∆ H = 1 × sin ⁡ ⁡ 2 π t T , and the pressure fluctuation at the outlet of the draft tube can be calculated using the formula ∆ P = 10000 ⁡ sin ⁡ ⁡ 10 π t 3 . When the wave period is 1.9 s, the wave frequency is 0.527 Hz. This pressure fluctuation is the boundary condition for the operation under the wave disturbance. The initial settings of the selected working conditions are shown in Table 1.

In this study, Fluent 2020 R2 was used in numerical simulation. The boundary conditions of the proposed model are as follows: the spiral case inlet is set as the mass flow inlet, and the draft tube outlet is set as the pressure outlet. Based on the coupling of pressure and velocity, the SIMPLEC algorithm is used for steady calculation, and the PISO algorithm is used for unsteady calculation. The second-order upwind format is used to discretize the momentum and turbulent kinetic energy dissipation rate. All over-current components are connected through an interface. The solid wall is set as a non-slip wall condition, and the residual convergence accuracy is set to 10⁻⁵.

We first assume that the fluid has a steady flow and the result of steady-flow calculation is the initial conditions of the model. The rotating speed of the runner is 760 r/min, and the fixed time step is set to 6.579×10⁻⁴ s (rotate 3° for each time step). The maximum number of iteration steps is set to 500, and the iterative convergence is stable in 50 steps. For the situation without wave disturbance, 15 rotor rotation periods are calculated. A total of 10 rotor rotation periods are calculated for analyzing the dynamics of fluid after being stable. Under the situation with the wave disturbance, two wave disturbance periods are calculated, and the results of the second wave disturbance period are analyzed.

In order to explore how the pressure pulsation varies with the flow field in the pump turbine, monitoring points are arranged every 45° in the circumferential direction of the axial cross section of the spiral case, every 30° in the vaneless space between the guide vane and the runner, and every 0.1 m in the center of the axial plane of the draft tube. The sketch of the positions of the monitoring points is shown in Figure 5.

Figure 6 shows the streamline in the guide vanes, runner, and draft tube under different operating conditions without wave disturbance. The velocity triangle corresponds to the position of the red dot in the streamline diagram, where U is the circumferential velocity, W is the relative velocity, and V is the absolute velocity. Figure 7 shows the vortex distribution in the draft tube under different operating conditions without wave disturbance.

As shown in Figure 6A and Figure 7A, under the operating condition of 0.3 Q_r, due to the small inlet flow and velocity, the inlet attack angle of the runner is large, and the water flow forms a high-speed water-retaining ring under the influence of the static and dynamic interference when it flows through the vaneless space. The generation of the high-speed water-retaining ring makes the fluid that should have flowed into the original runner rotate at high-speed flow into other runners. At this time, the circumferential velocity of the water flow on the inner wall of the draft tube is important, a rotating downward-flow state is formed on the wall of the draft tube, and a large amount of backflow is generated along the axial direction of the draft tube center. This flow state creates a high-intensity spiral vortex that almost fills the straight cone section of the draft tube. A high level of vortex in the draft tube may lead to hydraulic resonance. Therefore, the low-load operating condition should be avoided during the operation of the turbine.

As shown in Figure 6B and Figure 7B, under the operating condition of 0.7 Q_r, with the increase in the guide vane opening, the high-speed water ring in the vaneless space obviously weakens, and the maximum velocity decreases obviously, which is obviously improved compared with the low-flow operating condition. However, as the water flow at the inlet of the runner cannot maintain the best angle with the runner blades, the vortex and backflow remain on the water-facing surface of the runner blades. Although the water flow in the draft tube is inclined to the axial direction, it still maintains a large circumferential velocity component, and the flow field at the bottom of the straight cone section of the draft tube is chaotic. A significant vortex is observed, while the strength of the vortex is reduced.

As shown in Figure 6C and Figure 7C, under the operating load of Q_r, the water flow rate is increased, and the velocity of the fluid between the guide vane and the runner is uniform. There is no high-speed water ring in the vaneless space. In addition, the water enters the runner evenly and stably. The streamline in the draft tube is modeled with mainly axial flow, with a small circumferential velocity component. A small draft tube vortex area and low vortex strength are observed, which means that the unit runs smoothly under this operating condition.

Figure 8 shows pressure distribution in the pump turbine without wave disturbance under different operating conditions. The pressure in each working condition gradually decreases from the spiral case to the runner in the centripetal direction. In addition, the pressure gradient changes in the vaneless space. With the increase in the opening of the guide vane, the variation in pressure in the vaneless space gradually becomes stable. As the blades of the runner convert the kinetic energy of the water flow into the rotating mechanical energy of the runner, the pressure gradient in the runner is obvious, and the pressure on the water surface of the blade is generally higher than that on the back surface. The pressure on the waterfront surface of the blade is generally higher than that on the back surface, and the blade pressure gradually decreases from the inlet to the outlet, forming a negative pressure area at the blade outlet. The influence of the guide vane opening on the surface pressure of the blade is obviously less than that under the operating condition of a large opening. Under the operating condition of 0.3 Q_r, as the flow rate is low and the circumferential speed is high, the water flow clings to the draft tube wall in a rotating downward-flow state. The maximum pressure area can be observed on both sides of the draft tube inlet wall, and the straight tapered section of the draft tube presents an obvious conical pressure gradient. With the increase in the guide vane opening, the flow rate increases, the conical pressure gradient of the draft tube gradually stabilizes, and both sides of the draft tube inlet wall become the minimum pressure area.

Figure 9, Figure 10, and Figure 11 show the influence of wave disturbance on the internal flow pattern of the pump turbine under the operating condition of Q_r, allowing us to analyze the streamline, pressure, and vortex of the draft tube in channels with wave disturbance. The effect of wave disturbance on the water head is relatively small. Furthermore, the wave disturbance has no obvious influence on the flow direction, velocity, and vortex distribution of the draft tube. However, the pressure at different times is obviously different. The maximum pressure occurs at T/4 and the minimum pressure occurs at 3 T/4, which confirms the sinusoidal functions of wave disturbance. Therefore, the wave disturbance only affects the pressure in the pump turbine. Its influence on the streamline, draft tube vortex rope, and pressure distribution can be neglected.

To analyze the main characteristics of the pressure pulsation, two dimensionless numbers, the relative pressure pulsation amplitude Δ H ′ and the pressure pulsation coefficient C _p , are introduced to characterize the pressure pulsation (Xu et al., 2022). Δ H ′ = Δ H H = P i max − P i min ρ g H × 100 % , (16) C p = P i − P ¯ ρ g H × 100 % , (17) where ΔH is the peak value of pressure pulsation, H is the water head, P i m a x and P i m i n are the maximum and minimum pressures at point i, respectively, C _p is the dimensionless pressure pulsation coefficient, P i is the corresponding pressure at point i, and P ¯ is the average pressure over time.

Figure 12 shows the relative pressure pulsation at the monitoring points in the spiral case. For the condition without wave disturbance, the relative pressure pulsation amplitude increases with the increase in load, and the maximum value of Q_r is 1.40% under the rated working condition. For the condition with wave disturbance, the relative pressure pulsation amplitude of each working condition generally increases. However, the 0.7-Q_r working condition is affected by wave disturbance, and its relative pressure pulsation amplitude is almost the same as the rated working condition, which means that the relative pressure pulsation amplitude is related not only to the flow rate but also to the turbulent flow pattern.

Figure 13 shows the relative pressure pulsation amplitude at the monitoring points in the vaneless space. With the strong dynamic and static interference caused by the rotation of the runner blades, the relative pressure pulsation amplitude of the measuring points in the vaneless space is higher than that in other areas. For the condition without wave disturbance, the relative pressure pulsation amplitude increases with the increase in the load, and the maximum value is 16.55% under the rated working condition Q_r. For the condition with the wave disturbance, the relative pressure pulsation amplitude of each working condition generally increases, and the maximum value is 21.93% under the 0.7-Q_r working condition. However, due to the large relative pressure pulsation amplitude of the vaneless space under the large-flow working condition, the amplitude variation caused by wave disturbance is not significant under the rated working condition Q_r. Under the 0.7-Q_r working condition, part of the fluid flows with a high-speed water-retaining ring in the vaneless space, and the flow field is turbulent, which means that the pressure pulsation is obviously affected by wave disturbance.

Figures 14A, B show the relative pressure pulsation amplitude at the monitoring points in the draft tube for the conditions without and with wave disturbance, respectively. Compared with other areas, the relative pressure pulsation amplitude in this area is the smallest. For the situation without wave disturbance, the relative pressure pulsation amplitude in this area is related not only to the load but also to the strength of the eccentric vortex in the draft tube. The relative pressure pulsation amplitude has a positive correlation with the strength of the eccentric vortex in the draft tube, so the maximum relative amplitude of pressure pulsation is observed at the P4 measuring point with a value of 0.75% under the working condition of 0.7 Q_r. For the situation with wave disturbance, the relative amplitude of the pressure pulsation increases with the increase in load under different operating conditions, which shows that the eccentric vortex intensity of the draft tube is the main factor affecting the relative amplitude of the pressure pulsation for the case without wave disturbance.

As the frequency domain of the pressure pulsation at the monitoring points in the same area is highly similar, the representative monitoring points are selected for analyzing the frequency domain of the pressure pulsation, based on the internal flow characteristics of the pump turbine and the relative amplitude of the pressure pulsation. The W1 monitoring point is selected in the spiral case, and the Z1 monitoring point is selected in the vaneless space. Due to the influence of the vortex rope on the draft tube, the draft tube flow pattern changes are more complicated, and the relative amplitude of the pressure pulsation is more obviously affected by the wave disturbance. Three representative monitoring points, P1 at the draft tube inlet, P5 at the center of the straight cone section, and P9 at the bottom of the elbow section, are selected to analyze the pressure fluctuation in the frequency domain.

Figure 15 shows the relationship between the pressure pulsation factor and the frequency domain at monitoring point W1 in the spiral case area. f _n is the frequency conversion factor of the runner (approximately 12.7 Hz), and 0.04 f _n is the wave disturbance frequency (approximately 0.527 Hz). For the situation without wave disturbance, the pressure pulsation increases with the increase in load. This is because the frequency components are complicated and numerous, and there are more low-frequency and high-amplitude pressure pulsations under the operating condition of 0.7 Q_r (Case 2). For the situation with wave disturbance, the pressure pulsation has an additional wave disturbance frequency. However, the other main frequency characteristics did not change obviously. Thus, we only analyze the main frequency components under the wave disturbance (see Figure 15B). Here, the frequency of the wave disturbance is 0.04 f _n, which is caused by the wave pressure propagating from the downstream to the spiral case. Furthermore, the blade frequencies 7 f _n, 14 f _n, and 21 f _n and their multiples are mainly caused by static and dynamic interference. However, under different loads, the alternating hydraulic components acting on the runner blades are different, resulting in different amplitudes of pressure pulsation. In general, with the increase in the load, the frequency of the pressure pulsation increases. Moreover, 20 f _n is the number of guide vanes multiplied by the rotation frequency. This frequency type is mainly affected by the thickness of the guide vanes and the runner blades crowding out the water flow. Thus, the water flow pressure fluctuates periodically, which is particularly obvious under the light-load condition (Cases 1 and 4). Due to the high flow rate and the influence of the high-speed water-retaining ring on the vaneless space under the working condition of Cases 2 and 5, there are more vortices and secondary flows in the spiral case area, resulting in more high-amplitude pressure pulsations similar to 3.5 f _n.

Figure 16A shows the frequency domain results of pressure pulsation at monitoring point Z1 in the vaneless space, and Figure 16B shows the main frequency components of pressure pulsation under wave disturbance at this monitoring point. The pressure pulsation amplitude in the vaneless space is high, mostly due to the blade frequency and its frequency multiples, indicating that the pressure pulsation in the vaneless space is mainly caused by static and dynamic interference, and the higher the load, the higher the pressure pulsation amplitude caused by the static and dynamic interference. At 0.7 Q_r (Cases 2 and 5), there are many low-frequency and high-amplitude pressure pulsations; however, the other conditions do not appear, indicating that the occurrence of low-frequency and high-amplitude pressure pulsations in the vaneless space is related to the chaotic flow pattern. Due to the large pressure pulsation caused by the static and dynamic interference in the vaneless space, the main frequencies of the other two working conditions have not changed under the wave disturbance, except for the low load 0.3 Q_r (Cases 1 and 4), which is still 7 f _n.

Figure 17 shows the variation in the pressure pulsation coefficient to the frequency at the monitoring points P1, P5, and P9 in the draft tube. Due to the low-frequency and high-amplitude characteristics of the pressure pulsation in the draft tube, the pressure pulsation under wave disturbance is superimposed with the original pressure pulsation in the draft tube, resulting in more low-frequency and high-amplitude pressure pulsation. Therefore, the main frequency components for the two cases with and without wave disturbance are analyzed, respectively.

Figures 18A, C, E show the main frequency components of pressure pulsation in the draft tube without wave disturbance. There is a large amplitude of pressure pulsation of 0.11 f _n for all operating conditions. As a similar unsteady pressure phenomenon, many scholars have confirmed that the main frequency of the draft tube pressure pulsation of the pump turbine is approximately 0.11 f _n. Under the working conditions of 0.3 Q_r (Case 1) and 0.7 Q_r (Case 2), the upstream monitoring points P1 and P5 have a maximum amplitude of pressure pulsation at 0.23 f _n. However, under the working conditions of Q_r (Case 3), the pressure amplitude corresponds to 0.23 f _n, and the rotation frequency of the vortex rope is generally 1/3–1/5 of the rotation frequency, so 0.23 f _n is the main frequency of the draft tube vortex rope. The runner rotating frequency is 1 f _n, which is caused by the hydraulic imbalance in the runner components. This pressure pulsation has an impact on the whole draft tube. The closer to the runner, the higher the amplitude of pressure pulsation of 1 f _n, and the frequency of the maximum pressure amplitude at P1 point at the draft tube inlet is 1 f _n. For the frequency of 2 f _n, this type of frequency is formed by the vortex and secondary flow passing through the runner, and its corresponding amplitude of pressure pulsation is small.

Figures 18B, D, F show the main frequency components of pressure pulsation for the condition with wave disturbance. Under the working conditions of 0.3 Q_r (Case 4) and 0.7 Q_r (Case 5), the main frequency of the draft tube vortex rope is 0.23 f _n, and the pressure pulsation amplitude is still high. As the monitoring point P1 is close to the runner, it is influenced by the dynamic and static interference of the runner. Hence, it has a high-pressure amplitude at the rotating frequency of 1 f _n. Compared with the condition without wave disturbance, the main frequency of the pressure pulsation at each monitoring point is close to 0.527 Hz, which is the frequency of the wave disturbance. In addition, there are several high-amplitude pressure pulsations with a frequency doubling of 0.527 Hz, such as 0.08 f _n. In this case, the frequency is superimposed with the original low-frequency pressure pulsation of the draft tube. Therefore, the frequency spectrum of the pressure pulsation becomes complicated. The amplitude of the high-frequency pressure pulsation above 1 f _n shows no obvious variation, but the pressure amplitudes of 0.11 f _n and 0.23 f _n generally increase by 2–15 times. The wave disturbance has an impact on the low-frequency pressure pulsation of the draft tube.