Next Article in Journal
Effect of the Presence of Virus-like Particles on Bacterial Growth in Sunlit Surface and Dark Deep Ocean Environments in the Southern East China Sea
Previous Article in Journal
An Assessment of Hydroacoustic and Electric Fishing Data to Evaluate Long Term Spatial and Temporal Fish Population Change in the River Thames, UK
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
2
Power China Zhongnan Engineering Corporation Limited, Changsha 410014, China
*
Author to whom correspondence should be addressed.
Water 2021, 13(20), 2933; https://doi.org/10.3390/w13202933
Submission received: 23 September 2021 / Revised: 16 October 2021 / Accepted: 17 October 2021 / Published: 19 October 2021
(This article belongs to the Section Hydrogeology)

Abstract

:
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 verification. 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 verified in a field 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.

1. 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 simulation method is no longer applicable, so it is necessary to correct the simulation of the gas–liquid two-phase transient process.
The research of pipeline transient flow mainly focuses on liquid transient flow, gas transient flow and gas–liquid two-phase transient flow. (1) The current calculation methods for a liquid transient flow include the method of characteristics (MOC), the equivalent circuit method, the finite-difference (FD) and the finite-volume (FV) method. The MOC proposed by Wylie [1] and Chaudhry [2] is currently the most widely used method to solve the water hammer equation. It is mostly applied to a fixed grid in computer programming, but the nodes of the grid are not always on the characteristic lines, so linear or nonlinear interpolation is required on the space or time grid. Larock [3] applied the linear interpolation in the MOC method, while Holly and Preissmann [4] utilized the nonlinear interpolation. Zhao et al. [5,6] simulated the transient process of a pumped storage power station by the equivalent circuit method and reconciled the contradiction between simulation efficiency and accuracy through the novel pump-turbine model and space-time discretization. Chaudhry and Hussaini [7] solved the water hammer equation by MacCormak and Gabutti’s explicit FD methods, and they found that these second-order schemes required fewer computing nodes and less time for the same calculation accuracy compared with the MOC method. Ghidaoui [8], Guinot [9], Hwang [10], Zhao [11] and Sabbagh-Yazdi [12] formulated Godunov-type first-order or second-order explicit FV methods to simulate the water hammer, and all of the results show that the second-order Godunov-type method is more accurate. (2) The governing equations of liquid and gas transients in pipelines are similar, so the numerical methods mentioned above are also applicable to the simulation of gas transportation in pipelines. Wylie et al. [13] and Osiadacz [14] used the MOC method to simulate gas flows. Ke and Ti [15] used the equivalent circuit method to simulate isothermal gas flow in pipeline network, which has certain computational advantages compared with the traditional method. Thorley [16] used the FD method to simulate air flow and found that the FD method had high accuracy, but a higher accuracy can usually be achieved at the expense of increased computational labor. Kessal [17] reduced the computational time of sophisticated finite difference schemes by coupling two different types of finite difference schemes. Ibraheem [18] combined the TVD method with the upwind scheme to solve the problem of dynamic pressure wave and found that the whole process was stable, robust, and accurate. Zhou [19] simulated the transport of a fast transient in a short pipe through the TVD/Roe scheme and found that this scheme could capture and maintain the integrity of the wave fronts even after a long time. (3) Compared with other gas–liquid two-phase flow problems, the gas–liquid two-phase transient process in a pumped storage power station is a water filling or draining process in the pipeline. Wang et al. [20] established a mathematical model to solve the gas–liquid two-phase transient flow in undulating pipes by using the VOF model and the standard k-ε turbulence model, and carried out three-dimensional numerical simulation of the gas–liquid two-phase transient flow in the filling process of V-shaped undulating pipes with stagnant airbag. Zhang [21,22] solved the water filling process of a pipeline containing trapped air mass by using the MOC method and the method of lines. Liu [23] calculated the water filling process with gas at the end of the pipe by using the rigid plug elastic water column model. Wang [24,25] adopted the shock-fitting method to build comprehensive pipeline filling models and analyzed the transient flow of pipeline water filling in various layout forms.
Among the previously mentioned methods, the MOC method has the advantages of easy programming and high calculation speed, and it is more suitable for the calculation of unsteady flow in long invariable-area pipelines because of the limitations of a strict adherence to the time–step–space relationship. Though the FD and FV methods have higher calculation accuracy, they usually require more computing resources. In addition, none of them can solve the problem of gas–liquid two-phase transient flow. Though the three-dimensional numerical simulation can accurately simulate the gas–liquid two-phase problem, the calculation is time-consuming and it is difficult to set up the model of the 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.

2. Mathematical Model

2.1. 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]:
Continuity   equation : B 2 A M x + p t = 0
Momentum   equation : p x + α 2 A M t + p g sin β B 2 + f B 2 M | M | 2 D A 2 p = 0
where the relevant parameters are defined as follows: B is the wave velocity (m/s); A is the sectional area of the pipeline (m2); 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 (m2/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:
  • M = ρ v A , where v is flow velocity (m/s).
  • The wave velocity can be determined by using the state equation B = p / ρ = λ R T ; where ρ is the gas density (kg/m3), λ is the compressibility coefficient, R is the gas constant (m2/(s2·K)) and T is the absolute temperature (K).
  • 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.
  • The default values of the parameters are: B = 340 m/s, α = 1, f = 0.015, ρ = 1.205 kg/m3, p0 = 101,325 Pa.
The schematic of Preissmann’s four-point finite-difference implicit method with the temporal weighting coefficient θ is shown in Figure 1. Equation (3) represents the discretization of Preissmann’s scheme corresponding to Figure 1, where f denotes any function.
{ f ( x , t ) = θ 2 ( f i + 1 n + 1 + f i n + 1 ) + 1 θ 2 ( f i + 1 n + f i n ) f x = θ f i + 1 n + 1 f i n + 1 Δ x + ( 1 θ ) f i + 1 n f i n Δ x f t = f i + 1 n + 1 f i + 1 n + f i n + ! f i n 2 Δ t
Equations (1) and (2) are discretized according to Preissmann’s four-point finite-difference 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):
A 1 G , i · Δ p i + 1 + B 1 G , i · Δ M i + 1 = C 1 G , i · Δ p i + D 1 G , i Δ M i + F 1 G , i
A 2 G , i · Δ p i + 1 + B 2 G , i · Δ M i + 1 = C 2 G , i · Δ p i + D 2 G , i Δ M i + F 2 G , i
where:
  • A 1 G , i = 1 2 Δ t ,   B 1 G , i = B 2 θ A Δ x ,   C 1 G , i = 1 2 Δ t ,   D 1 G , 1 = B 2 θ A Δ x , F 1 G , i = 1 Δ x ( M i + 1 n M i n )
  • A 2 G , i = θ Δ x + g θ sin β 2 B 2 f B 2 ( M i + 1 n + M i n ) | M i + 1 n + M i n | 4 D A 2 ( p i + 1 n + p i n ) 2 , B 2 G , i = α 2 2 A Δ t + f B 2 | M i + 1 n + M i n | 2 D A 2 ( p i + 1 n + p i n )
  • C 2 G , i = θ Δ x g θ sin β 2 B 2 + f B 2 ( M i + 1 n + M i n ) | M i + 1 n + M i n | 4 D A 2 ( p i + 1 n + p i n ) 2
  • D 2 G , i = θ Δ x g θ sin β 2 B 2 + f B 2 ( M i + 1 n + M i n ) | M i + 1 n + M i n | 4 D A 2 ( p i + 1 n + p i n ) 2
  • F G , i = 1 Δ x ( p i + 1 n p i n ) g sin β 2 B 2 ( p i + 1 n + p i n ) f B 2 ( M i + 1 n + M i n ) | M i + 1 n + M i n | 4 D A 2 ( p i + 1 n + p i n )
For gas pipe sections solved by the implicit method, the following relationships exist among the internal calculated sections [28]:
Δ M i = E E G , i · Δ p i + F F F , i
Δ p i = L G , i · Δ p i + 1 + R G , i · Δ M i + 1 + N G , i
where:
  • L G , i = A 1 G , i C 1 G , i + D 1 G , i · E E G , i
  • R G , i = B 1 G , i C 1 G , i + D 1 G , i · E E G , i
  • N G , i = D 1 G , i · F F G , i + F 1 G , i C 1 G , i + D 1 G , i · E E G , i
  • E E G , i = A 1 G , i ( C 2 G , i + D 2 G , 2 · E E G , i 1 ) + A 2 G , i ( C 1 G , i + D 1 G , i · E E G , i 1 ) B 1 G , i ( C 2 G , i + D 2 G , i · E E G , i 1 ) B 2 G , i ( C 1 G , i + D 1 G , i · E E G , i 1 )
  • F F G , i = ( D 1 G , i · F F G , i 1 + F 1 G , i ) ( C 2 G , i + D 2 G , i · E E G , i 1 ) B 1 G , i ( C 2 G , i + D 2 G , i · E E G , i 1 ) B 2 G , i ( C 1 G , i + D 1 G , i · E E G , i 1 ) ( D 2 G , i · F F G , i 1 + F 2 G , i ) ( C 1 G , i + D 1 G , i · E E G , i 1 ) B 1 G , i ( C 2 G , i + D 2 G , i · E E G , i 1 ) B 2 G , i ( C 1 G , i + D 1 G , i · E E G , i 1 )
In the equation, E E i and F F i are the transfer coefficients and only depend on the E E i 1 and F F i 1 of the front section, while L i , R i , and N i depend on the transfer coefficients E E i and F F i . Through the boundary conditions of the first section surface, the E E 1 and F F 1 of the first node can be obtained, so the E E i , F F 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.

2.2. Basic Equation and Solution Method for Unsteady Liquid Flow

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 follows [1]:
Continuity   equation : Q A H x + H t + a 2 g A Q x + a 2 Q g A A x sin α · Q A = 0
Momentum   equation : g A 2 H x + Q Q x + A Q t + f Q | Q | 2 D = 0
where x is the distance (m); t is the time (s); Q is the discharge (m3/s); H is the pressure head (m); a is the wave speed (m/s); A is the cross-sectional area (m2); β is the pipe slope; g is the gravitational acceleration (m2/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):
A 1 L , i · Δ H i + 1 + B 1 L , i · Δ Q i + 1 =   C 1 L , i · Δ H i + D 1 L , i · Δ Q i + F 1 L , i
A 2 L , i · Δ H i + 1 + B 2 L , i · Δ Q i + 1 = C 2 L , i · Δ H i + D 2 L , i · Δ Q i + F 2 L , i
where:
  • A 1 L , i = 1 + θ Δ t Δ x ( Q i + 1 A i + 1 + Q i A i ) , B 1 L , i = a 2 g A i + 1 2 θ Δ t Δ x
  • C 1 L , i = 1 + θ Δ t Δ x ( Q i + 1 A i + 1 + Q i A i ) , D 1 L , i = a 2 g A i 2 θ Δ t Δ x
  • F 1 L , i = Δ t Δ x ( Q i + 1 A i + 1 + Q i A i ) ( H i + 1 H i ) a 2 g 2 Δ t Δ x ( Q i + 1 A i + 1 Q i A i ) a 2 g Δ t Δ x ( Q i + 1 A i + 1 2 + Q i A i 2 ) ( A i + 1 A i ) + Δ t sin α ( Q i + 1 A i + 1 + Q i A i )
  • A 2 L , i = g θ Δ t Δ x , B 2 L , i = 1 2 A i + 1 [ 1 + θ Δ t Δ x ( Q i + 1 A i + 1 + Q i A i ) ]
  • C 2 L , i = g θ Δ t Δ x , D 2 L , i = 1 2 A i [ 1 θ Δ t Δ x ( Q i + 1 A i + 1 + Q i A i ) ]
  • F 2 L , i = Δ t 2 Δ x ( Q i + 1 2 A i + 1 2 Q i 2 A i 2 ) g Δ t Δ x ( H i + 1 H i ) Δ t 4 ( f Q i | Q i | A i 2 D i + f Q i + 1 | Q i + 1 | A i + 1 2 D i + 1 )
For the liquid pipe section solved by the implicit scheme method, the internal calculation sections also have the following relationship:
Δ Q i = E E L , i · Δ H i + F F L , i
Δ H i = L L , i · Δ H i + 1 + R L , i · Δ Q i + 1 + N L , i
where:
  • L L , i = A 1 L , i C 1 L , i + D 1 L , i · E E L , i , R L , i = B 1 L , i C 1 L , i + D 1 L , i · E E L , i , N L , i = D 1 L , i · F F L , i + F 1 L , i C 1 L , i + D 1 L , i · E E L , i
  • E E L , i = A 1 L , i ( C 2 L , i + D 2 L , 2 · E E L , i 1 ) + A 2 L , i ( C 1 L , i + D 1 L , i · E E L , i 1 ) B 1 L , i ( C 2 L , i + D 2 L , i · E E L , i 1 ) B 2 L , i ( C 1 L , i + D 1 L , i · E E L , i 1 )
  • F F L , i = ( D 1 L , i · F F L , i 1 + F 1 L , i ) ( C 2 L , i + D 2 L , i · E E L , i 1 ) B 1 L , i ( C 2 L , i + D 2 L , i · E E L , i 1 ) B 2 L , i ( C 1 L , i + D 1 L , i · E E L , i 1 ) ( D 2 L , i · F F L , i 1 + F 2 L , i ) ( C 1 L , i + D 1 L , i · E E L , i 1 ) B 1 L , i ( C 2 L , i + D 2 L , i · E E L , i 1 ) B 2 L , i ( C 1 L , i + D 1 L , i · E E L , i 1 )
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 differential equations on two characteristic lines, as shown in Figure 2 and as expressed in Equations (14) and (15).
C + : Q P = C P B p · H P
C : Q P = C M + B M · H P
where:
  • B P = 1 ( C 0 C 2 ) / A + C 0 C 1 , B M = 1 ( C 0 + C 2 ) / A + C 0 C 3
  • C P = B P ( Q R ( C 0 + C 2 ) A + H R ) , C M = B M ( Q S ( C 0 + C 2 ) A + H S )
  • C 0 = a g , C 1 = f Δ t | Q R | 2 D A 2 , C 2 = 1 2 Δ t sin β , C 3 = f Δ t | Q S | 2 D A 2
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.

2.3. 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 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 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 Zn section, totaling Nk segments, where the length of Nk−1 segments is Δx and the length of a segment (i.e., the Z section to the Zk+1 section) is x G , Z = Δ x ( Z Z k / Z k + 1 Z k ) . 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 A G L , Z = A L G , Z , then Q G L , Z = Q L G , Z and the gas mass flow can be written as Equation (16):
M G L , Z = ρ G Q L G , Z
The relationship between the head of the liquid pressure measuring pipe and the gas pressure can be expressed as follows:
H L G , Z = Z Z D + ( p G L , Z p a ) ρ L g
The x G , Z 0 at the last time step has the following relationship with the x G , Z at the current time step:
x G L , Z x G L , Z 0 = ( M G L , Z + M G L , Z 0 ) ρ G Δ t 2 A G L , Z
Bringing x G , Z = Δ x ( Z Z k / Z k + 1 Z k ) and x G , Z 0 = Δ x ( Z 0 Z k / Z k + 1 Z k ) into Equation (18) and the variation of interface position of the gas–liquid coupling can be expressed as follows:
Δ Z = ( 2 M G L , Z 0 + Δ M G L , Z ) ρ G W 0 sin θ
where sin θ = Z k + 1 Z k Δ x and W 0 = Δ t 2 A G L , Z . 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:
Δ Q L G , Z = E E L G , Z · Δ H L G , Z + F F L G , Z
Linearizing Equations (16) and (17) can get Equations (21) and (22):
Δ M G L , Z = ρ G Δ Q L G , Z
Δ H L G , Z = Δ Z + 1 ρ L g Δ p G L , Z
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):
Δ M G L , Z ρ G = E E L G , Z [ ( 2 M G L , Z 0 + Δ M G L , Z ) ρ G W 0 sin θ + 1 ρ L g Δ p G L , Z ] + F F L G , Z
Simplify Equation (23) to get Equation (24):
Δ M G L , Z = E E L G , Z · Δ p G L , Z + F F G L , Z
where:
E E G L , Z = ρ G E E L G , Z ρ L g ( 1 E E L G , Z · W 0 sin θ ) ,   F F G L , Z = 2 E E L G , Z · W 0 sin γ · M G L , Z 0 + ρ G F F L G , Z ( 1 E E L G , Z · W 0 sin θ )

3. Boundary Conditions

3.1. 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:
C + : Q P = C P B M · H P
C : Q P = C M + B M · H S
The electric generator first-order equation can be expressed as follows:
n = n 0 + 0.1875 ( M t + M t 0 ) Δ t / G D 2
The unit parameter equation of the turbine can be expressed as follows:
Q P = Q 1 D 1 2 H P H s
n 1 = n D 1 / H P H s
n 1 = n D 1 / H P H s
The characteristic curve equation of the turbine can be expressed as follows:
Q 1 = f ( α , n 1 )
M 1 = g ( α , n 1 )
The guide vane specified equation can be expressed as follows:
α = α ( t )
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 (m3/s); M t denotes the turbine output,(N·m); M g denotes the generating torque (N·m); G D 2 denotes the rotary inertia (kg·m2); D 1 denotes the diameter of the turbine (m); n 1 denotes the unit speed (rpm); Q 1 denotes the unit discharge (m3/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.

3.2. 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.
  • 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:
Z U ( H P + Z D ) = 1 2 g ( φ A ) 2 Q P 2
The C equation behind the gate inlet can be expressed as follows:
Q P = C M + B M · H P
where Z U is the upper reservoir water level (m); Z D is the down reservoir water level, A is the area of pipeline (m2); φ 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 (m3/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 (1) the equation of the orifice outflow, (2) the C equation, (3) the equation of gas state, (4) the relationship equation between the piezometer head H f and the absolute pressure p in the tube, and (5) the equation of the venthole mass gas flow. These equations together include a total of five unknown values, namely, Q P 1 , H P , Q P 2 , V , and m ˙ .
The equation of the orifice outflow can be expressed as follows:
Z U ( H P + Z D ) = 1 2 g ( φ A ) 2 Q P 1 2
where Hp is the pressure head behind the gate (m); Q P 1 is the discharge behind the gate (m3/s).
The C equation behind the gate inlet can be expressed as follows:
Q P 2 = C M + B M · H P
where Q P 2 is the discharge behind airbag (m3/s).
The equation of the gas state can be expressed as follows:
p V = m M R T
where p is the gas pressure (m); V is the gas volume (m3); m is the gas mass (kg); M is the molar mass of gas (kg/mol); R is the gas constant (m2/(s2·K)); T is the absolute temperature (K).
The relationship equation between the piezometer head Hf and the absolute pressure p in the tube can be expressed as follows:
H P = p γ + ( Z K Z D ) H a
where Ha is the atmospheric pressure head (m); γ is the liquid specific weight (kg/m3); Zk 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:
m ˙ = C i n A i n 7 p a ρ a ( ( p p a ) 1.4286 ( p p a ) 1.714 ) , 0.528 p a < p < p a
where C i n is the venthole discharge coefficient; A i n is the area of venthole (m2); ρ a is the mass density of atmospheric air ρ a = p a M / R T a .
  • Critical flow in:
    m ˙ = C i n A i n 0.686 R T a p a , p 0.528 p a
  • Subsonic air flow out:
    m ˙ = C o u t A o u t 7 R T ( ( p p a ) 1.4286 ( p p a ) 1.714 ) , p a < p < p a 0.528
  • Critical flow out:
    m ˙ = C o u t A o u t 0.686 R T a p a , p a 0.528 p
3.
The inlet gate is fully closed.
When the intake gate is fully closed, Q P 1 = 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:
Δ p G , 0 = p G , 0 p G , 0 0 = 1 E E G , 0 ( M G , 0 M G , 0 0 F F G , 0 )

3.3. 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:
Δ Q L G , N = E E L G , N · Δ H L G , N + F F L G , N
where E E L G , N = C Q M , 0 , F F L G , N = Q C M , 0 0 + Q C M , 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:
E E G , N = E E G L , 0 , F F G , N = F F G L , 0
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., pG,0 = p0 and Δ p G , 0 = 0 By using Equation (7) the unknown quantity Δ M G of this boundary node can be obtained as follows:
Δ M G , 0 = F F G , 0

4. A Novel Analytical Model and Verification

4.1. 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 gas–liquid 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.

4.2. 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 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.
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 Figure 9 and Figure 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 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).
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 m3/s and then decreases to 0.0116 m3/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 gas–liquid 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.

5. Influence Mechanism of Gas–liquid Interface Coupling on Hydraulic Transient Process

5.1. Scenario1: Quitting Runaway Operations with Quick Closure of the Intake Gate

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.

5.1.1. 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 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.

5.1.2. 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.
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.

5.2. 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).
Comparing the P1 condition with the P2 condition, because their pumping processes are completely consistent, the change trend of characteristic parameters and the maximum pumping discharge are the same. After the pump failure, the unit speed decreases sharply. Since the guide vane is closed under the P1 working condition, the unit will not generate 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 m3/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 m3/s, 0.0264 m3/s, 0.0251 m3/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 m3/s, 0.052 m3/s, and 0.0494 m3/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.

6. 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.

Author Contributions

Conceptualization, C.L. and J.Y. (Jiandong Yang); methodology, T.P.; software, J.Y. (Jiebin Yang); writing—review and editing, Z.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was jointly supported by National Natural Science Foundation of China (No. 51879200, 51839008).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wylie, E.B.; Streeter, V.L. Fluid Transients in Systems; Prentice-Hall: Englewood Cliffs, NJ, USA, 1993. [Google Scholar]
  2. Chaudhry, H.M. Applied Hydraulic Transients; Van Nostrand Reinhold: New York, NY, USA, 1987. [Google Scholar]
  3. Larock, B.E.; Jeppson, R.W.; Watters, G.Z. Hydraulics of Pipeline Systems; Crc Press: Boca Raton, FL, USA, 1999. [Google Scholar]
  4. Holly, F.M.; Preissmann, A. Accurate Calculation of Transport in Two Dimensions. J. Hydraulic. Div. Proc. ASCE 1977, 103, 1259–1277. [Google Scholar] [CrossRef]
  5. Zhao, Z.G.; Yang, J.D.; Yang, W.J.; Hu, J.H.; Chen, M. A coordinated optimization framework for flexible operation of pumped storage hydropower system: Nonlinear modeling, strategy optimization and decision making. Energy Convers. Manag. 2019, 194, 75–93. [Google Scholar] [CrossRef]
  6. Zhao, Z.G.; Yang, J.D.; Chung, C.Y.; Yang, W.J.; He, X.H.; Chen, M. Performance enhancement of pumped storage units for system frequency support based on a novel small signal model. Energy 2021, 234, 121207. [Google Scholar] [CrossRef]
  7. Chaudhry, M.H.; Hussaini, M.Y. Second-Order Accurate Explicit Finite-Difference Schemes for Waterhammer Analysis. J. Fluids Eng. 1985, 107, 523–529. [Google Scholar] [CrossRef]
  8. Ghidaoui, M.S.; Zhao, M.; McInnis, D.A.; Axworthy, D.H. A Review of Water Hammer Theory and Practice. Appl. Mech. Rev. 2005, 58, 49–76. [Google Scholar] [CrossRef]
  9. Guinot, V. Riemann solvers for water hammer simulations using Godunov Method. Int. J. Numer. Methods Eng. 2000, 49, 851–870. [Google Scholar] [CrossRef]
  10. Hwang, Y.-H.; Chung, N.-M. A fast Godunov method for the water-hammer problem. Int. J. Numer. Methods Fluids 2002, 40, 799–819. [Google Scholar] [CrossRef]
  11. Zhao, M.; Ghidaoui, M.S. Godunov-type solutions for water hammer flows. J. Hydraul. Eng. 2004, 130, 341–348. [Google Scholar] [CrossRef]
  12. Sabbagh-Yazdi, S.R.; Mastorakis, N.E. Water Hammer Modeling by Godunov type Finite Volume Method. Int. J. Math. Comput. Simul. 2007, 1, 350–355. [Google Scholar]
  13. Streeter, V.L.; Wylie, E.B. Fluid Transients; McGraw-Hill: New York, NY, USA, 1978. [Google Scholar]
  14. Osiadacz, A. Simulation and Analysis of Gas Networks; Gulf Publishing Company: Houston, TX, USA, 1987. [Google Scholar]
  15. Ke, S.L.; Ti, H.C. Transient analysis of isothermal gas flow in pipeline network. Chem. Eng. J. 2000, 76, 169–177. [Google Scholar] [CrossRef]
  16. Thorley, A.R.D.; Tiley, C.H. Unsteady and transient flow of compressible fluids in pipelines—A review of theoretical and some experimental studies. Int. J. Heat Fluid Flow 1987, 8, 3–15. [Google Scholar] [CrossRef]
  17. Kessal, M. Simplified Numerical Simulation of Transients in Gas Networks. Chem. Eng. Res. Des. 2000, 78, 925–931. [Google Scholar] [CrossRef] [Green Version]
  18. Ibraheem, S.O.; Adewumi, M.A. On Total Variation Diminishing Schemes for Pressure Transients. J. Energy Resour. Technol. 1999, 121, 122–130. [Google Scholar] [CrossRef]
  19. Zhou, J.; Adewumi, M.A. Simulation of transients in natural gas pipelines using hybrid TVD schemes. Int. J. Numer. Methods Fluids 2000, 32, 407–437. [Google Scholar] [CrossRef]
  20. Wang, Q.; Zhang, X.; Zhang, K.; Yu, B. Gas-Liquid Two-Phase Flow Numerical Simulation of Undulating Pipelines. J. Beijing Inst. Petrochem. Technol. 2018, 26, 26–30. [Google Scholar]
  21. Zhang, Y.; Vairavamoorthy, K. Application of Method-of-Lines to Charging-Up Process in Pipelines with Entrapped Air. Tsinghua Sci. Technol. 2006, 11, 324–331. [Google Scholar] [CrossRef]
  22. Zhang, Y.; Vairavamoorthy, K. Transient Flow in Rapidly Filling Air-Entrapped Pipelines with Moving Boundaries. Tsinghua Sci. Technol. 2006, 11, 313–323. [Google Scholar] [CrossRef]
  23. Liu, D.; Zhou, L.; Karney, B.; Zhang, Q.; Ou, C. Rigid-plug elastic-water model for transient pipe flow with entrapped air pocket. J. Hydraul. Res. 2011, 49, 799–803. [Google Scholar] [CrossRef]
  24. Wang, L.; Wang, F.; Karney, B.; Malekpour, A. Numerical investigation of rapid filling in bypass pipelines. J. Hydraul. Res. 2017, 55, 647–656. [Google Scholar] [CrossRef]
  25. Wang, L.; Wang, F.-J.; Huang, J. Numerical investigation of filling transients in small-scale pipelines with submerged outlet. J. Hydrodyn. 2019, 31, 145–151. [Google Scholar] [CrossRef]
  26. Guo, W.; Yang, J.; Chen, J.; Teng, Y. Effect Mechanism of Penstock on Stability and Regulation Quality of Turbine Regulating System. Math. Probl. Eng. 2014, 2014, 349086. [Google Scholar] [CrossRef] [Green Version]
  27. Lyn, D.A.; Goodwin, P. Stability of a General Preissmann Scheme. J. Hydraul. Eng.-ASCE 1987, 113, 16–28. [Google Scholar] [CrossRef]
  28. Wang, C.; Yang, J.D. Water Hammer Simulation Using Explicit-Implicit Coupling Methods. J. Hydraul. Eng. 2015, 141, 04014086. [Google Scholar] [CrossRef]
  29. Vasconcelos, J.; Wright, S. Applications and Limitations of Single-Phase Models to the Description of the Rapid Filling Pipe Problem. J. Water Manag. Modeling 2005, 13, 377–408. [Google Scholar] [CrossRef]
Figure 1. Schematic diagram of Preissmann’s implicit method [27].
Figure 1. Schematic diagram of Preissmann’s implicit method [27].
Water 13 02933 g001
Figure 2. Characteristic lines in xt plane of MOC.
Figure 2. Characteristic lines in xt plane of MOC.
Water 13 02933 g002
Figure 3. Schematic diagram of pipeline segment.
Figure 3. Schematic diagram of pipeline segment.
Water 13 02933 g003
Figure 4. Schematic diagram of gas–liquid pipeline.
Figure 4. Schematic diagram of gas–liquid pipeline.
Water 13 02933 g004
Figure 5. Schematic diagram of intake gate closure. (a) The intake gate fully opened; (b) the intake gate closing; (c) the inlet gate fully closed.
Figure 5. Schematic diagram of intake gate closure. (a) The intake gate fully opened; (b) the intake gate closing; (c) the inlet gate fully closed.
Water 13 02933 g005
Figure 6. The complete simulation process of a transient process of a gas–liquid two-phase flow in a pumped storage power station.
Figure 6. The complete simulation process of a transient process of a gas–liquid two-phase flow in a pumped storage power station.
Water 13 02933 g006
Figure 7. Field diagram of the experimental platform.
Figure 7. Field diagram of the experimental platform.
Water 13 02933 g007
Figure 8. Schematic diagram of the experimental platform.
Figure 8. Schematic diagram of the experimental platform.
Water 13 02933 g008
Figure 9. The comparison of water head at each monitoring point under condition A. (a) Monitoring point 1; (b) monitoring point 2; (c) monitoring point 3; (d) monitoring point 4; (e) monitoring point 5; (f) monitoring point 6.
Figure 9. The comparison of water head at each monitoring point under condition A. (a) Monitoring point 1; (b) monitoring point 2; (c) monitoring point 3; (d) monitoring point 4; (e) monitoring point 5; (f) monitoring point 6.
Water 13 02933 g009
Figure 10. The comparison of the characteristic parameters under condition A. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Figure 10. The comparison of the characteristic parameters under condition A. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Water 13 02933 g010aWater 13 02933 g010b
Figure 11. The comparison of the characteristic parameters under condition B. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Figure 11. The comparison of the characteristic parameters under condition B. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Water 13 02933 g011
Figure 12. The comparison of results under different gate closing interval times. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Figure 12. The comparison of results under different gate closing interval times. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Water 13 02933 g012
Figure 13. The comparison of results under different gate closing interval times. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Figure 13. The comparison of results under different gate closing interval times. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Water 13 02933 g013
Figure 14. The result of transient process of pump failure after low-head startup. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Figure 14. The result of transient process of pump failure after low-head startup. (a) The pressure of spiral case; (b) the inlet pressure of draft tube; (c) the discharge of unit; (d) the rotational speed of unit.
Water 13 02933 g014aWater 13 02933 g014b
Table 1. The detail of pipeline system segments.
Table 1. The detail of pipeline system segments.
Pipe No.Length
(m)
Equivalent Diameter
(m)
Starting Point Elevation
(m)
End Point Elevation
(m)
Pipe No.Length
(m)
Equivalent Diameter
(m)
Starting Point Elevation
(m)
End Point Elevation
(m)
L12.2180.71513.54213.542L100.2760.151.1811.181
L22.1320.37713.54213.46L111.20960.151.1811.181
L37.97860.3513.3087.425L121.7960.2161.1120.694
L46.04730.357.4257.304L132.7770.30.6940.694
L58.08210.357.3041.349L141.4210.3680.6940.751
L61.4570.31.3491.232L1517.6240.4260.7512.296
L71.0070.21.2321.181L160.8010.4292.2962.296
L80.51130.1781.1811.181L172.2180.7352.2962.296
L90.3890.151.1811.181
Table 2. The comparison of extreme values of characteristic parameters.
Table 2. The comparison of extreme values of characteristic parameters.
Comparative ItemThe Valley of Spiral Case Pressure (cm)Maximum Draft Tube Inlet Pressure (cm)Unit Discharge Extremum (m3/s)Maximum Rotational Speed (rpm)
PositiveNegative
Experimental data849.48380.5970.0286−0.0556908.314
Numerical results777.33362.0250.0242−0.0538870.175
Relative error (%)8.494.8815.383.244.20
Table 3. The comparison of three simulation methods for the gas–liquid two-phase transient process.
Table 3. The comparison of three simulation methods for the gas–liquid two-phase transient process.
Comparison ItemsThe Novel ModelThe MOC MethodThree-Dimensional Numerical Simulation
Gas–liquid two-phase transient process?YesNoYes
Simulation timeMinutes/Days
Simulation accuracyHigh accuracy and small simulation step/High accuracy
Internal flow patternNoNoYes
Table 4. Detailed information of pump failure after low-head startup.
Table 4. Detailed information of pump failure after low-head startup.
Condition NumberStartup Water Level (m)Downstream Water Level (m)Guide Vane Opening Time (s)Pump Failure Water Level (s)Guide Vane State after Pump Failure
P14.03.85513.49Closes normally
P24.03.85513.49Fails to close
P34.03.85512.63Fails to close
P44.03.85511.78Fails to close
P55.23.85513.49Fails to close
P67.33.85513.49Fails to close
Table 5. Extreme values of characteristic parameters.
Table 5. Extreme values of characteristic parameters.
Condition NumberThe Valley of Spiral Case Pressure (cm)/Time (s)Maximum Draft Tube Inlet Pressure (cm)/Time (s)Unit Discharge ExtremumMaximum Rotational Speed (rpm) /Time (s)
Positive (m3/s)
/Time (s)
Negative (m3/s)
/Time (s)
P11028.76/45.31347.77/45.340.0146/46.01−0.0554/4.34−16.63/100
P21005.37/46.24351.96/45.260.028/48.28−0.0554/4.341010.71/52.67
P3961.06/39.31364.35/39.310.0264/42.47−0.0554/4.34921.52/47.05
P4901.06/36.52344.42/36.520.0251/39.84−0.0554/4.34872.79/44.65
P51005.22/43.00362.43/43.010.028/46.03−0.052/4.611009.71/50.41
P61004.76/37.47348.42/37.50.028/40.52−0.0494/5.611010.15/44.90
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, C.; Peng, T.; Yang, J.; Zhao, Z.; Yang, J. Transient Process of Pumped Storage System Coupling Gas–Liquid Interface: Novel Mathematical Model and Experimental Verification. Water 2021, 13, 2933. https://doi.org/10.3390/w13202933

AMA Style

Liu C, Peng T, Yang J, Zhao Z, Yang J. Transient Process of Pumped Storage System Coupling Gas–Liquid Interface: Novel Mathematical Model and Experimental Verification. Water. 2021; 13(20):2933. https://doi.org/10.3390/w13202933

Chicago/Turabian Style

Liu, Chengpeng, Tao Peng, Jiebin Yang, Zhigao Zhao, and Jiandong Yang. 2021. "Transient Process of Pumped Storage System Coupling Gas–Liquid Interface: Novel Mathematical Model and Experimental Verification" Water 13, no. 20: 2933. https://doi.org/10.3390/w13202933

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop