Simulation of unsteady flow in viscoelastic pipes

The paper focuses on modeling of unsteady flows in a hydraulic system built of a pressure tank, a plastic pipeline and a quick-closing valve. The influence of unsteady friction as well as the experimentally obtained creep function (necessary for modeling retarded strains) on simulation results was investigated. In addition, the effectiveness of the dimensionless parameter P known from the literature was analyzed, especially in the context of rejecting an unsteady term decision. Detailed investigations on the variability of terms describing unsteady friction have shown that one should still look for a dimensionless parameter with the help of which it will be possible to decide on the friction model before making the calculations. The quantitative analysis carried out showed that the use of unsteady hydraulic resistance in simulated runs brought simulated results closer to experimental results.


Introduction
Plastics are starting to play an increasingly important role in hydraulic systems each year. As of today, most of the hydraulic elements can be made of them [24]. Plastic pipes almost completely replaced metal pipes (steel, brass, etc.) in water supply systems. This was mainly due to their low price. They are produced from at least five different polymers: PP-polypropylene, PE-polyethylene, PVC-polyvinyl chloride, PBpolybutylene, ABS-acrylonitrile butadiene styrene. Each of the above-mentioned materials is characterized by different mechanical properties (Poisson's ratio, Young's modulus of elasticity, creep function, etc.), whose values strongly depend on the temperature. In these pipelines, as well as in metal, transient states may occur during the flow of liquids. From previous work carried out by the authors of this study [28][29][30][31], it turned out that the phenomena accompanying unsteady flows in metal pipelines such as: cavitation [1,7,14,21,35], unsteady friction [5,13,[25][26][27] or the interaction of the walls of the pipe with the flow (fluid structure interaction) [9,11,16] also plays an important role in plastic pipes. The phenomenon that is most responsible for flow damping is the viscoelasticity of the material of the pipe [2,10,[28][29][30][31]. Attempts to simulate the flow using a model for elastic pipes result in an unacceptable discrepancy between simulation and experimental results. Thermal phenomena [17,18] are usually neglected during the modeling of the examined transient flow (water hammer) in this paper (e.g., result of sudden blockage of flow). As proved in [2,10], the application of the calibration enables satisfactory model compliance to be obtained using the quasi-steady model of hydraulic resistance. The need to calibrate the creep function is mainly due to the lack of accurate experimental studies on their course in time (frequency) and temperature. In this paper, a modified unsteady flow model will be used, which properly takes into account frequency dependent hydraulic resistance. Unsteady flows in two experimental systems in HDPE pipelines were analyzed. Viscoelastic damping was modeled using the experimental creep function known from the literature [4]. According to paper [8], one can specify a certain dimensionless parameter P, which is dependent on internal diameter, pressure wave speed, friction factor, velocity and length of the pipe. Some authors [6,15,33] use this parameter to determine the use of friction model. In this work, a transient flow model will be used, dedicated to plastic pipes, discussed in detail in [28,30]. Analysis of the results obtained will demonstrate the need of using unsteady hydraulic resistance.

Water hammer in viscoelastic pipes
This section is devoted to introduce the equations describing water hammer in viscoelastic pipes and presentation of numerical solutions.

Mathematical model
Polymer pipelines do not respond according to Hook law when subjected to a certain instantaneous stress. Strain can be decomposed into a sum of instantaneous elastic strain e and a retarded strain r (see e.g., [3]) The retarded strain is a convolution integral of pressure and derivative of the creep function J which describes viscoelastic behavior of the pipe material (1) (t) = e + r (t). (2) where J i -creep compliance of the spring of the Kelvin-Voigt i-th element defined by J k = 1∕E k , E k -modulus of elasticity of the spring of i-th element, T i -the retardation time of the dashpot of i-th element.
Equation (2) can be written in the form [30] and its partial derivative with respect to time t as where w J (u) denotes creep weighting function and In unsteady flow in pipe, the instantaneous wall shear stress may be regarded as the sum of two components [32]: where -Darcy-Weisbach friction factor [−], w(t) -weighting function [−], -dynamic viscosity kg m s , utime, used in convolution integral [s].
The first component in (7), q , presents the quasi-steady wall shear stress and the second, u , is the additional contribution due to unsteadiness. Equation (7) relates the wall shear stress to the instantaneous average velocity and to the weighted past velocity changes.
Unsteady flow of liquid in viscoelastic pipe is represented by two one-dimensional hyperbolic partial differential equations. The momentum and continuity equations have the following form [31] where t-time

Numerical solution
There is no known analytical solution for system of hyperbolic partial differential equations (8). Therefore, there is a need to perform calculations using numerical methods. In this paper, a method of characteristics and finite difference method with classic constant rectangular grids is used, to avoid interpolation problems (Fig. 2). At the beginning, the initial conditions of the system will be discussed. In a pre-transient state, the mean velocity v 0 is constant on whole length of pipe x = L . The reservoir pressure (at x = 0 ) is constant and decreases linearly in the direction of the flow. At t = 0 , the valve closes suddenly ( x = L ) causing transient states, and the velocity changes from v 0 to 0 in one time step. After that, the pressure oscillation of fluid occurs until the full suppression. The final conditions on the entire pipe are v = 0 and p = p A (constant reservoir pressure). Now one can proceed to the numerical solution. The system of differential equations (8) was solved using method of characteristics. Thanks to this, the system of partial differential equations was transformed into two ordinary differential equations Using the finite difference method on characteristics grid ( Fig. 2), it was obtained where Δt-constant time step [s]. The term connected with strain rate is calculated as follows [30] (9) ± 1 c dp dt

with constants
Let now denote all the terms calculated for the previous time step as Finally, one gets a numerical solution for pressure and flow velocity The formulas for pressure and velocity at the pipe ends will now be derived. For the reservoir ( x = 0 ), the negative characteristics C − are used (Fig. 2). Let us assume that the pressure is constant ( p = p A ) and therefore, one needs the formula for velocity only. Using the second equation from (10), one obtains for the left boundary For the valve ( x = L ), the positive characteristics C + are used. Let us assume that the velocity is constant ( v Z = 0 ) so one needs the formula for pressure. Using the first equation from (10), one obtains Numerical calculations of time dependent of the wall shear stress u [second component in (7)] can be performed using the efficient numerical solution of this integral presented by Urbanowicz [27] where The constants n i and m i present in the above equations are associated with the effective weighting function and is the correction factor [25].

Experimental setup
To verify carried out numerical simulations, they were compared with experimental data. The experiments were carried out by D. Covas et al. at Imperial College London [3] and Evangelista et al. at the University of Cassino and Southern Lazio [6].
The first experimental reservoir-pipeline-valve system [3] consists of a 271.7 m HDPE pipe with an internal diameter 50.6 mm and 6.3 mm wall thickness. Form the analysis of observed experimental dynamical courses of pressure changes, the pressure wave speed in the tested system was determined at c = 395 m/s. The water temperature is T = 20 • C , on its basis, the density = 998.2 kg/m 3 and kinematic viscosity = 1 × 10 −6 m 2 /s were determined. The Poisson's ratio of the material from which the pipeline is made was P = 0.46 . Other input parameters for numerical calculations are summarized in Table 1. Parameters T k and J k of the Kelvin-Voigt models experimentally determined by Covas [4] have values presented in Table 2 (see Case 1, 3). The pipe wall constraint coefficient is calculated form the following formula [3] and for this test stand have a fixed value = 1.0647.
The second test stand was built to investigate the transient states in branched systems of HDPE pipes [6]. The system was inspired by irrigation system and is Y-shaped with three-way junction. The authors carried out preliminary experiments in single-pipe systems, which were used in this paper. The experimental system consists of a 203.3 m HDPE pipe with internal diameter 44 mm and wall thickness 3 mm. The average value of propagation of the pressure wave speed was determined at 368 m/s. The water temperature is T = 14.95 • C , on its basis, the density = 999.1 kg/m 3 and kinematic viscosity = 1.14 × 10 −6 m 2 /s were determined. The authors of the original reference did not mention the value of Poisson's ratio, so it was assumed the same value as Covas for HDPE pipes, i.e., P = 0.46 . The pipe wall constraint coefficient = 0.9372 was calculated using equation (18). To describe the creep function, the results of the Covas experiment were used. Parameters T k have been set as above. However, the J k parameters have been scaled based on the results of the Güney experiment [30]. Güney determined the parameters of the creep function at different temperatures. The scaling factor was assumed to be 0.75 which gave values presented in Table 2 (see Case 2, 4).

Comparison of experimental results and numerical simulations
To evaluate the effect of unsteady friction on the modeled pressure waves, numerical simulations were carried out for four different cases. Two of them concerned laminar flow and the other two were turbulent. In addition, the dimensionless parameter P defined by Ghidaoui et. al [8] were examined. The values of this parameter for the analyzed below the selected flows are shown in Table 1.
As described in paper [8], the physical meaning of the parameter P is the ratio of the diffusion time scale to the wave time scale in water hammer problems. When P ≫ 1 ,  the quasi-steady friction model is acceptable. When P ≈ 1 , it has been indicated that a unsteady friction model should be used. When P ≪ 1 , it does not belong to the water hammer problems, so similarly the quasi-steady friction model should be acceptable.
In  Figure 3 provides detailed comparisons between measured and calculated results corresponding to the quasisteady and unsteady friction model, which included the viscoelastic behavior of the pipe wall material. On their basis, the following conclusions can be formulated 1. The pressure peak on the first amplitude is better simulated using unsteady friction model.  On the other amplitudes, low pressures (runway valleys) are more accurately simulated with use of unsteady friction model, while high (amplitude peaks) using a quasi-steady friction model in low Reynolds number cases (Case 1 and 2). 3. Both models have a slight shifting over time. In the case of a quasi-steady friction model, the time between successive amplitudes is shortened, while in the unsteady friction model it becomes longer. 4. The qualitative analysis in laminar flow (Case 1 and 2) shows that the unsteady friction model better fits to the experimental data. 5. The qualitative analysis of turbulent flow (Case 3 and 4) shows the relative compatibility of both models. 6. The other differences observed in the comparative runs are the result of the neglecting of less significant phenomena accompanying the water hammer, as well as the result of experimental research in systems with many elbows, which are the result of increased hydraulic resistance. [23]. Figure 4 presents a comparison of dimensionless velocity ( v(t) = v(t)∕ max |v(t)| ) with a dimensionless quasi-steady wall shear stress ( ̂q(t) = q (t)∕ max | q (t)| ) and a comparison of dimensionless acceleration ( ̂v (t) =v(t)∕ max |v(t)| ) with dimensionless unsteady wall shear stress ( ̂u(t) = u (t)∕ max | u (t)| ) component in the middle section of the pipeline. Normalization was done in order to better present the studied phenomenon. One can observe the phase compatibility of the velocity with a quasi-steady wall shear stress component as well as phase

Fig. 4
Comparison of wall shear stresses for the middle of the pipe in a dimensionless form in Cases 2 (laminar) and 3 (turbulent) for unsteady friction model compatibility of the acceleration with a unsteady wall shear stress component. Thus, the total wall shear stress is a function dependent on velocity and acceleration. The above observation is consistent with the form of wall shear stress = f (v,v) being the sum of two components [12,19,20,22,34]: one proportional to the average speed and the second proportional to the acceleration. In addition, one can observe coverage of dimensionless velocity and quasi-steady wall shear stress plots in laminar case. In turbulent flow, this phenomenon does not occur due to the changing nature of the flow and the use of various models of quasi-steady wall shear stress in laminar and turbulent flows. On the basis of Fig. 5, the dominant role of the unsteady wall shear stress component on total wall shear stress can be observed for low Reynolds numbers (Case 2). On the other hand, for the large values of the Reynolds number (Case 3), the role of the quasi-steady and unsteady wall shear stress component is comparable. Similar observations also take place in Cases 1 and 4. The greater the Reynolds number will be, the smaller will be the impact of the unsteady friction term u . One of the assumptions of this paper was to examine whether using the P parameter one can make decisions about the choice of the use of a hydraulic resistance model. Therefore, the dynamics of changes of this parameter in different pipe cross-sections was also investigated. Selected plots of the parameter P for the case of the unsteady friction model are shown in Fig. 6. The above plots questioned the reasonableness of using the dynamic parameter P described by the Eq. (19), because it can be seen that whenever the flow velocity dropped and a change in its value occurred (as a result of a change in the flow direction), its singularities occurred. An increase in the P > 30 value in accordance with [6,15,33] should be considered a sufficient reason to reject and not designate unsteady term of friction. However, the attempt to reject the unsteady wall shear stress u from Eq. (7) in a dynamic manner, carried out as part of the simulation work, was unsuccessful. Generally, it is also known that the effect of unsteady friction increases with the decreasing Reynolds number, so that the behavior of the parameter P would be contradictory to the above.

Quantitative comparative analysis
In the previous section, the results of numerical simulation was presented. Calculated values of the initial parameter P (Table 1) suggests that the better results should be obtained using the quasi-steady friction model in Cases 1 and 2, while one should obtain better results using the unsteady friction model in Cases 3 and 4, but the qualitative analysis conducted does not give a definite answer which of the friction models gives better results. To this end, L 2 norm of errors was introduced (see e.g., [2]) and denoted as ‖ ⋅ ‖ 2 . Considering numerical p sim and measured experimental p exp values of pressure as a vector, the L 2 norm is defined here as where N p denotes the dimension of the vector p exp .
In the all considered cases, a smaller L 2 norm of errors was obtained in numerical simulations using unsteady friction model. At the same time, the computation time is higher than in the case of a quasi-steady friction model. The details are presented in Table 3.

Conclusions
The paper presents further analysis of the influence of a quasi-steady and unsteady friction model in viscoelastic pipes on simulated pressure runs. The model of water hammer in viscoelastic pipes was analyzed. Additional term describing the retarded deformation of the pipe wall was added to continuity equation. System of partial differential equations describing this type of flow was numerically solved using the method of characteristics and finite difference method. The research considered a qualitative and quantitative analysis of new comparisons of experimental results [6] with simulated ones. Unfortunately, the qualitative analysis does not give definite  answers as to the validity of the hydraulic loss model used on the basis of the P parameter.
In the absence of a clear answer, a quantitative analysis was used. The L 2 norm of errors was introduced to compare the simulated pressure waves and results of the experiment. A detailed analysis showed that the unsteady friction model minimizes L 2 norm of errors.
The research carried out shows that further work is needed to determine proper dimensionless coefficient used for initial and following estimations of the importance of selected friction model for numerical simulations.
In addition, the paper proves that experimental creep function is suitable for numerical simulations in various systems based on HDPE polymer pipes. The flow model analyzed in this work allows simulation of transient states occurring in plastic pipes. It is especially important in the design stage of the fluid systems: hydraulic, water supply, transmission, etc., because with it the systems are optimized in the selection of physical parameters for better operation and management.