Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 4, pp. 853-871, Warsaw 2007 SIMULATION OF TRANSIENT FLOWS IN A HYDRAULIC SYSTEM WITH A LONG LIQUID LINE Zbigniew Zarzycki Sylwester Kudźma Szczecin University of Technology, Department of Mechanics and Machine Elements, Poland e-mail: zbigniew.zarzycki@ps.pl; sylwester.kudzma@ps.pl Zygmunt Kudźma Michał Stosiak Wrocław University of Technology, Institute of Machine Design and Operation, Poland e-mail: zygmunt.kudzma@pwr.wroc.pl, michal.stosiak@pwr.wroc.pl Thepaperpresents theproblemofmodellingand simulationof transients phenomena in hydraulic systems with long liquid lines. The unsteady resistance model is used to describe the unsteady liquid pipe flow. The wall shear stress at thepipewall is expressedbymeansof the convolution of acceleration and a weighting function which depends on the (laminar or turbulent) character of the flow. The results of numerical simulation are presented for the waterhammer effect, which is caused by a sudden shift of the hydraulic directional control valve. The following cases of the system supply are considered: the first, with a constant delivery rate of the pump and the second, which additionally considers pulsation of the delivery of the pump. Computer simulations are compared with results of experiments. They are found to be very consistent in the case with the variable rate of the pump delivery taken into accent. Key words: unsteady pipe flow, transients, waterhammer, pulsation of pump 1. Introduction Drive and hydraulic control systems are often subjected to transient states, caused by dynamical excitation forces resulting either from sudden changes of the load of a motor or hydraulic actuator or else from changes of the flow 854 Z. Zarzycki et al. direction of the working liquid or speed caused by the control unit. It is often essential to precisely learn characteristics of dynamical runs in such states. It is very important in the design of automatic control systems or in the for analyzis of the strength of pipes and other elements of hydraulic systems. While analyzing transient states, particular attention should be paid to such cases when a hydraulic systemhas a long liquid line (Garbacik and Szew- czyk, 1995; Wylie and Streeter, 1978; Jelali and Kroll 2003), i.e. either when the length of the pipe is close to that of the pressurewavewhich is propagated in the system or when it is higher than that. Such a pipe is treated then as an element of a system with distributed parameters. Consequently, any changes of the flowpressure and rate are distributed along the pipes axiswith a limited speed in the form of progressive and reflected waves. In literature, the research of transient states concernmostly simple water- hammer cases (Wylie and Streeter, 1978; Ohmi et al., 1985), in which simple boundary conditions are used. It means that at one end of the hydraulic line the pressure is constant (reservoir) and at the other, the velocity of flow equals zero (suddenly closed valve). The present paper is an attempt at simulating transients caused by a sud- den change of settings of the hydraulic control valve, with taking into account at the same time the pulsation delivery rate resulting from kinematics of a positive-displacement pump.The results obtained from numerical simulations are compared and validated with the recorded runs of pressure changes on a specially built test stand. 2. Mathematical model of the flow The aim of the present paper is to present simulations of pressure transients of the considered system in cases of sudden jumps of pressure at the end of a long hydraulic line caused by e.g. operational overload. The fundamental component of the hydraulic system, shown in Fig.1, is the liquid long line which is treated as a distributed parameter element. The system is a subassembly, which is often found in many hydraulic systems of technical machinery and devices used e.g. inmining or ship building industry. 2.1. Fundamental equations The unsteady flow in liquid pipes is often represented by two 1D hyper- bolic partial differential equations. Linearized equations of momentum and Simulation of transient flows in a hydraulic system... 855 Fig. 1. The object of investigations continuity can be given by Jungowski (1976), Ohmi et al. (1985) and Zarzycki (1994) ρ0 ∂v ∂t + ∂p ∂z + 2 R τw =0 ∂p ∂t +ρ0c 2 0 ∂v ∂z =0 (2.1) where: t is time, z – distance along the pipe axis, v= v(z,t) – average value of velocity in a cross-section of the pipe, p= p(z,t) – average value of pressure in a cross-section of pipe, τw – shear stress in the pipe wall, ρ0 – density of the liquid (constant), R – inner radius of the pipe. The speed of the acoustic wave c0 in Eq. (2.1)2 takes into consideration compressibility of the liquid and elastic deformability of the pipewall, and can be given by the following relation (Wylie and Streeter, 1978) c0 = √ βc ρc 1 √ 1+2 βc E R g (2.2) where βc is the bulkmodulus of the liquid, E – Young’s modulus of the tube and g – thickness of the wall. In an unsteady flow in the pipe, the instantaneous stress τw may be re- garded as a sum of two components: the quasi-steady state shear component and the unsteady state shear component (Ohmi et al., 1985; Zarzycki, 2000) τw = 1 8 λρ0v|v|+ 2µ R t∫ 0 w(t−u)∂v ∂t (u) du (2.3) where: λ is the Darcy-Weisbach friction coefficient, w – weighting function, µ – dynamic viscosity, u – time used in the convolution integral. Thefirst component inEq. (2.3) presents the quasi-steady state of thewall shear stress, the second one is the additional contribution due to unsteadiness. 856 Z. Zarzycki et al. The second summand inEq. (2.3) relates thewall shear stress to the instanta- neous average velocity and to the weighted past velocity changes. The system of Eqs. (2.1) and (2.3) is closed because of p and v, as long as the weighting function w(t) is known for the wall shear stress at the pipe wall. 2.2. Weighting functions Zielke (1968) was first to give an analytical relation for theweighting func- tion w(t) for a laminar flow. He derived it from the analysis of a transient two-dimensional laminar flow. The Zielke model can be given by w(t̂)=    6∑ i=1 mit̂ i−2 2 for t̂¬ 0.02 5∑ i=1 exp(−nit̂) for t̂ > 0.02 (2.4) where mi is 0.28209, −1.25, 1.05778, 0.93750, 0.396696, −0.351563 and ni = −26.3744, −70.8493, −135.0198, −218.9216, −322.5544, respectively, and t̂ is the dimensionless time, defined by the following relation t̂= ν t R2 (2.5) Zielke’smodel requiresmuch computermemory and, therefore, itwasmodified by Trikha (1975) and Schoohl (1993) to improve its computational efficiency. Schoohl’s model can be given by w(t̂)= 5∑ i=1 miexp(−nit̂) (2.6) where: m1 = 1.051, m2 = 2.358, m3 = 9.021, m4 = 29.47, m5 = 79.55, n1 =26.65, n2 =100, n3 =669.6, n4 =6497, n5 =57990. In the case of the unsteady turbulent flow, theweighting function depends not only on the dimensional time but also on the Reynolds number. Vardy et al. (1993) andVardy and Brown (1996) derived amodel in which all viscosity effects were assumed to occur in the steady boundary layer (viscosity varied linearly across the outer annular shear layer). An approximated form of their weighting function model is wa = 1 2 √ πt̂ exp ( − t̂ C∗ ) (2.7) Simulation of transient flows in a hydraulic system... 857 where C∗= 12.86 Rek k= log10 15.29 Re0.0567 Zarzycki (1994, 2000) developed a weighting function model using a four- region discretization (four instead of two regions) for turbulent viscosity di- stribution. The model had a complex mathematical form and it was further approximated to a simpler form wapr =C 1√ t̂ Ren (2.8) where: C =0.299635, n=−0.005535. Zarzycki’s model yields the same results as Vardy’s model, but generates themmore quickly. InEq. (2.8), both time andReynold’s number (i.e. also the speed) are in the denominator, whichmakes calculations difficult. In order to eliminate these difficulties, Zarzycki and Kudźma (2004) and Kudźma (2005) presented a model similar to Schohl’s model for a laminar flow. Their model can be given by wN(t̂)= (c1Re c2 + c3) 6∑ i=1 Aiexp(−bit̂) (2.9) where: A1 = 152.3936, A2 = 414.8145, A3 = 328.2, A4 = 640.2165, A5 = 58.51351, A6 = 17.10735, b1 = 207569.7, b2 = 6316096, b3 = 1464649, b4 = 15512625, b5 = 17790.69, b6 = 477.887, c1 = −1,5125, c2 = 0.003264, c3 =2.55888. The value of critical Reynold’s number (between unsteady laminar and turbulent flow), which qualifies the application of an appropriate weighting function, can be calculated bymeans of the following empirical relation (Ohmi et al., 1985) Recn =800 √ Ω (2.10) Equation (2.10) can be used for an oscillatory flow, whereas for a pulsating flow Recn is (Ramaprian and Tu, 1980) Recn =2100 (2.11) where: Ω=ωR2/ν denotes the dimensionless frequency, ω=2π/T – dimen- sional frequency, T =4L/c0 – period of the waterhammer, L – length of the pipe. For further simulations of the hydraulic waterhammer effect, two models were adopted: model (2.6) for the laminar flow andmodel (2.9) for the turbu- lent flow. 858 Z. Zarzycki et al. 2.3. Method of characteristics. Computational codes The system of Eqs. (2.1) and (2.3) with the known weighting function presents a closed nonlinear system of differential-integral Voltera’s equations with a degenerated kernel. The system can be transformed into a pair of ordinary differential equations usingMOC– amethod of characteristics. As a result, we obtain (Zarzycki and Kudźma, 2004; Zarzycki and Kudźma, 2005) ±dp+ρ0c0dv+ 2τwc0 R dt=0 dz=±c0dt (2.12) The net for the method of characteristics is shown in Fig.2. Fig. 2. The net of characteristics Equations (2.12) were approximated using a differential scheme of the first order.As itwas proved (ChaudhryandHussaini, 1985) such an approximation gives satisfactory results provided that the time step ∆t remains small.Owing to that, a system of algebraic equations was created. The system enables cal- culation of cross-sectional mean values of the instantaneous pressure and flow rate: — for the internal nodal points of the net of characteristics pk+1,i = 1 2 [ (pk,i−1+pk,i+1)+ρ0c0(vk,i−1−vk,i+1)+ + 2∆z R (τw(k,i+1)− τw(k,i−1)) ] (2.13) vk+1,i = 1 2 { (vk,i−1+vk,i+1)+ 1 ρ0c0 [ (pk,i−1−pk,i+1)+ + 2∆z R (τw(k,i+1)+ τw(k,i−1)) ]} Simulation of transient flows in a hydraulic system... 859 — for the boundary nodal points of the net of characteristics: pk+1,1 = pk,2+ρ0c0 [ (vk+1,1−vk,2)+ 2∆t R τw(k,2) ] (2.14) vk+1,h+1 = vk,h+ 1 ρ0c0 (pk,h−pk+1,h+1)− 2∆t R τw(k,h) where i = 2,3, . . . ,h, k = 1,2, . . . ,m, m is the number of time steps, h – number of calculation sections along the hydraulic line. As it was mentioned in Section 2.1, the instantaneous shear stress τw can be given by a sum of two components τw(k,i) = τwq(k,i)+ τwn(k,i) (2.15) where τwq(k,i) = 1 8 ρ0λ(Rek,i)vk,i|vk,i| (2.16) τwn(k,i) = 2µ R [(vk,i−vk−1,i)w1,i+(vk−1,i−vk−2,i)w2,i+ . . .+ +(v2,i−v1,i)wk−1,i] The friction loss coefficient λ in Eq. (2.16)1 is expressed for the laminar flow by λ= 64 Re (2.17) And for the turbulent flow, from Prandtl’s formula, by 1√ λ =0.869ln(Re √ λ)−0.8 (2.18) The system of Eqs. (2.13)-(2.18) together with Eqs. (2.6), (2.9) and the boun- dary and initial conditions is the basis for creating an algorithmof calculations and then a computer program. Figure 3 presents the algorithm of calculations. 2.4. Verification of the model In order to compare the accuracy of unsteady and quasi-steady models of friction in relation to experimental data, simulations of a simplewaterhammer case (tank – long liquid line and cut-off valve) were conducted. 860 Z. Zarzycki et al. Fig. 3. Flow chart Simulation of transient flows in a hydraulic system... 861 The computed results were compared with experimental data reported by Holmboe and Rouleau (1967). They ran tests on a copper tube with radius 0.0127m and length 36.1 m 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.7 ·10−6m2/s. Themeasured sound speed was 1324m/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 midpoint of the line. From the above para- meters, it followed that it was a case of a laminar flow. It was determined in numerical calculations in which themodels of Zielke, Eq. (2.4), and Schoohla (2.6) were used. In addition, the calculation with the quasi-steadymodel only i.e., with Eqs. (2.16)1 and (2.17) was done as well. Results of simulations and experimental data are shown in Fig.4. It is clearly seen that the calculation using the weighting functions (changeable hydraulic resistance) is much closer to the experimental data. Therefore, in further calculations the weighting functions were used instead of the quasi- steady model. Fig. 4. Fluctuations of pressure at the midpoint of the line; 1 – experimental data, 2 – simulations with the unsteadymodel of friction, 3 – simulations with the quasi-steadymodel 2.5. Boundary conditions It is necessary to knowboundary conditions in order to be able to solve the systemofEqs. (2.1) and (2.4) using themethodof characteristics. Theanalysis involves pressure runs at the beginning and at the end of the long line (Fig.1) during a sudden shift of the hydraulic directional control valve in time t0. 862 Z. Zarzycki et al. At that moment, a sudden change in the pressure from p0 to p0+∆p takes place. Two cases are investigated. The first case assumes a constant rate of delivery of the positive-displacement pump, the second one takes into account pulsation of the delivery of the pump. These conditions can be expressed in the following way: — for a constant rate of delivery for z=0 Q(t)= const (2.19) and for z=L p= { p0 for 0¬ t< t0 p0+∆p for t­ t0 (2.20) — for a changeable rate of delivery for z=0 Q=Qm [ 1− 1 2 K=∞∑ K=1 δQK cos(ωKt) ] (2.21) where: ωK denotes pulsations of harmonic vibration of the pump, K – order of harmonics, δQK – relative amplitudes of harmonic vibration of the flow rate according to the literature data, Qm – mean theoretical efficiency. Relation (2.21) was obtained by Rohatynski (1968). The condition for z=L has the form expressed in Eq. (2.20). 3. The test stand and description of experimental investigations In order to validate the presented model and the method for simulation of the hydraulic waterhammer effect, some tests were carried out on a specially prepared test stand. A diagram of the hydraulic system of the test stand is presented in Fig.5. The central part of the system comprised a hydraulic li- ne. Two extensometer pressure converters (7), (9) of the working liquid were fixed to its two ends. The generated flow intensity through the axial-flowmul- tipiston pump with deflected disc (6) Z-PTOZ2-K1-100R1 was measured by flowmeter (13). At the end of the hydraulic line, hydraulic control valve (10) (4/2) was installed, whose function was to suddenly direct the liquid through throttle valve (11). In order to realises an increase in the system load an adju- stable throttle valvewas used.Toprotect the systemagainst an incidental and dangerous pressure increase, safety valve (8) was installed right at the pump. Simulation of transient flows in a hydraulic system... 863 Fig. 5. A scheme of the hydraulic system used in experimental investigations on the pressure wave velocity propagation: 1 – hydraulic constant delivery pump (pressure charging pump), 2 – safety valve, 3 – filter, 4 – throttle valve, 5 – vacuometer, 6 – hydraulic variable displacement pump, 7 – pressure transducer, 8 – safety valve, 9 – pressure transducer, 10 – directional control valve 4/2, 11 – throttle valve, 12 – check valve, 13 – flowmeter, 14 – water cooler, 15 – thermometer The unsteady state in the system was caused by the shifted hydraulic directional control valve directing the liquid flow through the throttle valve with higher hydraulic resistance d2 (Fig.1). The time of shift of the directional control valve was tz = 20ms. It was shorter than half of the hydraulic hammer time (tz