An accurate and efficient scheme involving unsteady friction for transient pipe flow

The following symbols are used in this paper:

a

wave speed

linearized coefficient matrix

A*

coefficient of original turbulent weighting function

B

coefficient matrix

B*

coefficient of original turbulent weighting function

C

coefficient matrix

Cr

Courant number

C*

Vardy's shear decay coefficient

D

pipe diameter

f

flux term

g

gravitational acceleration

H

piezometric head

H_L, H_R

piezometric head to left and right of interface

H_r

upstream reservoir pressure head

i

index for cell, reach

I_i

ith cell

J_Q

the head loss by the quasi-steady friction

J_u

the head loss due to unsteady friction

k

the unsteady friction coefficient in Brunone model

L

pipeline length

m_i

the weighting coefficient

the scaled universal weighting coefficient

n

index for time t

n_i

the weighting coefficient

the scaled universal weighting coefficient

N

number of cells along the pipeline

R

pipe radius

Re

Reynolds number

s

source term

t

time

u

flow variables H and V

U

cell mean value of u

intermediate flow variables in Runge–Kutta scheme

average cell value of u to the left and the right of interface at time nΔt

V

the average cross-sectional velocity

V_L, V_R

average cell velocity to the left and right of interface

V₀

initial steady pipe flow velocity

W (t)

a weighting function

x

distance along pipeline

Y_ai(t)

the weighting function

Δt

time increment

Δτ

4vΔt/D²

Δx

reach length

κ

coefficient of original turbulent weighting function

ρ

the density of the liquid

τ

the dimensionless time

τ_u

unsteady shear stress

ν

kinematic viscosity of the fluid

eigenvalues of

ɛ

pipe wall roughness

CFD

computational fluid dynamics

CFL

Courant–Friedrichs–Lewy

FDM

finite difference method

FVM

finite volume method

MINMOD

MINMOD slope limiter function

MOC

method of characteristics

MUSCL

monotone upstream-centered scheme for conservation laws

GTS

Godunov-type schemes

TVB

Trikha–Vardy–Brown

The complicated hydraulic transients often occur during the energy conversion process in the pipe system of hydropower station and pumped storage station. Dangerous water hammer events caused by some inappropriate operations of pump/turbine/valve in the water system and likely produce abnormally high pressures, which may induce pipe rupture and damage other hydraulic devices (Wylie et al. 1993). Therefore, the accurate and efficient numerical simulations of all the possible hydraulic events become more important for the proper design and safe operation of the real pipe systems.

Various numerical solution schemes, including method of characteristics (MOC) and finite difference method (FDM) (Chaudhry & Hussaini 1985), have been introduced in the traditional one-dimensional (1D) water hammer model to simulate the hydraulic transients for various pipe systems (Wylie et al. 1993). Among these methods, MOC proved to be popular among water hammer experts since it is simple, efficient, and accurate.

The real water pipe systems often contain various pipe sections with different materials and lengths. So, the MOC approach has to employ either interpolation or wavespeed adjustment in pipes since it is impossible to make the Courant number exactly equal to one in all pipes of a complex pipe system. Several significant and complicated interpolation techniques have been proposed in order to deal with the discretization problem associated with applying fixed-grid MOC to pipeline systems (Ghidaoui & Karney 1994; Karney & Ghidaoui 1997; Ghidaoui et al. 1998). Alternatively, Wylie & Streeter (1970) and Chaudhry & Hussaini (1985) solved the water hammer equations by FDM schemes. It was found that the second-order FDM schemes can get much better results than the first-order MOC; in particular, the implicit methods are stable for large time steps. However, implicit FDM schemes increase both the computation time and the storage requirement (Chaudhry & Hussaini 1985).

Moreover, the finite volume method (FVM) Godunov-type schemes (GTS) have been introduced for the solution of classic water hammer equations not involving the unsteady friction factor (Guinot 2000; Zhao & Ghidaoui 2004; Zhou et al. 2017, 2018a). The results show that the first-order FVM GTS is identical to the fixed-grid MOC scheme, and importantly, the second-order FVM GTS is more robust for simple water hammer events. However, it is necessary to further investigate the feasibility of GTS for more complicated hydraulic transient problems, including the incorporation of the unsteady friction.

Computational fluid dynamics (CFD) methods have also been applied to realize the three-dimensional (3D) flow simulations in the whole water system from the upstream and downstream reservoirs. Results show that the 3D CFD methods are feasible to accurately predict the transient pressures and energy transformation, and its significant advantage is to vividly reveal the physical process (Martins et al. 2014; Wang et al. 2016; Zhou et al. 2018b). However, it is extremely difficult to realize the fast dynamic simulating and monitoring by using the 3D CFD methods due to the extremely time-consuming and less efficient computations (Wang et al. 2016; Zhou et al. 2018b).

The main purpose of this work is to develop an accurate and efficient 1D numerical approach, which is significant for the dynamic prediction system monitoring all possible hydraulic transients in various real pipe systems, such as at the hydropower station and the pumped storage station. The second-order FVM GTS considering the unsteady friction factors was developed to realize the accurate and efficient simulations on the hydraulic transients, which is rarely involved in previous published work. One explicit-solution source item approach is proposed for the GTS to incorporate various forms of the existing unsteady friction models, including original accurate convolution-based models (Zielke model and Vardy–Brown model), the efficient simplified convolution-based models (Trikha–Vardy–Brown (TVB) model), and the Brunone instantaneous acceleration-based model. Results calculated by the proposed second-order GTS unsteady friction models are compared with published experimental data as well as predictions by the classical MOC scheme. The accuracy and efficiency of the proposed approach are discussed carefully.

The experimental pressures of Bergant et al. (2001) and Adamkowski & Lewandowski (2006) are used to verify the proposed models. The experimental system of Bergant et al. (2001) consists of an upstream constant pressure tank, an upward sloping straight pipeline with a 5.45% slope, a downstream ball valve, and a downstream pressure tank. Adamkowski & Lewandowski (2006) designed a similar experimental setup with higher upstream constant pressure and smaller pipe diameter. Figure 1 displays the experimental apparatus layout.

Initially, the pipe flow reaches a steady velocity V₀ by adjusting the downstream outlet pressures. All the experimental water hammer events of Bergant et al. (2001) and Adamkowski & Lewandowski (2006) were initiated by quickly closing the downstream ball valve. In the experiments of Bergant et al. (2001), upstream inlet pressure H_r = 32 m, pipeline length L = 37.23 m, pipe inner diameter D = 0.0221 m; wave speed a = 1,319 m/s; pipe wall thickness e = 1.63 mm and water temperature = 15.4 °C. In the experiments of Adamkowski & Lewandowski (2006), H_r = 127.47 m, L = 98.11 m, D = 0.016 m; a = 1,298.4 m/s; e = 1 mm and water temperature = 20 °C. The initial experimental conditions are summarized in Table 1.

In this work, the original convolution-based models, the simplified convolution-based model of Vardy & Brown (2004a, 2007), and the Brunone instantaneous acceleration-based model are all represented for the unsteady friction simulation of transient pipe flow.

The original convolution-based models are numerically accurate and nevertheless are computationally inefficient due to their continual involving all the past velocity. In order to improve computational efficiency, the Zielke model and Vardy–Brown model have been simplified by the introduction of the so-called effective weighting functions. Trikha (1975) carried out the first research to develop an effective method of solving the integral convolution. Later, Kagawa et al. (1983), Suzuki et al. (1991), Schohl (1993), Vardy & Brown (2004a, 2007), Vítkovský et al. (2004), Ghidaoui & Mansour (2002), and Urbanowicz & Zarzycki (2012) had improved the applicability range of effective weighting function as well as its degree of fit to the original weighting function.

The Trikha–Vardy–Brown (TVB) efficient simplified model (Vardy & Brown 2004a, 2007) with the ranges of 10⁻⁸ ≤Δτ < ∞ for laminar flow and 10⁻⁹ ≤Δτ ≤ 10⁻¹ for turbulent flow is introduced to simulate the unsteady friction of transient pipe flows since the possible range of Δτ is ranged from 10⁻⁷ to 10⁻⁵ according to the experimental parameters of Bergant et al. (2001) and Adamkowski & Lewandowski (2006), and the computational time step Δt (being associated with the grid length Δx, Δt=Cr·Δx/a).

For laminar flow, an efficient implementation of the Zielke weighting function can be realized by the effective weighting functions with {n_i, i= 1, …, 9} = {26.3744; 10²; 10^2.5; 10³; 10⁴; 10⁵; 10⁶;10⁷;10⁸} and {m_i, i= 1, …, 9} = {1; 2.1830; 2.7140; 7.5455; 39.0066; 106.8075; 359.0846; 1107.9295; 3540.6830}.

Similarly, for turbulent flow, Vardy & Brown (2004a, 2007) recommend the adoption of a scaling procedure proposed by Vitkovsky et al. (2004) to simplify the original Vardy–Brown weighting function, and presented the effective weighting functions with the scaled universal weighting coefficients (⁠⁠) and (⁠⁠) {⁠, i= 1, …, 17} = {10; 10^1.5; 10²; 10^2.5; 10³; 10^3.5; 10⁴; 10^4.5; 10⁵; 10^5.5; 10⁶; 10^6.5; 10⁷; 10^7.5; 10⁸; 10^8.5; 10⁹} and {⁠, i= 1, …, 17} = {9.06; −4.05; 12; 8.05; 22.7; 35.2; 65.9; 115; 206; 363; 664; 1,070; 2,060; 3,630; 6,640; 10,700; 26,200}.

and can be calculated by combining Equations (12) and (13). However, as shown in Equations (12) and (13), when Δt < Δx/a (Cr < 1) occurs in some complex pipe system with short pipe section or different pipe material, the space-line interpolation inevitably causes numerical damping, which is discussed later.

The estimation approach of and determines the accuracy order of the numerical scheme. In the first-order accuracy, and ⁠. Herein, the MUSCL-Hancock method is used to achieve second-order accuracy in space and time, while the MINMOD limiter is suggested to avoid the spurious oscillations. The details of the MUSCL-Hancock method and the MINMOD limiter can be found in Toro (2009).

For the upstream boundary (i = 0), coupling this Riemann invariant of Equation (17) with a head-flow boundary relation determines u_1/2(t) = (H_1/2, V_1/2), ⁠. Similarly, for the downstream boundary (i = N), coupling the Riemann invariant of Equation (18) with a head-flow boundary relation determines u_N+1/2(t) = (H_N+1/2, V_N+1/2), ⁠. Alternatively, virtual control volumes and I₀ adjacent to I₁, and virtual control volumes I_N+1 and I_N+2 adjacent to I_N, are used to realize the direct solution of the Riemann problem for the boundary unknowns U₁, U₂, U_N−1, and U_N. Here, ⁠, and ⁠, where and are calculated by combining the Riemann invariant with a head-flow boundary relation at time n.

The unsteady friction item is incorporated into the source terms s(u). Herein, a second-order Runge–Kutta discretization is introduced in the solution of the source terms s(u) of Equation (15), as the following explicit procedure.

The calculation procedure of the proposed algorithm is presented in Figure 3.

In the water hammer simulations of this work, when Δt = Δx/a, a little stability may occur due to the unsteady friction considered in the source terms s(u). Thus, the stability constraint is suggested, in which Cr near 0.9 is validated to get stable and considerably accurate results.

The specific purposes of this section are (1) validation of the second-order FVM GTS unsteady friction model by comparing the calculated and measured data; (2) exploration of the accuracy and efficiency of various unsteady friction models in the water hammer simulations by using the second-order FVM GTS; and (3) comparison of the stability and the accuracy of both the second-order FVM GTS and the MOC in the water hammer unsteady friction simulations by using different Courant number and grid number.

The original accurate convolution-based models (Zielke model and Vardy–Brown model) and the Brunone instantaneous acceleration-based model in both the proposed second-order FVM GTS and the MOC scheme are used to simulate the five experimental cases listed in Table 1. In the simulations, grid number N = 256; Cr = 0.9 and 1.0 are, respectively, used in the second-order FVM GTS and the MOC scheme. Figures 4 and 5 display the computed and observed pressure oscillation patterns of water hammer events in five cases. The experimental pressure was recorded by a pressure transducer at the pipe end.

The experimental Case 1 with initial velocity V₀ = 0.1 m/s and Re = 1870 is of laminar flow. So, Zielke model as the original accurate laminar-flow unsteady friction model is chosen to simulate the water hammer pressures in Case 1. Figure 4(a) shows that the quasi-steady friction model only predicts the first pressure peak, but seriously underestimates the pressure damping in the later pressure oscillations. The second-order GTS Zielke model can well reproduce the pressure attenuation. Moreover, the second-order GTS Brunone model seems to better predict the period of the whole pressure histories, but to produce a little excessive pressure damping.

Case 2 with Re = 3750 and Case 3 with Re = 5600 represent low Reynolds number turbulent flows. So, Vardy–Brown model as the original accurate turbulent-flow unsteady friction model is chosen to simulate the water hammer pressures in Cases 2 and 3. Compared with Case 1, the similar behavior and conclusion of both the quasi-steady friction model and the second-order GTS unsteady friction models can be found in Cases 2 and 3, as shown in Figure 4(b) and 4(c). In the cases of Reynolds number of a higher order of magnitude (10⁴) shown in Figure 4(d) and 4(e), the second-order GTS Vardy–Brown model and the Brunone model can still better reproduce the experimental pressure traces compared with the quasi-steady model.

To further validate the proposed second-order GTS, the classic MOC scheme is also introduced to simulate the experimental cases in Table 1. Figure 5 gives the pressure oscillations predicted by Zielke model, Vardy–Brown model, and the Brunone model in the MOC scheme. Comparisons of the results in Figures 4 and 5 indicate that the proposed second-order GTS unsteady friction models can produce the nearly identical pressure oscillations with the classic MOC unsteady friction models.

To quantificationally evaluate the proposed models, the root-mean-squared error (RMSE) and Nash–Sutcliffe efficiency coefficient (NSE, also known as R²) are used to identify the error magnitude. Table 2 shows simulation error values in numerical simulations of experimental Cases 1 and 3 (representing the laminar and turbulent flows, respectively). In theory, as the RMSE value decreases and the NSE is close to 1, the proposed model becomes more accurate. Results in Table 2 also show that (1) the proposed second-order GTS has the nearly identical accuracy with the classic MOC and (2) the unsteady models (Vardy–Brown model and Brunone model) can obtain more accurate results than the quasi-steady model.

Overall, the proposed second-order GTS unsteady friction models are capable of well reproducing the experimental pressure oscillations and accurately predicting the pressure damping of the water hammer pipe flow.

Some efficient simplified convolution-based models have been developed to approximate the original accurate convolution-based models (Zielke model and Vardy–Brown model). However, the simplified models have different application range associated with values of ⁠. The efficient TVB model (Vardy & Brown 2004a, 2007) with the ranges of 10⁻⁸ ≤ Δτ < ∞ for laminar flow and 10⁻⁹ ≤Δτ ≤ 10⁻¹ for turbulent flow is preliminarily chosen for the current simulations, since the range of is around from 10⁻⁷ to 10⁻⁵ in this work. Figure 6(a)–6(e) gives the numerical pressure oscillations of Cases 1–5 computed by the TVB model, which are compared with those calculated by the Zielke model and the Vardy–Brown model.

Figure 6(a)–6(e) indicates that the efficient TVB model has the identical accuracy with the original accurate Zielke model and the Vardy–Brown model in the simulations of water hammer events. Importantly, the TVB model is much more efficient than the original convolution-based models. For example, Lenovo ThinkPad T440 with Intel(R) Core(TM) i7-4510 U CPU and 8 GB RAM were used here to run the Fortran Codes of simulating five experimental cases. The computational time for different models in the numerical simulation of three cases is shown in Table 3. The TVB model takes only a few seconds to calculate the results in each case of Figure 4, whereas the Zielke model or Vardy–Brown Model consumes more than 1 h.

Consequently, the second-order GTS TVB model is accurate and efficient, and can represent the original accurate convolution-based models (Zielke model and Vardy–Brown model) to describe the unsteady friction in the simulations of transient pipe flow.

Courant number Cr and grid number N are the most important factors in determining the numerical accuracy. The experimental Case 1 was simulated by both the TVB model and the Brunone model in the solutions of both the proposed second-order GTS and the MOC scheme, in which N= 32 or 256 computational reaches and Cr = 0.1, 0.5, 0.9, and 1.0 (Cr = aΔt/Δx).

Figure 7 gives the pressure oscillations simulated by the MOC TVB model with N = 256 and 32, Cr = 0.1, 0.5, and 1.0. It can be found in Figure 7(a) that as Courant number decreases from 1.0 to 0.1, obvious numerical damping occurs after the first several periods of the pressure oscillations. In particular, when the coarse grid is used (N = 32 in Figure 7(b)), the numerical attenuation in the pressure results becomes more severe. The same behavior of numerical damping also arises in the simulation results of the MOC Brunone model, as shown in Figure 8(a) and 8(b). Actually, the numerical damping is attributed to the space-line interpolation of the fixed-grid MOC scheme. As shown in Equations (12) and (13), when Cr = aΔt/Δx < 1, the linear interpolation approach is used to calculate the variables at the MOC nodes at n time.

Figures 9 and 10 display the pressure results of Case 1 predicted by both the second-order GTS TVB model and the second-order GTS Brunone model with N = 256 and 32, Cr = 0.1, 0.5, and 0.9. It can be found in Figures 9(a) and 10(a) that the pressure oscillations of Cr = 0.1 and 0.5 are nearly consistent with that of Cr = 0.9 in the simulations of both the second-order GTS TVB model and the second-order GTS Brunone model. Even when coarse grid and small Courant number are used (in Figures 9(b) and 10(b)), the second-order GTS only produces slight numerical attenuation. Therefore, when Courant number is less than one, for the given Courant number and grid number, the second-order GTS is much more accurate than the MOC scheme.

Importantly, the second-order GTS is much more efficient than the MOC scheme. Table 4 presents the RMSE values and NSE coefficients of both MOC and the second-order GTS Brunone model when Cr = 0.5 with different grid numbers. It can be found that, for example, the second-order GTS with a sparse grid (N = 32) and small Courant number (Cr = 0.5) in Figure 10(b) can produce nearly consistent accuracy with the MOC scheme with the fine grid (N = 256) in Figure 8(a).

In all the numerical simulations of this paper, Lenovo ThinkPad T440 with Intel(R) Core(TM) i7-4510 U CPU and 8 GB RAM was used here to run the Fortran Codes of simulating all the cases. The computational time in the numerical simulations of Case 1 is shown in Table 5. When Courant number is less than one, for example Cr = 0.1 or 0.5 as shown in Table 5, for the same accuracy, the efficiency of the second-order GTS Brunone model is nearly 20 times that of the MOC Brunone model. More surprisingly, the second-order GTS TVB model is more than 100 times faster than the MOC TVB model.

Overall, for the water hammer simulations, the second-order GTS is accurate, efficient, and stable for Courant number less than one, even when a sparse grid is used. This is so important for the complicated pipe system, since it is difficult to ensure that Courant number is exactly equal to one, which inevitably causes the numerical dissipation in the traditional MOC scheme.

Numerical efficiency is also an important factor to evaluate the numerical model. Table 6 shows the computational time of Brunone model and TVB model in both the MOC scheme and the second-order GTS.

As shown in Table 6, for the given simulation conditions (Cr and N), the Brunone model is more efficient than the TVB model in both the MOC scheme and the second-order GTS. The main reason is that the computation of the weighting functions (Equation (9)) in the TVB model takes more time and storage than Brunone model.

Interestingly, for the given Cr and N, the second GTS TVB model is two times faster than the MOC TVB model, but the second GTS Brunone model is a little slower than the MOC Brunone model. The main reasons are (1) in the MOC scheme, each calculation of the pressure and velocity at i node involves twice computations of the unsteady friction items at i − 1 and i + 1 nodes, as shown in Equations (12) and (13); (2) in the Godunov scheme, each calculation of the pressure and velocity at i cell only involves once computation of the unsteady friction items at i cell, as shown in Equations (20) and (21); (3) the MUSCL-Hancock method of the Godunov scheme takes more computations than the Brunone formulation, but far less than the weighting function of the TVB model.

Importantly, as discussed in the previous section, as shown in Table 5, for the same accuracy, both the Brunone model and the TVB model in the second-order GTS are much faster than those in the MOC scheme.

The second-order FVM GTS unsteady friction models are developed to simulate the water hammer events in pipe systems. The explicit-solution source item approach is proposed for the GTS to easily incorporate various forms of the existing unsteady friction models, including original accurate convolution-based models (Zielke model and Vardy–Brown model), the efficient simplified convolution-based models (TVB model), and the Brunone instantaneous acceleration-based model. Results calculated by the proposed second-order GTS unsteady friction models are compared with published experimental data as well as predictions by the classical MOC scheme. The main conclusions are:

1. The explicit-solution source item approach is simple and feasible to incorporate various forms of the existing unsteady friction models in the water hammer solution of the second-order FVM GTS.

2. Both the experimental data and the results of the classic MOC scheme validate that the proposed second-order GTS unsteady friction models are capable of well reproducing the experimental pressure oscillations and accurately predicting the pressure damping of the water hammer pipe flow.

3. For the convolution-based unsteady friction models in the second-order Godunov scheme, the simplified TVB model has the identical accuracy with the original accurate Zielke model and the Vardy–Brown model in the current simulations of water hammer events. The proposed second GTS TVB model is much more efficient than the original convolution-based models.

4. The proposed second-order GTS unsteady friction models are accurate, efficient, and stable even for Courant number less than one. For the given Courant number and the same accuracy, the second-order GTS unsteady friction models are far more efficient than the MOC unsteady friction models.

5. The Brunone instantaneous acceleration-based model is computationally faster than the accurate convolution-based models in both the second-order GTS and the classic MOC scheme.

The writers gratefully acknowledge the financial support for this research from the Open Research Fund Program of State Key Laboratory of Water Resources and Hydropower Engineering Science (Wuhan University) (Grant No. 2016SDG01), National Natural Science Foundation of China (Grant Nos. 51839008 and 51679066), the Fundamental Research Funds for the Central Universities (Grant No. 2018B43114), and Fok Ying Tong Education Foundation (Grant No. 161068).

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