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
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, hf, 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 hf is decomposed into two terms.
where Ai and Bi = Vardy's parameters that depend on f×Re; U= average velocity in each cross-section; D= pipe inner diameter; ν = kinematic fluid viscosity.
where k3= empirical decay coefficient; U= average velocity in each cross-section; SGN= sign operator.
Linear-viscoelastic model
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
EXPERIMENTAL SET-UP
(a) Schematic configuration of the experimental facility and (b) example of collected data.
(a) Schematic configuration of the experimental facility and (b) example of collected data.
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 (Ed), which is approximately 1.5–2.0 times the static modulus of elasticity (E0) 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, J0, and then the modulus of pipe elasticity, E0, must be corrected, which originates the term dynamic modulus of elasticity, Ed. For non-plastic pipes, E0 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 (NKV) and of the relaxation times τk, calibration of the Jk 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 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, τk, can be estimated for each K-V element as follows: (i) the first element is equal to half of the valve closure time, tc, (τ1 = tc/2); (ii) the second element is equal to half the period of the pressure wave, 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 Jk (for the pre-set of τk) should be calibrated simultaneously with the wave speed, 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 Jk parameters for the range of celerity values should be done simultaneously with UF parameters, if any (Step 6).
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)
(a) Coefficient k3 variation with flow rate for different wave speeds; (b) OF variation with wave speed; (c) comparison of pressure head variation between numerical and experimental data for Q = 2.726 l/s and k3 = 0.238; (d) absolute error between numerical and experimental values.
(a) Coefficient k3 variation with flow rate for different wave speeds; (b) OF variation with wave speed; (c) comparison of pressure head variation between numerical and experimental data for Q = 2.726 l/s and k3 = 0.238; (d) absolute error between numerical and experimental values.
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 Jk (for a pre-set of τk) should be calibrated simultaneously with the wave speed, a.
(a) Creep function and OF variations with wave speed for Q = 2.726 l/s; (b) comparison of creep functions neglecting UF loss model (solid line) with those considering UF loss model (dashed line); (c) creep coefficients' variation with wave speed; (d) comparison of piezometric head variation between numerical and experimental data; (e) absolute error between numerical and experimental values.
(a) Creep function and OF variations with wave speed for Q = 2.726 l/s; (b) comparison of creep functions neglecting UF loss model (solid line) with those considering UF loss model (dashed line); (c) creep coefficients' variation with wave speed; (d) comparison of piezometric head variation between numerical and experimental data; (e) absolute error between numerical and experimental values.
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 k3, VE coefficients τk and Jk) 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.
(a) Artificial data vs. the best fitted solution and (b) average sensitivity of piezometric head with each parameter xi for the best solution (unsuccessful calibration of four different parameters).
(a) Artificial data vs. the best fitted solution and (b) average sensitivity of piezometric head with each parameter xi for the best solution (unsuccessful calibration of four different parameters).
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.