Finite-Volume Solutions to the Water-Hammer Equations in Conservation Form Incorporating Dynamic Friction Using the Godunov Scheme

In the preliminary design of hydraulic systems, practicing engineers use a general framework to design equipment to prevent and/or mitigate any excessive pressures caused by water hammer. The detailed design further attempts to refine the preliminary concept to elaborate the complete system response using various tools such as numerical modeling. In general, numerical simulation tools assume that the friction factor is considered steady or quasi-steady; this has the tendency to overestimate the water-hammer peak pressures because the primary mechanism that may significantly affect pressure waveforms is the unsteady friction. This results in oversizing of surge control equipment.

The origin of the theory of water hammer goes back to the contributions of Menabrea (1858), who published a short note on the calculation of water pressures (Anderson 1976). However, because of the mathematical rigor and significant theories presented, the papers of Michaud (1878, 1903), Allievi (1903, 1913, 1932), Schnyder (1932), and Jaeger (1933) have been the source of inspiration for all studies on water hammer. Water hammer in hydropower systems and its stability problems were reviewed by Lescovich (1967), Roche (1975), and Chaudhry (1987). Research by Chaudhry et al. (1985), Finnemore and Franzini (2002), and Azhdari Moghaddam (2004) provide implicit equations for analyzing the water-hammer wave. Chaudhry and Silvaaraya (1992) investigated the stability of water level oscillations during transient conditions in a surge tank. Kim (2010) applied a genetic algorithm (GA) coupled with an impulse response technique to derive the impedance functions for pipe systems equipped with a surge tank. Ramadan and Mustafa (2013) investigated the effects of different parameters such as the friction loss coefficient, tank dimensions, and total discharge on transient overpressures. Seck and Fuamba (2015) have developed the analytical equations to design cylindrical and conical surge tanks, and have provided user-friendly diagrams that are based on the equations developed. Their paper bridges the gap between the design concept and the detailed design phase for surge tanks but cannot be used as a general purpose tool for transient analysis because certain key parameters such as the dynamic friction are not incorporated. Guinot (2003) focused on the theoretical and practical implementation of the Godunov approach to simulate water hammer using finite volumes with steady friction that gives rise to a nonhyperbolic source term.

In general, water-hammer models assume that the friction factor is steady or quasi-steady. However, in transient flows with a constant pipe diameter, the one-dimensional velocity profile is distorted from its assumed steady state shape due to the large accelerations (varying across the cross section) that are imposed on the fluid particles. The equivalent friction then exceeds the value assumed at steady state. Brunone et al. (1991) and subsequently Bergant et al. (2001) decomposed the dynamic friction into the sum of a quasi-steady and unsteady parts. Bergant et al. (2008a, b) have experimentally investigated key parameters (unsteady friction, cavitation, fluid–structure interaction, viscoelastic behavior of the pipe wall, leakage, and blockage) that may affect the pressure waveform predicted by the classical theory of water hammer. Bousso and Fuamba (2013) have numerically represented the unsteady friction in transient two-phase flow with a hyperbolic source term. Shamloo et al. (2015) review a quasi-steady model and four unsteady friction models (Zielke 1968; Vardy and Brown 2003; Trikha 1975; Brunone et al. 1991) for transient pipe flow. The Vardy and Brown model is limited to smooth pipes while the Zielke model is developed for transient laminar flow. The Trikha (1975) model presents a simplification of the Zielke (1968) model. The Brunone model appears to be the only one developed for all types of transient flows and all roughness values. Shamloo and Mousavifard (2015) use a two-dimensional finite-difference method incorporating the $k\text{-}\omega$ turbulence model to study the dynamics of turbulence during different periods of water hammer in a polymeric pipe. However, comparison of their model with experimental tests reveal that the model underestimates the peak water-hammer pressures, which could pose a significant risk for pipe safety.

Although derived from the principles of conservation of mass and momentum, the water-hammer equations integrating dynamic friction are almost never expressed in conservation form. This is because the presence of the dynamic friction terms preclude presentation in the integral (conservation) form and also cause a source term to appear. The finite-volume method that was developed for hyperbolic equations in conservation form (e.g., the Euler equations in fluid dynamics) has to be adapted to tackle systems that are not in proper conservation form. Hyperbolic systems admit weak solutions in the form of discontinuities or shocks. The treatment of these discontinuities requires special treatment of the flux function to avoid spurious oscillations. A large body of literature is available on finite-volume schemes developed to handle the solution in the vicinity of these shocks. Among the most modern numerical schemes are total variation diminishing (TVD) (Zijlema and Wesseling 1998), essentially nonoscillatory (ENO) (Harten et al. 1987), and weighted essentially nonoscillatory (WENO) schemes (Liu et al. 1994; Jiang and Shu 1995). One limitation of the TVD schemes is unsatisfactory performance near extrema (Osher and Chakravarthy 1984; Titarev and Toro 2003; Park et al. 2010). Essentially nonoscillatory and WENO schemes have been developed to overcome this limitation and provide a better scheme that can handle both sharp interfaces and smooth gradients (Shu 1998). According to Gallerano and Cannata (2011) and Gallerano et al. (2012), ENO and WENO shock capturing higher-order schemes for the solution of hyperbolic systems can be considered as a further development of ideas that gave rise to the TVD schemes. The original contribution of this paper is to present the derivation of the water-hammer equations in conservation form incorporating dynamic friction. The hyperbolic system is then solved using a finite-volume scheme. The Godunov scheme was retained to obtain the numerical solutions to the present problem.

It is hoped that this paper will help practical engineers who are at the detailed design phase to further safely optimize any incorporated surge controls by assisting in the development and numerical simulation of the water-hammer equations with dynamic friction.

This paper focuses on implementation of the Godunov approach to one-dimensional hyperbolic systems of conservation laws that describe the phenomenon of water-hammer. It is first order accurate in both space and time. In the present context, given the simple geometry, the basic Godunov scheme is considered adequate by the authors. First, the computation of the fluxes is detailed followed by a description of the complete numerical algorithm. Finally, two case studies are investigated with comparisons being made between an analytical solution, the numerical solution with quasi-steady friction, the numerical solution with dynamic friction, and experimental results.

For a pipe of constant cross section, the one-dimensional flow equations presented in vector form (Guinot 2003) arewhere $t=\mathrm{time}$; $x$ = unit vector in the $x$-direction; $\mathbf{U}$ = flow variable vector; $\mathbf{F}$ = flux vector in the $x$-direction; $\mathbf{S}$ = source term vector; $\mu$ = mass of fluid per unit length of pipe; $Q_m$ = mass discharge; $A_0$ = cross-sectional area of pipe; $p$ = pressure; $\mathbf{V}$ = fluid velocity; and $f$ = friction coefficient dependant on the pipe roughness and the fluid viscosity.

For unsteady friction in transient flow with a constant pipe diameter, the friction factor $f$ is split into the sum of the quasi-steady $f_q$ and unsteady $f_u$ part i.e., $f=f_q+f_u$ as in the model by Brunone et al. (1991) and modified by Bergant et al. (2001)where $a$ = celerity of the pressure waves; $D$ = diameter of the pipe; and $C^{*}$ = Vardy’s shear decay constant depending on the flow regime (Reynolds number $\mathsf{R}$) expressed as

Eq. (2) may be rewritten in terms of the flow variables $\mu$ and $Q_m$ aswhere $\psi$ is expressed as

For turbulent flow, the Colebrook-White equation may be iteratively solved for the quasi-steady part $f_q$ of the fiction factor

In which $k_S$ = absolute roughness of the pipe.

Streeter and Wylie (1993) provide the following expression for the celerity (sonic velocity) $a$ in a circular pipe:where $E_V$ = bulk modulus of elasticity; $\rho$ = mass density; $e$ and $E$ = pipe wall thickness and Young’s modulus of elasticity for the pipe material, respectively.

Eq. (7a) may be rewritten

By taking into consideration the unsteady friction, Eq. (1) may be written as

In pipe flows, the flow velocity $V$ is negligible compared with $a$ (the celerity of the pressure waves). The authors also note that $2\mu /2\mu +\sqrt{C^{*}}D\approx 1$ because the term $\sqrt{C^{*}}D$ is negligible compared with $2\mu$. Eq. (8) may now be rewritten

In conservation form the governing Eq. (9) may be written, respectivelywhere $k=\sqrt{C^{*}}D/2\mu +\sqrt{C^{*}}D$

This is the general form of the water-hammer equations with unsteady friction.

In Eq. (11), $\mathbf{A}=\partial \mathbf{F}/\partial \mathbf{U}$ is the Jacobian matrix of $\mathbf{F}$ with respect to the matrix $\mathbf{U}$. In the expression for the matrix $\mathbf{A}$, $k{A}_{0}p/\mu \ll a^2$. The general Eq. (10) is equivalent to the Guinot Eq. (1) if unsteady friction is not considered. The inclusion of the unsteady friction results in the appearance of the quasi-linear term $\mathbf{R}\partial \mathbf{U}/\partial x$ in the hyperbolic system [Eq. (10)] thus destroying the integral form of the equations.

The Godunov method is a conservative numerical scheme for solving partial differential equations (PDE) in computational fluid dynamics. This conservative finite-volume method solves an exact or approximate Riemann problem (depending on the friction formulation) at each intercell boundary. It is first order accurate in both space and time. The next section describes the solution process for these equations with the Godunov method.

Guinot (2003) as well as Bousso and Fuamba (2013) have provided the procedure for solving PDE’s by Godunov’s method using the time splitting technique. Eq. (10) is solved in three steps. The first step uses the Guinot solution for the conservation part of the PDE [Eq. (10)]

The second step solves for the nonconservative term of Eq. (10)

The final step uses the Toro (2001) and Guinot (2003) treatment of the nonhyperbolic source term of Eq. (10)

The conservation component of Eqs. (10) and (12) is the same conservation component as Eq. (1) in Guinot (2003). For the internal cells, Guinot (2003) constructs the Riemann problem $\mathbf{U}(x,t^n)$ and presents the solution ${\mathbf{U}}_{i+1/2}^{n+1/2}$ of this Riemann problem in Eq. (15) computing the flux ${\mathbf{F}}_{i+1/2}^{n+1/2}$ between the time intervals $t^n$ and $t^{n+1}$ according to Eq. (16)where $p_{\mathrm{ref}}$ is a reference pressure at which the density ${\rho }_{\mathrm{ref}}$ is known. ${\mu }_{\mathrm{ref}}$ is calculated as ${\mu }_{\mathrm{ref}}=A_0{\rho }_{\mathrm{ref}}$.

The flux solutions for a prescribed pressure $p_b$ at the left-hand ($p_{b,L}$) and the right-hand ($p_{b,R}$) boundary are respectively given by

In Eqs. (17) and (18), $p_b$ represents a pressure to be prescribed at the boundary. The mass of fluid per unit length of pipe ${\mu }_{1/2}^{n+1/2}$ (at the left hand) and ${\mu }_{N+1/2}^{n+1/2}$ (at the right hand) is obtained from the prescribed pressure $p_b$ using

The fluxes for a prescribed discharge $Q_b$ at the left-hand and the right-hand boundary are respectively given bywhere ${\mu }_{1/2}^{n+1/2}$ in Eq. (20) and ${\mu }_{N+1/2}^{n+1/2}$ in Eq. (21) are obtained, respectively, from

For all cells, the balance is performed omitting the source term using

The intermediate solution ${\mathbf{U}}_i^{n+1,x}$ is taken as a starting point for the computation of the hyperbolic source term. This nonhyperbolic source term is solved according to the procedure detailed subsequently.

In Eq. (13), $\mathbf{R}=ka\psi \left[\begin{array}{cc}0& 0\\ -V& 1\end{array}\right]$. The two eigenvalues ${\lambda }^{(1)}$ and ${\lambda }^{(2)}$ of the matrix $\mathbf{R}$ are therefore

The matrix $\mathbf{K}$ of the $\mathbf{R}$ eigenvectors is given by

This yields the following generalized Riemann invariants

Noting that $d{Q}_m=\mu dV+Vd\mu$, the system Eqs. (28) and (29) may be rewritten as

Consider the Riemann problem as follows:

The solution ${\mathbf{U}}^h$ of this Riemann problem is given bywhere ${\mathbf{U}}_{i+1/2}^{*}=\frac{1}{A_0}\left[\begin{array}{c}{\mu }_{i+1}^{n+1,x}\\ A_0^2\frac{{\mu }_{i+1}^{n+1,x}}{{\mu }_i^{n+1,x}{\mu }_i^{n+1,x}}{Q}_{m_i}^{n+1,x}\end{array}\right]$where ${\mathbf{U}}_{i+1/2}^{*}=\frac{1}{A_0}\left[\begin{array}{c}{\mu }_i^{n+1,x}\\ A_0^2\frac{{\mu }_i^{n+1,x}}{{\mu }_{i+1}^{n+1,x}{\mu }_{i+1}^{n+1,x}}{Q}_{m_{i+1}}^{n+1,x}\end{array}\right]$

The value of ${\mathbf{U}}_i^{n+1,x}$ is obtained from

By substituting Eqs. (33) and (34) into Eq. (35), the following equation is obtained:

Taking the provisional solution ${\mathbf{U}}_i^{n+1,h}$ as the initial state, the solution ${\mathbf{U}}_i^{n+1}$ of the nonhyperbolic source term is presented by Toro (2001), Guinot (2003), and Bousso and Fuamba (2013) in the following form:where $\mathbf{S}({\mathbf{U}}_i^{n+1,h})$ is calculated by using ${\mathbf{U}}_i^{n+1,h}$ in Eq. (1) of the source term.

Because the scheme is explicit, the CFL condition for stability must be respected, i.e., the maximum permissible time step $\mathrm{\Delta }{t}_{\mathrm{Max}}$ as presented by Guinot (2003) is as follows:

The water-hammer wave attenuation ratio $\eta$ is computed from

As $t$ increases ($t\to \infty$) and assuming no pipe failure, the water-hammer wave is completely attenuated in the pipe.

Then, ${\eta }_i^{n+1}$ approaches 100% according to

The water-hammer equations in conservation form incorporating dynamic friction is developed and numerically solved using a finite-volume formulation. The computational algorithm based on the Godunov scheme for one-dimensional hyperbolic systems is presented in some detail. Introduction of the dynamic friction results in the appearance of a source term in the hyperbolic system of governing partial differential equations. Two case studies have been presented to compare and contrast the separate impacts of steady friction only and of dynamic friction on wave attenuation in time. The findings indicate that inclusion of the dynamic friction reduces the peak water-hammer pressures when compared with the standard quasi-steady assumption for the friction. The inclusion of the dynamic friction results in a better agreement between calculated and measured values. In these two case studies, grid size and time step were reduced to test for convergence; the numerical solution presented is the converged solution. Computational time for each of the cases studied took only 3 min on a PC (operating system: Windows 7, 64 bits; CPU: Intel(R) Xeon(R) CPU E5620 @ 2.40 GHz x 2 processors; mainboard: ASUS Z8NA-D6(C); memory: 12280 MB triple channel DDR3 ECC @ 1333 MHz). Thus several variants may be studied in a short period of time. The efficient algorithm introduced in this paper may be adapted, with some modifications, to incorporate higher-order schemes that may yield less computational effort. Because the principal objective of this paper was to refine the physics of the phenomenon, extensive numerical experimentation with different schemes was not undertaken. However, for more precision in engineering applications, other key parameters such as cavitation, fluid–structure interaction, viscoelastic behavior of the pipe wall, leakage, and blockage have to be properly modeled and incorporated into the formulation. This paper may be used as a basis for the study of real mass variations in a surge tank. It should also prove useful as a tool for the optimal design of surge controls.

The following symbols are used in this paper:

$\mathbf{A}$

Jacobian matrix of the flux $\mathbf{F}$ respecting to $\mathbf{U}$;

$A_0$

cross-sectional area of pipe;

$a$

celerity of the pressure waves;

$C^{*}$

Vardy’s shear decay constant;

$D$

diameter of the pipe;

$d$

(operator) differential;

$E$

Young’s modulus of elasticity for the pipe material;

$E_V$

bulk modulus of elasticity;

$e$

pipe wall thickness;

$\mathbf{F}$

vector flux in the $x$-direction;

$f$

friction factor;

$f_q$

quasi-steady part of the friction factor;

$f_u$

unsteady part of the friction factor;

$k_S$

absolute roughness of the pipe;

$N$

number of cells in the computational domain;

$p$

pressure;

$Q$

volume discharge;

$Q_m$

mass discharge;

$\mathbf{R}$

matrix for the hyperbolic source term;

$\mathsf{R}$

Reynolds number;

$\mathbf{S}$

vector source term;

$t$

time;

$\mathbf{U}$

vector variable;

$V$

fluid velocity;

$x$

unit vector in the $x$-direction;

$\mathrm{\Delta }t$

computational time step;

$\mathrm{\Delta }x$

cell size in the $x$-direction;

$\mu$

mass of fluid per unit length of pipe; and

$\rho$

mass density.

Subscripts

$b$

value to be prescribed at a boundary;

$b$, $L$

value to be prescribed at the left-hand boundary;

$b$, $R$

value to be prescribed at the right-hand boundary; and

$i$

cell number.

Superscripts

$n$

time level;

$n$, $h$

numerical solution obtained at the time level $n$ after the solution of the hyperbolic source term of the equation;

$n$, $x$

numerical solution obtained at the time level $n$ after the solution of the conservation part of the equation in the $x$-direction;

$n+1/2$

Average value between the time level $n$ and $n+1$;

$\infty$

time level if t gets larger; and

*

intermediate region in the solution of the Riemann problem.

The authors gratefully acknowledge the partial financial support of Natural Sciences and Engineering Research Council of Canada (NSERC).