Modelling water hammer in viscoelastic pipelines: short brief

The model of water hammer in viscoelastic pipelines is analyzed. An appropriate mathematical model of water hammer in polymer pipelines is presented. An additional term has been added to continuity equation to describe the retarded deformation of the pipe wall. The mechanical behavior of viscoelastic material is described by generalized Kelvin-Voigt model. The comparison of numerical simulation and experimental data from well known papers is presented. Short discussion about obtained results are given.


Introduction
Hydraulic transient analysis is important during the design of water supply systems. One of the most dangerous phenomena for pipeline systems is water hammer. Analysis of water hammer allows to choose an appropriate parameters of system like pipe material, wall thickness or surge protection devices. For many years, researches used only the classic water hammer model. The main assumption in this model is linear elastic behavior of pipe walls [1][2][3]. This theory is good for modelling systems build mainly from steel or brass.
In last 40 years, researchers started investigate a new water hammer model for pipes made from a new generation polymeric materials (see e.g. [7,10,11]). High resistant properties (mechanical, chemical, temperature) of polymers, like polyethylene or polyvinyl chloride, caused increasing popularity of this type of pipelines. These materials have viscoelastic behavior which influence on the pressure response of a pipe system.
Viscoelastic rheological behavior of polymers has been studied by many authors [4][5][6]. This behavior is characterized by a sum of instantaneous elastic strain and retarded strain. It makes a mechanical damping of the pressure wave. There are two main approaches to simulate the viscoelastic behavior of plastic pipe-walls.
The first approach assumes frequency-dependent wave speed. The viscoealstic behavior is modelled by a frequency-dependent creep function to calculate wave speed, instead of a constant elastic modulus [7][8][9][10].
The second, more popular, approach assumes a retarded viscoelastic term in continuity equation. The precursory paper were written by Gally et al. [11] and Rieutford and Blanchard [12]. Gally at al. presented the mathematical model and comparison of theoretical and experimental results. Covas et al. in [13] analyzed the dynamic effects of pipe wall viscoelasticity on hydraulic transients. This effects have been observed in data collected on the two experimental setups. In two-part paper [14,15], Covas with others widely described the viscoelasticity effects in polymer pipelines. The theoretical and experimental results were compared. The creep-function of the pipe material was experimentally determined. Bergant et al. [16,17] investigated parameters that may significantly affect water hammer wave attenuation, shape and timing. The authors considered five possible mechanism that affect pressure waveforms include unsteady friction, cavitation, fluid-structure interaction, viscoelastic behavior of the pipe-wall material, leakages and blockages. In [18], Soares et al. investigated the analysis of hydraulic transients in polyvinyl chloride pipes. The creep function of the PVC pipes was determined by using an inverse transient method. Different characterization of creep function of PVC in comparison with HDPE pipes was observed. Keramat et al. [19] considered water hammer model in viscoelastic pipes with Poisson's ratio as a function of time. K. Weinerowska-Bords discussed in [20] the main problems connected with applying the water hammer model in viscoelastic pipes. The aspects of estimation a group of parameters like wave speed, retardation time and creep compliance was analyzed.
The structure of the work is as follows. After Introduction, Section 2 presents mathematical model of water hammer in elastic pipes. Next, Section 3 contains efficient solution of convolution integral. Section 4 is devoted to present mathematical model of water hammer in viscoelastic pipelines. In Section 5, numerical simulation are conducted. In conclusions, the summary of obtained results is given.

Water hammer in elastic pipe -fundamental equations
Unsteady flow of liquid in elastic pipes (made of brass, steel etc.) is often represented by a two 1D hyperbolic partial differential equations. The momentum and continuity equations are [21][22][23][24][25][26][27]:  In unsteady flow in pipe, the instantaneous wall shear stress τ may be regarded as the sum of two components [28][29][30][31][32][33]: The first component in Equation (4) τ q presents the quasi-steady wall shear stress and the second τ u is the additional contribution due to unsteadiness. Equation (4) relates the wall shear stress to the instantaneous average velocity and to the weighted past velocity changes. Among methods, which enable to resolve system of above equations, particular attention should be paid to the method of characteristics (MOC), which perfectly interprets the essence of natural phenomena of transient flow, and at the same time is characterized by fast convergence, easiness to take into account various boundary conditions and the high accuracy of calculation results. With its help one can easily solve a system of partial differential equations of quasi-linear hyperbolic type (1)-(4).

Efficient solution of convolution integral
Zielke presented the weighting functions for laminar flow [28]:  (4)) can be conducted using the first-order differential approximation [26,28]: Recently Vardy-Brown [34] rightly pointed out an error in this methodology. To avoid it they suggest to calculate the τ u component using following equation: where ∆t -dimensionless time step [−]. Determination of the wall shear stress using above formulas ((6) or (7)) is very time consuming. Trikha [35] was the first who developed an effective method of solving the integral convolution. Next Kagawa et al. [36] and Schohl [37] have improved this method. But as all this solution fails to model correctly the wall shear stress, a new efficient solution were presented recently [38] that is fully consistent with corrected by Vardy-Brown (7) classic inefficient solution: where In above equation n i and m i are constants of efective weighting function and θ is a correction factor [38]. The new effective method (8) presented above requires that the weighting function have the form of finite sum of exponential expressions: The needed number of exponential terms that make up the final form have been discussed recently in some works [33,38]. And as it was pointed there for a case of water hammer two terms is just enough for proper simulation.
As they value depends on the dimensionless time step used in simulation -a proper weighting function terms n i and m i used in completed simulations will be presented with all other constants in section 5.

Water hammer in viscoelastic pipe -fundamental equations
Taking into account viscoelastic behavior of the pipe wall one can derive the continuity equation similarly as for the linear elastic pipe. The most important details of derivation are presented in the Appendix A.
Polymer pipeline do not respond according to Hook law when subjected to a certain instantaneous stress σ 0 , one notice an immediate-elastic response and a retarded-viscous response. Thus strain can be decomposed into a sum of instantaneous-elastic strain ε e and a retarded strain ε r (see e.g. [11,14,16,18]): The final form of continuity equation is very similar to the one used in elastic pipeline: For small strains, a combination of stresses that act independently in a system result in strains that can be added linearly. Thus, the total strain generated by a continuous application a stress σ(t) is [5,[13][14][15]: where J 0 -the instantaneous creep-compliance, J(u) -the creep function at u time. For all linear-elastic materials, creep compliance J 0 is assumed equal to the inverse modulus of elasticity For final solution we additionally assume: homogeneous and isotropic pipe material; linear viscoelastic behavior for small strains and constant Poisson's ratio ν (mechanical behavior is only dependent on a creep-function). Then the total circumferential strain ε = (D(t) − D 0 )/D 0 can be described by: where σ = α∆pD/2e -circumferential stress [P a], p(t) -pressure at time t, p 0 -initial pressure, D(t) and D 0 -inner diameters at time t and 0 respectively, e(t) and e 0 -the wall thicknesses at time t and 0 respectively, α(t) and α 0 -the pipe-wall constraint coefficient at time t and 0 respectively, J(t) -the creep-compliance function at time t.
The creep function which describes viscoelastic behaviour of the pipe material should be determined experimentally for polymer pipelines in independent mechanical tests. Using generalized Kelvin-Voigt model in this work (Fig. 1) viscoelastic solid are defined with creep function of the pipe wall in form [5] of finite sum: where J k -the creep-compliance of the spring of the Kelvin-Voigt k-element defined by  The set of differential equations (1) and (11) is solved by the Method of Characteristics (MOC). With the help of this method this system of equations can be transformed to pair of ordinary differential equations: Using for above equation finite differences method on characteristics grid (Fig. 2) one obtain: In above equations (16) where constants: Lets now all known terms calculated from previous numerical time step write as: In numerical procedure after determining these two values there is still need to determine the momentary retarded strain ε D,k which is used in the equation (17):

Simulation example
To verify results of simulation using elastic water hammer model and a modified viscoelastic one, they results have been compared to experimental results. The experiments have been carried out by D. Covas in London on a PE pipe-rig at Imperial College [14,45]. The experimental reservoir-pipeline Note that to get this coeficient a scalling procedure using an extended Vardy-Brown coefficient [39] and a universal weighting function assumption were made [26], also characteristic polymer pipe roughness size k s = 1.5 · 10 −6 were assumed) and m 1 = 4.4413, m 2 = 33.6391, n 1 = 72.591, n 2 = 2744.3 (for Re=1400 laminar case). Knowing Poisson ratio one can determine the pipe-wall constraint coefficient α. As there is at least four different equations: that can be used to determine α value for thick pipes that are anchored along its length. In the simulations carried out below examined will be not only presented models (elastic and viscoelastic) but also the impact of this parameter on the results (using constant c = 395m/swhich is equivalent with assuming α · J 0 = const).
The graphical representations of transient states are presented below on Fig. 3 and Fig. 4. From simulations presented in Fig. 3 and Fig. 4, we can formulate the following conclusions: (i) The linear elastic model (with unsteady friction factor) is not an appropriate model to simulate water hammer in pressurized pipes made from polymers. Even, the maximum pressure on first amplitude simulated by this model was overestimated. (ii) The pipe-wall constraint coefficient α has a slight influence on obtained results. The lower values of α: the higher max. pressures (on first and following pressures amplitudes) and the faster following amplitudes occurs (see Fig. 5).
The differences between numerical simulations of water hammer in viscoelastic pipes and experimental data can be explained by: (i) Simulations do not take into account the slope of the pipe and local resistance in the elbows.
(ii) An experimental creep function was used, which haved a selected J 0 coefficient, because there are problems in determining the creep function for small time t < 0.1 [s] [15].  (iii) During simulations a constant pressure on the reservoir were assumed. But the analysis of experimental pressure runs shows that this pressure after the valve closure linearly increases (unfortunately the authors of experiment did not explain what is the reason, despite showing this increas on Fig. 10 on page 65 of their most cited paper [15]).

Conclusions
The paper presents a numerical model used to calculate the transient flow (caused by rapid closing of the valve) in the polymeric pipes. From presented exemplary simulation tests and obtained from them results the following conclusions can be formulate: a) There is a great need to conduct new experimental research in a simple horizontal reservoirpipe-valve system, in which the pipes are made of a polymeric material. As it is still lack of the accurate results of analysis of water hammer in polymeric pressurized pipes. All known studies [11,[13][14][15] were mostly carried out in more or less inclined pipes. In them occurred an additional local resistance (in elbows or in curvature of the pipes caused by coiling them onto a cylinder) that caused problem with collected experimental data. There is a major problem with separation of their influence and their proper numerical modelling. This is particularly visible from presented in this paper exemplary results. b) Used by many authors procedure adjusting J coefficients (commonly referred to calibration [15,42,43]) that role is to minimize the simulation errors and fit simulation results to experimental one is in the view of the authors of this study unacceptable. Indeed, such a calibration can be performed for existing systems, but the priority should be to find an analytical solution by means of which, it will be possible to predict the unsteady behavior of the liquid at any time in polymeric pipes. To this end, new material research are necessary needed to properly describe the mechanical behavior of the polymeric materials used for pressure pipes. c) There is a need for additional numerical simulations to determine the significant of the unsteady friction in these type of flows. d) The real observed velocity pressure wave speed is during transient flow variable and can be calculated from following equation [8,9,43,44]: That's the reason why one should consider the attempt to modify the classical numerical solution based on method of characteristics built on a solid rectangular grid.