The current paper aims at the discussion of the uncertainties of hydraulic transient modelling using the classic method of characteristics and incorporating unsteady friction (UF) and pipe-wall viscoelasticity. A two-stage calibration procedure is used for estimating the model parameters: the first stage refers to the calibration of steady state conditions and the second stage to the transient state. This step-wise calibration procedure, although simple, consists of a set of guidelines useful for engineering practice when having pressure and flow data, instead of the blind simultaneous calibration of all known parameters. A data collection programme was carried out in an experimental facility, composed of high-density polyethylene pipes to collect transient pressure data for different initial discharges. Results obtained have shown that UF and viscoelasticity, despite being phenomena that have completely different natures, have similar effects in the transient pressure response, which the inverse transient solver is not capable of distinguishing.

## INTRODUCTION

The classic theory of water hammer is commonly used for design, as it reasonably well describes the maximum and minimum pressure variations. This approach typically assumes elastic behaviour of the pipe wall, steady state friction losses, one-phase flow, completely constrained pipe from axial movement, and no lateral in/out flows from the pipe. These assumptions are not always verified in practice. Examples of these are the energy dissipation during fast transient events, the mechanical damping in plastic pipes, and two-phase flow.

Several theoretical and practical studies have been carried out on hydraulic transients. Approaches which take into account one, or more, phenomena that affect transients have been developed, such as unsteady friction (UF) (Adamkowski & Lewandowski 2006; Duan *et al.* 2010; Storli & Nielsen 2011a; Mitosek & Szymkiewicz 2012) and viscoelasticity (Covas *et al.* 2005; Soares *et al.* 2008; Duan *et al.* 2010).

The use of mathematical models in hydraulic transient analysis helps to better understand the behaviour of a system as well as to predict extreme pressures under different operational conditions. The analysis of hydraulic transients in pressurized systems has become a common practice in engineering due to a general awareness of the economic losses and operational disturbances caused by transient events. Currently, there are new challenges for the development of more accurate and easy-to-use computational models to predict the maximum information on the behaviour of a hydraulic system (Ghidaoui *et al.* 2005).

The objective of this study is two-fold. First, it aims at showing and analysing transient data in an experimental facility composed of high-density polyethylene (HDPE) pipes. A data collection programme has been carried out in an experimental facility assembled at Instituto Superior Técnico of Lisbon, Portugal. Although Soares *et al.* (2012) have studied the dynamic effects of both cavitation and pipe-wall viscoelasticity in the same pipe-rig, the analysis of UF and pipe-wall viscoelasticity simultaneously has not been done. Second, it aims at showing the complexity of the hydraulic transient solver (TS) calibration when dealing with water hammer modelling incorporating UF and pipe-wall viscoelasticity. Although Duan *et al.* (2010) have quantitatively shown that in plastic pipes the role of UF is relevant only in the first phases of the transients, and that viscoelastic effects become increasingly more dominant over UF as time progresses, both effects have hardly been distinguished by inverse transient analysis. Results obtained in this study have shown that UF and viscoelasticity, despite being phenomena with a completely different nature, have similar effects in the transient pressure response, which the inverse transient solver (ITS) is not capable of distinguishing.

## UNCERTAINTIES IN TRANSIENT HYDRAULICS

An uncertainty can be defined as a lack of knowledge to deterministically or numerically describe or predict system behaviour or its characteristics. In transient hydraulics, examples of the unknown parameters are pipe roughness, wave speed, UF coefficients, creep coefficients, air cavity volumes or boundary condition characteristics (e.g., valve manoeuvre or pump rated conditions). Within the hydraulics literature, the parameter identification methods have tended to focus on the estimation of pipe friction parameters and of leak sizes (Zecchin *et al.* 2013). Many of these have focused on customized approaches for single pipe systems with either measured (Verde *et al.* 2007) or known boundary conditions (Wang *et al.* 2002). Calibration can be carried out by trial-and-error procedure or by optimization algorithms. The main difficulties in the calibration procedures are associated with uncertainties of different natures, namely topological and physical, hydraulic, data collection and the calibration method.

*Topological and physical* uncertainties include pipe characteristics (materials, lengths, diameters and axial constraints), system configuration (branches, valves, dead-ends, boundary conditions), air pockets or dissolved air. Discontinuities, like pipe diameter/material changes, tee-junctions, local losses or air pockets, introduce changes in the transient pressure signal: a leak creates a pressure drop, a dead-end or a closed-valve reflects an incident wave, and an air bubble creates a pressure drop followed by a pressure increase. These signals can be misleading, unless the topology of the system is extremely well known and each can be clearly identified in the signal. Another example is the reduction of the internal pipe diameter due to corrosion, tuberculation and scaling occurring (Jung & Karney 2008); a 10% decrease in the pipe diameter can increase steady state friction losses by 40% (Walski *et al.* 2001).

*Hydraulic uncertainties* are associated with the flow characteristics, the type of valve manoeuvre and wave speed/pipe creep. Uncertainties in the pipe characteristics and non-elastic behaviour of pipe walls have serious consequences in wave speed estimation. *Wave speed* is another challenging uncertainty in a pipeline. It is a function of fluid (water density, bulk coefficient, temperature, presence of air and solids) and pipe properties (pipe diameter, thickness and material, pipe constraints). Some of these conditions are poorly defined and uncertain (Jung & Karney 2008). For example, the accurate measurement of the free-air content in fluid is difficult; however, even a tiny amount of gas in liquid greatly reduces the propagation velocity of a pressure wave in a pipeline (Wylie & Streeter 1993). A roughly estimated wave speed can lead to the unsuccessful calibration of other parameters. A possible solution is the installation of pressure transducers at different locations and wave speed estimation based on ‘the transient travelling time between transducers’.

Examples of *data collection uncertainties* are pressure transducer location, number of sensors, acquisition rate, sensor accuracy, and signal noise. The pressure sensor needs to be next to the site where the transient is generated to facilitate calibration. Transients tend to be dissipated during travel in the pipe. Having several pressure sensors along the pipe is important to identify the location of multiple features and to better estimate wave speed. The sensor accuracy and the noise are also important to distinguish the pressure drop in the signal.

One of the most widely used methods to estimate the unknown parameters is the *inverse transient method*, also called inverse solvers. In recent years, there has been a significant interest in the application of the inverse transient approach for leak detection and calibration of water pipe system models (Liggett & Chen 1994; Nash & Karney 1999; Vitkovský *et al.* 2000b; Kapelan *et al.* 2003; Covas & Ramos 2010). Even though inverse transient techniques have been widely investigated, many challenges still remain.

## MATHEMATICAL MODEL

### Basic equations

*Q*= discharge;

*H*= piezometric head;

*g*

*=*gravity acceleration;

*S*

*=*pipe cross-sectional area;

*x*

*=*coordinate along the pipeline axis;

*t*

*=*time;

*a*

*=*elastic wave speed;

*h*= slope of the energy line.

_{f}The derivation of the governing equations of one-dimensional transient flow in pressurized conduits takes into account simplifying assumptions such as (i) a pseudo-uniform velocity profile, (ii) elastic rheological behaviour of the pipe wall, (iii) fluid being one-phase, homogenous and quasi-incompressible and (iv) the pipe being uniform and constrained. These equations are quasi-linear and hyperbolic, valid at intermediate sections of the pipe.

### Unsteady friction

In classical water hammer theory, the friction term, *h _{f}*, is calculated in the same manner as in the steady state regime. This approach is satisfactory for slow transients, in which velocity profiles do not significantly differ from the steady state ones (Zidouh 2009). For rapidly varying flows or higher pulsating frequencies, these approximations are inaccurate in the description of the damping and dispersion of the transient pressure wave (Covas 2003). In order to take into account UF losses and fluid inertial effects, the head loss per unit length

*h*is decomposed into two terms.

_{f}*Y*

_{1}and

*Y*

_{2}) developed for smooth turbulent flows:

where *A _{i}* and

*B*= Vardy's parameters that depend on

_{i}*f*

*×*

*Re*;

*U*

*=*average velocity in each cross-section;

*D*

*=*pipe inner diameter; ν = kinematic fluid viscosity.

*et al.*1991; Pezzinga 2000; Bergant

*et al.*2001; Ghidaoui

*et al.*2005; Storli & Nielsen 2011a, 2011b). The most widely used formulation is that of Vitkovský

*et al.*(2000a), which is a modification of Brunone's model:

where *k _{3}*

*=*empirical decay coefficient;

*U*

*=*average velocity in each cross-section;

*SGN*

*=*sign operator.

### Linear-viscoelastic model

*J(t)*. Decomposing total strain in the sum of an elastic strain, ɛ

*, and a retarded strain, ɛ*

_{e}*, the continuity equation yields:*

_{r}The third term of Equation (5) represents the retarded strain, and the instantaneous-elastic strain is included in the piezometric-head time derivative and in the elastic wave speed, *a*.

### Method of characteristics

*C*and

^{+}*C*), which can be solved after being expressed in terms of finite differences. The friction and the piezometric-head terms may be included in the equations and can handle complex systems of pipes. This method is not the only one that can be used to solve the equations, but it is increasingly used for the study of water hammer problems. The result of MOC's application to Equations (1) and (5) is the following set of equations:

^{−}*h*, and the retarded strain time derivative, ∂ɛ

_{f}*∂*

_{r}/*t*. Several mechanical models can be used, combining springs and dashpots connected in series or in parallel to numerically describe the viscoelastic behaviour of materials. PE pipes are viscoelastic solids, being described by the generalized Kelvin–Voigt model: where

*J*

_{0}= creep-compliance of the first spring (

*J*

_{0}= 1/

*E*

_{0}, with

*E*

_{0}as Young's modulus of pipe elasticity),

*J*= creep-compliance of the spring of the Kelvin–Voigt

_{k}*k*-element (

*J*= 1/

_{k}*E*, with

_{k}*E*as the modulus of elasticity of the

_{k}*k*-th spring element), τ

*=*

_{k}*μ*/

_{k}*E*the retardation time of the dashpot of element

_{k}*k*, with

*μ*as the corresponding viscosity, and

_{k}*N*= number of Kelvin–Voigt elements.

_{KV}## EXPERIMENTAL SET-UP

^{2}. The pipe was installed in two coils with 1 m of radius of curvature each. The HDPE pipe has a total length of 199 m. The installation has a total length of 203.37 m from the hydropneumatic vessel to the downstream valve. The downstream boundary condition of the experimental facility is an atmosphere-valve. This valve is a ball valve type, and is manually operated and closed and opened as fast as possible to simulate quasi-instantaneous valve manoeuvres. The discharge is controlled by another ball valve at the upstream of the tank. The configuration of the experimental facility and an example of the collected data set at the downstream end pressure transducer (T3) are depicted in Figure 1.

## CALIBRATION PROCEDURE

### Steady-state calibration (Steps 1–3)

The calibration of the steady state flow consists of the adjustment of the calculated piezometric heads to the measured values. This calibration depends on the initial discharge and, consequently, on the flow resistance equations. For example, for smooth pipes with high Reynolds numbers (Re >4,000), the von Kármán–Prandtl formulation may be used; for rough pipes and high Reynolds numbers, the Colebrook–White formula may be used; for the laminar regime, the Hagen–Poiseuille equation should be used (Step 1). If the selected flow resistance equation has parameters (e.g., pipe roughness), then these have to be calibrated (in Step 2). Usually, the continuous and the local head losses are incorporated in the equivalent pipe roughness. The last step of calibration of the steady state calibration (Step 3) focuses on the valve opening. The calibration of this parameter can be carried out by inverting the equations to obtain a reference value, or by a trial-and-error procedure.

### Transient state calibration (Steps 4–8)

#### Valve manoeuvre and wave speed (Step 4)

The calibration of the valve manoeuvre requires (i) the initial time of closure based on a preliminary analysis of measurements, (ii) the total time of closure (corresponding to the inflection point of the first pressure wave) and (iii) the shape of the manoeuvre's curve, which can be defined by scrutinizing the water hammer overpressure due to the valve closure.

#### Range of wave speed values (Step 5)

A range of wave speed values should be estimated using maximum and minimum values of dynamic modulus of elasticity, Poisson coefficient and axial constraints. Wave speed can be estimated by theoretical formulae (Korteweg expression). For plastic pipes the dynamic modulus of elasticity should be used (*E _{d}*), which is approximately 1.5–2.0 times the static modulus of elasticity (

*E*

_{0}) provided by manufacturers of HDPE pipes (Covas 2003), and 1.1–1.2 times for PVC pipes (Soares

*et al.*2008). This is because the first pressure peak due to valve closure depends on the creep-compliance of the first spring of the Kelvin–Voigt model,

*J*

_{0}, and then the modulus of pipe elasticity,

*E*

_{0}, must be corrected, which originates the term dynamic modulus of elasticity,

*E*. For non-plastic pipes,

_{d}*E*

_{0}can be used.

#### Using an elastic transient solver (Step 6)

The use of an elastic TS is more suitable for concrete and metal pipes since these pipe materials have an elastic rheological behaviour. The solver can incorporate either a steady state friction model or a UF model. In the former case, the wave speed is calibrated by a trial-and-error procedure or by using an optimization algorithm fitting the numerical results to the maximum and minimum observed pressures. In the latter case, the wave speed is calibrated together with the UF model parameters (e.g., decay coefficient) trying to fit both extreme pressures as well as the pressure wave damping and phase shift (Step 6).

#### Using a viscoelastic transient solver (Steps 6 and 7)

When the pipe is made of plastic, such as HDPE or PVC, a viscoelastic TS should be used. If a steady state friction loss model is used, the calibration model consists of determining the number of Kelvin–Voigt elements (*N _{KV}*) and of the relaxation times τ

*, calibration of the*

_{k}*J*parameters in a range of wave speed values, and selection of the best fitted solution. The number of Kelvin–Voigt (K-V) elements to consider in the mathematical model depends essentially on the type of pipe material. According to Covas (2003), for HDPE pipes three elements of the K-V model are required to obtain a good calibration, whereas for PVC pipes, Soares

_{k}*et al.*(2008) state that only one element is needed. Based on these studies, the K-V parameters, relaxation times and creep-compliance coefficients result in several combinations between them that fit the model results to measured data. In order to avoid undetermined solutions originated by the ITS, and based on the previous studies, relaxation times can be set by scrutinizing the pressure wave history observed. Relaxation times, τ

*, can be estimated for each K-V element as follows: (i) the first element is equal to half of the valve closure time,*

_{k}*t*, (τ

_{c}_{1}=

*t*/2); (ii) the second element is equal to half the period of the pressure wave,

_{c}*T*, (τ

_{2}=

*T*/2); and (iii) the third element is equal to one-third of the simulation time,

*t*, (τ

_{3}=

*t*/2) (Step 6). The creep function coefficient

*J*(for the pre-set of τ

_{k}*) should be calibrated simultaneously with the wave speed,*

_{k}*a*, by using an inverse transient method (Step 7). If instead of a steady state friction loss model an unsteady model is used then the calibration of the

*J*parameters for the range of celerity values should be done simultaneously with UF parameters, if any (Step 6).

_{k}## RESULTS

The mathematical model developed in this study has been coded in C ++ language. The TS, which is based on the MOC, was linked to the Levenberg–Marquardt algorithm in order to perform the step-wise calibration procedure presented herein. The inverse transient method obtained has been applied to calibrate the hydraulic model of the reservoir–pipe–valve system presented in Figure 1(a).

### Steady state calibration

In the developed mathematical model, an explicit formulation of the Colebrook–White equation was implemented. An estimation of the absolute equivalent roughness was done, obtaining the average value of 0.06 mm (Step 2). To fit the calculated piezometric head to the measured value, the valve opening was calibrated by a trial-and-error procedure (Step 3).

### Unsteady state calibration

In the current case, the valve manoeuvre (Step 4) was approximated to a tri-linear manoeuvre. The static modulus of elasticity for HDPE pipes varies between 0.7 and 1.0 GPa; accordingly, the estimated range of wave speed is between 273 and 371 m/s (Step 5). Calibration in the next steps was carried out for a wave speed range between 270 and 380 m/s with an increment of 10 m/s.

An elastic solver with UF modelling (UF solver), a viscoelastic solver with a steady state friction loss model (VE solver), and a viscoelastic solver with unsteady state friction loss model (VE + UF solver) were used. In the latter case, the UF formulation has parameters that have to be calibrated simultaneously with the wave speed (Step 6).

#### Elastic transient solver (UF solver)

*k*, is uncertain and needs to be calibrated simultaneously with the wave speed,

_{3}*a*. This calibration presupposes that previous Steps 1–4 have already been carried out. An ITS running the Levenberg–Marquardt optimization algorithm was used to calibrate

*k*for fixed values of

_{3}*a.*The procedure was repeated for each value within the range of expected wave speed values (i.e., from 270 to 390 m/s with a set of 10 by 10 m/s). The analysis of the results obtained (Figure 3) shows the following. (i) Although the optimum values occur for decreasing values of

*k*

_{3}coefficient (black line in Figure 3(a)) there is a tendency to almost linear increase of the

*k*

_{3}coefficient with increasing flow rate. This trend was expected since the higher the flow, the greater the forces of inertia and the flow resistance in unsteady flows, and hence the greater the value of the

*k*

_{3}coefficient. (ii) Theoretically, the minimum value of the objective function (OF) would match the optimum value of the empirical

*k*

_{3}coefficient that fits the numerical results to experimental data. For each test, optimum values (i.e., minimum values of the OF) occur for different wave speed values and, in this case, the empirical

*k*

_{3}coefficient increases as flow rate decreases (Figure 3(b)). (iii) OF values increase for each flow rate test, which would be expected since the higher the flow rate, the higher the pressure changes will be as well as the higher the differences between measured and calculated values. (iv) There is a reasonable fit between measured and calculated pressure peaks. However, Vitkovský

*et al.*’s (2000a) formulation smooths the pressure wave shape slightly but induces lower pressure damping than those observed in the experiments (Figures 3(c) and 3(d)).

#### Viscoelastic transient solver

The viscoelastic TS for describing hydraulic transients in plastic pipes was used with a steady state friction loss model (VE) as well as taking into account a UF loss model (VE + UF). For the sake of simplification, an unparameterized Vardy's (1992) UF formulation was used in the latter. For both cases, the creep function coefficients *J _{k}* (for a pre-set of τ

*) should be calibrated simultaneously with the wave speed,*

_{k}*a*.

*Q*= 0.499 l/s the component

_{0}*J*is null, which means that only two Kelvin–Voigt elements are required to represent the creep function for low flow rates in HDPE pipes (Figure 4(b)); (iii) whereas the value of

_{3}*J*rises almost linearly with the wave speed, and increasing values of

_{1}*J*and

_{2}*J*stabilize as the wave speed increases, which means that the coefficient

_{3}*J*

_{1}determined for τ

_{1}=

*t*/2 counterbalances the increase in the wave speed by delaying the pressure wave (Figure 4(c)) (i.e.,

_{c}*J*

_{1}delays the pressure wave whereas a higher wave speed diminishes the pressure wave period); (iv)

*J*

_{2}has a significant increase between wave speed values 270 and 280 m/s as well as small variations from 280 to 320 m/s; however, the determined creep coefficient

*J*

_{2}is practically constant for wave speed values higher than 320 m/s; (v) for wave speed values between 270 and 290 m/s, there is a marked decrease of

*J*

_{3}, whereas its value stabilizes around 320 m/s; (vi) numerical results fit to observed pressure heads for all tests and the error between them is negligible (Figures 4(d) and 4(e)); (vii) the difference between creep functions (VE and VE + UF) decreases as the discharge decreases (Figure 4(b)); (viii) the creep functions considering UF losses are always lower (around 2%) than those obtained neglecting UF. This is because the creep function describes artificially the viscoelasticity of the pipe wall, the energy dissipation due to UF losses (if they are ignored in the simulation), and other dissipative effects not taken into account by the hydraulic model (e.g., that the conduit is wound); however, higher creep functions have been obtained by implementing Vitkovský

*et al.*’s (2000a) model, since the formulation introduces a small artificial increase in the first pressure peak, which the calibrated creep function compensates by increasing its value.

## DISCUSSION

Inverse transient analysis attempts to estimate unknown parameters by using pressure data collected during the occurrence of simulated pressure surges. The parameter identification is an optimization problem in which the system's behaviour is simulated by a TS and the difference between observed and calculated variables is minimized by means of an optimization model – ITS. An ITS is an optimization algorithm that searches for the best-fitted solution by minimizing an OF defined by the average least-square errors between observed and calculated variables. Theoretically and ideally, when using an optimization algorithm all unknown parameters could be simultaneously calibrated instead of using a step-wise procedure as used here; however, there were two main reasons for the authors not having followed, or recommending, that approach.

The first reason is that when using real data, no matter how well calibrated the hydraulic solver is, the OF is never zero. First, there are always some uncertainties in the formulations used to describe considered phenomena (e.g., UF formulations, the linear viscoelasticity assumption); second, there are other dynamic effects that are neglected in the model (e.g., a small percentage of dissolved air); finally, the pressure signal (even after being filtered) always has inevitable noise resulting from other electrical interferences or mechanical vibrations not described by the mathematical model.

The second reason for not performing simultaneous calibration of all unknown parameters (i.e., valve closure time, valve closure manoeuvre, pipe-wall roughness, wave speed *a*, UF coefficient *k _{3}*, VE coefficients τ

*and*

_{k}*J*) is because the effect of these parameters overlaps and is not clearly distinguishable in terms of observed dissipation, delay and shape of transient pressure signal. Additionally, the inverse solver has different sensitivities to each parameter, and the result is a fast convergence to the parameters with higher sensitivities, and only afterwards to the other parameters. As a result, the final solution (i.e., the best fitted parameters) is one of many combinations of parameters that leads to minimum values of the OF (this solution can have an OF lower that the one corresponding to the true combination of parameters) as shown in Figures 3 and 4.

_{k}*k*and VE coefficients

_{3}*J*, as presented in Figure 5. The major difference in sensitivities is between the decay coefficient (|

_{k}*dH/dx*| ∼ 10) and the creep coefficients (|

_{i}*dH/dx*| ∼ 10

_{i}^{9}to 10

^{10}), whereas the other two parameters (leaks and pipe roughness) have intermediate sensitivities. This is one of the reasons why the ITS fails to accurately estimate the parameters. The second is that UF and VE, despite being phenomena with completely different natures, have similar effects in the transient pressure response, which the ITS is not capable of distinguishing. Similar results were obtained in this study (Figure 4), in which UF and pipe-wall viscoelasticity have similar effects in transient pressure response.

## CONCLUSIONS

This research work aimed at the discussion of the uncertainties of inverse transient analysis incorporating UF and pipe-wall viscoelasticity. A hydraulic TS incorporating UF formulations and two pipe-wall rheological behaviours was developed. The mathematical model was calibrated using data collected in a laboratory facility made of HDPE pipes. Due to the overlapping of dynamic effects, UF damping in the pressure transient signal is quite similar to the damping caused by the viscoelastic pipe-wall behaviour. Although recent studies have shown the relevance of both effects in transients, the results obtained in this study have shown that ITSs are not capable of distinguishing UF and viscoelasticity effects when their parameters are determined simultaneously. In this way, a step-wise calibration procedure of hydraulic TSs was suggested, instead of the blind simultaneous calibration of all known parameters with distinct uncertainties associated with water hammer modelling.