Investigation of Pumped Storage Hydropower Power-Off Transient Process Using 3 D Numerical Simulation Based on SP-VOF Hybrid Model

The transient characteristic of the power-off process is investigated due to its close relation to hydraulic facilities’ safety in a pumped storage hydropower (PSH). In this paper, power-off transient characteristics of a PSH station in pump mode was studied using a three-dimensional (3D) unsteady numerical method based on a single-phase and volume of fluid (SP-VOF) coupled model. The computational domain covered the entire flow system, including reservoirs, diversion tunnel, surge tank, pump-turbine unit, and tailrace tunnel. The fast changing flow fields and dynamic characteristic parameters, such as unit flow rate, runner rotate speed, pumping lift, and static pressure at measuring points were simulated, and agreed well with experimental results. During the power-off transient process, the PSH station underwent pump mode, braking mode, and turbine mode, with the dynamic characteristics and inner flow configurations changing significantly. Intense pressure fluctuation occurred in the region between the runner and guide vanes, and its frequency and amplitude were closely related to the runner’s rotation speed and pressure gradient, respectively. While the reversed flow rate of the PSH unit reached maximum, some parameters, such as static pressure, torque, and pumping lift would suddenly jump significantly, due to the water hammer effect. The moment these marked jumps occurred was commonly considered as the most dangerous moment during the power-off transient process, due to the blade passages being clogged by vortexes, and chaos pressure distribution on the blade surfaces. The results of this study confirm that 3D SP-VOF hybrid simulation is an effective method to reveal the hydraulic mechanism of the PSH transient process.


Introduction
To maintain the balance between increasing demands for energy and environment protection, more and more renewable energy resources, such as wind and solar power, have been exploited, even though their intermittent electric output will cause surges in power network [1].Besides, there is a larger peak-valley difference of electric load owing to the increasing net capacity [2].Pumped storage hydropower (PSH) stations, as the most efficient energy storage facilities [3], play an important role in the power grid, such as peaking load shifting, frequency and phase modulation, and as an emergency reserve.Thus, there has been strong interest in developing PSH facilities worldwide, particularly in China [4].To provide various grid services, PSH units have to switch between various operation conditions quickly and frequently.Some transient processes with dramatic changes in flow patterns will pose a threat to the stability and security of PSH plants [5].In addition, most pump-turbines have the unsteady S-shape region in synthetic characteristic curves [6], which can result in more risks than normal hydropower stations in transient processes.
The PSH transient process attracts much attention because it is a key factor for the safety of hydraulic facilities.Experimental method is currently considered as the most reliable way for characterizing the transient process [7].However, many transient processes cannot be evaluated experimentally, due to limitations such as safety, funding, or technical conditions [8].By contrast, numerical methods have the potential to overcome these limitations.The one-dimensional method of characteristics (MOC) is the most popular numerical method, and has been used extensively to investigate the transient process of hydraulic facilities.Thanapandi and Prasad [9] analyzed the transient characteristics of a centrifugal pump during startup and shutdown periods using the MOC.Afshar and Rohani [10] simulated the transient flow caused by load variation of hydroelectric power plants.Wan and Huang [11] presented an improved method, and obtained the complete characteristics and hydraulic transient of centrifugal pump.Zhang [12] researched the influence of water conveyance system layout on the hydraulic transients of pump-turbines load successive rejection in pumped storage station.Zhou [13] used an enhanced multi-objective gravitational search algorithm to optimize guide vane closing schemes of PSH unit.Although the MOC approach is useful to predict the transient characteristics, it could not show details of three-dimensional flow fields [14], such as vortexes, velocity vector, and streamlines, therefore, if involving turbo-machines, it usually relies on synthetic characteristic curves obtained using model experiments [15].
Nowadays, three-dimensional numerical simulations based on computational fluid dynamics (CFD) have been an important means of predicting hydraulic facilities performance under normal conditions.Yang and Wu [16] presented an innovative optimization method on runner blades in bulb turbine based on CFD analysis.Litfin [17] analyzed and designed a two-bladed wastewater pump, and Sedlar [18] researched the cavitation phenomena in a mixed-flow pump.Recently, this was also used to simulate the transient processes.Xia [19] simulated the 3D transient flow patterns in a corridor-shaped air-cushion surge chamber of a hydropower station.Nicolle [20] modeled the startup phase of a Francis turbine while it went from rest to speed under no-load condition.Zhou [21] simulated runaway transient of a propeller turbine model.A combined approach, the pipe system part using 1D-MOC method and the hydro unit part using 3D method, was used to research hydraulic turbine flows in transient operating regimes [22][23][24][25].Liu [26] simulated a part of the transient power failure process of a prototype pump-turbine.Li [27] used dynamic mesh to simulate the 3D dynamic turbulent flow of whole flow channels in start transition process and runaway process of bulb hydraulic turbine.Mao [28] investigated unsteady flow field in a reversible pump-turbine during the continuous load rejection by a 3D CFD analysis.Luo [29] carried out a 3D transient numerical simulation method to simulate the runaway process of tubular turbine through secondary development of CFX software and Fortran.However, due to the limitations of the computer capacity and the grid technology to deal with large solid boundary deformation, the boundary conditions are usually simplified so that the transient flow mechanism in the 3D and two-phase system has not been fully studied.This paper constructs a 3D transient hybrid simulation method that combines normal single-phase model with VOF model (SP-VOF), and verifies the model using a simple physical model calculation.Then, the SP-VOF hybrid model is applied to simulate the entire power-off transient process of the whole PSH flow system in pump mode, and reveal the varying law of dynamic characteristics parameters.

Governing Equations
The flow media in pump-turbine units and pipes is mainly regarded as liquid water, due to ignoring cavitation, and its continuity equation can be written as follows: and the conservation of momentum is described by where ρ is the density of fluid, t is time, u is velocity, p is the static pressure, = τ is the stress tensor, g is the gravitational body force, and F is external body forces.For water in single-phase model, F = 0.The stress tensor = τ is given by where µ is the molecular viscosity, I is the unit tensor, and the second term on the right hand side is the effect of volume dilation.
It is well known that the water will fluctuate in a surge tank during the transient process as a two-phase flow pattern.The volume of fluid (VOF) two-phase model was often used to research flow problems about pipes and surge tanks [30,31].For the VOF model [32], the form of governing equations is roughly the same as that of single-phase model, but there are two significant differences between them.One is that the stress tensor = τ in the VOF model is given by Equation (4).The other is that a single momentum equation is solved throughout the domain, and the resulting velocity field is shared among the phases.The momentum equation is dependent on the volume fractions of all phases through the properties ρ and µ.For example, if the water and air are represented by the subscripts 1 and 2 respectively, the density ρ will be given by Equation (5).All other properties are computed in this manner.
The tracking of the level between water and air is accomplished by the continuity equation for the volume fraction of each phase.For the phase of water, this equation takes the following form where .m 21 is the mass transfer from water to air, .m 12 is that from phase air to phase water, S α 1 is the source term, by default, S α 1 = 0. α 1 is the volume fraction of water in a cell, and the volume fraction of air, marked by α 2 , is calculated simply through α 2 = 1 − α 1 .

Turbulence Model
In consideration of many components and complex flow fields in pumped storage power stations, we require the turbulence model with good adaptability to close governing equations mentioned above.The realizable k − ε model [33], which is an improvement on the standard k − ε model, has been extensively validated for a wide range of flows, including the rotating homogeneous shear flows, the free flows with jets or mixed layers, the pipe flow and the separation flows.

3D Hybrid Simulation Method
After several attempts, the VOF model seemed a more feasible method to simulate the 3D and multiphase surge tank.However, normal single-phase model is the better choice for simulation of pipes and units flow field with the faster computation speed and more accurate results.With comprehensive considerations, we developed an innovative simulation method combining single phase and the VOF (SP-VOF) model.Specifically, the VOF model was only applied to surge tank and single-phase model was used in other regions, such as pipes and the unit.The two numerical models are solved by two independent CFD codes respectively and their properties data are transferred through the shared interface.If the static pressure p of the single-phase model and VOF model are respectively represented by the superscripts 1 and 2, and iteration step is represented by subscripts n.The transfer equation of static pressure from VOF model to single-phase model is written as follows: where β is the relaxation factor, and β = 0.9.All other properties, such as velocity and parameters of turbulence energy, are computed in this manner.
To verify the reliability and accuracy of the hybrid simulation method, a simplified model consisting of a surge tank and a portion of pipe was set up (Figure 1a), and one transient process with the same boundary conditions was simulated, which were defined as follows: (1) The tank outlet was opened to atmosphere, and its static pressure remained at zero.(2) The water flow rate at the pipe inlet was equal to 100 m 3 /s, and that at the pipe outlet varied with time (Figure 1b).The initial water volume in surge tank was initialized as Figure 1a.Finally, the flow rate across the interface of the hybrid method was compared with that of the single VOF method.
Flow rate values of the hybrid method agreed well with the VOF simulation results (Figure 1b), and the maximum relative error of flow rate was 1.32% when the time was equal to 5 s.In addition, the pumped storage power station was much bigger in volume than the surge tank, so the errors at other parts caused by the hybrid simulation method would be very small.

3D hybrid Simulation Method
After several attempts, the VOF model seemed a more feasible method to simulate the 3D and multiphase surge tank.However, normal single-phase model is the better choice for simulation of pipes and units flow field with the faster computation speed and more accurate results.With comprehensive considerations, we developed an innovative simulation method combining single phase and the VOF (SP-VOF) model.Specifically, the VOF model was only applied to surge tank and single-phase model was used in other regions, such as pipes and the unit.The two numerical models are solved by two independent CFD codes respectively and their properties data are transferred through the shared interface.If the static pressure p of the single-phase model and VOF model are respectively represented by the superscripts 1 and 2, and iteration step is represented by subscripts n.The transfer equation of static pressure from VOF model to single-phase model is written as follows:   where  is the relaxation factor, and  = 0.9.All other properties, such as velocity and parameters of turbulence energy, are computed in this manner.
To verify the reliability and accuracy of the hybrid simulation method, a simplified model consisting of a surge tank and a portion of pipe was set up (Figure 1a), and one transient process with the same boundary conditions was simulated, which were defined as follows: (1) The tank outlet was opened to atmosphere, and its static pressure remained at zero.(2) The water flow rate at the pipe inlet was equal to 100 m 3 /s, and that at the pipe outlet varied with time (Figure 1b).The initial water volume in surge tank was initialized as Figure 1a.Finally, the flow rate across the interface of the hybrid method was compared with that of the single VOF method.
Flow rate values of the hybrid method agreed well with the VOF simulation results (Figure 1b), and the maximum relative error of flow rate was 1.32% when the time was equal to 5 s.In addition, the pumped storage power station was much bigger in volume than the surge tank, so the errors at other parts caused by the hybrid simulation method would be very small.

Algorithm of Power-Off Transient Simulation
Algorithm routine of power-off transient simulation was implemented by a secondary developed and customized program basing on widely used CFD software, Fluent 6.3 (ANSYS, Canonsburg, PA, USA) and the algorithm routine diagram was described as follows (Figure 2):

Geometry Model and Initial Parameters
The pumped storage hydropower station involved in this simulation consists of an upstream reservoir, a diversion tunnel, surge tank, pump-turbine units, and so on, whose detailed structure was displayed in

Grid Technology
According to physical characteristics of different regions, we made different meshing schemes.Narrow prism grids were meshed in pipes, and tetrahedral grids were meshed at other places.We also performed the grid refinement for some places, such as the region near pipe wall, the runner and the guide vanes.Varying speed sliding mesh technology was used to transfer flow field parameters between runner zone and adjacent part passages, meanwhile, mesh rotating of guide vanes was resolved by dynamic mesh technology [34].

Equations Discretization and Boundary Condition
The governing equations of the hybrid SP-VOF method were discrete with finite volume

Geometry Model and Initial Parameters
The pumped storage hydropower station involved in this simulation consists of an upstream reservoir, a diversion tunnel, surge tank, pump-turbine units, and so on, whose detailed structure was displayed in Figure 3

Geometry Model and Initial Parameters
The pumped storage hydropower station involved in this simulation consists of an upstream reservoir, a diversion tunnel, surge tank, pump-turbine units, and so on, whose detailed structure was displayed in Figure 3

Grid Technology
According to physical characteristics of different regions, we made different meshing schemes.Narrow prism grids were meshed in pipes, and tetrahedral grids were meshed at other places.We also performed the grid refinement for some places, such as the region near pipe wall, the runner and the guide vanes.Varying speed sliding mesh technology was used to transfer flow field parameters between runner zone and adjacent part passages, meanwhile, mesh rotating of guide vanes was resolved by dynamic mesh technology [34].

Equations Discretization and Boundary Condition
The governing equations of the hybrid SP-VOF method were discrete with finite volume

Grid Technology
According to physical characteristics of different regions, we made different meshing schemes.Narrow prism grids were meshed in pipes, and tetrahedral grids were meshed at other places.We also performed the grid refinement for some places, such as the region near pipe wall, the runner and the guide vanes.Varying speed sliding mesh technology was used to transfer flow field parameters Energies 2018, 11, 1020 6 of 16 between runner zone and adjacent part passages, meanwhile, mesh rotating of guide vanes was resolved by dynamic mesh technology [34].

Equations Discretization and Boundary Condition
The governing equations of the hybrid SP-VOF method were discrete with finite volume method, and the first order implicit format was used for time item.For the single-phase model applied in pipes and unit, the second order format was used for pressure item, the second order upwind format for convection item, and Semi-implicit method for pressure linked equations-consistent (SIMPLEC) algorithm was chosen to obtain the coupling solution from velocity and pressure equations.For the VOF model in the surge tank, the PRESTO!Format was used for pressure item, the first order upwind format for momentum item and turbulent kinetic energy, and PISO algorithm for velocity-pressure coupling solution.
Due to small change of reservoir water level in the transient process, their inlet or outlet pressure value was assumed to be constant, according to their water depth, and the surge tank outlet was opened to atmosphere, which had a relative pressure that was equal to zero.We also set the roughness values of different passage walls: 2.55 mm for surge tank, diversion system, and tailrace system; 0.06 mm in the casing; 0.0016 mm on runner blades and vanes.In addition, runner rotate speed was calculated by Equation ( 8), and closing law of guide vanes was given by Equation (9).
where ω is the runner rotate speed, M is resultant torque of the rotational system, ∆t is time step size, and ∆γ is the change amount of guide vane opening.

Data Monitoring and Processing
In this study, for the assessment of this transient process, and some dynamic parameters, such as runner rotate speed n, pumping lift H, unit flow rate Q and torque M, are obtained and normalized and divided by corresponding parameters values under the initial condition.However, the surge tank flow rate Q st is normalized by initial unit flow, because its initial value is regarded as zero.In order to monitor variable regularities of static pressure, three measuring points were set in the pump-turbine passage.Point 1 was at the inlet section center of spiral case, point 2 was at the half height of trailing edge of guide vane, and point 3 was at the center of draft tube inlet.

Test Results of Grid and Time Dependency
To eliminate influences caused by computational grids, three different meshing schemes are employed in the CFD simulation: the corresponding mesh grid size is about 2.5 million, 3.5 million, and 4.5 million respectively.The two parameters of pre-calculated results, namely unit rotate speed and static pressure of spiral casing inlet, were plotted in Figure 4. Though the two parameters' overall varying tendency under the three meshing schemes is similar, the local details still retain slightly different, especially under the intense fluctuation region.Obviously, the curves of scheme 2 are closer to that of scheme 3 than that of scheme 1.Thus, scheme 2 was chosen based on the available computing power and subdivision scheme for each part of flow passage, as was listed in the Table 1.Finally, Figure 5 shows some details of the mesh in different regions of the PSH model.Similarly, different time step sizes also have influences on the unsteady numerical results of the transient parameters.Three different time step sizes were used to pre-compute: 0.02 s, 0.01 s, and 0.005 s, respectively.The corresponding results of runner rotating speed were described in Figure 6, and showed that the shorter time step led to the slower decreasing rotating speed and smaller runaway speed value.To obtain more detailed transient information, the 0.005 s time step size was adopted finally.Similarly, different time step sizes also have influences on the unsteady numerical results of the transient parameters.Three different time step sizes were used to pre-compute: 0.02 s, 0.01 s, and 0.005 s, respectively.The corresponding results of runner rotating speed were described in Figure 6, and showed that the shorter time step led to the slower decreasing rotating speed and smaller runaway speed value.To obtain more detailed transient information, the 0.005 s time step size was adopted finally.Similarly, different time step sizes also have influences on the unsteady numerical results of the transient parameters.Three different time step sizes were used to pre-compute: 0.02 s, 0.01 s, and 0.005 s, respectively.The corresponding results of runner rotating speed were described in Figure 6, and showed that the shorter time step led to the slower decreasing rotating speed and smaller runaway speed value.To obtain more detailed transient information, the 0.005 s time step size was adopted finally.

The Varying Law of Dynamic Characteristics Parameters
Once the transient process began, the PSH station would undergo pump mode, braking mode, and turbine mode in sequence (Figure 7).In the beginning period, unit flow rate Q, pumping lift H and torque M declined quickly with decreasing runner rotate speed n after power failure.Nevertheless, the surge tank flow rate Qst rose rapidly, and supplied water to the diversion tunnel due to its pressure decreasing.As time arrived at 6.4 s, unit flow rate Q was equal to 0, which indicated the end of pump mode.Meanwhile, pumping lift H and torque M reached the minimum values, and from then on, the unit was running in braking mode for 11.5 s.During this time, reversed rate Q, torque M, and pumping lift H had been increasing until t = 11.6 s, and then decreased slowly with guide vanes closed.It was noted that there were sharp turning points on the curves of Q, Qst, M , and H at t = 11.6 s because of the water hammer effect, and more detailed discussion will be conducted later, combined with flow field characteristics.After t = 17.8 s, the unit entered into turbine mode marked as the unit rotating speed direction changed, and all dynamic parameters became stable gradually, except the small changes at t = 22.95 s for the variable closing law of guide vanes.At last, runner rotate speed got the maximum reversed value 43.8 r/min when the torque M was equal to approximately zero, and the lift H was lower than the initial one, of which the difference was approximate to the hydraulic loss of pipes in initial pump mode.
In this simulation of power-off transient process, surge tank always supplied water to diversion tunnel.Before t = 10.1 s, the water flow rate from the surge tank had been increasing until reaching the maximum value (112% of the unit flow rate).Then, the discharge of the surge tank decreased gradually due to the decline of the water level in the surge tank, and increased diversion tunnel pressure.

The Varying Law of Dynamic Characteristics Parameters
Once the transient process began, the PSH station would undergo pump mode, braking mode, and turbine mode in sequence (Figure 7).In the beginning period, unit flow rate Q, pumping lift H and torque M declined quickly with decreasing runner rotate speed n after power failure.Nevertheless, the surge tank flow rate Q st rose rapidly, and supplied water to the diversion tunnel due to its pressure decreasing.As time arrived at 6.4 s, unit flow rate Q was equal to 0, which indicated the end of pump mode.Meanwhile, pumping lift H and torque M reached the minimum values, and from then on, the unit was running in braking mode for 11.5 s.During this time, reversed flow rate Q, torque M, and pumping lift H had been increasing until t = 11.6 s, and then decreased slowly with guide vanes closed.It was noted that there were sharp turning points on the curves of Q, Q st , M, and H at t = 11.6 s because of the water hammer effect, and more detailed discussion will be conducted later, combined with flow field characteristics.After t = 17.8 s, the unit entered into turbine mode marked as the unit rotating speed direction changed, and all dynamic parameters became stable gradually, except the small changes at t = 22.95 s for the variable closing law of guide vanes.At last, runner rotate speed got the maximum reversed value 43.8 r/min when the torque M was equal to approximately zero, and the lift H was lower than the initial one, of which the difference was approximate to the hydraulic loss of pipes in initial pump mode.
In this simulation of power-off transient process, surge tank always supplied water to diversion tunnel.Before t = 10.1 s, the water flow rate from the surge tank had been increasing until reaching the maximum value (112% of the unit flow rate).Then, the discharge of the surge tank decreased gradually due to the decline of the water level in the surge tank, and increased diversion tunnel pressure.
The dynamic loads on the runner affect the stability of the pump-turbine unit directly.The axial force F z and radial force F r acting on the runner were plotted in Figure 7.The absolute value of axial force F z had been decreasing in pump mode and fluctuating in braking mode until the force direction changed at t = 11.6 s, when the reversed unit flow rate reached maximum value.Then, the axial force F z jumped to the maximum reverse value of nearly 10 6 N suddenly, concurrently with the pumping lift and torque.By comparison, the variation of radial force F r was smaller, but fluctuated at a higher frequency in the end stage of pump mode when the unit flow rate went down to the minimum.The dynamic loads on the runner affect the stability of the pump-turbine unit directly.The axial force Fz and radial force Fr acting on the runner were plotted in Figure 7.The absolute value of axial force Fz had been decreasing in pump mode and fluctuating in braking mode until the force direction changed at t = 11.6 s, when the reversed unit flow rate reached maximum value.Then, the axial force Fz jumped to the maximum reverse value of nearly 10 6 N suddenly, concurrently with the pumping lift and torque.By comparison, the variation of radial force Fr was smaller, but fluctuated at a higher frequency in the end stage of pump mode when the unit flow rate went down to the minimum.

Comparisons between Numerical and Experimental Results
To verify the reliability and the accuracy of the power-off transient simulation, some numerical values, namely runner rotating speed and static of measuring points (Figure 8), were compared with prototype test results (Figure 9).The overall trends of numerical results agreed well with experimental ones, particularly under the pump and turbine conditions located on both sides of the figure.However, the rotating speed derived from numerical approach decreased a bit faster than experimental results, which may be due to that the moment of inertia value used in the numerical simulation was not exactly the same as the actual value.In addition, the pressure fluctuation characteristics between the two methods had significant discrepancy, especially when pump-turbine unit was operated in the braking mode.Compared with experimental data, the numerical results displayed the higher pressure amplitude at

Comparisons between Numerical and Experimental Results
To verify the reliability and the accuracy of the power-off transient simulation, some numerical values, namely runner rotating speed and static pressure of measuring points (Figure 8), were compared with prototype test results (Figure 9).The dynamic loads on the runner affect the stability of the pump-turbine unit directly.The axial force Fz and radial force Fr acting on the runner were plotted in Figure 7.The absolute value of axial force Fz had been decreasing in pump mode and fluctuating in braking mode until the force direction changed at t = 11.6 s, when the reversed unit flow rate reached maximum value.Then, the axial force Fz jumped to the maximum reverse value of nearly 10 6 N suddenly, concurrently with the pumping lift and torque.By comparison, the variation of radial force Fr was smaller, but fluctuated at a higher frequency in the end stage of pump mode when the unit flow rate went down to the minimum.

Comparisons between Numerical and Experimental Results
To verify the reliability and the accuracy of the power-off transient simulation, some numerical values, namely runner rotating speed and static pressure of measuring points (Figure 8), were compared with prototype test results (Figure 9).The overall trends of numerical results agreed well with experimental ones, particularly under the pump and turbine conditions located on both sides of the figure.However, the rotating speed derived from numerical approach decreased a bit faster than experimental results, which may be due to that the moment of inertia value used in the numerical simulation was not exactly the same as the actual value.In addition, the pressure fluctuation characteristics between the two methods had significant discrepancy, especially when pump-turbine unit was operated in the braking mode.Compared with experimental data, the numerical results displayed the higher pressure amplitude at The overall trends of numerical results agreed well with experimental ones, particularly under the pump and turbine conditions located on both sides of the figure.However, the rotating speed derived from numerical approach decreased a bit faster than experimental results, which may be due to that the moment of inertia value used in the numerical simulation was not exactly the same as the actual value.In addition, the pressure fluctuation characteristics between the two methods had significant discrepancy, especially when pump-turbine unit was operated in the braking mode.Compared with experimental data, the numerical results displayed the higher pressure amplitude at point 2, but lower amplitude at point 1 and point 3, besides the remarkable pressure jump at t = 11.6 s.The three possible explanations for their differences could be deduced as follows: (1) they have different acquisition frequency; (2) some kinds of vortexes or cavitation could not be fully simulated due to the limitations of the using numerical model and grid size, so that their contributions to the pressure fluctuation could not be reflected; (3) the measured pressure data from field test may include the vibration signals of the unit.point 2, but lower amplitude at point 1 and point 3, besides the remarkable pressure jump at t = 11.6 s.The three possible explanations for their differences could be deduced as follows: (1) they have different acquisition frequency; (2) some kinds of vortexes or cavitation could not be fully simulated due to the limitations of the using numerical model and grid size, so that their contributions to the pressure fluctuation could not be reflected; (3) the measured pressure data from field test may include the vibration signals of the unit.Then, the more details of the static pressure variation obtained by the numerical method would be described and explained.The static pressure of point 1 at the spiral casing had been declining from 2.6 MPa to 1.9 MPa near the end of pump mode, due to pumping energy reduction.Once entering the braking mode, it showed an uptrend and ascended to 2.63 MPa suddenly, when the maximum reversed flow rate appeared.Since the guide vanes' closing speed changed, it stepped to a lower level at t = 23 s.In pump mode, the variation law of static pressure at point 2 was same as that of point 1, but presented strong pulsation in braking mode, of which the frequency was related with the runner rotate speed, and the maximum amplitude, equal to about 1 MPa, was recorded when the maximum reversed flow rate appeared, too.In the whole process, the variation range of static pressure at point 2 was also much bigger than that of other two points, decreasing from 2 MPa to 0.35 MPa.The curve of point 3 showed that static pressure in draft tube had an opposite variation trend towards that in the spiral casing, which caused the decrease of pumping lift H.In summary, the static pressure during power-off transient process had a complex variation, not only related with position, but also affected by dynamic parameters, such as runner rotate speed, unit flow rate, and guide vanes opening.So the visualization of transient flow field can help to increase our understanding of hydrodynamic mechanism.

Transient Flow Fields
Obviously, the flow field in pumped storage power station passage will take on different features during power-off transient process, which is the intrinsic factor leading to variations of dynamic parameters.In the following, transient flow fields of different parts at some critical moments are plotted and visualized for analysis.
Figure 10 displayed flow fields between spiral casing and guide vanes at different moments.The initial condition had a perfect performance with water flowing upstream smoothly (Figure 10a).When pump rotational speed dropped, unit flow rate and static pressure decreased during the pump mode.When the unit just changed its flow direction at t = 6.5 s, the water in the region flowed slowly and disorderly, with minimum static pressure (Figure 10b).Then, in braking mode at t = 12 s, the water flowed downstream smoothly in the region, and static pressure in the spiral casing Then, the more details of the static pressure variation obtained by the numerical method would be described and explained.The static pressure of point 1 at the spiral casing had been declining from 2.6 MPa to 1.9 MPa near the end of pump mode, due to pumping energy reduction.Once entering the braking mode, it showed an uptrend and ascended to 2.63 MPa suddenly, when the maximum reversed flow rate appeared.Since the guide vanes' closing speed changed, it stepped to a lower level at t = 23 s.In pump mode, the variation law of static pressure at point 2 was same as that of point 1, but presented strong pulsation in braking mode, of which the frequency was related with the runner rotate speed, and the maximum amplitude, equal to about 1 MPa, was recorded when the maximum reversed flow rate appeared, too.In the whole process, the variation range of static pressure at point 2 was also much bigger than that of other two points, decreasing from 2 MPa to 0.35 MPa.The curve of point 3 showed that static pressure in draft tube had an opposite variation trend towards that in the spiral casing, which caused the decrease of pumping lift H.In summary, the static pressure during power-off transient process had a complex variation, not only related with position, but also affected by dynamic parameters, such as runner rotate speed, unit flow rate, and guide vanes opening.So the visualization of transient flow field can help to increase our understanding of hydrodynamic mechanism.

Transient Flow Fields
Obviously, the flow field in pumped storage power station passage will take on different features during power-off transient process, which is the intrinsic factor leading to variations of dynamic parameters.In the following, transient flow fields of different parts at some critical moments are plotted and visualized for analysis.
Figure 10 displayed flow fields between spiral casing and guide vanes at different moments.The initial condition had a perfect performance with water flowing upstream smoothly (Figure 10a).When pump rotational speed dropped, unit flow rate and static pressure decreased during the pump mode.When the unit just changed its flow direction at t = 6.5 s, the water in the region flowed slowly and disorderly, with minimum static pressure (Figure 10b).Then, in braking mode at t = 12 s, the water flowed downstream smoothly in the region, and static pressure in the spiral casing achieved a peak value when the reverse discharge reached maximum value (Figure 10c).After the guide vanes were closed completely, the flow pattern in the channels between stay vanes became disorder and the circumfluent flow with small velocity between the stay vanes and guide vanes emerged, due to the previous flow inertia (Figure 10d).achieved a peak value when the reverse discharge reached maximum value (Figure 10c).After the guide vanes were closed completely, the flow pattern in the channels between stay vanes became disorder and the circumfluent flow with small velocity between the stay vanes and guide vanes emerged, due to the previous flow inertia (Figure 10d).The runner is the most important part of a pump-turbine unit [35], which not only decides steady state performance, but also has great influence on transient characteristics.Obviously, there was no vortex or flow separation in the runner region under the beginning condition (Figure 11a).Once power-off happened, flow configuration in the runner would undergo very complex and unsteady change.When the runner rotate speed reduced 40% at t = 6.5 s, low speed vortex flow appeared at the inlet of runner due to the loss of pumping function for water (Figure 11b).Since then, water changed flow direction, but direction of runner rotate speed remained the same as initial one.As a result, there were jets towards the front of blade, leading to high static pressure on front blade and negative pressure on the middle part.Vortex flow also became stronger, not only behind the leading edge of the blade, but also near the trail of guide vanes (Figure 11c), which could be the typical flow characteristics in braking mode.In addition, the amplitude of pressure fluctuation at point 2 was determined by varying rotor-stator interactions between runner blades and guide vanes.With guide vanes closed, as well as the variations of both reversed flow rate and rotating speed, the angle between jets flowing from guide vanes and blade chord line, the so-called angle of attack, changed from acute angle to right angle, and continued to increase from right angle to obtuse angle, observed from Figure 11c-e.During the first stage, the vector vertical component of jets towards blade strengthened gradually, so that the drag torque exhibited corresponding increase.However, when the angle of attack reached nearly 90°, the maximum value of pressure on the blade convex face occurred, with a large vortex formed among the inlet of blade passages (Figure 11d), which caused water flow to be clogged suddenly, and led to water hammer resulting in sharp turning points at t = 11.8 s on curves plotted in Figure 7.While entering into the turbine mode, runner was rotating reversely, pushed by water.At t = 33 s, although the complete closure of guide vanes broke the connection between runner and spiral casing, the runner still maintained runaway condition for a short time with no difference of static pressure on both side of runner blades (Figure 11f).Then, the unit rotating speed would continue to decline because of the effects of drag torque, until it was forced to brake.The runner is the most important part of a pump-turbine unit [35], which not only decides steady state performance, but also has great influence on transient characteristics.Obviously, there was no vortex or flow separation in the runner region under the beginning condition (Figure 11a).Once power-off happened, flow configuration in the runner would undergo very complex and unsteady change.When the runner rotate speed reduced 40% at t = 6.5 s, low speed vortex flow appeared at the inlet of runner due to the loss of pumping function for water (Figure 11b).Since then, water changed flow direction, but direction of runner rotate speed remained the same as initial one.As a result, there were jets towards the front of blade, leading to high static pressure on front blade and negative pressure on the middle part.Vortex flow also became stronger, not only behind the leading edge of the blade, but also near the trail of guide vanes (Figure 11c), which could be the typical flow characteristics in braking mode.In addition, the amplitude of pressure fluctuation at point 2 was determined by varying rotor-stator interactions between runner blades and guide vanes.With guide vanes closed, as well as the variations of both reversed flow rate and rotating speed, the angle between jets flowing from guide vanes and blade chord line, the so-called angle of attack, changed from acute angle to right angle, and continued to increase from right angle to obtuse angle, observed from Figure 11c-e.During the first stage, the vector vertical component of jets towards blade strengthened gradually, so that the drag torque exhibited corresponding increase.However, when the angle of attack reached nearly 90 • , the maximum value of pressure on the blade convex face occurred, with a large vortex formed among the inlet of blade passages (Figure 11d), which caused water flow to be clogged suddenly, and led to water hammer resulting in sharp turning points at t = 11.8 s on curves plotted in Figure 7.While entering into the turbine mode, runner was rotating reversely, pushed by water.At t = 33 s, although the complete closure of guide vanes broke the connection between runner and spiral casing, the runner still maintained runaway condition for a short time with no difference of static pressure on both side of runner blades (Figure 11f).Then, the unit rotating speed would continue to decline because of the effects of drag torque, until it was forced to brake.
Though geometric structure of draft tube is simple, its transient flow pattern changed significantly since it connected with the runner directly.According to Figure 12a, initial flow fields in the draft tube were harmonious before the transient process happened.With the runner lifting capacity decreasing, static pressure of the draft tube inlet had been increasing in pump mode, due to the effect of water hammer formed by the flow rate difference between the inlet and outlet of the runner.After entering into braking mode, the water would flow into the draft tube via the outside lane of the inlet, and flow back to the runner through the inlet center.When reversed flow rate reached maximum near t = 12 s, the flow regime became very disordered (Figure 12d), and the water flowed into draft tube choppily with 9 branched vortexes induced by the runner, and the static pressure in the entrance center of draft tube is lowest and asymmetrical.When the unit entered turbine mode, the flow became much gentler than that in braking mode.After t = 30 s, unit flux was equal to zero again, but the flow in draft tube was driven by the runner and whirled, and static pressure took on layered distribution, due to gravity effect (Figure 12f).Though geometric structure of draft tube is simple, its transient flow pattern changed significantly since it connected with the runner directly.According to Figure 12a, initial flow fields in the draft tube were harmonious before the transient process happened.With the runner lifting capacity decreasing, static pressure of the draft tube inlet had been increasing in pump mode, due to the effect of water hammer formed by the flow rate difference between the inlet and outlet of the runner.After entering into braking mode, the water would flow into the draft tube via the outside lane of the inlet, and flow back to the runner through the inlet center.When reversed flow rate reached maximum near t = 12 s, the flow regime became very disordered (Figure 12d), and the water flowed into draft tube choppily with 9 branched vortexes induced by the runner, and the static pressure in the entrance center of draft tube is lowest and asymmetrical.When the unit entered turbine mode, the flow became much gentler than that in braking mode.After t = 30 s, unit flux was equal to zero again, but the flow in draft tube was driven by the runner and whirled, and static pressure took on layered distribution, due to gravity effect (Figure 12f).

Conclusions
The 3D SP-VOF hybrid model is an effective method to study the transient process of the whole PSH station flow system in pump mode.During the power-off process, the hydraulic performance of the pump-turbine changed significantly as well as its inner flow configurations while undergoing pump mode, braking mode, and turbine mode, sequentially.Especially in the braking mode, the flow is the most unstable with the great pressure pulsation.The flow was blocked by the strong vortexes in the runner blade channels at t = 11.6 s under the 40% of the rated rotate speed and 34% of

Conclusions
The 3D SP-VOF hybrid model is an effective method to study the transient process of the whole PSH station flow system in pump mode.During the power-off process, the hydraulic performance of the pump-turbine changed significantly as well as its inner flow configurations while undergoing pump mode, braking mode, and turbine mode, sequentially.Especially in the braking mode, the flow is the most unstable with the great pressure pulsation.The flow was blocked by the strong vortexes in the runner blade channels at t = 11.6 s under the 40% of the rated rotate speed and 34% of the initial guide blade opening, which induced a notable water hammer and resulted in the sudden changes of the dynamic parameters.The axial force on the runner is not larger than the initial value during the transient process, while the radial force will fluctuate irregularly in braking mode.The pressure fluctuation between guide vanes and runner is the most intense, whose frequency is determined by runner rotate speed and amplitude, and is related to the dynamic rotor-stator interactions between the guide vanes and runner blades.
Future studies should establish a high-precision and compressible 3D numerical model that considers the effects of the water density change on the results during the transient process [36,37], and the calculation conditions should be increased so that the transient process of power-off transient process in PSH could be more precisely and fully investigated.

Figure 1 .
Figure 1.(a) Simplified pipe and surge tank model; (b) The flow rate across the shared interface obtained by two different methods and the flow rate of pipe outlet varied with the time.

Figure 1 . 16 Figure 2 .
Figure 1.(a) Simplified pipe and surge tank model; (b) The flow rate across the shared interface obtained by two different methods and the flow rate of pipe outlet varied with the time.

Figure 3 .
The length of diversion tunnel, penstock, and tailrace tunnel is 1260 m, 190 m, and 355 m, respectively, and the corresponding diameter is 9 m, 5.6 m, and 10 m.The specific speed of pump-turbine unit ns = 182 m•kW, and its runner outlet diameter D = 3.57 m, total unit moment of inertia J = 4,825,000 kg•m 2 .In addition, the pump-turbine unit has 20 stay vanes, 20 guide vanes and 9 runner blades.Some important parameters of initial condition are as follows: pumping lift 0 205m  H , unit flow rate Q0 = 139 m 3 /s, runner rotating speed n0 = 250 r/min, the initial torque M0 = −11.7 × 10 6 N, guide vane opening angle γ = 23.75°(75% of fully opening), and input power P0 = 306 kW.

Figure 3 .
Figure 3.The whole flow system of the pumped storage power station.
. The length of diversion tunnel, penstock, and tailrace tunnel is 1260 m, 190 m, and 355 m, respectively, and the corresponding diameter is 9 m, 5.6 m, and 10 m.The specific speed of pump-turbine unit ns = 182 m•kW, and its runner outlet diameter D = 3.57 m, total unit moment of inertia J = 4,825,000 kg•m 2 .In addition, the pump-turbine unit has 20 stay vanes, 20 guide vanes and 9 runner blades.Some important parameters of initial condition are as follows: pumping lift 0 205m  H , unit flow rate Q0 = 139 m 3 /s, runner rotating speed n0 = 250 r/min, the initial torque M0 = −11.7 × 10 6 N, guide vane opening angle γ = 23.75°(75% of fully opening), and input power P0 = 306 kW.

Figure 3 .
Figure 3.The whole flow system of the pumped storage power station.

Figure 3 .
Figure 3.The whole flow system of the pumped storage power station.

Figure 4 .
Figure 4. Pre-calculated results of different meshing schemes.

Figure 5 .
Figure 5. (a) Local section of pipeline grid; (b) Grid near the guide vane in the initial time; (c) Grid near the guide vane at full closed time.

Figure 4 .
Figure 4. Pre-calculated results of different meshing schemes.

Figure 4 .
Figure 4. Pre-calculated results of different meshing schemes.

Figure 5 .
Figure 5. (a) Local section of pipeline grid; (b) Grid near the guide vane in the initial time; (c) Grid near the guide vane at full closed time.

Figure 5 .
Figure 5. (a) Local section of pipeline grid; (b) Grid near the guide vane in the initial time; (c) Grid near the guide vane at full closed time.

Figure 6 .
Figure 6.Pre-calculated results of different time step schemes.

Figure 6 .
Figure 6.Pre-calculated results of different time step schemes.

Figure 7 .
Figure 7.The left Y-axis is shared by normalized value of dynamic parameters variations with time; the right Y-axis is shared by axial and radial force on runner vs time.

Figure 8 .
Figure 8. Measuring points of static pressure in the passage.

Figure 7 .
Figure 7.The left Y-axis is shared by normalized value of dynamic parameters variations with time; the right Y-axis is shared by axial and radial force on runner vs. time.

Figure 7 .
Figure 7.The left Y-axis is shared by normalized value of dynamic parameters variations with time; the right Y-axis is shared by axial and radial force on runner vs time.

Figure 8 .
Figure 8. Measuring points of static pressure in the passage.

Figure 8 .
Figure 8. Measuring points of static pressure in the passage.

Figure 9 .
Figure 9. Curves of experiment and simulation results.

Figure 9 .
Figure 9. Curves of experiment and simulation results.

Table 1 .
Mesh subdivision schemes of the PSH station.

Table 1 .
Mesh subdivision schemes of the PSH station.

Table 1 .
Mesh subdivision schemes of the PSH station.