Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 53, 2, pp. 295-305, Warsaw 2015 DOI: 10.15632/jtam-pl.53.2.295 IMPROVED LUMPING FRICTION MODEL FOR LIQUID PIPE FLOW Kamil Urbanowicz, Zbigniew Zarzycki West Pomeranian University of Technology, Department of Mechanical Engineering and Mechatronics, Szczecin, Poland e-mail: kamil.urbanowicz@zut.edu.pl; zbigniew.zarzycki@zut.edu.pl Normally, during one-dimensional pipe flow, the friction terms are calculatedwith the use of a numericalmethod (for exampleMOC–method of characteristics) at every computational node along the pipe and at every time step. This procedure tends to increase the compu- tational effort greatly. A considerable increase in computational speed can be archived by calculating the frequency-dependent friction at the end of the pipe only. To avoid possible problems (no damping at closed walls, underestimate damping on high impedance compo- nents) the frequency-dependent friction term is calculated from the flowwaves.The lumping frictionmodel in this work is based on amodificated Schohl convolution integral solution. In addition, thework examined the impact of using of simplified effectiveweighting function on the obtained results of numerical simulations. Themodifiedmethod in conjunctionwith the use of simplified weighting function allow determination of real-time estimate of the basic parameters representing the fluid flow in complex hydraulic systems, water supply, etc. Keywords: transient pipe flow, frequency-dependent friction, effective numerical solution, water hammer, cavitation 1. Introduction In a transient liquid pipe flow velocity and pressure are changing with time. In fluid systems, transient occurs during starting-off or stopping of a pump, closing or opening a valve, or changes in tank levels. Often transient flow conditions persist as oscillating pressure and velocity waves for some time after the initial event that caused it. A water hammer is an example of that condition where the compressibility of the fluid has the dominant effect after a sudden closure of the valve. A transient flow can also result in significant transient pressures that may exceed the design limit of pipes and fittings, especially when the designer calculated it from a verywell and most known Joukovsky formula. Detailed understanding of the transient flow is then very important for safe operating of fluid systems. Thehereby study refers tonumericalmodellingof unsteadyflow inpressureconduits. It takes into account the friction loss coefficient as a total of quasi-steady and time-variable coefficient. The time-variable coefficient is calculated with the use of a new effective solution of convolution integral (Urbanowicz andZarzycki, 2012). The study comprises detailed comparison of numerical results with the use of the conventional method (hydraulic resistance calculated in all computer nodes as a total of two components) and the lumped method (where the resistance coefficient is determined in a simplified way – discussed in detail in the present study), see also Johnston (2006). The results presented in the study showwhether the simplifiedmethodwill be consistent enough to be successfully employed in practice. 296 K. Urbanowicz, Z. Zarzycki 2. Modelling of the coefficient of friction losses During unsteady flow of a liquid in pressure conduits, the value of the friction factor λ depends on the averaged flow velocity λ=λq+ 32 |Re|v S (2.1) whereλq is the quasi-steady friction coefficient, Re –Reynolds number, v – average current flow velocity. On the other hand, the sum S is determined from the below correlation (Urbanowicz and Zarzycki, 2012) S = j ∑ i=1 [yi(t)Ai+ηBi(v(t+∆t)−vt)+(1−η)Ci(vt−v(t−∆t))] ︸ ︷︷ ︸ yi(t+∆t) (2.2) where η is the corrective coefficient (for details, see Urbanowicz and Zarzycki, 2012), yi(t) – a time-variable parameterwhich is a value of the above total calculated for the previous time step, Ai,Bi andCi – constants which are more widely discussed later, v – average flow velocity in a proper time stop. In the laminar flow, the quasi-steady resistance coefficient (the subscript ql) does not depend on the roughness of the walls and is determined from the below correlation λql = 64 |Re| (2.3) On the other hand, in the turbulent flow (hence the subscript qt), the impact of the internal geometrical structure of the walls must not be neglected and it has influence on the resistance coefficient value (Colebrook, 1939) 1 √ λqt =−2log ( 2.51 Re √ λqt + ε/D 3.7 ) (2.4) where ε/D is the relative roughness of the pipeline internal walls. The above empirical Colebrook equation may be solved by numerical methods with high accuracy. There are also analytical approximations of the above formula which determines the friction coefficient in turbulent flows. The best approximation form was given by Goudar and Sonnad (2008) (for details see Appendix A). The research carried out recently has shown that a sufficient condition for correct (as ac- curately as required) modelling of hydraulic resistances occurring in course of unsteady states in conduits is constructing the sum S of just two terms. Constants Ai, Bi and Ci occurring in formula (2.2) for this total are calculated from the following correlations Ai =e −ni∆t̂ Bi = mi ∆t̂ni (1−Ai) Ci =AiBi (2.5) The above correlation shows that the constantsAi,Bi andCi depend on the dimensionless time step (constant in the method of characteristics) as well as coefficients mi and ni describing the so-called effective weighting function ∆t̂=∆t v R2 weffective = j ∑ i=1 mie −ni∆t̂ (2.6) Improved lumping friction model for liquid pipe flow 297 The integral of the aforementioned weighting functionmultiplied by the derivative of the liquid local velocity is defined as the so-called convolution integral. It is a result of a reverse Laplace transform of the correlation describing the impedance of the hydraulic system presented by Brown (1962). This integral describes the unsteady component of the wall shear stress. It was presented for the first time by Zielke (1968). It takes the following form τu = 2µ R t∫ 0 w(t−u) ∂v ∂t (u) du (2.7) Rquation (2.7)may obviously beused also for the unsteady turbulent flow (Zarzycki, 1997, 2000; Vardy and Brown, 2003, 2004) if the weighting function is known. The classic solution of the integral (2.7) employed in the method of characteristics is inef- fective because the calculation of the momentary value τu in the conduit cross section being analysed requires knowledge and use of the entire flow history. Trikha (1975) showed, however, that this solution can be presented in an effective formwhich allows avoiding the need to store and use the data referring to the average velocity changes occurring in time from the moment of occurrence of the unsteadiness in a specific point of the conduit. In the case of the Trikha method, it is sufficient to know the last two average flow velocities occurring in the current time moment at a specific pipe cross section. Subsequent researchers have perfected the Trikha solu- tion. There are known effective solutions byKagawa et al. (1983), Schohl (1993) or Urbanowicz and Zarzycki (2012). The recently carried out benchmark tests of these effective solutions have shown that (Urbanowicz-Zarzycki, 2012): • the Trikha solution was derived with too many simplifications and it is not suitable for the correct modelling of resistances, • the solution by Kagawa and others has its classical counterpart in the solution according to Zielke (1968), • theSchohl solutionhas its classical counterpart inZielke-Vardy-Brown solution (Vardyand Brown, 2010) on condition that the weighting function is expressed with a large number of exponential terms (which increases the time of calculation), • the Urbanowicz-Zarzycki solution is an effective counterpart of the Zielke-Vardy-Brown classical solution but it does not require expressing the effective weighting function with a large number of exponential terms. Further in the present study, a solution in form of (2.2) presented at Conference of Fluid Me- chanics (KKMP2012) held in Gliwice will be used (Urbanowicz-Zarzycki, 2012). Formula (2.1) shows that when lower velocities are dealt with, bigger is the influence of flow unsteadiness on themomentary value of the resistance coefficient. It can be stated that for large Reynolds numbers (Re> 106), the coefficient λu can be entirely negligible since λu ≪λq. 3. Research results The correctmodelling of hydraulic resistances occurring in hydraulic, water supply pipeline and heating systems in unsteady flow conditions is still not commonly employed in the design stage using specialised software. Most of the common computer programmes continue to be based on an ordinary quasi-steady approachwhich shouldnot beused in the case of unsteadyflows.There is no doubt that the developers of the commercial software dedicated to the modelling of flows are not interested in the correctmodelling of resistances because of themathematical complexity of the known correct solution. Therefore, this study has been dedicated to the comparison of 298 K. Urbanowicz, Z. Zarzycki the results of numerical research obtained from simplified effective methods of the modelling of the hydraulic resistances. The analysis of literary sources shows that there are two approaches which can possibly be adaptedwhile determining hydraulic resistances with the use of theMOC method of characteristics: • Standard Effective Method (SM) – the resistances are determined for all analysed cross sections (mostly, there are 16 cross-sections being analysed as the division into such a number ensures that the results obtained are sufficiently accurate) and all time steps with the use of effective solutions of the convolution integral which are based on the effective weighting function, which have been expressed in the previously presented form (2.6)2; • LumpedFrictionMethod (LFM) – the unsteady hydraulic resistances are determined only in the boundary cross-sections, while only the quasi-steady resistances are determined in the internal cross-sections (Johnston, 2006). According to the study by Johnston (2006), in the last of the aforementioned approaches, the momentary values of velocities in boundary nodes (Fig. 1) should be corrected in the following manner vRB,c = 1 2 ( vRB,r + pRB,r ρc ) vLB,c = 1 2 ( vLB,r+ pLB,r ρc ) (3.1) where vRB,c and vLB,c are subsequently, velocities corrected at the right and the left boundary, vRB,r and vLB,r – velocities calculated in the classical way on the right and the left boundary, ρ – fluid density, c – pressurewave propagation velocity, pRB,r and pLB,r – pressures at the right and the left boundary of the system, respectively. Fig. 1. Method of characteristic grid The corrected velocity values are used in numerical calculations of the unsteady corrected hydraulic resistance coefficient λu, Eqs. (2.1) and (2.2). The below study examines in detail not only the influence of themethod of calculation (the SM standardmethod or the LFM lumped frictionmethod) but also the correlation between the amount of terms constituting the effective function of weight and the quality of the results ob- tained. The detailed examination comprises forms of the effective weighting functions composed of the exponential terms 11, 4, 3, 2 – ver. I and 2 – ver. II (for detailed information about the coefficients, see Appendix B). The percentage errors of the adjustment of the above effective weighting functions with respect to Zielke (1968) classical weight are presented in Fig. 2. The results of thenumerical tests carried outwith theuse of theabovemethodsSMandLFM as well as various forms of the effective weighting function will be compared against the results of experimental tests of pressure change trends in course of hydraulic impact of water hammer. The experimental research was registered byAdamkowski and Lewandowski (2006, 2012) at the test facility (Fig. 3) situated at the Institute of Fluid-Flow Machinery of the Polish Academy Improved lumping friction model for liquid pipe flow 299 Fig. 2. Effective weighting function errors of Sciences (PAN) in Gdańsk. This study focused on four experimental courses obtained at the following initial velocities: • Re=1111, v0 =0.066m/s (Adamkowski and Lewandowski, 2006) • Re=6899, v0 =0.47m/s (Adamkowski and Lewandowski, 2012) • Re=25040, v0 =1.7m/s (Adamkowski and Lewandowski, 2012) • Re=40358, v0 =2.74m/s (Adamkowski and Lewandowski, 2012) Fig. 3. Layout of the test stand In the flow inwhich thephenomenonof liquid column separation occurs, commonly knownas vapour cavitation, pressurepulsations canbemuchgreater than thosewhichmight be calculated from the commonly used Joukovsky formula ∆p= ρc∆v (3.2) This is related to the occurrence of sudden accelerations of fluid at the moment when the mo- mentary vapour cavitation areas are closed (Simpson andWylie, 1991). Therefore, the turbulent courses during which liquid column separation took place play such an important role in this study. The simulation of unsteady flows with cavitation in the present study employed the use of the cavitation model according to Adamkowski and Lewandowski (2009) in which the artificial dampening taking place in the classical column separation model has been eliminated. The numerical courses obtained as a result of the simulation have been compared against the experimental courses by means of a simple quantitative method which has been used to 300 K. Urbanowicz, Z. Zarzycki determine time correctness – the parameter Et (times of occurrence of pressure maximums within subsequent amplitudes of the course of pressure – Fig. 4) Et= 1 k k ∑ i=1 ∣ ∣ ∣ tis− tie tie ∣ ∣ ∣ ·100% (3.3) where tis are – subsequent times for which a local pressure maximum occurred, tie – the same times in the experimental course, k – number of pressure amplitudes being analysed. Fig. 4. Analyzed pressures peaks The value correctness is defined by the parameter Ep (the analysis covered values of maxi- mum pressures on subsequent amplitudes – Fig. 4) Ep= 1 k k ∑ i=1 ∣ ∣ ∣ pis−pie pie ∣ ∣ ∣ ·100% (3.4) where pis are maximum values of simulated pressures on subsequent pressure amplitudes, pie – maximum values of pressures on subsequent amplitudes recorded during the experiment. The Tables 1-4 include all obtained results of the quantitative analysis. Furthermore, the smallest errors in these tables havebeenmarkedbold face tomake themattract special attention. Tabel 1.Errors parameters (for laminar flow v0 =0.0066m/s) Type of Standard unsteady Lumped unsteady weighting friction – SM friction – LFM function Et Ep Et Ep 2 terms – ver. I 0.693314 0.318876 0.865383 0.292364 2 terms – ver. II 0.693426 0.560728 0.88126 0.430243 3 terms 0.703865 0.409423 0.88045 0.465801 4 terms 0.702829 0.411929 0.875647 0.514488 11 terms 0.702829 0.352031 0.863143 0.30647 3.1. Laminar flow results Theobtained results of thenumerical tests for laminarflow show(Table 1) that the simulated changes of pressure by the SM standardmethod and the LFMmethodwith unsteady resistance Improved lumping friction model for liquid pipe flow 301 Tabel 2.Errors parameters (for turbulent flow v0 =0.45m/s) Type of Standard unsteady Lumped unsteady weighting friction – SM friction – LFM function Et Ep Et Ep 2 terms – ver. I 1.700704 17.92573 1.010862 15.49663 2 terms – ver. II 2.859451 17.6086 2.333399 21.32954 3 terms 1.977994 16.17788 1.177909 15.0219 4 terms 2.463139 17.34088 1.784986 10.85285 11 terms 2.13715 15.41202 1.628331 13.00306 Tabel 3.Errors parameters (for turbulent flow v0 =1.7m/s) Type of Standard unsteady Lumped unsteady weighting friction – SM friction – LFM function Et Ep Et Ep function Et Ep Et Ep 2 terms – ver. I 2.979111 3.342316 2.269701 5.367631 2 terms – ver. II 2.696019 4.498614 1.85463 6.653467 3 terms 2.733686 4.344046 2.761818 5.027982 4 terms 2.329251 4.200977 2.600806 3.812173 11 terms 2.860539 3.923151 2.419137 6.48246 Tabel 4.Errors parameters (for turbulent flow v0 =2.76m/s) Type of Standard unsteady Lumped unsteady weighting friction – SM friction – LFM function Et Ep Et Ep 2 terms – ver. I 2.123979 2.606999 2.126099 6.09562 2 terms – ver. II 2.478371 3.962234 2.292239 3.336577 3 terms 2.309905 3.883457 2.177465 4.336589 4 terms 2.20718 3.934879 2.740817 1.863973 11 terms 2.20718 3.835296 1.87071 4.136423 lumped in boundarynodes (lumped frictionmodel) featured similar accuracy. The average value (computedbasedon the results obtained for all five effective functions ofweigh underanalysis) of the parameterEp defining the degree of conformity of the simulated pressure peaks in reference to the experimental results employing the use of the SM and LFM methods was very close to EpSM = 0.41% and EpLFM = 0.40%. The time coincidence of the pressure peaks determined by means of the parameter Et, the second parameter calculated in the quantitative method, is slightly worse for the average results obtained bymeans of the LFMmethod (EtSM =0.7% and EtLFM =0.87%). What is interesting is that the best result of the parameterEp is obtainedwith the use of the LFM method and an effective function of weigh composed of just two terms (2 terms – ver. I). This confirms the fact that, in the case of laminar flows, it is possible to adjust the effective weighting function in such away that the obtained resultswill coincidewith those obtainedwith the use of very accurate effective functions ofweight (composed of a large number of exponential terms and therefore requiringmanymore calculations in order to estimate themomentary local resistances). 302 K. Urbanowicz, Z. Zarzycki 3.2. Turbulent flow results a) v0 =0.45m/s The average results calculated from Table 2 shows that the time coincidence of pressure peaks (parametr Et) occurring while using the SM standard method is insignificantly worse (EtSM =2.23%) than in the case of the LFMmethod (EtLFM =1.59%). Evenmore interesting results have been obtained for the Ep parameter which is more important from the practical point of view. As amatter of fact, it can be seen that a similar error has been obtained by using the SM method with a very accurate function of weight (11 terms) Ep(SM−11terms) =15.4% as by using the LFMmethod and a very inaccurate function of weight consisting of just two terms Ep(LFM−2terms ver.I) =15.5%. The course bearing the smallest error is the one simulated by means of the LFM method with the use of a weighting function composed of 4 exponential termsEp(LFM−4terms) =10.85% (Fig. 5). It is also worth noticing that the analysis of the first two results obtained by means of the LFMmethod with the use of weighting functions composed of just two terms shows that the influence of the course of the effective weighting function is crucial; therefore, the span of the result Ep(LFM−2terms ver.II)−Ep(LFM−2terms ver.I) =5.83% obtained between these two examinations is rather large. At the present stage of the research, it is, however, difficult to determine how to optimally select the shape of the effective weighting function and onwhat to focus at the stage of determi- ning the coefficients which describe it, so as to ensure that, with a small number of terms, the simplified numerical results are represented with an error similar to that of the accurate results. Fig. 5. Pressure variation results at the valve (v0 =0.45m/s, 4 terms eff. weighting function) b) v0 =1.7m/s In the analysed case, slightly better average results determining the time coincidence of the occurrence of the pressure course maximum values have been obtained with the use of the simplified LFM method: EtSM = 2.72% and EtLFM = 2.38%. The average values of the parameterEp for the two analysedmethods have shown that the SM standardmethod is better at simulating maximum pressures for higher Reynolds numbers: EpSM = 4.06% and EpLFM = 5.47%. Itmust not be left unnoticed that also in this case a very good results has been obtained bymeansof the simplifiedLFMmethodand theweighting function composedof four exponential terms: Ep(LFM−4terms) = 3.81%. The result obtained is lower than the average results, which Improved lumping friction model for liquid pipe flow 303 confirms that if the effective weighting function is selected properly, it is possible to obtain a very good course of pressure which is coincident with the experiment by using the simplified LFMmethod. c) v0 =2.76m/s The last case of the turbulent flow to be analysed concerns the flow with an initial velocity of v0 =2.76m/s. The average values of the parametersEp andEt for the two analysedmethods are:EpSM =3.64%,EpLFM =3.95%,EtSM =2.26%andEtLFM =2.24%.So, theaverage results are comparable. This is confirmed by the fact noticed in the course of the analysis of formula (2.1)which becomes more accurate as the average flow velocity increases and the influence of the unsteady resistance coefficient λu decreases. The best results of the parameterEp are again obtainedwith theuse of aweighting function composedof 4 exponential termsEp(LFM−4terms) = 1.86% (Fig. 6), which finally confirms the possibility of commonly using the simplified LFM method while modelling the unsteady hydraulic resistances in pressure conduits. Fig. 6. Pressure variation results at the valve (v0 =2.76m/s, 4 terms eff. weighting function) 4. Conclusions The presented study compares the research results obtained with the use of the two methods beinganalysed. In thefirst of thesemethods, SM(the standardmethod), thehydraulic resistance coefficient has been determined in all computational nodes of the methods of characteristics with the use of formula (2.1). In the second method, the so-called lumpedmethod (LFM), the resistance coefficient has been determinedwith the use of the same formula (2.1), but the sumS found in formula (2.2) is calculated only in the case of boundary nodes while in all internal nodes S =0 is implied (that is, as a matter of fact, the hydraulic resistances are determined in these internal nodes in a quasi-steady way). The results obtained in the study have shown that the simplified LFMmethod can be used for the modelling of transient pressure courses with a high accuracy. As a matter of fact, using this method, features a degree of accuracy which is close to that obtained when using amethod causing much higher computer workload. The results also show that it is possible to achieve very good results in laminar flows by using a very simple effective weighting function composed of just two exponential terms. In the case when the LFMmethod is used for turbulent flow, the best results would be obtained for an effective weighting function composed of four terms. On the other hand, in the SM standardmethod, very good results are obtained in the modelling of turbulent flows when simplest weighting function composed of two terms is used. The results obtained for a flow qualified as turbulent and featuring a small Reynolds number of Re=6899 show too high error level (Ep≈ 15%). It is difficult to clearly determine the cause 304 K. Urbanowicz, Z. Zarzycki of such a significant error. It is possible that this is caused by the fact that the Adamkowski and Lewandowski model used for the modelling of sizes of cavitation areas does not take into account the influence of gas cavitation or the nature of unsteady resistances occurring in the course of the transition flow (2320