Efficient approach toward the application of the Godunov method to hydraulic transients

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.

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.

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

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:

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:

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

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

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:

MINMOD:

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

Test model description

Reservoir boundary condition at upstream:

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

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.

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.

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

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.

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 m³/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⁻⁵ 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.

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.

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.

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

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.

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.

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.

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.

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.

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

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 m³ 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.

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.

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.

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.

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.

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