A finite-volume second-order Godunov-type scheme (GTS) combining the unsteady friction model (UFM) is introduced to simulate free surface flow in pipelines. The exact solution to the Riemann problem calculates the mass and momentum fluxes while considering the Brunone unsteady friction factor. One simple boundary treatment with double virtual cells is proposed to ensure the whole computation domain with second-order accuracy. Results of various transient free-surface flows achieved by the proposed models are compared with exact solution, experimental data, the four-point implicit Preissmann scheme solution, as well as predictions by the classic Method of Characteristics (MOC). Results show that the proposed second-order GTS UFMs are accurate, efficient, and stable even for Courant numbers less than one and sparse grid. The four-point implicit Preissmann scheme may produce severe numerical attenuation in the case of large time steps and unsuitable weighting factors, while the MOC scheme may produce severe numerical attenuation in the case of a low Courant number and could not maintain mass conservation. The numerical simulations considering the unsteady friction factor are closer to the measured water depth variations. The effect of unsteady friction becomes more important as the initial water depth difference increases significantly.

  • A finite-volume second-order Godunov-type scheme combining the unsteady friction model is introduced to simulate free surface flow in pipelines.

  • One simple boundary treatment with double virtual cells is proposed to ensure the whole computation domain with second-order accuracy.

  • The accuracy and stability of various numerical models (including the proposed model, MOC scheme, and Preissmann scheme) are investigated.

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

S0

slope of the pipe

Sf

pipeline friction

Ss

steady friction

Su

unsteady friction

v

water velocity

nm

Manning roughness coefficient

k

Brunone friction coefficient

C*

Vardy's shear decay coefficient

Re

Reynolds number

MINMOD

MINMOD slope limiter function

UL, UR

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

FL, FR

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

SL, SR

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

AL, AR

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

JS

head losses per unit length due to steady friction

JU

head losses per unit length due to unsteady friction

t

time

x

distance

Δx

space step

Δt

time step

Δτ

4vΔt/D2

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

Governing equations

The continuity and momentum equations of one-dimensional free surface flow can be written as (Chaudhry 1987):
(1)
where the variable U, the flux F, and the source term S can be written as
(2)
where A is the cross-sectional area of water; Q is the flow rate; is the average pressure of the water column on the cross-sectional area; ρ is the density of liquid; g is the acceleration of gravity; S0 is the slope of the pipe; and Sf is the channel friction.
For a circular channel (Figure 1), the cross-sectional area A and the item can be expressed as
(3)
(4)
where d is the diameter of the channel; θ is the central angle of the water surface; and y is the water depth.

Brunone unsteady friction model

The friction item contains the steady friction and unsteady friction, which can be expressed as
(5)
where Ss is steady friction and Su is unsteady friction items. Based on the Manning formula, the steady friction can be defined as
(6)
where n is the pipe Manning roughness coefficient and R is the hydraulic radius.
The improved Brunone formulation proposes the following formula for calculating the unsteady friction term (Bergant et al. 2001; Vítkovský et al. 2006),
(7)
in which sign(v) = (+1 for v ≥ 0 or −1 for v < 0). The Brunone friction coefficient k can be predicted by analytically using Vardy's shear decay coefficient C*,
(8)
where C* depends on the flow regime. For Reynolds number Re < 2,320, the flow state is laminar flow:
(9)
For Reynolds number Re > 2,320, the flow state is turbulent flow:
(10)

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.

Figure 1

Definition of variables in circular cross-sections.

Figure 1

Definition of variables in circular cross-sections.

Close modal
Figure 2

Discretization grids.

Figure 2

Discretization grids.

Close modal
By integrating Equation (1) from control surface i 1/2 to control surface i + 1/2 in the x direction, it can be obtained,
(11)
where the superscripts n and n + 1 denote time t and t + Δt, respectively.

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.

First step: initial data reconstruction
(12)
Second step: evolution
(13)

HLL approximate Riemann solver

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 (UL and UR) by two infinitely thin waves. Figure 3 illustrates this approximation.

Figure 3

Principle of the HLL Riemann solver.

Figure 3

Principle of the HLL Riemann solver.

Close modal
Combined with the conservation law, the flux Fi+1/2 in the HLL Riemann solver is given.
(14)
where SL and SR are the wave speed of the left and right waves; FL and FR are the fluxes on the left and right sides of the interface i + 1/2, respectively. UL and UR are the variables on the left and right sides of the interface i + 1/2.
(15)
In this method, ML and MR are obtained by assuming that the star region and the respective constant state are separated by shock waves. The left or right wave is a rarefaction wave if A*AK (K = L, R), it is suggested to replace Mk with the gravity wave celerity cK (K = L, R), otherwise the wave is a shock wave. Thus, Mk can be expressed as
(16)
where A* is the estimated value of the star area. cK (K = L, R) is the velocity of the gravity wave.
As for the estimation of A*, Toro (2001) has proposed various methods. In this paper, A* is given by solving the Riemann problem for the linearized hyperbolic system ∂U/∂t + ∂F(U)/∂x = 0. The two eigenvalues for are and . The application of the Rankine–Hugoniot condition across the two eigenvalues waves gives the following relation:
(17)
For a circular channel (Figure 1), the gravity wave celerity c can be expressed as
(18)
where ϕ = ∫(c/A)dA. For a circular cross-section, ϕ is given by
(19)

Double virtual boundary approach

Equation (11) is used for the calculation of all computational domains, including upstream and downstream boundaries F1/2 and FN+1/2. For each Fi+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.

Figure 4

Control volume and virtual boundary. (a) Leon's virtual boundary and (b) double virtual boundary.

Figure 4

Control volume and virtual boundary. (a) Leon's virtual boundary and (b) double virtual boundary.

Close modal
In order to achieve the second-order accuracy of the boundaries, the double virtual cells method is introduced. The parameters of ghost cells are the same as those of the boundary, so that the introduction of ghost cells does not introduce other disturbances while improving the accuracy. The double virtual cells method is given by
(20)
All variables in the computational domain −1 to N + 2 are known at the previous time, and a second-order linear reconstruction is performed for elements 1 to N. According to the Riemann invariant theory, the variables on the boundary can be obtained. Taking the subcritical flow as an example, the incoming flow or water depth at the upstream boundary is known, and the remaining variables can be obtained according to the following equation.
(21)

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.

Incorporation of source terms

The second-order format of the Fi−1/2 and Fi+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.

First step:
(22)
Second step:
(23)
Last step:
(24)

Stability constraints

The stability constraint includes the Courant–Friedrichs–Lewy (CFL) criterion for convection and the constraint on the source term. The CFL constraint is given by
(25)
It can also be seen that the allowable time step of the source term (Δtmax,S) is (some details of these derivations are presented in the Supplementary Material)
(26)
Finally, the maximum permissible time step is given by:
(27)

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.

Case 1: comparison of the proposed model and MOC scheme

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 m3/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 5

Schematic diagram of Case 1.

Figure 5

Schematic diagram of Case 1.

Close modal

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 6

Water depth at x= 2.5 m for Case 1.

Figure 6

Water depth at x= 2.5 m for Case 1.

Close modal
Figure 7

Mass traces for Case 1.

Figure 7

Mass traces for Case 1.

Close modal

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.

Figure 8

Water depth profile for Case 1. (a) t = 36 s and (b) t= 197 s.

Figure 8

Water depth profile for Case 1. (a) t = 36 s and (b) t= 197 s.

Close modal

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: comparison of the proposed model and the four-point implicit Preissmann scheme

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 m3/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 m3/s in the system, and the discharge is decreased linearly to 0.0 m3/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 9

Water depth at x = 1,500 m for Case 2.

Figure 9

Water depth at x = 1,500 m for Case 2.

Close modal

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.

Figure 10

Water depth profile for Case 2. (a) Preissmann four-point implicit scheme and (b) the proposed model.

Figure 10

Water depth profile for Case 2. (a) Preissmann four-point implicit scheme and (b) the proposed model.

Close modal

Case 3: comparison with the experimental data

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−6 m2/s. The initial base flow depth of the system is 0.0768 m, and the base flow is 0.004984 m3/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.

Figure 11

Inflow hydrograph for Case 3.

Figure 11

Inflow hydrograph for Case 3.

Close modal
Figure 12

Experimental installation of Ackers and Harrison.

Figure 12

Experimental installation of Ackers and Harrison.

Close modal

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.

Figure 13

Water depth at x= 8.66 m and x = 77.94 m for Case 3. (a) x= 8.66 m and (b) x= 77.94 m.

Figure 13

Water depth at x= 8.66 m and x = 77.94 m for Case 3. (a) x= 8.66 m and (b) x= 77.94 m.

Close modal

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.

Table 1

Error values when Cr = 0.3

Locationx = 8.66 m
x = 77.94 m
ModelGTS-SFMGTS-UFMGTS-SFMGTS-UFM
Relative error (%) 4.23 4.02 1.37 1.04 
Mean squared error 0.313 0.302 0.0004 0.0002 
Locationx = 8.66 m
x = 77.94 m
ModelGTS-SFMGTS-UFMGTS-SFMGTS-UFM
Relative error (%) 4.23 4.02 1.37 1.04 
Mean squared error 0.313 0.302 0.0004 0.0002 

Case 4: effect of unsteady friction on the free surface flow

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−6 m2/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 105 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.

Figure 14

Water depth for Case 4 at x= 990 m.

Figure 14

Water depth for Case 4 at x= 990 m.

Close modal

Influence of initial water depth

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.

Figure 15

Simulation results of unsteady friction for Case 4 at x= 990 m.

Figure 15

Simulation results of unsteady friction for Case 4 at x= 990 m.

Close modal

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.

Table 2

Water depth variation and relative water depth variation caused by unsteady friction

Initial upstream water depth (m)Equilibrium depth y_s (m)PeaksWater depth (SFM) y1 (m)Water depth (UFM) y2 (m)Water depth variation by UF y1 − y2 (m)Relative water depth variation by UF (y1 − y2)/y_s
13 7.87 First 9.58 8.99 0.60 7.58% 
Second 8.32 6.44 1.89 23.96% 
Third 8.66 7.25 1.42 17.99% 
10 6.62 First 6.30 5.91 0.38 5.81% 
Second 6.59 5.41 1.18 17.87% 
Third 7.55 6.53 1.02 15.48% 
5.10 First 4.48 4.28 0.19 3.81% 
Second 5.81 5.50 0.31 6.05% 
Third 5.44 4.80 0.65 12.69% 
Initial upstream water depth (m)Equilibrium depth y_s (m)PeaksWater depth (SFM) y1 (m)Water depth (UFM) y2 (m)Water depth variation by UF y1 − y2 (m)Relative water depth variation by UF (y1 − y2)/y_s
13 7.87 First 9.58 8.99 0.60 7.58% 
Second 8.32 6.44 1.89 23.96% 
Third 8.66 7.25 1.42 17.99% 
10 6.62 First 6.30 5.91 0.38 5.81% 
Second 6.59 5.41 1.18 17.87% 
Third 7.55 6.53 1.02 15.48% 
5.10 First 4.48 4.28 0.19 3.81% 
Second 5.81 5.50 0.31 6.05% 
Third 5.44 4.80 0.65 12.69% 

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 16

Time variation curves of unsteady friction for Case 4 at x= 990 m. (a) Water depth variation by UF and (b) relative water depth variation by UF.

Figure 16

Time variation curves of unsteady friction for Case 4 at x= 990 m. (a) Water depth variation by UF and (b) relative water depth variation by UF.

Close modal

Influence of courant number

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.

Figure 17

Water depth profile for Case 4 at x = 10 m. (a) GTS-UFM and (b) MOC-UFM.

Figure 17

Water depth profile for Case 4 at x = 10 m. (a) GTS-UFM and (b) MOC-UFM.

Close modal

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.

Ackers
P.
&
Harrison
A. J. M.
1964
Attenuation of flood waves in part-full pipes
.
Proc. Inst. Civ. Eng.
28
,
361
382
.
Alcrudo
F.
,
Garcia-Navarro
P.
&
Saviron
J.-M.
1992
Flux-difference splitting for 1D open channel flow equations
.
Int. J. Numer. Methods Fluids
14
,
1009
1018
.
Bazin
H. E.
1865
Recherches experimentales sur lecoulement deleau dans les canaux decouverts. In Memories presents par diverssavants al Academie des Sciences (Paris), p. 19 (in French)
.
Bergant
A.
,
Simpson
A. R.
&
Vitkovský
J. P.
2001
Developments in unsteady pipe flow friction modelling
.
J. Hydraul. Res.
39
(
3
),
249
257
.
Brunone
B.
,
Golia
U. M.
&
Greco
M.
1991
Some remarks on the momentum equations for fast transients
. In
Proc., Int. Meeting on Hydraulic Transients with Column Separation (9th and Last Round Table of the IAHR Group)
.
International Association for Hydro-Environment Engineering and Research
,
Valencia
,
Spain
, pp.
201
209
.
Chaudhry
M. H.
1987
Applied Hydraulic Transients
.
Springer
,
Berlin
.
Chow
V. T.
1959
Open Channel Hydraulic
.
McGraw-Hill Publishing
,
New York
.
Daily
J. W.
,
Hankey
W. L.
,
Olive
R. W.
&
Jordaan
J. M.
1956
Resistance coefficients for accelerated and decelerated flows through smooth tubes and orifices
.
Trans. ASME
78
(
7
),
1071
1077
.
Duan
H.-F.
,
Pan
B.
,
Wang
M.
,
Chen
L.
,
Zheng
F.
&
Zhang
Y.
2020
State-of-the-art review on the transient flow modeling and utilization for urban water supply system (UWSS) management
.
J. Water Supply Res. Technol.
69
(
8
),
858
893
.
https://doi.org/10.2166/aqua.2020.048
.
French
R. H.
1985
Open Channel Hydraulic
.
McGraw-Hill Publishing
,
New York
.
Fujihara
M.
&
Borthwick
A. G. L.
2000
Godunov-type solution of curvilinear shallow water equations
.
J. Hydraul. Eng.
126
(
11
),
827
836
.
Ghidaoui
M. S.
&
Mansour
S.
2002
Efficient treatment of Vardy-Brown unsteady shear in pipe transients
.
J. Hydraul. Eng.
128
(
1
),
102
112
.
Henderson
F. M.
1966
Open Channel Hydraulic
.
McGraw-Hill Publishing
,
New York
.
Hirsch
C.
1990
Numerical Computation of Internal and External Flows
, Vols.
1 and 2
.
Wiley
,
New York
.
Leon
A. S.
,
Ghidaoui
M. S.
,
Arthur
R. S.
&
Marcelo
H. G.
2006
Godunov-type solutions for transient flows in sewer
.
J. Hydraul. Eng.
132
(
8
),
800
813
.
Lim
H. S.
2018
Open channel flow friction factor: logarithmic law
.
J. Coastal Res.
34
(
1
),
229
237
.
Negharchi
S. M.
&
Shafaghat
R.
2021
The numerical development of MOC for analyzing the inclined pipelines using the experimental network of Babol Noshirvani University as a case study
.
J. Water Supply Res. Technol.
70
(
5
),
741
756
.
https://doi.org/10.2166/aqua.2021.002
.
Qiang
Z.
,
Fengchuan
Z.
,
Yuchen
Y.
&
Zhaoyu
D.
2019
Diagnostic function analysis of the logarithmic law in open channel turbulence
.
J. Tsinghua Univ: Nat. Sci. Ed. R.
59
(
12
),
999
1005
.
Retsinis
E.
,
Daskalaki
E.
&
Papanicolaou
P.
2020
Dynamic flood wave routing in prismatic channels with hydraulic and hydrologic methods
.
J. Water Supply Res. Technol.
69
(
3
),
276
287
.
https://doi.org/10.2166/aqua.2019.091
.
Sart
C.
,
Baume
J. P.
,
Malaterre
P. O.
&
Guinot
V.
2010
Adaptation of Preissmann's scheme for transcritical open channel flows
.
J. Hydraul. Res.
48
(
4
),
428
440
.
Straub
L. G.
,
Silberman
E.
&
Nelson
H. C.
1958
Open channel flow at small Reynolds numbers
.
Trans. ASCE
123
,
685
706
.
Toro
E. F.
2001
Shock-Capturing Methods for Free-Surface Shallow Flows
.
Wiley
,
Chichester
,
UK
.
Trikha
A. K.
1975
An efficient method for simulating frequency-dependent friction in transient liquid flow
.
J. Fluids Eng.
97
(
1
),
97
105
.
https://doi.org/10.1115/1.3447224
.
Urbanowicz
K.
2018
Fast and accurate modelling of frictional transient pipe flow
.
Z. Angew. Math. Mech.
98
(
5
),
802
823
.
https://doi.org/10.1002/zamm.201600246
.
Vardy
A. E.
&
Brown
J. M. B.
1995
Transient, turbulent, smooth pipe friction
.
J. Hydraul. Res.
33
(
4
),
435
456
.
https://doi.org/10.1080/00221689509498654
.
Vardy
A. E.
&
Brown
J. M. B.
2003
Transient turbulent friction in smooth pipe flows
.
J. Sound Vib.
259
(
5
),
1011
1036
.
https://doi.org/10.1006/jsvi.2002.5160
.
Vardy
A. E.
&
Brown
J. M. B.
2004
Transient turbulent friction in fully rough pipe flows
.
J. Sound Vib.
270
(
1–2
),
233
257
.
https://doi.org/10.1016/S0022-460X(03)00492-9
.
Varwick
F.
1945
Zur Fließformel für offene Künstliche Gerinne
.
Dresden University
,
Dresden
,
Germany
(in German)
.
Vítkovský
J. P.
,
Bergant
A.
,
Simpson
A. R.
&
Lambert
M. F.
2006
Systematic evaluation of one-dimensional unsteady friction models in simple pipelines
.
J. Hydraul. Eng.
132
(
7
),
696
708
.
Wan
W.
&
Huang
W.
2018
Water hammer simulation of a series pipe system using the MacCormack time marching scheme
.
Acta Mech.
229
,
3143
3160
.
https://doi.org/10.1007/s00707-018-2179-2
.
Wang
Y.
,
Zhang
C.
,
Li
Z.
,
Sun
B.
&
Zhou
H.
2019
Applicability of Preissmann box scheme for calculation of transcritical flow in pipes
.
Water Supply
19
(
5
),
1429
1437
.
https://doi.org/10.2166/ws.2019.010
.
Yen
B. C.
2001
Hydraulics of Sewer Systems. Stormwater Collection Systems Design Handbook
.
McGraw-Hill
,
New York
.
Zhang
B.
,
Wan
W.
&
Fan
L.
2019
Investigation into complex boundary solutions of water filling process in pipeline systems
.
J. Water
11
(
4
),
641
.
https://doi.org/10.3390/w11040641
.
Zhao
M.
&
Ghidaoui
M. S.
2004
Godunov-type solutions for water hammer flows
.
J. Hydraul. Eng.
130
(
4
),
341
348
.
Zhou
L.
,
Li
Y.
,
Karney
B.
,
Cheng
Y.
&
Liu
D.
2021a
Godunov-type solutions for transient pipe flow implicitly incorporating Brunone unsteady friction
.
J. Hydraul. Eng.
147
(
7
),
04021021
.
Zhou
L.
,
Cao
Y.
,
Karney
B.
,
Vasconcelos
J. G.
,
Liu
D.
&
Wang
P.
2021b
Unsteady friction in transient vertical-pipe flow with trapped air
.
J. Hydraul. Res.
59
(
5
),
820
834
.
doi:10.1080/00221686.2020.1844808
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Supplementary data