Uncertainties of inverse transient modelling with unsteady friction and pipe-wall viscoelasticity

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.

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.

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.

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_f is decomposed into two terms.

where A_i and B_i = Vardy's parameters that depend on f×Re; U= average velocity in each cross-section; D= pipe inner diameter; ν = kinematic fluid viscosity.

where k₃= empirical decay coefficient; U= average velocity in each cross-section; SGN= sign operator.

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.

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.

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.

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₀) 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₀, and then the modulus of pipe elasticity, E₀, must be corrected, which originates the term dynamic modulus of elasticity, E_d. For non-plastic pipes, E₀ can be used.

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

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 τ_k, calibration of the J_k 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, t_c, (τ₁ = t_c/2); (ii) the second element is equal to half the period of the pressure wave, T, (τ₂ = T/2); and (iii) the third element is equal to one-third of the simulation time, t, (τ₃ = t/2) (Step 6). The creep function coefficient J_k (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 J_k parameters for the range of celerity values should be done simultaneously with UF parameters, if any (Step 6).

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

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

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

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 τ_k) should be calibrated simultaneously with the wave speed, a.

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₃, VE coefficients τ_k and J_k) 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.

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.