Abstract

A robust prediction system should monitor all possible hydraulic transients, which is significant for the appropriate and safe operation of pipe systems. A second-order finite volume method (FVM) Godunov-type scheme (GTS) considering unsteady friction factors is introduced to simulate hydraulic transients, which was rarely involved in previous work. One explicit-solution source item approach developed in this work is crucial for the proposed GTS to easily incorporate various forms of the existing unsteady friction models, including original convolution-based models (Zielke model and Vardy–Brown model), simplified convolution-based model (Trikha–Vardy–Brown (TVB) model), and Brunone instantaneous acceleration-based model. Results achieved by the proposed models are compared with experimental data as well as predictions by the classic Method of Characteristics (MOC). Results show that the MOC scheme may produce severe numerical attenuation in the case of a low Courant number. The proposed second-order GTS unsteady friction models are accurate, efficient, and stable even for Courant numbers less than one and sparse grid, and only need much less grid number and computation time to reach the same numerical accuracy. The TVB convolution-based model and Brunone model in the second-order GTS are suggested for further applications in hydraulic transients due to their high accuracy and efficiency.

HIGHLIGHTS

  • A Godunov-type scheme involving unsteady pipe friction is developed.

  • An explicit-solution source item approach is introduced to incorporate unsteady friction.

  • The numerical models are validated by the experimental data.

  • The proposed models are accurate, efficient, and stable for transient pipe flow.

  • A Brunone unsteady friction model is more efficient than the convolution-based model.

NOTATION

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

  •  
  • HL, HR

    piezometric head to left and right of interface

  •  
  • Hr

    upstream reservoir pressure head

  •  
  • i

    index for cell, reach

  •  
  • Ii

    ith cell

  •  
  • JQ

    the head loss by the quasi-steady friction

  •  
  • Ju

    the head loss due to unsteady friction

  •  
  • k

    the unsteady friction coefficient in Brunone model

  •  
  • L

    pipeline length

  •  
  • mi

    the weighting coefficient

  •  
  • the scaled universal weighting coefficient

  •  
  • n

    index for time t

  •  
  • ni

    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

  •  
  • VL, VR

    average cell velocity to the left and right of interface

  •  
  • V0

    initial steady pipe flow velocity

  •  
  • W (t)

    a weighting function

  •  
  • x

    distance along pipeline

  •  
  • Yai(t)

    the weighting function

  •  
  • Δt

    time increment

  •  
  • Δτ

    4vΔt/D2

  •  
  • Δ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

ABBREVIATIONS

     
  • 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

INTRODUCTION

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.

MATERIALS AND METHODS

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.

Figure 1

Schematic diagram of the experimental setup in Bergant et al. (2001) and Adamkowski & Lewandowski (2006).

Figure 1

Schematic diagram of the experimental setup in Bergant et al. (2001) and Adamkowski & Lewandowski (2006).

Initially, the pipe flow reaches a steady velocity V0 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 Hr = 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), Hr = 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.

Table 1

Conditions for the experimental cases

Case no.V0 (m/s)ReFlow regimeData source
0.1 1,870 Laminar flow Bergant et al. (2001)  
0.2 3,750 Turbulent flow 
0.3 5,600 Turbulent flow 
0.631 10,634 Turbulent flow Adamkowski & Lewandowski (2006)  
0.94 15,843 Turbulent flow 
Case no.V0 (m/s)ReFlow regimeData source
0.1 1,870 Laminar flow Bergant et al. (2001)  
0.2 3,750 Turbulent flow 
0.3 5,600 Turbulent flow 
0.631 10,634 Turbulent flow Adamkowski & Lewandowski (2006)  
0.94 15,843 Turbulent flow 

WATER HAMMER EQUATIONS INVOLVING UNSTEADY FRICTION

The classic water hammer equations considering the friction factors can be written as (Wylie et al. 1993):
formula
(1)
formula
(2)
where H is the piezometric head; V is the average flow velocity; a is the wave speed; g is the gravitational acceleration; x is the distance; t is the time; D is the pipe diameter; JQ is the head loss caused by the quasi-steady friction which is determined by using the Hagen Poiseuille law and Colebrook–White formula (Adamkowski & Lewandowski 2006); Ju is the head loss due to unsteady friction.

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.

Original convolution-based unsteady friction models

Zielke (1968) developed pioneering work on the original convolution-based model by considering the convolution of previous fluid accelerations and a weighting function to express the unsteady friction items for laminar flow as follows (Zielke 1968):
formula
(3)
Zielke (1968) proposed the staggered-grid MOC scheme for the unsteady friction solution by the first-order approximation of convolution integral. Here, in the rectangular-grid MOC scheme solution, the unsteady friction can still be calculated by the first-order approximation of convolution integral:
formula
(4)
The effects of the past velocity can be described by using the Zielke weighting function :
formula
(5)
formula
(6)
formula
(7)
where j and k are multiples of the time step Δt; ν is the kinematic viscosity; τ is the dimensionless time; and coefficients of weighting function W(τ) {ni, i= 1, …, 5} = {− 26.3744, −70.8493, −135.0198, −218.9216, −322.5544} and {mi, i= 1, …, 6} = {0.282095, −1.25, 1.057855, 0.937500, 0.396696, −0.351563}.
Subsequently, for the turbulent flow, Vardy & Brown (1995, 1996, 2003, 2004b) developed the original turbulent weighting functions related to the Reynolds number and time (Vardy & Brown 1995, 1996, 2003, 2004b):
formula
(8)
where and , for smooth pipes and and for flows in rough pipes. The ratio ɛ/D is called the relative roughness.

Efficient simplified convolution-based model

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−8 ≤Δτ < ∞ for laminar flow and 10−9 ≤Δτ ≤ 10−1 for turbulent flow is introduced to simulate the unsteady friction of transient pipe flows since the possible range of Δτ is ranged from 10−7 to 10−5 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).

The TVB method of evaluating the unsteady shear stress τu is
formula
(9)

For laminar flow, an efficient implementation of the Zielke weighting function can be realized by the effective weighting functions with {ni, i= 1, …, 9} = {26.3744; 102; 102.5; 103; 104; 105; 106;107;108} and {mi, 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; 101.5; 102; 102.5; 103; 103.5; 104; 104.5; 105; 105.5; 106; 106.5; 107; 107.5; 108; 108.5; 109} 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}.

Brunone instantaneous acceleration-based unsteady friction model

Brunone et al. (1991) presented the unsteady friction part with the instantaneous local acceleration and instantaneous convective acceleration. Original Brunone formulation was improved as shown below (Bergant et al. 2001; Vítkovský et al. 2006):
formula
(10)
in which for V0, and for V< 0. The Brunone friction coefficient k can be predicted either empirically by trial and error or analytically using Vardy's shear decay coefficient C*:
formula
(11)
where C* depends on the flow regime. For laminar flow, C* = 0.00476; for turbulent flow, presented by Vardy & Brown (1995), are used in this work.

NUMERICAL SOLUTION BY USING SECOND-ORDER FVM

The first-order fixed-grid MOC scheme is widely used to solve the water hammer equations with the unsteady friction factors (Bergant et al. 2001; Vítkovský et al. 2006). In this work, since the effects of Courant number (Cr = a·Δtx) are considered, the space-line interpolation fixed-grid MOC scheme can give along C+ and C characteristic lines, as follows (Wylie et al. 1993):
formula
(12)
formula
(13)

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.

Alternatively, a robust second-order FVM GTS is presented for the water hammer solution considering the unsteady friction factors. The water hammer equations can be written in the form of a Riemann problem:
formula
(14)
where ; ; and .

Godunov-type scheme for water hammer equations

Figure 2 displays the FVM grid system, in which the pipeline is divided into N reaches by the fixed-grid length Δx. For the ith control volume, the integration of Equation (14) between control interfaces i − 1/2 and i+ 1/2 yields
formula
(15)
where Ui is the average value of u within [i − 1/2, i + 1/2]; the superscripts n and n + 1 indicate the t and t + Δt time levels, respectively; fi+1/2 and fi−1/2 are the mass and momentum fluxes at the control interfaces i − 1/2 and i+ 1/2, which is determined by solving a local Riemann problem at each cell interface (Zhao & Ghidaoui 2004; Toro 2009; Zhou et al. 2017, 2018a).
Figure 2

Grid system in the FVM Godunov scheme.

Figure 2

Grid system in the FVM Godunov scheme.

Applying Rankine–Hugoniot conditions where the eigenvalues and , the fluxes at i+ 1/2 for all internal nodes and for t[tn, tn+1] can be calculated by the following equation:
formula
(16)
in which = A; = average value of u to the left of interface i+ 1/2 at time n; and = average value of u to the right of interface i+ 1/2 at time n.

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

In the Godunov scheme, the Rankine–Hugoniot condition across each wave of speed gives the following relations:
formula
(17)
formula
(18)

For the upstream boundary (i = 0), coupling this Riemann invariant of Equation (17) with a head-flow boundary relation determines u1/2(t) = (H1/2, V1/2), . Similarly, for the downstream boundary (i = N), coupling the Riemann invariant of Equation (18) with a head-flow boundary relation determines uN+1/2(t) = (HN+1/2, VN+1/2), . Alternatively, virtual control volumes and I0 adjacent to I1, and virtual control volumes IN+1 and IN+2 adjacent to IN, are used to realize the direct solution of the Riemann problem for the boundary unknowns U1, U2, UN−1, and UN. Here, , and , where and are calculated by combining the Riemann invariant with a head-flow boundary relation at time n.

Incorporation of unsteady friction terms

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.

First step:
formula
(19)
Second step:
formula
(20)
where . Note: is updated by using , and is only related to the past variables rather than .
Last step:
formula
(21)
where . Note: is updated by using , and is still calculated by the past variables .

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

Figure 3

Procedure of the proposed algorithm.

Figure 3

Procedure of the proposed algorithm.

The Courant–Friedrichs–Lewy (CFL) criterion is satisfied in advance (because a is constant) by taking small enough for given Δx. However, the second-order Runge–Kutta procedure also has the stability constraint:
formula
(22)
Thus, the permissible time step for the source term (Δtmax,s) is given by:
formula
(23)
Finally, the maximum permissible time step is given by:
formula
(24)

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.

RESULTS AND DISCUSSION

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.

Validations of original convolution-based model and Brunone model in 2nd GTS

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.

Figure 4

Comparison of the measured data and the results calculated in the second-order GTS: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; and (e) Case 5.

Figure 4

Comparison of the measured data and the results calculated in the second-order GTS: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; and (e) Case 5.

Figure 5

Comparison of the measured data and the results calculated in the MOC scheme: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; and (e) Case 5.

Figure 5

Comparison of the measured data and the results calculated in the MOC scheme: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; and (e) Case 5.

The experimental Case 1 with initial velocity V0 = 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 (104) 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 R2) 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.

Table 2

Simulation error values for the experimental cases

Case no.ModelRMSE (m)NSE or R2
Case 1 MOC Quasi-steady 2.990 0.911 
Zielke 2.623 0.932 
Brunone 2.351 0.945 
2nd GTS Quasi-steady 2.990 0.911 
Zielke 2.574 0.934 
Brunone 2.435 0.941 
Case 3 MOC Quasi-steady 4.020 0.980 
Vardy–Brown 2.885 0.990 
Brunone 2.877 0.990 
2nd GTS Quasi-steady 4.020 0.980 
Vardy–Brown 2.883 0.990 
Brunone 2.832 0.990 
Case no.ModelRMSE (m)NSE or R2
Case 1 MOC Quasi-steady 2.990 0.911 
Zielke 2.623 0.932 
Brunone 2.351 0.945 
2nd GTS Quasi-steady 2.990 0.911 
Zielke 2.574 0.934 
Brunone 2.435 0.941 
Case 3 MOC Quasi-steady 4.020 0.980 
Vardy–Brown 2.885 0.990 
Brunone 2.877 0.990 
2nd GTS Quasi-steady 4.020 0.980 
Vardy–Brown 2.883 0.990 
Brunone 2.832 0.990 

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.

Evaluation of the second-order GTS TVB model

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−8 ≤ Δτ < ∞ for laminar flow and 10−9 ≤Δτ ≤ 10−1 for turbulent flow is preliminarily chosen for the current simulations, since the range of is around from 10−7 to 10−5 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

Comparison of the results calculated by the 2nd GTS Vardy–Brown model and the 2nd GTS TVB model: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; and (e) Case 5.

Figure 6

Comparison of the results calculated by the 2nd GTS Vardy–Brown model and the 2nd GTS TVB model: (a) Case 1; (b) Case 2; (c) Case 3; (d) Case 4; and (e) Case 5.

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.

Table 3

Computational time of the Brunone model and the TVB model in the schemes of the MOC and the second-order GTS at the same conditions

Model – Computation timeCase 1Case 2Case 3Case4Case5
2nd GTS Zielke or Vardy–Brown model-T1 (s) 6,375 6,680 6,044 7,450 8,163 
2nd GTS TVB model-T2 (s) 8.594 8.764 7.262 9.615 10.465 
Model – Computation timeCase 1Case 2Case 3Case4Case5
2nd GTS Zielke or Vardy–Brown model-T1 (s) 6,375 6,680 6,044 7,450 8,163 
2nd GTS TVB model-T2 (s) 8.594 8.764 7.262 9.615 10.465 

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.

Stability of the second-order GTS and the MOC scheme for unsteady friction models

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Δtx).

The MOC scheme

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Δtx < 1, the linear interpolation approach is used to calculate the variables at the MOC nodes at n time.

Figure 7

Influences of Courant number and grid number in the MOC TVB model: (a) N= 256; (b) N= 32.

Figure 7

Influences of Courant number and grid number in the MOC TVB model: (a) N= 256; (b) N= 32.

Figure 8

Influences of Courant number and grid number in the MOC Brunone model: (a) N= 256; (b) N= 32.

Figure 8

Influences of Courant number and grid number in the MOC Brunone model: (a) N= 256; (b) N= 32.

The second-order Godunov-type scheme

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.

Figure 9

Influences of Courant number and grid number in the second-order GTS TVB model: (a) N= 256; (b) N= 32.

Figure 9

Influences of Courant number and grid number in the second-order GTS TVB model: (a) N= 256; (b) N= 32.

Figure 10

Influences of Courant number and grid number in the second-order GTS Brunone model: (a) N= 256; (b) N= 32.

Figure 10

Influences of Courant number and grid number in the second-order GTS Brunone model: (a) N= 256; (b) N= 32.

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

Table 4

Error values when Cr = 0.5

ModelCrNRMSENSE or R2
MOC Brunone 0.5 32 4.589 0.792 
256 2.115 0.956 
Second-order GTS Brunone 0.5 32 2.268 0.949 
256 0.189 1.000 
ModelCrNRMSENSE or R2
MOC Brunone 0.5 32 4.589 0.792 
256 2.115 0.956 
Second-order GTS Brunone 0.5 32 2.268 0.949 
256 0.189 1.000 

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.

Table 5

Computational time of the MOC scheme and the second-order GTS for the given numerical accuracy

ModelMOC scheme
Second-order GTS
T9/T10
Case 1T9 (s)Case 1T10 (s)
Brunone model N= 256, Cr= 0.5 0.7344 N= 32, Cr= 0.5 0.0313 23.5 
N= 256, Cr= 0.1 3.5471 N= 32, Cr= 0.1 0.1412 25.1 
TVB model N= 256, Cr= 0.5 32.1412 N= 32, Cr= 0.5 0.3125 102.9 
N= 256, Cr= 0.1 162.1875 N= 32, Cr= 0.1 1.3594 119.3 
ModelMOC scheme
Second-order GTS
T9/T10
Case 1T9 (s)Case 1T10 (s)
Brunone model N= 256, Cr= 0.5 0.7344 N= 32, Cr= 0.5 0.0313 23.5 
N= 256, Cr= 0.1 3.5471 N= 32, Cr= 0.1 0.1412 25.1 
TVB model N= 256, Cr= 0.5 32.1412 N= 32, Cr= 0.5 0.3125 102.9 
N= 256, Cr= 0.1 162.1875 N= 32, Cr= 0.1 1.3594 119.3 

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.

Efficiency of the second-order GTS and the MOC scheme for unsteady friction models

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.

Table 6

Computational time of the Brunone model and the TVB model in the schemes of the MOC and the second-order GTS at the same conditions

Case 1N= 256
N= 32
Cr= 0.5Cr= 0.1Cr= 0.5Cr= 0.1
MOC Brunone model-T5 (s) 0.7344 3.5471 0.0156 0.0681 
MOC TVB model-T6 (s) 32.1412 162.1875 0.8125 3.1251 
2nd GTS Brunone model-T7 (s) 1.6412 8.3282 0.0313 0.1412 
2nd GTS TVB model-T8 (s) 12.0781 63.6702 0.3125 1.3594 
T6/T5 43.8 45.7 52.1 45.9 
T8/T7 7.4 7.6 10.0 9.6 
T7/T5 2.2 2.3 2.0 2.1 
T6/T8 2.7 2.5 2.6 2.3 
Case 1N= 256
N= 32
Cr= 0.5Cr= 0.1Cr= 0.5Cr= 0.1
MOC Brunone model-T5 (s) 0.7344 3.5471 0.0156 0.0681 
MOC TVB model-T6 (s) 32.1412 162.1875 0.8125 3.1251 
2nd GTS Brunone model-T7 (s) 1.6412 8.3282 0.0313 0.1412 
2nd GTS TVB model-T8 (s) 12.0781 63.6702 0.3125 1.3594 
T6/T5 43.8 45.7 52.1 45.9 
T8/T7 7.4 7.6 10.0 9.6 
T7/T5 2.2 2.3 2.0 2.1 
T6/T8 2.7 2.5 2.6 2.3 

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.

CONCLUSIONS

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.

ACKNOWLEDGEMENTS

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

DATA AVAILABILITY STATEMENT

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

REFERENCES

Adamkowski
A.
Lewandowski
M.
2006
Experimental examination of unsteady friction models for transient pipe flow simulation
.
J. Fluids Eng.
128
(
6
),
1351
1363
.
https://doi.org/10.1115/1.2354521
Bergant
A.
Simpson
A. R.
Vitkovský
J. P.
2001
Developments in unsteady pipe flow friction modelling
.
J. Hydraul. Res.
39
(
3
),
249
257
.
https://doi.org/10.1080/00221680109499828
Brunone
B.
Golia
U. M.
Greco
M.
1991
Some remarks on the momentum equations for fast transients
. In
Hydraulic Transients with Column Separation (9th and Last Round Table of the IAHR Group)
.
IAHR
,
Valencia
,
Spain
, pp.
201
209
.
Chaudhry
M. H.
Hussaini
M. Y.
1985
Second-order accurate explicit finite-difference schemes for waterhammer analysis
.
J. Fluids Eng.
107
(
4
),
523
529
.
https://doi.org/10.1115/1.3242524
Ghidaoui
M. S.
Karney
B. W.
1994
Equivalent differential equations in fixed-grid characteristics method
.
J. Hydraul. Eng.
120
(
10
),
1159
1175
.
https://doi.org/10.1061/(ASCE)0733-9429(1994)120:10(1159)
Ghidaoui
M. S.
Mansour
S.
2002
Efficient treatment of the Vardy-Brown unsteady shear in pipe transients
.
J. Hydraul. Eng.
128
(
1
),
102
112
.
https://doi.org/10.1061/(ASCE)0733-9429(2002)128:1(102)
Ghidaoui
M. S.
Karney
B. W.
McInnis
D. A.
1998
Energy estimates for discretization errors in waterhammer problems
.
J. Hydraul. Eng.
124
(
4
),
384
393
.
https://doi.org/10.1061/(ASCE)0733-9429(1998)124:4(384)
Guinot
V.
2000
Riemann solvers for water hammer simulations by Godunov method
.
Int. J. Numer. Methods Eng.
49
(
7
),
851
870
.
https://doi.org/10.1002/1097-0207(20001110)49:7&lt;851::AID-NME978 > 3.0.CO;2-%23
Kagawa
T.
Lee
I.
Kitagawa
A.
Takenaka
T.
1983
High speed and accurate computing method of frequency-dependent friction in laminar pipe flow for characteristic method
.
Trans. Jpn. Soc. Mech. Eng., Ser. A
49
(
447
),
2638
2644
.
Karney
B. W.
Ghidaoui
M. S.
1997
Flexible discretization algorithm for fixed-grid MOC in pipelines
.
J. Hydraul. Eng.
123
(
11
),
1004
1011
.
https://doi.org/10.1061/(ASCE)0733-9429(1997)123:11(1004)
Martins
N. M. C.
Carrico
N. J. G.
Ramos
H. M.
Covas
D. I. C.
2014
Velocity-distribution in pressurized pipe flow using CFD: accuracy and mesh analysis
.
Comput. Fluids
105
,
218
230
.
https://doi.org/10.1016/j.compfluid.2014.09.031
Schohl
G. A.
1993
Improved approximate method for simulating frequency-dependent friction in transient laminar flow
.
J. Fluids Eng.
115
(
3
),
420
424
.
https://doi.org/10.1115/1.2910155
Suzuki
K.
Taketomi
T.
Sato
S.
1991
Improving Zielke's method of simulating frequency-dependent friction in laminar liquid pipe flow
.
J. Fluids Eng.
113
,
569
573
.
https://doi.org/10.1115/1.2926516
Toro
E. F.
2009
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
, 3rd edn.
Springer-Verlag
,
Berlin, Heidelberg
.
Trikha
A. K.
1975
An efficient method for simulating frequency-dependent friction in transient liquid flow
.
J. Fluids Eng.
97
,
97
105
.
https://doi.org/10.1115/1.3447224
Urbanowicz
K.
Zarzycki
Z.
2012
New efficient approximation of weighting functions for simulations of unsteady friction losses in liquid pipe flow
.
J. Theor. Appl. Mech.
50
(
2
),
487
508
.
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.
1996
On turbulent, unsteady, smooth-pipe flow
. In
Proc. Int. Conf. on Pressure Surges and Fluid Transients
.
BHR Group
,
Harrogate
,
UK
, pp.
289
311
.
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.
2004a
Efficient approximation of unsteady friction weighting functions
.
J. Hydraul. Eng.
130
(
11
),
1097
1107
.
https://doi.org/10.1061/(ASCE)0733-9429(2004)130:11(1097)
Vardy
A. E.
Brown
J. M. B.
2004b
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
Vardy
A. E.
Brown
J. M. B.
2007
Approximation of turbulent wall shear stresses in highly transient pipe flows
.
J. Hydraul. Eng.
133
(
11
),
1219
1228
.
https://doi.org/10.1061/(ASCE)0733-9429(2007)133:11(1219)
Vitkovsky
J.
Stephens
M.
Bergant
A.
Martin
L.
Simpson
A.
2004
Efficient and accurate calculation of Zielke and Vardy-Brown unsteady friction in pipe transients
. In
Proc., 9th Int. Conf. on Pressure Surges
.
BHR Group
,
Cranfield
,
UK
, pp.
405
419
.
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
.
https://doi.org/10.1061/(ASCE)0733-9429(2006)132:7(696)
Wang
H.
Zhou
L.
Liu
D. Y.
Karney
B.
Wang
P.
Xia
L.
Ma
J. J.
Xu
C.
2016
CFD approach for column separation in water pipelines
.
J. Hydraul. Eng.
142
(
10
),
11
.
https://doi.org/10.1061/(ASCE)HY.1943-7900.0001171
Wylie
E. B.
Streeter
V. L.
1970
Network system transient calculations by implicit method
. In
45th Annual Meeting of the Society of Petroleum Engineers of AIME, Paper Number 2963
.
Houston
.
Wylie
E. B.
Streeter
V. L.
Suo
L.
1993
Fluid Transients in Systems
.
Prentice Hall
,
Englewood Cliffs, NJ
.
Zhao
M.
Ghidaoui
M.
2004
Godunov-type solutions for water hammer flows
.
J. Hydraul. Eng.
130
(
4
),
341
348
.
https://doi.org/10.1061/(ASCE)0733-9429(2004)130:4(341)
Zhou
L.
Wang
H.
Liu
D. Y.
Ma
J. J.
Wang
P.
Xia
L.
2017
A second-order finite volume method for pipe flow with water column separation
.
J. Hydro-Environ. Res.
17
,
47
55
.
https://doi.org/10.1016/j.jher.2016.11.004
Zhou
L.
Wang
H.
Bergant
A.
Tijsseling
A. S.
Liu
D.
Guo
S.
2018a
Godunov-type solutions with discrete gas cavity model for transient cavitating pipe flow
.
J. Hydraul. Eng.
144
(
5
),
04018017
.
https://doi.org/10.1061/(ASCE)HY.1943-7900.0001463
Zhou
L.
Wang
H.
Karney
B.
Liu
D.
Wang
P.
Guo
S.
2018b
Dynamic behavior of entrapped air pocket in a water filling pipeline
.
J. Hydraul. Eng.
144
(
8
),
04018045
.
https://doi.org/10.1061/(ASCE)HY.1943-7900.0001491
Zielke
W.
1968
Frequency-dependent friction in transient pipe flow
.
J. Basic Eng.
90
(
1
),
109
115
.
https://doi.org/10.1115/1.3605049
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/).