Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 50, 2, pp. 487-508, Warsaw 2012 50th Anniversary of JTAM NEW EFFICIENT APPROXIMATION OF WEIGHTING FUNCTIONS FOR SIMULATIONS OF UNSTEADY FRICTION LOSSES IN LIQUID PIPE FLOW Kamil Urbanowicz Zbigniew Zarzycki West Pomeranian University of Technology, Faculty of Mechanical Engineering and Mechatronics, Szczecin, Poland; e-mail: kamil.urbanowicz@zut.edu.pl; zbigniew.zarzycki@.zut.edu.pl Most papers dealing with calculations and simulations of the unsteady liquid pipe flow are based on the assumption that the formula for quasi- steady friction (Darcy-Weisbach formula) can be applied. In the case of fast changes, like fast transients e.g. water hammer, it fails. In this work, the wall shear stress is presented as a sum of quasi-steady and unsteady component. The unsteady component of thewall shear stress ismodeled as a convolution of local fluid acceleration and aweighting function. The originalweighting functionhas usually avery complicated structure, and what is more, makes impossible to carry out an efficient simulation of dynamical runs. In this paper, in order to enable efficient calculation of the unsteadycomponentof thewall shear stress, newweighting functions are presented as sums of exponential components. Key words: transient flow, weighting function, friction losses 1. Introduction Most paper dealing with numerical calculations and simulations of the unste- ady liquid pipe flow are based on the assumption that the formula for quasi- steady friction can be applied. However it is only correct in the case of slow changes in the fluid velocity field in the pipe cross-section, i.e. either for small accelerations or for low frequencies It fails in the case of simulation of a fast- changing flow, e.g. in the case of water hammer simulation, because received results significantly differ from the results of experimental studies. Earlier models of unsteady friction losses were based on instantaneous va- lues of velocity and acceleration (Daily et al., 1956; Cartens and Roller, 1959; 488 K. Urbanowicz, Z. Zarzycki Safwat and Polder, 1973; Shuy and Apelt, 1987; Brunone et al., 1991; Vit- kovsky et al., 2000; Bughazem and Anderson, 2000; Bergant et al., 2001). Currently, models based on the history of the flow are commonly used. The forerunner in this group of models was Zielke (1968), who presented the in- stantaneouswall shear stress τu in formof an integral convolution of themean local acceleration of the liquid and a weighting function w(t) τu = 2µ R t∫ 0 ∂v ∂t (u)w(t−u) du (1.1) where: µ – dynamic viscosity, R – inner radius of pipe, v – instantaneous mean flow velocity, t – time, u – time used in the convolution integral, w(t) – weighting function. This dependence is correct for any changes in the average velocity of flow in the pipe cross-section and relates to laminar flow.Thismodel requires time- consuming numerical calculations due to continuous referring to the history of the flow velocity. Therefore, it has been simplified by the introduction of the so-called effective weighting function. Trikha (1975) was first to present the effective expression of the Zielke weighting model, but this relationship has a limited range of applications. Then, new computing models, based on the approximation of the Zielke weighting function were created byKagawa et al. (1983), Suzuki et al. (1991) and Schohl (1993). In the case of transient turbulent flow, inmost of the scientific papers, the models of unsteady friction losses are based on the two-dimensional Reynolds equation and the Boussinesqu hypothesis. In addition, the experimental data for the turbulent viscosity coefficient distribution in the pipe cross section (in its various layers) are used. In the literature one can find so-called two region models (Vardy et al., 1993; Vardy andBrown, 1995, 1996; Popov, 1982; Brown et al., 1969), three-regionmodels (Brown, 1969) and four-regionmodels (Ohmi et al., 1985; Zarzycki, 1994, 2000). Similarly, as it was in the case of laminar flow, the relation describing the instantaneous wall shear stress in form of Eq. (1.1) canalsobeused forunsteadyturbulentflow,but thenweighting functions w(t) are determined on the basis of the abovementionedmulti-regionmodels. In this case, expressions which describe the weighting functions effectively are represented by the followingmodels: Vardy andBrown (2003, 2004), Zarzycki and Kudźma (2004). In the turbulent flow, weighting functions depend not only on dimensionless time (as in the case of laminar flow), but also on the Reynolds number and relative roughness height. Domains of dimensionless time and Reynolds numbers of presented weighting functions do not always correspondtopractical applications.Accordingly, thisworkextended the scope New efficient approximation of weighting functions... 489 of practical applications of new weighting functions. In the case of laminar flow, the most strict model (Zielke, 1968) was used as a base for the process of approximation (a new function is valid in the range 10−9 ¬ t̂¬∞). In the case of turbulent flow, the commonly used weighting functions: by Zarzycki (2000) andbyVardyandBrown (2003)weremodified to the range ofReynolds numbers 2320¬Re¬ 108 and dimensionless time 10−9 ¬ t̂¬∞. Using themethod of characteristics (MOC) for solution of the equations of motion and continuity, there are a few examples of the application of the new weighting function to thewaterhammer phenomenon.The results of numerical calculations are compared with the results of experimental studies using the experimental data by Holmboe and Rouleau (1967) and by Adamkowski and Lewandowski (2004). The comparison confirmed a good agreement. 2. Governing equations The unsteady flow, accompanying the water hammer effect, may be described by a set of the following partial differential equations (Ohmi et al., 1985; Zarzycki, 1994, 2000): — equation of continuity ∂p ∂t + c2ρ ∂v ∂x =0 (2.1) — equation of motion ∂v ∂t + 1 ρ ∂p ∂x +g sinγ+ 2τ ρR =0 (2.2) where: v = v(x,t) – average velocity of the liquid in the pipe cross-section, p= p(x,t) – average pressure in the pipe cross-section, τ – wall shear stress, ρ – density of the fluid, g – gravitational acceleration, γ – inclination angle of the hydraulic line, c – velocity of the pressure wave propagation, t – time, x – axial location along the pipe. Amongmethods, which enable one to resolve the system of the above equ- ations, particular attention should be paid to the method of characteristics (MOC), which perfectly interprets the essence of natural phenomena of tran- sient flow, and at the same time is characterized by fast convergence, taking easily to take into account various boundary conditions and high accuracy of calculation results. With its help, one can easily solve a system of partial differential equations of the quasi-linear hyperbolic type, Eqs. (2.1) and (2.2). The solution is to find the equivalent to four ordinary differential equations, 490 K. Urbanowicz, Z. Zarzycki which are then solved using themethod of finite differences.Approximation to first-order differential schemes gives satisfactory results. Most computer pro- grams used in computational simulations of transients in pipe systems have implemented simple computational algorithms, which adopt a quasi-steady hydraulic resistance. Novelty of thiswork is to include decription of the unste- ady hydraulic resistance of the transient flow in a pipeline, for which it is easy to conduct an efficient simulation for both laminar and turbulent flow. 3. Equations of unsteady hydraulic resistance Zielke analyzed the relationship presented in the work of Brown describing the impedance of a hydraulic line versus frequency. Application of the inverse Laplace transform brought him to the following expression (Zielke, 1968) τ(t)= ρv|v| 8 λ ︸ ︷︷ ︸ τq + 2µ R t∫ 0 w(t−u) ∂v ∂t (u) du ︸ ︷︷ ︸ τn (3.1) where: w(t−u) – weighting function; t – actual time in numerical simulation; u – integral variable having dimension of time. The first component τq of equation (3.1) represents the quasi-steady amo- untofwall shear stress.Theother one, τu, describes the impactof theunsteady effect of flow on the wall shear stress. 3.1. The laminar flow Zielke (1968) presented the weighting functions for laminar flow in the following form w(t̂)=            6 ∑ i=1 mit̂ (i−2)/2 for t̂¬ 0.02 5 ∑ i=1 e−nit̂ for t̂ > 0.02 (3.2) where: t̂ = νt/R2 – diensionless time, coefficients mi and ni successively take the following values: mi = 0.282095, −1.25, 1.057855, 0.9375, 0.396696, −0.351563; ni =26.3744, 70.8493, 135.0198, 218.9216, 322.5544. New efficient approximation of weighting functions... 491 Numerical calculation of the time-dependent component of the wall shear stress τu (second component in Eq. (3.1)) can be made using the first-order differential approximation (Zarzycki and Kudźma, 2004; Zielke, 1968) τn = 2µ R k−1 ∑ j=1 (vi,j+1−vi,j)w ( (k− j)∆t̂− ∆t̂ 2 ) = 2µ R k−1 ∑ j=1 (vi,k−j+1−vi,k−j)w ( j∆t̂− ∆t̂ 2 ) (3.3) where: j –number of computational time step changing from1 to k for k­ 2, ∆t̂= ν∆t/R2. Determination of the wall shear stress using above formula (3.3) is very time consuming. Trikha (1975)was first to develop an effective method of so- lving the integral convolution (later Kagawa et al. (1983), Suzuki et al. (1991) and Schohl (1993) had improved that method). In this paper, the effective solution of Kagawa et al. (1983) is used τu = 2µ R k∑ i=1 ( yi(t)e −ni∆t̂+mie −ni∆t̂/2[v(t+∆t)−v(t)] ︸ ︷︷ ︸ yi(t+∆t) ) (3.4) Thismethod, however, requires that theweighting function has to bewrit- ten as a finite sum of exponential expressions w(t̂)= k ∑ i=1 mie −nit̂ (3.5) Number of exponential terms that make up the final form of the effective weighting function, affects the range of its applicability as well as its degree of fit to the original function (according to Zielke). Over the past 35 years, many authors have presented their effective weighting functions for the case of laminar flow (Kagawa et al., 1983; Schohl, 1993; Trikha, 1975; Vardy and Brown, 2004; Vitkovský et al., 2004). For the ranges of their applicability and their estimated coefficients – see Tables A1 and A2 in Appendix A. The course of the laminar weighting function is shown in Fig.1. 3.2. Turbulent flow BothVardyandBrown (1995, 1996, 2003, 2004), andZarzycki (1994, 2000) suggest that the relation for theunsteadypartof thewall shear stresspresented 492 K. Urbanowicz, Z. Zarzycki Fig. 1. Comparison of laminar function runs; (a) log-linear scale, (b) log-log scale byZielke, Eq. (3.1), is also true for the turbulent flow.However, theweighting function in this case has different shapedue to its dependence on theReynolds number. The original turbulent weighting functions are: — Zarzycki’s weighting function (Zarzycki, 2000) w(t̂,Re)=C 1√ t̂ Ren (3.6) where: C =0.299635, n=−0.005535. —Vardy’s and Brown’s weighting function (Vardy and Brown, 2003, 2004) w(t̂,Re)= A∗e−B ∗t̂ √ t̂ (3.7) where: A∗ = √ 1/(4π) and B∗ = Reκ/12.86, κ = log(15.29/Re0.0567) for smooth pipes and A∗ =0.0103 √ Re(ε/D)0.39 and B∗ =0.352Re(ε/D)0.41 for flows in rough pipes. The ratio ε/D is called the relative roughness (in one of the recent works by Vardy and Brown, there is a proposition to set the parameters A∗ and B∗ from more complicated equations, shown in Vardy and Brown (2007)). But they were not suitable for efficient simulations using expression (3.3) for wall shear stress calculations. From computational point of view it is very important to be able to conduct effective (in terms of high accuracy, fast working numerical scheme and less memory usage) simulations of transients both in laminar and turbulent flow. Therefore, in literature, one can find, so-called, effective expressions which are approximations of original models of turbulent weighting functions. For ranges of their applicability and their estimated coefficients – see Tables A3 and A4 in Appendix A. New efficient approximation of weighting functions... 493 4. Developing effective weighting functions 4.1. Weighting function for laminar flow In view of the need to extend the applicability of the effective function of weighting which was noted by, among others, Vardy and Brown (2004, 2007), a new model which is the sum of exponential expressions will be presented further. The new model is consistent with the original function by Zielke in the following range of applicability 10−9 ¬ t̂ < ∞. The final form of the new function consists of 26 exponential expressions. Since the value of Zielke we- ighting function for the dimensionless time t̂ > 0.02must be determined from the following formula (3.2) w(t̂)= 5∑ i=1 e−nit̂ (4.1) where: n1 = 26.3744, n2 = 70.8493, n3 = 135.0198, n4 = 218.9216, n5 =322.5544. In this work these first five exponential expressions were kept unchanged (in order to receive perfect fit of the newweighting function in this interval of dimensionless time). For the mapping of the interval 10−9 ¬ t̂ < 0.02 it was decided to add extra similar components. It was assumed that a very good accuracy would be received by describing each dimensionless time interval 10n−1 ¬ t̂ < 10n with three exponential expressions (except for the range 10−3 ¬ t̂ < 0.02, which also describes the three formulas – so the worst match of the newweighting functionwas expected). In each of those intervals, matching was carried out to uniformly distributed 1000 points (in log scale), which were the results obtained by using the original function of weighting according to Zielke. The effective weighting function coefficients were determined by using the LSQNONLIN function, which is a module of MATLAB. In this function the Levenberg-Marquardt algorithm is implemented, consideredas one of themost effective among the minimization algorithms. It combines the linear approxi- mation away fromtheminimumand square approximation near theminimum, so that it is specialized to problems of the method of least squares. Following the procedure outlined above, all coefficients of the new effective weighting function for laminar flow were determined wapr(t̂)= 26∑ i=1 mie −nit̂ (4.2) 494 K. Urbanowicz, Z. Zarzycki where: m1 = 1, m2 = 1, m3 = 1, m4 = 1, m5 = 1, m6 = 2.141, m7 = 4.544, m8 = 7.566, m9 = 11.299, m10 = 16.531, m11 = 24.794, m12 = 36.229, m13 = 52.576, m14 =78.150, m15 = 113.873, m16 = 165.353, m17 = 247.915, m18 = 369.561, m19 = 546.456, m20 = 818.871, m21 = 1209.771, m22 = = 1770.756, m23 = 2651.257, m24 = 3968.686, m25 = 5789.566, m26 = = 8949.468, n1 = 26.3744, n2 = 70.8493, n3 = 135.0198, n4 = 218.9216, n5 = 322.5544, n6 = 499.148, n7 = 1072.543, n8 = 2663.013, n9 = 6566.001, n10 =15410.459, n11 =35414.779, n12 =80188.189, n13 =177078.960, n14 = = 388697.936, n15 = 850530.325, n16 = 1835847.582, n17 = 3977177.832, n18 = 8721494.927, n19 = 19120835.527, n20 = 42098544.558, n21 = = 92940512.285, n22 = 203458923.000, n23 = 445270063.893, n24 = =985067938.878, n25 =2166385706.058, n26 =4766167206.672. 4.2. Weightig functions for turbulent flow Theneweffectiveweighting function for turbulentflowwill bebased on the Zarzycki original model ((Zarzycki, 1994, 2000; Zarzycki andKudźma, 2004). The form of a new function is selected so that it is easy to scale it. Scaling will depend in detail on multiplying the estimated coefficients (assuming a constant value of the Reynolds number Re during the transient flow) by a function of the Reynolds number (which will be fixed prior to simulation) (Vitkovský et al., 2004). It will ensure a good fit for the standard functions (the originalweighting function) of the effective function throughout the range of applicability. Theweighting function byZarzycki is a function of dimensionless time and the Reynolds number Re (3.6). Following the way presented by Vitkovský et al. (2004), it can be written wapr(t̂,Re)= CRen√ t̂ ≈ k∑ i=1 mie −nit̂ (4.3) Then, dividing the above dependence by CRen, the weighting function has form w∗apr = 1√ t̂ ≈ k ∑ i=1 mi CRen e−nit̂ w∗apr ≈ Re−n C k ∑ i=1 mie −nit̂ w∗apr ≈ Re0.005535 0.299635 k ∑ i=1 mie −nit̂ w∗apr ≈ k ∑ i=1 m∗ie −nit̂ (4.4) where: m∗i – universal coefficients for the effective weighting function of tur- bulent flow. New efficient approximation of weighting functions... 495 This procedure shows how the universal coefficients of the effective weigh- ting function for the turbulent flow are determined. The current form of the effective weighting function is determined with their help as a result of the reverse scaling. In this paper, the universal coefficients m∗i were determined by adjusting the effective weighting function to the original by Zarzycki for the constant Reynolds number Re=10000 m∗i = mi(Re=10000) C Re0.005535 = mi(Re=10000) 0.299635 100000.005535 (4.5) Themaking use of universal coefficientsmust ensure that the number of Re= 10000 are derived values, which were initially estimated, i.e. mi(Re=10000). Therefore, thedesired effectiveweighting function for turbulentflowmusthave the following form wapr(t̂,Re)≈ C Re−n k ∑ i=1 m∗ie −nit̂ = 0.299635 Re0.005535 k ∑ i=1 m∗ie −nit̂ (4.6) Assume that the scope of the new effective weighting function for turbulent flowmustbeas follows: 10n−1 ¬ t̂ < 10n tobe fully suited for theuse inpresent technical issues. Then it was found that for each interval 10n−1 ¬ t̂ < 10n (for n = −9 to n = 3) for a good estimate of the coefficients characterizing this function, just two power terms (mie −nit̂) are efficient. Hence, for the entire range of applicability in the field of dimensionless time 10−9 ¬ t̂ < 103 24 exponential expressions were finally received. The sum of these expressions is theneweffectiveweighting function for turbulentflow. Inorder to determine such a large number of factors (in total 48 coefficients) it was assumed that the successive coefficients of determination in stepsmust be followed. For each step of estimation, therewere twonew exponential terms added (4 consecutive coefficients) to the searched function (power series). The searching procedure was started for the dimensionless time period 102 ¬ t̂ < 103. In order to accurately estimate the coefficients in each interval of analy- sis 10n−1 ¬ t̂ < 10n, the estimation was based on equally distributed 1000 points (similar like when setting the new effective function for laminar flow), which were the results obtained from using the original weighting function by Zarzycki. Following the above procedure, all the estimated coefficients of the new effective weighting function for turbulent flow were determined wapr(t̂,Re)= C Re−n 24∑ i=1 m∗ie −nit̂ = 0.299635 Re0.005535 24∑ i=1 m∗ie −nit̂ (4.7) 496 K. Urbanowicz, Z. Zarzycki where: m∗1 = 0.06054, m ∗ 2 = 0.09698, m ∗ 3 = 0.17971, m ∗ 4 = 0.31240, m∗5 = 0.56562, m ∗ 6 = 0.98348, m ∗ 7 = 1.77243, m ∗ 8 = 3.08626, m ∗ 9 = 5.57348, m∗10 = 9.7254, m ∗ 11 = 17.591, m ∗ 12 = 30.723, m ∗ 13 = 55.603, m ∗ 14 = = 97.138, m∗15 = 175.825, m ∗ 16 = 307.176, m ∗ 17 = 551.342, m ∗ 18 = 954.362, m∗19 = 1727.71, m ∗ 20 = 3171.2, m ∗ 21 = 5899.4, m ∗ 22 = 11013, m ∗ 23 = 19923, m∗24 = 37929, n1 = 0.000671, n2 = 0.00838, n3 = 0.04504, n4 = 0.1790, n5 =0.6457, n6 =2.159, n7 =7.088, n8 =22.563, n9 =72.215, n10 =227.12, n11 = 723.19, n12 = 2270.23, n13 = 7226.1, n14 = 22686.2, n15 = 72226.7, n16 =226796, n17 =720015, n18 =2234661, n19 =7050737, n20 =22553627, n21 =74840660, n22 =253286747, n23 =856109205, n24 =2893640000. The range of applicability of the new effective weighting functions, presen- ted in the last two subsections, (laminar flow 10−9 ¬ t̂ < ∞, turbulent flow 10−9 ¬ t̂ < 103 and 2300 ¬ Re ¬ 108) virtually guarantees a very accurate assessment of hydraulic resistance which is useful in simulation of the unste- ady cavitating flow in pipelines. However, one can move beyond this range, when it is necessary to substantially thicken the grid of characteristics (trac- king changes of flow parameters in a number of cross-sections in the analyzed pressure line). Then, in the numerical analysis, the changes in a very short time t̂ < 10−9 can be taken into account (see example in Appendix B). Therefore, there are no decisive experimental results supporting one of the two (considered in this paper) original weighting functions for turbulent flow (byVardyandBrownandbyZarzycki).There is alreadyan effectiveweighting function presented by Vitkovský et al. (2004) (for details see Appendix A) based on the model by Vardy and Brown (3.7), which is characterized by a goodfit,buta small rangeof applicability.There is thereforenoneed topresent a completely new function. It is simply enough to extend the applicability of the existing one by finding new extra exponential terms. Tominimize the impact of the last termof theVitkovský et al. (2004) func- tion on quality ofmappingof the expanded function to the original function, it was also replaced by a new one (new estimated values of the coefficients: m∗10 and n∗10). Moreover, in the estimation process it was assumed that the lower range of applicability of the new function must be t̂=10−9. For this purpo- se, it was necessary to find seven new exponential expressions. The procedure used to achieve this goal was almost identical to that used in the estimating of the weighting function for the turbulent flow based on the original weigh- ting function by Zarzycki. Thus, below is only presented the final form of this function wapr(t̂,Re)= 16∑ i=1 A∗m∗ie −(n∗ i +B∗)t̂ (4.8) New efficient approximation of weighting functions... 497 where: A∗ and B∗ are parameters known from Vardy’s and Brown’s equ- ation (3.7) and m∗1 = 5.03362, m ∗ 2 = 6.4876, m ∗ 3 = 10.7735, m ∗ 4 = 19.904, m∗5 = 37.4754, m ∗ 6 = 70.7117, m ∗ 7 = 133.460, m ∗ 8 = 251.933, m ∗ 9 = 476.597, m∗10 =902.22, m ∗ 11 =1602.04, m ∗ 12 =2894.84, m ∗ 13 =5085.55, m ∗ 14 =9190.11, m∗15 = 16118.6, m ∗ 16 = 29117.3, n ∗ 1 = 4.78793, n ∗ 2 = 51.0897, n ∗ 3 = 210.868, n∗4 = 765.03, n ∗ 5 = 2731.01, n ∗ 6 = 9731.44, n ∗ 7 = 34668.5, n ∗ 8 = 123511, n∗9 =440374,n ∗ 10 =1578229,n ∗ 11 =5481659,n ∗ 12 =18255921,n ∗ 13 =59753474, n∗14 =192067361, n ∗ 15 =616415963, n ∗ 16 =1945566788. The scope of applicability of all new weighting functions presented above (4.2), (4.7) and (4.8) can be any further extended by adding the new expo- nential expressions. 4.3. Comparison of the newweighting functionswith their original coun- terparts and the best known effective features from the literature The following figures are presented in order to show a comparison of the new, proposed in the previous section, weighting functions with their ineffi- cient counterparts and the most precise effective functions known from the literature. In addition, in order to demonstrate the degree of matching obta- ined by the weighting functions with their prototypes, results of quantitative analysis are graphically presented. Fig. 2. Comparison of weighting functions for laminar flow: (a) log-log graph, (b) relative percentage error As aparameter determining quantitatively the degree ofmatching the new efficient weighting function with its original counterpart, the relative percen- tage error determined from the following relationship was incorporated Rerror = wapr−w w 100% (4.9) 498 K. Urbanowicz, Z. Zarzycki Fig. 3. Comparison of weighting functions (original Zarzycki’s and its known effective counterparts) for turbulent flow (Re=104): (a) log-log graph, (b) relative percentage error Fig. 4. Comparison of weighting functions (original Vardy and Brown and its known effective counterparts) for turbulent flow (Re=104): (a) log-log graph, (b) relative percentage error Comparisons of theweighting functions for turbulent flow for awider range of theReynolds numbers are not shown in thiswork. Indeed, the use of scaling procedurespresentedbyVitkovský et al. (2004)means that for otherReynolds numbers (from the range of applicability) the errors remain the same. 5. Numerical results In order to compare the accuracy of unsteady (with the use of original and universal weighting functions) and quasi-steady models of friction in relation New efficient approximation of weighting functions... 499 to experimental data, simulations of a simple waterhammer case (tank – long liquid line and cut-off valve) were conducted. Thecomputed results (themethodof characteristics basedona rectangular gridwasused;hydrayulic linewasdividedon N =16elements)were compared with the experimental data reported by Holmboe and Rouleau (1967), and Adamkowski and Lewandowski (2004). 5.1. Holmboe and Rouleau experiment Holmboe and Rouleau (1967) ran their tests on a copper tube with radius 0.0127m and length 36.09m connected upstream to a tank which was main- tained at a constant pressure by the compressed air. The liquid used in the experiment was an oil having viscosity 39.6710−6m2/s. The measured speed of the pressure wave was 1324.36m/s, and the initial flow velocity 0.128m/s (Re = 82). The downstream valve was rapidly closed in the pipe line during flow. Pressure fluctuation was measured at the endpoint of the line (near the valve). From the above parameters, it followed that it was a case of laminar flow. It was determined in numerical calculations, in which themodels of Ziel- ke (3.2) and the new effective laminar model (4.2) were used. In addition, the calculation with the usage of quasi-steady model was conducted. The results of simulations and experimental data are shown in Fig.5. Fig. 5. Fluctuations of pressure at the endpoint of the line (Re=82) – Holmboe and Rouleau experiment It is clearly seen that the calculation using unsteady friction models is much closer to the experimental data. Therefore, in further calculations, the unsteady friction models were used instead of the quasi-steady formula. Ad- ditionally, it is significant that the differences between the simulation results 500 K. Urbanowicz, Z. Zarzycki according to the classical Zielke model (3.2) and the new proposed function (4.2) (efficient way of unsteady wall shear stress computation) are hardly di- stinguishable. 5.2. Adamkowski and Lewandowski experiment Adamkowski and Lewandowski (2004) conducted experiments at a test rig specially designed in order to investigate the unsteady pipe flows. Its main component was a 98.11m long copper pipe with internal diameter of 0.016m and wall thickness of 0.001m. The pipe was rigidly mounted to the ground using bearings in order to minimize its vibrations. A quick-closing spring dri- ven ball valve was installed at one end of the pipe. The initial conditions were defined by the high pressure reservoir pressure head and the initial flow velo- city in the pipeline. During the tests, temperature of water was 22.6◦C, and the kinematic viscosity coefficient for these conditions was 9.493 ·10−7m2/s. Three runs were selected for the purpose of verifying new unsteady friction models. The initial flow velocities were 0.066m/s (Re = 1100), 0.631m/s (Re = 10600) and 0.927m/s (Re = 15600). The results of simulations and experimental data are shown in Figs. 6, 7 and 8. Fig. 6. Fluctuations of pressure at the endpoint of the line (Re=1100) – Adamkowski and Lewandowski experiment Fromtheabove graphs, it is clear that theuse of thenewefficientweighting functions (4.2), (4.7) and (4.8) ensures compatibility of results with those obtained from the original models by Zielke (3.2), Vardy andBrown (3.7) and Zarzycki (3.6). New efficient approximation of weighting functions... 501 Fig. 7. Fluctuations of pressure at the endpoint of the line (Re=10600) – Adamkowski and Lewandowski experiment Fig. 8. Fluctuations of pressure at the endpoint of the line (Re=15600) – Adamkowski and Lewandowski experiment 6. Summary The main drawback of the original models (according to Zielke, Vardy and Brown, and Zarzycki) describing the unsteady hydraulic resistance is their inefficiency! In each successive time step they require more resources of the computer memory (this is due to augmentations in quantity of elements that make up the sum creating the solution to the integral convolutions (3.3), as well as due to enlargement of the matrix which stores information about past velocity changes v(t)). Often, after the time when the entire memory of the computer is used, it comes to forced interruption of performance simulation. Hence the original weighting functions thatmake the obtained results coincide with the results of experimental studies (Vardy and Brown, 1995; Zarzycki, 1994), even today, with tremendous development of computer technology, fit for simulations of very short transients only. 502 K. Urbanowicz, Z. Zarzycki The known efficient solution of the integral convolution (3.4) for both la- minar and turbulent flow (Kagawa et al., 1983; Schohl, 1993; Zarzycki and Kudźma, 2004) allowed, in the case of using the weighting function as a finite sum of exponential terms ∑k i=1mie −nit̂, a significant reduction in computa- tion time and reduced demand for operational memory. Importantly, the time of calculation with the increasing number of time steps, increased almost line- arly, as shown in Fig.9 (calculations were carried out on a standard desktop computer –Fujitsu Siemens-Intel Core 2QuadCPUQ6600, 2.4GHz, 2048MB RAM). The impact of number of weighting function terms on time of calcu- lation was also investigated. The two cases are presented in Fig.9: 10 terms (Eq. (3.4) and Vitkovský et al. laminar function (2004)) and 26 terms (Eqs. (3.4) and (4.2)). The increase of time consumption for efficient cases is very small in relation to the original model (Eqs. (3.2) and (3.3)). Fig. 9. Time of numerical calculations depending on the number of time steps Numerous simulation tests carried out by the authors revealed that the nu- merical results usingmodelswith small ranges of applicability (e.g., byTrikha) often deviate from the experimental results (as a result of going beyond the scope of applicability of the used weighting function). Professional software, which may in future be based on weighting functions presented in this work, must prevent crossing outside the range of applicability of these functions – by informing the user about the need to change numerical parameters (e.g. the number of simultaneously observed pipe cross sections). The presented newweighting functions characterized by the extended ran- ges of applicability in the domain of dimensionless time are suitable for the efficientmodeling of transients. The studies show that the results of numerical simulations, in which a newweighting functions were used, overlap with those that use the original models. New efficient approximation of weighting functions... 503 In the next stage of research, the usefulness of the newweighting functions presented in this paper in simulation of transients with cavitation should be explored. Appendix A TableA1.Ranges of applicability of the effective laminarweighting functions Model by Trikha Schohl Kagawa Vitkovský Vardy and (1975) (1993) et al. (1983) et al. (2004) Brown (2004) range of 7.41 ·10−5 ¬ 1.26 ·10−5 ¬ 6.31 ·10−6 ¬ t̂¬∞ 10 −8 ¬ applicab. t̂¬ 10 ¬ t̂¬ 1 ¬ t̂¬∞ Table A2.Estimated coefficients of the effective laminar weighting functions Trikha Schohl Kagawa Vitkovský Vardy and i (1975) (1993) et al. (1983) et al. (2004) Brown (2004) mi ni mi ni mi ni mi ni mi ni 1 1 26.4 1.051 26.65 1 26.3744 1 26.3744 1 A 2 8.1 200 2.358 100 1.16725 72.8033 1.09301 72.044 2.1830 102 3 40 8000 9.021 669.6 2.20064 187.424 1.82206 166.931 2.714 102.5 4 29.47 6497 3.92861 536.626 3.34085 435.932 7.5455 103 5 79.55 57990 6.78788 1570.60 5.89377 1229.74 39.0066 104 6 11.6761 4618.13 10.2835 3584.84 106.8075 105 7 20.0612 13601.1 17.9006 10621.7 359.0847 106 8 34.4541 40082.5 31.1516 31757 1107.9295 107 9 59.1642 118153 54.4168 95563.7 3540.683 108 10 101.59 348316 99.4360 293268 A=26.3744 Table A3. Ranges of applicability of the effective turbulent weighting func- tions Model by Zarzycki and Kudźma Vitkovský Vardy and Kudźma (2004) (2005) et al. (2004) Brown (2004) range of 10−7 ¬ t̂¬ 10−3 10−5 ¬ t̂¬ 10 10−6 ¬ t̂¬ 10−1 10−9 ¬ t̂¬ 10−1 applicab. 504 K. Urbanowicz, Z. Zarzycki TableA4.Estimated coefficients of the effective turbulentweighting functions Based on Zarzycki’s Based on Vardy and Brown original expression (3.6) original expression (3.7) Zarzycki and Kudźma Vitkovský Vardy and i Kudźma (1975)1 (2005)2 et al. (2004) Brown (2004) wapr(t̂,Re)= k ∑ i=1 C∗m∗ie −nit̂ wapr(t̂,Re)= k ∑ i=1 A∗m∗ie −(n∗ i +B∗)t̂ m∗i ni m ∗ i ni m ∗ i n ∗ i m ∗ i n ∗ i 1 17.10735 477.887 0.2137 0.09834 5.03362 4.78793 9.06 10 2 58.51351 17790.69 1.568 8.44 6.4876 51.0897 −4.05 101.5 3 152.3936 207569.7 2.799 88.02 10.7735 210.868 12 102 4 328.2 1464649 5.527 480.5 19.904 765.03 8.05 102.5 5 414.8145 6316096 10.76 2162 37.4754 2731.01 22.7 103 6 640.2165 15512625 18.99 8425 70.7117 9731.44 35.2 103.5 7 33.26 29250 133.46 34668.5 65.9 104 8 60.73 96940 251.933 123511 115 104.5 9 476.597 440374 206 105 10 932.86 1590300 365 105.5 11 651 106 12 1150 106.5 13 2060 107 14 3630 107.5 15 6640 108 16 10700 108.5 17 26200 109 1 C∗ =(−1.5125Re0.003264+2.55888); 2 C∗ =(−13.27813Re0.000391+14.27658) A∗ and B∗ are parameters known fromVardy’s and Brown’s equation (3.7) Appendix B In the numerical analysis of transients, one can sometimes gobeyond the scope of applicability of the weighting functions (thus committing a serious error in determining the hydraulic resistance). This has not been previously discussed in the literature on the research of such states. Below, there are two theoretical cases explaining the importance of exten- ding the scope of applicability of the weighting functions (by adding new exponential expressions) to the lower range of dimensionless time t̂=10−9. Case I L = 25m, ν = 10−6m2/s, N = 100, c = 1225m/s, R = 0.13m – we will have New efficient approximation of weighting functions... 505 ∆t̂= Lν NcR2 = 25 ·10−6 100 ·1225 ·0.132 =1.2076 ·10−8 In the case of a rectangular grid for the first computational step (j=1) using classical formula (3.3) τu(i,k) = 2µ R k−1 ∑ j=1 (v(i,k−j+1)−v(i,k−j))w ( j∆t̂− ∆t̂ 2 ) as well as using effective formula (Vitkovský et al., 2004) τu(i,k) = 2µ R z ∑ j=1 [ yj(i,k−1)e −nj∆t̂+mje −nj ∆t̂ 2 (v(i,k)−v(i,k−1)] ] ︸ ︷︷ ︸ yj(i,k) the value of the weighting function is needed for the argument equal to ∆t̂/2. Itmeans that for this case, the functionhave to be in the range of applicability beginning at ∆t̂/2= 6.038 ·10−9. The values of different weighting functions for this argument are shown in Table B1. TableB1.Values of the laminarweighting function for thedimensionless time ∆t̂/2=6.038 ·10−9 Function Zielke New function Vardy and Vitkovský Kagawa (3.2) (4.2) Brown (2004) et al. (2004) et al. (1983) Value of weighting 3629.103 3629.157 3494.923 226.123 241.764 function [–] Relative 0 0.0015 −3.6973 −93.7692 −93.3382 error [%] Case II If: L=1000m, ν =10−5m2/s,N =16, c=1280m/s,R=0.008m, then ∆t̂= Lν NcR2 = 1000 ·10−5 16 ·1280 ·0.0082 =7.6 ·10−3 and ∆t̂/2=3.8 ·10−3. In this case, there is no need to usemany exponential expressions – just 5 is enough (for both turbulent and laminar flow) in order to properly simulate the hydraulic resistance. 506 K. Urbanowicz, Z. Zarzycki References 1. Adamkowski A., Lewandowski M., 2004, Experimental examination of unsteady friction models for transient pipe flow simulation, 9th Int. conf. on Pressure Surges BHR 2004, Proc., II, 421-437 2. Bergant A., Simpson A., Vtkovský J., 2001, Developments in unste- ady pipe flow friction modelling, Journal of Hydraulic Research, IAHR, 39, 3, 249-257 3. BrownF.T., 1969,A quasimethod of characteristicswith application to fluid lines with frequency dependent wall shear and heat transfer,Trans. ASME, J. Basic Engineering, 91, 2, 217-227 4. Brown F.T., Margolis D.L., ShahR.P., 1969, Small-amplitude frequency behaviour of fluid lines with turbulent flow, Trans. ASME, J. Basic Engine- ering, 678-693 5. Brunone B., Golia U.M., Greco M., 1991, Some remarks on the momen- tum equations for fast transients,Proceedings of the International Meeting on Hydraulic Transients with Column Separation, 9thRoundTable, IAHR,Valen- cia, Spain, 201-209 6. Bughazem M.B., Anderson A., 2000, Investigation of an unsteady friction model for waterhammer and column separation, Proceedings of the 8th Inter- national Conference on Pressure Surges, BHR Group Conf. Series Publ., The Hague, The Netherlands, 483-498 7. Carstens M.R., Roller J.E., 1959, Boundary-shear stress in unsteady tur- bulent pipe flow,ASCE Journal of the Hydraulics Division, 85, HY2, 67-81 8. Dailey J.W., HankeyW.L., OliveR.W., Jordaan J.M., 1956,Resistan- ce coefficient for accelerated and decelerated flows through smooth tubes and orifices, J. Basic Engineering, 78, July, 1071-1077 9. Holmboe E.L., Rouleau W.T., 1967, The effect of viscous shear on tran- sients in liquid lines, Journ. of Basic Eng., Trans. ASME, s. D,89, 11, 174-180 10. Kagawa T., Lee I., Kitagawa A., Takenaka T., 1983, High speed and accurate computing method of frequency-dependent friction in laminar pipe flow for characteristicsmethod, Trans. Jpn. Soc. Mech. Eng., Ser. A, 49, 447, 2638-2644 [in Japanese] 11. Kudźma S., 2005, Modeling and simulation dynamical runs in closed condu- its of hydraulics systems using unsteady friction model, PhD Thesis, Szczecin University of Technology [in Polish] 12. Ohmi M., Kyomen S., Usui T., 1985,Numerical analysis of transient flow in a liquid line,Bulletin of JSME, 28, 239, 799-806 New efficient approximation of weighting functions... 507 13. Popov D.N., 1982,Unsteady-State Hydromechanical Processes, Moscow,Ma- shinostroenie [in Russia] 14. Safwat H.H., Polder J. van Den, 1973, Experimental and analytic data correlation study of water column separation, Journal of Fluids Engineering, March, 91-97 15. Schohl G.A., 1993, Improved approximate method for simulating frequency – dependent friction in transient laminar flow, Journal of Fluids Engineering, Trans. ASME, 115, September, 420-424 16. ShuyE.B.,ApletC.J., 1987,Experimental studies ofunsteadyshear stress in turbulent flow in smooth roundpipes,Conf. onHydraulics inCivil Engineering, Melbourne, Australia, 137-141 17. Suzuki K., Taketomi T., Sato S., 1991, Improving Zielke’s method of si- mulating frequency – dependent friction in laminar liquid pipe flow, Journal of Fluids Engineering, Trans. ASME, 113, December, 569-573 18. Trikha A.K., 1975, An efficient method for simulating frequency-dependent friction in transient liquid flow, Journal of Fluids Engineering, Trans. ASME, March, 97-105 19. Vardy A.E., Brown J.M.B., 1995, Transient, turbulent, smooth pipe fric- tion, J. Hydraul. Res., 33, 435-456 20. Vardy A.E., Brown J.M.B., 1996, On turbulent unsteady, smooth – pipe friction,Proc. of the 7th International Conf. on Pressure Surges – BHRGroup, Harrogate, United Kingdom, 289-311 21. Vardy A.E., Brown J.M.B., 2003, Transient turbulent friction in smooth pipe flows, Journal of Sound and Vibration, 259, 5, 1011-1036 22. VardyA.E.,BrownJ.M.B., 2004a,Efficient approximationof unsteady fric- tion weighting functions, Journal of Hydraulic Engineering, ASCE, 130, 11, 1097-1107 23. Vardy A.E., Brown J.M.B., 2004b, Transient turbulent friction in fully ro- ugh pipe flows, Journal of Sound and Vibration, 270, 233-257 24. Vardy A., Brown J., 2007, Approximation of turbulent wall shear stres- ses in highly transient pipe flows, Journal of Hydraulic Engineering, 133, 11, 1219-1228 25. Vardy A.E., HwangK.L., Brown J.M.B., 1993,Aweighting functionmo- del of transient turbulent pipe friction, J. Hydraul. Res., 31, 4, 533-548 26. Vitkovsky J., Lambert M., Simpson A., Bergant A., 2000, Advances in unsteady friction modelling in transient pipe flow, Proceedings of the 8th International Conference on Pressure Surges, BHR Group Conf. Series Publ., The Hague, The Netherlands, 471-482 508 K. Urbanowicz, Z. Zarzycki 27. Vtkovský J.P., Stephens M.L., Bergant A., Simpson A.R., Lambert M.F., 2004,Efficient and accurate calculation of Zilke andVardy-Brownunste- ady friction inpipe transients,9th International Conference onPressure Surges, Chester, United Kingdom, 405-419 28. ZarzyckiZ., 1994,AHydraulicResistance’s ofUnsteadyLiquidFlow inPipes, Published by Technical University of Szczecin, No 516, Szczecin [in Polish] 29. ZarzyckiZ., 2000,Onweighting function forwall shear stressduringunsteady turbulent flow,Proc. of 8th International Conference on Pressure Sergues, The Hague, The Netherlands, BHRGroup Conference Series, 39, 529-534 30. Zarzycki Z., Kudźma S., 2004, Simulations of transient turbulent flow in liquid lines using time – dependent frictional losses, Proceedings of the 9th International Conference on Pressure Surges, BHR Group, Chester UK, 439-455 31. ZielkeW., 1968, Frequency-dependent friction in transient pipe flow, Journal of ASME, 90, March, 109-115 Nowe efektywne funkcje wagowe umożliwiające symulację niestacjonarnych oporów hydraulicznych podczas przepływów cieczy w przewodach ciśnieniowych Streszczenie Wdużej części prac, dotyczącychobliczania oraz symulacji niestacjonarnegoprze- pływucieczywprzewodachciśnieniowych, oporywyznacza się jedynie zdobrze znanej formuły Darcy-Weisbacha. Niestety w szybkozmiennych przepływach, takich jak np. uderzeniehydrauliczne, założeniequasi-stacjonarnościoporówjestpoważnymbłędem. Wponiższej pracynaprężenie stycznena ścianceprzewoduwyznacza się z sumyquasi- stacjonarnego składnika oraz niestacjonarnego. Niestacjonarny składnik naprężenia stycznegomodelowany jest jako splot lokalnego przyśpieszenia cieczy oraz funkcji wa- gi. Klasyczne funkcje wagi mają skomplikowaną postać, co uniemożliwia efektywne wyznaczanie oporów. W tej pracy przedstawione zostały nowe efektywne postacie funkcji wagowych o rozszerzonych zakresach stosowalności. Manuscript received June 6, 2011; accepted for print September 21, 2011