Abstract
The main contribution of the paper is to incorporate pipe-wall viscoelastic and unsteady friction in the derivation of the water-hammer solutions of non-conservative hyperbolic systems with conserved quantities as variables. The system is solved using the Godunov finite volume scheme to obtain numerical solutions. This results in the appearance of a new term in the mass conservation equation of the classical governing system. This new numerical algorithm implements the Godunov approach to one-dimensional hyperbolic systems of conservation laws on a finite volume stencil. The viscoelastic pipe-wall response in the mass conservation part of the source term has been modeled using generalized Kelvin–Voigt theory. For the momentum part of the source term a fast, robust and accurate numerical scheme linked to the Lambert W-function for calculating the friction factor has been used. A case study has been used to illustrate the influence of the various formulations; a comparison between the classical solution, the numerical solution including quasi-steady friction, the numerical solution incorporating the viscoelastic effects, and measurements are presented. The inclusion of viscoelastic effects results in better agreement between the measured and solved values.
HIGHLIGHTS
Incorporate the viscoelastic pipe and unsteady friction in water-hammer formulations.
Present the solutions of non-conservative hyperbolic systems with conserved quantities as variables.
Adapt the finite volume approach, which was developed for hyperbolic equations to tackle systems that are not in ‘proper’ conservation form.
Give a new numerical algorithm in finite volume for study several variants.
NOMENCLATURE
The following notation is used in this contribution.
Variables and operators
Celerity of the pressure waves
Jacobian matrix of F respecting to U
Cross-sectional pipe area
Dimensionless parameter
(operator) Differential
Diameter of the pipe
Pipe wall width
Young's modulus of the pipe material
Modulus of the spring of
-element elasticity
Bulk modulus of elasticity
Vector flux in the x-direction
Quasi-steady part of the friction factor
Gravitational acceleration
- H
Piezometric head
Creep-compliance of the spring of Kelvin–Voigt
-element
Absolute roughness of pipe
Cells number in the computational domain
Pressure
Volume discharge
Mass discharge
Pipe radius
Reynolds number
- S
Vector source term
- t
Time
- U
Vector variable
Longitudinal fluid velocity
- x
Unit vector in the x-direction
Calculation time step
Maximum tolerable calculation time step
Maximum tolerable calculation time step for stability of the discretized source term
Cell size in x-direction
Stress
Total strain
Elastic strain
Retarded viscoelastic strain
Constant eigenvalue of the Jacobian matrix
Mass of fluid per unit length of pipe
Mass density
Retardation time of the dashpot of
-element- ν
Poisson's ratio
- νr
Transverse velocity
- νR
Radial velocity
Solution obtained at n after the solution of the conservation part in the
-direction
Average value between n and
Sub-script
INTRODUCTION
In the detailed design of pipeline systems, water-hammer analysis is very important to select pipe characteristics and to specify surge control to prevent and/or mitigate any excessive pressures. The detailed design further attempts to elaborate the complete system response using various tools such as numerical modeling. Adamkowski & Lewandowski (2006) and Shamloo et al. (2015) reviewed a quasi-steady and four unsteady friction formulations (Zielke 1968; Trikha 1975; Brunone et al. 1991; Vardy & Brown 2003) for transient pipe flow. The Trikha (1975) model gives the Zielke (1968) model with assumptions that simplify the scheme. Bousso & Fuamba (2013) have numerically represented the unsteady friction in transient two-phase flow. Other unsteady friction alternatives have been proposed by other authors (Zidouh et al. 2009; Weinerowska-Bords 2015). Their contribution bridges the gap between the detailed phase and the general concept but cannot be used as a full purpose tool for transient analysis since certain key parameters such as the effects of the viscoelastic pipe-wall are not incorporated.
Bergant et al. (2008a, 2008b) have experimentally studied key criteria that influence the pressure profile predicted by the water-hammer phenomenon. In general, numerical simulation tools assume that the friction factor is considered steady, quasi-steady, or dynamic; this has the tendency to overestimate the water-hammer peak pressures since the primary mechanism that may significantly affect pressure waveforms is the viscoelastic of the non-steel pipe-wall. In materials with viscous and elastic properties, the mechanical damping of the pressure response is larger than the damping phenomenon due to the fluid friction during the event of transient occurrences (Gally et al. 1979; Rieutord & Blanchard 1979; Rieutord 1982; Franke 1983; Güney 1983; Suo & Wylie 1990; Covas et al. 2004, 2005). Ferrante & Capponi (2018) have compared viscoelastic models with a different number of parameters for transient simulations.
The Godunov method is a conservative numerical scheme for solving PDE in CFD. This well-known conservative finite-volume approach computes a relative or exact Riemann solver at each cell. It is first-order accuracy in both time and space. This framework has been used in several papers to design first-order quasi-linear numerical methods in one-dimensional non-conservative hyperbolic systems, in the finite volume framework; see for example, Muñoz-Ruiz & Parés (2007), Guinot (2002a, 2002b, 2004), Bousso & Fuamba (2013), and Seck et al. (2017). Guinot (2003) focused on the scheme of the Godunov method to solve the classical water-hammer equations using finite volumes with steady friction. Pal et al. (2020) compared the first- and second-order FVMs based on the Godunov scheme with the method of characteristics (MOCs) and experimental water-hammer measurements available in the literature. Seck et al. (2017) developed the derivation of water-hammer formulations’ hyperbolic systems of conservation laws with unsteady friction, which are not applicable for viscoelastic pipes because their non-negligible effects have not been taken into account.
Derived from the principles of conservation of momentum and mass, the water-hammer equations integrating unsteady friction and viscoelastic pipe-wall are expressed in hyperbolic systems of conservation laws with conserved quantities as variables. The presence of the viscoelastic pipe-wall terms preclude presentation in the integral (conservation) form and cause a source term to appear. The 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.
The innovative components of this contribution are to present the water-hammer derivation equations with the effects of viscoelastic pipe-wall and unsteady friction, and adapt the finite volume approach with conserved variables, for hyperbolic equations of conservation laws. The aim of this contribution is to provide a new user-friendly complete numerical algorithm in finite volume for studying several variants or design water-hammer, in a short period, of the series of different characteristics of viscoelastic pipes in the distribution network that can be of benefit to engineers.
METHODOLOGY
This paper focuses on implementation of the Godunov approach to the hyperbolic systems in one-dimensional conservation laws that describe the phenomenon of water-hammer integrating the viscoelastic pipe-wall and unsteady friction. The solution of a hyperbolic system of conservation laws by Godunov-type algorithms comprises six steps:
- (i)
discretization of space in finite volumes,
- (ii)
construction of generalised Riemann problems (GRP) at the cell,
- (iii)
conversion of the GRP into equivalent Riemann problems (ERP),
- (iv)
solution of the ERP,
- (v)
balance over the cells,
- (vi)
incorporation of the effect of the source terms.
The computation of the fluxes is detailed followed by a description of the complete numerical algorithm. A case study is examined with comparisons being made between the numerical solution with quasi-steady friction, the numerical solution with viscoelastic pipe-wall, a classical solution, and the experimental results.
Derivation from the conservation of mass and momentum
is expressed as:
, 
is the flow conserved variable vector, 
the flux vector, 
the source term vector, D the inner pipe diameter (m), 
the Vardy's shear decay constant related to the flow regime, x the unit vector in x-direction, t the time (s), 
the mass of liquid per unit length of pipe (kg/m), 
the mass discharge (kg/s), 
the fluid density (kg/m3), 
the cross-sectional pipe area (m2), V the longitudinal fluid velocity (m/s), 
wave speed (m/s), 
the steady state friction coefficient, and p the pressure (Pa) where the formula will be given later (Equation (21)).
In the definitions and basic notion of conservation laws, the source term 
can be a function of time t, space coordinate x and flow variable vector U, combined or separately. In this form (Equation (1)), where only the conservation part remains on the left, the term 
on the right is considered as being part of the source term. This source term thus contains an additional term 
. 
is referred to as a hyperbolic source term, if it incorporates the derivative 
. Systems that can be written under the form of Equation (1) are also called systems of conservation laws with source term, or hyperbolic systems of conservation laws, or balance laws (Parés 2006). The incorporation of the unsteady friction results in the appearance of the quasilinear term 
in the hyperbolic system (Equation (1)) thus destroying the classical integral form of the equations (Seck et al. 2018).
in Equation (1) may be obtained, respectively, by the Colebrook–White formula and the Hagen–Poiseuille formula which are well-known and documented worldwide. If the flow is turbulent, the momentum part of the source term contains an implicit expression 
for the steady friction factor. Clamond (2009) has proposed a robust, accurate, and fast numerical scheme for calculating the Colebrook–White formula that is more efficient than the result regarding the Lambert W-function, also sometimes called the Omega (Corless et al. 1996, 1997) or simple approach such as the Haaland (1983) equation:
is equivalent absolute roughness 
divided by the pipe diameter D, 
and 
the Reynolds number. The 
-function is linked to the Lambert 
-function as:
-function has been replaced with Lambert 
-function for the sake of simplicity.The viscoelastic pipe-wall inclusion
is divided into the sum of the retarded strain 
and elastic strain 
(Covas et al. 2004, 2005; Weinerowska-Bords 2006; Bergant et al.2008a, 2008b):
the fluid bulk modulus of elasticity (Pa), 
the pipe material Young's modulus (Pa), and 
is the constant that takes into account constraints and pipe cross-section dimensions; it is defined as (Streeter & Wylie 1993) for a thin-wall pipe fixed along its length.
The viscoelastic influences only the continuity equation (Covas et al. 2004, 2005; Weinerowska-Bords 2006; Bergant et al.2008a, 2008b).
is the transverse velocity (m/s) and r the radial coordinate. By integrating that continuity Equation (9) along the pipe cross-sectional area,
is the radial velocity at pipe-wall (m/s), R the pipe radius. Rewritten as:
and pipe radial expansion is given by (Güney 1983; Duan et al. 2010; Chaudhry 2014):
From a physical point of view, the third term in Equation (16) is the flux in the radial direction. That component is a function of conserved vector variable 
, and time derivative of retarded strain 
, combined, but does not contain functions of any of derivatives of 
.
contained in the source term can be neglected. Equation (18) may be rewritten as:
in the mass conservation part of the source term 
. To close the vector Equation (19), the wave speed 
that also relate the variations of 
to those of the pressure p is introduced as:
is the (constant) Jacobian matrix of the flux vector in x-direction with respect to the flow conserved variable vector. Since the homogeneous part of the system is linear, the two (constant) eigenvalues 
and 
of 
are:
The next section describes the solution process for the system.
Numerical solution
:
The first-order time splitting is used. The viscoelastic model (based on the generalized Kelvin–Voigt application) in the mass conservation part of the source term is solved using the modified Covas et al. (2005) numerical scheme. For the momentum part of the source term, the Colebrook–White equation is solved using a robust, fast, and accurate numerical scheme due to Clamond (2009).
Solutions at the internal cells and at the boundaries
of Equation (19) is Equation (24). It is the same conservation component as in classical Guinot's (2003) equation. For the internal cells, Guinot (2003) constructs the Riemann problem 
and presents the solution 
of this Riemann problem in Equation (26) computing the flux 
between the time intervals 
and 
according to Equation (27).
according to Equation (21); 
is a reference pressure at which the density 
is known. 
is calculated as 
.
boundary at the left side 
and the right side 
are, respectively, given by:
represents a pressure to be determined at the boundary. 
(at the left) and 
(at the right) is obtained from the prescribed pressure 
using:
at the left and at the right side boundary are, respectively, given by:
in Equation (31) and 
in Equation (32) are obtained, respectively, from:
Discretization of the source term
obtained at the end of the Godunov step as a starting point, the final solution 
of the source term at the end of the time step is presented by Toro (2001), Guinot (2003), Bousso & Fuamba (2013), and Seck et al. (2017) in the following form:
is the source term value, determined by using 
.
at time t must be known to calculate the time derivative of the retarded strain 
. The generalized Kelvin–Voigt model for a material with viscous and elastic properties, i.e., Equation (41) represented in Figure 1, is generally needed to characterize the function 
of the creep-compliance (Aklonis et al. 1972; Gally et al. 1979; Suo & Wylie 1990; Covas et al. 2005):
is the creep-compliance of the initial spring (Pa−1), 
the creep-compliance of the spring (Pa−1) of Kelvin–Voigt 
-element, 
the modulus of elasticity (Pa) of the spring of 
-element, 
the retardation time (s) of the dashpot of 
-element.
Generalized Kelvin–Voigt model of viscoelastic behavior of pipe hoop.
and 
were calibrated with the creep-compliance test details. The retarded viscoelastic strain 
is computed as the addition of these parameters of individual Kelvin–Voigt element k as:
and 
, respectively, are defined by:
is the gravitational acceleration and H the piezometric head (m). If the pipe material is isotropic and homogeneous, the Poisson's ratio is constant and viscoelastic for small strains is linear; function 
is equivalent to the stress 
(Pa) for the pipe subjected to the internal pressure expressed as (Streeter & Wylie 1993):
.Computational time step
(Guinot 2003) is as follows:
is needed. It is given by:
Resolution algorithm
The user-friendly proposed algorithm is presented in Figure 2.
Schematic of the proposed algorithm.
RESULTS AND DISCUSSION
A case study has been simulated. Data from the experimental test conducted by Covas et al. (2004) have been compared with the numerical results.
Experimental setup
The results of experimental studies conducted by Covas et al. (2004) obtained from a PE pipe-rig at Imperial College and considered to be of good quality are compared with the mathematical model. The experimental setup is composed of a 277 m long PE pipe with a thickness of 6.25 mm and an inner diameter of 50.6 mm that connects to an upstream air vessel. Initial steady flow (Q0 = 1.008 l/s) are followed by a transient occurrence initiated by the sudden closure of a globe valve at the downstream end. 
= 1,000 kg/m3 and g = 9.81 m/s2. The initial wave celerity a0 was evaluated at 395 m/s. The pressure waves were computed and compared with the corresponding experimental values at location #1 (distance x= 273 m from upstream), location #5 (distance x= 117.5 m from upstream) and location #8 (distance x= 199 m from upstream). The values of the creep coefficients 
for the retardation times 
are specified as follows. A full description of the collected data and the experimental system may be found in Covas et al. (2004).
J0 = 0.7.10−9 Pa−1
J1 = 0.1394.10−9 Pa−1;
= 0.05 s;J2 = 0.0062.10−9 Pa−1;
= 0.5 s;J3 = 0.1148.10−9Pa−1;
= 1.5 s;J4 = 0.3425.10−9Pa−1;
= 5 s;
= 10 s;
Results and discussion
The water-hammer equations incorporating viscoelastic pipe-wall were solved using the Godunov scheme detailed above. In keeping with the objectives of this paper, results were obtained for the analytical solution, the numerical solution with quasi-steady friction, and the numerical solution including the tuned viscoelastic model using a time step of 
. The experimental measurements are compared with these results. Good accord between measured and computed viscoelastic profiles has been obtained (Figures 3–5).
Pressure profiles (Δt = 0.0008 s, N = 100 cells) at location #1.
Pressure profiles (Δt = 0.0008 s, N = 100 cells) at location #8.
Pressure profiles (Δt = 0.0008 s, N = 100 cells) at location #5.
Comparison between the results indicates that the inclusion of viscoelasticity has a serious impact on the pressure profiles. The quasi-steady friction model does not reproduce the progression of the aspect of the pressure waves very well. They also display phase shift errors. The quasi-steady friction model undervalues the propagation and damping predicted by the (physically) more accurate viscoelastic model which tends to attenuate the magnitude of the overpressures more rapidly. Due to damping, the waves will be of decreasing amplitude until the final equilibrium pressure is reached. In the present context given the simple geometry and non-varying fluid properties (the fluid is slightly compressible), the basic Godunov scheme is appropriate for the analysis of water-hammer phenomenon. The wave profile of the viscoelastic line is smooth and the wave fronts are smeared which justify the presence of a rarefaction wave. That front smearing is caused by the exponential term contained in the mass conservation equation (Equation (40)). The Godunov scheme introduces a non-negligible diffusion. This explains the use of a large number of calculation cells to obtain a quality solution similar to that of a higher order scheme that requires few computing cells. More cells and therefore more computational effort is required for one-order scheme. But, for a robust and fast one-order scheme that runs in a short time with many cells, the use of a higher order scheme is not justified, because it will only improve the solution just weakly.
CONCLUSIONS
The water-hammer equations, in conservation form, combining viscoelastic pipe-wall and unsteady friction, are developed and numerically solved using a finite volume formulation for engineering applications. A robust, simple, accurate, and fast computational algorithm based on the Godunov scheme for one-dimensional hyperbolic systems is presented in some detail. Incorporation of the viscoelastic pipe-wall results in the appearance of a new term in the mass conservation part of the source term in the hyperbolic system of governing partial differential equations. From a physical point of view, this source term in mass conservation part is the flux in the radial direction. A case study has been presented to compare the separate impacts on wave attenuation with time. The findings indicate, as expected, inclusion of viscoelasticity reduces the peak water-hammer pressures. The model developed appears to be able to satisfactorily predict transient pressures in viscoelastic pipes and can be a useful tool for qualitative analysis.
These equations, as expected, admit ‘weak’ solutions with hydraulic jumps being the discontinuities. Using local first-order approximations to handle the spurious behavior near shocks does work without the expense of additional computational time. Computational time for the case studied took only 1 minute and 28 secs on a PC. Thus, for engineering applications, several variants for the design of series of different characteristic viscoelastic pipes in distribution networks may be studied in a short period of time with a spatial discretization not necessarily uniform.
The algorithm introduced in this paper could be used however, of course with some modifications, in the high-order schemes. This first-order proposed algorithm, in the context of hydraulic transients in pressurised pipelines incorporating viscoelastic pipe-wall and unsteady friction, is an interesting start to achieve this aim.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
%20cropped.png?versionId=7142)