Abstract

The proposed study investigated the applicability of the finite volume method (FVM) based on the Godunov scheme to transient water hammer with shock front simulation, in which intermediate fluxes were computed using either first-order or second-order Riemann solvers. Finite volume (FV) schemes are known to conserve mass and momentum and produce the efficient and accurate realization of shock waves. The second-order solution of the Godunov scheme requires an efficient slope or a flux limiter for error minimization and time optimization. The study examined a range of limiters and found that the MINMOD limiter is the best for modeling water hammer in terms of computational time and accuracy. The first- and second-order FVMs were compared with the method of characteristics (MOCs) and experimental water hammer measurements available in the literature. Both the FV methods accurately predicted the numerical and experimental results. Parallelization of the second-order FVM reduced the computational time similar to that of first-order. Thus, the study presented a faster and more accurate FVM which is comparable to that of MOC in terms of computational time and precision, therefore it is a good substitute for the MOC. The proposed study also investigated the implementation of a more complex convolution-based unsteady friction model in the FVM to capture real pressure dissipation. The comparison with experimental data proved that the first-order FV scheme with the convolution integral method is highly accurate for computing unsteady friction for sudden valve closures.

HIGHLIGHTS

  • Error estimation of different flux limiters and finding the best flux limiter for water hammer modeling.

  • Establishing that first-order FVM is acceptable for water hammer modeling.

  • Using parallelization routines in FVM.

  • Implementing convolution-based unsteady friction method.

  • Methodology along with complete derivations is provided for better understanding of FVM water hammer modeling.

INTRODUCTION

Water is one of the principal requirements of human life, and water pipeline distribution systems are one of the most important aspects of civic operations. Often, pipeline distribution systems suffer from hydraulic transients caused by sudden valve closures, sudden valve openings and pump tripping. These transients, if not properly contained, can cause devastation to property and lives. Understanding hydraulic transients and associated modeling of water hammer pressures have been prevalent topics of research for several decades. Numerous numerical modeling methods have been adopted, among which the method of characteristics (MOCs) has become the most popular due to its simplicity, understandability, accuracy and computational efficiency. The MOC also offers some limitations.

The convective term of the governing equation is neglected in the MOC (Guinot 2000). Furthermore, other basic disadvantages of the MOC scheme are as follows: the MOC scheme poses risk of artificially interpolating the wave celerity and threats to the solution with infusing artificial damping in cases of Courant number less than one (Ghidaoui & Karney 1994; Zhang et al. 2018a, 2018b; Li et al. 2019). According to Zhang et al. (2018a, 2018b), results are slightly affected due to neglecting the convective acceleration term in MOC, which can be significant in multidimensional complex pipeline networks. In addition, according to Afshar & Rohani (2008), the explicit MOC is cumbersome if more than one device is located between two pipes. The uniform celerity approximation (Guinot 2000) for the MOC scheme is an ideal case which is unquestionable when pipe and fluid characteristics are constant, but it may not be conserved for some specific types of problems such as cavitation resulting from the water hammer and in the cases where the impedance varies along the pipe, e.g., for different wall thickness. As a result, it is better to search for other numerical schemes that can solve the governing equations effectively and efficiently apart from the extensively used MOC schemes. Therefore, researchers have started exploring alternative methods for solving discontinuity-based hydraulic transients, especially water hammer-type problems. Some of the alternate methods include FVM, Lattice Boltzmann method (LBM), method of lines (MOLs) and time-splitting method. Budinsky (2016) modeled one-dimensional water hammer using the LBM on an adaptive grid, which is independent of the Courant number. The LBM simulation was validated using the MOC for water hammer resulting from pump tripping. The LBM is computationally fast and equally accurate to that of MOC for low-frequency water hammer waves; however, LBM becomes unstable and inaccurate for high-frequency water hammer waves. Chen et al. (2015) used the MOL to approximate the water hammer Partial Differential Equations (PDEs) by a system of Ordinary Differential Equations (ODEs). In this way, Chen et al. (2015) used a system of ODEs as an optimal control to mitigate water hammer. The MOL has stability and accuracy issues. Niroomandi et al. (2012) investigated the time-splitting projection method (TSPM) to model transient flow in elastic pipes using the quasi-steady model of static friction and Brunone's model of unsteady friction and evaluated its accuracy by L2 relative error norm thereby proving the applicability of TSPM for compressible pipe-type problems. The TSPM is mainly applicable for incompressible flows and poses challenges while applying for compressible flows. In view of the limitations of above-mentioned alternate methods, the FVM is applied and analyzed in this paper for its accuracy and speed.

In this direction, pioneering research was carried out by Guinot (2000), in which a simple pipeline reservoir system was solved using the finite volume (FV)-based Godunov scheme without frictional head losses. The FV approach is well known in modeling gas dynamics problems because of its sharp resolution producing capability, conservation of mass and momentum, and applicability to non-uniform unstructured grids. Hwang & Chung (2002) utilized the FVM for water hammer without neglecting the advective term and using the conservative form with density as an unknown instead of the pressure head. Szydlowski (2002) used a second-order accurate FV scheme for one-dimensional transient flow in a single pipeline using a Riemann solver with the Roe linearization technique to analyze the water hammer pressure. Zhao & Ghidaoui (2004) further investigated the application of FV technique based on mass and momentum conservation principles for water hammer modeling. They considered an MOC-type system of two pipes connected in series, valve, valve junction and a reservoir for the study and solved the numerical flux by exact solution of the Riemann problem using both first- and second-order Godunov schemes. Zhao & Ghidaoui (2004) concluded that the first-order Godunov scheme showed high numerical dissipation with reduced accuracy, but using the second-order Godunov scheme with the MINMOD limiter increased the accuracy and efficiency of the model to a significant extent. Yazdi et al. (2007) used a similar approach of second-order FV Godunov to develop a one-dimensional coupled model based on Riemann solution with the analogous type of arrangement and compared the results with analytical solutions. They concluded that the second-order scheme is in good correspondence with the analytical method at Courant number () equal to one, but for a lesser Courant number it shows traces of numerical dissipation. Later, Leon et al. (2008) continued their work on the problem and introduced a FV shock capturing scheme for both one- and two-phase flows. They compared the results between the FV model and the available MOC for single-phase flow and concluded that the MOC is suitable for , but for the efficiency and accuracy of the shock capturing FV scheme proved to be significantly better. Daude et al. (2018) studied vaporous cavitation-induced column separation using one-dimensional FV technique and compared the simulations with the published experimental data.

The aims of this paper are (i) to provide detailed formulations of first- and second-order FVM, (ii) to find the accurate flux limiter for second-order FVM, (iii) to check whether first-order FVM is adequate for water hammer, (iv) to investigate parallel computing implementation in the FVM modeling of water hammer for improving computational speed and (v) to implement a more accurate weighting function-based unsteady friction model in the FVM. The organization of the paper is divided into the following sections. The first section explains the methodology of the finite volume method (FVM) for solving the governing equations, first-order solution, second-order solution, examination of the range of flux limiters, initial and boundary conditions, and parameters for study, validation and discussion. The methodology also covered implementation details of a more complex convolution-based unsteady friction model. The Results and discussion section demonstrates successful validation of the proposed first-order and second-order FV models by comparing them with the MOCs on a hypothetical pipeline for water hammer. The final part of the Results and discussion demonstrates the successful validation of the proposed FV models by testing with experimental datasets available in the literature as well as the accuracy of convolution-based unsteady friction model in dissipating pressure. Lastly, a summary of the important findings are presented.

METHODOLOGY

Governing equations

According to Guinot (2000), the governing equation for water hammer can be expressed in the conservative form as
formula
(1)
formula
where is the vector of flow variables, is the flux, is the cross-sectional area of the pipe, is the friction factor based on the pipe roughness and the fluid viscosity, P is the pressure, is the mass discharge, u is the fluid velocity, is the mass of fluid per unit length of pipe and = . Variables and Qm are conservative and those can be converted to nonconservative pressure and density form as required using the following equation:
formula
(2)
where , is the reference fluid density and c is the pressure wave celerity.
Now expressing the governing equation in its characteristic form
formula
(3)
where
formula

Equation (3) is a hyperbolic system of conservation laws, in which B is the Jacobian matrix of F in Equation (1) with respect to U. If U is derivable, then this characteristic form will be continuous with respect to space.

Guinot (2000) applied the well-known Godunov scheme to calculate the flux at the interface as follows:
formula
(4)
where i is the number of spatial cells, n is the number of temporal cell, is the time step, is the grid size, and and are the fluxes at cell interfaces. Guinot (2000) considered only the hyperbolic term of the equation and solved F(U) for in between to cell. The source term S was incorporated into the equation as shown in Equation (4).
Riemann's concept of wave propagation was applied to find the fluxes, and , at the cell interfaces. The eigenvalues of the Jacobian matrix represent the celerities of the pressure waves. A disturbance travels along the given characteristic lines throughout the domain as shown in Figure 1. Therefore, the solution of the Riemann problem consists of two contact discontinuities with travel velocities and . The solution of the intermediate state is expressed as
formula
(5)
where and represent the two contact discontinuities separating those three regions (Guinot 2000).
Figure 1

Solution of the Riemann problem in the physical space (top) and in the phase space (bottom).

Figure 1

Solution of the Riemann problem in the physical space (top) and in the phase space (bottom).

In Figure 1, the represents the intermediate state of the two solutions at left and right. Guinot (2000) suggested solving the Riemann problem at the interface by using a number of exact or approximate solvers. Further, he suggested the application of the Rankine Hugoniot condition for capturing the shock-wave characteristics due to the jump condition.

The application of the FVM needs descritization of the computational domain in space and time. Figure 2 shows the typical discretization adopted in the FVM, in which computational cells are used instead of computational nodes as in the finite difference method. The FVM does not store values at grid points, whereas it stores average values over the intervals. Therefore, variables of interest are average values over the computational cell which are closely linked to the integral form of the governing equations. The variables (vectors) are subscripted and superscripted. The superscript n and subscript i refer to the nth time level and the ith index of the cell, respectively. The subscripts and refer to interfaces between i and , and i and , respectively. Superscript refers to the time average of a variable between n and time levels. For example, the vector denotes the average value of U between n and time levels and at the interface between i and cells. For the one-dimensional problem such as here, the length of the cell is considered as .

Figure 2

Discretization of computational domain in the FVM.

Figure 2

Discretization of computational domain in the FVM.

First-order Godunov method

In the first-order framework, the average value of over each cell i at a time level n is assumed to be known and the objective is to find the average value at the next time level. The actual distribution of U over the cell is not known and it does not need to be continuous within the cell. A first step consists in ‘guessing’ the distribution of the variable U within the cell i at the time level n. This step is also called the reconstruction step; in the original Godunov scheme, the profile in a given cell is assumed to be equal to the average of U over the cell.

Discretization of first-order equation:

According to Guinot (2000), mass and momentum equations can be written as follows:
formula
(6)
formula
(7)
formula
(8)
and where the source term .

Here, the friction term can be further written as where .

Where is the density of water, f is the friction factor, A is the area and D is the diameter of the pipe.

Characteristic lines:

A negative characteristic line can be given as follows:
formula
(9)
So,
formula
(10)
Similarly, the positive characteristic equation can be written as follows:
formula
(11)
So,
formula
(12)
Now
formula
(13)
So
formula
(14)
Let
formula
formula
(15)
formula
(16)
formula
(17)
formula
(18)
and
formula
(19)
formula
(20)
formula
(21)
So,
formula
(22)
and
formula
(23)
Solving we get
formula
(24)
formula
(25)
So, the equation
formula
(26)
formula
(27)
Godunov scheme:
formula
(28)
in which
formula
Integrating we get
formula
(29)
Solving the value at a new time step can be expressed as follows:
formula
(30)
And,
formula
(31)
Similarly,
formula
(32)
formula
(33)
Assuming , we get the final expression for pressure and discharge as follows:
formula
(34)
Similarly, the next equation can be solved and simplified as follows:
formula
(35)
Stability condition for source term:
In Equation (3), the source term S can be descritized as follows:
formula
(36)

Here, is the flow variable at a new time level without the source term, whereas is the flow variable incorporated with the source term.

For explicit problem, it is subjected to a stability constraint; therefore, the problem should satisfy
formula
Using this in Equation (4)
formula
So
formula
Hence,
formula

So, maximum permissible time step, If the grids are not uniform, then and for uniform grid .

Second-order Godunov method

The second-order solution only differs in the solution technique at the interfaces. The first-order Godunov method assumes that each cell was bearing a single constant value and the solution was obtained accordingly as shown in Figure 3. But, for the case of second-order Godunov scheme, the reconstruction of cell values was adopted. For the MUSCL-Hancock (Monotonic Upwind Scheme for Conservation Laws) scheme, the reconstruction is taken in the form of a straight line with some specific slope. Here, the cell value is assumed as a linear value instead of a constant value as in the case of first order, and the solution is carried out to find the flux through the cell interface, i + 1/2. To keep the cell values within the cell limit, a suitable limiter function is required as follows:

Figure 3

(a) Reconstruction for first-order Godunov scheme using the constant value for each cell and (b) reconstruction for second-order Godunov scheme using the piecewise linear function for each cell.

Figure 3

(a) Reconstruction for first-order Godunov scheme using the constant value for each cell and (b) reconstruction for second-order Godunov scheme using the piecewise linear function for each cell.

The following section presents a brief description of different types of limiter functions.
formula
(37)
Superbee:
formula
formula
formula
MC (monotonized central difference limiter):
formula
Van Leer:
formula
formula
Van Albada:
formula
formula
Double MINMOD:
formula

MINMOD:

In the proposed study, the MINMOD limiter has been used since this has given the best result as shown in the comparative study presented in the section ‘Results and discussion’.
formula
formula
So,
formula
(38)
formula
(39)
Second stage, i.e. corrector stage is given as follows:
formula
(40)
formula
(41)

Those two values are used as the left and right cell value of each interface while solving the Godunov flux.

Error estimation of slope limiter:

The above-mentioned different flux limiters may contain numerical errors and the most accurate flux limiter has to be ascertained for water hammer modeling. For the estimation of error extent of each flux limiter method, two statistical loss functions have been used in this study. Those are known as L1 and L2 norms. L1 loss function is also known as least absolute deviations (LADs), and L2 loss function is known as least square (LS) error.

The L1 function is extensively used for error reduction in terms of summation of absolute difference of the true value and the estimated value.
formula
(42)
In addition, the relative error can be explained as follows:
formula
(43)
whereas the L2 function is for the reduction of error in terms of the summation of squared difference of the true value and the estimated value.
formula
(44)
And relative error can be explained as follows:
formula
(45)
where X represents any numerical parameter and its subscript denotes the nature of the parameter.

The derivations presented so far are summarized here. Assuming a constant speed of sound and steady-state friction along the pipe as the source term, the conservative expressions for pressure and velocities at the interfaces were developed for the Godunov first-order FVM and the second-order MUSCL FVM method. In the MUSCL scheme, the Godunov method extended with piecewise linear approximations of each cell, which results in a central difference scheme that is second-order accurate in space. The computations are carried out in fully explicit approach which offers reduced complexity and reduced computational time without an iterative solver. Further, the details of different flux limiters of MUSCL and the method for error estimation of the flux limiters were described.

Initial and boundary conditions

Initial and boundary conditions need to be specified for flow simulations. In this study, the initial condition is taken as a steady state with a constant discharge flowing through the pipeline and a constant inlet pressure. Figure 4 shows the schematic of the model explained here. A constant inlet pressure at the upstream and the valve boundary at the downstream were used. In addition, negative and positive characteristics have been used at the upstream and downstream boundaries, respectively.

Figure 4

Model schematic for a simple pipeline with a large reservoir and the valve system.

Figure 4

Model schematic for a simple pipeline with a large reservoir and the valve system.

Test model description

Reservoir boundary condition at upstream:

Water level in the reservoir is assumed to be constant during the unsteady flow. The negative characteristic equation at the upstream reservoir can be written as
formula
(46)
Static pressure after entrance losses at the upstream reservoir is . At left boundary . Hence, incorporating pressure boundary condition and simplifying for discharge, Equation (46) becomes
formula
(47)
Valve boundary condition at downstream:
The valve is used to control water flow by partially or fully opening and closing. The valve regulates pressure drop of fluid flow. The implementation of the valve boundary condition is given as follows (Chaudhry 2014):
formula
(48)
In Equation (48), is flow coefficient, steady-state flowrate, and steady-state pressure loss across the valve. For unsteady flow, the equation can be written as follows, where subscript ‘p’ represents unsteady terms:
formula
(49)
Using Equations (48) and (49) and applying some modifications, we get
formula
(50)
where,
formula

For full valve opening = 1, and this value of s is decreased from 1 to 0 during the valve closure.

Flow chart

A flow chart of the second-order Godunov method for water hammer modeling is presented in Figure 5 for clarity. This method follows the MUSCL-Hancock approach which uses a Roe-type Riemann invariant solver. The MUSCL-Hancock scheme provides highly accurate numerical solutions which are shock free and without discontinuities (van Leer 1979). In the following flowchart, the second-order solution is shown only with the MINMOD limiter, but it is replaceable by any other limiter as presented in the previous section.

Figure 5

Flowchart for the FV scheme using second-order Godonov method with the MINMOD limiter.

Figure 5

Flowchart for the FV scheme using second-order Godonov method with the MINMOD limiter.

Unsteady friction computaton using the convolution theory

It is well known that steady-state friction is inadequate for simulating transient flows except for the first peak. To model unsteady friction, instantaneous acceleration-based (IAB) and convolution-based methods are available. The IAB methods require the selection of correct decay coefficients following the methods given in Pezzinga (2000) and Vardy & Brown (1996). The IAB model was already implemented in the FVM by Yazdi et al. (2007) and Seck et al. (2017). In this paper, to account accurate dissipation of pressure, wall friction during the transient state was modeled using a weighting function model as proposed by Vardy & Brown (2007) based on two-region viscosity distribution.

The wall shear stress is a sum of steady and unsteady friction terms as given as follows:
formula
(51)
where is the total wall shear stress, is the steady-state wall shear stress and is the unsteady-state wall shear stress.
The steady and unsteady wall shear stresses are modeled as
formula
(52)
formula
(53)
where is the dynamic viscosity of fluid, is the cross-sectional average velocity, is the diameter of pipe, w is the weighting function, t is the current instant of time and is the integration variable.
The weighting function can be calculated in the form of numerical integration of convolution integral. The numerical integration is carried by setting up intervals of length and evaluating the weighting function and local accelerations at the midpoints of the each integral. For convenience, we will use the same discretization for t and , i.e. so that .
formula
(54)
Because and are evenly spaced = , where k and j refer to current and previous time indexes.
formula
(55)
According to Vardy & Brown (2007), the weighting function can be determined as follows:
formula
(56)
where and are given by Vardy & Brown (2007) for Reynold number in the range of . These coefficients depend on the frozen viscosity Reynolds number and characteristic roughness. Empirical equations to determine the cofficients are given as follows:
formula
(57)
formula
(58)

In which, is the characteristic roughness ratio, is the diameter of the pipe and is the frozen viscosity Reynold number.

RESULTS AND DISCUSSION

The proposed FVM numerical model for simulating unsteady pipeflows was validated using the MOC and experimental data available in the literature. In these validations, downstream valve closures as well as upstream valve opening scenarios were simulated using the FVM.

Comparative study with the MOC

For this validation, a simple pipeline with a reservoir and a valve model system as shown in Figure 3 was considered. Pipe characteristics such as internal diameter, D and length, L were taken as 0.063 and 10,000 m, respectively. Initial steady-state discharge, , i.e. prior to valve closure was taken as 0.0021 m3/s and constant upstream reservoir head, was taken as 102 m. Steel pipe data were given as wave celerity 1,405 m/s, wall thickness 0.005 mm, roughness height 2.54 × 10−5 m and friction factor 0.0235. The water and pipeline rigidity are assumed for the solution. The domain was specially discretized into 100 segments. Total time for simulation was taken for 100 s after the valve closure. A constant time step of 0.01 s was chosen for the computations. A valve closure time of 1 s was used for simulating the water hammer effect. The valve characteristic s was varied linearly from 1.0 to 0.0 during the valve closure simulation. The valve started closing at t = 2 s and completely closed at t = 3 s. Computational specifications of this simulation are also given in Table 1 for better understanding.

Table 1

Detail of numerical parameters for the validation study

ParameterNotationUnitValue
Pipe length  10,000 
Pipe diameter  cm 6.3 
Initial discharge  m3/s 0.0021 
Initial head  102 
Wave celerity  m/s 1,405 
Fluid density  kg/m3 1,000 
Total time of computation  100 
Time step size  0.01 
Spatial step size  100 
Valve closure time  
Friction factor   0.0235 
Courant number   1.0 
ParameterNotationUnitValue
Pipe length  10,000 
Pipe diameter  cm 6.3 
Initial discharge  m3/s 0.0021 
Initial head  102 
Wave celerity  m/s 1,405 
Fluid density  kg/m3 1,000 
Total time of computation  100 
Time step size  0.01 
Spatial step size  100 
Valve closure time  
Friction factor   0.0235 
Courant number   1.0 

Investigation of flux limiters for water hammer

The second-order solution of the Godunov scheme requires an efficient slope or flux limiter for error minimization and time optimization. For this purpose, various limiters were employed and corresponding results were compared for computational time and for erroneous deviation. Time series length taken in this study was calculated by , which is the time required for one complete wave travel cycle. The results of the study were presented in Table 2.

Table 2

Computed error norms of different limiters

Limiter usedNorm L1Norm L2Relative error norm L1Relative error norm L2Computational time (s)
Superbee 162.64 16.16 3.6  0.91 
MC 156.10 14.65 3.5  0.85 
Van Leer 154.71 14.34 3.4  0.86 
Van Albada 152.51 14.09 3.3  0.87 
Double MINMOD 155.59 14.63 3.4  0.88 
MINMOD 147.68 13.08 3.2  0.89 
Limiter usedNorm L1Norm L2Relative error norm L1Relative error norm L2Computational time (s)
Superbee 162.64 16.16 3.6  0.91 
MC 156.10 14.65 3.5  0.85 
Van Leer 154.71 14.34 3.4  0.86 
Van Albada 152.51 14.09 3.3  0.87 
Double MINMOD 155.59 14.63 3.4  0.88 
MINMOD 147.68 13.08 3.2  0.89 

Bold values refer to the minimum values.

From Table 2, it is evident that although the MC limiter is showing the least computational time, the MINMOD limiter is showing the least error among the limiters studied. The ascending order of computational time for the limiters can be given as follows: Time of MC < Time of Van Leer < Time of Van Albada < Time of Double MINMOD < Time of MINMOD < Time of Superbee, whereas increasing order of errors (Norm L1 and Norm L2) is given as: error of MINMOD < error of Van Albada < error of Van Leer < error of Double MINMOD < error of MC < error of Superbee. From the analysis of computational time and the accuracy point of view, the MINMOD limiter is found to be the best. Therefore, it can be concluded that the application of MINMOD will be beneficial and the same is adopted for continuing further research on the problem.

Comparison with the MOC

Figure 6(a) and 6(b) show the validation of first- and second-order FVM results with that of MOC. The MOC scheme presented in Chaudhry (2014) and Wylie & Streeter (1993) was used in these simulations. In Figure 6(b), the second-order computation was done using the MINMOD limiter since it was proved to be the best from the study discussed in the previous section. For this test, the Courant number was kept at 1.0. It is evident from Figure 6 that the proposed method is in good agreement with the well-established MOC scheme. The advantages of the FVM are it works for lesser Courant numbers without significant numerical dissipation, whereas the MOC completely fails for Courant number less than 1.0. This particular characteristic of FVM helps in solving such problems where the wave speed varies quite a bit along the length of pipe due to change in pipe thickness caused by wear and tear, deposit thickness and corrosion. In addition, the FVM is useful where wave speed changes along the pipes caused by the presence of air cavities (Daude et al. 2018).

Figure 6

(a) Comparison of FVM first-order scheme with the MOC. (b) Comparison of MUSCL-Hancock scheme with MOC.

Figure 6

(a) Comparison of FVM first-order scheme with the MOC. (b) Comparison of MUSCL-Hancock scheme with MOC.

From the comparison, it is evident that the proposed scheme is a good alternative for solving hydraulic transients in a simple pipe and large reservoir environment. The plot shows the dissipation of pressure head due to steady friction. Also, the linear increase and decrease in the water pressure head can be attributed to linepack pressure as described in Hanmaiahgari et al. (2019).

Codes for MOC and FVMs were programmed in MATLAB. The computations were carried out using an Intel Pentium (R) quad core CPU 3.1 GHz processor with 12 GB RAM. This validation shows an excellent agreement between pressures computed by the MOC and Godunov first- and second-order schemes. The computational time of the first-order scheme is comparable to that of MOC, whereas the second-order scheme requires 15 times that of MOC (refer Table 3). Hence, in the case of linear (serial) programming, the first-order FV method shows more applicability over the second-order scheme in terms of time consumption. But the FVM is well known for its structured algorithms and suitability for parallel computing. The implementation of the MATLAB parallel computing toolbox (MATLAB and Parallel Computing Toolbox User's Guide Release 2018) reduced the computational time of the second-order FVM comparable to the first-order FVM and computational time of the first-order FVM to that of MOC. This demonstrated that parallel computing improves the applicability of FVM for real-time monitoring applications. The results also demonstrate that first-order FVM results are acceptable for modeling water hammer resulting from sudden valve closures at downstream and sudden pump shutdowns at upstream as tested in this paper.

Table 3

Time comparison between MOC and FVM

MethodElapsed time (s)
Method of characteristics 0.21 
Godunov first-order method 0.37 
Godunov second-order method 2.73 
MethodElapsed time (s)
Method of characteristics 0.21 
Godunov first-order method 0.37 
Godunov second-order method 2.73 

Validation with experimental results

The proposed FVM model was compared with experimental results presented in Adamkowski & Lewandowski (2006), Meniconi et al. (2014), Silva-Araya & Chaudhry (2001), and Pezzinga (1999, 2009). The results of the comparisons are presented in the following section. It has been observed that both first- and second-order Godunov schemes give exact results for all four cases. Hence, FVM denotes both the first-order and second-order in the figures for space optimization.

Comparison with Adamkowski & Lewandowski (2006): Adamkowski & Lewandowski (2006) conducted several tests for transient pipe flow with a wide range of Reynold numbers. They have used a specially prepared test rig at IMP PAN in Gdansk and study the water hammer phenomenon. It consisted of a 98.11 m-long copper pipe of diameter 0.016 m and wall thickness of 0.001 m. For quick closure, a spring ball valve was used at the end and valve closure time was ensured below 0.003 s that is 2% of the half period of propagation of wave back and forth. The unsteady pressures were measured at four locations at L/4 spacing. Those transducers' properties and experimental conditions are given in Tables 4 and 5, respectively.

Table 4

Transducer properties Adamkowski & Lewandowski (2006) experiments

PropertyRange
Measuring range 0–4 Mpa 
Transmitted frequency band 0–2 kHz 
Precision class 0.2% 
PropertyRange
Measuring range 0–4 Mpa 
Transmitted frequency band 0–2 kHz 
Precision class 0.2% 
Table 5
PropertyValue
Initial flow velocity, V0  
Initial discharge,   
Reynolds number, Re 15,800 
Water bulk modulus,  2.1 109 Pa 
Density of water,   
Pipe line poisson coefficient,  0.35 
Density,  
Young's modulus,   
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,  , 127.497 m 
Pipe length,   
Friction factor,   
PropertyValue
Initial flow velocity, V0  
Initial discharge,   
Reynolds number, Re 15,800 
Water bulk modulus,  2.1 109 Pa 
Density of water,   
Pipe line poisson coefficient,  0.35 
Density,  
Young's modulus,   
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,  , 127.497 m 
Pipe length,   
Friction factor,   

A turbine flowmeter of range was used for measuring the discharge. A steady pressure head was maintained at 1,250 kPa at inlet reservoir. Temperature and kinematic viscosity of water in the experiment were recorded as 22.6 °C and , respectively. The experimental time duration of pressure measurement was 3.78 s. A data logger of frequency 2.5 kHz was used for capturing the data. More details about the experimental setup can be found at Adamkowski & Lewandowski (2006).

The comparison between the FVM and experimental pressures at the valve is shown in Figure 7. The first-order FVM and second-order FV MUSCL-Hancock are giving the same results. It can be concluded that for the first cycle, the model has predicted both the positive and negative peak pressures as well as the rate of dissipation correctly; however, for longer time periods, the model overestimates the pressure peaks because of lesser dissipation with the steady-state friction. The experimental valve closure is fast, but not instantaneous resulting in a small drift between the model and experimental data.

Comparison with Meniconi et al. (2014): Meniconi et al. (2014) experimented with a high Reynold number unsteady flow without any pressure surge control. They have used a single main pipeline and studied the sudden pump shutdown effect on the pressures. Details of the model can be found in Meniconi et al. (2014). The experimental conditions are summarized in Table 6.

Table 6

Model properties (Meniconi et al. (2014))

PropertyValue
Initial flow velocity, V0  
Initial discharge,   
Reynolds number, Re 240,446.39 
Density of water,   
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,   
Pipe length,   
Friction factor,   
Pipe material Steel pipe with the wall thickness of 6.5 mm 
PropertyValue
Initial flow velocity, V0  
Initial discharge,   
Reynolds number, Re 240,446.39 
Density of water,   
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,   
Pipe length,   
Friction factor,   
Pipe material Steel pipe with the wall thickness of 6.5 mm 

The same experiment was simulated with the proposed first-order and MUSCL-Hancock second-order Godunov FV scheme as shown in Figure 8. The comparison of the simulated data with the experimental data accurately represents the same initial peak pressure of 400 m as well as down surge pressure of 100 m along with the proper damping rate.

Figure 8

Validation of FVM with the experiment of Meniconi et al. (2014).

Figure 8

Validation of FVM with the experiment of Meniconi et al. (2014).

Comparison with Silva-Araya & Chaudhry (2001): Silva-Araya & Chaudhry (2001) developed an energy dissipation model for computation of unsteady friction losses for fully rough pipe and compared it experimentally. Their experimental data were used to validate and compare the presented FV model. A brief description of the experimental model is provided here. The reader can refer to Silva-Araya & Chaudhry (2001) for more details.

In this model, a constant head reservoir open to the atmosphere is connected by a horizontal pipe with a motorized butterfly valve at the other end. The pipeline is made of steel for strength considerations. The pressure transient was produced by the slow closure of the butterfly valve at the downstream end. Three equally spaced pressure transducers were used for measuring the pressure transients. Valve closure time was taken as 1.085 s. The model properties are summarized in Table 7.

Table 7
PropertyValue
Initial flow velocity, V0 2.035 m/s 
Initial discharge,   
Reynolds number, Re  
Density of water,  998 kg/m3 
Pipe line length, L 32 m 
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,   
Friction factor,   
PropertyValue
Initial flow velocity, V0 2.035 m/s 
Initial discharge,   
Reynolds number, Re  
Density of water,  998 kg/m3 
Pipe line length, L 32 m 
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,   
Friction factor,   

The comparison of experimental and simulated pressure is presented in Figure 9, which shows that the current FVM model predicts the peak pressure perfectly, but the peak occurs slightly prior to the experiment peak occurence. This may be due to the small error in the time measurement of the experimental pressure. The FVM results agree well with the first peak and its time as measured by Silva-Araya & Chaudhry (2001).

Figure 9

Validation of FVM with the slow valve closure experiment of Silva-Araya & Chaudhry (2001).

Figure 9

Validation of FVM with the slow valve closure experiment of Silva-Araya & Chaudhry (2001).

Comparison with Pezzinga (1999, 2009): Pezzinga (1999) used a pipeline network for generating unsteady flow transients. The pipe system was a centrifugal pump fed network with a zinc-plated steel pipe. A 1.0 m3 tank was used at the inlet of the main pipe. Electromagnetic flowmeter and pressure transducers were used for measuring flow and pressures at different locations of the pipe, respectively. The details of transducer properties are given in Table 8.

Table 8

Transducer properties Pezzinga (1999, 2009) experiments

PropertyRange
Measuring range 0–10 bar 
Transmitted frequency band 200 Hz 
Precision class 0.5% 
PropertyRange
Measuring range 0–10 bar 
Transmitted frequency band 200 Hz 
Precision class 0.5% 

A number of experiments were conducted by Pezzinga (1999) and for each experiment a valve closure time of 0.04 s was adopted. For the current study, a particular case of those experiments, i.e., Re = 9,250, was considered to validate the proposed FVM schemes. The experimental particulars of the model are presented in Table 9. The experimental pressure at the valve is compared with the proposed FVM model as shown in Figure 10.

Table 9

Model properties Pezzinga (1999, 2009)

PropertyValue
Initial flow velocity,   
Initial discharge,   
Reynolds number, Re  
Density of water,   
Young modulus,   
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,  37.74 m 
Pipe length,   
Friction factor,   
PropertyValue
Initial flow velocity,   
Initial discharge,   
Reynolds number, Re  
Density of water,   
Young modulus,   
Internal pipe diameter and area, D and A  and  
Pipe wall thickness,   
Celerity,   
Initial reservoir head,  37.74 m 
Pipe length,   
Friction factor,   
Figure 10

Validation of FVM (first and second order) with Pezzinga (1999, 2009).

Figure 10

Validation of FVM (first and second order) with Pezzinga (1999, 2009).

Figure 10 shows that the maximum and the minimum pressures, as well as the damping rate, are almost similar in simulations and the experiment. This validates the current numerical model and ensures its applicability to practical pipeline network problems. The small drift between the model and the experiment is due to fast but not instantaneous closure of the valve in the experiment.

Unsteady friction using the convolution intregral method

An unsteady friction model based on convolution theory as proposed by Vardy & Brown (2007) was implemented in the proposed FVM. This particular Vardy & Brown (2007) was selected for its applicability for full range of roughness sizes and wide Reynolds numbers. A brief description of the convolution-based unsteady model was presented in the section ‘Unsteady friction computaton using the convolution theory’. The proposed investigation demonstrates the implementabiliy of the convolution-based unsteady friction model in the FVM for practical purposes. The experimental data (Adamkowski & Lewandowski 2006) for sudden valve closure was chosen for the validation. A brief description of the experiment has already been presented in the section ‘Validation with experiment results’ (refer Table 6 and Figure 8). More specifics of experiments can be obtained from Adamkowski & Lewandowski (2006). The convolution method was compared with IAB models to evaluate the accuracy and computational time. Various IAB models proposed by Brunone et al. (1995), Vítkovský et al. (2006) and Pezzinga (2000) were considered. Vardy & Brown's (1996) decay coefficient as a function of Reynolds number and with a value of 0.0144 was implemented in models by Brunone et al. (1995) and Vítkovský et al. (2006). Pezzinga (2000) proposed design plots for the decay coefficient, which is a function of discharge and a dimensionless longitudinal step. The estimated decay coefficient 0.0182 is implemented in the Pezzinga (2000) model. Finally, the comparison between numerical and experimental pressures for 3.75 s of simulation is presented in Figure 11. The analysis of Figure 11 shows that IAB models were not able to adequately capture the shape of the observed pressure signal. However, the convolution method has accurately predicted the damping as well as the shape. Regarding the computational time of the experimental simulation, the IAB models required only 10 s, whereas the convolution method required 3 min. The computational time of the convolution-based method has been minimized by following the methodology given in Vitkovsky et al. (2004, 2006). This simulation illustrates that Vardy & Brown's (2007) method is highly accurate for computing unsteady friction during sudden valve closures and can be effectively implemented in the FVM-based dynamic model.

Figure 11

Comparison between FVM with the convolution method, IAB methods and experiment.

Figure 11

Comparison between FVM with the convolution method, IAB methods and experiment.

CONCLUSIONS

The FVM is becoming an emerging way to explore many engineering problems with ease and rectify many boundaries that cannot be overcome by other methods. The state-of-the-art method is known for very low numerical diffusion, efficient handling of shocks and discontinuities in the simulations. In this direction, the proposed study investigated the accuracy and computational times of first-order and second-order FVM schemes for modeling of unsteady flows in pipelines. The detailed methodology of first- and second-order FVM schemes for modeling the water hammer problem was presented in this paper. The second-order solution of the Godunov scheme requires an efficient slope or flux limiter for error minimization and time optimization. The study examined a range of limiters and found that the MINMOD limiter is the best for modeling water hammer in terms of computational time and accuracy.

The first- and second-order FVMs were compared with the MOC and experimental water hammer pressure measurements available in the literature. Both the FV methods accurately predicted the numerical and experimental results. However, the second-order FVM requires a larger computational time. The parallel computing of first-order and second-order FVM is implemented to reduce the computational time. The parallel computing of the second-order scheme reduced its computational time to the level of the first-order scheme and similarly computational time of first-order FVM reduced to that of MOC. This demonstrated that parallel computing enables the applicability of FVM for real-time monitoring applications. The proposed study also demonstrates that first-order FVM results are acceptable for modeling water hammer resulting from sudden valve closures at downstream and sudden pump shutdowns at upstream. Thus, the study presented a faster and more accurate FVM which is comparable to that of MOC in terms of computational time and precision, therefore it is a good substitute for MOC.

One of the novelties of the proposed study is the implementation of a more complex weighting function-based unsteady friction model in the FVM to capture real pressure dissipation. The comparison with experimental data proved that the convolution integral method is highly accurate for computing pressure damping during sudden valve closures.

DATA AVAILABILITY STATEMENT

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

REFERENCES

Adamkowski
A.
Lewandowski
L.
2006
Experimental examination of unsteady friction models for transient pipe flow simulation
.
J. Fluids Eng.
128
,
1351
1363
.
doi:10.1115/1.2354521
.
Afshar
M. H.
Rohani
M.
2008
Water hammer simulation by implicit method of characteristic
.
Int. J. Press. Vessels Pip.
85
(
12
),
851
859
.
Brunone
B.
Golia
U. M.
Greco
M.
1995
The effects of two dimensionality on pipe transients modeling
.
J. of Hydraul. Eng., ASCE
,
121
(
12
),
906
912
.
Chaudhry
M. H.
2014
Applied Hydraulic Transients
, 3rd edn.
Springer
,
New York
.
Chen
T.
Ren
Z.
Xu
C.
Loxton
R.
2015
Optimal boundary control for water hammer suppression in fluid transmission pipelines
.
Comput. Math. Appl.
69
(
94
),
275
290
.
Ghidaoui
M. S.
Karney
B. W.
1994
Equivalent differential equations in fixed grid characteristics methods
.
J. Hydraul. Eng.
120
,
1159
1175
.
Hanmaiahgari
P. R.
Kottam
R. R.
Kaushik
M.
2019
Estimation and examination of linepack pressures in long water pipelines
.
Sādhanā
44
(
101
),
2019
.
Hwang
Y. H.
Chung
N. M.
2002
A fast Godunov method for the water-hammer problem
.
Int. J. Numer. Methods Fluids
40
(
6
),
799
819
.
Leon
A. S.
Ghidaoui
M. S.
Schmidt
A. R.
Garcia
M. H.
2008
Application of Godunov-type schemes to transient mixed flows
.
J. Hydraul. Res.
47
(
2
),
147
156
.
Li
X.
Tang
X.
Zhu
M.
Shi
X.
2019
1D-3D coupling investigation of hydraulic transient for power-supply failure of centrifugal pump-pipe system
.
J. Hydroinf.
21
(
5
),
708
726
.
https://doi.org/10.2166/hydro.2019.122
.
MATLAB and Parallel Computing Toolbox User's Guide Release
2018
The MathWorks, Inc., Natick, MA
.
Meniconi
S.
Duan
H. F.
Brunone
B.
Ghidaoui
M. S.
Lee
P. J.
Ferrante
M.
2014
Further developments in rapidly decelerating turbulent pipe flow modeling
.
J. Hydrual. Eng.
140
(
7
),
04014028
.
https://doi.org/10.1061/(ASCE)HY.1943-7900.0000880
.
Niroomandi
A.
Borghei
M.
Bohluly
A.
2012
Implementation of time splitting projection method in water hammer modeling in deformable pipes
.
Int. J. Press. Vessels Pip.
98
,
30
42
.
Pezzinga
G.
1999
Quasi-2D model for unsteady flow in pipe networks
.
J. Hydrual. Eng.
125
(
7
),
676
685
.
Pezzinga
G.
2009
Local balance unsteady friction model
.
J. Hydrual. Eng.
135
(
1
),
45
56
.
doi:10.1061/ASCE0733-9429
.
Seck
A.
Fuamba
M.
Kahawita
R.
2017
Finite-volume solutions to the water-hammer equations in conservation form incorporating dynamic friction using the Godunov scheme
.
J. Hydraul. Eng.
143
(
9
),
04017029
.
doi:10.1061/(ASCE)HY.1943-7900.0001333
.
Silva-Araya
W. F.
Chaudhry
M. H.
2001
Unsteady friction in rough pipes
.
J. Hydrual. Eng.
127
(
7
),
607
618
.
Szydlowski
M.
2002
Finite Volume Method for Water Hammer Simulation
.
TiASWiK'02. In: Proceedings of International Scientific and Technical Conference on Technology, Automation and Control of Wastewater and Drinking Water Systems (K. Konarczak & D. Trawicki, eds). TiASWIK'02. Gdańsk-Sobieszewo, June 19–21, 2002. Gdańsk Technol. Univ., Gdańsk, pp. 159–165
.
Vardy
A. E.
Brown
J. M. B.
1996
On turbulent, unsteady, smooth pipe friction
. In:
Proceedings of 7th International Conference on Pressure Surges and Fluid Transients in Pipelines and Open Channels
.
BHR Group
,
Harrogate
,
UK
, pp.
289
311
.
Vardy
A. E.
Brown
J. M. B.
2007
Approximation of turbulent wall shear stresses in highly transient pipe flows
.
J. Hydrual. Eng.
133
(
11
),
1219
1228
.
doi:10.1061/ASCE0733-94292007133:111219
.
Vítkovský
J.
Stephens
M.
Bergant
A.
Lambert
M. F.
Simpson
A. R.
2004
Efficient and accurate calculation of Zielke and Vardy-Brown unsteady friction in pipe transients
. In
9th International Conference on Pressure Surges
,
Chester, UK
.
Vítkovský
J. P.
Stephens
M.
Bergant
A.
Simpson
A.
Lambert
M. F.
2006
Numerical error in weighing function-based unsteady friction models for pipe transients
.
J. Hydrual. Eng.
132
(
7
),
709
721
.
Wylie
E. B.
Streeter
V. L.
1993
Fluid Transients in Systems
.
Prentice-Hall
,
Englewood Cliffs
.
Yazdi
S. R. S.
Mastorakis
N. E.
Abbasi
A.
2007
Water hammer modeling by Godunov type finite volume method
.
Int. J. Math. Comp. Sim.
350
355
.
Zhang
C.
Zecchin
A. C.
Lambert
M. F.
Gong
J.
Simpson
A. R.
2018a
Multi-stage parameter-constraining inverse transient analysis for pipeline condition assessment
.
J. Hydroinf.
20
(
2
),
281
300
.
https://doi.org/10.2166/hydro.2018.154
.
Zhao
M.
Ghidaoui
M. S.
2004
Godunov-type solutions for water hammer flows
.
J. Hydrual. Eng.
130
(
4
),
341
348
.