Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90 © 2020 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2019.6324	Original Scientific Paper
Received for review: 2019-09-11 Received revised form: 2019-11-07 Accepted for publication: 2019-12-03
Transient Liquid Flow in Plastic Pipes
Kamil Urbanowicz1* - Huan-Feng Duan2 - Anton Bergant34 1 West Pomeranian University of Technology, Department of Mechanical Engineering and Mechatronics, Poland 2 The Hong Kong Polytechnic University, Department of Civil and Environment Engineering, China
3 Litostroj Power d.o.o., Slovenia 4 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia
Nowadays, plastic pipes play a key role in fluid conveyance, for example, in urban water supply systems. Dynamic analysis of designed water pipe networks requires the use of numerical methods that allow for solving basic equations that describe transient flows occurring in plastic pipes. In this paper, the main equations were formulated with pressure (p) and velocity (v) as the principal variables. A novel simplified retarded strain solution is used to properly model the pipe-wall viscoelasticity effect during transient flow process. Unsteady friction is calculated with the use of a filtered weighting function (with three exponential terms). The proposed numerical method enables the efficient modelling of transient flow in plastic pressurized pipes. Extensive simulations are performed and compared with experimental results known from three different European research centres (London, Cassino, and Lyon), with the aim of demonstrating the impacts of plastic pipe properties and frequency-dependent hydraulic resistance on transient pipe flow oscillations. The effectiveness of the proposed method for determining the creep functions of plastic pipes is also examined and discussed.
Keywords: plastic pipes, viscoelasticity, water hammer, transient flow, method of characteristics, unsteady friction Highlights
•	A proper mathematical model is derived based on the principal variables of pressure and velocity.
•	An efficient numerical solution of transient flow model in plastic pipes is used.
•	A simplified, effective numerical convolution integral describing the retarded strain is presented.
•	A detailed qualitative and quantitative analysis for the results of comparative tests is carried out.
•	A discussion on the influence of unsteady friction on transient flows is presented.
0 INTRODUCTION
The work of Tison, written in 1958, [1] has been considered to be foundational regarding the dynamic behaviour of the viscoelastic pipe wall and its influence on the transient flow behaviour associated with a rapid valve opening. In 1962, Hardung [2] described the physical phenomena that come into play in the arterial system as a consequence of heart action. He studied the influence of internal wall friction and liquid viscosity on damping of pressure wave velocity. From the linearized dynamic equations of viscous incompressible fluid flow in viscoelastic tubes, Martin [3] calculated frequency response of specific tube flow system in 1969. Lorenz and Zeller [4] analysed wave propagation in viscoelastic pipes, to find the analytical solution based on an analogy with a piston problem. Their main conclusion was that the shock waves could not be formed in viscoelastic pipes like those in elastic ones. Instead, a fully dispersed wave is formed, which has a width of several times of the tube diameter. In these early scientific studies, the main focuses was in analytical analysis.
A milestone to modern modelling was attributed to deriving general equations of the transient flow in pipes with viscoelastic properties by Rieutord and
Blanchard in 1972 [5]. In 1977, Güney, in his PhD thesis [6], presented the first working numerical method (based on finite difference method) that included a laminar unsteady friction [7] and was able to calculate the pressures and mean velocities of flow in plastic pipes during water hammer event including column separation. For cavitation modelling, Güney successfully applied the simplest discrete vapour cavity model. Meanwhile, Rieutord and Blanchard [8] presented a theoretical study in a one-dimensional form on the effect of viscoelastic material properties on the pipe-conveying process. Their analysis results indicated that, in addition to exponential attenuation of wavefronts, the disturbance is transformed into a diffusion front. Limmer and Meißner [9] incorporated the complex creep compliance into the transient flow models and then derived the wave speed and damping factor for an oscillating pressure wave propagating in a thin-walled viscoelastic pipeline. Franke and Seyler [10] utilized an impulse response method to calculate water hammer.
In Italy, in the mid-1980s, the idea of reducing the unsteady flow oscillations in pressure pipelines by inserting an in-line polymer section, which is characterized with low pressure wave speed, was examined by Ghilardi and Paoletti [11]. This
*Corr. Author's Address: Department of Mechanical Engineering and Mechatronics, West Pomeranian University of Technology Szczecin, Piastow 19, Szczecin 70-310, Poland, kamil.urbanowicz@zut.edu.pl
77
Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
solution was then further studied theoretically and experimentally by Pezzinga and Scandura [12]. The authors used a Kelvin-Voigt model with one element. The coefficients of creep function were calibrated based on experimental results. The additional plastic pipe in their tests significantly reduced unsteady-flow oscillations. Tijsseling [13] presented an extended literature review about FSI (fluid-structure interaction) in steel and plastic pipes.
At the beginning of this millennium, Covas et al. investigated transient flows in a laboratory high-density polyethylene (HDPE) pipe with a length of 277 m [14] to [16]. A new simplified model for simulating viscoelastic transient pipe flows was implemented into the numerical scheme of the method of characteristics (MOC), which allows a fast calculation of pressure runs in viscoelastic conduits. Thereafter, Duan et al. [17] argued that an inaccurate representation of unsteady friction or turbulence in one-dimensional (1D) model may lead to non-physical calibration of viscoelastic parameters and that the influence mechanism of pipe wall viscoelasticity is totally different from friction effects (steady and unsteady) in terms of both amplitude damping and phase shifting of transient pressures. Their results have been validated by the experimental test data of Brunone and Berni [18]. Similar results regarding the contribution of unsteady friction and viscoelasticity were also obtained by Bergant et al. [19]. The fluid-structure interaction (FSI) in plastic pipes was studied by Keramat et al. [20]. The authors found that damping in transient flow may be attributed to friction effects (steady and unsteady), valve resistance, gas (air) in liquid (dissolved and trapped), pipe-wall viscoelasticity, as well as fluid-structure interaction and other non-elastic behaviours of supports. As many of these effects are unknown and have not been properly included in the transient solver, the KelvinVoigt model not only represents viscoelasticity but also all the other effects. As a result, the calibrated results of viscoelastic parameters for such complex situations (i.e., various influence factors) may become non-physical.
Pezzinga et al. [21] used a microgenetic algorithm to calibrate creep compliance function coefficient. Calibration of Kelvin-Voigt models with 2, 3, 5, and 7 parameters, respectively, proved the substantial independence of the elastic modulus and the dependence of the retardation time on the pipe period (i.e., the pipe length). In other recent studies, e.g., [22] to [24], it is concluded that the one-element Kelvin-Voigt model is reasonably accurate to achieve satisfactory numerical simulation for typical transient
flows in common plastic pipes. Kodura [25] showed that the characteristic of butterfly valve closure has a significant influence on water hammer in polyethylene (PE) pipes for the valve closure time larger than 25 % of the wave characteristic period. For shorter closure time, good agreement could be obtained between the model and experimental test in terms of maximum pressures. More recently, Ferrante and Capponi [26] introduced a generalized Maxwell model, based on fractional derivatives, to better represent the viscoelastic behaviour of plastic pipes during transient flow process. Through different tests for HDPE and oriented polyvinyl chloride (PVC-O), they concluded that the proposed model performs slightly better than the well-known generalized Kelvin-Voigt model. Today, not only pipes but also many hydraulic devices can be made of different plastics, as needed [27]. As a result, various previous studies on unsteady flows in plastic pipes have demonstrated the importance and necessity of this research subject. Note that in plastic pipes during water hammer process, the same phenomena may occur that are known for the flow in metal pipes, such as cavitation, frequency-dependent friction and fluid-structure interaction (FSI). However, in plastic pipes, the main factor responsible for transient damping is the time-dependent viscoelastic material behaviour, which is the main difference from the elastic pipes studied previously.
The objective of this work is to present a mathematical model and its effective numerical solution method that allows the simulation of the phenomenon of transient liquid flows in plastic pipes with sufficient accuracy and efficiency. The model will consider two important factors: unsteady friction of pipe flows and viscoelasticity of retarded strain. These two factors described by convolutional integrals require effective solutions to accelerate the calculation process. The integral describing unsteady friction was effectively solved by Urbanowicz [28], using the suggestions of Vardy and Brown [29]; in contrast, the integral for the viscoelastic nature of retarded strain has been previously solved in a complicated manner [14] and [30]. In this study, it will be solved in a much simpler way, using Schohl's effective solution [31], to further improve its effectiveness. An attempt is made to assess, in both qualitative and quantitative manners, the impacts of pipe-wall viscoelasticity and unsteady friction on transient flows in plastic pipes. The situations of laminar, transitional, and turbulent flows in pipes are tested and analysed.
78
Urbanowicz, K. - Duan, H.-F. - Bergant, A.
Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
1 MATHEMATICAL MODEL
The hyperbolic partial differential equations of motion and continuity, which are described by the basic parameters representing flow: p pressure and v averaged velocity in the cross-sectional area of the pipeline, including the pipe wall friction in a horizontal pipe, can be defined in the following form
[5]:
1 dp i dv : „j dP,
- +--+ ^
pc dt dx
(u)
w.,, du = 0 dt J(t-u)
p
dv dp 32p 16p r dv(u)
. (1)
- + — + dt dx Dl
d2 0 dt
w^-u )du = 0
Among the different numerical methods enabling resolving the system of above equations, particular attention should be paid to the MOC scheme, which perfectly interprets the essence of physical phenomena of transient flow, and at the same time is characterized by fast convergence, ease of incorporating various boundary conditions, together with the high accuracy of calculation results. The detailed numerical solution may be found in recent papers [30], [32] and [33].
1.1 Simplified Numerical Solution of Retarded Strain
In the case of flow in plastic pipes, the equation of continuity has the same form as for metal pipes. A polymer pipeline does not respond according to Hook's law when it is subjected to a certain instantaneous stress a. It consists of an immediate-elastic response and a retarded-viscous response. Thus the strain can be decomposed into a sum of instantaneous-elastic strain ee and a retarded strain er [6] and [16].
The partial time derivative of retarded strain is a convolution integral of pressure and derivative of the creep function J that describes viscoelastic behaviour of the pipe-wall material as:
deT(, ) S \ dp(u )
—— = — J —— • w
dt
2 o dt
J (t-u)
du.
(2)
An efficient numerical solution scheme to the above type of convolution integral has been developed by Schohl [31]:
de
(t+At)
dt
z
-e
At Rt,
S
2 At
1 - e
Y(l+Al)
(P(t+At ) P(t ))
(3)
where y, is the time-dependent recursion coefficient, of which the starting value (before transient) is equal to zero. Re-ordering the equation yields:
ds	kk
r{l+Al)-=p^YX-lL((^r¡(Y, (4)
dt
where: Y = — e ' 2
Rt' and X. =-J ' 2 At
1 - e
Let us now define the time-dependent term G(i) in Eq. (4):
i=l
(5)
Then, the final form of the simplified equation that describes the derivative of retarded strain is as follows,
ds
(t+At )
dt
'' P(t+At )
X - G,
(t )•
(6)
Comparing with an existing solution that was introduced by Covas et al. [14] and [16], this novel solution simplifies the numerical algorithm for analysis of transient events in engineering polymer pipes. It may be very useful for the various pipe engineering systems: design of new pressure systems based on plastic pipes; the modernization of existing plastic pressurized pipe systems; and design strategies to protect systems against large transient loads.
1.2 Numerical Solution of Wall Shear Stress
The pipe wall friction effect plays a significant role in modelling transient events in elastic pipes. In viscoelastic pipes, some authors note the importance and impact of unsteady friction on transient flow process [6], [16], [17] and [33], while the vast majority try to simplify the model to a minimum and model unsteady friction by calibrating coefficients in the creep function. The action of the latter means that the calibrated function is no longer a material function describing the strength properties of the pipe material. This function, in a simplified way, describes together all the phenomena that accompany the transient flow (including but not limited to: unsteady friction resistance, the impact of FSI and the phenomenon of energy dissipation resulting from retarded strain of pipe walls) [20]. In this work, due to the fact that only experimentally obtained creep functions are used, there is a need to model resistance in a conventional way, i.e., by including the unsteady component tu contained in the formula for describing the wall shear stress. Urbanowicz [28] has shown that there is a lack
0
Transient Liquid Flow in Plastic Pipes
79
Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
of compliance of the effective friction weighting function with original classic weighting function for small time scales. The corrected final effective solution of convolution integral representing the additional contribution due to unsteadiness becomes [28]:
Pfv{
(t+At) I (t+At ) I
(t+At )
d h
8
Ay
[V(t+At )- V(t)}
M]C [v(( )--At )]
(7)
where
n =
f4
J 0
(u ~)du
f A0'WeffXU )dU
y/(t+At)
A = e-
B = [1 - A-],
n At
C = AB-
(8)
The correction factor n in the above solution is a quotient of the integral of the classic weighting function wc[ass. (for laminar flow, Zielke's function [7], and for turbulent flow, Vardy-Brown's function [34]) and the integral of effective weighting function weJf (approximated by a finite sum of exponential terms):
"eff.
=2>
(9)
which is just used to approximate the conventional solution. It is a reminder that integration range of both integrals are the same from 0 to first dimensionless time A. The mt and nt coefficients are functions of the dimensionless time step At. Their values can be determined with use of the analytical formulas presented in [35]. The range of applicability of the filtered weighting function is from Af up to 103 Af.
observed experimentally (EXv) and simulated by using unsteady (Sv UF) and quasi-steady hydraulic resistance ($>,qsf), and the time are determined by the following formulas:
P = P(^f = ZÍO-^ and Í = t / (4L / c). (10)
Ap	pcv0
For quantitative analysis, it is essential to define a minimum of two coefficients, whose role is to mathematically determine the compliance of the simulated pressure histories with respect to the experimental tests. To the authors' knowledge, only few studies in the literature [24], [28], [36] and [37] in the field of water hammer flows are concerned on such quantitative analysis. Urbanowicz's analysis carried out in [28] was based on pressure histories in dimensional form. However, that analysis method may not be applicable to achieve the purpose of this study. Specifically, for laminar flows, the pressure pulsation from the value of the final pressure to which the transient state ends (equal to the pressure inside the reservoir) is small. Also, the pressure drop along the length of the pipe becomes lower as the Reynolds number is smaller. This is particularly noticeable in laminar flows with very small Reynolds numbers. In preliminary tests, the values of quantitative parameters determining the pressure compliance were calculated using the method presented in [28], which however resulted in very small values for laminar runs. Hence, in this work, to avoid the influence of the Reynolds number on this analysis, the quantitative analysis has been carried out in a dimensionless form. For this purpose, a subprogram was defined to search for maximum and minimum values and their occurrence times (calculated from the beginning of the analysed transient state). For demonstration, Fig. 1 illustrates the procedure of "collecting" maximal and minimal pressures, and their relevant times (circles).
2 COMPARISON OF NUMERICAL SOLUTION WITH EXPERIMENTAL RESULTS
To assess the effects of pipe-wall viscoelasticity and unsteady friction on the transient flows, extensive comparative studies have been carried out by researchers. Experimental results obtained by the researchers from European scientific centres (London (UK), Cassino (Italy) and Lyon (France)) are first used to validate the numerical model described in the previous chapter of the paper. For convenience, the comparisons of all the results are presented in dimensionless form. Then the pressure changes
Fig. 1. Selecting the maximal and minimal pressures
and their times
80	Urbanowicz, K. - Duan, H.-F. - Bergant, A.
Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
On the one hand, the coefficient determining the compliance of the maximum and minimum simulated pressures is calculated as follows:
s :
e=-
Ps,i - P e Pel
:
•100%,
(11)
where psi is simulated maximal or minimal dimensionless pressure, and pei is experimentally recorded maximal or minimal dimensionless pressure. On the other hand, the coefficient that determines the time compliance of subsequent simulated amplitudes was calculated as follows:
Z k
E, =-
ts,i - te,i te,i
k - 1
•100%,
(12)
where ts,i is the simulated dimensionless time ^of maximal or minimal dimensionless pressure and ^ is the experimentally recorded dimensionless time of maximal or minimal dimensionless pressure.
The degree of simulation compatibility and accuracy increases with decreasing values of the above coefficients (more accurate, matching of the simulated pressure runs with respect to the known experimental runs). Tables 1, 2 and 3 (presented later in this paper) summarize the quantitative analysis results of the Ep and Et parameters for all test cases in this work.
The wave speeds were estimated from the empirically observed durations of the first pressure amplitude. Then the instantaneous creep coefficient J0 can be estimated by the formula:
Jo =-
1
pac
1
K, a
(13)
D
where S = —Ç is expanded pipe constraint
e
coefficient. For Covas and Evangelista test stands, the pipe constraint coefficient is [3] to [5]:
i-1 (1 + V
D
D + e
(1 -vi )
and for Güney test stand [6] and [7]:
e ] 2e 2 i - e
1 + l DJ I1 -D
(14)
(15)
In all experiments, the pressure inside the supply reservoir increased with different extents during the water hammer event. The largest increase was observed in the experimental runs carried out by the Evangelista et al. [38] in Casinno (Italy). A relatively smaller increase was recorded during tests carried out
by Covas [16] and [39] in London (United Kingdom). Finally, the smallest increase in pressure, inside the supply reservoir, was noted in tests carried out by Guney [6] and [40] in Lyon (France). The main reason for the increase in pressure might be attributed to the poor synchronization of the rapid closure of the valve with the shut-off of the working pump. In the pre-transient state the reservoir pressure pri=0 is the sum of the theoretical pressure loss over the pipe length ApL and the pressure value recorded in the valve cross-section pvt= 0:
Pr ,t=0 = Pv,t=0 +
fLp v2 D 2
= Pvt=0 + APl
(16)
where Darcy-Weisbach friction factor for laminar flow is f= 64 /Re, and for turbulent flow (assuming
smooth pipe walls) f =	. With use of the
above-calculated reservoir pressure at the start of transients prt=0 , and known final pressure inside the reservoir pFinai , the value of the pressure increase in the cross-section at the reservoir Apr can be estimated approximately by:
so that
APr = PFinal - Pr ,t=0
PFinal =APr +APl + Pv,t=
(17)
(18)
For all Covas cases with known overpressures, there is no information about valve pressure pVtt=0 in steady flow prior to transient. Based on the residual information available from their works [16] and [39], the values of pvt= 0 for Covas test stand were assumed to achieve the simulations (see Table 1). For all test systems, the information needed to conduct simulations (pipe length, inner diameter, wall thickness, the Poisson's ratio, etc.) is collected in Tables 1, 2 and 3.
The creep function used for Covas test stand simulations (water temperature T = 20 °C), which was needed to calculate the retarded strain, was archived from experimental test [16], by giving the following creep and retardation time coefficients: J0 = 0.674-10-? Pa-i, J = 0.1394-10-9 Pa-i, J2 = 0.0062 10-9 Pa-i, J3 = 0.114810-9 Pa-i, J4 = 0.3425 10-9 Pa-1, J5 = 0.0928 10-9 Pa-1 and Rt0 = 0 s, Rt1 = 0.05 s, Rt2 = 0.5 s, Rt3 = 1.5 s, Rt4 = 5 s, Rt5 = 10 s. The parameters of creep compliance functions vary strongly with temperature. Especially for Evangelista et al. test stand (water temperature T = 15 °C), the above function has to be scaled. The scaling factor was assumed, as in work [41] to
2
0'
2
Transient Liquid Flow in Plastic Pipes
81
Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
be 0.75 which gives the following values for the simulations: J0 = 0.513-10-9 Pa-1, J = 0.1046-10-9 Pa-i, J2 = 0.0046 10-9 Pa-1, J3 = 0.0861 10-9 Pa-1, J4 = 0.2569 10-9 Pa-1, J5 = 0.0696 10-9 Pa-1 and Rt1 = 0.05 s, Rt2 = 0.5 s, Rt3 = 1.5 s, Rt4 = 5 s, Rt5 = 10
s. The creep function for Guney test stand was experimentally defined for different temperatures. Original values received by Guney are required to be filtered according to the filtration procedure proposed by Urbanowicz and Firkowski [33]. The obtained
Table 1. Covas test stand details and quantitative results
Case	vo	Reo [-]	pv.t= o [Pa]	Apl [Pa]	Apr [Pa]	Unsteady friction (UF) [%]		Quasi-steady friction (QSF) [%]		AEp	AEt [%]
	o [m/s]					AEp.UF	AEt.UF	AEp.QSF	AEt.QSF	[%]	
CO1	O.O2S	1417	4.7S5e5	O.O1e4	O.14e4	5.45	O.45	S2.6S	1.7O	27.2S	1.25
CO2	O.124	6274	4.4O5e5	O.15e4	O.45e4	2.OS	O.97	21.4S	1.66	19.4O	O.69
COS	O.249	12599	4.4O5e5	O.5Oe4	1.OOe4	2.77	1.2O	15.S6	1.S7	1S.O9	O.67
CO4	O.S7S	1SS74	4.4O5e5	1.O1e4	1.S9e4	2.O1	O.45	14.94	1.5S	12.9S	1.OS
CO5	O.5OO	25SOO	4.2S9e5	1.65e4	2.OOe4	S.9S	O.29	14.4O	1.2O	1O.47	O.91
CO6	O.746	S774S	S.425e5	S.S9e4	S.11e4	S.O4	O.5S	12.O5	1.SO	9.O1	O.72
CO7	O.S7O	44O22	S.425e5	4.4Se4	S.SOe4	S.41	O.46	11.21	1.SS	7.SO	O.S7
COS	O.995	5OS47	4.7S5e5	O.O1e4	O.14e4	S.7O	O.56	11.5S	1.17	7.SS	O.61
L =	271.7 m; D	= O.O5O6 m; e = O.OO6S m; vP			= O.46 ;	= 1.O65; S =	S.551; T= 2O	"C; p = 99S.1 kg/mS; v = 1O-6 m2/s; c =			4OO m/s
Table 2. Evangelista test stand details and quantitative results
Case	vo [m/s]	Reo [-]	pv.t= o [Pa]	Apl [Pa]	Apr [Pa]	Unsteady friction (UF) [%]		Quasi-steady friction (QSF) [%]		AEp	AEt [%]
						AEp.UF	AEt.UF	AEp.QSF	AEt.QSF	p [%]	
EO1	O.O7O	27O2	4.19Se5	496.5	S95O	7.SO	S.2S	S4.9S	4.S4	27.1S	1.O6
EO2	O.O9S	S7SS	4.1SSe5	S94.5	516O	4.47	S.15	SS.O9	4.S6	2S.62	1.21
EOS	O.164	6SSO	4.12Se5	22O2	S5OO	S.67	2.76	27.6O	4.19	2S.9S	1.4S
EO4	O.SS5	129SO	S.992e5	76S6	166OO	2.SO	2.7S	2O.6S	4.1O	17.SS	1.S2
EO5	O.657	25S5S	S.66Se5	2.5e4	16SOO	5.15	2.S6	14.94	4.S7	9.79	1.51
EO6	O.9S9	SS172	S.2S5e5	5.11e4	4S4OO	5.16	2.75	12.42	4.19	7.26	1.44
EO7	1.S15	5O754	2.7OOe5	S.414e4	67S6O	5.16	2.22	1O.6O	S.76	5.44	1.54
EOS	1.645	6S491	2.O45e5	1.245e5	9SOOO	5.7S	2.S9	S.61	S.SS	2.SS	1.49
L = 203.3 m; D= 0.044 m; e =0.003 m; vP = O.46 ; £= 0.937; S = 13.745; T = 15 "C; p = 999.1 kg/mS; v = 1.14-1O-6 m2/s; c = 365 m/s
Table 3. Güney test stand details and quantitative results
Case	vo [m/s]	Reo [-]	pv.t= o [Pa]	Apl [Pa]	Apr [Pa]	Unsteady friction (UF) [%]		Quasi-steady friction (QSF) [%]		AEp	AEt [%]
						AEp.UF	AEt.UF	AEp. QSF	AEt. QSF	p [%]	
GO1	O.49O	17422	1O1S25	S42S	475O	S.54	S.7S	14.2S	5.94	1O.74	2.21
GO2	O.55O	2565O	1O1S25	S9O6	5OOO	S.9S	2.77	19.57	4.15	1O.59	1.SS
GOS	O.57O	SO245	1O1S25	4O19	1OOOO	S.S1	1.S6	19.7S	1.S7	1O.92	O.51
GO4	O.55O	S1646	1O1S25	S695	1OOOO	S.4O	2.S5	11.91	S.54	S.51	1.19
GO5	O.56O	S451S	1O1S25	S74S	SOOO	9.S9	2.56	19.S9	S.77	9.5O	1.21
GO6	O.S2O	5O5S6	1O1S25	7296	1OOOO	S.S4	2.S2	16.S2	2.94	7.9S	O.62
L = 43.1 m; D= 0.0416 m; e =0.0042 m; vP = 0.38 ; £= 0.97; S = 9.61
Table 4. Values of temperature dependent Güney's research coefficients
Case	T ["C]	c [m/s]	p [kg/mS]	v [m2/s]	K, [Pa]	At [s]	Jo [Pa-1]	Ji [Pa-1]	J2 [Pa-1]	Rti [s]	Rt2 [s]
GO1	1S.S	S1O	999.S	1.171O-6	2.141O9	O.OO4S	1.OS5 1O-9	O.6S71O-9	O.S711O-9	O.O166	1.747
GO2	25	265	997.1	O.S92 1O-6	2.241O9	O.OO51	1.44O1O-9	1.O461O-9	1.2S71O-9	O.O222	1.S64
GOS	S1	247	995.S	O.7S41O-6	2.271O9	O.OO55	1.66S1O-9	1.S971O-9	1.62S1O-9	O.O221	1.S22
GO4	S5	24O	994.1	O.72S1O-6	2.2S51O9	O.OO56	1.772 1O-9	1.797 1O-9	2.S491O-9	O.O265	2.S92
GO5 GO6	SS.5	22O	992.6	O.675 1O-6	2.295 1O9	O.OO61	2.1211O-9	2.O971O-9	S.57O1O-9	O.OS47	S.O77
82
Urbanowicz, K. - Duan, H.-F. - Bergant, A.
Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
values and other temperature-dependent parameters details are collected in Table 4.
In all MOC-based simulations, a constant number N = 32 of reaches was used. Such a pipe division gave the following time steps: for Covas test stand At = 0.0212 s, for Evangelista test stand At = 0.0174 s and for Guney test stand, detailed values can be found in Table 4. A rapid-closing valve (tc[ < 2L/c) at the downstream end of the pipes was applied to all scenarios (Covas, Evangelista and Guney test stands) to generate water hammer flows in the system. In all of the three test stands, two distinct flow cases were selected (marked with gray in Tables 1 to 3) and investigated in detail.
2.1 Research Based on Covas Experimental Results
The results of the selected simulation cases C01 (laminar flow with initial: velocity v0 = 0.028 m/s and Reynolds number Re0 = 1410) and C05 (turbulent flow with initial: velocity v0 = 0.5 m/s and Reynolds number
Re0 = 25300) in relation to the Covas experimental data are plotted in Figs. 2 and 3, respectively.
It can be seen from Fig. 2a that the simulation runs realized for laminar flow (C01), taking into account in a simplified "step way" that the recorded pressure increases inside the supply reservoir, the model results present a high accuracy. Fig. 2b presents the comparison of the same experimental data with simulation results, but with implementing only the quasi-steady hydraulic resistance t = tq in the simulation. It can be observed that in this case the modelled damping (amplitude attenuation) was much smaller, which results in overstated simulation runs (Sv,qSF). Moreover, comparison in Fig. 2 reveals that a change in the method of modelling the wall shear stress could greatly affect the shift of simulated pressure runs - the use of unsteady friction delays the moment of occurrence of subsequent amplitudes. This means that unsteady friction indirectly affects the change (during the transients) of the pressure wave speed.
a)
b)
Fig. 2. Comparison of simulated and experimental results for Covas C01 (laminar flow); a) unsteady friction model, and b) quasi-steady friction model
a)
b)
Fig. 3. Comparison of simulated and experimental results for Covas C05 (turbulent flow); a) unsteady friction model, and b) quasi-steady friction model
Transient Liquid Flow in Plastic Pipes
83
Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
From the analysis of Fig. 3 obtained for the C05 for turbulent flow scenario, most of the results are similar to those just discussed for laminar case C01. Of the interesting noticeable differences between Figs. 2a and 3a showing the normalized Covas cases analysed in detail in this work (laminar C01 and turbulent C05), the shape of the first amplitude deserves additional attention. In the laminar flow case (C01), it can be found that the maximum value occurred in the initial period, i.e., just after the valve was closed, there was a significant pressure drop at the peak of the first amplitude (see Fig. 4). A similar tendency was also noted from the observation of experimental results. For the turbulent flow case in Fig. 3, it can be seen that there was only a slight decrease in pressure to the minimum value at the top of first amplitude, which maintained a constant value for a long time (normalized times from 0.17 up to 0.35), and, in the final phase, there was a slight increase in pressure. In Fig. 3, the pressure runs at the first amplitude top is almost flat. For clarity, an enlarged plot for this first amplitude is shown in Fig. 4 for a detailed comparison. It is suspected that the above change in pressure behaviour at the top of the first amplitudes is the result of a change in the relationship between two simulated phenomena: friction and retarded strain (both described by convolutional integrals).
Fig. 4. Auxiliary drawing - the top of the first pressure amplitude
From the quantitative results collected in Table 1 regarding experimental tests carried out by Covas, it can be concluded that:
• pressure runs simulated with applying unsteady friction presented a much better accuracy than those in which the resistances were modelled with a quasi-steady component only. This is because for all test cases (C01 up to C08), an inclusion of unsteady friction may reduce the values of Ep and Et coefficients;
•	the average values of Ep and Et parameters in the test runs, taking into account unsteady resistances, were Epu UF, meanC = 3.29 % and Et, UF,meanC = 0.62 %, respectively. However, for test runs simulated with quasi-steady resistances
Ep,QSF,meanC = 16.77 % and Et,QSF,meanC = 1.47
%, respectively;
•	the highest Ep value (means the worst simulation fit) was obtained for the laminar flow case. For the model with unsteady friction, the errors were: Ep,uf,co\ = 5.45 % and E,ufc01 = 0.45 %, while in the model with quasi-steady resistance: Ep,qsf,C01 = 32.68 % and EtqskC01 = 1.7 %. The possible reason of these poor results can be attributed to a large pressure signal noise (before filtering) noticed here for laminar results (the same problem also affected other tests with a low initial flow rate (Re < 6500) carried out by Evangelista et al.; see Table 2).
2.2 Research Based on the Experimental Results of Evangelista et al.
Two cases for the Evangelista test stand are analysed herein, in which the main difference is the initial velocity of the flowing liquid (Reynolds number). In the first case (E01) (transitional flow case), the initial flow velocity was Vq, e0\ = 0.070 m/s, and in the second analysed case (E08) (turbulent case) it was vo,Eo8 = 1 645 m/s. The initial velocity significantly affects not only the Darcy friction factor and unsteady friction calculation but also the pressure pulsations appearing after the rapid valve closures. Meanwhile, what follows from equations presented in the first chapter, the simulated retarded strains are also affected (which are usually treated by researchers as a source of additional damping of flow pulsations). Experimental results (especially for low flow rates) were characterized by high noise. For this purpose, the MatLab "filtfilt" filter was used, which performs zero-phase digital filtering by processing the input data in both the forward and reverse directions. After filtering the data in the forward direction, the "filtfilt" function reverses the filtered sequence and runs it back through the filter. The result has the following characteristics: zero phase distortion; a filter transfer function equal to the squared magnitude of the original filter transfer function. The results of selected simulation tests (E01 and E08) are presented in Figs. 5 and 6, respectively.
From Fig. 5a, it can be seen that the degree of matching of simulation results (using unsteady friction) to experimental results seems not perfect. This difference may be attributed to the following factors.
84
Urbanowicz, K. - Duan, H.-F. - Bergant, A.
Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
First, the Reynolds number Re = 2702 indicates that it is a transitional flow between laminar and turbulent. In these flows, according to experimental research related estimation of friction coefficients [42], there is a sharp change in the value of this coefficient. In the simulation, the friction coefficient and the unsteady hydraulic resistances were calculated in a
simplified manner, assuming the laminar nature of the flow f= 64/Re and mu n, representing the weighting function were used for laminar flow). Secondly, the experimental results were initially characterized by high noise, as well as visible effects of reflected waves (which may be the result of narrowing at pressure sensor connection points, sharp knees, intense
a)
b)
Fig. 5. Comparison of simulated and experimental results for Evangelista E01 (transitional flow); a) unsteady friction model, b) quasi-steady friction model
0.6 0.4 0.2
t 0 -0.2 -0.4 -0.6
n		-sv,uf '
		--EX v
		
i	1 \ l	
1
a)
2 3 il-1
b)
Fig. 6. Comparison of simulated and experimental results for Evangelista E08 (turbulent flow); a) unsteady friction model, b) quasi-steady friction model
a)
b)
Fig. 7. Enlargement of: a) top of the first amplitude, and b) first and second valley of the course of transient pressure changes
Transient Liquid Flow in Plastic Pipes	85
Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
vibration of the pipeline, etc.) on pressure courses, referring to the sharp increase at the middle of the first amplitude (Fig. 7a) as well as peculiar increase visible on the first and during subsequent valleys of the pressure course (Fig. 7b).
The quality of matching of the simulation and experimental runs may be considered satisfactory. From the results in Table 2, describing the quantitative results of Evangelista scenarios, it is evident that transient test runs modelled with additional unsteady friction could match better with experiment than those modelled using the quasi-steady friction only. Along with the increase in the Reynolds number (from case E02 represented with the Reynolds number Re = 3783), a significant decrease in the significance of unsteady friction on the maximum and minimums values of pressure histories was noticed in this system, and evidenced by the decrease of such difference of pressure compliance parameters:
AEp = Ep,QSF - Ep,UF .	(19)
For example, in the case E08 in which the highest flow rate was recorded, the difference AEp was only AEpE08 = 2 83 %. However, the unsteady friction should not be neglected even for large Reynolds numbers, because it affects the time compliance of simulated transient runs, by referring to the similar difference value of:
AEt = Et, QSF - Et, UF ,	(20)
for different Reynolds numbers as shown in Tables 1 and 2.
2.3 Research Based on Güney Experimental Results
For Güney test stand six cases of turbulent flow were analysed, in which the main difference was the water temperature. Only two of these tested cases are selected for detailed discussion in this paper. It is known that temperature may significantly affect the mechanical properties of the plastic pipe material. With such influence, a large change in the creep function may occur during the transient tests (Table 4). The creep function derivative is the main function affecting the modeled damping of pressure oscillation, and in this simulation it was filtered according to filtration procedure proposed by Urbanowicz and Firkowski [33]. The numerical simulation results of the experimental tests are provided in Figs. 8 and 9.
From the results of the test G01 for relatively low temperature T = 13.8 °C, it is noticed that the model including unsteady friction (Fig. 8a) may match with
high accuracy the top of the first amplitude (Fig. 10), and even the peak of pressure on the top of this first amplitude could be reflected in the simulation process. This confirms the acceptable accuracy for the used mathematical model. As concluded earlier in this work, the results with better matching accuracy were obtained after applying unsteady friction (Fig. 8a) than using the quasi-steady resistance only (Fig. 8b). From the results in Fig. 9, it is also observed that an increase in temperature in the same system may cause a decrease in the number of appearing amplitudes in the same time interval. For example, when the temperature was 13.8 °C within two seconds after the rapid closing of the valve, there were three pressure amplitudes; while for the temperature of 38.5 °C, only two amplitudes are noticed at the same time.
This is mainly due to the following two factors. The first is the change of the wave speed c, from cG01 = 310 m/s for temperature T = 13.8 °C to cG06 = 220 m/s for temperature T = 38.5 °C. The second factor may be attributed to the change of the creep function coefficients J and tt (see Table 4). The creep parameters describing the properties of the material may strongly depend on the prevailing flow temperature (as seen in Fig. 11 illustrating the effect of temperature on creep functions).
In the test case with relatively high temperature (Fig. 9), a large divergence of the modelled peak was observed just after the valve was closed (Fig. 10). It is difficult to clearly indicate the key reason for this. In this example, it is also difficult to confirm the positive effect of unsteady friction by analysing these results only (Fig. 9). Normalized pressure runs show only a slight advantage of the unsteady model over the quasi-steady one. Therefore, it is necessary to conduct a quantitative analysis on such results so as to determine the key factors affecting the modelling accuracy.
From the results in Table 3 that were obtained for quantitative analysis for the Guney experimental system, the average value of the obtained coefficients illustrating model compliance was much higher (Ep.uF.meanG = 7.24 % and Et,uF,meanG = 2.52 %) than for the other two quantified systems. This may be due to two aspects as follows. Firstly, these results were tested and collected a relatively long time ago (i.e., in 1977), which forced the use of pressure sensors with much larger ranges of measurement errors than those currently used. Secondly, the process of digitizing these results, carried out using the charts attached to the PhD thesis [6], could have also contributed to the erroneous reading of real experimental results. These charts were often characterized by casual information on the value of the ordinate axes. Even though these
86
Urbanowicz, K. - Duan, H.-F. - Bergant, A.
Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
a)
b)
Fig. 8. Comparison of simulated and experimental results for Guney G01 (T = 13.8 °C, Re = 17422); a) unsteady friction model, b) quasi-steady friction model
a)
b)
Fig. 9. Comparison of simulated and experimental results for Guney G06 (T = 38.5 °C, Re = 50536); a) unsteady friction model, b) quasi-steady friction model
0 0.1 0.2 0.3 0.4 0.5 0.6
t H
Fig. 10. Enlargement of top of the first amplitudes
Fig. 11. Courses of the creep functions used in the simulations
differences in this system, it can still be confirmed that the use of unsteady friction may result in a significant decrease of quantitative parameters, and thus the substantial increase of the simulation accuracy.
Based on all quantitative analyses of the results for these three considered test systems, it can be concluded that the simulations performed using the
experimentally defined creep functions by including additional unsteady friction effect may provide more accurate modelling results for transient flows in plastic pipes. However, it is also revealed that the lack of experimental data regarding the temperature effect on viscoelastic creep function is usually an issue that requires solutions in practical applications, because
87
Transient Liquid Flow in Plastic Pipes
Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
from this study the temperature is found to be an important factor affecting transient flow behaviour in plastic pipes. Furthermore, the proposed parameters of quantitative analysis in Eq. (11) and Eq. (12) in this study can be useful measures for characterizing and evaluating the importance of unsteady friction effect to transient flows with the existence of the pipe-wall viscoelasticity effect in plastic pipes.
3 DEFINING THE EFFECTIVENESS OF NEW PROPOSED METHOD
Further verification of the effectiveness of a computer program is performed herein for the newly developed simple mathematical formulas in this study. As an illustrative example, the laminar Covas C01 run was chosen (details of this case can be found in Table 1), in which the effect of discretization was examined. Computation times told obtained using the original calculation algorithm [30] (complicated solution of retarded strain convolution integral, and unfiltered long exponential weighting function consisting 26 terms is used), were compared with calculation times obtained with the help of new procedures discussed in this work tnew.
Table 5. Calculation details
N [-]	16	32	64	128	256	512
told [s]	0.30	1.82	11.65	157	1152	10272
tnew [s]	0.11	0.52	3.79	66	487	4377
Nodes	469 x17	940 x33	1882 x65	3766 x129	7535 x257	15073 x513
At [e-6]	66.3	33.2	16.6	8.29	4.15	2.07
The tested reservoir-pipe-valve (R-P-V) system was divided into a predefined number of reaches N in accordance with Table 5. An increase in discretization corresponds to the corresponding intensity of the calculation (growing number of computational nodes). The calculation time was measured using a MatLab functions "tic" and "toc". The "tic" command was introduced in the software after specifying the boundary and initial conditions, just before the loop function in which the transient values are calculated. The "toc" command was introduced in the computer program just after the loop function. Simulations were repeated using selected discretization for both calculation algorithms with at least five times. The table summarizes the minimal calculation times obtained. Calculations were carried out on a typical laptop for personal use. It can be seen from Table 5 that the new procedure reduces the calculation time by
more than two times (or even three times), in the case of the discretization commonly used in practice, i.e., N < 64.
4 CONCLUSIONS
The main results and findings of this paper can be summarized as follows:
(1)	The comparative studies have demonstrated the significant role of unsteady friction effects in plastic pipes during rapid transient events. These effects could usually be masked in many previous research studies during their creep function calibrations.
(2)	There is a strong need to conduct detailed experimental research on polymers used for pressure pipes, aiming to determine the exact temperature-dependent creep functions. Such tests cannot be carried out only on individual dynamical mechanical thermal analyser (DMTA) devices, because the frequency of forcing deformations of the tested material may turn out to be too low compared to those factors that affect pipes in real conditions, i.e., systems in which water hammer occurs. A detailed analysis of these creep function courses is required immediately after applying the step load.
(3)	A new experimental research is needed in straight and horizontal R-P-V (reservoir-pipe-valve) system, in which the local resistances can be minimized. This is because local resistances can significantly affect the flow, as shown by Stosiak et al. [43]. To this end, it is necessary to design pressure sensors working exactly inside the shut-off valve and to use the optimal rounding radius of the plastic pipe in the input cross-section located on the wall of the pressure reservoir.
5 ACKNOWLEDGEMENTS
The authors would like to thank Professor Angelo Leopardi from University of Cassino and Southern Lazio for providing the experimental results to one of the tested scenarios. The co-author Duan would like to thank the Hong Kong Research Grants Council (RGC) for the support on his research work (15201017 and 15200719). Finally, Bergant gratefully acknowledges the support of the Slovenian Research Agency (ARRS) conducted through the research project L2-1825 and the programme P2-0162.
88
Urbanowicz, K. - Duan, H.-F. - Bergant, A.
Strojniski vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
6 NOMENCLATURES
c pressure wave speed, [m/s] D inner diameter of the pipe, [m] e thickness of the pipe wall, [m] E0 Young modulus of pipe material, [Pa] f Darcy-Weisbach friction factor, [-] JçD compliance creep function, [Pa1] J0 instantaneous creep coefficient, [Pa1] J creep-compliance of the ith spring of the KelvinVoigt model, [Pa1] K bulk modulus of liquid, [Pa] N number of reaches, [-] p liquid pressure, [Pa] Rt retardation times, [s] t time, [s]
tc[ valve closing time, [s] T temperature, [°C] u time used in convolution integral, [s] v liquid flow velocity, [m/s] W(t-u) friction weighting function,[-] wJ(t-u) creep weighting function, [Pa^s1] x distance along the pipe axis, [m] se instantaneous elastic strain, [-] sr retarded strain, [-] H dynamic viscosity of liquid, [Pa-s] v kinematic viscosity, [m2/s] g pipe support coefficient, [-] S expanded pipe support coefficient, [-] p liquid density, [kg/m3] a circumferential stress, [Pa] t wall shear stress, [Pa] At constant time step in the MOC, [s]
7 REFERENCES
[1]	Tison, G. (1958). Le mouvement non permanent succédant à l'ouverture d'une vanne sur une conduite en polyethylene (The unsteady flow produced by the opening of a valve on a polyethylene pipe). Becetel, communication N°3, Ghent University, Ghent. (in French)
[2]	Hardung, V. (1962). Propagation of pulse waves in visco-elastic tubings. Handbook of Physiology, (Field, J., ed.) p. 107135, American Physiological Society, Washington.
[3]	Martin, T. (1969). Pulsierende Strömung in visko-elastischen Leitungssystemen. Ingenieur-Archiv, vol. 37, no. 5, p. 315325, D0l:10.1007/BF00532579. (in German)
[4]	Lorenz, J., Zeller, H. (1972) An analogous treatment of waves propagation in liquid filled elastic tubes and gaz filled rigid tubes. Proceedings of the 1st International Conference on Pressure Surges, Paper B5.
[5]	Rieutord, E., Blanchard, A. (1972). Influence d'un comportement viscoélastique de la conduite dans le phénomène du coup de bélier (Influence of a viscoelastic pipe behavior in the phenomenon of water hammer). Reports of
the Academy of Sciences, Paris, vol. 274, p. 1963-1966. (in French)
[6]	Güney, M.S. (1977). Contribution à l'étude du phénomène de coup de bélier en conduite viscoélastique (Contribution to the study of the water hammer phenomenon in viscoelastic pipe), (Phd Thesis). Claude Bernard University, Lyon. (in French)
[7]	Zielke, W. (1968). Frequency-dependent friction in transient pipe flow. Journal of Basic Engineering, vol. 90, no. 1, p. 109115, D0l:10.1115/1.3605049.
[8]	Rieutord, E., Blanchard, A. (1979). Pulsating viscoelastic pipe flow - water-hammer. Journal of Hydraulic Research, vol. 17, no. 3, , p. 217-229, D0I:10.1080/00221687909499585.
[9]	Limmer, J., Meißner, E. (1974). Druckwellen in visko-elstischen Leitungen. In Mitteilungen des Institutes für Hydraulik und Gewässerkunde, Technische Universität München, München, vol. 14. (in German)
[10]	Franke, P.G., Seyler, F. (1983) Computation of unsteady pipe flow with respect to visco-elastic material properties. Journal of Hydraulic Research, vol. 21, no. 5, p. 345-353, D0I:10.1080/00221688309499456.
[11]	Ghilardi, P., Paoletti, A. (1987). Viscoelastic parameters for the simulation of hydraulic transients in polymeric pipes. Excerpta, vol. 2, p. 51-62.
[12]	Pezzinga, G., Scandura, P. (1995). Unsteady flow in installations with polymeric additional pipe. Journal of Hydraulic Engineering, vol. 121, no. 11, p. 802-811, D0I:10.1061/(ASCE)0733-9429(1995)121:11(802).
[13]	Tijsseling, A.S. (1996). Fluid-structure interaction in liquid-filled pipe systems: a review. Journal of Fluids and Structures, vol. 10, no. 2, p. 109-146, D0I:10.1006/jfls.1996.0009.
[14]	Covas, D., Stoianov, I., Ramos, H., Graham, N., Maksimovic, C., Butler, D. (2004). Water hammer in pressurized polyethylene pipes: conceptual model and experimental analysis. Urban Water Journal, vol. 1, no. 2, p. 177-197, D0I:10.1080/15730 620412331289977.
[15]	Covas, D., Stoianov, I., Ramos, H., Graham, N., Maksimovic, C. (2004). The dynamic effect of pipe-wall viscoelasticity in hydraulic transients. Part I - experimental analysis and creep characterization. Journal of Hydraulic Research, vol. 42, no. 5, p. 516-530, D0I:10.1080/00221686.2004.9641221.
[16]	Covas, D., Stoianov, I., Ramos, H., Graham, N., Maksimovic, C. (2005). The dynamic effect of pipe-wall viscoelasticity in hydraulic transients. Part II - model development, calibration and verification. Journal of Hydraulic Research, vol. 43, no. 1, p. 56-70, D0I:10.1080/00221680509500111.
[17]	Duan, H.-F., Ghidaoui, M., Lee, P.J., Tung, Y.-K. (2010). Unsteady friction and visco-elasticity in pipe fluid transients. Journal of Hydraulic Research, vol. 48, no. 3, p. 354-362, D0I:10.1080/00221681003726247.
[18]	Brunone, B., Berni, A. (2010). Wall shear stress in transient turbulent pipe flow by local velocity measurement. ASCE Journal of Hydraulic Engineering, vol. 136, no. 10, p. 716-726, D0I:10.1080/00221681003726247.
[19]	Bergant, A., Hou, Q., Keramat, A., Tijsseling, A.S. (2013). Waterhammer tests in a long PVC pipeline with short steel end sections. Journal of Hydraulic Structures, vol. 1, no. 1, p. 2334, D0I:10.22055/JHS.2013.10069.
Transient Liquid Flow in Plastic Pipes
89
Strojniški vestnik - Journal of Mechanical Engineering 66(2020)2, 77-90
[20]	Keramat, A., Tijsseling, A.S., Hou, Q., Ahmadi, A. (2012). Fluid-structure interaction with pipe-wall viscoelasticity during waterhammer. Journal of Fluids and Structures, vol. 28, p. 434-455, DOI:10.1016/jjflmdstracts.2011.11.001.
[21]	Pezzinga, G., Brunone, B., Meniconi, S. (2016). Relevance of pipe period on Kelvin-Voigt viscoelastic parameters: 1D and 2D inverse transient analysis. Journal of Hydraulic Engineering, vol. 142, no. 12, DOI:10.1061/(ASCE)HY.1943-7900.0001216.
[22]	Meniconi, S., Brunone, B., Ferrante, M., Massari, C. (2014). Energy dissipation and pressure decay during transients in viscoelastic pipes with an in-line valve. Journal of Fluids and Structures, vol. 45, p. 235-249, DOI:10.1016/j. jfluidstructs.2013.12.013.
[23]	Weinerowska-Bords, K. (2015). Alternative approach to convolution term of viscoelasticity in equations of unsteady pipe flow. Journal of Fluids Engineering, vol. 137, no. 5, art. ID: 054501, DOI:10.1115/1.4029573.
[24]	Ferrante, M., Capponi, C. (2017). Comparison of viscoelastic models with a different number of parameters for transient simulations. Journal of Hydroinformatics, vol. 20, no. 1, p. 1-17, DOI:10.2166/hydro.2017.116.
[25]	Kodura, A. (2016). An analysis of the impact of valve closure time on the course of water hammer. Archives of HydroEngineering and Environmental Mechanics, vol. 63, no. 1, p. 35-45, DOI:10.1515/heem-2016-0003.
[26]	Ferrante, M., Capponi, C. (2017). Viscoelastic models for the simulation of transients in polymeric pipes. Journal of Hydraulic Research, vol. 55, no. 5, p. 599-612, DOI:10.1080/ 00221686.2017.1354935.
[27]	Stryczek, J., Banas, M., Krawczyk, J., Marciniak, L., Stryczek, P. (2017). The fluid power elements and systems made of plastics. Procedia Engineering, vol. 176, p. 600-609, DOI:10.1016/j.proeng.2017.02.303.
[28]	Urbanowicz, K. (2018). Fast and accurate modelling of frictional transient pipe flow. Zeitschrift für Angewandte Mathematik und Mechanik, vol. 98, no. 5, p. 802-823, DOI:10.1002/zamm.201600246.
[29]	Vardy, A.E., Brown, J.M.B. (2010). Evaluation of unsteady wall shear stress by Zielke's method. Journal of Hydraulic Engineering, vol. 136, p. 453-456, DOI:10.1061/(ASCE) HY.1943-7900.0000192.
[30]	Urbanowicz, K., Firkowski, M., Zarzycki, Z. (2016). Modelling water hammer in viscoelastic pipelines: short brief. Journal of Physics: Conference Series, vol. 760, no. 1, art. ID: 012037, DOI:10.1088/1742-6596/760/1/012037.
[31]	Schohl, G.A. (1993). Improved approximate method for simulating frequency - dependent friction in transient laminar flow. Journal of Fluids Engineering, vol. 115, no. 3, p. 420424, DOI:10.1115/1.2910155.
[32]	Urbanowicz, K., Flrkowskl, M. (2018) Modelling water hammer with quasi-steady and unsteady friction in viscoelastic pipelines. Dynamical Systems in Applications, Awrejcewicz J. (ed.) Springer Proceedings in Mathematics & Statistics, vol. 249, p. 385-399, Springer, Cham, DOI:10.1007/978-3-319-96601-4_35.
[33]	Urbanowicz, K., Firkowski, M. (2018). Effect of creep compliance derivative in modeling water hammer in viscoelastic pipes. Proceedings of 13th International Conference Pressure Surges, p. 305-324.
[34]	Vardy, A.E., Brown, J.M.B. (2003). Transient turbulent friction in smooth pipe flows. Journal of Sound and Vibration, vol. 259, no. 5, p. 1011-1036, DOI:10.1006/jsvi.2002.5160.
[35]	Urbanowicz, K. (2017). Analytical expressions for effective weighting functions used during simulations of water hammer. Journal of Theoretical and Applied Mechanics, vol. 55, no. 3, p. 1029-1040, DOI:10.15632/jtam-pl.55.3.1029.
[36]	Adamkowski, A., Lewandowski, M. (2006). Experimental examination of unsteady friction models for transient pipe flow simulation. Journal of Fluids Engineering, vol. 128, no. 6, p. 1351-1363, DOI:10.1115/1.2354521.
[37]	Bertaglia, G., loriatti, M., Valiani, A., Dumbser, M., Caleffi, V. (2018). Numerical methods for hydraulic transients in viscoelastic pipes. Journal of Fluids and Structures, vol. 81, p. 230254, DOI:10.1016/j.jfluidstructs.2018.05.004.
[38]	Evangelista, S., Leopardi, A., Pignatelli, R., de Marinis, G. (2015). Hydraulic transients in viscoelastic branched pipelines. Journal of Hydraulic Engineering, vol. 141, no. 8, art. lD 04015016, DOI:10.1061/(ASCE)HY.1943-7900.0001030.
[39]	Covas, D. (2003). Inverse Transient Analysis for Leak Detection and Calibration of Water Pipe Systems Modelling Special Dynamic Effects. PhD Thesis, University of London, London.
[40]	Gally, M., Guney, M., Rieutord, E. (1979). An investigation of pressure transients in viscoelastic pipes. ASME Journal of Fluids Engineering, vol. 101, no. 4, p. 495-499, DOI:10.1115/1.3449017.
[41]	Firkowski, M., Urbanowicz, K., Duan, H.F. (2019). Simulation of unsteady flow in viscoelastic pipes. SN Applied Sciences, vol. 1, no. 519, DOI:10.1007/s42452-019-0524-2.
[42]	Swanson, C.J., Julian, B., Ihas, G.G., Donnelly, R.J. (2002). Pipe flow measurements over a wide range of Reynolds numbers using liquid helium and various gases. Journal of Fluid Mechanics, vol. 461, p. 51-60, DOI:10.1017/ S0022112002008595.
[43]	Stosiak, M., Zawislak, M., Nishta, B. (2018). Studies of resistances of natural liquid flow in helical and curved pipes. Polish Maritime Research, vol. 25, no. 3, p. 123-130, DOI:10.2478/pomr-2018-0103.
90
Urbanowicz, K. - Duan, H.-F. - Bergant, A.