© Strojni{ki vestnik 49(2003)3,137-149 © Journal of Mechanical Engineering 49(2003)3,137-149 ISSN 0039-2480 ISSN 0039-2480 UDK 621.227:519.87:004.43 UDC 621.227:519.87:004.43 Izvirni znanstveni ~lanek (1.01) Original scientific paper (1.01) Matemati~no modeliranje hidravli~nega ovna Mathematical Modelling of a Hydraulic Ram Pump System Veljko Filipan, Zdravko Virag, Anton Bergant V prispevku podajamo matematični model hidravličnega ovna (HO), vgrajenega v cevni sistem. Obravnavamo poenostavljeni delovni postopek. Pretok v sistemu z ovnom je neustaljen, zato podajamo enačbe neustaljenega toka v ceveh in metodo karakteristik (MK). Podan je podrobni oris matematičnega modeliranja elementov sistema HO. Pogonsko in dobavno cev modeliramo s pomočjo pravokotne mreže karakteristik, pogonski in dobavni zbiralnik pa sta robna pogoja z nespremenljivo gladino vode. HO modeliramo kot robni pogoj, ki popisuje delovanje ovna v področjih delovnega postopka. Ta robni pogoj je popisan z 11 enačbami in 11 odvisnimi spremenljivkami, rešujemo ga z metodo iteracij v časovnih korakih računskega postopka. Postavljen model je programiran na računalniku. V podani računski simulaciji smo uporabili podatke iz literature, časovni potek tlaka in pretočne hitrosti so prikazani v grafični obliki. © 2003 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: ovni hidravlični, sistemi črpalni, udari tlačni, metode karakteristik) In this paper we present the mathematical modelling of a hydraulic ram rump (HRP) system and an explanation of the simplified working cycle of its action. The flow in the HRP system is unsteady; therefore, unsteady pipe-flow equations and the method of characteristics (MOC) are given. The mathematical modelling of particular components of the HRP system is explained in detail. The drive and delivery pipes are modelled by a fixed grid MOC, and the supply and delivery tanks are the boundary conditions - open reservoirs with a constant water level. The HRP is modelled as a unit boundary condition describing the physics of its action. This boundary condition consists of 11 equations with 11 dependent variables and is solved by an iterative procedure for each time step of the computational run. The derived model is programmed for a digital computer. A computer simulation using input data from the literature is made and the results are presented in the form of pressure and velocity vs. time graphs. © 2003 Journal of Mechanical Engineering. All rights reserved. (Keywords: hydraulic rams, pumping systems, water hammer, method of characteristics - MOC) 0 UVOD Sistem s hidravličnim ovnom (HO) je preprost in obsega dve cevi, dva ventila in tlačni kotel. Hidravlični oven izkorišča energijo vodnega toka pri relativno majhni pogonski višini za črpanje ustrezne pretočne količine na višjo dobavno raven. HO obratuje brez dodatnih zunanjih pogonskih naprav (električnih ali mehaničnih), prav tako ni potreben dodatni vir energije. Preprosta izvedba in zanemarljivi obratovalni stroški sta dva ključna elementa za uporabo sistema s hidravličnim ovnom v moderni inženirski praksi, več ko 200 let po njegovem izumu. V začetnem obdobju razvoja od leta 1797 (leto, ko je J. M. Montgolfier izumil HO) do konca 19. stoletja so številni raziskovalci neuspešno poskušali postaviti teoretični model delovanja HO. Žukovskij je leta 1898 [1] postavil teoretični model vodnega udara, po 0 INTRODUCTION The hydraulic ram pump (HRP) system is a very simple device that consists of two pipes, two valves and a pressure vessel. The momentum of the water flow under a relatively low drive head is used to pump a small portion of the flow to a considerably higher delivery head. The HRP operates without any other external (electrical or mechanical) driving device, and with no other external energy input. The simplicity of its structure and the negligible operating costs are the two main reasons why the HRP is currently very interesting, although it was invented more than 200 years ago. There have been many unsuccessful theoreti-cal attempts to explain the HRP’s action, from its in-vention in 1797 (by J. M. Montgolfier) until the end of the 19th century, when the water hammer – the main use of the HRP’s action – was first explained in detail gfin^OtJJIMISCSD 03-3 stran 137 | ^BSSITIMIGC Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling katerem deluje hidravlični oven. Delovanje hidravličnega ovna je bilo zadovoljivo popisano na podlagi sklopljene teoretične in eksperimentalne analize HO v prvi polovici 20. stoletja. Najpomembnejši izvajalci sklopljenih raziskav v tem obdobju so bili: O’Brien in Gosline [2], Lansford in Dugan [3] ter Krol [4]. Omeniti moramo tudi Eytelweina [5] in Calverta [6], ki sta izvedla obširne eksperimentalne raziskave HO. Na podlagi analize časovnega poteka tlaka sta O’Brien in Gosline [2] leta 1933 prva uspešno popisala načelo črpanja z uporabo HO. Kakorkoli že, raziskovalca sta v svojem modelu uporabila teorijo elastičnega vodnega udara samo v črpalnem taktu. Preostale takte delovnega postopka HO sta popisala s teorijo togega vodnega udara. Podoben model z malenkostnimi izboljšavami sta razvila Lansford in Dugan [3]. Krol [4] je uporabil teorijo togega udara za popis celotnega delovnega postopka HO. Obravnavani avtorji [2] do [4] so v matematičnih modelih uporabili statične karakteristike ventilov in cevi. Poenostavljen model togega vodnega udara je razvil tudi Iversen [7]. Basfeld in Müller [8] sta predstavila podoben model ob uporabi izmerjenih koeficientov pretočne hitrosti. Chiang [9] je uporabil impedančno metodo za popis delovanja in optimizacijo HO. Rennie in Bunt [10] sta razvila grafoanalitični model, ki sloni na grafoanalitični rešitvi enačb neustaljenega toka v cevi za praktično celoten delovni postopek HO. Rezultati izračuna so se dobro ujemali z rezultati meritev. Sam postopek grafičnega reševanja v ravnini piezometrična višina (h) - pretočna hitrost (u), ki je fizikalno jasen, pa je težaven in zamuden. Avorja sta objavila tudi rezultate svojih obsežnih eksperimentalnih raziskav [11]. Razvoj na koncu 20. stoletja je bil usmerjen v izdelavo preprostih izvedb HO, namenjenih za proizvodnjo v skromno opremljenih delovnih obratih dežel v razvoju. Levji delež raziskav so opravili na Univerzi v Warwicku v obdobju med 1985 in 1997. Na tej univerzi so razvili niz izvedb hidravličnih ovnov in postavili poenostavljeni algebraični model [12], da bi čim bolj razširili uporabo HO v nerazvitih državah. V letih 1982 do 1984 so tudi na Tehnični univerzi v Delftu potekale obsežne laboratorijske raziskave, ki so obsegale preizkus 12 hidravličnih ovnov iz industrije. Rezultate laboratorijskih raziskav je objavil Tacke [13]. V tem poročilu Tacke podaja matematični model, ki sloni na teoriji togega vodnega udara z izjemo črpalnega takta. Črpalni takt je popisan s teorijo elastičnega udara. Leta 1999 so Najm, Azoury & Piasecki [14] razvili numerični model HO, ki sloni na metodi karakteristik. Model zajema vse sestavne dele naprave, razviti program pa omogoča modeliranje vseh pretočnih režimov HO. V ^BSfiTTMlliC | stran 138 by Zhukovsky in 1898 [1]. Combined theoretical and experimental investigations in the first half of the 20th century led to an adequate explanation of the HRP’s action. In a combined approach O’Brien & Gosline [2], Lansford & Dugan [3] and Krol [4] made the most important contributions. Eytelwein [5] and Calvert [6] should also be mentioned because of their exhaustive and extensive experimental investigations of the HRP. In 1933 using experimental data of pressure vs. time O’Brien & Gosline [2] were the first to successfully explain the HRP’s pumping action. However, in their model they used the elastic water hammer theory only for the pumping period. For other parts of the HRP’s working cycle they considered using the rigid column theory. A similar approach, with some improvements, was reported by Lansford & Dugan [3]. Krol [4] used the rigid column theory for the whole cycle of the HRP’s action. All these authors [2] to [4] used the measured static characteristics of valves and pipes in their mathematical models. A simplified rigid column theory model was also derived by Iversen [7]. Basfeld & Müller [8] introduced a similar model using experimentally obtained velocity coefficients. Chiang [9] used a model for the simulation and optimisation of a HRP based on an impedance method. Rennie & Bunt [10] derived a graphic-analytical model and they used unsteady pipe-flow equations and their graphic-analytical solution for almost the whole of the HRP cycle. The results of the model agreed well with their experimental results, but the graphical procedure in the plane of the piezometric head (h) and the velocity (u), although physically clear, is very cumbersome. The authors presented comprehensive experimental results with all the necessary details in [11]. Developments at the end of the 20th century were directed towards simplified solutions of the HRP so that they could be manufactured in the workshops of developing countries. Great efforts were made at the University of Warwick, where they ran a project from 1985 to 1997. They developed several designs of HRP and introduced a simplified algebraic model [12], all with the aim of enabling a wider application of the HRP in developing countries. An extensive investigation with a large number of laboratory tests on twelve commercially available hydraulic rams was also carried out at the Delft University of Technology in the period 1982 to 1984. The results of the laboratory investigations were reported by Tacke in 1988 [13]. In this report Tacke developed a mathematical model based on rigid column theory, except for the pumping period. This period was considered using the pressure-wave theory. In 1999 Najm, Azoury & Piasecki [14] developed a numerical model of the HRP based on the method of characteristics. The model encompasses all the parts of the device and the derived program accounts for all possible flow situations of the HRP. This is very similar to the model of the HRP system developed in a doctoral Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling pričujočem prispevku podajamo podoben model, ki ga je v doktorskem delu razvil prvi avtor [15]. Model sloni na teoriji elastičnega vodnega udara in zajema vse elemente HO. Obravnavani model je programiran na računalniku, v numerični simulaciji pa smo uporabili podatke Rennieja in Bunta [11]. Dobljeni rezultati izračuna so primerjani z rezultati meritev iz literature. 1 DELOVNI POSTOPEK HIDRAVLIČNEGA OVNA Slika 1 prikazuje značilni sistem z vgrajenim HO. Sistem sestavljajo: HO oziroma črpalka (1), pogonska cev (2), dobavna cev (3), pogonski zbiralnik (4) in dobavni zbiralnik (5). Črpalko (1) sestavljajo: okrov (1a), vzbujevalni ventil (1b), dobavni ventil (1c) in tlačni kotel (1d) z zračno blazino. Oven izrablja spremembo kinetične energije ustrezno velike količine vode (pogonski pretok, Q), ki se pretaka iz pogonskega zbiralnika (pogonska višina, H), pri hitri zapori vzbujevalnega ventila. Narastek tlaka omogoči črpanje vode (dobavni pretok, q) v dobavni zbiralnik (dobavna višina, h). Obratovanje HO sloni na izmeničnem delovanju obeh ventilov. Vzbujevalni ventil se odpre z uporabo stalno delujoče sile uteži ali vzmeti. Obravnavani ventil avtomatično zapre hidrodinamična sila delujoča na ventil. Dobavni ventil je preprosta izvedba povratnega ventila, ki odpira in zapira na podlagi razlike tlakov. V prvem približku lahko delovni postopek HO razdelimo na tri značilne delovne takte, kakor je prikazano na sliki 2. thesis by the first author [15] and outlined in this paper. The model is based on elastic water-hammer theory and takes into account all the parts of the HRP system. The HRP unit is modelled as a complex boundary condition and a significant extension is made in the impulse valve modelling. The model is programmed for a digital computer and a simulation using the experimental data of Rennie and Bunt [11] is made. The obtained results are compared with experimental data from the literature. 1 THE WORKING CYCLE OF THE HRP A typical HRP system is shown in Fig 1. It consists of a HRP unit or pump (1), drive pipe (2), delivery pipe (3), supply tank (4) and delivery tank (5). The pump (1) consists of a pump body or casing (1a), impulse valve (1b), delivery valve (1c) and a pressure vessel (1d) partially filled with air. The HRP operates by extracting the energy from a large quantity of water (drive flow, Q), which flows from a supply tank (driv-ing head, H), by rapid closure of the impulse valve. A small quantity of water (delivery flow, q) is then pumped to a delivery tank (delivery head, h). The operation is due to the automatic action of two valves that work intermittently. The impulse valve is open due to a con-tinuously acting force (dead weight or spring). The valve is closed automatically under the dynamic force of flowing water. The delivery valve is a simple non-return valve that is opened and closed under a pres-sure difference. In a simplified manner the HRP working cycle can be divided into three different periods of action, as depicted in Fig. 2. 3 t l~L_5. q 1a 1c 1b] Sl 1. Značilni sistem s HO (1-črpalka: 1a-okrov, 1b-vzbujevalni ventil, 1c-dobavni ventil, 1d-tlačni kotel; 2-pogonska cev; 3-dobavna cev; 4-pogonski zbiralnik; 5-dobavni zbiralnik) Fig 1. A typical HRP system (1-pump: 1a-casing, 1b-impulse valve, 1c-delivery valve, 1d-pressure vessel; 2-drive pipe; 3-delivery pipe; 4-supply tank; 5-delivery tank) ^vmskmsmm 03-3 stran 139 | ^BSSITIMIGC Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling u u u -u Sl. 2. Poenostavljen delovni postopek HO Fig. 2. The simplified working cycle of an HRP V pospeševalnem taktu (T) teče voda skozi odprt vzbujevalni ventil, povratni ventil je zaprt. Voda v pogonski cevi je pospešena od u=0 do u=u (faza 1) in nato na u=u (faza 2). V prvi fazi je vzbujevalni ventil polno odprt, v drugi fazi pa se zapira zaradi hidrodinamične sile delujoče na ventil. V sklepnem delu pospeševalnega takta se vzbujevalni ventil hitro zapre, pretok vode se zaustavi, zaustavitev pa povzroči vodni udar z visokim tlakom. Tlak vodnega udara v okrovu črpalke je večji od tlaka v tlačnem kotlu. Posledično se nato odpre dobavni ventil, vodni tok se usmeri proti tlačnemu kotlu in naprej v dobavni cevovod, kjer se na koncu izteka v dobavni zbiralnik. V drugem taktu postopka, imenovanem črpalni takt (T), voda teče skozi odprt dobavni ventil, medtem ko je vzbujevalni ventil zaprt. Pretočna hitrost v sistemu se korakoma zmanjšuje glede na dinamiko potovanja udarnih valov v pogonski cevi. Tlak v okrovu črpalke se postopoma zmanjša na kritični tlak, pri katerem se zapre dobavni ventil. Temu sledi odprtje vzbujevalnega ventila z uporabo sile uteži ali vzmeti in pojav povratnega toka v pogonski cevi (voda se pretaka iz okrova črpalke v pogonsko cev proti pogonskemu zbiralniku). Pretok v pogonski cevi postopoma pojema od u=-u na u=0, obravnavani tretji takt imenujemo sunek (Tr). Sunek konča delovni proces hidravličnega ovna, ki se nato samodejno ponavlja s periodo T. 2 ENAČBE NEUSTALJENEGA TOKA V CEVI IN METODA KARAKTERISTIK Tok v sistemu s HO lahko popišemo z enačbami neustaljenega toka v cevi [16]. To sta kontinuitetna enačba: In the acceleration period (Ta) water discharges through the opened impulse valve while the delivery valve is closed. Water in the drive pipe is accelerated from u=0 to u=uv (phase 1) and after that to u=uc (phase 2). In phase 1 the impulse valve is fully opened and in phase 2 it is closing due to the dynamic force of the water. At the end of the acceleration period the impulse valve is closed rapidly, the flow is stopped, and a water hammer with high pressure is created as a consequence. The induced pressure in the pump body is higher than the pressure in the air chamber. As a consequence, the delivery valve opens, enabling water to flow to the pres-sure vessel and further through the delivery pipe to the delivery tank. In this second part of the cycle, called the pumping period (Tp), water continues to flow through the opened delivery valve while the impulse valve is closed. The flow velocity diminishes stepwise after each returning of the subsequent pressure waves that are travelling along the drive pipe. Gradually the pressure in the pump body drops below the critical pressure, causing the delivery valve to close. After this the impulse valve opens due to active force and negative flow in the drive pipe (from the pump body to the supply tank) occurs. The flow is gradually decelerated from u=-ur to u=0 in this part of action, called the recoil period (Tr). The end of the recoil period is the beginning of the second cycle and the action of the HRP contin-ues automatically with a period of Tc. 2 THE UNSTEADY-PIPE-FLOW EQUATIONS AND THE METHOD OF CHARACTERISTICS The flow in a HRP system can be described using unsteady-pipe-flow equations [16]. These are the continuity equation: (1) VH^tTPsDDIK in gibalna enačba: stran 140 dh/dt + udh/dx + u sinß + c2/gdu/dx = 0 and the momentum equation: du/dt + udu/dx + gdh/dx + X\u\u/(2D) = 0 (2), Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling kjer so: h - piezometrična višina, t - čas, u - pretočna hitrost, x - koordinata, ß - strmina cevovoda, c -hitrost vala, g - zemeljski pospešek, X - Darcy-Weisbachov koeficient trenja in D - premer cevi. Pri izpeljavi enačb (1) in (2) smo prevzeli naslednje postavke: enakomerna porazdelitev tlaka, hitrost je povprečena po prerezu cevi, tlačni tok v cevi, cev je okroglega prereza, elastične napetosti v steni cevi, Newtonska tekočina. Par kvazilinearnih parcialnih diferencialnih enačb (1) in (2) nima analitične rešitve, za reševanje obravnavanih enčb so razvite števile metode ([16] in [17]). V prispevku bomo uporabili metodo karakteristik (MK). Sprememba enačb (1) in (2) po metodi karakteristik da štiri navadne diferencialne enačbe, dve združljivostni enačbi: in which h is the piezometric head, t is the time, u is the flow velocity, x is the distance, b is the pipe slope, c is the wave speed, g is the gravitational acceleration, l is the Darcy-Weisbach friction factor and D is the pipe diameter. Equations (1) and (2) are derived using the following assumptions: the distribution of pressure is uniform, the velocity is averaged over the pipe’s cross section, the pipe is full with fluid, the pipe is of circular cross section, the pipe strains are in the elastic domain, the fluid is Newtonian. The pair of quasi-linear partial differential Equa-tions (1) and (2) has no analytical solution and a number of methods for their solution have been developed ([16] and [17]). The method of characteristics (MOC) is used in this paper. The MOC transformation of equa-tions (1) and (2) results in four ordinary differential equations, two compatibility equations: ± g/c dh/dt + du/dt ± g/c u sinß + X \u\ u/(2D) = 0 (3) in dve enačbi karakteristik: and two characteristic equations: dx/dt = u ± c (4). Vsaka združljivostna enačba (3) je veljavna samo vzdolž ustrezajoče karakteristične enačbe (4). Člen g/c u sinß lahko zanemarimo, ker je hitrost vala c mnogo večja od pretočne hitrosti u. Enačbi integriramo z metodo končnih razlik ([16] in [17]). Rezultat integracije je sistem algebraičnih enačb. Te enačbe podajamo v poglavju 3 skupaj z enačbami robnih pogojev v sistemu s hidravličnim ovnom. 3 MATEMATIČNO MODELIRANJE SISTEMA S HIDRAVLIČNIM OVNOM Elementi sistema s HO z uporabljenimi označbami so prikazani na sliki 3. Each compatibility Equation (3) is only valid along the corresponding characteristic Equation (4). The term g/c u sinb in Eqn. (3) can be neglected be-cause the wave speed c is much higher than the flow velocity u. The equations are integrated using a dis-crete finite-difference method ([16] and [17]). As a result a set of algebraic equations is obtained. These equations are presented in Chapter 3 together with the boundary conditions of the HRP model. 3 MATHEMATICAL MODEL OF THE HRP SYSTEM The components of the HRP system with symbols used in model are shown in Fig 3. Sl. 3. Elementi modela sistema s HO Fig. 3. The components of HRP system model gfin^OtJJlMlSCSD 03-3 stran 141 | ^BSSITIMIGC Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling Modeliranje pogonske in dobavne cevi V numeričnem modelu [13] razdelimo pogonsko cev na izbrano število cevnih odsekov I (NE=1 do I) enake dolžine Dx. Dobavna cev je razdeljena na ustrezno število cevnih odsekov NEMAX-(I+1) (NE=J do NEMAX-1). Vsak cevni odsek je omejen z dvema vozliščema, ki ju oštevilčimo z NC=1 do NCMAX. V numeričnem izračunu uporabimo pravokotno mrežo MK. Časovni korak Dt izberemo na podlagi Courant-Friedrich-Lewyjevega kriterija stabilnosti: Modelling the Drive and Delivery Pipes In the numerical model [15], the drive pipe is divided into a selected number of pipe reaches I (NE=1 to I) of equal length Dx. The delivery pipe is divided into a corresponding number of pipe reaches NEMAX-(I+1) (NE=J to NEMAX-1). Each pipe reach is bounded by two nodes numbered NC=1 to NCMAX. A fixed grid of the MOC is used for the numerical computation. The time step Dt is selected according to the stability criterion of Courant-Friedrich-Lewy: Dt\jPr-> P 1 P2 , , t0 R/ O1 O2 S ¦ 1 O Dx Dx 0 Sl. 4. Ravnina x-t Fig. 4. The x-t plane VH^tTPsDDIK gfin^(5ül[Fi!]DaGC] | stran 142 x Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling rt0+At t0+Dt t0 NEMAX-1 NC=NCMAX-2 NC=NCMAX-1 Sl. 5. Pogonski (a) in dobavni (b) zbiralnik Fig. 5. Supply (a) and delivery (b) tanks v odprtem zbiralniku (H=konst. in h=konst, sl. 5). Pretočni hitrosti uP2 in uP izračunamo iz ustreznih enačb (8) in (7). Enotni robni pogoj za HO Hidravlični oven (črpalka) (sl. 1), sestavljajo ga okrov, vzbujevalni in dobavni ventil ter tlačni kotel, modeliramo kot enotni robni pogoj. Nespremenljivi tlak (višino) na izstopnem delu vzbujevalnega ventila v računski točki NC=NCMAX modeliramo kot zbiralnik s stalno gladino vode (zIS =konst, sl. 3). Enačbe, ki popišejo druge dele obravnavanega robnega pogoja, so podane kakor sledi. (1) Na navzdolnjem koncu pogonske cevi, ki je pritrjena na okrov ovna (NC=N), lahko izrazimo zdruljvostno enačbo C+ karakteristike v naslednji obliki: reservoirs with a constant water level (H=const. and h=const., see Fig. 5). The velocities uP2 and uP1 are cal-culated from the corresponding Equations (8) and (7). The HRP Unit Boundary Condition The HRP unit (pump) (see Fig. 1) that con-sists of a pump body, impulse and delivery valves, and a pressure vessel, is modelled as a unit bound-ary condition. Constant pressure (head) at the outlet of the impulse valve at the computational node NC=NCMAX is modelled as an air-open outlet reservoir with a constant water level (zIS=const., see Fig. 3). The equations that describe the other parts of the HRP boundary condition are explained below. (1) At the downstream-end section of the drive pipe coupled to the HRP pump body (NC=N) the compatibility equation of the C+ characteristic can be expressed as: ob upoštevanju R=NI po sliki 3. Neznanki v enačbi (11) sta višina v okrovu ovna N (hN) in pretočna hitrost na navzdolnjem koncu pogonske cevi (u I) v računskem času t0+At (2) Na navzgornjem koncu dobavne cevi, ki je pritrjena na tlačni kotel (NC=NK), lahko izrazimo združljivostno enačbo C" karakteristike v naslednji obliki: g/cR (hN-hR) + u+I - uR + FR (u+I+uR) = 0 (11) with R=NI according to Fig. 3. The unknowns in Eq. (11) are the head in the pump body (hN) and the velocity at the downstream end of the drive pipe (u+I) for the new time step t0+Dt. (2) At the upstream-end section of the deliv-ery pipe coupled to the pressure vessel (NC=NK) the compatibility equation of the C- characteristic can be expressed as: ob upoštevanju S=NJ po sliki 3. Neznanki v enačbi (12) sta višina na priključku v tlačni kotel (hNK) in pretočna hitrost na navzgornjem koncu dobavne cevi (u J) v računskem času t0+At. (3) Enačbo politrope za zračno blazino v kotlu izrazimo v obliki: -g/cS (hNK-hS) + u-J - uS + FS (u-J+uS) = 0 (12) with S=NJ according to Fig. 3. The unknowns in Eq. (12) are the head just below the pressure vessel (hNK) and the velocity at the upstream-end section of the delivery pipe (u-J) for the new time step t0+Dt. (3) The polytropic equation for air in the pres-sure vessel can be expressed as: p V n = p V n = konst. p0 p0 (13), | IgfinHŽslbJlIMlIgiCšD I stran 143 glTMDDC t Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling kjer so: n eksponent politrope, pp in V sta neznana absolutni tlak in prostornina zraka v tlačnem kotlu v računskem času t+Dt, absolutni tlak p 0 in prostornina zraka V sta izračunana v predhodnem računskem času t. p0 (4) Absolutno piezometrično višino v vozlišču NK izrazimo z enačbo : where n is a polytropic exponent, pp and Vp are the unknown absolute pressure and the volume of air in the pressure vessel for the new time step t0+Dt, while pp0 and Vp0 are the known absolute pressure and the volume of air for the previous time step t0. (4) The absolute piezometric head at the node NK is expressed as : kjer so r - gostota tekočine, zTP in zIS - gladini vode (sl. 3). (5) Model dobavnega ventila upošteva izmerjene koeficiente tlačnih izgub na ventilu x =f(uK). Ventil je odprt, ko je razlika tlakov med vozliščema N in NK pozitivna. V nasprotnem primeru, tj. pri negativni razliki tlakov, je ventil zaprt. Pretočno hitrost skozi ventil popišemo takole: (6) Kontinuitetno enačbo za računsko vozlišče NK zapišemo v obliki: hNK = pp + r g (zTP-zIS) (14), in which r is the fluid density, and zTP and zIS are the water-level elevations (see Fig. 3). (5) The delivery valve is modelled using the measured valve head-loss coefficient xDV=f(uK). If the pressure difference between points N and NK is positive, the valve is opened. On the other hand, if the pressure difference is negative, the valve is closed. The velocity through the valve is modelled as fol-lows: )0.5 za/for hN > h za/for^ 0 za/for sPV = 0 (20) ob upoštevanju NM=NCMAX po sliki 3. (8) Za rešitev sistema enačb (11) do (20) potrebujemo še kontinuitetno enačbo, zapisano za vozlišče N: with NM=NCMAX according to Fig. 3. (8) One additional equation necessary for solving Eqs. (11) to (20) is the continuity equation for node N: u+I AI - uM AM - uK AK = 0 MM KK (21). Sistem 11 nelinearnih algebraičnih enačb (11) do (21) z 11 neznankami (h, hNK, u , u , uK, uM, s, vPV, aPV, p in V) rešimo z iterativno metodo. Uporabimo prerezno metodo [18], ki ustreza fiziki delovanja HO. Obravnavana metoda je preprosta in robustna (stabilna in konvergentna) [15]. Pri zagonu HO v času t=0. s postavimo hidrostatične tlačne razmere v sistemu. Dobavni ventil je zaprt, vzbujevalni ventil pa je polno odprt (nenadno odprtje). Temu sledi pospeševanje vodne mase v pogonski cevi po sliki 2. 4 RAČUNALNIŠKI PROGRAM IN REZULTATI IZRAČUNA Obravnavani matematični model programiramo v FORTRANu. Uporabimo standardno pravokotno mrežo MK za reševanje enačb neustaljenega toka v cevi. Razvili smo izvirni podprogram za HO, kjer v vsakem časovnem koraku izračunamo 11 neznank v enačbah (11) do (21) z uporabo iteracije. Za izračun tlakov (piezometričnih višin) in pretočnih hitrosti v računskih točkah uporabimo standardne podprograme. Izvedli smo računalniško simulacijo delovanja HO ob uporabi eksperimentalnih podatkov, dobljenih v prispevku Rennieja in Bunta ([10] in [11]). Glavni podatki so naslednji: pogonska višina //=1,52 m, dobavna višina Ä=32 m (ta višina je bila krmiljena z ventilom, vgrajenim v dobavni cevi), dolžina pogonske cevi L=12,2 m, premer pogonske cevi D=0,036 m, dolžina dobavne cevi 1=1,36 m, premer The system of 11 non-linear algebraic Equa-tions (11) to (21) with 11 unknowns (hN, hNK, u+I, u-J, uK, uM, sPV, vPV, aPV, pp and Vp) is solved by an iterative method. The bisection method [18] based on the phys-ics of the HRP action is found to be a simple and robust one (stable and converging) [15]. For the HRP start-up at time t=0 it can be as-sumed that the water is under hydrostatic pressure in the system. The delivery valve is closed and the impulse valve is fully opened (it opens instantane-ously). Then the water in the drive pipe starts to accelerate according to Fig. 2. 4 THE COMPUTER PROGRAM AND THE SIMULATION RESULTS The described mathematical model is programmed in the FORTRAN language. A standard fixed-grid MOC procedure for solving the unsteady pipe flow is used. An original subroutine for modelling the HRP is developed, in which an iterative procedure is used for each time step to predict 11 unknowns from Eqs. (11) to (21). Standard subroutines are used for the computation of pressures (piezometric heads) and the velocities at other computational sections. As an example, the computer simulation was carried out using the experimental pump system from the paper of Rennie & Bunt ([10] and [11]). The most important data of this HRP system are as follows: drive head H=1.52 m, delivery head h=32 m (the head was adjusted by a control valve in the delivery pipe), length of drive pipe L=12.2 m, diameter of drive pipe D=0.036 m, length of delivery pipe L=1.36 m, diameter stran 145 Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling dobavne cevi D=0,012 m, največji pomik vzbujevalnega ventila sPVm =0,00381 m, potrebna sila za odprtje vzbujevalnega ventila F =2,68 N, prečni prerez vzbujevalnega ventila APV=0,001134 m2, vplivna masa gibajočih delov vzbujevalnega ventila m =0,273 kg, prostornina tlačnega kotla VT =0,0083 m3, stalnica tlačnega kotla pp0 Vp0n =380 Nm8/5, eksponent politrope n=1,2. Izmerjeni koeficienti dobavnega in vzbujevalnega ventila so približni, kakor sledi ([10], [11] in [15]): xPV=0,00761 s 138, C =0,000216 s-1,824 in x =200=konst. Ocenjena hitrost tlačnega vala v ceveh znaša c=1390 m/s [10]. Pogonsko cev smo razdelili na 36 cevnih odsekov dolžine Dx*0,339 m, dobavno cev pa na 4 odseke enake dolžine. Časovni korak v izračunu smo določili s pomočjo enačbe (5), dobimo vrednost Dt*0.000242 sekunde. Rezultati izračuna postavljenega modela so primerjani z rezultati izračuna in meritev Rennie in Bunta ([10] in [11]). Rezultati se dobro ujemajo, kar of delivery pipe D=0.012 m, maximum displacement of impulse valve sPVm =0.00381 m, force on impulse valve F =2.68 N, impulse valve cross-sectional area APV=0.001134 m2, effective mass of impulse-valve moving parts m =0.273 kg, the volume of the pressure vessel VT =0.0083 m3, the pressure vessel constant pp0 Vn =380 Nm8/5, the polytropic exponent n=1.2. The measured coefficients of the impulse and delivery valves are approximated as follows ([10], [11] and [15]): x =0.00761 s-1.38, C =0.000216 s-1.824 and x =200=const. The wave speed in the pipes was estimated as c=1390 m/s [10]. The drive pipe was divided into 36 pipe reaches of length Dx*0.339 m and the delivery pipe into 4 reaches of the same length. The time increment for the computational run was estimated according to Eq. (5) as Dt*0.000242 seconds. The simulation results of the presented model are compared with the computational and experimental results of Rennie & Bunt ([10] and [11]), and good agreement, especially in the case of the HRP 1.25 1 0.75 ^*^ ^^ yrtyr !_¦¦ ¦ ¦ ¦ ¦ ¦ ^ ¦ ^*r. R.&BU NT M. NT C. 0.5 L 0.25 0 0 1 02 03 04 05 06 07 08 09 01 00 h, m Sl. 7. Perioda T v odvisnosti od dobavne višine h (meritev - R.&BUNT M - in izračun - R.&BUNT C -c Rennie in Bunta; izračun s postavljenim modelom - MODEL) Fig. 7. The cycle period T vs. delivery head h (measured - R.&BUNT M - and computed - R.&BUNT C - by Rennie & Bunt; computed by described model - MODEL) 0.8 0.6 0.4 >^ .„^ : \ ¦ Xl ^ 1 , 1 W lPJ j\ .ill/ .# ¦« Mill I« ,JJJF • 8.6 8.8 9.2 9.4 9.6 Sl. 8. Pretočna hitrost na navzdolnjem koncu pogonske cevi (točka +I) Fig. 8. Velocity at the downstream end of the drive pipe (section +I) 9.8 t,s grin^sfcflMISDSD ^BSfiTTMlliC | stran 146 0.2 0 -0.2 9 Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling 5 4 A JUUVAA-V* „*m JJJJJJJJJJ IJXLU.UJJJJJJJJjUJOJJW -JWMMMM^WWWWl ^ t,s Sl. 9. Tlak na navzdolnjem koncu pogonske cevi, spojene z okrovom HO (vozlišče N) Fig. 9. Pressure at the downstream end of the drive pipe and pump body (node N) posebno velja za potek karakteristik HO [15]. Slika 7 prikazuje primerjavo rezultatov izračuna periode (T) v širokem pasu dobavnih višin (h), dobljene z uporabo postavljenega modela (MODEL) ter Rennie in Buntovega modela (R.&BUNT C) z rezultati meritev Rennie in Bunta (R.&BUNT M). Iz rezultatov izluščimo načelo delovanja in zapletene hidravlične interakcije v samem HO. Kot primer na slikah 8 in 9 podajamo rezultate izračuna časovnenega poteka pretočne hitrosti in tlaka na navzdolnjem koncu pogonske cevi za en delovni proces hidravličnega ovna (dobavna višina h=32 m). S slik 8 in 9 razberemo, da tlačne in hitrostne motnje v črpalnem taktu ovna vplivajo na celotni delovni postopek. Računalniški program računa v vsakem delovnem postopku naslednje karakteristike sistema s HO: pogonski pretok (Q), dobavni pretok (q), periodo postopka (T) in izkoristek (h). Izkoristek sistema s HO (h) je definiran kot Rankinov izkoristek: characteristics, is obtained [15]. Fig. 7 shows a comparison between the computational results for the cycle period (Tc) in a wide range of delivery heads (h) obtained by the described model (MODEL) and the model of Rennie and Bunt (R.&BUNT C), and the results of the measurements (R.&BUNT M). The results show the operational principle and the complex hydraulic interactions occurring within the HRP. As an example, the simulation results of velocity and pressure versus time traces at the downstream end of the drive pipe are depicted in Figs. 8 and 9, respectively, for one complete cycle of the HRP action (delivery head h=32 m). Figs. 8 and 9 show that the pressure and velocity disturbances occurring in a pumping period influence the whole cycle of operation. The computer program calculates the following HRP system characteristics: drive flow (Q), delivery flow (q), cycle period (Tc) and efficiency (h) for each successive cycle. The efficiency of the HRP system (h) is defined as the Rankine efficiency: h = (h/H - 1) q/Q (22), kjer sta: H - pogonska višina in h - dobavna višina (sl. 1). V občutljivostni analizi smo spremljali vpliv mase vzbujevalnega ventila (v obravnavanem primeru sile F, ki deluje na ventil z utežjo). Rezultati izračuna na v sliki 10 prikazujejo odvisnost karakteristik hidravličnega ovna Q, q, T in h od sile F. Ti rezultati so zelo uporabni v dejanskem primeru, ko moramo ustrezno nastaviti parametre sistema s HO. V splošnem zmanjšanje sile F vpliva na zmanjšanje pretočne količine (Q in q) in periode delovnega postopka (T), poveča pa se izkoristek (h). V primeru, ko mora sistem s HO dobaviti večje količine vode, povečamo silo F , s tem pa zmanjšamo izkoristek. Povedati moramo, da je pas nastavitve sile F omejen za vsak posamezni sistem s HO. where H is the driving head and h is the delivery head (see Fig. 1). The sensitivity analysis with respect to the mass of the impulse valve (e.g. the force Fv in present case where the valve is loaded by a dead weight) was made. The simulation results of the analysis are shown in Fig. 10 as dependencies of the hydraulic ram pump characteristics Q, q, Tc and h vs. force Fv. These results are very useful in a real case where the HRP system parameters should be adequately adjusted. Generally speaking, lower Fv values result in lower flow rates (Q and q) and a lower cycle period (Tc), whereas the efficiency (h) is higher. In the case when the HRP system has to deliver more water, the adjustment must be made for higher forces Fv, but then the efficiency is lower. It must be pointed out that there is only a limited range of Fv values for each particular HRP system. gfin^OtJJIMISCSD 03-3 stran 147 \ ^BSSITIMIGC 3 2 1 0 8.6 8.8 9 9.2 9.4 9.6 9.8 Filipan ^" _ Q, kg/c. --------A------ q, g/c. /• y Tc, s \^ / — -0 — Eta, - ./' / - // / - ^ X^ / - ~~- -L i ^ a'' y a - s' ®>^r y' - / 5 y^ ~~ e - ^c }— /' _ M __ w^^"^ ,--'"" ~~C x^ _ d , _-i ~--¦*""' ^ ^ . i i i 01234 5 Fv,N Sl. 10. Vpliv F na karakteristike hidravličnega ovna Q, q, T in h Fig. 10. The influence of Fv on hydraulic ram pump characteristics Q, q Tc and h 5 SKLEPI V prispevku podajamo matematični model hidravličnega ovna (HO), vgrajenega v cevni sistem. Obravnavamo poenostavljeni delovni postopek. Pretok v sistemu z ovnom je neustaljen, zato podajamo enačbe neustaljenega toka v ceveh in metodo karakteristik (MK). Podan je podrobni oris matematičnega modeliranja elementov sistema HO. Pogonsko in dobavno cev modeliramo s pravokotno mrežo MK, pogonski in dobavni zbiralnik pa sta robna pogoja s stalno gladino vode. HO modeliramo kot robni pogoj, ki popisuje delovanje ovna v področjih delovnega postopka. Ta robni pogoj je popisan z 11 enačbami in odvisnimi spremenljivkami, rešujemo ga s pomočjo metode iteracij v časovnih korakih računskega postopka. Postavljen model je programiran na računalniku. V podani računski simulaciji smo uporabili podatke iz literature, časovni potek tlaka in pretočne hitrosti so prikazani v grafični obliki. Postavljeni model omogoča bolj natančno analizo delovanja hidravličnega ovna in izluščitev zapletenih hidravličnih interakcij med elementi, ki so vgrajeni v ovnu. V občutljivostni analizi smo spremljali vpliv sile, delujoče na vzbujevalni ventil. 5 CONCLUSIONS Mathematical modelling of the HRP system is presented and the physics of a simplified working cycle of its action is explained. The flow in the HRP system is unsteady, therefore, unsteady-pipe-flow equations and the method of characteristics (MOC) are presented. The mathematical modelling of particular components of the HRP system is explained in detail. The drive and delivery pipes are modelled by the fixed grid MOC, and the supply and delivery tanks are the boundary conditions – open reservoirs with constant water levels. The HRP is modelled as a unit boundary condition describing the physics of its action. This boundary condition consists of 11 non-linear algebraic equations with 11 unknowns and is solved by an iterative procedure for each time step of the computational run. The derived model is programmed for a digital computer. A computer simulation using a pump system from the literature was made and the results are presented in the form of pressure and velocity vs. time graphs. The presented model enables a more detailed analysis of the operational principle of the HRP and gives insight into the complex hydraulic interactions occurring within the HRP. The sensitivity analysis presents the influence of the force acting on the impulse valve. 6 LITERATURA 6 REFERENCES [1] Joukowsky, N. (1900) Über den hydraulischen Stoss in Wasserleitungsröhren. Memoires de l’Academie Imperiale des Sciences de St.-Petersbourg, Series 8, Classe Physico-Mathematique, Vol. 9, No. 5. [2] O’Brien, M. P, J. E. Gosline (1933) The hydraulic ram. University of California Publication in Engineering, Vol. 3, No. 1, 1-58. [3] Lansford, W. M., W. G. Dugan (1941) An analytical and experimental study of the hydraulic ram. University of Illinois Bulletin, Vol 38, No. 22, 1-70. [4] Krol, J. (1951) The automatic hydraulic ram. Proc. Inst. Mech. Eng., Vol. 165, 53-65. VBgfFMK stran 148 Filipan V., Virag Z., Bergant A.: Matemati~no modeliranje - Mathematical Modelling [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] Eytelwein, JA. (1805) Bemerkungen über die Wirkung und vorteilhafte Anwendung des Stoßhebers (Belier hydrauliques). Realschulbuchhandlung, Berlin. Calvert, N. G. (1957) The hydraulic ram. The Engineer, Vol. 203, April 19, 597-600. Iversen, H. W. (1975) An analysis of the hydraulic ram. Trans. ASME, Journal of Fluids Engineering, Vol. 97, No. 2, 191-196. Basfeld, M., E.-A. Müller (1984) The hydraulic ram. Forsch. Ing.-Wes., Band 50, Nr. 5, 141-147. Chiang, Y.-C. (1984) Simulation and optimisation of transient oscillations, flow, and sound in complex piping system. Dissertation, University of Wisconsin, Madison. Rennie, L. C, E. A., Bunt (1981)The automatic hydraulic ram. The South African Mech. Eng, Vol. 31, No. 10, 258-273, No. 11, 286-311. Rennie, L. C, E. A. Bunt (1990) The automatic hydraulic ram - experimental results. Proc. Inst. Mech. Eng., Vol. 204, Part A: Journal of Power and Energy, 23-31. Thomas, T H. (1994) Algebraic modelling of the behaviour of hydraulic ram pump, Working Paper WP41, University of Warwick, Coventry. Tacke, J.H.P.M. (1988) Hydraulic rams - a comparative investigation. Communications on Hydraulic and Geotechnical Engineering, Delft University of Technology, Faculty of Civil Engineering, Report 88-1. Najm, H.N., PH. Azoury, M. Piasecki (1999) Hydraulic ram analysis: a new look at an old problem. Proc. Inst. Mech. Eng, Part A: Journal of Power and Energy, Vol. 213, Nr. 2, 127-141. Filipan, V (1998) Hydrodynamic modelling of power hydraulic ram. Dissertation, University of Zagreb, Zagreb, (in Croatian). Fox, J. A. (1989) Transient flow in pipes, open channels and sewers. Ellis Horwood Ltd., Chichester. Wylie, E. B., V.L. Streeter (1993) Fluid transients in systems. Prentice Hall, Englewood Cliffs. Press, W H, S.A. Teukolsky, W.T. Vettterling, B.P Flannery (1994) Numerical recipes in Fortran, The art of scientific computing. 2nd Ed., Cambridge University Press, Cambridge. Naslovi avtorjev: dr. Veljko Filipan Fakulteta za kemijsko inženirstvo in tehnologijo Univerza v Zagrebu Savska 16 HR-10000 Zagreb, Hrvaška veljko.filipan@fkit.hr Authors’ Addresses: Dr. Veljko Filipan Faculty of Chemical Eng. and Technology University of Zagreb Savska 16 HR-10000 Zagreb, Croatia veljko.filipan@fkit.hr prof dr. Zdravko Virag Fakulteta za strojništvo in ladjedelništvo Univerza v Zagrebu Ivana Lučiča 5 HR-10000 Zagreb, Hrvaška zdravko.virag@fsb.hr Prof Dr. Zdravko Virag Faculty of Mechanical Eng. and Naval Architecture University of Zagreb Ivana Lučiča 5 HR-10000 Zagreb, Croatia zdravko.virag@fsb.hr dr. Anton Bergant Litostroj E.I. d.o.o. Litostrojska 40 SI-1000 Ljubljana anton.bergant@litostroj-ei.si Dr. Anton Bergant Litostroj E.I. Ltd. Litostrojska 40 SI-1000 Ljubljana, Slovenia anton.bergant@litostroj-ei.si Prejeto: Received: 14.1.2003 Sprejeto: Accepted: 29.5.2003 Odprt za diskusijo: 1 leto Open for discussion: 1 year