Investigate the full characteristic of a centrifugal pump-as-turbine (PAT) in turbine and reverse pump modes

Abstract

Large centrifugal pumps can generate electricity in reverse, reducing the power supply pressure of the power grid. However, S region may occurs during the operation and regulation of the pump as a steam turbine (PAT). The S region can cause unstable unit startup, not only causing hydraulic system oscillation, but also accompanied by strong vibration and noise. Understanding the characteristic curve of the S region is of great significance for improving unit stability and efficiency. Previously, the main research object in S area was small generator units. This article conducted unsteady numerical simulation and experimental research on a large centrifugal pump. Analysed the internal flow rate, blade load, entropy generation, and pressure pulsation under three special operating conditions: the optimal operating condition, runaway operating condition, and braking operating condition in the S region. In addition, a comparative analysis was conducted on three pairs of operating conditions: S region turbine mode and reverse pump mode, with opposite flow direction, similar flow rate, and the same speed. Research has found that in the S region of large centrifugal pumps, the flow chaos in the turbine mode mainly occurs in the impeller section, and the maximum pressure occurs in the stationary blade section. The flow chaos under reverse pump operation mainly occurs in the stationary blade section, and the maximum pressure occurs in the impeller section.

KEYWORDS:

• Pump as turbine
• S region
• turbine mode
• reverse pump mode
• entropy production

1. Introduction

With the rapid improvement of scientific level, the proportion of renewable energy in the world energy structure is becoming higher and higher (Amiri et al., 2016; Varun et al., 2009). Due to the intermittence and randomness of wind and light energy (Bazilian et al., 2012; Favrel et al., 2018), the large scale grid connection of power plants has led to increased operational pressure on the power grid. Existing large centrifugal pump units can generate electricity by reversing operation of turbine, which can provide better feasibility and economy for peak shaving and valley filling of power grid (Paish, 2002). Compared with direct use of turbine, PAT has the advantages of extensive application and mature manufacturing process. In addition, there are more types of pumps with greater working flow and head ranges than hydraulic turbines (Sinagra et al., 2017; Williams, 1996). There has been a large amount of theoretical, experimental and numerical research on PAT. The early investigation attempts of PAT began in the 1930s. In recent decades, it has attracted much attention. Raman's experimental studies have demonstrated that pump operation in turbine mode requires a higher head and flow rate than in pump mode. However, centrifugal pumps can operate just as well as water turbines without any mechanical problems (Raman et al., 2013). It is found that the optimum efficiency of the pump in reverse operation is not lower than that in positive rotation, and the head and flow rate of the optimum efficiency point are higher (Nautiyal et al., 2010). Orchard and Klos (2009) provided the application and advantages of PAT and elaborated on its application in the water industry (Orchard & Klos, 2009). Huang et al. (2017) carried out flow analysis according to the matching principle between rotor and volute (Huang et al., 2017). Performance prediction methods considering model loss have been gradually used to explore PAT performance in recent years. This method has been successful in pumps (Derakhshan & Nourbakhsh, 2008; Liu et al., 2019). In addition, the polynomial fitting method is also used to calculate the performance of PAT (Fecarotta et al., 2016; Jain et al., 2015).

There is hump problem in pumping mode of pump turbine (Liu et al., 2021), which will have S region problem in turbine mode. Similarly, the same problem exists in PAT (Amblard et al., 1985). However, PAT and pumped storage units have the same problem. The potential dangers of S characteristics include: (1) When the unit operates at a small opening, the S characteristics can cause instability in the unit startup and make it difficult to connect to the grid (Hu et al., 2019); (2) When the unit operates at large opening, S characteristic curve will cause severe water hammer pressure and serious vibration and pressure pulsation during load rejection, which will seriously affect the operation stability of the unit and may cause safety problems (Hu et al., 2018).

There are already many research foundations and developments for the S characteristics research of pump turbine. In this process, the research method is based on the evolution from prototype observation to model testing, from one-dimensional simulation to three-dimensional computational fluid dynamics simulation. The research scale is from global to local, and the research object is from system to individual. With the continuous improvement of turbulence model in engineering application and the continuous development of science and technology, the internal flow mechanism of S characteristic of pump turbine is gradually revealed. The relevant research progress can be divided into two parts according to the technical route: model test and numerical simulation. Among them, model experimental research needs to rely on advanced testing methods and high-precision model experimental platforms to achieve flow regime observation, velocity measurement, and pressure measurement. Hasmatuchi et al. (2011) based on a static model test rig, visualised the rotating stall phenomenon of the pump turbine in the S region by injecting bubbles, and quantitatively analysed the pressure pulsation characteristics of measuring points such as the vaneless region and the volute inlet (Hasmatuchi et al., 2011). Zeng et al. (2017) discussed the impact of S characteristics on runaway stability and how to improve the water hammer pressure under load rejection conditions by controlling the guide vane closing law (Zeng et al., 2017). Walseth et al. (2016) tested low specific speed pump turbine in runaway condition and predicted the dynamic characteristics of the unit after runaway condition through one dimensional transition process calculation (Walseth et al., 2016). Botero et al. (2014) proposed a non-intrusive diagnostic method for rotating stall based on the stationary pulsating signal with certain periodicity (Botero et al., 2014). This method only needs to measure the vibration signals of guide vane, which has the advantages of simple operation and applicability to practical power plants. As a supplement to the experiment, CFD numerical simulation can help to obtain the details of the flow field inside the unit, thereby establishing the relationship between the instability of the unit's external characteristics and the flow characteristics. The application of CFD to the study of S characteristics includes flow patterns, inlet vortices, pressure pulsation propagation, rotor force distribution, and rotational stall propagation. Through detailed calculations, Xia et al. (2017) pointed out that the counter-current vortical flow structure at the runner inlet of the pump turbine is related to the instability when it operates in S region (Xia et al., 2017). Lenarcic and Gehrer (2019) explained the characteristics of S region through numerical simulation and model experiment, and proposed methods to improve the operation stability of unit such as asynchronous guide vane, variable speed operation and reducing static pressure at runner outlet. Jacquet et al. (2016) simulated the external characteristics of model pump turbine at medium specific speed and constant opening, and analysed the flow velocity distribution and pressure fluctuation at the runner inlet (Jacquet et al., 2016). The results showed that using the SAS-SST turbulence model can obtain relatively close results to model tests. He also pointed out that the S characteristic curve is related to the backflow occurring in the vaneless region, which is caused by the flow separation at the blade inlet side. Based on numerical simulation, Liu et al. (2016) quantitatively analysed the reflux phenomenon of pump turbine under off-design conditions, and found that guide vane opening has a significant effect on the spatial distribution of reflux area, and there are differences in the spatial distribution of the backflow phenomenon under the runaway point and hydraulic turbine braking conditions (Liu et al., 2016). Zhang et al. (2021) investigated the hydrodynamic characteristic of a single-stage centrifugal pump with inlet inducer and outlet Radial Guided Vanes (RGVs) influenced by the clocking effect for the first time (Zhang et al., 2021). Staubli et al. (2010) demonstrated through numerical simulation that when the runner operates near the runaway point, there will be partial reverse flow back to the flow area in the runner flow path, and these fluctuating flows can account for approximately 50% of the mainstream flow (Staubli et al., 2010). He also pointed out that this flow instability is caused by local vortices at the inlet side of the blade flow path. Aiming at the phenomenon of flow vortices and draft tube vortex bands in the flow passage of a pump turbine that deviates from the design conditions, Kim et al. (2019) pointed out through numerical simulation that these vortex occurrence regions have significant blockage effects and reflux phenomena, and there are also significant pressure fluctuations in the regions where the vortices are relatively significant (Kim et al., 2019). Widmer et al. (2011) found that stationary vortex formation and rotating stall have initially the same physical cause, but it depends on the mean convective acceleration within the guide vane channels, whether the vortex formations will rotate or not. Lu et al. (2021) studied the pressure fluctuation characteristics of a pump turbine under runaway operating conditions in turbine mode and found that the characteristic frequency of pressure fluctuation in the flow channel is caused by the interaction between vortices moving in the bladeless and runner regions and the runner blades. According to the experimental results, Zhang et al. (2020) found that in turbine mode, guide vane opening had the greatest influence on pressure fluctuation. This effect is weaker in turbo braking mode and even weaker in reverse pump mode. Flow separation and pressure fluctuation under different flow rates of turbine mode and reverse pump mode are compared. It is found that pressure fluctuation in bladeless space is larger in turbine braking mode. At a low flow rate, the separation in the leafless space becomes very obvious, while the separation in the flow passage weakens, resulting in the decrease of pressure fluctuation in the leafless space. Under reverse pumping condition, the pressure fluctuation in the bladeless space is relatively small due to the obvious symmetrical separation in the flow passage. Y. Zhang and Wu (2017) pointed out that in generating mode, rotating stall initiates at runaway and is fully developed at low discharges in turbine brake mode. At deep reverse pump mode, rotating stall also exists but is not a dominant mechanism.

It can be seen that on the Q₁₁-n₁₁ curve, the conventional turbine guide vane opening line slightly bends downward at a high unit speed, and the included angle formed by the runaway line (T = 0) is large, as shown in the dashed line in Figure 1 (Zuo et al., 2016). Generally, the diameter of the runner of a pump turbine or a PAT is larger than that of a conventional turbine, which makes the centrifugal force acting on the water flow in the runner greater. The guide vane opening line exhibits a significant curvature at high unit speeds, and the angle formed by the runaway line is smaller, as shown by the solid line in Figure 1. The slope of the guide vane opening line of the pump turbine in the turbine braking area (The direction of rotation is the same as the turbine mode and the flow rate is positive) is greater than zero, indicating that when operating in this area, the unit flow Q₁₁ decreases, and the unit rotational speed n₁₁ also decreases, mainly because the water flow in the braking area hinders the rotation of the runner. When the unit flow rate is close to 0, the rotational speed is still relatively high, and the water flow may flow out in the opposite direction due to the greater centrifugal force when entering the runner. At this point, the pump turbine enters the reverse pump mode (an area where the direction of rotation is the same as the turbine mode and the flow rate is positive negative). In reverse pump mode, the flow rate increases as the speed increases. In general, the guide vane opening line presents an inverse ‘S’ region in turbine mode and reverse pump area. At this time, a single unit speed corresponds to several units of flow. This means that when the pump turbine operates in these areas, the operating point may jump, switching from the turbine operating condition to an unstable turbine braking condition or even a reverse pump operating condition. When the head is low and the governor is not involved in automatic adjustment of the unit's no-load opening, The switching of unit operating conditions is represented by the oscillations of unit speed. The lower the head, the smaller the speed at which oscillation occurs, the greater the amplitude of oscillation, and the larger the range of entering the unstable zone; When the governor participates in automatic regulation, the hydraulic inertia of the unit affects the interaction between the governor and the S characteristic, making it difficult for the unit to connect to the grid or unable to achieve no-load stability or even trip after load rejection. In addition, when the unit is connected to the grid, if it passes through the S characteristic area or its vicinity, in order to maintain the synchronous speed of the unit and the power grid unchanged, the unit will absorb a large amount of power from the system. If the absorbed power exceeds the specified value, it will cause the reverse power protection action of the unit to trip. In addition, due to the influence of the runner S characteristics, the water flow in the pressure steel pipe will cause the same effect as the guide vane closing quickly, which leads to the rapid reduction of the flow into the unit, thus forming the peak pressure of the water hammer in the pressure steel pipe. Therefore, studying the S region characteristic curve to reduce the water hammer pressure and pulsating pressure during the transition process has always been the goal of researchers from domestic and foreign water turbine research and development teams and research institutions. Due to the rarity of reverse power generation in large pump units, this article investigates the energy loss, blade load, and pressure pulsation in the reverse S region of a large centrifugal pump, laying a theoretical foundation for the stable and efficient operation of reverse power generation in large vertical centrifugal pumps in the future.

Figure 1. S characteristic curve of pump turbine.

2. Research object and calculation model

2.1. Research object

The basic parameters of the unit are shown in the Table 1 and the structural diagram is shown in Figure 2.

Figure 2. (a) Structural drawing of large vaned-voluted centrifugal pump (b) Unit impeller inlet and outlet diagram.

Table 1. Basic parameters of the unit.

Name	Parameter
Number of impeller blades	9
Number of stay vanes	12
Impeller inlet diameter	0.28m
Impeller outlet diameter	0.56m
Rated speed (pump mode)	1000r/min
Q of BEP point (pump mode)	0.237 kg/s
H of BEP point (pump mode)	47.33m
η of BEP point (pump mode)	90.1%

2.2. Model test

The schematic of the model test rig is shown in Figure 3. Table 2 shows the main parameters of the test rig. Experimental measurements were carried out using the closed loop water circuit to get realistic conditions. Water from the cavitation tank is pumped to the reservation tank and sequentially flow down to the flow stabilisation tanks. The test bench site is shown in Figure 4, and test instrument parameters are listed in Table 3. Before the experiment, we calibrated various parameter sensors to ensure the accuracy of the experimental measurement. All instruments and metres were calibrated using standard gauges with valid qualification certificates, or directly sent to the national metrology department for regular verification. In the actual experimental process, we take the average of multiple measurements to minimise random errors. The experimental PXI system is used for data acquisition and testing, it can provide high-performance modular instrumentation and a rich set of I/Os with dedicated synchronisation and major software functions for fluid machinery measurement applications. The PXI system acquires measurement data from various sensors and then transmits it to the computer for processing. The flow of the inlet and outlet is divided by the section area to obtain the average flow velocity thus the dynamic pressure head. The flow measurement uses an inner diameter of 300 mm electromagnetic flowmeter. The electromagnetic flowmeter uses standard metre method for on line calibration. The speed is measured by an encoder connected to the generator shaft, which is installed directly on the shaft end of the dynamometer motor. When the unit rotational speed is in to steady state, the measured rotational speed frequency is calibrated with a calibrated high precision frequency metre, or the speed is verified with a stroboscope. Energy test is conducted to determine the relationship between head H (m), unit speed n₁₁, unit flow rate Q₁₁ and unit torque M₁₁ under various operating conditions. Several operating conditions located in the S zone are selected to test on the model machine. During the test, the fluctuation value of the rotational speed is Δn≤±0.2%. Each parameter is calculated according to the following formula, the meanings of each symbol in the formula are shown in the Table 4:

• (1) Head H (m) (1) $\begin{array}{r}H=\frac{\mathrm{\Delta }P}{\rho g}+\frac{V_s^2-V_d^2}{2g}\end{array}$(1)

The reverse pump mode can be replaced by the inlet and outlet in the formula

• (2) Unit speed n₁₁ (r/min) (2) $\begin{array}{r}{n}_{11}=\frac{n\times D}{H^{0.5}}\end{array}$(2)

• (3) Unit flow Q₁₁ (m³/s) (3) $\begin{array}{r}{Q}_{11}=\frac{Q}{D^2\times H^{0.5}}\end{array}$(3)

• (4) Unit torque M₁₁ (N·m) (4) $\begin{array}{r}{M}_{11}=\frac{M}{D^3\times H}\end{array}$(4)

2.3. Numerical simulation

At present, there are three methods to deal with turbulence numerical calculation: direct numerical simulation (DNS) method, large eddy simulation (LES) method and Reynolds mean N-S equation (RANS) method. In summary, because of the large size of the unit, LES or DNS method consumes too much resources and is difficult to calculate. The core of Reynolds average method is not to solve the instantaneous Navier-Stokes equation directly, but to find a way to solve the time-homogenized Reynolds equation. This can not only avoid the problem of large calculation, but also achieve good results for practical application of engineering. Therefore, this paper chooses RANS method which is more suitable for engineering calculation. Due to the different starting points of model processing, turbulence model theory can be divided into two categories: Reynolds stress model and vortex viscous closed model (abbreviated as vortex viscous model). Due to the limitation of calculation conditions, the Reynolds stress model has a large amount of computation and its application scope is limited. Therefore, the vortex viscosity model is widely used in engineering turbulence problems. The vortex viscosity model was proposed by Boussinesq (1877) imitating the idea of molecular viscosity, namely, let Reynolds stress be: (5) $\begin{array}{rl}{v}_t& =C_{\mu }{k}^2/ϵ,v_t=C_{\mu }\frac{k}{\omega },v_t=C_{\mu }k\tau ,\\ v_t& =C_{\mu }\frac{q^2}{\omega },v_t=C_{\mu }\sqrt{kl}\end{array}$(5) (6) $\begin{array}{rl}\overline{u_i{u}_j}& =-v_t\left(U_{i,j}+U_{j,i}+\frac{2}{3}{U}_{k,k}{\delta }_{ij}\right)+\frac{2}{3}k{\delta }_{ij}\end{array}$(6) where is the turbulent kinetic energy, v_t is called the eddy viscosity coefficient. This is the first reference eddy viscosity mode proposed, which assumes that Reynolds stress is linear with the average velocity strain rate. When the average velocity strain rate is determined, the six reynolds stresses can be completely determined by determining only one eddy viscosity coefficient v_t. different eddy viscosity modes can be obtained, such as ‘k-ε’, ‘k-ω’, ‘k-τ’, ‘q-ω’ and ‘k-l’. In order to make the governing equation closed, how many additional turbulent quantities are introduced, how many additional differential equations must be solved at the same time. According to the number of additional differential distances required to be solved, the vortex viscous modes can be generally divided into three categories: zero-equation mode and half-equation mode, one-equation mode and two-equation mode. The ‘k-ω’ model chosen in this paper is one of the most widely used turbulence models.

Figure 3. Schematic diagram of test device.

Figure 4. Experimental device.

Table 2. Parameters of the test device.

Maximum test head	120 m
Dynamometer motor power	600 kW
Maximum speed of dynamometer motor	1900 r/min
Maximum flow of test bench	1.20 m^3/s
Tail water pressure	−85.0∼250.0 kPa
Reference pressure (absolute pressure)	0 Pa

Table 3. Test instrument parameters.

Number	Measured physical quantity	Name	Accuracy
1	Flow	Electromagnetic flowmeter	0.18%
2	Rotation speed	Rotary encoder	0.02%
3	Head	Differential pressure sensor	0.05%
4	Moment	Load cell	0.05%
5	Moment	Excitation amplifier	0.02%
6	Pressure pulsation	Dynamic pressure sensor	1%
7	Pressure pulsation	PCB signal amplifier	1%
8	Atmospheric pressure	Digital pressure gauge	0.1%

Table 4. List of symbols.

Parameter	Meaning
ΔP	Static pressure difference between pump inlet and outlet pressure measuring sections, measured by differential pressure sensor (Pa)
Q	Test flow (m^3/s)
Ad	Outlet area (m)
As	Inlet area (m)
Vd	Average flow velocity at the outlet pressure measuring section (m/s), Vd  = Q/Ad
Vs	Average flow velocity of inlet pressure measuring section (m/s), Vs  = Q/As
ρ	Test water density (kg/m^3)
G	Acceleration of gravity (m/s^2)
N	Impeller speed (r/min)
D	Pump impeller outlet diameter (m)
H	head (m)
M	Shaft torque (N·m)

This model greatly improves the treatment of low Reynolds number flows in the near-wall region. It no longer requires the construction of nonlinear attenuation function, and the requirement of the first layer of grid is no longer located in the log-law layer. Therefore, the applicability of turbulence model is improved, and flow separation phenomena with adverse pressure gradient can be predicted more ideally. Currently available ‘k-ω’ models mainly include Wilcox ‘k-ω’ model, Baseline ‘k-ω’ model and SST ‘k-ω’ model. The following will mainly introduce SST ‘k-ω’ model:

Dr. Menter from NASA Ames Research Center proposed the famous SST ‘k-ω’ in the turbulence model, a cross diffusion term has been added after the omega equation, which happens to be the difference between the ‘k-ω’ and ‘k-ε’ models (Menter, 1992). Dr. Menter cleverly multiplies a mixing function on this cross diffusion term, which can easily control the transformation of the turbulence model: when the mixing function is 1, the omega equation is still the omega equation, which is used for solving near the boundary layer, while when the mixing function is 0, the omega equation is transformed into the epsilon equation, which is used for solving the mainstream, of course, a mixing zone is used between the two Domain for a smoother transition.

In the near wall region, the value of the mixing function is equal to 1, which is equivalent to the standard ‘k-ω’ Model.

In the far wall region, the value of the mixing function is equal to 0, which is equivalent to the standard ‘k-ε’ Model. (7) $\begin{array}{rl}{v}_t& =\frac{\rho a_{1}k}{max\left(a_1,\omega ,s{F}_2\right)}\end{array}$(7) (8) $\begin{array}{rl}{F}_2& =\tan \left({\mathrm{\Phi }}_2^2\right)\end{array}$(8) (9) $\begin{array}{rl}{\mathrm{\Phi }}_1& =max\left[2\frac{\sqrt{k}}{{\text{β}}^{\mathrm{\prime }}\omega y},\frac{500v}{y^{2\omega }}\right]\end{array}$(9)

In the formula: a₁ is the empirical coefficient; F₂ is a mixed function.

The overall calculation domain is shown in Figure 5 and the calculation settings are shown in the Table 5. In turbine mode, the volute is set as the mass flow inlet, and the draft tube is set as the average static pressure outlet, with the relative pressure of 0 Pa. The grid convergence is checked by using GCI criterion based on Richardson extrapolation method, parameter φ takes the head H of the unit, as shown in Table 6. The final grid number is 7.63 million, as shown in Table 7. The unsteady simulation calculation is carried out on the basis of 1000 step steady simulation. The impeller is iterated 180 steps per revolution, a total of 20 cycles are calculated, and the total iteration steps are 4600 steps.

Figure 5. Calculation domain.

Table 5. Numerical simulation setting.

Name	Set
Turbulence model	SST k-ω model
Frame Change/Mixing Model	Transient Rotor Stator
Medium	Liquid water at 25℃
Turbulence intensity	5%
Convergence criteria	RMS Residuals <1.0×10^−5
Grid Type	Hexa

Table 6. Grid independence check.

N1	N2	N3	φ1	φ2	φ3	GCI32fine	GCI21fine
7630000	4080000	1920000	33.40	33.93	33.53	3.08%	5.17%

Table 7. Number of grids.

Name	Number of grids	Number of nodes
Stay vane	3370204	584953
Impeller	4030277	699111
Volute	160591	28836
Inlet	71335	12508
Total	7,632,407	1,325,408

3. S region characteristic curve

3.1. Calculation results and errors

After unsteady calculation, the S region curve is obtained as shown in Figures 6 and 7. The naming and data of the calculated 22 operating points are shown in Table 8. Since the only variable in both n₁₁ and Q₁₁ is H, the error bar of n₁₁ drawn in Figure 8 also represents the error level of Q₁₁. It can be seen from the figure that the test data is in good agreement with the simulated data.

Figure 6. n₁₁–Q₁₁ curve of S region.

Figure 7. n₁₁–M₁₁ curve of S region.

Figure 8. Error bar of n₁₁.

Table 8. Simulated calculation of operating point data.

Operating Point	Q(m^3/s)	H(m)	n(r/min)	n11	Q11	M11
O3	−0.05323	12.1718	650.2	52.18283	−194.609	−1123.51
O2	−0.03754	13.6267	650.17	49.3162	−129.713	−757.174
O1	−0.01394	14.4173	650.17	47.94496	−46.8279	−431.624
A1	0.00289	13.8584	650.17	48.9022	9.902053	−339.252
A2	0.0087	13.893	650.27	48.84878	29.77181	−293.287
A3	0.01816	13.6007	650.27	49.37091	62.80861	−244.811
A4	0.02849	13.4692	650.19	49.60522	99.01602	−192.832
A5	0.0353	13.4692	650.21	49.60675	122.6839	−149.919
A6	0.0436	13.2438	650.17	50.02403	152.8143	−96.7633
A7	0.0525	13.0484	650.24	50.40262	185.3807	−38.2791
A8	0.06051	13.0461	650.26	50.40861	213.6834	−16.5976
A9	0.06378	13.0409	650.27	50.41943	225.2759	41.98358
A10	0.08588	13.3235	650.09	49.86805	300.1006	215.8843
A11	0.11314	14.0439	649.6	48.53558	385.0847	461.0633
A12	0.13168	14.8117	649.62	47.26231	436.4166	645.8628
B13	0.14188	14.1181	599.59	44.68114	481.6338	889.7474
C14	0.14874	13.9063	549.67	41.27188	508.7518	1052.176
D15	0.15312	13.4669	499.82	38.13624	532.2088	1199.677
E16	0.15739	12.8116	399.91	31.28374	560.8664	1380.03
F17	0.15947	13.4307	299.9	22.91317	555.0264	1482.854
G18	0.15952	14.1577	199.97	14.88081	540.7577	1542.187
H19	0.15759	13.9718	99.88	7.481874	537.7574	1615.029

3.2. Analysis of pressure fluctuation in S region curve

Based on the overall flow characteristics of the S region, we analysed the maximum amplitude values at different measuring points and different operating conditions. Through unsteady calculation of 22 operating conditions of the pump S region, we analysed the maximum pressure conditions of 22 operating conditions. After FFT changes, we obtained a pressure fluctuation nephogram, with the horizontal axis f_z defined as the blade passage frequency (97.5 Hz) at 650 r/min, f_z= n*Z/60, n is the speed (r/min) and Z is the number of impeller blades. In addition, the pressure fluctuation coefficient C_p is defined to represent the magnitude of the pressure fluctuation: (10) $\begin{array}{rl}{C}_p& =\frac{\mathrm{\Delta }p}{\frac{1}{2}\rho u_2^2}\end{array}$(10) (11) $\begin{array}{rl}{u}_2& =\frac{n\pi {\text{r}}_2}{30}\end{array}$(11) Δp is the detection pressure minus the average pressure; ρ is the density of the fluid; u₂ is the circumferential velocity at the impeller outlet, r₂ is the impeller outlet radius, and n is the impeller rotational speed. The overall results are shown in Figure 9:

Figure 9. Pressure analysis diagram for different conditions at same location.

From the pressure fluctuation coefficient of each flow area, the highest pressure conditions are H19 and G18 operating conditions of turbine mode. This may be due to the small unit speed of H19 and G18 operating conditions and the high unit flow resulting in the increase of pressure caused by liquid impinging on the overflow components. A11 and A12 of turbine mode have higher pressure, which may be caused by high flow rate and high speed. The pressure of O1, O2 and O3 of the reverse pump mode is also slightly higher, which may be caused by mismatch of liquid flow direction and impeller rotation direction, indicating that the operation stability of the reverse water pump is poor. From the frequency domain diagram of each flow area, it can be seen that the main frequency of most operating conditions from A1 to 12 in the stay vane and volute basically occurs at the position of one time the blade frequency. As the speed decreases from 650r/min, the main frequency position of other low-speed operating conditions gradually approaches 0. In the impeller outlet and bladeless section, the main frequency basically appears around 1.5 times the blade frequency, and the position of the main frequency is not obvious at low speeds. Among them, A6 and A9 are relatively special, with their main frequencies approaching 0.5 times the blade frequency. In summary, G18 and H19 have the highest pressure, but their maximum amplitude is relatively small. The pressure and pulsation amplitudes of turbine modes A11, A12, and reverse pump modes O1, O2, O3 are at high levels. It indicates that when the unit speed Q₁₁and unit flow n₁₁ of the unit are both high or low at the same time, the unit pressure will increase significantly and appear unstable state. Therefore, during the actual operation of the unit, the speed and flow rate should be avoided from being too high or too low at the same time to prevent vibration and other problems caused by too high pressure pulsation.

3.3. Flow analysis under special operating conditions of S region

After a certain analysis of the overall stability of S region, the optimal operating point E16, runaway operating point A9 and braking operating point A1 are extracted from the 22 calculated points for analysis (as shown in Figure 10 and Table 9). In numerical simulation, the Entropy Production Rate (EPR) can be used as a method to characterise the Flow Energy Dissipation (FED), which can be divided into E_loss and E_disp. When using the RANS method, it can be calculated as time-average term and pulsation term (Herwig & Kock, 2006): (12) $\begin{array}{rl}{E}^{\mathrm{\prime }}& =E_{\text{loss}}+E_{\text{disp}}\\ & =T\left[\right(S_{\overline{PC}})+(S_{P{C}^{\mathrm{\prime }}}\left)\right]+T\left[\right(S_{\overline{PD}})+(S_{P{D}^{\mathrm{\prime }}}\left)\right]\end{array}$(12) where $S_{\overline{PC}}$ and $S_{P{C}^{\mathrm{\prime }}}$ are the subterm of entropy production caused by loss term, $S_{\overline{PD}}$ and $S_{P{D}^{\mathrm{\prime }}}$ are the subterm of entropy production caused by dissipation term, T is the initial temperature: (13) $\begin{array}{rl}{S}_{\overline{PC}}& =\frac{{\lambda }_t}{{\overline{T}}^2}\left[{\left(\frac{\mathrm{\partial }\overline{T}}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{T}}{\mathrm{\partial }y}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{T}}{\mathrm{\partial }z}\right)}^2\right]\end{array}$(13) (14) $\begin{array}{rl}{S}_{P{C}^{\mathrm{\prime }}}& =\frac{{\lambda }_t}{{\overline{T}}^2}\left[\overline{{\left(\frac{\mathrm{\partial }{T}^{\mathrm{\prime }}}{\mathrm{\partial }x}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{T}^{\mathrm{\prime }}}{\mathrm{\partial }y}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{T}^{\mathrm{\prime }}}{\mathrm{\partial }z}\right)}^2}\right]\end{array}$(14) (15) $\begin{array}{rl}{S}_{\overline{PD}}& =\frac{\mu }{\overline{T}}[2\left\{{\left(\frac{\mathrm{\partial }\overline{u}}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{v}}{\mathrm{\partial }y}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{w}}{\mathrm{\partial }z}\right)}^2\right\}\\ & +{\left(\frac{\mathrm{\partial }\overline{u}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\overline{v}}{\mathrm{\partial }x}\right)}^2+{\left(\frac{\mathrm{\partial }\overline{u}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }\overline{w}}{\mathrm{\partial }x}\right)}^2\\ & +{\left(\frac{\mathrm{\partial }\overline{v}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }\overline{w}}{\mathrm{\partial }y}\right)}^2]\end{array}$(15) (16) $\begin{array}{rl}{S}_{P{D}^{\mathrm{\prime }}}& =\frac{\mu }{\overline{T}}[2\left\{\overline{{\left(\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }x}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{v}^{\mathrm{\prime }}}{\mathrm{\partial }y}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{w}^{\mathrm{\prime }}}{\mathrm{\partial }z}\right)}^2}\right\}\\ & +\overline{{\left(\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{v}^{\mathrm{\prime }}}{\mathrm{\partial }x}\right)}^2}+\overline{{\left(\frac{\mathrm{\partial }{u}^{\mathrm{\prime }}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }{w}^{\mathrm{\prime }}}{\mathrm{\partial }x}\right)}^2}\\ & +\overline{{\left(\frac{\mathrm{\partial }{v}^{\mathrm{\prime }}}{\mathrm{\partial }z}+\frac{\mathrm{\partial }{w}^{\mathrm{\prime }}}{\mathrm{\partial }y}\right)}^2}]\end{array}$(16) where x, y, z are coordinate components. The pulsation term $S_{P{C}^{\mathrm{\prime }}}$ can be empirically approximated as $S_{P{C}^{\mathrm{\prime }}}=A_t{S}_{\overline{PC}}$. $A_t$ is an empirical constant equal to 100. The pulsation term $S_{P{D}^{\mathrm{\prime }}}$ can be empirically approximated as $S_{P{D}^{\mathrm{\prime }}}=\rho ϵ/\overline{T}$, where $ϵ$ is the turbulence dissipation rate. Due to the small temperature changes during operation, the proportions of $S_{\overline{PC}}$ and $S_{P{C}^{\mathrm{\prime }}}$ are very small, so they are ignored. In addition, the proportion of $S_{\overline{PD}}$ on the wall is very large, while the proportion in the flow channel is very small. Since we mainly consider the energy changes in the flow channel, after removing the distribution of $S_{\overline{PD}}$ on the wall, we only consider the influence of the $S_{P{D}^{\mathrm{\prime }}}$ part, which accounts for a very large proportion. In order to better compare the distribution and size of $E^{\mathrm{\prime }}$, we dimensionless it as ‘e’, as shown in the formula, Q_r is the flow rate, P_r is the shaft power (Li et al., 2022). (17) $\begin{array}{r}e=E^{\mathrm{\prime }}\frac{Q_r}{P_r}\sqrt{\frac{H_r}{\text{g}}}\end{array}$(17) In addition, the dimensionless pressure parameter $C_r$ and velocity dimensionless coefficient $v^{\ast }$ are defined as shown in the formula: (18) $\begin{array}{r}{C}_r=\frac{p-p_{ref}}{\rho gH}\end{array}$(18) where p is the pressure and p_ref is the pressure at the inlet/outlet boundary of the draft tube. (19) $\begin{array}{r}{v}^{\ast }=\frac{v}{\omega r_2}\end{array}$(19) where $\omega$ is the angular velocity of the impeller, r₂ is the outer radius of the impeller. The velocity vector, streamline, pressure, and EPR distribution of the impeller are shown in Figure 11:

Figure 10. Schematic diagram of special operating conditions points.

Figure 11. Internal flow state at special operating points.

Table 9. Data of special operating condition points.

Operation point	Q(m^3/s)	H(m)	n(r/min)	n11	Q11	M11
A1	0.00289	13.8584	650.17	48.9022	9.902053	−339.252
A9	0.06378	13.0409	650.27	50.41943	225.2759	41.98358
E16	0.15739	12.8116	399.91	31.28374	560.8664	1380.03

In the optimal operating point, no obvious vortical flow can be seen in the impeller stay vane at the optimum operating point from the streamline and speed diagram. It can also be found in the graph of the EPR that the EPR is very small at the optimum state and the overall flow is uniform and stable.

With the decrease of unit flow rate, the unit is running out of control. We can see a certain degree of flow disorder in the runner passage at this time in the streamline diagram. In the speed diagram, we can see that the speed at the junction of stay vane and impeller increases significantly, which may be caused by insufficient flow due to the decrease of unit flow rate. Although the velocity in the stay vane area decreases significantly, the overall flow in the stay vane area is uniform and stable with less vortical flow. In the graph of the EPR, the partial EPR at the outer edge of the impeller is very high, almost in a chaotic state outside the impeller, while the EPR in the stay vane is not obvious. This may be due to the fact that a certain flow still enters the impeller when the flow rate decreases. At this time, the impeller speed is insufficient and the flow is insufficient, which leads to the chaotic flow state at the junction of the guide van. A circular flow with higher velocity formed in the vaneless space makes it difficult for the fluid to flow into the impeller, and leads to vortical flow, we can also see that high e values are mainly concentrated in the vaneless space, which indicates that the flow in the vaneless space is very unstable.

Under the condition of braking operating point, the vortical flow in the stay vane area increases obviously, which leads to an obvious return flow between the vanes of the impeller. The vaneless area of the impeller is full of vortical flow. From the velocity vector diagram, it can be seen that the speed of the connecting part between the stay vane and the blade is very large and most of it is return flow. The fluid speed in the impeller is significantly reduced and the flow is very insufficient. The increase of the entropic energy mainly occurs in the corner part of the outer edge of the impeller and the vaneless area of the stay vane. This may be due to the low flow rate, the absence of a full flow path when the fluid passes through the stay vane, the increase in the EPR caused by the impact of the liquid against the blade, and the continued impact of the fluid from the stay vane to the impeller, resulting in a slight increase in the partial EPR at the outer edge of the blade.

In summary, in each typical unstable working condition, the flow state of the bladeless region of the impeller is the most chaotic and its energy loss is the largest, followed by the outer edge of the blade.

In addition, for the above three special operating points, we conducted a comparative analysis of their blade loads, and selected the inter flow surface of the blade and the parts close to the upper crown and lower ring respectively. It was found that under different operating conditions, the load conditions of the runner blade would also change. The results are shown in Figure 12:

Figure 12. Blade load at E16, A9, A1.

We can see that the load variation range of blade inlet and outlet under optimum conditions is much smaller than that under unstable conditions. It can be seen that as the flow rate decreases, the overall pressure distribution at the inlet of the blade changes slightly, while the pressure at the working surface and back of the outlet side increases to a certain extent. Under unstable operating conditions, the pressure at the middle of the blade is slightly higher, while the maximum pressure of the blade at the outlet side under unstable and braking conditions basically does not change, but the minimum pressure at the outlet side under unstable operating conditions is slightly smaller than other operating conditions. In general, under stable operating conditions, there is a constant pressure difference between the pressure surface and the suction surface of the blade, and the pressure difference gradually increases from the inlet side to the outlet side of the impeller, which is beneficial for the impeller to work on the fluid.

Under runaway operation, there is a constant pressure difference from the outlet edge to the middle of the blade, but due to changes in the inflow angle, the position of the water impinging on the blade changes. At the same time, unstable vortices and backflow in the runner flow channel cause changes in the pressure distribution between the pressure surface and the suction surface of the blade when it is close to the outlet edge of the blade, which can lead to a serious decrease in the power performance of the impeller.

Under braking operation conditions, there is only a certain pressure difference at the inlet side of the blade. The pressure difference between the working surface in the middle of the blade and the back side is basically the same. At this time, the pressure difference is basically 0, so the work done by the runner is 0. Also, because the pressure difference on both sides of the blade head and tail, it turns the blade to work on the water flow. In summary, under braking conditions, the impeller basically does not work, and the confusion of pressure difference can lead to flow separation and vortex generation, indicating that the instability of the flow field is the greatest. In the pressure fluctuation analysis section of three special working conditions, we obtained the results shown in Figure 13:

Figure 13. Pressure of different parts under the same working condition.

After analysing the pressure magnitude of each part of the three special working conditions, we can find that the pressure fluctuation amplitude of the stay vane has always been at a high level in A1 and A9. This is because the flow is fully fixed at the optimal operating point and runaway operating point, which plays a guiding role and alleviates the impact of the fluid, resulting in high pressure. At the braking operation point, the unit flow rate is very small, and the fluid in each flow passage component is very small, so the overall pressure is low.

3.4. Comparative analysis between turbine mode and reverse pump mode of S region

In order to better understand the characteristics of S region, we have selected three pairs of operating points (as shown in Table 10) in the turbine mode and the reverse pump mode with opposite flow directions, similar flow, and the same rotational speed for analysis.

Table 10. Turbine mode and reverse pump mode points.

Turbine mode	Q(m^3/s)	H(m)	n(r/min)	n11	Q11	M11
A3	0.01816	13.6007	650.27	49.37091	62.80861	−244.811
A5	0.0353	13.4692	650.21	49.60675	122.6839	−149.919
A7	0.0525	13.0484	650.24	50.40262	185.3807	−38.2791
Pump mode	Q(m^3/s)	H(m)	n(r/min)	n11	Q11	M11
O1	−0.01394	14.4173	650.17	47.94496	−46.8279	−431.624
O2	−0.03754	13.6267	650.17	49.3162	−129.713	−757.174
O3	−0.05323	12.1718	650.2	52.18283	−194.609	−1123.51

Considering the change process from the braking of the turbine mode to the reverse pump mode, we first extracted the blade load conditions at the operating points O1, O2, O3 of the reverse pump mode, as shown in Figure 14:

Figure 14. Blade load diagram under three conditions of reverse pump mode.

In the blade load diagram of the reverse pump mode, we can see that the blade pressure gradually increases as the flow rate increases, with the rotational speed unchanged. At the inlet edge of the blade, the pressure on the working surface is significantly greater than that on the back. After flowing into the impeller, the fluid impinges on the pressure surface of the leading edge of the blade, forming a high pressure on the pressure surface, while the other side of the blade does not receive the impact, resulting in a small change in pressure, so the leading edge of the blade appears the phenomenon of load jumping. From the inlet side to the middle of the blade, the pressure on both sides of the blade is very small. In the part near the outlet edge of the blade, the pressure gradually increases, and the pressure difference also increases significantly. In addition, we conducted a comparative analysis of blade load and impeller internal flow at these three pairs of approximate operating conditions, as shown in Figure 15:

Figure 15. Comparative analysis of blade loads in turbine mode and reverse pump mode.

After comparing three pairs of operating points with the same rotational speed, similar flow rate, and different flow directions between the turbine mode and the reverse pump mode, we can see that the variation amplitude of blade pressure in the reverse pump mode is around 200,000 Pa, significantly greater than that in the turbine mode. In the turbine mode, there is significant pressure drop and pressure difference reversal at the blade outlet, while in the corresponding reverse pump mode, there is no significant decrease in pressure near the blade outlet, and the pressure reversal is also very small; The same point is that regardless of the positive or negative flow rate, there is always a situation where the pressure on both sides of the blade approaches the area from the blade outlet to the middle of the blade, resulting in a very small pressure difference. This may seriously weaken the power generation ability of the blade, lead to poor flow conditions of the runner, generate backflow vortices, and so on. On the basis of blade load analysis, we extracted the internal flow situation, pressure magnitude, and EPR distribution at corresponding operating points for comparative analysis, as shown in Figures 16 and 17:

Figure 16. Comparative analysis of internal flow status in turbine mode and reverse pump mode.

Figure 17. Comparative analysis of pressure magnitude in turbine mode and reverse pump mode.

In the internal flow analysis of the turbine model, we found that the speed chaos mainly occurs in the blade fillet at the outer edge inlet of the impeller. In the EPR diagram, it can also be seen that the position where the entropy increases mainly exists in the impeller blade fillet. In the pressure analysis, it can be seen that the part with the highest pressure in turbine model is the stay vane, and the pressure in the remaining part is very small, which can be almost ignored. This may be caused by the pressure inversion between the pressure surface and the suction surface of the blade analysed previously; Compared to the turbine mode, the flow state of the reverse pump mode is more chaotic. In the reverse pump mode, the velocity chaos mainly occurs in the circumferential portion of the stay vane outer edge outlet. We can also see that the entropy increases significantly near the stay vane outlet fillet. In the pressure coefficient analysis, we found that the region with the highest pressure in the reverse pump mode is the impeller outlet portion, and next is the impeller blade free zone. This is due to the fact that the working fluid of the reverse pump flows into the stay vane, where the stay vane is impacted by the flow, which plays a role in guiding and reducing circulation, resulting in an increase in entropy production. After the fluid flowing into the impeller from the stay vane, the fluid impinges on the rotating blade leads to pressure instability due to changes in the flow direction at the inlet and outlet. In conclusion, the shape of the impeller blade is crucial to the transition process of the unit. In order to ensure the stability of the unit in the process of pumped storage, the shape of the impeller blade can be improved and optimised.

4. Conclusion

In this paper, the S region characteristics of a large vertical centrifugal pump under reversing operation are simulated numerically and the following conclusions are obtained:

1. Through the analysis of S region, we found that the unit at high unit speed, low unit flow condition of the pressure and pulsation is larger; The pressure is maximum under the condition of low unit speed and high unit flow, but the pulsation amplitude is small. In the three special operating conditions, the best operating conditions have the strongest power capacity. Because the flow rate is too small, the blade does not work and there is basically no pressure change at the braking operating point. Due to the large volume and weight of the unit, it is necessary to pay attention to the compression performance, structural strength and operation stability of the unit under braking conditions in the actual operation process to ensure safety.

2. After comparing the approximate operating conditions of the two modes of the unit, we found that in the turbine mode, the entropy inside the impeller slightly increases, and the pressure fluctuations inside the stationary blades are greater, resulting in more stable unit operation. In reverse pump mode, the impeller pressure and the EPR of the stationary blades is high. In practical applications, it is necessary to avoid entering the reverse water pump area as much as possible, reduce energy loss, and improve the lifespan and efficiency of the unit.

3. It is found that the EPR did not increase significantly under the optimal operating point, and the overall flow was stable; During runaway operating point, there is a very significant entropy increase from the stay vane outlet to the impeller inlet, which is caused by high flow and low rotational speed under this operating condition; Under braking operating point, entropy increase generation mainly occurs in the fillet of the stay vane. In the comparison between the hydraulic turbine mode and the reverse pump mode, there was no significant increase in entropy at the three operating points of the hydraulic turbine mode; At a certain speed, the entropy of the reverse pump mode gradually increases with the increase of flow rate, especially in the rounded part of the impeller.

In summary, this article analyses the reasons and mechanisms for the rare reverse S region operation of large vertical centrifugal pumps from the characteristics of load, energy, and pressure of different flow under different operating conditions. However, there are still certain limitations in this article. In future work, it is also necessary to consider comparative analysis of the operation of units at different scales in the S region, as well as the cavitation and erosion situation of units under high instability conditions.

Acknowledgements

The authors would like to acknowledge the financial support of National Natural Science Foundation of China.

Disclosure statement

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

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [grant numbers U22A20238; 52079142].