## 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*-directionQuasi-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*-directionCalculation time step

Maximum tolerable calculation time step

Maximum tolerable calculation time step for stability of the discretized source term

Cell size in

*x*-directionStress

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

*et al.*2017) are:where 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/m^{3}), the cross-sectional pipe area (m^{2}), *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).

*W*-function, also sometimes called the Omega (Corless

*et al.*1996, 1997) or simple approach such as the Haaland (1983) equation:where is equivalent absolute roughness divided by the pipe diameter

*D*, and the Reynolds number. The -function is linked to the Lambert -function as:

### The viscoelastic pipe-wall inclusion

*et al.*2004, 2005; Weinerowska-Bords 2006; Bergant

*et al.*2008a, 2008b):

*e*is the pipe wall width, 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.in which,

*ν*is the Poisson's ratio.

The viscoelastic influences only the continuity equation (Covas *et al.* 2004, 2005; Weinerowska-Bords 2006; Bergant *et al.*2008a, 2008b).

*r*the radial coordinate. By integrating that continuity Equation (9) along the pipe cross-sectional area,where is the radial velocity at pipe-wall (m/s),

*R*the pipe radius. Rewritten as:

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

*et al.*1979; Rieutord & Blanchard 1979; Rieutord 1982; Franke 1983; Güney 1983; Suo & Wylie 1990; Covas

*et al.*2004, 2005). In fact, the primary mechanism that may seriously influence pressure waveforms is the viscoelastic of the non-steel pipe-wall (Bergant

*et al.*2008a, 2008b). Then, specifically for viscoelastic materials, the additional term quasilinear term contained in the source term can be neglected. Equation (18) may be rewritten as:

*p*is introduced as:where it is assumed that the wave speed depends on fluid compressibility, pipe physical characteristics, and external constraints.

*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

*et al.*(2017), the first step uses the Guinot (2003) solution for the mass conservation part of the PDE omitting the right-hand part :The second step uses Seck

*et al.*(2017), Toro (2001), Guinot (2003), and Bousso & Fuamba (2013) treatment of the source term in (Equation (19)),

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

### Discretization of the source term

*et al.*(2017) in the following form:where is the source term value, determined by using .

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

*k*, Covas

*et al.*(2005) proposed Equations (43)–(45), which are numerical approximations:where the functions and , respectively, are defined by:where 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

### Resolution algorithm

The user-friendly proposed algorithm is presented in Figure 2.

## 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 (*Q _{0}* = 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/m

^{3}and g = 9.81 m/s

^{2}. The initial wave celerity

*a*was evaluated at 395 m/s. The pressure waves were computed and compared with the corresponding experimental values at location #1 (distance

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

J

_{0}= 0.7.10^{−9}Pa^{−1}J

_{1}= 0.1394.10^{−9}Pa^{−1}; = 0.05 s;J

_{2}= 0.0062.10^{−9}Pa^{−1}; = 0.5 s;J

_{3}= 0.1148.10^{−9}Pa^{−1}; = 1.5 s;J

_{4}= 0.3425.10^{−9}Pa^{−1}; = 5 s;

_{5}= 0.0928.10

^{−9}Pa

^{−1}; = 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).

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.