Godunov-type solutions for free surface transient flow in pipeline incorporating unsteady friction

U

vector of flow variables

F

flux vector

S

vector containing source terms

A

cross-sectional area of water

Q

flow rate

average pressure of the water column on the cross-sectional area

ρ

density of liquid

g

the acceleration of gravity

S₀

slope of the pipe

S_f

pipeline friction

S_s

steady friction

S_u

unsteady friction

v

water velocity

n_m

Manning roughness coefficient

k

Brunone friction coefficient

C*

Vardy's shear decay coefficient

Re

Reynolds number

MINMOD

MINMOD slope limiter function

U_L, U_R

vector of flow variables that to left and right of interface i − 1/2, respectively

F_L, F_R

flux vector that to left and right of interface i − 1/2, respectively

S_L, S_R

wave velocities that to left and right of interface i − 1/2, respectively

A_L, A_R

cross-sectional area of water that to left and right of interface i − 1/2, respectively

θ

central angle of water surface

y

water depth

c

gravity wave celerity

average gravity wave celerity

g

gravitational acceleration

Cr

Courant number

N

number of grids

J_S

head losses per unit length due to steady friction

J_U

head losses per unit length due to unsteady friction

t

time

x

distance

Δx

space step

Δt

time step

Δτ

4vΔt/D²

τ_u

unsteady shear stress

ν

kinematic viscosity of the fluid

Free surface flow may frequently occur in various pipe systems, including sewer pipelines, power plants, and urban supply pipelines. Efficient and accurate free surface flow models are essential factors for real-time control of the pipeline systems to reproduce the rapidly changed hydraulics.

As pioneers of the field of free surface flow, Chow, Henderson, French, and Chaudhry have made great contributions to the research and propagation of free surface flow (Chow 1959; Henderson 1966; French 1985; Chaudhry 1987). The numerical models range from explicit numerical schemes to implicit numerical schemes (Yen 2001; Zhang et al. 2019; Retsinis et al. 2020). Most of the models that solve these equations use the implicit four-point Preissmann scheme when the formation of transients is secondary to issues such as conveyance capacity (Leon et al. 2006). Although the Preissmann scheme is unconditionally stable, the numerical dissipation caused by the time steps and weighting factors and damping limits its accuracy. The Preissmann scheme was also found not only to be unsuitable for transcritical flows but also has complex modeling, large storage requirement, and long simulation time when solving for branches or junctions in drainage systems (Sart et al. 2010; Wang et al. 2019). On the other hand, most free surface flow models are mainly designed to study the formation of hydraulic transients and adopt a scheme based on the characteristic method (MOC) (Leon et al. 2006). The MOC solution is easy to solve complicated boundary conditions. However, when the Courant number is less than one, spaced interpolations are inevitable, which often leads to numerical damping of the water depth at crest and the responding time delay of water surface line. Importantly, the MOC scheme cannot guarantee the conservation of water mass. Therefore, it is vital to develop an alternative solution scheme that can easily handle boundary conditions but provide robust and higher resolution.

The finite-volume Godunov-type scheme (GTS) is widely used because of its ability to conserve mass and to provide a high resolution of discontinuous points without spurious oscillations (Hirsch 1990). Alcrudo et al. (1992) proposed an upwind finite difference scheme based on flux difference splitting for solving the equations of open channel hydraulics. A limited second-order correction is added, making the resulting scheme robust and precise in calculating open channel flow. Fujihara & Borthwick (2000) presented an accurate Godunov-type numerical model for two-dimensional conservative hyperbolic shallow water equations. The results show that the model can accurately simulate a flow discontinuity between subcritical and supercritical conditions. Zhao & Ghidaoui (2004) proposed the second-order GTS to solve pressurized flow and pointed out that compared with MOC, the second-order GTS calculation is more accurate with the same Courant number, but the method is actually first-order accurate at the boundary. Leon et al. (2006) improved the GTS model to realize the second-order solution at the boundaries.

In general, the existing GTS model only considers the steady friction in the numerical simulations of the free surface transient pipe flow. However, in fast transient pipe flows, the steady friction approximation cannot accurately model pressure decay beyond the first pressure wave period (Daily et al. 1956; Wan & Huang 2018; Duan et al. 2020; Negharchi & Shafaghat 2021). Subsequently, various unsteady friction models (UFMs) have been developed to simulate pressure dissipation (Zielke 1968; Trikha 1975; Brunone et al. 1991; Suzuki et al. 1991; Vardy & Brown 1995; Bergant et al. 2001; Ghidaoui & Mansour 2002; Vardy & Brown 2003; Vardy & Brown 2004; Vítkovský et al. 2006; Urbanowicz 2018; Zhou et al. 2021a, 2021b). One of the most commonly employed UFMs is the Brunone UFM. Brunone et al. (1991) related pipe friction to the local instantaneous acceleration and convective acceleration. Vítkovský et al. (2006) and Bergant et al. (2001) improved the original Brunone model by extending the sign of the convective term, resulting in an improved model that is computationally simple and overlaps more completely with experimental results.

Zhou et al. (2021a, 2021b) have presented the GST model incorporating the unsteady friction factor for the transient pressurized pipe flows. So far, the unsteady friction factor is rarely considered in the free surface pipe flows. However, the previous work of the free surface pipe flows has demonstrated that the effect of friction is associated with the flow intensity, and different formulas are usually used to calculate the friction coefficient for different flow regimes (Lim 2018). Some measured data (Bazin 1865; Varwick 1945; Straub et al. 1958) prove that the open channel flow is laminar flow when the Reynolds number is less than 2,320, and in most cases, the flow conditions are found to be smooth turbulent. As the Reynolds number increases, the friction factor of open channel flow also decreases, similar to that of pipe flow. This is called smooth turbulent flow (Lim 2018; Qiang et al. 2019). The effect of unsteady friction on the free surface pipe flow is not clear.

Thus, the purposes of this paper are (1) to develop and apply the second-order GTS method to calculate the free surface flow considering Brunone unsteady friction and to compare it with existing models and experimental results and (2) to numerically explore the effects of unsteady friction on the free surface flow. The double virtual cells boundary treatment is proposed to ensure the whole computation domain with second-order accuracy. A sensitivity analysis of the numerical solution to changes of computational parameters (initial water depth, Courant number) is performed.

The pipe x-axis is discretized into N + 4 grids of fixed length Δx. The t-axis is discretized into grids of time interval Δt. For the free surface flow, each time step Δt is different. The ith control volume is centered at node i and extends from i− 1/2 to i + 1/2. Figure 2 shows a partial discrete grid diagram of the computational model.

In this paper, the MUSCL-Hancock scheme is introduced to reconstruct the GTS with second-order accuracy, and the MINMOD slope limiter is used to avoid spurious oscillations. The linear reconstruction steps of the MUSCL-Hancock method are as follows.

As a classical Riemann solver, HLL has been adopted in various fields for its simple calculation and reliable results. In this approach, the solution of the Riemann problem is approximated by an intermediate star region U_*, of constant state separated from the left and right states (U_L and U_R) by two infinitely thin waves. Figure 3 illustrates this approximation.

Equation (11) is used for the calculation of all computational domains, including upstream and downstream boundaries F_1/2 and F_N+1/2. For each F_i+1/2, the linear reconstruction using the MUSCL-Hancock scheme requires data from both the left and right sides of the interface. However, the control volume only includes cells 1 to N. Therefore, the second-order GTS proposed by Zhao & Ghidaoui (2004) has only first-order accuracy in solving the water hammer problem at the boundary. Leon et al. (2006) proposed to add a virtual cell on each side of the boundary, as shown in Figure 4(a), to achieve the second-order accuracy, but the fluxes of both boundaries still need to be handled separately. This section focuses on the simplicity of the second-order accuracy (both internal and boundary cells) of the double virtual cells method for free surface flow modeling in this paper. By adding two virtual grids on each side of the upstream and downstream boundaries, the second-order linear reconstruction can be performed to unify the boundary grids to achieve the second-order accuracy.

For a circular cross-section channel, if the inflow is to be prescribed at upstream boundary and and are unknown variables associated with θ. Therefore, there is only one unknown variable θ in the above equation, which can be solved by iteration, and then and can be found directly. The downstream boundary can also be found in the same way.

The second-order format of the F_i−1/2 and F_i+1/2 flux terms has been found in the previous section. In order to advance the solution from n to n + 1, time integration was required. In the case of friction, the source term S is introduced into the solution by using a second-order Runge–Kutta discretization for time segmentation, and the following explicit procedure is obtained.

This section has several objectives: (1) to validate the proposed model by comparing the MOC scheme; (2) to compare the results calculated by the proposed model and the four-point implicit Preissmann scheme; (3) to demonstrate the accuracy of GTS-UFM in modeling unsteady flow and to explore the effects of unsteady friction on the free surface flow; and (4) to conduct the parameter analysis, including initial water depth and Courant number Cr.

As shown in Figure 5, the simple system consists of a 1,000 m long frictionless pipe with 15 m inner diameter, and a gate in the middle of the pipeline. The initial moment is a stationary state with a water depth of 10 m upstream of the gate and 3 m downstream of the gate. The inflow of the upper boundary and outflow of the lower boundary are 0 m³/s, so the total mass in the sewer is constant with time. The change of water depth after the instantaneous full opening of the gate is simulated. Since friction is not considered in this case, it can be assumed that any attenuation generated in the simulation results is caused by numerical dissipation.

Figure 6 compares the water depth variation curves at x = 2.5 m for grid number N = 200 and Courant number Cr = 0.9 between the proposed model and the existing model (including the conventional boundary treatment proposed by Leon et al. (2006) and MOC). The results show that the simulation results of both GTS schemes gradually converge to the exact equilibrium level as time advances, while the simulation results of MOC are much smaller than the exact equilibrium level. The reason is that, when the Courant number is less than one, the interpolation of MOC will result in wave smoothing (damping). Damping reduces the peak and artificially delays the time when surcharging occurs. For the same Courant number, the GTS with the double virtual cells boundary method and the single virtual cell method proposed by Leon et al. (2006) have almost the identical accuracy, thus verifying the correctness and feasibility of the double virtual cells boundary method. It is also worth pointing out that the boundary treatment method of adding double virtual cells can realize the unification of the simulation of the whole calculation domain, simplify the calculation and programming steps, and thus improve the programming efficiency. In this case, the gate in the middle of the pipe suddenly opens and no water flows into or out of the system, so the total mass in the pipe is invariant in the transient event. As shown in Figure 7, the results show that the second-order GTS proposed in this paper can maintain mass conservation, and the same conclusion can be reached when steady and unsteady friction are considered, while the MOC method leads to the loss of water mass. In addition, when Cr = 0.9, the total water loss can reach 10%. Therefore, theoretically, the double virtual boundary treatment method and the second-order Godunov scheme proposed in this paper are fully self-consistent, obeying mass conservation, and are significantly better than the simulation results using the MOC method with the same parameters.

Figure 8 shows the water depth variation along the pipe at t= 36 s and t = 197 s, the calculation results of the proposed method are compared with the existing models and exact solutions. It can be seen from Figure 8 that the double virtual boundary treatment proposed in this work is more accurate than the computational results of MOC. In particular, when t = 197 s, the water depth variation of the GTS is almost consistent with the exact solution, while the simulated curve of MOC significantly lags behind the exact solution (about 100 s). An accurate reproduction of the jump speed in the sewer system is very significant because it determines the occurrence time of the surcharge.

In general, when the Courant number is less than 1.0, the MOC needs to be interpolated, which will lead to wave smoothing (damping) and wave front diffusion. The diffusion makes the wave reach the boundary earlier, while damping reduces the peak of the wave and artificially delays the time of overload occurrence. Therefore, MOC is not well suitable for real-time control of sewer system, and GTS proposed in this paper can well avoid the above two problems. Compared with Leon et al. (2006), the boundary treatment method proposed in this paper is simpler to calculate.

Case 2 considers a 1,500 m long trapezoidal channel. The channel has a prismatic cross-section with a bottom width of 10 m, a longitudinal slope of 0.0002, a side slope of 0.5, and a Manning roughness coefficient of 0.016. The channel delivers 40 m³/s of water to the hydroelectric scheme, and the flow in the channel is uniform during normal operation. The channel is open to the upstream reservoir, and when the plant is shut down, it is filled with stationary water at the reservoir level. The initial flow rate is 40 m³/s in the system, and the discharge is decreased linearly to 0.0 m³/s in 2 min.

The downstream transient water depth predicted by the proposed model are compared with those predicted by the four-point implicit Preissmann scheme and the MOC scheme, as shown in Figure 9. Ten reaches are used for three programs. In the four-point implicit Preissmann scheme, time step is 20 s and the weighting factor is Ψ = 0.6. As discussed in Case 1, the proposed model can well reproduce the results of exact solution, but MOC leads to serious calculation errors due to a large mass loss. It can be found in Figure 9 that transient water depth predicted the four-point implicit Preissmann scheme is closer to that of the proposed model, compared with the result of the MOC scheme. This means the four-point implicit Preissmann scheme is more accurate than the MOC scheme.

Figure 10(a) presents results from the four-point implicit Preissmann scheme analysis of Case 2. The effect of the time steps is investigated. The increased dissipation was obvious with the larger time steps in the four-point implicit Preissmann scheme, although the four-point implicit Preissmann scheme was stable. In addition, the influence of the weighting factor Ψ is also investigated. For this case, a value of Ψ = 0.6 is recommended. The numerical damping with Ψ = 1.0 is large. Therefore, the accuracy of the four-point implicit Preissmann scheme is limited by the time steps and weighting factors. However, for the method proposed in this paper, as shown in Figure 10(b), the number of reaches remains constant and the increase in dissipation is not significant as the Cr decreases. In addition, the method in this paper is as easy to operate as MOC when dealing with complex boundary conditions such as branches or nodes. In summary, the scheme proposed in this paper is more suitable for realistic open channel transient simulations.

To evaluate the proposed model, the free surface transient flow experiment of Ackers & Harrison (1964) is considered. The computational results of the second-order GTS-UFM and GTS-SFM are both compared with the experimental results of Ackers & Harrison (1964).

In the experiments of Ackers & Harrison (1964), the pipe length is 304.8 m with 0.3048 m inner diameter and 0.001 slope, and the Manning roughness coefficient is 0.0115. The Brunone unsteady friction is taken into account in this case, and the kinematic viscosity of the fluid is taken as 0.8937 × 10⁻⁶ m²/s. The initial base flow depth of the system is 0.0768 m, and the base flow is 0.004984 m³/s. The upstream boundary flow varied with time as shown in Figure 11. In this case, transients occurred in the free surface pipeline as the upstream inflow changed. In addition, the downstream boundary condition is assumed to have a free outlet, which is consistent with the laboratory setup. The whole installation is shown in Figure 12.

The model uses N = 300 cells and Cr = 0.3. Figure 13 shows the calculated and experimental values at 8.66 and 77.94 m from the upstream. As shown in Figure 13, the proposed model could accurately simulate the case results, including crest, periods, and decay, under the above case conditions. Compared with the results of GTS-SFM, the calculated results of GTS-UFM are closer to the experimental results. After adding the Brunone unsteady friction, the dispersion of the simulated data increases and is closer to the experimental value, that is, the real-time control accuracy of GTS-UFM is higher than that of GTS-SFM.

Table 1 presents the relative error (RE) values and the mean squared error (MRE) values of both GTS-SFM and GTS-UFM when Cr = 0.3 with the same grid number. Theoretically, the proposed models become more accurate as the RE values and the MRE values decrease. As shown in Table 1, the model developed in this paper is more accurate than the model that only considers steady friction.

Based on Case 1, this case considers friction factors, in which the Manning roughness coefficient n is 0.015 and the fluid viscosity is 0.893710⁻⁶ m²/s.

As shown in Figure 14, when the grid number N= 50 and the Courant number Cr = 0.9, the smoothing models (SMM), SFM, and UFM of the second-order GTS proposed in this paper are used to obtain the water depth versus time plots for x= 990 m, respectively. Figure 14 shows that the simulation results of the three models are almost identical at initial moment, but there are large differences in the simulation results of subsequent water depth cycles in the free water transient flow. The simulation results for the UFM decayed faster than the SMM and SFM and took the shortest time to reach the exact equilibrium level, with the entire system reaching steady state at 7,110.285 s. However, when the SMM and SFM simulate this system flow, the flow does not fully reach steady state even when calculated to 10⁵ s. This is because the Brunone UFM considered the real roughness of the actual engineering pipeline and the unsteady friction caused by the transient acceleration and convective acceleration. The smooth model and steady friction model may differ greatly in the time dimension when predicting the transient of free surface flow.

The initial flow velocity in Case 4 is 0.0 m/s. As the gate in the middle of the pipeline is opened instantaneously, the water at the high water depth flows downstream and drives the water in the whole pipeline to flow back and forth until the water depth reaches to a constant value. Different initial water depths could lead to different flow velocities and transient process. In order to analyze the effect of unsteady friction on the free surface transient flow at different initial water depths, the numerical simulations keep the initial downstream water depth = 3 m and upstream water depth = 13, 10, and 7 m.

Figure 15 shows the variation of water depth with time considering the UFM under different initial upstream water depth. Obviously, the equilibrium depth increases as the initial upstream water depth increases.

To analyze the effect of unsteady friction in simulating the unsteady flow on the free surface, the water depth difference of the results calculated by the steady friction model and the UFM under different initial upstream water depth. The water depth peaks during the first three oscillation cycles are given in Table 2. When the upstream water depth is 13 m, the maximum water depth variation caused by unsteady friction can reach 1.89 m, accounting for 23.96% of the system equilibrium depth. When the upstream initial water depth is 7 m, the maximum water depth variation caused by unsteady friction is 0.6475 m, accounting for 12.7% of the system equilibrium depth. These demonstrate that the unsteady friction factor has important effects on the water depth oscillations.

Figure 16(a) and 16(b) presents water depth variation and relative water depth variation caused by unsteady friction during the transient event. In Figure 16(a), the effect of unsteady friction increases first and then decreases gradually. This is because initial flow velocity is 0.0 m/s; when the valve is suddenly opened, the flow accelerates (velocity increases) and the effect of unsteady friction becomes important due to drastic flow variation; as the whole flow gradually tends to the equilibrium state, the effect of the unsteady friction decreases. In Figure 16(b), the relative water depth variation caused by the unsteady friction during the transient event is larger with increasing initial upstream water depth until 1,200 s, but after 1,200 s it becomes smaller with the higher initial upstream head. This is because the effect of the unsteady friction is larger in the beginning stage of the transient event, while it decreases as the system stabilizes with the advance of time. The relative water depth variation caused by unsteady friction during the transient event includes the effect of the equilibrium depth, and increases as the equilibrium depth becomes smaller in the case that the difference of unsteady friction is not large.

Figure 17 displays the depth results of Case 4 with N = 50 and Cr = 0.3, 0.5, and 0.9, predicted by the second-order GTS-UFM and the space-line interpolation MOC-UFM scheme. As expected, the numerical dissipation of both the MOC-UFM and the second-order GTS-UFM in this paper becomes more severe as Cr decreases. For MOC-UFM, with the variation of Cr, the numerical simulation results mainly have the following changes: (1) the numerical dissipation increases with the decrease of Cr, (2) as Cr decreases, the system cannot maintain mass conservation, especially when Cr is taken as 0.5–0.9, the mass loss increases sharply, (3) as the Cr decreases, the time phase deviation becomes larger. While GTS-UFM is used to simulate the free surface flow, with the decrease of Cr, there is no significant difference in the results of water depth simulation for each water depth cycle except for the maximum and minimum values with certain numerical dissipation. Moreover, the GTS-UFM does not have phase deviation and can also maintain mass conservation. Therefore, GTS-UFM is better in realizing the real-time control simulation of sewer system.

Traditionally, the MOC scheme or the four-point implicit Preissmann scheme considering the steady friction models has often been used to simulate free surface flow. In this paper, the second-order GTS combining the improved Brunone UFM is developed to simulate free surface transient events. The double virtual boundary cell method is proposed to deal with the boundary calculation problem. Results calculated by the GTS-UFM are compared with existing models as well as experimental data, the main conclusions are as follows:

1. The double virtual boundary cell method is proposed to establish a finite-volume second-order GTS for open channel flow. By adding two virtual cells at the boundary, the unified calculation of the whole computational domain is realized. Compared with the calculation results of MOC, the model established in this paper satisfies the mass conservation law and is more accurate. In contrast, the MOC method cannot guarantee mass conservation.

2. The accuracy of the four-point implicit Preissmann scheme is limited by the time steps and weighting factors. The four-point implicit Preissmann scheme was also found not only to be unsuitable for transcritical flows but also has complex modeling, large storage requirement, and long simulation time when solving for branches or junctions in drainage systems. Therefore, compared with the four-point implicit Preissmann scheme, the method proposed in this paper is more stable and easier to operate when dealing with a real network.

3. The numerical dissipation of both the MOC and the method proposed in this paper becomes more and more severe as the Cr decreases. However, when Cr changes to a very small level, the simulation results of the method proposed in this paper are still accurate, while MOC will have a large mass loss. Therefore, compared with MOC, the method proposed in this paper is more robust.

4. The second-order GTS considering the improved Brunone unsteady friction accurately predicts the measured water depth histories and water depth peaks of free surface transient pipe flow in the systems.

5. For free surface pipeline systems, the Brunone friction increases with the increase of initial water depth difference and initial flow velocity. On the other hand, the unsteady friction cannot be neglected even in the case of a small initial water depth difference.

This work was supported by the National Natural Science Foundation of China [Grant Nos. 51839008 and 51679066]; Fok Ying Tong Education Foundation [Grant No. 161068].

All relevant data are included in the paper or its Supplementary Information.