Transient Process of Pumped Storage System Coupling Gas–Liquid Interface: Novel Mathematical Model and Experimental Veriﬁcation

: The traditional calculation method for a transient process has high accuracy when the pipeline only contains liquid, but when the pipeline contains both gas and liquid the accuracy is greatly reduced. The coupling characteristics of gas–liquid interface movement in hydraulic transient processes are not clear due to the lack of high-precision mathematical model and experimental veriﬁcation. This paper proposes a novel mathematical model of a gas–liquid pipeline system in a hydropower station based on Preissman’s implicit difference scheme and the method of characteristics. The solving mechanism of the transient process of gas–liquid movement was developed on the gas–liquid interface tracking method. Subsequently, the models proposed in this paper were applied in two typical scenarios of a gas–liquid transient process in a hydropower system, and their accuracy were veriﬁed in a ﬁeld experiment. The comparison results showed that the novel model could accurately capture the movement of the gas–liquid interface, and the average relative error of the characteristic parameter was about 7.2%. Under the load rejection condition, the change speed of characteristic parameters was positively correlated with the pipeline slope. Under the pump failure after low-head startup condition, the maximum pumping discharge was negatively correlated with startup water level and the maximum reversal discharge and speed were positively correlated with the pump failure water level. Compared with the conventional method, the proposed model has advantages in solving the complex transient process coupling gas–liquid. It has potential value in applications such as the safe operation of hydropower stations, the transient process of water diversion projects and in urban pipe network operation.


Introduction
With the low-carbon transition of the global energy, the operation mode of hydropower has changed, and its flexible regulation ability has attracted more and more attention. In particular, pumped storage, which can function in two-way modes of the turbine and pump, is an important means to support the development of renewable energy. The frequent start and stop and the conversion of operating conditions of pumped storage increase the possibility of extreme operating conditions. The study of the hydropower transient process is a key way to ensure the safe and stable operation of the unit. As a large number of pumped storage power stations are put into operation, it can be found that the pipeline no longer only contains liquid during the transient process, but may contain both gas and liquid. The current transient process simulation method of a pumped storage power station can simulate relatively accurate results when the pipeline only contains liquid. When the pipeline has a gas-liquid two-phase flow, the current transient process simu-whole pipeline system of a pumped storage power station. For the pipeline water filling simulation, the boundary conditions of pumps-turbines and intake gate are not involved, hence its results cannot be referenced directly.
To determine the transient process of a gas-liquid two-phase flow in a pipeline, this paper proposes a novel mathematical model of the transient process of a gas-liquid two-phase flow in a pipeline based on Preissman's implicit difference scheme and the MOC. The recurrence equation of the gas-liquid coupling surface is derived based on the interface tracing method and the treatment method of corresponding boundary conditions of gas and liquid pipelines are given. The applicability and accuracy of numerical simulation as proposed in this work have been verified by comparing with experimental results. Finally, two typical scenarios of a gas-liquid transient process under different conditions are analyzed.

Basic Equation and Solution Method for Unsteady Gas Flow
For the case of an unsteady flow in a gas pipeline, the basic equations include the momentum equation and the continuity equation; their expressions can be written as follows [26]: where the relevant parameters are defined as follows: B is the wave velocity (m/s); A is the sectional area of the pipeline (m 2 ); M is the mass flow (kg/s); x is the position along the axis of the pipeline (m); p is the absolute pressure of gas (m); t is the time (s); α is the inertia factor; g is the gravitational acceleration (m 2 /s); θ is the included angle between the axis of the pipeline and the horizontal plane; f is the Darcy-Weisbach coefficient of friction resistance; D is the diameter of the pipeline (m), D = 2A/S, where S is the pipeline section perimeter (m). It should be noted that: where v is flow velocity (m/s).

2.
The wave velocity can be determined by using the state equation B = p/ρ = √ λRT; where ρ is the gas density (kg/m 3 ), λ is the compressibility coefficient, R is the gas constant (m 2 /(s 2 ·K)) and T is the absolute temperature (K). 3.
When the axis rises up along the direction of +x, θ, the included angle between the axis of the pipeline and the horizontal plane, takes a positive value.
Equations (1) and (2) are discretized according to Preissmann's four-point finitedifference implicit scheme, and the mass flow and gas pressure are expressed in increments. The discrete control equations are obtained as shown in Equations (4) and (5): where: For gas pipe sections solved by the implicit method, the following relationships exist among the internal calculated sections [28]: where: In the equation, EE i and FF i are the transfer coefficients and only depend on the EE i−1 and FF i−1 of the front section, while L i , R i , and N i depend on the transfer coefficients EE i and FF i . Through the boundary conditions of the first section surface, the EE 1 and FF 1 of the first node can be obtained, so the EE i , FF i , L i , R i and N i of each section can be obtained from Equations (6) and (7). This process is called forward scanning. The value of ∆p n and ∆M n can be obtained by incorporating Equation (6) and the boundary condition equation of the end node. The parameter ∆p n−1 is obtained from Equation (7) and subsequently ∆M n−1 is obtained from Equation (6). Following such recursions, the pressure of gas and mass flow of all the sections can be obtained. This process is called backward scanning. Considering the elasticity of the water body and pipe wall, the basic equation of  unsteady flow under pressure includes a continuity equation and a momentum equation;  their expressions can be written as

Basic Equation and Solution Method for Unsteady Liquid Flow
where x is the distance (m); t is the time (s); Q is the discharge (m 3 /s); H is the pressure head (m); a is the wave speed (m/s); A is the cross-sectional area (m 2 ); β is the pipe slope; g is the gravitational acceleration (m 2 /s); f is the Darcy-Weisbach friction factor; and D is the inner diameter of the pipe (m), D = 2A/S, where S is the pipeline section perimeter (m). When the pipeline is prismatic, ∂A/∂x = 0. For the liquid section in the gas-liquid pipeline, Equations (8) and (9) are discretized according to the four-point implicit scheme difference scheme in Preissman space, and the discharge and pressure head are expressed in increments. The discrete control equations are obtained as shown in Equations (10) and (11): where: For the liquid pipe section solved by the implicit scheme method, the internal calculation sections also have the following relationship: where: The solution method of the liquid section in a gas-liquid pipeline is the same as that of a gas pipeline, so it is not repeated here.
For liquid pipelines not in contact with gas pipelines, the MOC approach transforms the quasi-linear hyperbolic partial differential equation into two groups of ordinary dif-Water 2021, 13, 2933 6 of 23 ferential equations on two characteristic lines, as shown in Figure 2 and as expressed in Equations (14) and (15). where: Water 2021, 13, x FOR PEER REVIEW 7 of 26 where: The parameters HP and QP denote the unknown pressure head and unknown discharge at point P at the time t + Δt, respectively, and the parameters QR, QS, HR, and HS denote the pressure head and discharge at points R and S, respectively, at the time t. Their values can be interpolated from adjacent grid nodes. Therefore, at time t + Δt, CP, BP, CM, and BM can be obtained from Equations (14) and (15). From these two equations, the unknown pressure head and discharge at any node in the pipeline can be obtained.

Gas-liquid Coupling Interface and Recursive Equations
The interface tracking method is used to calculate and judge the position of the gasliquid coupling surface at every moment in the simulation calculation and establish corresponding control equations for different flow regions. For each pipe in the calculation model, it is divided into N ( ) x a t Δ = Δ pipe segments and N + 1 sections. If the elevation of the first and end sections of a pipeline are above the water level, the pipeline is a gas pipeline. If the elevation of the first and end sections of a pipeline are below the water level, the pipeline is a liquid pipeline. If the water level is between the elevation of the first and end sections of a pipeline, the pipeline is a gas-liquid pipeline. The schematic diagram of pipeline segmentation is shown in Figure 3. In the Figure 3, yellow represents gas and blue represents liquid. If it is a gas pipeline or a liquid pipeline, the pipeline division remains unchanged and is solved according to their respective control equations. The gas pipeline is solved by the implicit difference method and the liquid pipeline is solved by the MOC method. The point boundary between pipelines is calculated according to the The parameters H P and Q P denote the unknown pressure head and unknown discharge at point P at the time t + ∆t, respectively, and the parameters Q R , Q S , H R , and H S denote the pressure head and discharge at points R and S, respectively, at the time t. Their values can be interpolated from adjacent grid nodes. Therefore, at time t + ∆t, C P , B P , C M , and B M can be obtained from Equations (14) and (15). From these two equations, the unknown pressure head and discharge at any node in the pipeline can be obtained.

Gas-liquid Coupling Interface and Recursive Equations
The interface tracking method is used to calculate and judge the position of the gas-liquid coupling surface at every moment in the simulation calculation and establish corresponding control equations for different flow regions. For each pipe in the calculation model, it is divided into N (∆x = a∆t) pipe segments and N + 1 sections. If the elevation of the first and end sections of a pipeline are above the water level, the pipeline is a gas pipeline. If the elevation of the first and end sections of a pipeline are below the water level, the pipeline is a liquid pipeline. If the water level is between the elevation of the first and end sections of a pipeline, the pipeline is a gas-liquid pipeline. The schematic diagram of pipeline segmentation is shown in Figure 3. In the Figure 3, yellow represents gas and blue represents liquid. If it is a gas pipeline or a liquid pipeline, the pipeline division remains unchanged and is solved according to their respective control equations. The gas pipeline is solved by the implicit difference method and the liquid pipeline is solved by the MOC method. The point boundary between pipelines is calculated according to the boundary equation.
In case of a gas-liquid pipeline, the position of the gas-liquid coupling surface shall be judged according to the gas-liquid coupling surface calculated at the previous time and the elevation of each section. When (Z − Z k )(Z − Z k+1 ) < 0, the gas-liquid coupling surface is located between the k section and k + 1 section of the pipeline. The liquid pipe section is the section from the Z 0 section to the Z section, totaling k + 1 segments, where the length of k segments is ∆x and the length of a segment (i.e., the Z k section to the Z section) is x L,Z = ∆x(Z − Z k /Z k+1 − Z k ). The gas pipe section is the section from the Z section to the Z n section, totaling N−k segments, where the length of N−k−1 segments is ∆x and the Water 2021, 13, 2933 7 of 23 length of a segment (i.e., the Z section to the Z k+1 section) is The schematic diagram of a gas-liquid pipeline is shown in Figure 4. In the Figure 4, yellow represents gas and blue represents liquid. In case of a gas-liquid pipeline, the position of the gas-liquid coupling surface shall be judged according to the gas-liquid coupling surface calculated at the previous time and the elevation of each section. When ( )( ) surface is located between the k section and k+1 section of the pipeline. The liquid pipe section is the section from the Z0 section to the Z section, totaling k + 1 segments, where the length of k segments is Δx and the length of a segment (i.e., the Zk section to the Z section) is L, The gas pipe section is the section from the Z section to the Zn section, totaling N-k segments, where the length of N-k-1 segments is Δx and the length of a segment (i.e., the Z section to the Zk+1 section) is The schematic diagram of a gas-liquid pipeline is shown in Figure 4. In the Figure 4, yellow represents gas and blue represents liquid. The recurrence equation is established by taking the gas-liquid contact surface as the first boundary of the gas pipe section, the gas-liquid pipe and the gas pipeline contact surface or the atmospheric section as the final boundary. At the gas-liquid coupling interface, the liquid moves at the same speed as the gas. When LG Z Q Q = and the gas mass flow can be written as Equation (16): The relationship between the head of the liquid pressure measuring pipe and the gas pressure can be expressed as follows:  In case of a gas-liquid pipeline, the position of the gas-liquid coupling surface shall be judged according to the gas-liquid coupling surface calculated at the previous time and the elevation of each section. When ( )( ) surface is located between the k section and k+1 section of the pipeline. The liquid pipe section is the section from the Z0 section to the Z section, totaling k + 1 segments, where the length of k segments is Δx and the length of a segment (i.e., the Zk section to the Z section) is L, The gas pipe section is the section from the Z section to the Zn section, totaling N-k segments, where the length of N-k-1 segments is Δx and the length of a segment (i.e., the Z section to the Zk+1 section) is The schematic diagram of a gas-liquid pipeline is shown in Figure 4. In the Figure 4, yellow represents gas and blue represents liquid. The recurrence equation is established by taking the gas-liquid contact surface as the first boundary of the gas pipe section, the gas-liquid pipe and the gas pipeline contact surface or the atmospheric section as the final boundary. At the gas-liquid coupling interface, the liquid moves at the same speed as the gas. When LG Z Q Q = and the gas mass flow can be written as Equation (16): The relationship between the head of the liquid pressure measuring pipe and the gas pressure can be expressed as follows: The recurrence equation is established by taking the gas-liquid contact surface as the first boundary of the gas pipe section, the gas-liquid pipe and the gas pipeline contact surface or the atmospheric section as the final boundary. At the gas-liquid coupling interface, the liquid moves at the same speed as the gas. When A GL,Z = A LG,Z , then Q GL,Z = Q LG,Z and the gas mass flow can be written as Equation (16): The relationship between the head of the liquid pressure measuring pipe and the gas pressure can be expressed as follows: The x 0 G,Z at the last time step has the following relationship with the x G,Z at the current time step: Bringing (18) and the variation of interface position of the gas-liquid coupling can be expressed as follows: Water 2021, 13, 2933 . Superscript 0 denotes the value at the last time step.
The gas-liquid contact surface equation of the liquid part can be expressed as follows: Linearizing Equations (16) and (17) can get Equations (21) and (22): Bringing Equations (19), (21), and (22) into Equation (20) and the recursive equation of the gas-liquid coupling surface of the gas part can be obtained, as shown in Equation (23): Simplify Equation (23) to get Equation (24): where:

Unit Load Rejection
The governing equations of the pump turbine load rejection include (1) the boundary equation of the pipelines linked to the machine set, (2) the electric generator equation, (3) the equations for the machine set unit parameters, (4) the equation for the turbine characteristics curve, and (5) the speed governor equation. These equations together include a total of nine unknown values, namely, H P , H S , Q P , Q 1 , n 1 , M 1 , n, M t , and y.
Equation C + at the turbine inlet and equation C − at the turbine outlet can be expressed as follows: The electric generator first-order equation can be expressed as follows: The unit parameter equation of the turbine can be expressed as follows: The characteristic curve equation of the turbine can be expressed as follows: Water 2021, 13, 2933 9 of 23 The guide vane specified equation can be expressed as follows: where H P and H S denote the pressure head of the turbine inlet and outlet side, respectively, (m); Q P denotes the discharge of turbine (m 3 /s); M t denotes the turbine output,(N·m); M g denotes the generating torque (N·m); GD 2 denotes the rotary inertia (kg·m 2 ); D 1 denotes the diameter of the turbine (m); n 1 denotes the unit speed (rpm); Q 1 denotes the unit discharge (m 3 /s); M 1 denotes the unit torque (N·m); n denotes the revolving speed (rpm); denotes the opening of the guide vane at a given value. Subscript 0 denotes the value at the last time step.

The Intake Gate Closure
The rapid closing of the upstream intake gate is a complex gas-liquid two-phase flow process. To accurately simulate this process, the process of gate closure is divided into three stages, as shown in Figure 5. In the Figure 5, yellow represents gas, blue represents liquid and red represents the intake gate. 1. The intake gate is fully opened.
When the intake gate is fully opened, the governing equations include the equation of orifice outflow and the C − equation.
The equation of orifice outflow can be expressed as follows: The C − equation behind the gate inlet can be expressed as follows: where U Z is the upper reservoir water level (m); D Z is the down reservoir water level, A is the area of pipeline (m 2 ); φ is the discharge coefficient of gate orifice discharge, P H is the pressure head of the inlet (m); P Q is the discharge of the inlet (m 3 /s).
2. The intake gate is closing.
When the intake gate starts to close, the air vent will start to intake air. Assuming that the airbag only develops downstream, the governing equations of the gate closure include Where Hp is the pressure head behind the gate (m); 1 P Q is the discharge behind the gate (m 3 /s).
The C − equation behind the gate inlet can be expressed as follows:

1.
The intake gate is fully opened.
When the intake gate is fully opened, the governing equations include the equation of orifice outflow and the C − equation.
The equation of orifice outflow can be expressed as follows: The C − equation behind the gate inlet can be expressed as follows: where Z U is the upper reservoir water level (m); Z D is the down reservoir water level, A is the area of pipeline (m 2 ); ϕ is the discharge coefficient of gate orifice discharge, H P is the pressure head of the inlet (m); Q P is the discharge of the inlet (m 3 /s).

2.
The intake gate is closing.
When the intake gate starts to close, the air vent will start to intake air. Assuming that the airbag only develops downstream, the governing equations of the gate closure include

m.
The equation of the orifice outflow can be expressed as follows: where H p is the pressure head behind the gate (m); Q P1 is the discharge behind the gate (m 3 /s). The C − equation behind the gate inlet can be expressed as follows: where Q P2 is the discharge behind airbag (m 3 /s). The equation of the gas state can be expressed as follows: where p is the gas pressure (m); V is the gas volume (m 3 ); m is the gas mass (kg); M is the molar mass of gas (kg/mol); R is the gas constant (m 2 /(s 2 ·K)); T is the absolute temperature (K). The relationship equation between the piezometer head H f and the absolute pressure p in the tube can be expressed as follows: where H a is the atmospheric pressure head (m); γ is the liquid specific weight (kg/m 3 ); Z k is the bottom elevation of the venthole (m). Depending on the air inlet and outlet venthole velocity, the equation of the venthole mass gas flow can be divided into the following four cases: • Subsonic air flow in: where C in is the venthole discharge coefficient; A in is the area of venthole (m 2 ); ρ a is the mass density of atmospheric air ρ a = p a M/RT a .
The inlet gate is fully closed.
When the intake gate is fully closed, Q P1 = 0 and the gas-liquid coupling interface gradually moves downstream. The vent hole is taken as the end section of the gas pipeline (or gas-liquid pipeline) and the mass flow of the section can be calculated from Equations (38)-(41). The air pressure change value can be calculated from the end section equation of the forward scanning and can be expressed as follows:

Connecting Pipes
1. Connecting pipes between the liquid pipeline and the gas-liquid pipeline The junctions of pipes in a series meet the continuous flow condition. Meanwhile, the pressure of the sections just before and after the junction should be treated the same, if the local loss of head is neglected. By linearizing the C + equations, the transfer equation of the connecting pipe between liquid pipe and gas-liquid pipe can be expressed as follows: where EE LG,N = C QM,0 , FF LG,N = −Q 0 CM,0 + Q CM,0 . Superscript 0 denotes the value at the last time step.
2. Connecting pipes between the gas pipeline and the gas-liquid pipeline Assuming that the areas of the series points are same, the transfer equation at the series points can be expressed as follows: 3. Contact surface of the gas pipeline with the atmosphere The pressure of the contact surface connecting the gas pipeline with the atmosphere is always atmospheric pressure, i.e., p G,0 = p 0 and ∆p G,0 = 0 By using Equation (7) the unknown quantity ∆M G of this boundary node can be obtained as follows:

Analysis of Numerical Simulation
In accordance with the mathematical model developed in Section 2 and the boundary conditions from Section 3, the simulation of a transient process of a gas-liquid two-phase flow in a pumped storage power station can be carried out. The steps involved in the computational procedure are as follows: (1) Divide the water conveyance system into several pipelines, and each pipeline is divided into several sections. The space steps and the time steps can be denoted as ∆x and ∆t, respectively, where ∆x = a∆t. (2) According to the water level at time t = t 0 and the elevation of the head and end sections of each pipeline, the pipeline is classified into a gas, liquid, or gas-liquid pipeline, and the position of the gas-liquid coupling surface is determined. (3) On the basis of the initial condition and boundary condition at time t = t 0 , the H i and Q i at time t + ∆t of each section can be determined by using Equation (14) and (15). (4) The contact surface between the gasliquid pipeline and the liquid pipeline is used as the first section to perform the forward scanning to the liquid section of the gas-liquid pipeline until the gas-liquid coupling surface, and then the gas-liquid coupling surface is used as the first section to perform the forward scanning to the gas section of the gas-liquid pipeline and gas pipeline until the atmosphere. (5) The ∆M N and ∆p N of the end section of the gas pipeline can be determined by using the boundary conditions of the atmosphere. From the end section of the gas pipeline to the gas-liquid coupling surface, the backward scanning for gas is performed by using Equations (6) and (7) and the M i and p i at time t + ∆t of each section can be determined. Then, the ∆Z can be determined by using Equation (18). (6) From the gas-liquid coupling surface to the contact surface between the gas-liquid pipeline and the liquid pipeline, the backward scanning for liquid is performed by using Equations (12) and (13) and the H i and Q i at time t + ∆t of each section can be determined. (7) When the variable value on each pipeline section is determined at time t + ∆t, the variable value on each pipeline section at time t + 2∆t can be calculated, and so on until the required time. The complete simulation process is shown in Figure 6.

Model Verification
To verify the accuracy of the simulation method and simulation program for the transient process of a gas-liquid two-phase flow in a pumped storage power station based on the gas-liquid interface coupling proposed in this paper, two typical working conditions were selected for a comparison between the simulation and experiment results. The field diagram of an experimental platform is shown in Figure 6 and the two typical working conditions are: (A) When the upstream tank water level is 14.5 m and downstream tank water level is 3.85 m, the unit rejects its load during stable operation, and the unit guide vane and ball valve fail to close. Meanwhile, the intake gate is closed within 1.5 s. (B) When the water level in the upstream pipeline is flush with the downstream reservoir, the unit starts under the pumping condition, and the unit guide vane opens within 5 s. After 32 s, the pump rejects its load and the guide vane fails to close.

Model Verification
To verify the accuracy of the simulation method and simulation program for the transient process of a gas-liquid two-phase flow in a pumped storage power station based on the gas-liquid interface coupling proposed in this paper, two typical working conditions were selected for a comparison between the simulation and experiment results. The field diagram of an experimental platform is shown in Figure 6 and the two typical working conditions are: (A) When the upstream tank water level is 14.5 m and downstream tank water level is 3.85 m, the unit rejects its load during stable operation, and the unit guide vane and ball valve fail to close. Meanwhile, the intake gate is closed within 1.5 s. (B) When the water level in the upstream pipeline is flush with the downstream reservoir, the unit starts under the pumping condition, and the unit guide vane opens within 5 s. After 32 s, the pump rejects its load and the guide vane fails to close.
The experimental platform consists of five parts: upstream tank, downstream tank, pump turbine unit, penstocks, and measurement system. The field diagram of experimental platform is shown in Figure 7. The upstream and downstream tanks adjust the water level by adjusting the height of the water-stabilizing slot. When the water level is higher than the water-stabilizing slot, excess water will be discharged from the drainage pipe in the middle of the water-stabilizing slot. The upstream water tank is an open tank. By adjusting the height of the water-stabilizing slot, the water level of the upstream water can be stabilized between 14.25 m and 14.75 m. An electric gate capable of closing quickly was installed at the inlet of the upstream pipeline, and three air vents with a diameter of 5 cm were set after the inlet. The air vent allows air to enter the pipeline after the intake gate is closed. The downstream tank is also an open tank, but the water-stabilizing slot of the downstream tank cannot be adjusted and the downstream water level under steady The experimental platform consists of five parts: upstream tank, downstream tank, pump turbine unit, penstocks, and measurement system. The field diagram of experimental platform is shown in Figure 7. The upstream and downstream tanks adjust the water level by adjusting the height of the water-stabilizing slot. When the water level is higher than the water-stabilizing slot, excess water will be discharged from the drainage pipe in the middle of the water-stabilizing slot. The upstream water tank is an open tank. By adjusting the height of the water-stabilizing slot, the water level of the upstream water can be stabilized between 14.25 m and 14.75 m. An electric gate capable of closing quickly was installed at the inlet of the upstream pipeline, and three air vents with a diameter of 5 cm were set after the inlet. The air vent allows air to enter the pipeline after the intake gate is closed. The downstream tank is also an open tank, but the water-stabilizing slot of the downstream tank cannot be adjusted and the downstream water level under steady state is 3.85 m. The pump turbine runner diameter is 0.28 m, rated speed is 1000 r/min, and installation elevation is 1.18 m. The penstocks consist of an upstream pipe and a downstream pipe. The total length of upstream piping is about 31.3 m and that of downstream piping is about 26.6 m. Detailed information of each pipeline segment is shown in Table 1. Six pressure sensors were arranged in the upstream pipe of the experimental platform to measure the change of water head which can reflect the change of water level. Pressure sensors were also installed at the inlet of the spiral case and draft tube to measure the changes of spiral case inlet pressure and draft tube inlet pressure. An electromagnetic flowmeter was installed in front of the spiral case to measure the discharge of the unit and a speed sensor was installed on the unit to measure the unit speed. The schematic diagram of the experimental platform is shown in Figure 8.     stream piping is about 26.6 m. Detailed information of each pipeline segment is shown in Table 1. Six pressure sensors were arranged in the upstream pipe of the experimental platform to measure the change of water head which can reflect the change of water level. Pressure sensors were also installed at the inlet of the spiral case and draft tube to measure the changes of spiral case inlet pressure and draft tube inlet pressure. An electromagnetic flowmeter was installed in front of the spiral case to measure the discharge of the unit and a speed sensor was installed on the unit to measure the unit speed. The schematic diagram of the experimental platform is shown in Figure 8.    The simulation region was the whole experimental platform, and the calculation conditions were the same as the experimental ones. The comparison between the numerical simulation and experimental data of the water head at each monitoring point of the pipeline and the characteristic parameters in working condition A are shown in Figures 9 and 10.
The experimental data show that when the unit load rejection and guide vane of the unit fails to close, closing the upstream intake gate can prevent the unit from entering the runaway operations. After the intake gate is closed, the water in the upstream reservoir will no longer enter the pipeline, and the air enters the pipeline through the air vent. Meanwhile, with the continuous operation of the unit, the upstream water level is getting lower, and the working head of the unit is decreasing, resulting in the continuous reduction of unit speed and discharge. The change rate of the water head at each measuring point of the pipeline, and the inlet pressure of the spiral case of the unit, the unit discharge and the unit speed depend on the layout of the pipeline where the gas-liquid coupling interface is located. When the pipeline is horizontal, the change rate decreases, and when the pipeline is inclined, the change rate increases. By comparing the change trend of the pressure at each monitoring point with the pressure at the inlet of the spiral case, it can be seen that the pressure of the spiral case can reflect the change of the gas-liquid coupling surface in the upstream pipeline. The simulation region was the whole experimental platform, and the calculation conditions were the same as the experimental ones. The comparison between the numerical simulation and experimental data of the water head at each monitoring point of the pipeline and the characteristic parameters in working condition A are shown in Figures 9 and  10.
(e) (f) The experimental data show that when the unit load rejection and guide vane of the unit fails to close, closing the upstream intake gate can prevent the unit from entering the runaway operations. After the intake gate is closed, the water in the upstream reservoir will no longer enter the pipeline, and the air enters the pipeline through the air vent. Meanwhile, with the continuous operation of the unit, the upstream water level is getting level is slowly decreasing, but in the simulation, the water level will drop sharply at the end of the gate closing. As a result, the duration of fluctuation increases.
To further verify the accuracy of the solution method proposed in this paper, the condition B was simulated. The comparison of experimental data and simulation results of characteristic parameters under condition B are shown in Figure 11 and the comparison of extreme values of characteristic parameters is shown in  The experimental data show that after the guide vane is opened, the discharge of the unit rises rapidly and then decreases with the increase of water level in the upstream pipeline after reaching the maximum value. The pressure of the spiral case decreases rapidly within seconds after starting and then rises with the rise of water level in the upstream pipeline. The inlet pressure of the draft tube decreases due to the increase of discharge and then rises slightly due to the decrease of discharge. Because the variation trend of upstream water level depends on the arrangement form of the pipeline, the variation trend of the spiral case pressure and discharge also depends on the arrangement form of the pipeline. The larger the slope of the pipeline is, where the upstream water level is located, the larger the variation rate of the spiral case pressure and discharge is.
After the pump failure, because the guide vane fails to close, the unit discharge decreases sharply after the pump failure and then the positive discharge appears. The positive discharge increases to 0.0286 m 3 /s and then decreases to 0.0116 m 3 /s due to the increase of the unit speed. Finally, the unit discharge gradually decreases with the decrease of upstream water level. Similarly, the unit negative speed decreases rapidly and the positive speed appears due to the upstream water level. After the unit positive speed rapidly rises to 908.314 rpm, the unit positive speed slowly decreases with the decrease of upstream water level. The pressure of the spiral case decreases sharply to 849.48 cm after a power failure, then rises, and finally decreases continuously with the water level decreasing, while the draft tube inlet pressure increases sharply to 380.59 cm, then decreases, and finally remains stable.
In the experiment, the pump starts with the ball valve closed, which results in a high pressure in the spiral case after the pump is started and before the ball valve is opened. After the guide vane and ball valve are opened, the case pressure decreases quickly. This process cannot be simulated in the calculations, and hence the experimental results are quite different from the simulation results within seconds of startup.
Comparing the numerical results with the experimental data, the numerical simulation can accurately simulate the change process of each parameter. From Table 1 it can be The calculation results and experimental data show that the simulation results of the characteristic parameters and water head of each monitoring point are consistent with the change process of the experimental data. Comparing with the experimental data, the speed is different from the actual result. The main reason for these differences is that the characteristic curve cannot reflect the actual operation of the unit. In this condition, the guide vane fails to close and the intake gate closes after the unit rejection load, the unit reaches the runaway condition and finally enters the "S" area of the pump-turbine. In the "S" area, the characteristic curves are often quite different from the actual unit, so there are certain errors between the calculated results and the experimental results.
When the gas-liquid coupling interface is in a horizontal pipe, there are certain errors between the numerical simulation and the experiment. In the experiment, the values of spiral case pressure, unit discharge and unit speed are slowly reduced, but the corresponding values remain unchanged in the simulation. This is because the interface tracing method is based on the assumption that the water surface is always perpendicular to the axis of the pipe and that the water surface is at the center of the pipeline [29]. When the gas-liquid coupling interface is located in a horizontal pipe, the numerical simulation will assume that the water level remains constant while the water level remains at the center of the pipe. However, in the experiment, the water level always keeps horizontal in the horizontal pipeline, and slowly reduces from the top to the bottom of the pipeline. In addition, during the closing phase of the intake gate, the fluctuation of the characteristic parameters in the simulation has a longer duration than in the experiments. The main reason is that the mathematical model of the gate closing cannot represent the actual closing condition of the gate. During the closing phase of the gate, the actual upstream water level is slowly decreasing, but in the simulation, the water level will drop sharply at the end of the gate closing. As a result, the duration of fluctuation increases.
To further verify the accuracy of the solution method proposed in this paper, the condition B was simulated. The comparison of experimental data and simulation results of characteristic parameters under condition B are shown in Figure 11 and the comparison of extreme values of characteristic parameters is shown in Table 2 (note that speed and discharge under pump condition are negative).  Figure 12 show that: The adoption of different closing intervals does not affect the overall change trend of the characteristic parameters, and still depends on the layout of upstream pipelines. Closing the intake gate 5 s after the unit load rejection can reduce the time for the unit to quit runaway condition and help the unit return to the safe state faster. After the interval exceeds 10 s or more, with the extension of the interval, the longer the unit stays in the runaway state of high speed, the greater the harm to the unit. Therefore, closing the intake gate at an appropriate interval time after the load rejection of the unit can shorten the time for the unit to quit runaway condition, so as to ensure the safety of the unit.

Analysis of Closing Time of Intake Gate
The intake gate is closed after 50 s after load rejection of the unit, and the closing time of the intake gate is given as 1.5 s, 5 s, 10 s and 15 s respectively. It is calculated by using the program in this paper and the calculation results are shown in Figure13.  The experimental data show that after the guide vane is opened, the discharge of the unit rises rapidly and then decreases with the increase of water level in the upstream pipeline after reaching the maximum value. The pressure of the spiral case decreases rapidly within seconds after starting and then rises with the rise of water level in the upstream pipeline. The inlet pressure of the draft tube decreases due to the increase of discharge and then rises slightly due to the decrease of discharge. Because the variation trend of upstream water level depends on the arrangement form of the pipeline, the variation trend of the spiral case pressure and discharge also depends on the arrangement form of the pipeline. The larger the slope of the pipeline is, where the upstream water level is located, the larger the variation rate of the spiral case pressure and discharge is.
After the pump failure, because the guide vane fails to close, the unit discharge decreases sharply after the pump failure and then the positive discharge appears. The positive discharge increases to 0.0286 m 3 /s and then decreases to 0.0116 m 3 /s due to the increase of the unit speed. Finally, the unit discharge gradually decreases with the decrease of upstream water level. Similarly, the unit negative speed decreases rapidly and the positive speed appears due to the upstream water level. After the unit positive speed rapidly rises to 908.314 rpm, the unit positive speed slowly decreases with the decrease of upstream water level. The pressure of the spiral case decreases sharply to 849.48 cm after a power failure, then rises, and finally decreases continuously with the water level decreasing, while the draft tube inlet pressure increases sharply to 380.59 cm, then decreases, and finally remains stable.
In the experiment, the pump starts with the ball valve closed, which results in a high pressure in the spiral case after the pump is started and before the ball valve is opened. After the guide vane and ball valve are opened, the case pressure decreases quickly. This process cannot be simulated in the calculations, and hence the experimental results are quite different from the simulation results within seconds of startup.
Comparing the numerical results with the experimental data, the numerical simulation can accurately simulate the change process of each parameter. From Table 1 it can be seen that the extreme values of each parameter are basically in agreement with the measured values, the relative error of the minimum pressure of the spiral case is 8.49%, the relative error of the maximum inlet pressure of the draft tube is 4.88%, the relative error of the maximum positive discharge of the unit is 15.38%, the relative error of the minimum negative discharge of the unit is 3.24%, and the relative error of the maximum reverse speed of the unit after pump failure is 4.2%. The main reason for the difference of discharge after reversal is that the characteristic curve used for the calculations cannot represent the running state of the pump-turbine in the experiment.
In conclusion, the solution method and program for the transient process of a gasliquid two-phase flow in a pumped storage power station based on a gas-liquid interface coupling proposed in this paper can accurately simulate the position of water surface in the pipeline and the change trend of characteristic parameters during the transient process. The comparison of three simulation methods for the gas-liquid two-phase transient process is shown in Table 3. Table 3. The comparison of three simulation methods for the gas-liquid two-phase transient process.

Comparison Items
The In case of unit load rejection and guide vane closing failure, closing the intake gate can prevent the unit from entering runaway operations, so as to ensure the operational safety of the unit. Therefore, it is very important to select the appropriate gate closing interval time and gate closing time. Based on the calculation model of the experimental platform, this paper calculates and analyzes the effect of gate closing interval time and gate closing time on the characteristic parameters.

Analysis of Interval Time of Intake Gate Closing
The intake gate shall be closed after 0 s, 5 s, 10 s, 20 s and 40 s, respectively, after the unit load rejection. The calculation is carried out by using the program in this paper, and the calculation results are shown in Figure 12.  Figure 13 shows that the change of gate closing time cannot change the overall change trend of the unit quitting the runaway process. The state of water flow passing through the intake gate becomes gentle due to the extension of gate closing time, and the amplitude of the parameter oscillation decreases. Meanwhile, the time of the gas-liquid coupling surface from the upper flat section to the inclined pipe section increases with the increase of the closing time of the intake gate, and the time of the unit at high speed also increases, which is not conducive to the safe operation of the unit. Therefore, when the intake gate is closed, it should be closed in a short time.

Scenario2: Pump Failure after Low-Head Startup
When the startup mode of the first unit is under the pumping condition, the water storage requirements of the upper reservoir can be reduced and the construction period of the project can be shortened. However, since a low-head startup under pump condition is a complex gas-liquid two-phase process, the previous one-dimensional numerical calculation cannot accurately simulate it, and the time required for three-dimensional numerical calculation is too long to systematically analyze the complex operating conditions after the pump startup at the low-head. Therefore, the unit stability is not clear under a pump failure condition after a low-head startup.
To study the stability of the unit under pump failure after a low-head startup condition and the influence mechanism of the gas-liquid coupling interface on the hydraulic transient process, based on the mathematical model, solution method and calculation program proposed in this paper, four different pump failures after low-head startup conditions were selected for analysis. The detailed information of the four conditions is shown  Figure 12 show that: The adoption of different closing intervals does not affect the overall change trend of the characteristic parameters, and still depends on the layout of upstream pipelines. Closing the intake gate 5 s after the unit load rejection can reduce the time for the unit to quit runaway condition and help the unit return to the safe state faster. After the interval exceeds 10 s or more, with the extension of the interval, the longer the unit stays in the runaway state of high speed, the greater the harm to the unit. Therefore, closing the intake gate at an appropriate interval time after the load rejection of the unit can shorten the time for the unit to quit runaway condition, so as to ensure the safety of the unit.

Analysis of Closing Time of Intake Gate
The intake gate is closed after 50 s after load rejection of the unit, and the closing time of the intake gate is given as 1.5 s, 5 s, 10 s and 15 s respectively. It is calculated by using the program in this paper and the calculation results are shown in Figure 13.  Table 4 and the results are shown in Figure 14 and   Figure 13 shows that the change of gate closing time cannot change the overall change trend of the unit quitting the runaway process. The state of water flow passing through the intake gate becomes gentle due to the extension of gate closing time, and the amplitude of the parameter oscillation decreases. Meanwhile, the time of the gas-liquid coupling surface from the upper flat section to the inclined pipe section increases with the increase of the closing time of the intake gate, and the time of the unit at high speed also increases, which is not conducive to the safe operation of the unit. Therefore, when the intake gate is closed, it should be closed in a short time.

Scenario2: Pump Failure after Low-Head Startup
When the startup mode of the first unit is under the pumping condition, the water storage requirements of the upper reservoir can be reduced and the construction period of the project can be shortened. However, since a low-head startup under pump condition is a complex gas-liquid two-phase process, the previous one-dimensional numerical calculation cannot accurately simulate it, and the time required for three-dimensional numerical calculation is too long to systematically analyze the complex operating conditions after the pump startup at the low-head. Therefore, the unit stability is not clear under a pump failure condition after a low-head startup.
To study the stability of the unit under pump failure after a low-head startup condition and the influence mechanism of the gas-liquid coupling interface on the hydraulic transient process, based on the mathematical model, solution method and calculation program proposed in this paper, four different pump failures after low-head startup conditions were selected for analysis. The detailed information of the four conditions is shown in Table 4 and the results are shown in Figure 14 and Table 5 (note that the speed and discharge under pump condition are negative).      a positive speed under the action of water flow, but will slowly decrease. P1 and P2 conditions generate a positive discharge under the action of high upstream water level, but the maximum positive discharge under P1 is less than that under P2 and the maximum discharge occurs earlier than under P2. Due to the continuous closure of the guide vane, the discharge under the P1 condition gradually decreases to 0 m 3 /s after the maximum positive discharge. The spiral case pressure has a downward fluctuation trend after the pump failure. The valley of spiral case pressure under P1 condition is 1028.76 cm. After the pressure rises, the spiral case pressure under the P1 condition does not decrease but oscillates continuously. The inlet pressure of the draft tube has an upward fluctuation trend. The maximum pressure under the P1 condition is 347.77 cm, which is less than 351.96 cm under the P2 condition.
Comparing the P2, P3 and P4 conditions, with the decrease of the pump failure water level, the pump failure time of the unit is advanced. After a power failure, the unit speed and discharge will be positive. When the water level of pump failure is 13.49 m, 12.63 m, and 11.78 m, the maximum reverse discharge is 0.028 m 3 /s, 0.0264 m 3 /s, 0.0251 m 3 /s, and the maximum reverse speed is 1010.71 rpm, 921.52 rpm, and 872.79 rpm, respectively. With the decrease of the water level, the maximum reversal flow discharge and reversal speed also decrease. Therefore, it can be considered that the maximum reverse discharge and reverse speed are positively correlated with the water level of the power failure.
Comparing P2, P5, and P6 conditions, the spiral case pressure and draft tube inlet pressure under the P5 and P6 conditions generate large numerical oscillations during the opening of the guide vane in the simulation. After the guide vane has fully opened, the variation trend of characteristic parameters under each working condition is consistent. When the initial water level of pump startup is 4.0 m, 5.2 m, and 7.3 m, the maximum pumping discharge is 0.0554 m 3 /s, 0.052 m 3 /s, and 0.0494 m 3 /s, respectively. With the increase of the startup water level, the maximum pumping discharge decreases. Therefore, it can be considered that the maximum pumping discharge is negative correlated with the startup water level. Due to the consistency of the pump failure water level, the change trend of characteristic parameters after the pump failure is also consistent, and the extreme values are basically the same.
In general, after the pump startup, the change trend of the spiral case pressure and the discharge depend on the pipeline layout above the startup water level. After the pump failure, if the guide vanes are closed, the characteristic parameters can be maintained within the safety limits. If the guide vane fails to close, the speed and discharge of the unit will rapidly reverse under the action of high upstream water level to generate positive discharge and speed, and then will slowly decrease with the decrease of upstream water level. The value of maximum positive discharge and the positive speed depend on the pump failure water level.

Conclusions
Aiming at the limitations of traditional methods, which cannot accurately simulate gas-liquid coupling transient process of pumped storage system, a novel mathematical model was proposed based on Preissman's implicit difference scheme and the method of characteristics. Thereafter, the solving mechanism of the transient process of gas-liquid movement was developed on the gas-liquid interface tracking method. The proposed model was applied in two typical scenarios of the gas-liquid transient process in a hydropower system, i.e., the unit quitting runaway operations by the intake gate quick closure and a pump failure after the startup at low water head, and the field experiments were carried out accordingly. The major conclusions can be summarized as follows: (1) The novel mathematical model could be used to accurately simulate the position of the gas-liquid coupling surface and the characteristic parameters in a gas-liquid two-phase transient process. The average relative error was about 7.2%. (2) The influence mechanism of the gas-liquid coupling interface on the hydraulic transient process depended on the pipeline layout. The change speed of the spiral case pressure, unit discharge and unit speed were positively correlated with the slope of the pipeline where the gas-liquid coupling surface was located. (3) Under the quitting of runaway operations with a quick closure of the intake gate condition, the gate closing interval time and gate closing time did not affect the change trend of parameters, but the appropriate gate closing interval time could reduce the time for quitting the runaway condition. Under the pump failure after low-head startup condition, the maximum pumping discharge was negatively correlated with startup water level. The maximum reversal discharge and the maximum reversal speed after pump failure were positively correlated with the pump failure water level. Meanwhile, the state of guide vanes after a pump failure had a great influence on the reverse discharge and the reverse speed.
The proposed model could simulate the transient process of a pumped storage system coupling gas-liquid interface from a one-dimension perspective, it has potential application value in the safe operation of hydropower stations, the transient process of water diversion projects and in urban pipe network operation.