Research on Hydraulic Characteristics in Diversion Pipelines under a Load Rejection Process of a PSH Station

Transient analysis in diversion pipelines should be performed to ensure the safety of a hydropower system. After the establishment of a three-dimensional (3D) geometric model from the part upstream reservoir to the diversion pipeline end in a pumped storage hydropower (PSH) station, the hydraulic characteristics of the diversion system were solved by Reynold average Navier–Stokes (RANS) equations based on a volume of fluid (VOF) method under the condition of simultaneous load rejection of two units. The variations of the water level in the surge tank, the pressure at the pipeline end, and the velocity on the different pipeline sections with time were obtained through the calculation. The numerical results showed that the water level changing in the surge tank simulated by VOF was consistent with the field test data. These results also showed that a self-excited spiral flow occurs in the pipeline when the flow at the end of the pipeline was reduced to zero and its intensity decreased with the flow energy exhaustion. The discovery of the self-excited spiral flow in the study may provide a new explanation for the pressure wave attenuation mechanism.


Introduction
The increasing popularity of renewable energy in power systems requires large-scale energy storage technology to compensate for their intermittency [1].A pumped storage hydropower (PSH) station is an ideal peaking and emergency appliance because of its flexible operation and quick response.It is not only suitable to use as the peaking load shifting, but also has several performances and dynamic benefits such as frequency and phase modulation, emergency reserves, black startup, and load adjustment, among others [2,3].In the PSH station and diversion power station, a long pipeline is always connected to the front of the unit spiral case and the tail of draft tube [4].The flow rate of the unit varies with the working condition, which leads to dramatic rise and fall of the internal pressure in the pipeline, also known as the phenomenon of water hammer, and finally induces instability of the unit [5][6][7].Therefore, it is necessary to investigate the changing process of the water flow in a pipeline to ensure the safety of the hydropower systems [8].
Nowadays, there are three methods to research the transient process: the one-dimensional method of characteristics (MOC), the experimental method, and three-dimensional (3D) numerical simulations based on computational fluid dynamics (CFD).Among these approaches, the experimental method [9] is considered as the most reliable way to investigate the transient process.However, this method is restricted by factors such as funding, security, and test control conditions, so the advantages of this method could not be fully harnessed.The MOC approach can perform an effective prediction of water hammer characteristics, but is inadequate to predict detailed dynamic 3D characteristics such as unit variation regularity, inner flow configuration, and dynamic instability of the PSH unit [10].Currently, in order to study the transient behavior on mechanisms and the flow field of turbines, 3D CFD simulation [11] received extensive attention, as it is very useful to predict flow configuration of hydraulic machinery under normal conditions.Furthermore, the spiral flow can be found through 3D numerical method and the attenuation law of pressure wave obtained by 3D calculation is closer to the actual situation, which is helpful to analyze the transient flow characteristics of pressure wave superposition or wave peak not in the first period.
On the basis of the technical and economic requirements of power stations' construction, most PSH stations adopt one conduit system shared by multiple units, which will inevitably lead to hydraulic disturbances between units during the transient process.Chen et al. [12] performed research on the hydraulic turbulence of a high-head PSH station in China.They found that the stable area was different from one of the small load disturbances, and if the parameters of the governor were optimum for load disturbance, the remaining machines were stabilized under hydraulic disturbance.Zhang et al. [13] conducted some experiments based on a hydropower model, which combines partial diversion tunnels with tailrace tunnels, and interactions between the air and water phases were observed in the combined diversion tunnel, which would significantly alter flow dynamics.With the increasing attention to pipeline safety, the study on flow regime in pipelines has become a relevant research topic.Cheng et al. [14] applied a 3D-CFD code with the volume of fluid (VOF) model to simulate the complex free surface-pressurized transient flow and the combined transient flow under a load rejection condition.Wang et al. [15] combined the MOC and the method of implicit (MOI) to simulate an unsteady flow in a pipeline and hydropower transient processes, then the effectiveness of the coupled method was verified by simulating the water hammer in a uniform and variable area duct, and water level fluctuation in a surge tank.Xia et al. [16] conducted numerical simulations through 1D-MOC and 3D methods, respectively, of a long distance water pipe with an air-cushion surge tank, and achieved relatively consistent results.Zhou et al. [17] considered 2D and 3D VOF models to compute the valve opening process in pressure pipelines, and the results were consistent with those of the 1D calculation.Simultaneously, the results verified the advantages of the VOF model in the simulation of a gas-liquid two phase flow.Jiang et al. [18] analyzed a multi-valve protection of the water hammer in a long distance pipeline based on a 1D-MOC method in the transient process for a water supply project.Dutta et al. [19] studied the influence of Reynolds number on flow separation and reattachment in 90 • pipe blends, and the flow separation could be clearly visualized for bend with a low curvature ratio, as well as the secondary motion, which was clearly induced by the movement of fluid from the inner to outer wall of the bend, leading to flow separation.Modesto et al. [20] studied a pump as turbines (PATs) operating state in water networks from the perspective of the control valve maneuver and over speed effect, and their research characterizes the water hammer phenomenon in the design of PAT systems, emphasizing the transient events that can occur during a normal operation.Jing et al. [21] studied the pressure variations during transients in pipelines and clarified the effects of an unsteady wall shear stress on the propagation of pressure waves through the analytical and experimental method.For PSH stations, ball valves are widely used, which play an important role in controlling runaway accidents, thus some studies involving ball valve movement were reported.Wen et al. [22] performed steady calculations in different openings of the ball valve and simulated the water hammer led by linear closing of the ball valve by both 3D-CFD and 1D-MOC; they found that the unstable flow fields within the valve are essential causes of pressure and force fluctuations.Ferreira et al. [23] investigated the behavior of ball valves under steady and unsteady conditions through the experimental method, and found that the valve geometry, closure percentage, and flow regime have the largest effect on the response of the valve under the steady condition, while the valve effective closure time and pipe length would affect the unsteady characteristics of the ball valve.Martins et al. [24] analyzed the hydraulic transient flows in pressurized pipes using the CFD method, Water 2019, 11, 44 3 of 13 and the calculated velocity profiles showed two regions when the valve closure was described by a hyperbolic time-domain function.
The above studies have achieved good results in the flow simulation of diversion pipelines during the transient process, but they have not dealt with the spiral flow phenomenon in a pipeline when the load rejection process occurred.Swirl flows have a wide range of applications in various engineering areas such as chemical and mechanical mixing and separating devices, chemical reactors, combustion chambers, turbo machinery, rocketry, fusion reactors, and pollution control devices [25][26][27][28].However, public literature rarely reported the study on the spiral flow in the diversion pipeline during the transient process.In addition, the attenuation mechanism of the water hammer pressure wave in a pipeline is also the focus of some of the literature mentioned above.The attenuation rate of the pressure wave measured by the 1D method is always slower than that of the experimental results, and a good agreement cannot be achieved.There are also some scholars who use some improved models to measure the pressure wave to explain the reason that the 1D calculation results do not agree with the experimental results [29].Compared with the traditional methods, this study provides further access by considering the influence of the self-excited spiral flow on the pressure wave attenuation.On the basis of the previous research on the pipeline transient process, a 3D numerical simulation method with a VOF two-phase flow model was carried out to study and analyze the 3D internal flow characteristics of the water diversion system with an impedance surge tank during the load rejection transient process of the PSH station, and the numerical data were compared to the results from 1D-MOC.

Governing Equations
In order to accurately observe the change of water level in the surge tank, the VOF model was applied to compute the load rejection transient process of a PSH station diversion system in this paper.The continuity and momentum equations are given below: where ρ is the fluid density; t is the time; x i and x j are the coordinate components in different directions, x,j = 1,2,3; u i and u j are the velocity components in different directions, x,j = 1,2,3; p is the pressure; and µ and µ t represent the kinetic viscosity and turbulent viscosity, respectively.
ε , where C µ is a constant.

Turbulence Model
In this paper, the standard k − ε turbulence model [30] was used to numerically calculate the flow field.The turbulence kinetic energy k and specific dissipation rate ε were calculated using the following equations: where G k is the production of turbulence kinetic energy; µ t is turbulent viscosity; and the empirical constant values of the above formula are C 1ε = 1.44,C 2ε = 1.92,C µ = 0.09,σ k = 1.0, and σ ε = 1.3.
Water 2019, 11,44 For the surge chamber, the VOF model seemed to be a more feasible method to simulate the 3D surge chamber, so the VOF model was only applied to the surge chamber.The VOF method [31,32] is a surface tracking technique applied to a fixed Eulerian mesh where the Reynold average Navier-Stokes (RANS) equations, which describe the motion of the flow, have to be solved separately.The method relies on the fact that two or more fluids (or phases) are not interpenetrating.For each additional phase added to the model, a variable is introduced-the volume fraction of the phase in the computational cell.In each control volume, the volume fractions of all phases sum to a unity, that is, ∑ n i=1 α i = 1 in each cell, where i represents the phase and n is the number of phases.For air-water two phase flow, the VOF method defines α w and α a as the volume fraction of the water phase and air phase, respectively, in the same computation cell.Then, α w should satisfy the following equations: In each cell, when α w = 1, the cell contains only water; when α w = 0, the cell contains only air.For cells that contain the air-water interface, 0 < α w < 1.
From the definition of the VOF method, the fluid density in each cell is as follows: Similarly, the expression for the kinetic viscosity is as follows: For the volume fraction equation of the VOF model in this paper, the explicit scheme was used for time discretization, and the first-order upwind was applied to obtain the face fluxes for all cells.

Mesh Independence Analysis
In order to analyze the flow pattern of a water diversion system under the transient process in detail, after considering the complex geometries of the computational domain, structured hexahedral elements were chosen for tubes and unstructured tetrahedron elements were adopted for the zones except for the tube.The variation curves of the water level in the surge tank under different mesh schemes are plotted in Figure 1.Considering the accuracy and available computer resources, 2.2 million elements were used to perform the analysis and the grid element numbers of each part were shown in Table 1.Additionally, different time-step sizes of 0.1 s, 0.01 s, and 0.001 s have been used to try to simulate the transient process.After tests, the time-step size of 0.01 s was small enough to ensure the convergence.In this paper, the workstation equipped with Intel Xeon E5-1650 v3 CPU with 3.5 GHz basic frequency and 128.0 GB RAM was used to carry out the calculation.The typical CPU time required to implement the 3D simulation for the 2.2 million grid system was about seven days on this configured computer.For the 1D simulation, about 15 min was required by running the program in MATLAB software.

Equations Discretization and Boundary Conditions
In this study, ANSYS Fluent 16.0 software (Ansys2016, ANSYS, Canonsburg, PA, USA) was used to implement the calculation.Fluent is a CFD software for simulating and analyzing fluid flow problems in complex geometric regions.It provides mesh adaptive characteristics to enable users to obtain flow field solutions with high accuracy.At the same time, it provides a user programming interface, which can be secondarily developed according to solving needs.
The calculation was divided into two parts: steady and unsteady simulation.Firstly, the steady state of the flow field was achieved by steady calculation.The finite volume method was adopted to discretize the governing equations in space [33,34].A central difference scheme was used for the pressure item; a first-order upwind format was adopted to calculate the velocity, turbulent kinetic energy, and turbulent viscosity coefficient term, and a semi-implicit method for pressure linked equations-consistent (SIMPLEC) algorithm was applied to a velocity-pressure coupling solution.Then, the unsteady simulation was chosen for the load rejection transient process calculation.The control equations are discrete with the finite volume method.The pressure staggering option (PRESTO) format was employed for the pressure item; a first-order upwind format was adopted for the velocity, the turbulent kinetic energy, and the turbulent viscosity coefficient item; and a pressure implicit with splitting of operators (PISO) method was applied for the coupling solution of the velocity and pressure equations [35].
The steady calculation under design conditions was performed to provide an initial flow field for an unsteady simulation.The boundary conditions were set as follows: if initial pressure is added to the pipeline entrance directly, large errors will occur because of the large diameter of the water diversion pipe.So it is necessary to add an appropriate reservoir simplified model to the entrance according to the actual situation of the power plant.Velocity outlet was set at the diversion system outlet and the change of the flow rate at the outlet can be assumed to be a substitute for the load rejection action of the unit.Additionally, the variation of the flow rate was loaded into the end of the diversion pipeline by a user defined function (UDF) program.A wall boundary condition was applied to the walls, and roughness values were taken into account because of the route loss of the pipeline.

Equations Discretization and Boundary Conditions
In this study, ANSYS Fluent 16.0 software (Ansys2016, ANSYS, Canonsburg, PA, USA) was used to implement the calculation.Fluent is a CFD software for simulating and analyzing fluid flow problems in complex geometric regions.It provides mesh adaptive characteristics to enable users to obtain flow field solutions with high accuracy.At the same time, it provides a user programming interface, which can be secondarily developed according to solving needs.
The calculation was divided into two parts: steady and unsteady simulation.Firstly, the steady state of the flow field was achieved by steady calculation.The finite volume method was adopted to discretize the governing equations in space [33,34].A central difference scheme was used for the pressure item; a first-order upwind format was adopted to calculate the velocity, turbulent kinetic energy, and turbulent viscosity coefficient term, and a semi-implicit method for pressure linked equations-consistent (SIMPLEC) algorithm was applied to a velocity-pressure coupling solution.Then, the unsteady simulation was chosen for the load rejection transient process calculation.The control equations are discrete with the finite volume method.The pressure staggering option (PRESTO) format was employed for the pressure item; a first-order upwind format was adopted for the velocity, the turbulent kinetic energy, and the turbulent viscosity coefficient item; and a pressure implicit with splitting of operators (PISO) method was applied for the coupling solution of the velocity and pressure equations [35].
The steady calculation under design conditions was performed to provide an initial flow field for an unsteady simulation.The boundary conditions were set as follows: if initial pressure is added to the pipeline entrance directly, large errors will occur because of the large diameter of the water diversion pipe.So it is necessary to add an appropriate reservoir simplified model to the entrance according to the actual situation of the power plant.Velocity outlet was set at the diversion system outlet and the change of the flow rate at the outlet can be assumed to be a substitute for the load rejection action of the unit.Additionally, the variation of the flow rate was loaded into the end of the diversion pipeline by a user defined function (UDF) program.A wall boundary condition was applied to the walls, and roughness values were taken into account because of the route loss of the pipeline.

Case Model
Taking the diversion system of a PSH station as an example, a load rejection transient process was numerically simulated, and the outlet flow rate of the water diversion pipeline was linearly changed to zero within 30 s.The schematic diagram of the water diversion system is displayed in Figure 2.
The relevant parameters are as follows [36]: the normal storage water level and the dead water level of the upstream reservoir are 308 m and 291 m, respectively; the downstream reservoir's normal storage water level and the dead water level are 104 m and 96 m, respectively; and the ultimate dead water level is 91 m.The water diversion system [10] adopts the arrangement of one tunnel shared by two units and the whole length is 1425.39m~1499.70 m.The whole system consists of the inlet and outlet of an upper reservoir, an upper diversion adit, a surge chamber, a diversion shaft, a down diversion adit, a diversion bifurcation, and a steel lining diversion branch.The diversion tunnel diameter is 9.00 m and the branch pipe diameter is 5.60 m.The wall surface roughness of the upper diversion adit, the surge chamber, the diversion shaft, and the sown diversion adit were set to 0.015, and that of the diversion bifurcation and the steel lining diversion branch was set to 0.01.

Case Model
Taking the diversion system of a PSH station as an example, a load rejection transient process was numerically simulated, and the outlet flow rate of the water diversion pipeline was linearly changed to zero within 30 seconds.The schematic diagram of the water diversion system is displayed in Figure 2. The relevant parameters are as follows [36]: the normal storage water level and the dead water level of the upstream reservoir are 308 m and 291 m, respectively; the downstream reservoir's normal storage water level and the dead water level are 104 m and 96 m, respectively; and the ultimate dead water level is 91 m.The water diversion system [10] adopts the arrangement of one tunnel shared by two units and the whole length is 1425.39m~1499.70 m.The whole system consists of the inlet and outlet of an upper reservoir, an upper diversion adit, a surge chamber, a diversion shaft, a down diversion adit, a diversion bifurcation, and a steel lining diversion branch.The diversion tunnel diameter is 9.00 m and the branch pipe diameter is 5.60 m.The wall surface roughness of the upper diversion adit, the surge chamber, the diversion shaft, and the sown diversion adit were set to 0.015, and that of the diversion bifurcation and the steel lining diversion branch was set to 0.01.

Water Level Fluctuation in Surge Chamber
For a long diversion-type hydropower station with a surge tank, a special unsteady hydraulic phenomenon in the surge tank (such as long period, large amplitude oscillation, and slow attenuation of wave) would affect the transient performance of the hydropower unit [37].The water level fluctuation in the surge tank during the simultaneous load rejection of the two units in the PSH station with different time steps and different grid numbers after 3D numerical simulation is presented in Figure 3, as well as the comparison between experiment results and 1D-MOC results.

Water Level Fluctuation in Surge Chamber
For a long diversion-type hydropower station with a surge tank, a special unsteady hydraulic phenomenon in the surge tank (such as long period, large amplitude oscillation, and slow attenuation of wave) would affect the transient performance of the hydropower unit [37].The water level fluctuation in the surge tank during the simultaneous load rejection of the two units in the PSH station with different time steps and different grid numbers after 3D numerical simulation is presented in Figure 3, as well as the comparison between experiment results and 1D-MOC results.From Figure 3, the water level fluctuation was basically the same with different time steps and grid numbers.The overall results of the 3D numerical simulation were smaller than those of the test results when the water level reached its maximum value, but the curves were basically the same when the water level achieved its minimum value.However, the 3D numerical results were closer to the test value when the water level reached its maximum and minimum value, compared with the 1D-MOC results.
From Figure 3, the water level fluctuation was basically the same with different time steps and grid numbers.The overall results of the 3D numerical simulation were smaller than those of the test results when the water level reached its maximum value, but the curves were basically the same when the water level achieved its minimum value.However, the 3D numerical results were closer to the test value when the water level reached its maximum and minimum value, compared with the 1D-MOC results.

The Pressure Fluctuation at the End of Pipeline
The initial value of the pressure at the end of the pipeline was 305 m before load rejection.When the outlet mass flow rate of the diversion pipeline was reduced to zero, the pressure change at the end of the pipeline was worthy of more attention.The amplitude of the pressure curve using the 3D calculation was smaller than that of the 1D calculation, and the attenuation period of the 3D calculation was less than that of the 1D calculation, whereas the overall changing trend was the same (Figure 4).The maximum value of the 3D calculation exceeded the stable pressure of 6.2% and the maximum value of 1D calculation exceeded the steady pressure of 3.9%, while the 3D calculation of the lowest value was lower than the stable pressure of 2.3%, and 1D results were below its stable pressure of 3.3%.The sudden increase of the pressure at 30 s in the graph was because of a sudden decrease to zero of the man-made hypothesis.Overall, the complex flow pattern in the pipeline should be the main reason for the difference between the 1D and 3D calculation results.

The Pressure Fluctuation at the End of Pipeline
The initial value of the pressure at the end of the pipeline was 305 m before load rejection.When the outlet mass flow rate of the diversion pipeline was reduced to zero, the pressure change at the end of the pipeline was worthy of more attention.The amplitude of the pressure curve using the 3D calculation was smaller than that of the 1D calculation, and the attenuation period of the 3D calculation was less than that of the 1D calculation, whereas the overall changing trend was the same (Figure 4).The maximum value of the 3D calculation exceeded the stable pressure of 6.2% and the maximum value of 1D calculation exceeded the steady pressure of 3.9%, while the 3D calculation of the lowest value was lower than the stable pressure of 2.3%, and 1D results were below its stable pressure of 3.3%.The sudden increase of the pressure at 30 s in the graph was because of a sudden decrease to zero of the man-made hypothesis.Overall, the complex flow pattern in the pipeline should be the main reason for the difference between the 1D and 3D calculation results.

Flow Velocity Distribution at Different Cross Sections in the Pipeline
Figure 5 demonstrates the distribution of an axial and a circumferential velocity at different times and different positions in the pipeline.Where the first section, the second section, the third section, the fourth section, and the fifth section of each figure have a distance of 1 m, 60 m, 160 m, 210 m, and 400 m from the end of the diversion pipeline, and the first, second, and third sections are located between the bifurcation and the valve, and the fourth and fifth sections are located between the bifurcation and the surge tank.When the flow rate at the end of the pipeline was reduced to zero (Figure 4a), the axial velocity did not change much in the pipeline on the downstream side of the bifurcation tube, and its value was small, whereas the pipeline on the upstream side of the bifurcation pipe had a big fluctuation because of the influence of the water level variation of the surge tank; the circumferential velocity distribution showed an approximate pipe center or one-quarter radius axis  located between the bifurcation and the valve, and the fourth and fifth sections are located between the bifurcation and the surge tank.When the flow rate at the end of the pipeline was reduced to zero (Figure 4a), the axial velocity did not change much in the pipeline on the downstream side of the bifurcation tube, and its value was small, whereas the pipeline on the upstream side of the bifurcation pipe had a big fluctuation because of the influence of the water level variation of the surge tank; the circumferential velocity distribution showed an approximate pipe center or one-quarter radius axis symmetric distribution.When the circumferential velocity reached the maximum at a three-quarter radius, it declined quickly as a result of the friction of the pipe wall; the axial velocity and circumferential velocity decreased with time because of the water energy losses.

Flow Velocity Distribution at Different Cross Sections in the Pipeline
Figure 5 demonstrates the distribution of an axial and a circumferential velocity at different times and different positions in the pipeline.Where the first section, the second section, the third section, the fourth section, and the fifth section of each figure have a distance of 1 m, 60 m, 160 m, 210 m, and 400 m from the end of the diversion pipeline, and the first, second, and third sections are located between the bifurcation and the valve, and the fourth and fifth sections are located between the bifurcation and the surge tank.When the flow rate at the end of the pipeline was reduced to zero (Figure 4a), the axial velocity did not change much in the pipeline on the downstream side of the bifurcation tube, and its value was small, whereas the pipeline on the upstream side of the bifurcation pipe had a big fluctuation because of the influence of the water level variation of the surge tank; the circumferential velocity distribution showed an approximate pipe center or one-quarter radius axis symmetric distribution.When the circumferential velocity reached the maximum at a three-quarter radius, it declined quickly as a result of the friction of the pipe wall; the axial velocity and circumferential velocity decreased with time because of the water energy losses.
Water 2018, 10, x FOR PEER REVIEW 9 of 14

Flow Regime Analysis in the Bifurcated Pipeline
The internal flow pattern in the bifurcated pipeline is provided in Figure 6 and Figure 7.The spiral flow in the upstream main pipeline is more obvious than that in the downstream branch pipe.For the two branches in the downstream, some obvious rotational flows were observed in branch pipe A, while branch pipe B had no obvious similar flow, because of the different connection angles for the two branches with the main pipeline (Figure 6).However, from Figure 6, when the flow rate

Flow Regime Analysis in the Bifurcated Pipeline
The internal flow pattern in the bifurcated pipeline is provided in Figures 6 and 7.The spiral flow in the upstream main pipeline is more obvious than that in the downstream branch pipe.For the two branches in the downstream, some obvious rotational flows were observed in branch pipe A, while branch pipe B had no obvious similar flow, because of the different connection angles for the two branches with the main pipeline (Figure 6).However, from Figure 6, when the flow rate at the end of the diversion pipe was reduced to zero, that is, t = 30 s, the flow patterns in the two pipes were almost the same swirling flow.When t = 50 s, the original weak swirling flow structure in branch B was destroyed by the centrifugal force of the spiral flow in the main pipeline and branch A, and a vortex was formed when some water flows in from one side of branch B and some water flows out from the other side.After load rejection, as time increases, the water flow produced shear deformation due to the centrifugal force and the irregular boundaries when flowing through the elbow, bifurcation pipe, and other flow passage components, which induced a complex pipeline spiral linear flow, as shown in Figure 8.The strength of shear deformation decreased with time, then the particle flow spiral route length become longer to accelerate the energy attenuation of the water flow.After load rejection, as time increases, the water flow produced shear deformation due to the centrifugal force and the irregular boundaries when flowing through the elbow, bifurcation pipe, and other flow passage components, which induced a complex pipeline spiral linear flow, as shown in Figure 8.The strength of shear deformation decreased with time, then the particle flow spiral route length become longer to accelerate the energy attenuation of the water flow.

Conclusions
The 3D numerical simulation of the long water diversion system with a surge tank during the load rejection process of the PSH station was carried out using a VOF two-phase flow model, and it can intuitively reflect the change of water level in a surge chamber at any time.The numerical results of the highest and lowest surge water level showed good agreement with the experimental data.
The variation pattern of the pressure at the end of the water diversion pipeline was obviously different between the 1D-MOC results and the 3D-CFD numerical simulation results.The pressure fluctuation rule obtained by the 3D numerical simulation method was more consistent with the actual situation in attenuation amplitude and attenuation period because the 3D numerical simulation method can simulate the actual flow state in the pipeline more realistically.The complex spiral flow in the water diversion pipeline was the internal cause of the difference between the 3D and 1D calculations.When the flux at the end of the water diversion pipe was reduced to zero, the upstream side of the whole pipeline system was affected by the fluctuation of the surge tank water level, and the axial velocity of the upstream bifurcation pipeline was larger than the downstream.The circumferential velocity was reduced quickly because of the friction of the pipe wall.The existence of spiral flow makes the actual travel distance of the flow in the pipe longer, so that the friction distance between the flow and the pipe wall is longer.Because of the circumferential velocity of water flow, the friction between the water flow in the pipe and the pipe wall was intensified, which accelerates the attenuation of water flow energy.
Future studies could also consider the ball valve movement during the transient process, which may also have an impact on the results.Furthermore, the calculation conditions should be increased for different pipeline layout modes.In the future work, the whole flow system including the pumpturbine unit should be established, so that more precise and comprehensive flow characteristics can be investigated during the transient process of the complex PSH station.
Author Contributions: D.Z. suggested the methodology and performed critical revision of the article; H.C. performed the numerical simulation, analyzed the numerical results, and wrote the paper; S.C. designed the numerical methods and performed the calculations.

Conclusions
The 3D numerical simulation of the long water diversion system with a surge tank during the load rejection process of the PSH station was carried out using a VOF two-phase flow model, and it can intuitively reflect the change of water level in a surge chamber at any time.The numerical results of the highest and lowest surge water level showed good agreement with the experimental data.
The variation pattern of the pressure at the end of the water diversion pipeline was obviously different between the 1D-MOC results and the 3D-CFD numerical simulation results.The pressure fluctuation rule obtained by the 3D numerical simulation method was more consistent with the actual situation in attenuation amplitude and attenuation period because the 3D numerical simulation method can simulate the actual flow state in the pipeline more realistically.The complex spiral flow in the water diversion pipeline was the internal cause of the difference between the 3D and 1D calculations.When the flux at the end of the water diversion pipe was reduced to zero, the upstream side of the whole pipeline system was affected by the fluctuation of the surge tank water level, and the axial velocity of the upstream bifurcation pipeline was larger than the downstream.The circumferential velocity was reduced quickly because of the friction of the pipe wall.The existence of spiral flow makes the actual travel distance of the flow in the pipe longer, so that the friction distance between the flow and the pipe wall is longer.Because of the circumferential velocity of water flow, the friction between the water flow in the pipe and the pipe wall was intensified, which accelerates the attenuation of water flow energy.
Future studies could also consider the ball valve movement during the transient process, which may also have an impact on the results.Furthermore, the calculation conditions should be increased for different pipeline layout modes.In the future work, the whole flow system including the pump-turbine unit should be established, so that more precise and comprehensive flow characteristics can be investigated during the transient process of the complex PSH station.and the Fundamental Research Funds for the Central Universities (2016B41514).The support of Hohai University and Wuhan University, China is also gratefully acknowledged.

Conflicts of Interest:
All the authors declare that there is no conflict of interests regarding the publication of this paper.All the authors do not have a direct financial relation with the commercial identities mentioned in this paper that might lead to a conflict of interests for any of the authors.

Figure 1 .
Figure 1.Test results of different mesh schemes and time step sizes.

Figure 2 .
Figure 2. The schematic diagram of the diversion system of a pumped storage hydropower (PSH) station.

Figure 2 .
Figure 2. The schematic diagram of the diversion system of a pumped storage hydropower (PSH) station.

Figure 3 .
Figure 3. Water fluctuation in a surge tank (where Δt is time step).

Figure 3 .
Figure 3. Water fluctuation in a surge tank (where ∆t is time step).

Water 2018 ,
10, x FOR PEER REVIEW 8 of 14

Figure 4 .
Figure 4.The curve of static pressure at the end of pipeline.CFD-computational fluid dynamics; MOC-method of characteristics.

Figure 4 .
Figure 4.The curve of static pressure at the end of pipeline.CFD-computational fluid dynamics; MOC-method of characteristics.

Figure 5
Figure 5 demonstrates the distribution of an axial and a circumferential velocity at different times and different positions in the pipeline.Where the first section, the second section, the third section, the fourth section, and the fifth section of each figure have a distance of 1 m, 60 m, 160 m, 210 m, and 400 m from the end of the diversion pipeline, and the first, second, and third sections are

Water 2018 ,Figure 4 .
Figure 4.The curve of static pressure at the end of pipeline.CFD-computational fluid dynamics; MOC-method of characteristics.

Figure 5 .
Figure 5. Flow velocity distribution of each section at different times: (a) the axial velocity of different sections when t = 50 s, (b) the circumferential velocity of different sections when t = 50 s, (c) the axial velocity of different sections when t = 90 s, (d) the circumferential velocity of different sections when t = 90 s, (e) the axial velocity of different sections when t = 130 s, (f) the circumferential velocity of different sections when t = 130 s.

Figure 5 .
Figure 5. Flow velocity distribution of each section at different times: (a) the axial velocity of different sections when t = 50 s, (b) the circumferential velocity of different sections when t = 50 s, (c) the axial velocity of different sections when t = 90 s, (d) the circumferential velocity of different sections when t = 90 s, (e) the axial velocity of different sections when t = 130 s, (f) the circumferential velocity of different sections when t = 130 s.

Figure 7 .
Figure 7. Water streamline in a bifurcation pipe at t = 30 s.

Figure 7 .
Figure 7. Water streamline in a bifurcation pipe at t = 30 s.

Figure 7 .
Figure 7. Water streamline in a bifurcation pipe at t = 30 s.After load rejection, as time increases, the water flow produced shear deformation due to the centrifugal force and the irregular boundaries when flowing through the elbow, bifurcation pipe, and other flow passage components, which induced a complex pipeline spiral linear flow, as shown in Figure8.The strength of shear deformation decreased with time, then the particle flow spiral route length become longer to accelerate the energy attenuation of the water flow.

:
The study was supported by the Key Program of National Science Foundation of China (No. 51839008; No.51806058; No. 51579080; No.51339005).The study was also supported by the Open Research Fund Program of State Key Laboratory of Water Resource and Hydropower Engineering Science (2015SDG04)

Table 1 .
Mesh number (N) and type of each part.
Figure 1.Test results of different mesh schemes and time step sizes.

Table 1 .
Mesh number (N) and type of each part.Passages Surge tank Bifurcated pipe Diversion pipe Other region Total