Modeling of the Free-Surface-Pressurized Flow of a Hydropower System with a Flat Ceiling Tail Tunnel

Abstract

For a water diversion hydropower system with a flat ceiling tail tunnel with high elevation, during transient states with relatively low tail water levels, free-surface-pressurized flow inevitably appears and its transient characteristics have obvious effects on the system’s operating stability. Using Newton–Raphson linearization in the characteristic implicit format for modeling of the free-surface-pressurized flow in the tail tunnel, the mathematical models for necessary boundary conditions were derived and linear algebraic equations with a band coefficient matrix were grouped for further transient simulation. Then, a unified mathematical model was established for hydraulic transient analysis of the hydropower system with free-surface-pressurized flow. Combined with experimental research and numerical simulation, the wave speed for the free-surface-pressurized flow was experimentally analyzed for further correctness in the unified model, and by comparative analysis the hydraulic characteristics of the free-surface-pressurized flow in the flat ceiling tail tunnel were investigated. It was found that the derived mathematical model can basically represent water behaviors in the water-surface-pressurized flow, the wave speed for the mixed water-surface-pressurized flow can be set to approximately 50m/s, and with this correctness the numerical results are in good agreement with the experimental results. Therefore, the obtained mathematical model combined with an experimental wave speed or a reference wave speed of 50 m/s for the free-surface-pressurized flow is preferable during the design stage of the hydropower system.

1. Introduction

For the development of a water diversion hydropower station with an underground powerhouse, a practical type of tail tunnel is introduced into the layout design of a large-scale tail system. In this tail system, in order to reduce the tail tunnel’s excavation and rock mass stability, the diversion tunnel that originally served during construction is redesigned and reconstructed as the downstream part of the tail tunnel of the hydropower system (see Figure 1).

As shown in Figure 1, the tail tunnel of this hydropower system mainly includes a lower tail tunnel directly connecting with the tail branches and a higher flat ceiling tail tunnel, which was originally part of a diversion tunnel during construction, combined with a mid-connecting tunnel with a large reverse slope. In addition, in order to mitigate the pressure oscillations along the tail tunnel, three air vents are appropriately set at the crown of the tail tunnel. During normal operation, two typical flow patterns exist along the tail tunnel system, namely:

• Pressurized flow. If the tail water level is higher than the top elevation of the tail tunnel outlet, the flat ceiling tail tunnel is always in a typical pressurized state.
• Free-surface-pressurized flow. If the tail water level is obviously lower than the top elevation of the tail tunnel outlet, the flat ceiling tail tunnel is always in free channel flow under both steady and transient states, and a unique interface between pressurized flow and free surface flow is located along the connecting tunnel. Particularly as the tail water level is slightly lower than the top elevation of tail tunnel outlet, the interface between the pressurized flow and the free surface flow is located along the flat ceiling tail tunnel under steady states, and possible mixed free-surface-pressurized flow will inevitably happen under transient states, with one or more air masses existing along the crown of the tail tunnel in some cases.

For the pressurized flow, the method of characteristics for pressurized systems is commonly used for detailed hydraulic transients and reveals the inherent transient characteristics of the tail system [1]. For free-surface-pressurized flow, the flow pattern in the tail tunnel is widely varied, including free surface flow, pressurized flow, and air–water two-phase flow. Particularly when the tail water level is slightly lower than the top elevation of the tail tunnel outlet, there are transient flow problems with entrapped air and trapped air mass in the mixed free-surface-pressurized flow regime. It is difficult to accurately simulate and clearly understand the transient characteristics of the tail system under these complex flow patterns or to clarify their effect on the hydropower system’s operating stability.

Therefore, more research studies have been devoted to the profound analysis and exact modeling of the free-surface-pressurized flow based on experimental research and numerical simulation, leading to many achievements. Free-surface-pressurized flow is free surface flow in which the conduit is pressurized during the transient state, and often occurs in sewers and in the conduits of hydroelectric power plants or pumped storage projects [2]. Based on an overall review and analysis, it was pointed out that the free-surface-pressurized flow in storm water systems is difficult to capture in modeling, and rapid pipe filling or emptying of water mains and sewer systems is accompanied by transitions between free surface and pressurized flow regimes and subatmospheric unsteady flow [3,4]. Focusing on the free-surface-pressurized flow, based on different experimental setups and further data analysis, the obtained experimental results were not only used for verification of numerical models but also revealed some obvious phenomena involved in air–water interactions; for example, the air near the pipe crown may pressurize and lead to consequent intense pressure oscillations [5,6,7,8]. Furthermore, considering the obvious air–water interaction in the free-surface-pressurized flow, the effect of air cavity intrusion into horizontal and inclined pipes on flow behavior was also experimentally investigated, enabling different degrees of ventilation by various orifices [9,10,11]. Considering the typical flat ceiling tail system in Figure 1 in particular, for the possible free-surface-pressurized flow, the types of air–water interactions observed in the combined diversion tunnel are divided into single air pocket motion, multiple air pocket motion, interfacial instability, and negligible interactions [12].

With a clear understanding of the free-surface-pressurized flow provided by experimental research, recent research studies have emphasized that it is more important to perform the exact numerical simulation for the free-surface-pressurized flow in different water systems, and in most cases the given prototype systems are too difficult or uneconomical to be modeled in the lab. For numerical simulation of the free surface flow, a family of well-balanced, semi-implicit numerical schemes was proposed and further proved to be reliable for solving engineering problems [13]. In simple hydraulic systems containing a channel, considering the partial free surface and partial pressurized flow, fundamental numerical simulation models were derived and presented by means of the control volume method or other methods [14,15]. More practically, based on the traditional Preissmann model [2], some numerical simulation models were constructed for drainage networks, including a model based on an introduced virtual slot on the crown of the pipe to treat a separated gas–liquid flow [16,17]; a model capable of simulating transient flows in closed conduits, ranging from free surface flows to mixed flows and fully pressurized flows [18]; the storm water management model (SWMM), which has wide practical applications [19]; a model based on the first-order Roe’s scheme within the framework of finite volume methods [20]; and a discrete model with a four-point linear implicit format [21]. From the viewpoint of engineering applications, considering the possible transient mixed free-surface-pressurized flow in tailrace tunnels of large hydropower stations, a new method was proposed—namely, the implicit method of characteristics based on an implicit finite difference scheme [22]—and was used for detailed characteristic analysis of the free-surface-pressurized flow for all kinds of tunnel network topologies, including large fluctuation computation, hydraulic disturbance analysis, and small disturbance analysis, mostly for changing top-altitude tail tunnels [23,24,25]. Particularly, based on the assumption of a rigid incompressible water column and a compressible air bubble, a numerical model was derived to simulate pressure fluctuation, void fraction, air–water flow rate, and water velocity in a closed conduit [26]. Focusing on the special free-surface-pressurized flow with a clear and regularly moving interface region in the changing top-altitude tail tunnel, a three-dimensional computational fluid dynamics research code with the volume of fluid (VOF) model was applied [27]. In addition, for some other specified water systems with possible free-surface-pressurized flow, the numerical methods mainly include the computational fluid dynamics (CFD) models [28], the high-precision discontinuous Galerkin finite element method [29], a model used to simulate vertical water level fluctuations with coupled liquid and gas phases [30], and a mathematical model based on dam break theory and bore flow theory [31]. For the irrigation distribution system, a fully implicit time scheme for free-surface-pressurized water flow was developed and validated using a typical test system and prototype experiment [32]. Most recently, an explicit smoothed particle hydrodynamics model for incompressible fluid was presented to simulate flow in conduits during transitions between free surface and pressurized flow [33]. A semi-implicit numerical model with a linear solver was proposed for mixed free surface and pressurized flow in hydraulic systems [34]. A novel 1D–2D coupled model was presented for accurate simulation of transient flow hydrodynamics in urban drainage systems, which was further confirmed by some test cases using comparative analysis with other existing models [35]. Using comparative analysis of all the aforementioned numerical models, most of them are presented for different draining networks and specified waterworks, while for modeling of the free-surface-pressurized flow in the hydropower system in Figure 1, the characteristic implicit method [22,23,24,25] can be appropriately used for the hydraulic transient analysis, with further mathematical modeling of necessary algorithm solutions and boundary conditions.

This paper aims to numerically model the free-surface-pressurized flow in a hydropower system with a flat ceiling tail tunnel, along with further hydraulic characteristics analysis. Considering that other state-of-the-art models are basically used for typical water systems and considering the difficulty in introducing these into the numerical simulation for an entire hydropower system with various and complex boundaries, the characteristic implicit format, which is an improved slot model, is preferred. Therefore, based on the Newton–Raphson linearization of the basic equations for the characteristic implicit method, the corresponding mathematical models for necessary boundary conditions are built according to the characteristic implicit format, and then the mathematical model for the hydraulic characteristics in the connecting tunnel and flat ceiling tail tunnel, which is presented by linear algebraic equations with a band coefficient matrix, is constructed with appropriate algorithm methods. Next, combined with the method of characteristics for pressurized pipelines, the mathematical model of a downstream surge tank and surge unit’s motion equation, and their detailed hydraulic characteristics, a unified mathematical model is established for hydraulic transient analysis of the given hydropower systems. For the tail tunnel system with possible free-surface-pressurized flow, experimental research is preferred to investigate its complex hydraulic characteristics. Based on the experimental setup in the lab and further data analysis, particularly on the wave speed for the free-surface-pressurized flow, the corresponding wave speed for the computation of Preissmann slot in the unified model is corrected together with sensitivity analysis, and then the detailed hydraulic characteristics of the free-surface-pressurized flow in the flat ceiling tail tunnel are further revealed, accompanied by comparative analysis with experimental data.

2. Material and Methods

2.1. Governing Equations

For both the free surface flow and pressurized flow during transient states, the general governing equations include a continuity equation and momentum equation:

$$B\frac{\partial h}{\partial t}+\frac{\partial Q}{\partial x}=0$$

(1)

$$\frac{\partial Q}{\partial t}+\frac{2Q}{A}\frac{\partial Q}{\partial x}+gA\frac{\partial h}{\partial x}-\frac{Q^2}{A^2}\frac{\partial A}{\partial x}=gA(i-J_f)$$

(2)

where Q is the transient volumetric discharge at different sections, h is the water depth, A is the cross-sectional area, B is the width of the water surface, i is the bottom slope of the open channel, and J_f is the hydraulic gradient calculated from the Manning formula.

2.2. Characteristic Implicit Format

Considering that the commonly used Preissmann slot model [2] may result in divergent computation, in order to find a reasonable convergent format, two coefficients are introduced:

$$l^{\pm }=-\frac{Q}{A}\pm \sqrt{gA/B},c^{\pm }=\frac{Q}{A}\pm \sqrt{gA/B}$$

(3)

After multiplying Equation (1) by $l^{\pm }$ and then combining Equation (2) with the introduction of $c^{\pm }$, two characteristic equations are obtained.

$$B{l}^{+}\frac{\partial h}{\partial t}+gA\frac{\partial h}{\partial x}+\frac{\partial Q}{\partial t}+c^{+}\frac{\partial Q}{\partial x}=gA(i-J_f)+\frac{Q^2}{A^2}\frac{\partial A}{\partial x}$$

(4)

$$B{l}^{-}\frac{\partial h}{\partial t}+gA\frac{\partial h}{\partial x}+\frac{\partial Q}{\partial t}+c^{-}\frac{\partial Q}{\partial x}=gA(i-J_f)+\frac{Q^2}{A^2}\frac{\partial A}{\partial x}$$

(5)

For Equations (4) and (5) at a given section j and reference time step n, the forward difference is applied for ∂/∂t items and different difference methods for ∂/∂x items. The gravity item, friction item, and ∂A/∂x item are calculated on time step (n + 1), and the others are on time step n. Finally, the characteristic implicit format [22] is presented with the two basic difference equations listed below:

$$a_1{h}_{j-1}^{}+b_1{Q}_{j-1}^{}+c_1{h}_j^{}+d{}_{{}_1}Q{}_j^{}=u_1$$

(6)

$$a_2{h}_j+b_2{Q}_j^{}+c_2{h}_{j+1}^{}+d{}_{2}Q{}_{j+1}^{}=u_2$$

(7)

where subscripts j − 1, j, j + 1 are the numbers of three neighboring sections; coefficients a_i, b_i, c_i, d_i (i = 1, 2), and right items u_i (i = 1, 2) are determined by the known parameters and variables of the relevant sections.

For frequent and exact numerical computation, after the incremental expressions of water depth and volumetric discharge of each section are introduced by linearization with the Newton–Raphson method, Equations (6) and (7) are rewritten as below:

$$a_{1j}\Delta h_{j-1}^{}+b_{1j}\Delta Q_{j-1}^{}+c_{1j}\Delta h_j^{}+d_{1j}\Delta Q_j^{}=e_{1j}$$

(8)

$$a_{2j}\Delta h_j^{}+b_{2j}\Delta Q_j^{}+c_{2j}\Delta h_{j+1}^{}+d{}_{2j}\Delta Q_{j+1}^{}=e_{2j}$$

(9)

where Δh and ΔQ are the water depth increment and volumetric discharge increment, respectively; and the second subscript j in coefficients a_ij, b_ij c_ij, d_ij, and right term e_ij (i = 1, 2) is the calculated section.

Based on the characteristic implicit format, Equations (8) and (9)—the total linear equations—can be obtained to describe the transient characteristics for a simple tail tunnel with m total sections, which are grouped as a matrix form as below:

AX = E

(10)

where A is the 2 m × 2 m coefficient matrix presented as a band matrix; X is a 2 m column vector, X = (Δh₁, ΔQ₁, Δh₂, ΔQ₂, …, Δh_j, ΔQ_j, …, Δh_m, ΔQ_m)^T; and E is also a 2 m column vector, E = (e₁₁, e₂₁, e₁₂, e₂₂, …, e_1j, e_2j, …, e_1m, e_2m)^T.

2.3. Boundary Conditions

For a typical water diversion system containing an open channel, free flow tunnel, changing top-altitude tail tunnel, or tail tunnel with possible free-surface-pressurized flow, the commonly used boundary conditions comprise an inlet section, serial sections with or without overflow weirs, bifurcations, and gate shafts. Considering the hydropower system in Figure 1, the free-surface-pressurized flow will inevitably appear along the connecting tunnel, which will have a relatively large reverse slope, together with its downstream flat ceiling tail tunnel, while its upstream tail tunnel will always be in a pressurized state. This is because the corresponding boundaries relating to the free-surface-pressurized flow mainly involve the inlet section of the connecting tunnel, which will have a relatively large reverse slope, the gate shaft at the mid-section of flat ceiling tail tunnel, the outlet of the flat ceiling tail tunnel, and the series sections. Here, based on the built characteristic implicit format (Equations (8) and (9)), the mathematical models of the above boundaries are derived.

2.3.1. Inlet Section

At the inlet section of the tail tunnel system with possible free-surface-pressurized flow (shown in Figure 2), the upstream side is a lower tail tunnel that is always in a pressurized state. The last section of the upstream tunnel is k and the first section of the downstream tunnel is 1, and then based on the method of characteristics, the C⁺ equation for section k is:

$$C^{+}:\begin{array}{cc}& H_{Pk}=C_P-B_P{Q}_{Pk}\end{array}$$

(11)

where C_P and B_P are constants calculated from the piezometric head and volumetric discharge of the neighboring section at time t − △t; H_Pk and Q_Pk are the piezometric head and volumetric discharge of section k at time t.

According to the volumetric discharge and head balance conditions at the inlet section, $Q_{Pk}=Q_1$and $H_{Pk}=h_1+{\nabla }_1$, in which H₁ and Q₁ are the piezometric head and volumetric discharge of section 1 at time t, ∇₁ is the bottom elevation at the inlet section, and then Equation (11) is reorganized into the following form:

$$F_1=h_1+B_P{Q}_1+{\nabla }_1-C_P=0$$

(12)

Equation (12) is linearized by using the Newton–Raphson method with the introduction of the water depth increment Δh and volumetric discharge increment ΔQ. Then, the required first-row elements of band matrix A and the corresponding right-hand item can be obtained, which describe the hydraulic characteristics of the inlet section.

$$\left[\begin{array}{cccc}1& B_P& & \\ a_{21}& b_{21}& c_{21}& d_{21}\end{array}\right]\left\{\begin{array}{c}\Delta h_1\\ \Delta Q_1\\ \Delta h_2\\ \Delta Q_2\end{array}\right\}=\left\{\begin{array}{c}-F_1\\ e_{21}\end{array}\right\}$$

(13)

2.3.2. Gate Shaft in Mid-Section of the Flat Ceiling Tail Tunnel

For a relatively long flat ceiling tail tunnel, a tail gate shaft is often set at its mid-section. Under free-surface-pressurized flow, the water level at this section is sometimes below the crown of the bottom tunnel and is just one part of a free surface channel, while sometimes the water level varies in the gate well and works as a surge tank. In Figure 3, at the gate shaft section, the neighboring sections are j–j and (j+1)–(j+1), and the controlling equations include:

$$H_j=H_{j+1}=H_P$$

(14)

$$H_P=Z_{PS}+R_S\left|Q_{PS0}\right|Q_{PS}^{}$$

(15)

$$Q_j=Q_{j+1}=Q_{PS}$$

(16)

$$Z_{PS}=Z_{PS0}+\frac{\left(Q_{PS}+Q_{PS0}\right)\Delta t}{2A_S}$$

(17)

where H_P is piezometric head at the bottom tunnel of gate shaft; H_j and H_j+1 are piezometric heads at sections j–j and (j+1)–(j+1); Q_j and Q_j+1 are volumetric discharges at sections j–j and (j+1)–(j+1); Z_PS and Z_PS0 are the instantaneous and initial water levels in the gate shaft, respectively; Q_PS and Q_PS0 are the instantaneous and initial volumetric discharge flowing into gate shaft, respectively; R_S is the head loss coefficient for water flowing into or out of the gate shaft; A_S is the effective area of gate shaft.

If we substitute Equation (17) into Equation (15) and then let$S_T=\frac{\Delta t}{2A_S}+R_S\left|Q_{PS0}\right|$, $C_T=\frac{\Delta t}{2A_S}{Q}_{PS0}+H_{PS0}$, the linear characteristic equation for the gate shaft is derived:

$$H_P=S_T{Q}_{PS}+C_T$$

(18)

We substitute Equation (16) into Equation (18), and then yield

$$H_P=S_T(Q_j-Q_{j+1})+C_T$$

(19)

According to head balance condition, $H_P=h_j+{\nabla }_j,H_P=h_{j+1}+{\nabla }_{j+1}$, in which ${\nabla }_j$ and ${\nabla }_{j+1}$are bottom elevations at sections j–j and (j+1)–(j+1), combined with ${\nabla }_j$=${\nabla }_{j+1}$, we substitute these conditions into Equation (19) and yield:

$$F_2=h_j-S_T{Q}_j+S_T{Q}_{j+1}+{\nabla }_j-C_T=0$$

(20)

$$F_3=-S_T{Q}_j+h_{j+1}+S_T{Q}_{j+1}+{\nabla }_{j+1}-C_T=0$$

(21)

Equations (20) and (21) are linearized by using the Newton–Raphson method. With the introduction of the water depth increment and volumetric discharge increment at each section, the required elements in 2j and 2j+1 rows of band matrix A and the corresponding right-hand items can also be obtained, which describe the hydraulic characteristics of the gate shaft section.

$$\left[\begin{array}{cccccccc}{a}_{1j}& b_{1j}& c_{1j}& d_{1j}& & & & \\ & & 1& -S_T& 0& S_T& & \\ & & 0& -S_T& 1& S_T& & \\ & & & & a_{2,j+1}& b_{2,j+1}& c_{2,j+1}& d_{2,j+1}\end{array}\right]\left\{\begin{array}{c}\Delta h_j\\ \Delta Q_j\\ \Delta h_{j+1}\\ \Delta Q_{j+1}\end{array}\right\}=\left\{\begin{array}{c}{e}_{1j}\\ -F_2\\ -F_3\\ e_{2,j+1}\end{array}\right\}$$

(22)

Equation (22) is used for the gate shaft section when the water level varies in the gate shaft, while if the water level at this section is below the crown of the bottom tunnel, the gate shaft section is simplified into a series section with the volumetric discharge and head balance conditions $Q_j=Q_{j+1}$and $h_j=h_{j+1}$, respectively. Then, Equation (22) is transferred into

$$\left[\begin{array}{cccccccc}{a}_{1j}& b_{1j}& c_{1j}& d_{1j}& & & & \\ & & 1& 0& -1& 0& & \\ & & 0& 1& 0& -1& & \\ & & & & a_{2,j+1}& b_{2,j+1}& c_{2,j+1}& d_{2,j+1}\end{array}\right]\left\{\begin{array}{c}\Delta h_j\\ \Delta Q_j\\ \Delta h_{j+1}\\ \Delta Q_{j+1}\end{array}\right\}=\left\{\begin{array}{c}{e}_{1j}\\ -h_j+h_{j+1}\\ -Q_j+Q_{j+1}\\ e_{2,j+1}\end{array}\right\}$$

(23)

2.3.3. Outlet of the Flat Ceiling Tail Tunnel

Figure 4 is the outlet of the flat ceiling tail tunnel with the tail water level, Z_w, and the outlet’s bottom elevation, ${\nabla }_m$. This boundary meets the following equation

$$F_4=h_m+{\nabla }_m-Z_w-\xi \frac{Q_m\left|Q_m\right|}{2g{A}_m{}^2}=0$$

(24)

where ξ is the minor head loss coefficient at the outlet, A_m is the sectional area of the outlet, and the subscript m represents the outlet section of the flat ceiling tail tunnel.

After linearizing Equation (24) with the Newton–Raphson method and introducing the incremental representation of corresponding parameters, the 2 m elements row in the band matrix A and the corresponding right-hand items can be modified according to Equation (25).

$$\left[\begin{array}{cccc}{a}_{1,m}& b_{1,m}& c_{1,m}& d_{1,m}\\ & & 1& -\xi \frac{\left|Q_m\right|}{g{A}_m{}^2}\end{array}\right]\left\{\begin{array}{c}\Delta h_{m-1}\\ \Delta Q_{m-1}\\ \Delta h_m\\ \Delta Q_m\end{array}\right\}=\left\{\begin{array}{c}{e}_{1,m}\\ -F_4\end{array}\right\}$$

(25)

2.4. Analysis Model of Transient Flow in System

Based on the characteristic implicit format (Equations (8) and (9)), as well as the corresponding initial and boundary conditions, a complete band matrix and further integrated mathematical model are derived to simulate the free-surface-pressurized flow in the flat ceiling tail tunnel. Combined with the method of characteristics for pressurized pipelines, the mathematical model of a downstream surge tank and surge unit’s motion equation, and their detailed hydraulic characteristics [1,2], a unified mathematical model is established for hydraulic transient analysis of the hydropower systems with possible free-surface-pressurized flow along the tail tunnel.

3. Results and Discussions

3.1. Experimental Research

3.1.1. Experiment Description

Before detailed numerical simulation and analysis based on the presented characteristic implicit format and corresponding boundary conditions, a complete experimental investigation is conducted for a hydropower system with a flat ceiling tail tunnel, as shown in Figure 1. This hydropower system consists of two parallel water diversion penstocks and tail branches, a downstream surge tank at the tail bifurcation with a diameter of 42.0 m, a D-shaped lower tail tunnel with a length of 626.80 m and B × H = 14.5 m × 18.0 m, a connecting tunnel with a length of 130.0 m and B × H = 17.5 m × 22.0 m, and a flat ceiling tail tunnel with a length of 485.83 m and B × H = 17.5 m × 22.0 m, in which B and H are the width and height of the D-shaped section. Figure 5 is the detailed longitudinal layout of the tail system with prototype sizes, including the length and elevation in m, and along the tail tunnel there are three air vents. The experimental setup is built with the model length scale λ_L = 60.0 and two pressure transducers are installed along the flat ceiling tail tunnel, as shown in Figure 5. One transducer is at the bottom of the combined section with the original diversion tunnel, and another is at the bottom of tail gate shaft section. The prototype length between these two transducers is 286.82 m.

As mentioned above, during normal operation or transient states, as the tail water level is slightly lower than the top elevation of the tail tunnel outlet, the interface between the pressurized flow and free surface flow is located along the flat ceiling tail tunnel under steady states, and possible free-surface-pressurized flow or even mixed flow will inevitably appear under transient states, with one or more air masses existing along the crown of the tail tunnel. Therefore, in the experimental research, the above worst cases with mixed free-surface-pressurized flow are underlined, and the details for these cases include:

C1: Reservoir’s water level is 825.0 m and tail water level is 595.83 m. Two units are in rated operation and have load rejection, with wickets normal closing at the same time.

C2: Reservoir’s water level is 825.0 m and tail water level is 595.83 m. One unit is in rated operation and another unit has load acceptance from idle state to rated operation, with wickets in normal open state.

The aforementioned two cases have a tail water level 595.83 m, which is slightly less than the top elevation of the tunnel’s outlet of 596.0 m. These cases are used for experimental research after the water levels are transferred to model test levels by referencing the elevation of the tail tunnel system, and in the following numerical simulation, they are also defined as the computation cases.

3.1.2. Wave Speed Analysis along the Flat Ceiling Tail Tunnel

Based on the experimental setup and according to the hydropower system in Figure 5, the aforementioned cases, case 1 and case 2, are simulated in the model hydropower system. With the measured data from the installed pressure transducers at two monitoring sections and the connected acquisition system, the corresponding dynamic data for the prototype can easily be obtained with reference to the model scale, including the pressure scale λ_H = λ_L = 60.0 and time scale λ_t= $\sqrt{{\lambda }_L}$ = 7.746. Figure 6 and Figure 7 give the time histories of the piezometric head at monitoring sections A and B under these two cases, together with the specified time as the first peak value at which pressure oscillations appear.

It can be revealed that because the tail water level of 595.83 m is slightly lower than the top elevation of the tunnel’s outlet of 596.0 m, during load rejection and load acceptance, the mixed free-surface-pressurized flow inevitably appears along the flat ceiling tail tunnel. The flow phenomena in case C1 were carefully observed, showing that at the time of load rejection, the water level along the flat ceiling tail tunnel horizontally falls down with crown void of the tunnel. In the 1st period with rising water level, the aerated flow enters into two air vents and the tail gate shaft, causing observable air–water interaction along the flat ceiling tail tunnel. In the 2nd period with rising water level, the observed phenomena are similar to the 1st period, with decayed pressure oscillation. Even in the 3rd period with rising water level, there is still a little aerated flow entering into the two air vents and the tail gate shaft. Until the 4th period with rising water level, no upwelling aerated flow is observed and the water level along the flat ceiling tail tunnel tends to be stable with decayed oscillation. The observed phenomena are in agreement with the dynamic curves at sections A and B in Figure 6.

Based on theoretical analysis and experimental observation, it is known that the pressure oscillation at the gate shaft (section B) always lags behind that at the combining section (section A), and the wave speed for the free-surface-pressurized flow, which is defined as a_f, can approximately reflect the pressure propagation characteristics from section A to section B, particularly in the process of the free-surface-pressurized flow, so according to the prototype distance between sections A and B (286.82 m), the approximate wave speed a_f for the free-surface-pressurized flow can be calculated and analyzed, which is listed in

Table 1. Wave speed analysis for the free-surface-pressurized flow.

Typical Cases	Reference Period	Peak Time (s)		Propagating Time (s)	Wave Speed, af (m/s)
		Section A	Section B
C1: Load rejection	1st pressure increase period	105.1	111.9	6.8	42.2
	2nd pressure increase period	266.9	272.8	5.9	48.6
C2: Load acceptance	1st pressure increase period	——	——	——	——
	2nd pressure increase period	153.6	158.7	5.1	56.2

.

It can be seen from Table 1 that for case C2, at the time of load acceptance, the free surface flow in the flat ceiling tail tunnel is pressurized steadily from upstream to downstream, without any obvious pressure oscillation in the 1st rising period. Therefore, the 2nd rising period, which involves entrapped air and resulting air–water interaction, is considered for wave speed analysis. For both load rejection and load acceptance cases, the calculated wave speed a_f in the free-surface-pressurized flow is from 40 to 60 m/s, and its approximate value can be set to a_f = 50 m/s, which can partly describe the propagation characteristics of the free-surface-pressurized flow.

3.2. Numerical Simulation by Using the Characteristic Implicit Method

3.2.1. Effect of Wave Speed on Transient Process

Based on the unified mathematical model established for hydraulic transient analysis of hydropower systems with possible free-surface-pressurized flow along the tail tunnel, further numerical computation and analysis can be conducted. Because the set air vents have a small sectional area and can be looked at as piezometric tubes, their effects are not considered in the numerical simulation. In the algebraic solution process with designed iterations, it is necessary to decide the wave speed a_f for the free-surface-pressurized flow in the tail tunnel, which is used to calculate the width of the Preissmann slot, $B_P=\raisebox{1ex}{$gA$}\!\left/ \!\raisebox{-1ex}{$a_f^2$}\right.$. Traditionally, through reference to pressurized flow, empirical data is often used for the wave speed a_f, which is relatively large and may result in an inevitable error in terms of simulation results, particularly for pressure oscillations. Hence, before detailed hydraulic transient computation and analysis, sensitivity analysis of the wave speed a_f or simulation of the free-surface-pressurized flow is carried out. Figure 8 gives the piezometric head at the combining section A for two computing cases, in which three different wave speeds a_f of 25 m/s, 50 m/s, and 100 m/s are introduced into the unified mathematical model, respectively.

It can be clearly found that for both load rejection and load acceptance, the calculated wave speed a_f for the free-surface-pressurized flow has an evident effect on the water behavior in the flat ceiling tail tunnel, basically on the pressure oscillation during the free-surface-pressurized flow; with the increase of wave speed a_f, the maximum pressure varies greatly, while the minimum pressure varies less. For computing case C1, because the maximum pressure is controlled by the maximum oscillation pressure in the 1st rising pressure period, the deviation of wave speed a_f will result in inaccurate evaluation of maximum pressure in the given sections, while for computing case C2, with the increase of wave speed a_f, the increasing maximum pressure in the 2nd rising pressure period tends to be greater than that in the 1st pressure rising period, leading to misinterpretation of the maximum pressure and its occurrence time. In summary, to accurately evaluate the hydraulic characteristics of the free-surface-pressurized flow along the flat ceiling tail tunnel, the premise is to take a relatively exact wave speed a_f for detailed numerical computation.

3.2.2. Comparative Analysis with Experimental Results

Based on the corresponding experimental results and the aforementioned analysis of the wave speed a_f of the free-surface-pressurized flow, the wave speed a_f is set to 50 m/s, and further numerical computation and analysis are implemented for two computing cases, C1 and C2. The obtained dynamic curves of the water level in the downstream surge tank and the piezometric head at combining section A are given in Figure 9 and Figure 10. In Figure 9 and Figure 10, the corresponding curves obtained from experimental research are also drawn for comparative analysis. The analysis of maximum and minimum values of water levels in the surge tank and piezometric head at combining section A for two computing cases C1 and C2 is shown in

Table 2. Comparative analysis of maximum (max.) and minimum (min.) values of two cases.

Typical Cases	Transient Variables		Experimental Data	Numerical Data	Error
C1: Load rejection	Water level in surge tank (m)	Max.	602.32 (134.9)	603.29 (126.8)	0.97
		Min.	585.33 (43.9)	585.49 (41.8)	0.16
	Piezometric head in section A (m)	Max.	605.10 (105.1)	605.18 (103.0)	0.08
		Min.	590.69 (62.0)	590.72 (57.9)	0.03
C2: Load acceptance	Water level in surge tank (m)	Max.	602.98 (50.3)	602.25 (46.8)	−0.73
		Min.	593.17 (126.0)	593.42 (120.6)	0.25
	Piezometric head in section A (m)	Max.	598.22(153.6)	598.09(160.3)	−0.13
		Min.	594.52(125.2)	594.51(128.6)	−0.01

, in which the data in parentheses is the occurrence time in s for the corresponding maximum or minimum value, and error = numerical data − experimental data.

As can be observed in Figure 9 and Figure 10 and Table 2, the numerical results for both water levels in the surge tank and piezometric heads in combining section A are in good agreement with the corresponding experimental results, including the maximum or minimum values and their corresponding occurrence times. Besides the approximately equal leading oscillation period, all the differences between numerical results and experimental results at any time are acceptable, and particularly the difference between the maximum and minimum data is less than 1.0 m. Most importantly, during the free-surface-pressurized flow, the detailed water behavior along the flat ceiling tail tunnel, characterized by the induced pressure oscillation with maximum amplitude, oscillation times, and decay rate, is highly similar. Therefore, with a relatively exact wave speed a_f for detailed numerical computation, the established unified model can accurately show typical water behavior for water-surface-pressurized flow along a flat ceiling tail tunnel, and can reveal the effects on the hydropower system’s dynamic characteristics.

4. Conclusions

In some water diversion hydropower systems with an underground powerhouse, the tail tunnel system includes a flat ceiling tail tunnel with high elevation, connected to the upstream lower tail tunnel via a connecting tunnel with a large reverse slope. During normal operation and in transient states, as the tail water level is lower than the top elevation of the tunnel’s outlet, free-surface-pressurized flow or even mixed flow will inevitably appear, and the complex transient characteristics have obvious effects on the hydropower system’s operating stability. Therefore, focusing on mathematical modeling and further numerical analysis of the free-surface-pressurized flow in this tail system, accompanied by experimental research and conducted studies, the obtained conclusions are as follows:

• Based on the characteristic implicit method for modeling of the free-surface-pressurized flow in the tail tunnel together with Newton–Raphson linearization, the linear algebraic equations with a band coefficient matrix are constructed, with the introduction of necessary boundary conditions for transient simulation of the free-surface-pressurized flow. Then, a unified mathematical model is established for hydraulic transient analysis of the given hydropower systems. This unified model can accurately reveal typical water behaviors in the water-surface-pressurized flow.
• With the built experimental setup in the lab and further data analysis, considering the dynamic curves of the piezometric head at two typical reference sections along the flat ceiling tail tunnel, the wave speed a_f for the free-surface-pressurized flow is experimentally analyzed, which is used for the correctness in the unified model. It is found that the wave speed a_f for the mixed water-surface-pressurized flow in the flat ceiling tail tunnel is close to 50 m/s.
• After the sensitivity analysis of wave speed a_f in the free-surface-pressurized flow, the detailed hydraulic characteristics of the free-surface-pressurized flow in the flat ceiling tail tunnel are further investigated and then confirmed by comparative analysis with experimental data. With appropriate correctness of wave speed a_f, the numerical results are in good agreement with the experimental results.

Therefore, for some hydropower stations, during the design stage of a tail tunnel system with possible free-surface-pressurized flow, experimental research is preferred to investigate the complex hydraulic characteristics, including wave speed evaluation. Then, by introduction of the obtained mathematical model combined with an experimental wave speed or reference wave speed 50 m/s for the free-surface-pressurized flow, detailed hydraulic transient computation and analysis under various cases for the entire hydropower system can be smoothly conducted.