## 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

*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}*i*th 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*_{0}initial steady pipe flow velocity

*W*(*t*)a weighting function

*x*distance along pipeline

*Y*(_{ai}*t*)the weighting function

- Δ
*t*time increment

- Δ
*τ*4

*v*Δ*t*/*D*^{2} - Δ
*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

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

Initially, the pipe flow reaches a steady velocity *V*_{0} 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*= 127.47 m,

_{r}*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.

Case no. . | V_{0} (m/s)
. | Re
. | Flow regime . | Data source . |
---|---|---|---|---|

1 | 0.1 | 1,870 | Laminar flow | Bergant et al. (2001) |

2 | 0.2 | 3,750 | Turbulent flow | |

3 | 0.3 | 5,600 | Turbulent flow | |

4 | 0.631 | 10,634 | Turbulent flow | Adamkowski & Lewandowski (2006) |

5 | 0.94 | 15,843 | Turbulent flow |

Case no. . | V_{0} (m/s)
. | Re
. | Flow regime . | Data source . |
---|---|---|---|---|

1 | 0.1 | 1,870 | Laminar flow | Bergant et al. (2001) |

2 | 0.2 | 3,750 | Turbulent flow | |

3 | 0.3 | 5,600 | Turbulent flow | |

4 | 0.631 | 10,634 | Turbulent flow | Adamkowski & Lewandowski (2006) |

5 | 0.94 | 15,843 | Turbulent flow |

## WATER HAMMER EQUATIONS INVOLVING UNSTEADY FRICTION

*et al.*1993):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;

*J*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);

_{Q}*J*is the head loss due to unsteady friction.

_{u}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

*j*and

*k*are multiples of the time step Δ

*t*;

*ν*is the kinematic viscosity;

*τ*is the dimensionless time; and coefficients of weighting function

*W*(

*τ*) {

*n*,

_{i}*i*

*=*1, …, 5} = {− 26.3744, −70.8493, −135.0198, −218.9216, −322.5544} and {

*m*,

_{i}*i*

*=*1, …, 6} = {0.282095, −1.25, 1.057855, 0.937500, 0.396696, −0.351563}.

*ɛ*/

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

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

^{2}; 10

^{2.5}; 10

^{3}; 10

^{4}; 10

^{5}; 10

^{6};10

^{7};10

^{8}} 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^{2}; 10^{2.5}; 10^{3}; 10^{3.5}; 10^{4}; 10^{4.5}; 10^{5}; 10^{5.5}; 10^{6}; 10^{6.5}; 10^{7}; 10^{7.5}; 10^{8}; 10^{8.5}; 10^{9}} 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

*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):in which for

*V*

*≥*

*0*, 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**: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

*et al.*2001; Vítkovský

*et al.*2006). In this work, since the effects of Courant number (

*Cr*=

*a*·Δ

*t*/Δ

*x*) are considered, the space-line interpolation fixed-grid MOC scheme can give along

*C*

^{+}and

*C*

^{−}characteristic lines, as follows (Wylie

*et al.*1993):

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.

### Godunov-type scheme for water hammer equations

*N*reaches by the fixed-grid length Δ

*x*. For the

*i*th control volume, the integration of Equation (14) between control interfaces

*i*− 1/2 and

*i*

*+*1/2 yieldswhere

**U**

*is the average value of*

_{i}**u**within [

*i*− 1/2,

*i*+ 1/2]; the superscripts

*n*and

*n*+ 1 indicate the

*t*and

*t*+ Δ

*t*time levels, respectively;

**f**

_{i}_{+1/2}and

**f**

_{i−}_{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).

*i*

*+*1/2 for all internal nodes and for

*t*[

*t*,

^{n}*t*

^{n+}^{1}] can be calculated by the following equation: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).

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*

_{0}adjacent to

*I*

_{1}, and virtual control volumes

*I*

_{N+}_{1}and

*I*

_{N+}_{2}adjacent to

*I*, are used to realize the direct solution of the Riemann problem for the boundary unknowns

_{N}**U**

_{1},

**U**

_{2},

**U**

_{N−}_{1}, and

**U**

*. Here, , and , where and are calculated by combining the Riemann invariant with a head-flow boundary relation at time*

_{N}*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.

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.

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

The experimental Case 1 with initial velocity *V*_{0} = 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^{4}) 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*^{2}) 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.

Case no. . | Model . | RMSE (m) . | NSE or R^{2}
. | |
---|---|---|---|---|

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. . | Model . | RMSE (m) . | NSE or R^{2}
. | |
---|---|---|---|---|

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

Model – Computation time . | Case 1 . | Case 2 . | Case 3 . | Case4 . | Case5 . |
---|---|---|---|---|---|

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 time . | Case 1 . | Case 2 . | Case 3 . | Case4 . | Case5 . |
---|---|---|---|---|---|

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*Δ*t*/Δ*x*).

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

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

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

Model . | Cr
. | N
. | RMSE . | NSE or R^{2}
. |
---|---|---|---|---|

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 |

Model . | Cr
. | N
. | RMSE . | NSE or R^{2}
. |
---|---|---|---|---|

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.

Model . | MOC scheme . | Second-order GTS . | T9/T10 . | ||
---|---|---|---|---|---|

Case 1 . | T9 (s) . | Case 1 . | T10 (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 |

Model . | MOC scheme . | Second-order GTS . | T9/T10 . | ||
---|---|---|---|---|---|

Case 1 . | T9 (s) . | Case 1 . | T10 (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.

Case 1 . | N= 256. | N= 32. | ||
---|---|---|---|---|

Cr= 0.5
. | Cr= 0.1
. | Cr= 0.5
. | Cr= 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 1 . | N= 256. | N= 32. | ||
---|---|---|---|---|

Cr= 0.5
. | Cr= 0.1
. | Cr= 0.5
. | Cr= 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:

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.

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.

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.

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.

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.