Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 55, 3, pp. 1029-1040, Warsaw 2017 DOI: 10.15632/jtam-pl.55.3.1029 ANALYTICAL EXPRESSIONS FOR EFFECTIVE WEIGHTING FUNCTIONS USED DURING SIMULATIONS OF WATER HAMMER Kamil Urbanowicz West Pomeranian University of Technology, Szczecin, Poland e-mail: kamil.urbanowicz@zut.edu.pl For some time, work has been underway aimed at significant simplification of themodelling of hydraulic resistanceoccurring in thewaterhammerwhilemaintaininganacceptable error. This type of resistance is modelled using a convolution integral, among others, from local acceleration of a liquid and a certainweighting function.The recently completedwork shows that during efficient calculations of the convolution integral, the effectiveweighting function used does not have to be characterised by large convergencewith a classical function (accor- ding to Zielke during laminar flow and to Vardy-Brown during turbulent flow). However, it must be a sumof at least two or three exponential expressions so that the final results of the simulation could be considered as satisfactory. In this work, it has been decided to present certain analytical formulas using which it will be possible to determine the coefficients of simplified effective weighting functions in a simple direct way. Keywords: unsteady flow,water hammer, convolution integral, frequency-dependent friction 1. Introduction Unsteady flows occur in hydraulic systems, water supply systems, heating systems, thermal- -hydraulic systems (cooling cores of nuclear power plants), etc., during start-up, braking or failure. Propermodelling of flows of liquids under pressure in such systems remains a significant challenge. Among the key issues widely discussed in new publications on the subject, special emphasis is placed on the correctmodelling: of the time-varying hydraulic resistance (Vardy and Brown, 2003; Zarzycki et al., 2011; Reddy et al., 2012), cavitation (Zarzycki and Urbanowicz, 2006; Adamkowski and Lewandowski, 2009, 2012; Bergant et al., 2006; Karadžić et al., 2014; Soares et al., 2015), the interaction between the liquid andwalls of the conduit (Keramat et al., 2012; Henclik, 2015; Zanganeh et al., 2015), the viscoelastic phenomenon that occurs during the flow in a piping made of engineering polymers (Weinerowska-Bords, 2015; Soares et al., 2012; Keramat et al., 2013; Pezzinga et al., 2014; Urbanowicz et al., 2016). Taking into account all of the above phenomena while simulating unsteady flows has seemed impossible until recently. But now, thanks to the work carried out by Keramat and Tijsseling (2012) presented at the international conference on the analysis and damping of pressure surges associated with the phenomenon of water hammer (BHR Pressure Surges, Lisbon), we know that it is possible. In this generalmodel, aswell as inmany others having a simplified design, themethod ofmodelling of the time-varying hydraulic resistance has a very large impact on pressure runs. During acceleration, deceleration or as a result of rapid suppression of the fluid flow resulting fromquick valve closing (the so called water hammer effect occurs then), the friction of the fluid against pipe walls, as well as the internal friction between its elements, has a significant impact on transient flow parameters. It was already noted byRoiti, Helmholtz, Stearn andGromeka in thefirst studies concerningunsteadyfluidflow inpipes about150years ago.Arapiddevelopment of numerical methods in the 1950th and the 1960th, particularly development of the method 1030 K. Urbanowicz of characteristics being commonly used to date, induced further studies, the main objective of whichwas to properly describe the friction occurring during the flow in amathematical manner. At present, the models enabling the pipe wall shear stress to be simulated can be divided into two groups. The first group is simple models in which the stress is directly proportional to amomentary local and convective acceleration of the fluid. Themodel developed by a group of researchersunder the leadershipofDaily (1956), inwhich the stressdependedonlyonmomentary local acceleration of the fluid and a certain constant coefficient, is considered a prototype. The above model was developed with time by other researchers as Carstens and Roller, Safwat and Polder and Shuy and Apelt. A significant adjustment was introduced by Brunone et al. (1991), additionally making the stress conditional on momentary convective acceleration. Vı́tkovský et al. (2000) introduced a sign next to the convective derivative, while Laurerio and Ramos (2003) made a final adjustment of that model consisting in the splitting of the single constant coefficient k into two new coefficients kt and kx, which are to be found next to adequate velocity derivatives τw(t)= λqρv|v| 8 + ρD 8 ( kt ∂v ∂t +kxc |v| v ∣ ∣ ∣ ∂v ∂x ∣ ∣ ∣ ) (1.1) where: λq is the quasi-steady friction coefficient, ρ – liquid density, kt and kx – empirical coef- ficients, D – pipe inner diameter, v – velocity, t – time, c – pressure wave speed, x – distance along the pipe. Ramos then numerically proved the impact of particular expressions of this solution on the phase shifts and the speed of pressure wave damping, whereas Reddy et al. (2012), based on knownexperimental results, presentedamethodconsisting in the empirical selection of constants whencalculating the coefficients kt andkx that are tobe found infinal solution (1.1). Themodels of the above group are limited due to the need for empirical determination of the coefficients kt and kx. There are no papers thatwould showdetails of their numerical implementation; besides, they are characterized by a limited qualitative compatibility of pressure course being modelled (Adamkowski and Lewandowski, 2006), which is their major disadvantage. The second group of models consists of theoretical models being based on the so called convolution integral. The author of their prototype was Zielke (1968) who postulated τw(t)= 4µ R v+ 2µ R t∫ 0 w(t−u) ∂v ∂t (u) du (1.2) where: µ is the dynamic viscosity coefficient, R – pipe inner radius,w(t) – weighting function. The convolutional integral, being a product of the weight functionw(t) and themomentary value of fluid acceleration, is the inverse Laplace transform from the expression describing the impedance of a hydraulic line. In laminar flows, this impedance is being calculated froma simple analytical formula introduced by Brown, whereas in the turbulent ones it has a very complex analytical and empirical form (empirical because an empirical distribution of the coefficient of turbulent viscosity in the pipe cross-section is needed to resolve it), the derivation of which was reached at the same time by Zarzycki (1997, 2000) and Vardy and Brown (1996, 2003, 2004). The first numerical resolution of this integral being suitable for implementation in the me- thod of characteristics was already shown in thework byZielke (1968). Unfortunately, itwas not suitable for effective calculations, therefore a few years later Trikha (1975) presented another numerical procedure based on a three term weighting function being constructed from expo- nential terms. Unfortunately, also Trikha made too many simplifications, thus – with time – other authors presented their revised versions of that procedure (Kagawa et al., 1983; Schohl, 1993). Recently, Vardy and Brown (2010) noticed and corrected a significant error in the ori- ginal procedure according to Zielke, consisting in approximation instead of integration of the Analytical expressions for effective weighting functions... 1031 weight function in the dimensionless time interval from 0 to ∆t̂, but they did not present a re- vised effective calculation procedure. Such a procedure was, however, presented by Urbanowicz (2015). It is also worth emphasising that Zarzycki (1997, 2000) as well as Vardy and Brown (1996, 2003, 2004) proved that solution in form of equation (1.2) may be also used for the modelling of turbulent unsteady flows, provided that an adequate weight function was going to be used. Because all effective solutions are based on the weighting functions being a finite sum of exponential terms, the authors of effective numerical solutions frequently showed new forms of those functions in their papers referring to themodelling of laminar flow (Trikha 1975, Kagawa et al., 1983; Schohl, 1993; Vı́tkovský et al., 2004) or the turbulent one (Vı́tkovský et al., 2004; Zarzycki et al., 2011). Up to this day, themost accurate functions representedwith an extended range of use were presented by Urbanowicz and Zarzycki (2012). They are very useful in all cases that require a complete weighting function, for example in themodeling of one-directional accelerated or decelerated flows. The coefficients describing the effective weighting functions in turbulent flow are closely dependent on the Reynolds number and the internal roughness of the pipe walls. For correct determination, the classical scaling procedure developed byVı́tkovský et al. (2004) can be used, or the universal procedure (Urbanowicz et al., 2012). The advantage of them is providing the shape of the effective weighting function compatible with the shape of the classical laminarweighting function presented byZielke (1968) for the critical Reynolds number. In this paper, analytical formulas that enable coefficients describing simplified effective we- ighting functions composed of two or three terms to be determined are presented.Themethod is responsible for offloading computermemoryandaccelerating the iterative computational process without losing accuracy. The examplary results of simulation tests presented confirmhigh com- patibility of the simulated courses (with the use of the weighting function with limited ranges and the same being characterized by a simple structure) with the experimental ones. 2. New idea Recently completed studies have shown that unsteady flows can be modelled accurately using simplified effective weighting functions consisting of only two k=2 or three k=3 exponential expressions (Urbanowicz, 2015; Urbanowicz and Zarzycki, 2015) w(t) = k ∑ i=1 mie −nit̂ (2.1) wheremi andni are coefficients of effectiveweighting function, t̂ is thedimensionless time.These expressions are combined with the new improvedmethod for calculating shear stress. Functions used in the studies are characterised by a limited yet essential range of application (Fig. 1). The lower end of this range in the general case is set equal to the dimensionless time step∆t̂, and the upper end to themultiplicity there of 103∆t̂. The aforementioned time step in numerical calculations is calculated individually for all pressure conduits that retain their shape stability using the formula ∆t̂= L f ν cR2 (2.2) where:L is conduit length,ν –kinematic viscosity of liquid, f –numberof analysed cross-sections. Calculation of thecoefficientsmi andnidescribingthepreviouslyanalysed simplifiedeffective weighting functions required the use of a numerical method developed in 2012 (Urbanowicz, 2012). Elimination in this work of the numerical procedure mentioned above at the stage of 1032 K. Urbanowicz Fig. 1. Significant range of weighting functions determining the weighting function coefficients reduces time needed for numerical computation, facilitating the modelling of unsteady flows to perform a simple verification of the effectiveness of the method presented in the work of Urbanowicz and Zarzycki (2015) and enabling simple implementation of this method in the existing commercial software by introducing a variable hydraulic resistance coefficient λ(t+∆t) =λq,(t+∆t) (2.3) + 16ν R|v(t+∆t)|v(t+∆t) j ∑ i=1 [yi(t)Ai+ηBi(v(t+∆t)−v(t))+(1−η)Ci(v(t)−v(t−∆t))] ︸ ︷︷ ︸ yi(t+∆t) ︸ ︷︷ ︸ λu,(t+∆t) In the equation above, constantsAi,Bi andCi depend only on coefficientsmi andni describing the effective weighting function and the dimensionless time step Ai = e −ni∆t̂ Bi = mi ni∆t̂ (1−Ai) Ci =AiBi (2.4) 3. Analytical approximate solution Difficulties in widespread use of the simplified methodology presented in the work of Urbano- wicz and Zarzycki (2015), arising from the need to use the numerical procedure (Urbanowicz, 2012), may discourage practical use of the solutions discussed. Therefore, to further simplify themodelling of unsteady resistance, analytical solutions will be presented that can help one to accurately calculate coefficients mi and ni of simplified effective weighting functions as a func- tion of the dimensionless time step ∆t̂. The studies carried out previously (Urbanowicz, 2015; Urbanowicz andZarzycki, 2015) indicate that the relative percentage error of effective weighting functions should not exceed 30% for the two-expression functions or 10% for three-expression functions and that the optimal range of applicability of these functions should be from ∆t̂ to 103∆t̂. As can be seen in equation (2.2), the dimensionless time step ∆t̂ which is the starting point for the applicability of effective weighting functions assumes different values depending on the properties of flowing liquid, the conduit and the numerical method. Based on the analysis of practical and theoretical examples, the possible range of its variability can be specified using the domain of∆t̂∈ [10−10;10−1]. To determine the analytical function describing variation of the values of coefficients of two-expression effective weighting functions, it is necessary to first identify the set of values of these coefficients. For this purpose, themethod discussed in 2012 was used (Urbanowicz, 2012). Analytical expressions for effective weighting functions... 1033 When searching for an analytical solution for the effective two-expression functions, 93 sets of coefficients were determined (m1,m2,n1,n2) for the range from 10 −10 every 100.1 to 10−0.8. It is worth noting that for ∆t̂= 10−0.8, the values of these coefficients coincided with the values known from the classical weighting function for laminar flow (i.e. weighting according to Zielke (1968)), i.e.m1 =1,m2 =1,n1 =26.3744,n2 =70.8493. In the case of effective three-expression functions, 89 sets of coefficientsweredetermined (m1,m2,m3,n1,n2,n3) for the range from10 −10 every 100.1 to 10−1.2. The difference in the number of these sets resulted from the fact that in this case, already for ∆t̂ = 10−1.2, the values of these coefficients coincided with the values known from the classical weighting function for laminar flow, i.e. m1 = 1, m2 = 1, m3 = 1, n1 = 26.3744, n2 = 70.8493, n3 = 135.0198. Knowing all the above values, the next step was to adopt appropriate forms of analytic functions, which would accurately describe variability of these coefficients as a function of the dimensionless time step. Analysis of the variability of individual coefficients and the tests performed with other forms revealed that in the case of two-expression functions, their coefficients can be described using the formula mi,ni = 3 ∑ i=1 Ai∆t̂ Bi +C (3.1) in the range of their linearity on a log-log graph (Fig. 2a) (interval for ∆t̂ from 10−10 to 10−4, exceptionally for n1 to 10 −5). And the formula mi,ni = 4 ∑ i=1 Die −Ei∆t̂+F (3.2) for the range of their non-linearity on a log-log graph (Fig. 2b) (interval for ∆t̂ from 10−4 to∞, exceptionally forn1 from 10 −5 to∞). On the graphs presented below (Figs. 2 and 3), the abbreviation “sol.”means that these are the coefficients calculated using thepresentedanalytical formulas. Fig. 2. Compatibility of the analytical solution – two-expression functions To find definitive values of the coefficient of the functions adopted above, i.e. A1, . . . ,A3, B1, . . . ,B3 andC, andD1, . . . ,D4,E1, . . . ,E4 andF, the Curve Fitting Toolboxmodule imple- mented inMATLABwas used. The values of the estimated final coefficients are summarised in Table 1 and 2. Analysis of the variability of individual coefficients m1, . . . ,m3 and n1, . . . ,n3 representing three-expression functions showed that the forms of analytical functions, which were assumed 1034 K. Urbanowicz Fig. 3. Compatibility of the analytical solution – three-expression functions Table 1.Coefficients of the analytical solution of effective two-expression functions for the range of small dimensionless time steps (Eq. (3.1)) m1 m2 n1 n2 Coeff. Interval Interval Interval Interval [10−10;10−4] [10−10;10−4] [10−10;10−5] [10−10;10−4] A1 0.03234 0.1963 0.001476 0.09021 A2 48.35 2.88 0.1203 0.382 A3 9.717 −0.2661 526.7 223.1 B1 −0.5 −0.5 −1 −1 B2 0.5437 3.575 −0.5 −0.4592 B3 3.85 5.276 0.5567 0.2615 C −1.318 −0.2351 6.091 0 Table 2.Coefficients of the analytical solution of effective two-expression functions for the range of large dimensionless time steps (Eq. (3.2)) m1 m2 n1 n2 Coeff. Interval Interval Interval Interval (10−4;∞) (10−4;∞) (10−5;∞) (10−4;∞) D1 0.1480 2.214 9.317 56.56 D2 0.3227 4.155 87 136.5 D3 0.8039 7.929 188.1 396.7 D4 2.458 20.485 477.43 1903.3 E1 188.8 62.02 4459 79.71 E2 1316 386.6 29320 489.6 E3 5728 2191 104300 2880 E4 19270 12570 290500 15760 F 1 1 26.3744 70.8493 in the case discussed above, would also work. Thus, in range of its linearity on a log-log graph (Fig. 3a), the function sought has the form mi,ni = 4 ∑ i=1 Ai∆t̂ Bi (3.3) whereas for the range of non-linearity (Fig. 3b), we can describe it using the form exactly the same as in equation (3.2). In the case of the aforementioned analytical solution describing Analytical expressions for effective weighting functions... 1035 coefficients of the effective three-expression weighting functions, much greater volatility of the dimensionless time was noted at which the functions describing individual coefficients correctly passed from the linear form (in log-log scale) to the non-linear form. Specific times of transition from one form to another are shown in Tables 3 and 4. Table 3. Coefficients of the analytical solution of effective three-expression functions for the range of small dimensionless time steps (Eq. (3.3)) m1 m2 m3 n1 n2 n3 Coeff. Interval Interval Interval Interval Interval Interval [10−10;10−4] [10−10;10−4] [10−10;10−4] [10−10;10−5] [10−10;10−4.4] [10−10;10−4.2] A1 0.02239 0.06549 0.2336 0.0009749 0.02208 0.3037 A2 −1.123 −0.1334 11.52 0.09783 0.1233 0.1641 A3 34.85 −2.54 −11.62 6.215 11.55 5.039 A4 2.114e+06 2559 7.868 887.8 2025 1.011e+04 B1 −0.5 −0.5 −0.5 −1 −1 −1 B2 0 0 0 −0.5 −0.5 −0.5 B3 0.5138 0.2948 0.0002657 0.001247 0.001441 −0.07303 B4 1.789 2.894 3.297 0.5838 0.6193 0.6172 Table 4. Coefficients of the analytical solution of effective three-expression functions for the range of large dimensionless time steps (Eq. (3.2)) m1 m2 m3 n1 n2 n3 Coeff. Interval Interval Interval Interval Interval Interval (10−4;∞) (10−4;∞) (10−4;∞) (10−5;∞) (10−4.4;∞) (10−4.2;∞) D1 0.02449 0.8285 3.272 1.16 26.05 216 D2 0.06897 1.547 6.819 25.91 71.93 729.2 D3 0.2359 2.776 13.42 96.44 263.8 2522 D4 1.8429 5.9004 22.9793 251.6091 1427 12006.2 E1 246 190.8 83.86 2939 314.5 140.2 E2 995.2 907.7 645.4 1.792e+04 2054 969.4 E3 4787 4112 3779 6.098e+04 1.09e+04 5460 E4 1.696e+04 1.608e+04 1.895e+04 2e+05 4.32e+04 2.803e+04 F 1 1 1 26.3744 70.8493 135.0198 The maximum values of relative percentage errors, which are represented by the effective weighting functions determined using the above analytical formulas, calculated with reference to the classical function according toZielke (1968) are illustrated in the chart below (Fig. 4).The graph shows that themaximumerror for thedimensionless time stepof∆t̂≈ 10−4 systematically decreases until reaching zero. Achieving the zero value is equivalent to overlaping of coefficients calculated using the analytical method with coefficients from the classical laminar weighting function according to Zielke. With the use of the analytical formulas presented in this Section, it is possible only to determine coefficients that describe effective laminar functions. In a situation where there is turbulent flow, these coefficients have to be rescaled in accordance with the procedure described in (Vı́tkovský et al., 2004; Urbanowicz et al., 2012) for, as we know, the form of a classical turbulent weighting function according to Vardy and Brown (2007) is highly dependent on the Reynolds number. 1036 K. Urbanowicz Fig. 4. Maximum relative percentage errors of weighting functions designated analytically 4. Example calculation results To examine the impact of this new effective weighting function procedure presented in the previous Section, comparative studies for purewater hammer have beenmade. Basic continuity (4.1)1 andmomentum (4.1)2 equations describing unsteady flow in a horizontal pipe are ∂p ∂t +ρc2 ∂v ∂x =0 ∂p ∂x +ρ ∂v ∂t + 2 R τw =0 (4.1) where: p is pressure, v –mean velocity in pipe cross-section. To derive the above equations, the following assumptions are made: flow in the pipe is assumedas one-dimensional and thevelocity distributionuniformover thepipe cross-section; the pipewalls and the fluid are assumed as linearly elastic (stress proportional to strain). Equations (4.1) have been solved using the well-known method of characteristics. In this paper, the results of comparisons for two significant simulated and experimentally obtained pressure runs are presented. The experimental data have been obtained in a copper pipeline at the IMP in Gdańsk by Adamkowski and Lewandowski (2006) and previously publi- shed. All the details of the experimental test rig and the numerical procedures input data are presented in Table 5. Table 5.Test rig details and input data for simulations L=98.11m, ρ=997.65kg/m3,D=0.016m, ν =9.493 ·10−7m2/s, f =32, e=0.001m, c=1300m/s v0 [m/s] Re0 [–] pr [Pa] Type of flow 0.066 1112 1.265 ·106 laminar 0.94 15843 1.264 ·106 turbulent In the numerical analyses beingmade, the dimensionless time step amounted to ∆t̂=∆t ν R2 =3.5 ·10−5 where:∆t=∆x/c=0.0024s and∆x=L/f =3.066m. For the above dimensionless step, coefficients of optimum simplified effective weighting func- tions have been determined with the use of the procedure presented in this paper, see Table 6. The coefficients for turbulent tests required the re-scaling. Thedetails referring to the scaling procedurewere discussed in the papers byVı́tkovský et al. (2004) andUrbanowicz et al. (2012). Owingto the fact that thepipewalls areassumedtoberough(k=0.0000015 [m]), the coefficients Analytical expressions for effective weighting functions... 1037 Table 6.Calculated weighting function coefficients L=98.11m, ρ=997.65kg/m3,D=0.016m, ν =9.493 ·10−7m2/s, f =32, e=0.001m, c=1300m/s m1 m2 m3 n1 n2 n3 type of flow – no terms 4.333 32.954 – 70.45 2636 – laminar – 2 terms 2.864 10.816 39.43 52.92 666.9 8738 laminar – 3 terms 4.364 33.195 – 503.59 3069 – turbulent – 2 terms 2.885 10.895 39.72 486.05 1100 9171 turbulent – 3 terms mi and ni are scaled; the coefficients of exact weighting function according toVardy andBrown (2007) are used for scaling. The coefficients of effective weighting function with extended range of applicability (26 terms), which with high accuracy corresponds to the classical weighting function according to Zielke (numerical calculations in this paper were also made using this function), were previously discussed in the paper by Urbanowicz and Zarzycki (2012). The results of numerical tests obtained are illustrated in Figs. 5 and 6. Fig. 5. Results for laminar pipe flow (Re=1112) Themain conclusions from the comparisons presented above are as follows: a) Simplified modelling with a new weighting function constructed with only two exponen- tial terms is responsible for a gentle phase shift in the course being modelled. They are particularly visible in the final phase of flow deceleration (Figs. 5b and 6b). However, for the needs of engineering practice, the obtained results can be considered sufficient. b) Application of anewthree termweighting function, theapplicability rangeofwhich strictly dependson thehydraulic systemanalysed aswell as on thenumerical density of grid on the pipe length, allowed obtaining numerical results qualitatively compatible with the exact results of numerical tests (obtained using the exact extended 26 term weighting function and the full convolution based on the classical weighting function). 1038 K. Urbanowicz Fig. 6. Results for turbulent pipe flow (Re=15843) c) Phase shifts between the experimental results and the numerical ones shown inFigs. 5b, 5d and6b, 6d canbe explainedbyagentle variation in the speed of pressurewave propagation during recording of experimental courses. This variation of results may result from the impact of non-dissolved gases (air) found in the experimental system. The simulation tests performed clearly show the impact of unsteady friction on the courses obtained as a result of examining the water hammer effect. The applied and very simplified weighting functions present themselves perfectly against the results obtained using only the quasi-steadymodel of friction. Thus, it is possible to safely recommend the presented procedure for engineers who are involved in protection of hydraulic systems against negative effects of the water hammer. 5. Summary The analytical solutions presented in the paper allow one to quickly determine simplified forms of effective weighting functions composed of two or three exponential expressions. These corre- lations could be used in a simple manner by applying the instantaneous resistance coefficient (Eq. (2.3)) in commercial and custom computer programs used for the modelling of unsteady flows of liquids in conduits under pressure. The biggest problem associated with implementing the presented solution is the need to introduce into the program many constants estimated in this paper, which describe individual solutions. Another issue which the future user of the pre- sented formulas should pay attention to, is the right choice of themethod of the characteristics grid. With the range of application of effective weighting functions simplified in this manner, the number of computing sections should not be higher than f=50, because for this value, the instantaneous hydraulic resistance calculated is a function of velocity changes occurring in the last five periods of the water hammer. Analytical expressions for effective weighting functions... 1039 References 1. AdamkowskiA.,LewandowskiM., 2006,Experimental examinationofunsteady frictionmodels for transient pipe flow simulation, Journal of Fluids Engineering, ASME, 128, 6, 1351-1363 2. Adamkowski A., Lewandowski M., 2009, A new method for numerical prediction of liquid column separation accompanying hydraulic transients in pipelines, Journal of Fluids Engineering, ASME, 131, 7, paper 071302 3. Adamkowski A., Lewandowski M., 2012, Investigation of hydraulic transients in a pipeline with column separation, Journal of Hydraulic Engineering, ASCE, 138, 11, 935-944 4. Bergant A., Simpson A.R., Tijsseling A.S., 2006,Water hammer with column separation: A historical review, Journal of Fluids and Structures, 22, 2006, 135-171 5. BrunoneB.,GoliaU.M.,GrecoM., 1991, Some remarks on themomentumequations for fast transients,Proceedings of International Meeting on Hydraulic Transients with Column Separation, 9th Round Table, IAHR, Valencia, Spain, 201-209 6. Dailey J.W., Hankey W.L., Olive R.W., Jordaan J.M., 1956, Resistance coefficient for accelerated and decelerated flows through smooth tubes and orifices, Journal of Basic Engineering, ASME, 78, 1071-1077 7. Henclik S., 2015, A numerical approach to the standard model of water hammer with fluid- structure interaction, Journal of Theoretical and Applied Mechanics, 53, 3, 543-555 8. KagawaT., Lee I.,KitagawaA.,TakenakaT., 1983,High speedandaccurate computingme- thod of frequency-dependent friction in laminar pipe flow for characteristicsmethod (in Japanese), Transactions of the Japan Society of Mechanical Engineers, Part A, 49, 447, 2638-2644 9. Karadžić U., Bulatović V., Bergant A., 2014, Valve-induced water hammer and column separation in a pipeline apparatus, Strojnǐski vestnik – Journal of Mechanical Engineering, 60, 11, 742-754 10. Keramat A., Tijsseling A.S., Hou Q., Ahmadi A., 2012, Fluid-structure interaction with pipe-wall viscoelasticity during water hammer, Journal of Fluids and Structures, 28, 1, 434-456 11. Keramat A., Kolahi A.G., Ahmadi A., 2013, Waterhammer modelling of viscoelastic pipes with a time-dependent Poisson’s ratio, Journal of Fluids and Structures, 43, November, 164-178 12. Keramat A., Tijsseling A.S., 2012,Waterhammer with column separation, fluid-structure in- teraction and unsteady friction in a viscoelastic pipe,Proceedings of 11th International Conference on Pressure Surges, BHRGroup, Lisbon, Portugal, October 24-26, 443-460 13. Loureiro D., Ramos H., 2003, A modified formulation for estimating the dissipative effect of 1-d transient pipe flow, [In:] Pumps, Electromechanical Devices and Systems Applied to Urban Management, Cabrera E. et al. (Edit.), A.A. Baklema Publishers, The Netherlands, 2, 755-763 14. Pezzinga G., Brunone B., Cannizzaro D., Ferrante M., Meniconi S., Berni A., 2014, Two-dimensional features of viscoelastic models of pipe transients, Journal of Hydraulic Engine- ering, 140, 8, Art. No. 04014036 15. Reddy H.P., Silva-Araya W.F., Chaudhry M.H., 2012, Estimation of decay coefficients for unsteady friction for instantaneous, acceleration-basedmodels, Journal of Hydraulic Engineering, 138, 260-271 16. Schohl G.A., 1993, Improved approximatemethod for simulating frequency – dependent friction in transient laminar flow, Journal of Fluids Engineering, ASME, 115, September, 420-424 17. Soares A.K., Martins N.M.C., Covas D.I.C., 2012, Transient vaporous cavitation in visco- elastic pipes, Journal of Hydraulic Research, 50, 2, 228-235 18. Soares A.K., Martins N., Covas D.I.C., 2015, Investigation of transient vaporous cavitation: experimental and numerical analyses,Procedia Engineering, 119, 235-242 19. Trikha A.K., 1975, An efficient method for simulating frequency-dependent friction in transient liquid flow, Journal of Fluids Engineering, ASME, 97, 1, 97-105 20. Urbanowicz K., 2012, New approximation of unsteady friction weighting functions, Proceedings of 11th International Conference on Pressure Surges, Lisbon, Portugal, October 24-26, 477-492 1040 K. Urbanowicz 21. Urbanowicz K., 2015, Simplemodelling of unsteady friction factor,Proceedings of 12th Interna- tional Conference on Pressure Surges, Dublin, Ireland, 18-20 November, 113-130. 22. Urbanowicz K., Firkowski M., Zarzycki Z., 2016, Modelling water hammer in viscoelastic pipelines: short brief, Journal of Physics: Conference Series, 760, paper 012037 23. UrbanowiczK., ZarzyckiZ., 2012,Newefficient approximationofweighting functions for simu- lations of unsteady friction losses in liquid pipe flow,Journal of Theoretical andAppliedMechanics, 50, 2, 487-508 24. Urbanowicz K., Zarzycki Z., 2015, Improved lumping friction model for liquid pipe flow, Journal of Theoretical and Applied Mechanics, 53, 2, 295-305 25. Urbanowicz K., Zarzycki Z., Kudźma S., 2012, Universal weighting function in modeling transient cavitating pipe flow, Journal of Theoretical and Applied Mechanics, 50, 4, 889-902 26. Vardy A.E., Brown J.M.B., 1996, On turbulent unsteady, smooth – pipe friction, Proceedings of 7th International Conference on Pressure Surges, BHR Group, Harrogate, United Kingdom, 289-311 27. Vardy A.E., Brown J.M.B., 2003, Transient turbulent friction in smooth pipe flows, Journal of Sound and Vibration, 259, 5, 1011-1036 28. VardyA.E., Brown J.M.B., 2004,Transient turbulent friction in fully roughpipe flows, Journal of Sound and Vibration, 270, 233-257 29. Vardy A.E., Brown J.M.B., 2007, Approximation of turbulent wall shear stresses in highly transient pipe flows, Journal of Hydraulic Engineering, 133, 11, 1219-1228 30. Vardy A.E., Brown J.M.B., 2010, Evaluation of unsteadywall shear stress by Zielke’smethod, Journal of Hydraulic Engineering, 136, 453-456 31. Vı́tkovský J., Lambert M., Simpson A., Bergant A., 2000, Advances in unsteady friction modeling in transient pipe flow, Proceedings of 8th International Conference on Pressure Surges, BHRGroup, Hague, Holland, 471-482 32. Vı́tkovský J.P., StephensM.L.,BergantA., SimpsonA.R.,LambertM.F., 2004,Efficient and accurate calculation of Zielke and Vardy–Brown unsteady friction in pipe transients, Proce- edings of 9th International Conference on Pressure Surges, BHRGroup,Chester,UnitedKingdom, 405-419 33. Weinerowska-Bords K., 2015, Alternative approach to convolution term of viscoelasticity in equations of unsteady pipe flow, Journal of Fluids Engineering, ASME, 137, 5, paper 054501 34. Zanganeh R., Ahmadi A., Keramat A., 2015, Fluid-structure interaction with viscoelastic supports during waterhammer in a pipeline, Journal of Fluids and Structures, 54, April, 215-234 35. Zarzycki Z., Kudźma S., Urbanowicz K., 2011, Improvedmethod for simulating transients of turbulent pipe flow, Journal of Theoretical and Applied Mechanics, 49, 1, 135-158 36. Zarzycki Z., 1997,Hydraulic resistance of unsteady turbulent liquid flow in pipes,Proceedings of 3rd International Conference on Water Pipeline Systems, BHRGroup, Hague, Holland, 163-178 37. Zarzycki Z., 2000, On weighting function for wall shear stress during unsteady turbulent flow, Proceedings of 8th International Conference on Pressure Surges, BHRGroup, Hague, Holland, 39, 529-534 38. Zarzycki Z., Urbanowicz K., 2006,Modelling of transient flow during water hammer conside- ring cavitation in pressure pipes (in Polish),Chemical and Process Engineering, 27, 3, 915-933 39. Zielke W., 1968, Frequency-dependent friction in transient pipe flow, Journal of Basic Engine- ering, ASME, 90, 1, 109-115 Manuscript received November 23, 2015; accepted for print April 7, 2017