Budapest University of Technology and Economics
Department of Hydrodynamic Systems
Dynamic Behaviour of Hydraulic Drives by Csaba H˝os, MSc
Supervisor: László Kullmann, MSc, PhD
Submitted in partial fulfilment of the requirements for the degree of Doctor of Philosophy in Mechanical Engineering in the Pattantyús Ábrahám Géza School of PhD Studies of the Budapest University of Technology and Economics 2005
c
Csaba H˝os, 2005
Nyilatkozat Alulírott H˝os Csaba kijelentem, hogy ezt a doktori értekezést magam készítettem és abban csak a megadott forrásokat használtam fel. Minden olyan részt, amelyet szó szerint vagy azonos tartalomban, de átfogalmazva más forrásból átvettem, egyértelm˝uen, a forrás megadásával megjelöltem.
Budapest, 2005. július 13.
A dolgozat bírálatai és a védésr˝ol készült jegyz˝okönyv a kés˝obbiekben a Budapesti M˝uszaki és Gazdaságtudományi Egyetem Gépészmérnöki karának Dékáni Hivatalában érhet o˝ k el.
iii
Dynamic Behaviour of Hydraulic Drives PhD Thesis by Csaba H˝os, MSc supervisor: László Kullmann, MSc, PhD Abstract This thesis studies dynamic phenomenon in hydraulic drives, which can occur either by design (e.g. sudden load change) or due to instabilities (e.g. valve chatter). The pressure transients in a linear actuator circuit are studied experimentally. The measured values are compared to those ones obtained by numerical simulation. For 1D simulation purposes novel mathematical models for transmission lines and for pilot-operated pressure relief valves are developed. By using the techniques of nonlinear dynamics the short, chamfered type valve seat is compared with the long, conical valve seat in terms of dynamic behaviour of poppet valves. The third part deals with the stability issues of the PD position controlling of a linear actuator circuit and it is shown how the control parameters and the digital sampling interval influence the accuracy and stability of the positioning.
Hidraulikus hajtások dinamikus viselkedése Doktori disszertáció készítette: H˝os Csaba, okleveles gépészmérnök konzulens: Dr. Kullmann László, egyetemi docens Kivonat A dolgozatban hidraulikus rendszerek m˝uködése közben fellép o˝ tranziens jelenségeket vizsgálunk, amelyek lehetnek üzemszer˝uek (pl. munkahenger terhelésváltás) vagy káros rezgések, instabilitások (pl. szelep rezgések). Egy hidraulikus emelo˝ berendezésen vizsgáljuk kísérleti úton a rendszer egyes pontjaiban a tranziensek hatására kialakuló nyomásváltozásokat. A mérések eredményeit összevetjük numerikus szimulációval kapott értékekkel. A szimulációhoz továbbfejlesztjük a hidraulika csövek és az el o˝ vezérelt nyomáshatároló szelep dinamikus modelljét. A dolgozat második részében a kúpos szelepek stabilitási kérdéseivel foglalkozunk és nemlineáris dinamikai módszerekkel összehasonlítjuk a rövid, letört ülék˝u és a hosszú, kúpos ülék˝u szelepek dinamikus viselkedését. Végül, egy hidraulikus munkahenger pozicionálási feladatán keresztül bemutatjuk a PD szabályozás paramétereinek és a digitális mintavételezés hatását a szabályozás stabilitására.
Köszönetnyilvánítás Ellentétben a tudományos publikációkkal, a doktori disszertáció címlapján csak egy szerz o˝ neve szerepelhet, noha a dolgozat számos ember közös munkájának eredményeit tartalmazza. Így - más lehet o˝ ség hiányában - kísérletet teszek felsorolni mindazokat, akik id˝ot és fáradtságot nem kímélve segítették e dolgozat létrejöttét. Mindenekel˝ott köszönöm feleségemnek, Orsolyának, hogy házasságunk els o˝ éveiben elnézte, hogy szabadido˝ met nem vele töltöm, hanem gyakran hétvégén is a jegyzeteimet böngészem vagy a számítógép képerny o˝ jét bámulom. Hálás vagyok szüleimnek azért a biztos családi háttérért, ami leheto˝ vé tette, hogy a tanulásra, majd késo˝ bb a tudományra koncentráljak. Köszönöm témavezet˝omnek, Kullmann Lászlónak, hogy a megfelelo˝ pillanatban mindig átlendített a holtponton, egy-egy jó ötlettel vagy észrevétellel. To˝ le a puszta szakmánál sokkal többet tanultam; precíz és alapos gondolkodást, ˝ vezetett be a keleti sz˝onyegek szövési és festési technikáinak rejtelmeibe is. Szintudományos becsületességet - és O tén mély hálával tartozom szobatársamnak, Pandula Zoltánnak, aki türelmesen segített megkeresni a programhibákat és ha kellett, szó nélkül újratelepítette az operációs rendszert a számítógépemre (miközben türelmesen hallgatta a “de-éncsak” kezdet˝u magyarázkodásomat). A laboratóriumi mérések során szintén felbecsülhetetlen segítséget kaptam T o˝ le. Ezúton szeretném megköszönni Halász Gábor tanszékvezet o˝ úrnak a lehet˝oséget, hogy valós mérnöki munkákban is részt vehettem és Bálint Lajosnak, hogy hatalmas gyakorlati tapasztalatával a labormunkák során annyit segített. Hálás vagyok Vajna Zoltán akadémikus úrnak és Paál Györgynek a dolgozat alapos átolvasásáért és hogy gondolkodásmódban annyit tanulhattam T˝olük. A Hidrodinamikai Rendszerek Tanszék (volt Vízgépek Tanszék) minden dolgozójának köszönöm azt az oldott családias légkört, amiben öröm volt dolgozni. Külön köszönet illeti Stépán Gábor akadémikus urat, aki megtisztelt közvetlenségével és bizalmával. Lelkesedése mindig átragadt rám is és ez a kritikus pillanatokban sokat segített. Neki köszönhetem érdekl o˝ désemet a nemlineáris dinamika iránt és els˝o hosszabb külföldi tanulmányutamat is, ami során felbecsülhetetlen tapasztalatokkal gazdagodtam. Ezúton köszönöm Bende Margit tanárno˝ nek azt a bizalmat, amivel másodéves hallgató koromban - gyenge teljesítményem ellenére - megajándékozott. Nélküle nem lett volna bátorságom nekivágni doktori tanulmányaimnak. Hálás vagyok Stépán Gábornak és Gáspár Tibornak, hogy méréseket végezhettem a M˝uszaki Mechanikai Tanszék dinamikus rezgésmér˝o m˝uszerével. Egyetemi tanulmányaim alatt nagy hatással voltak rám az Áramlástan Tanszéken hallgatott el o˝ adásaim, ezeket szeretném megköszönni Lajos Tamásnak, Kristóf Gergelynek, Vad Jánosnak és Koscsó Gábornak.
Acknowledgements Part of this work was performed while the author worked at the University of Bristol, Department of Engineering Mathematics as a research fellow under the supervision of Alan Champneys. I gratefully acknowledge his support and patience while guiding me through the field of nonlinear dynamics and 72-column-long f77 codes. I enjoyed working with him and all the people at the Nonlinear Dynamics Group in Bristol very much.
v
Summary Although the electronic drive technology superseded hydraulic drives in several industrial applications in the last decade, industrial hydraulics is still popular and wide-spread in many cases such as e.g. machining tools, fatigue-testing machines, airplane industry or mining. Hydraulic drives exhibit many advantageous properties, e.g. they are simple, large force (torque) can be exerted, the velocity of the drive can be set simply and continuously, the hydraulic fluid has excellent lubricant properties and overload protection is simple. The role of numerical simulation increases nowadays in the design of hydraulic elements and systems. On the one hand, by calculating and analysing the full three-dimensional flow field the characteristics of system elements can significantly be improved. On the other hand, with the help of one-dimensional transmission line models and lumped parameter element models, the transient response of an entire hydraulic system can be predicted. In this thesis, fast and accurate techniques are developed for one-dimensional hydraulic transmission line modelling and for calculating the unsteady friction factor. A dynamic model for the pilot-operated pressure relief valve is presented, in which the fluid forces are calculated by simple, closed formulae fine-tuned by Computational Fluid Dynamics (CFD) simulations. The capabilities of these models are verified with the help of system simulation; the pressure transients induced by sudden load change and piston impingement against the cylinder top are studied both experimentally and numerically. The second part of the thesis is mainly theoretical and unfolds the interesting phenomenon emerged during the modelling and the experimental work. While studying the pilot-operated pressure relief valve, it was found that the poppet valve body in the pilot stage is poorly damped and during relief valve operation, characteristic peaks appear in the acceleration spectra measured on the pilot housing. These undesired vibrations not only increase the wear of the valve elements but also increase the noise emission. To overcome these vibration problems, long conical valve seats are suggested. On one hand, the performance of the valve is easier to predict in such cases. On the other hand, the long, conical valve seat is more advantageous from the dynamical behaviour point of view. This latter statement is proved with stability analysis, in which it is shown that the long, conical valve seat exhibit an order of magnitude larger relative damping compared to the short, chamfered type valve seat. The third part of the thesis deals with the stability issues of the PID control of a linear actuator circuit. The behaviour of the control algorithm is studied in the case of two different control valves, namely two-way control valve and proportional valve. It is well known from control theory that if the proportional gain P increases so does the accuracy of the position control. However, beyond a critical value of the proportional gain self-excited oscillations appear. Theoretically, increasing the differential gain D may stabilise the system in such cases but our experiments show that the increase in the differential gain may have an opposite effect. In this part, the main challenge is to perform stability analysis on a system, whose governing equations are not continuous or not continuously differentiable. It is shown that by employing the geometric singular perturbation theory and by separating the ’fast’ and ’slow’ variables, one can draw a qualitative picture of how the slight compressibility of the hydraulic fluid leads to self-excited oscillations. By ’smoothing’ the discontinuities, closed-form relationships are given for the critical parameter values corresponding to the onset of oscillations. Finally, it is shown how the delay due to digital control influences the stability of the control. The theoretical results are verified by experimental ones.
vi
Összefoglalás Bár az utóbbi évtizedben az elektromos hajtástechnika számos területen kiszorította az ipari gyakorlatból a hidraulikus hajtástechnikát, napjainkban még mindig igen népszer˝uek és széles körben elterjedtek a hidraulikus hajtások, pl. szerszámgépek, fárasztógépek, repülo˝ gépek vagy a bányászat területén. A hidraulikus hajtások számos el o˝ nyös tulajdonsággal bírnak, így például egyszer˝uek, igen nagy ero˝ t vagy nyomatékot képesek kifejteni nagyfokú merevség mellett, a hajtás sebessége (vagy a szögsebessége) egyszer˝uen és fokozatmentesen állítható, a munkavégz o˝ közeg egyben kiváló kenési tulajdonságokkal is rendelkezik, egyszer˝u a túlterhelés elleni védelem stb. Napjainkban a numerikus szimuláció egyre nagyobb szerephez jut a hidraulikus rendszerek és rendszerelemek tervezése során. Egyrészr˝ol a teljes háromdimenziós áramlási tér kiszámításával jelento˝ sen javíthatók az egyes elemek jellemz˝oi, másrészr˝ol egydimenziós cs˝omodellek és koncentrált paraméter˝u elem modellek segítségével teljes hidraulikus hajtások tranziens átállási folyamatai számíthatók. A dolgozat els o˝ részében egydimenziós tranziens szimulációhoz kifejlesztett gyors és pontos módszereket mutatunk be a hidraulika csövekben lejátszódó hullámjelenségek és az instacionárius cs˝osúrlódási tényez˝o számítására. Bemutatjuk a kísérletekben alkalmazott elo˝ vezérelt nyomáshatároló szelep matematikai modelljét, melyben az áramlástani ero˝ ket olyan egyszer˝u, zárt alakú összefüggésekkel számítjuk ki, melyekben egyes paramétereket (pl. átfolyási tényezo˝ ) háromdimenziós áramlástani szoftverrel (CFD, Computational Fluid Dynamics) határoztuk meg. A matematikai modellek alkalmasságát rendszerszimulációval bizonyítjuk; a kísérleti berendezésen a hirtelen terhelésváltás, ill. munkahenger felütközés hatására kialakuló nyomáslengéseket vizsgáljuk mérésekkel és numerikus számítással. A dolgozat második része alapveto˝ en elméleti jelleg˝u és az els˝o részben leírt modellezési, ill. kísérleti munka közben felmerült kérdésekben mélyül el. A nyomáshatároló szelep modellezése közben azt találtuk, hogy az el o˝ vezérl˝o fokozat kúpos szeleptestje igen gyengén csillapított. Valóban, a szelep m˝uködése során jellemz o˝ frekvenciájú rezgéseket tapasztaltunk a szelepházon. Ez nemcsak a szeleptest és a szelepülék fokozott igénybevétele miatt káros, hanem megnöveli a berendezés által kibocsátott zajszintet is. Hosszú, kúpos ülék alkalmazását javasoljuk, így egyrészr o˝ l pontosabban és megbízhatóbban méretezheto˝ a szelep (a méretezési összefüggéseket megadjuk a dolgozatban), másrészr o˝ l a hosszú, kúpos ülék dinamikailag is kedvezo˝ bb. Ez utóbbi állítást stabilitásvizsgálattal támasztjuk alá, melyben kimutatjuk, hogy a hosszú, kúpos ülék esetén már nagyságrendekkel kisebb csillapítási tényez o˝ vel is elérhet˝o a stabil üzemállapot a rövid, letört ülékhez képest. A dolgozat harmadik része a kísérleti hidraulikus emelo˝ berendezés szabályozásának stabilitási kérdéseivel foglalkozik. Az iparban leggyakrabban elterjedt P(I)D szabályozási algoritmust vizsgáljuk különböz o˝ munkahenger szabályzó szelepek esetén. A szabályozáselméletben közismert, hogy az arányos tag (P) növelése javítja a beállási pontosságot, de egy kritikus érték felett destabilizálja a rendszert. Ilyenkor a csillapítás (D tag) növelése javíthat ugyan a helyzeten, de gyakran éppen ellentétes hatást vált ki a differenciáló tag növelése. Ebben a részben a leíró differenciálegyenletek folytonosságának és folytonos differenciálhatóságának hiánya jelenti a kihívást; ez a fejezet alapkutatás jelleg˝u. Bemutatjuk a szinguláris geometriai perturbációelmélet alkalmazását hidraulikai rendszerekre, melynek segítségével a “gyors” és “lassú” mozgások szétválasztásán keresztül min o˝ ségileg megérthet˝o, hogy a hidraulika folyadék kismérték˝u összenyomhatósága hogyan vezet az öngerjesztett lengésekhez és a teljes rendszer m˝uködésképtelenségéhez. A szakadások ’kisimításával’ zárt alakú összefüggéseket adunk a kritikus arányos- és differenciáló tag értékekre. Végül bemutatjuk, hogy a digitális mintavételezés okozta ido˝ késés hogyan hat a szabályzás stabilitására és az elméleti számításokat mérésekkel támasztjuk alá.
vii
Contents 1 Introduction
1
1.1
Background - History and future trends . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.2
Hydraulic research at the Department of Hydrodynamic Systems . . . . . . . . . . . . . . . . . . . . . .
3
1.3
Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2 Experimental set-up and mathematical modelling of the system elements 2.1
Introduction and outline of the chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.2
Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2.1
Elements of the circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.2.2
Data acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2.3
Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.3
Hydraulic cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.4
Pilot-operated pressure relief valve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
2.4.1
Pressures and flow rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
2.4.2
Fluid forces on the valve bodies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
2.4.3
Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
2.4.4
Static performance curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
Transmission line modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
2.5.1
Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
2.5.2
Friction modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
2.5.3
Numerical techniques for transmission line modelling . . . . . . . . . . . . . . . . . . . . . . .
31
2.5.4
Comparison and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
2.5.5
Handling pipe junctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
2.6
Numerical simulation of hydraulic transients in the experimental rig . . . . . . . . . . . . . . . . . . . .
42
2.7
Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
46
2.5
3
6
Analysis of the flow around a poppet valve in the case of long conical seat
49
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
3.1.1
Outline of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
3.2
The laminar flow through a pilot valve with long conical seat . . . . . . . . . . . . . . . . . . . . . . . .
50
3.3
Pressure force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
3.4
Viscous force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
3.4.1
Estimation 1. - Equivalent cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
3.4.2
Estimation 2. - Poiseuille flow between conical walls . . . . . . . . . . . . . . . . . . . . . . . .
57
viii
3.5
3.4.3
Estimation 3. - 2D flow in a diverging channel . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
3.4.4
Comparison of the models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
4 An analysis of valve chatter 4.1 4.2
4.3
4.4
4.5
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
4.1.1
Outline of this chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
Poppet valve with chamfered-type seat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
4.2.1
Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
4.2.2
Linear analysis with no spring precompression . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
4.2.3
Two-dimensional centre manifold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
68
Poppet valve with long conical seat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
4.3.1
Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
4.3.2
Linear analysis with no spring precompression . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
4.3.3
Two-dimensional centre manifold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
Further studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
4.4.1
Comparison in terms of stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
4.4.2
Effect of spring precompression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
4.4.3
Impact dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
Measurements and summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
5 Stability of controlled hydraulic systems 5.1
64
80
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
5.1.1
Outline of the chapter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
81
5.1.2
Some preliminary remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
5.2
The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
5.3
Relay control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
85
5.3.1
Application of singular perturbation theory for relay-like control . . . . . . . . . . . . . . . . . .
86
5.3.2
Asymptotically smooth analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
87
5.3.3
Piecewise continuation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
88
Proportional control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
5.4.1
Theoretical analysis of analogue control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
5.4.2
Numerical results for analogue control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
91
5.4.3
Effect of digital control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
92
5.4
5.5
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6 Papers presenting the results of this thesis
101
7 Bibliography
102
“The first principle is that you must not fool yourself - and you are the easiest person to fool. So you have to be very careful about that. After you’ve not fooled yourself, it’s easy not to fool other scientists. You just have to be honest in a conventional way after that. ... One example of the principle is this: If you’ve made up your mind to test a theory, or you want to explain some idea, you should always decide to publish it whichever comes out. If we only publish results of a certain kind, we can make arguments look good. We must publish both kinds of results.” Richard P. Feynman
x
Chapter 1
Introduction 1.1 Background - History and future trends The history of hydraulic power starts surprisingly long ago; a water mill - the first hydraulic motor - is depicted on a mosaic at the Great Palace in Byzantium, dating from the early fifth century. The first record of a water mill goes back even further, to around 100 BC, and the origins may indeed have been much earlier [52]. One of the earliest applications which combined both the high power capability and the fine precision of hydraulics was a system to elevate guns of the battleship Virginia in 1906. Later, the development of small, light and efficient power units was accelerated by the use of hydraulic power in aircraft applications; the first commercial application was on the Boeing 247D twin engine transport in the early 1930s. Nowadays, hydraulic systems have many applications in modern industry ranging from automotive industry through heavy industry to aerospace applications. The main advantages of industrial hydraulics are • exertion of large forces (torques) with small dimensions (large power density), • the force exerted automatically adapts to the load, • the motion can be started with full load (no need for gradual load increase), • good controllability, i.e. velocity, force, torque, etc. can be easily and continuously set, • simple overload prevention, • ability to realize both quick and extremely slow motions and • advantageous labour safety properties, e.g. no danger of electric shock. However, some of the unbeatable characteristics of former days, e.g. excellent dynamics and ease of control nowadays can be achieved with the state of art control techniques. Moreover, bad energetic efficiency of the resistance controlled drives caused already some hydraulic drive application to fall in hand of electric drives (for details, see [68]). To keep the market of hydraulic drives for the future, Backe in [10] suggests the following potential trends of fluid power technology. 1. The increased use of scientific methods of optimisation of stiffness and deformation characteristics of metallic and non-metallic materials, for the optimisation of flow within components and the minimisation of noise emission. 2. The increased use of simulation programs to reduce layout, testing-time and price of components and systems. 3. The use of novel materials and designs in order to improve the wear characteristic and tribological properties, in order to reduce construction effort and price. 1
1. Introduction
4. Protecting the environment through: (a) the use of energy efficient components and systems, (b) further developments in sealing and connection technology and (c) development and use of biodegradable pressure media. 5. Improving the static and dynamic characteristics of fluid power drives through: (a) the use of digital electronic signal circuits, (b) the integration of device-specific electronics into the components and (c) increased use of modern controlling concepts. Let us explain a few of the above issues in detail concerning this thesis. There is an impressive development in both the speed and storage capacity of computers 1. The consequences in the field of engineering are threefold. First, such detailed simulations are possible nowadays that would have taken unacceptably long before. On the one hand, the flow inside basic elements (valves, pumps, etc.) with complex geometry can be computed with CFD (Computational Fluid Dynamics) techniques, giving valuable information such as detailed pressure and velocity field, separation/cavitation locations or even complete performance curves can be predicted. On the other hand, systems of such elements can be built up in a suitable computer simulation package allowing the engineer to answer typical ’what if’-type questions - taking into account the dynamic characteristics of all systems elements simultaneously. Both techniques reduce the layout and the testing time - thus the price. Secondly, fast computers allow us to gather an enormous (sometimes overwhelming) amount of data with experimental techniques like PIV or LDA. These techniques help us understand the flow structure and give a direct comparison to e.g. CFD simulations. Thirdly, high CPU speed makes it possible to adopt complex control laws which need a massive amount of calculations - with the same (or even quicker) sampling frequency. Beside cost, efficiency and other traditional criteria, noise emission became more and more important in the last decade when rating hydraulic systems and/or elements. Typical noise sources are pumps, motors and valves. Pumps and motors used in fluid power systems are invariably of the positive displacement type. These pumps and motors do not produce a steady flow, but a pressure pulsation, usually called flow ripple. This ripple originates from the compressibility of the oil; as each ’pumped volume’ gets connected to the high-pressure discharge side of the pump, there will be a momentary backflow to equalise the pressure. This gives rise to large pressure peaks and this mechanism is responsible for noise generation in pumps (similar phenomena occur in the case of motors). Valve noise can arise from instability, turbulence and cavitation. Instability generates considerable noise (sometimes referred to as ’valve chatter’) and is due to the problems of valve design and/or dynamic interaction between moving valve elements and alternating pressure level. However, there is a slight confusion of the terms ’instability’ and ’resonance’. The main difference between the two phenomena is that resonance is caused by external excitation of an element close to its eigenfrequency. Should the excitation disappear, the heavy vibration vanishes also. However, in the case of instability, there is no ’need’ for external excitation (e.g. from the pump), thus this phenomena is also called self-excited vibration. Turbulence noise 2 is almost certainly present in most control elements (valves) but it is usually negligible in relative terms. The most common noise source in the case of valves is cavitation, in which case bubbles are formed as the pressure drops while the fluid accelerates through narrow gaps. As the bubbles collapse, short duration but intense pressure peaks are generated which results in a high frequency noise. When considering hydraulic systems, their control has received considerable attention because they provide highperformance motion control. Typical tasks are position, velocity and force control. However, the dynamic characteristics of hydraulic elements are complex, highly nonlinear and often even discontinuous, which makes the control design a 1 Moore’s
law - doubling of data density and processing speed in every approx. 18 months - still holds:
http://www.intel.com/research/silicon/mooreslaw.htm. 2 Noise emission can be present without turbulence, e.g. due to periodic vortex formation. However, the term ’turbulent noise’ is used in the hydraulic community for such noise sources.
2
1.2. Hydraulic research at the Department of Hydrodynamic Systems
challenging task. The main nonlinearities include the compressibility of the fluid, complex flow properties of the valves and friction forces arising in hydraulic actuators. Moreover, the operating conditions and the unpredictable disturbances may vary in a complicated fashion and usually parameters are also not known precisely, e.g. temperature-dependent viscosity. An ideal controller would be robust against disturbances and uncertainties - and give simultaneously the best performance. Beside traditional linear control techniques (mostly PID control), several new control strategies were tested in the last decade, such as sliding mode control, fuzzy control, neuro-control or adaptive (self-tuning) control. An excellent overview of these techniques is given in [40]. However, 90% of the industrial controllers is still fixed-gain PID control, thus it is essential to better understand and improve this technique. The main problem with theoretical development is the lack of mathematical theory for non-differentiable or discontinuous differential equations. Consider for example a simple piston controlled by (say) a P controller via a critically centred servovalve. This system is not linearisable at the equilibrium state because of the different pressure-terms in the two spool directions; in one direction the tank pump pressure is connected to one of the piston chambers, in the other direction the tank pressure is connected to the same chamber. Moreover, if the control valve happens to be an on-off valve, the system is not even continuous at he equilibrium state3 . It is very interesting that recently, on-off type valves are becoming popular again as if several on-off valves are connected in parallel with different sizes (Q N , 2QN , 4QN , etc.), flow rate can be controlled precisely with switching on only the appropriate ones. This technique is sometimes called ’digital hydraulic control’ and its main advantage is that (say) 5 on-off valves may still be cheaper than one servovalve, not to mention the clear advantages in contamination sensitivity and maintenance costs. Finally, let us mention that the presence of time delay due to the sampling frequency of the control gives rise to completely new problems; delay equations may have surprising behaviour, e.g. such parameter variations (notably increasing damping), which stabilise a continuous system, may destabilise a delay system. As it can be seen, hydraulic systems exhibit a surprisingly wide range of interesting engineering challenges both practical and theoretical ones. Due to the strong competition between hydraulic and electronic drives, a profound knowledge of solid and fluid mechanics, control systems, nonlinear dynamics, up-to-date mathematical techniques and practical skills are needed for designers to be able to build suitable systems for various industrial areas and maintain the market for hydraulics.
1.2 Hydraulic research at the Department of Hydrodynamic Systems The predecessor of the Department of Hydrodynamic Systems, Department for Hydraulics and Hydraulic Machines was founded in 1899 by Donát Bánki at the Budapest University of Technology and Economics,whose name was Royal Joseph Technical University at that time. The first recorded study on hydraulic transients traces back to 1956, in which the protection of a pump-pipeline system was investigated from overpressures generated by pump start and shut-down. During the decades, the analysis of hydraulic transients has remained one of the most intensely researched area at the Department, which was (and is still) motivated by the frequent industrial applications. The challenge of such computations lay in the need for simultaneous solution of partial differential equations (describing the wave propagation in the long transmission lines), ordinary differential equations (unsteady pump models, check valves, etc.) and algebraic equations (e.g. continuity equation in the pipe junctions). The advantages of digital computer simulation were early recognised; the senior researchers of the Department are still telling horror tales about the ’debugging’ of punch cards used in transient calculations. Thanks to the recent advances in computer sciences (hardware and software), nowadays a calculation including a few tens of pumps and a few hundred (long) transmission lines is also possible and finished within a few days (and the most of this time is the data pre- and post-processing). 3 Moreover,
in such cases the equilibrium does not even exist in the traditional sense.
3
1. Introduction
The first recorded oil hydraulic research, which was the stability analysis of a servomechanism was performed by Kullmann in 1971, [48]. An experimental study of a hydraulic pressure relief valve was conducted in 1976. This work was followed by several, mainly experimental studies on relief valves and on the (unsteady) loadability of hydraulic hoses. In his PhD thesis [20], Bálint applied the technique for pump system transient simulation for the fully nonlinear numerical simulation of a linear actuator circuit. With this technique, it was possible to predict the whole starting (or stopping) process contrary to the previously applied techniques, which were mainly linearising the system and then solving it analytically. The work focused mainly on the effect of the length of the transmission lines and its modelling issues and concluded that is it essential to take into account the compressibility of the hydraulic fluid and the elasticity of the transmission line in order to predict the pressure histories accurately.
1.3 Thesis overview The main objectives of this work were 1. to build up a computer simulation of an existing hydraulic test rig, perform measurements and computer simulations and compare the results, 2. based on the experiments, identify the discrepancies between the computer model and the experimental results and develop novel ways of modelling with improved performance and 3. focus on the nonlinear phenomena present in the system, understand the underlying mechanisms and find techniques to eliminate them if they worsen the system’s performance. Each problem dealt with in this thesis emerged while working with the experimental apparatus described in detail in §2.2. This simple hydraulic linear actuator circuit consisting of only four elements (gear pump, hydraulic cylinder, cylinder control valve, pressure relief valve) plus pipes yielded a surprisingly wide range of interesting phenomena - from theoretical mathematics and hydrodynamics through nonlinear dynamics to measuring techniques. The first part of this thesis, §2 describes the elements of the computer simulation. The two issues of the system, where the corresponding literature was found to be insufficient was the transmission line modelling and the dynamic modelling of the pilot operated pressure relief valve. In the case of the transmission line modelling, the following two problems emerged. 1. First, the traditional way of distributed parameter modelling is the method of characteristics, which requires a fixed time step. There is obviously the possibility of interpolating between the time steps but in this case the accuracy and simplicity - the main strengths - of the technique are corrupted. Meanwhile, changing the integration step size in the case of the ordinary differential equations describing other elements (e.g. the piston or the relief valve) is essential for an efficient numerical code. §2.5.3 describes some numerical schemes, which allow the variation of integration time step while calculating the 1D flow in hydraulic transmission lines. These results are also presented in [4]. 2. Another issue concerning to transmission line modelling is the insufficient damping experienced when comparing measured and calculated hydraulic transients. There are several unsteady friction models available in the literature, of which the two most wide-spread are Zielke’s model and the Brunone-Vítkovsky model. However, the implementation of the first one is computationally expensive while the literature reports stability problems when using the second one with an explicit scheme. §2.5.2 describes two numerical techniques overcoming these difficulties and allowing their implementation in an easy and robust way. The corresponding publications are [4] and [6]. 4
1.3. Thesis overview
Concerning the dynamic modelling of the pilot operated pressure relief valve, little work was found about the dynamic issues of the pilot stage. Our aim was to establish equations for describing the fluid forces on the valve body as a function of inlet pressure and valve opening. The main challenge of the model development described in §2.4 was not the elaborate mathematical derivation of the equations (actually, it was mainly basic fluid mechanics and some algebra) but the parameter identification. The author’s experience was that it was hard to identify the various parameters and match them with experimental results. Typically, a slight change in one of the parameters (say, some inner orifice radius hard to measure manually or some damping coefficient) may cause strongly unrealistic computed performance (e.g. characteristic curves). The results are presented in [5]. For more than 30 years, valve chattering has been a problem for valve designers, which is still not properly understood and mainly empirical guidelines are given on how to avoid it. It is also known that poppet valves with chamferedtype seat have a stronger tendency to oscillate than those ones with long, conical seats. Interestingly, similar phenomena were found in the test rig under analysis, namely that there was a characteristic peak in the spectra of the acceleration signal measured on the relief valve housing, which occurred only when the relief valve worked. This phenomena initiated §4, where an extremely simple model is presented for valve chatter analysis. In this chapter, both chamfered-type seated and long, conical seated poppet valves are analysed. The study focuses on the nonlinear dynamics of these valves and calculates critical damping parameters for stable discharge. Although the models are strongly simplified and several physical phenomena are omitted, the results clearly show that poppet valves with long, conical seats are indeed dynamically more stable (i.e. smaller relative damping is needed for stable operation). The formulae needed for the valve chatter analysis in §4 are the relationships between the flow rate, pressure drop and fluid forces on a poppet valve with long, conical seat, where laminar flow develops. Surprisingly, little literature was found concerning this problem, notably the damping forces on such a moving valve body. §3 describes the calculation of flow rate and pressure distribution in the case of such a poppet valve configuration. Pressure forces and viscous forces are also calculated and compared to CFD results. Note that although the results show a limited agreement when compared to CFD simulations, the author of this thesis is not aware of any similar formulae. The last section, §5 deals with the positioning problem of a simple hydraulic linear actuator system with the help of PD control. The first part shows how to perform a stability analysis in the case of an on/off valve. The reason of analysing this highly unusual control is that recently there is an effort in hydraulics to substitute expensive and sensitive servovalves with cheap and robust on/off valves connected in parallel. §5.3 describes some mathematical techniques adequate to cope with the discontinuous nature of on/off valves and shows how they can contribute to control design for such systems. The second part of §5 describes what kind of possible instabilities can be encountered if one uses PD control together with a proportional valve for piston position control. Again, stability analysis is performed in §5.4 and the resulting analytical formulae show how the control parameters influence the system’s performance. Finally, in §5.4.3 an interesting phenomenon is studied. During the experimental verification of the results of §5.4.1, the controller parameters influenced the system’s performance contrariwise as expected; notably the damping parameter D destabilised the system although theory suggested that increasing D helps the controller to overcome oscillations. It was finally found that the delay due to digital control has a fundamental effect on the stability of the control. The maximal allowable proportional gain P as a function of the sampling (and controlling) interval is computed to avoid instability and the theoretical results are verified by experiments. The results of this section are also presented in [3], [2] and [1].
5
Chapter 2
Experimental set-up and mathematical modelling of the system elements 2.1 Introduction and outline of the chapter This chapter describes the elements of the experimental vertical linear actuator circuit, the mathematical modelling of its elements, the issues corresponding to the numerical solution (computer simulation) and finally presents and compares the results of the numerical solution and the experiments. §2.2 describes the experimental set-up, the data acquisition process and the identification of the fluid properties (density, viscosity and bulk modulus). The main emphasis is then on the mathematical modelling of the system elements, especially on the issues corresponding to the pilot-operated pressure relief valve in §2.4 and the transmission line modelling in §2.5. Based on the corresponding literature (notably [20]), the model of the hydraulic cylinder was straightforward and is given in §2.3. The only specific issue corresponding to cylinder modelling was the description of the friction force, which was experimentally determined and then the approach of Bisztray-Balku presented in [14] was employed for an analytical description. After the literature survey on the modelling of pilot operated relief valves, it was found that (a) due to the large variety of their construction, no general model can be developed (as e.g. for hydraulic cylinders), (b) little attention was devoted to the dynamic operation of poppet valves (notably unsteady thrust estimation and damping forces) and (c) the problem of ’valve chatter’ is often encountered but is not clearly understood. §2.4 describes the model development for the particular valve used in the experiments, with a special emphasis on the thrust estimation, which, in the case of the pilot stage, is verified with CFD calculations. The model is quasi-steady in terms of flow forces; the thrust on the valve bodies obtained under steady-state assumptions is assumed to hold under dynamic operation as well. The result of §2.4 is a set of ordinary differential equations, which are capable of describing the dynamic behaviour of the pressure relief valve in a realistic way. Traditionally, the method of characteristics (MOC) is used for 1D transmission line modelling together with steady-state friction terms (Darcy’s friction factor). However, due to the fixed time step of the MOC, it is cumbersome to achieve compatibility with the varying (and usually much smaller) step-size needed by the ordinary differential equation solvers used for integrating the equations describing e.g. the cylinder or valve behaviour. Developing variable step size techniques for transmission line modelling is thus essential for an efficient integration. §2.5.3 presents some numerical schemes, which allow the transmission line solver to vary the time steps, yet provide an accurate solution with small extra effort. To be specific, the 1D continuity and equation of motion are first discretized in space with (a) three-point and (b) 6
2.2. Experimental setup
five-point central difference scheme, (c) finite element method and (d) Chebyshev collocation technique. The resulting system of ordinary differential equations are then integrated using several time marching schemes; (a) implicit Euler scheme, (b) second-order backward differentiation formula and (c) fourth- and (d) fifth-order explicit first stage, singly diagonally implicit Runge-Kutta schemes. The pairs of schemes are compared in terms of accuracy and efficiency and it is shown that higher-order schemes are advantageous. Another - from the viewpoint of practical application - important question in the field of transmission line modelling is the estimation of unsteady pipe friction. It is well-known that the damping effects in the real pipeline systems are usually larger than that one predicted by simulation. This is due to the fact that steady-state friction modelling is used although the process itself (pressure wave propagation) may be highly unsteady1. Unsteady friction models are available in the literature but they are either computationally expensive to implement or suffer from (numerical) stability problems. §2.5.2 presents two techniques of implementing the two most wide-spread unsteady friction models in a quick and robust way. Finally, in §2.6 it is demonstrated that the previously described models are capable of predicting the dynamic behaviour of the experimental hydraulic circuit presented in §2.2. The mathematical models developed in the previous sections are built together in a suitable programming framework and then calculations are performed to predict the unsteady behaviour of the experimental rig. The test cases were sudden load change and piston impingement against the cylinder top; the measured and calculated pressure histories are compared in §2.6.
2.2 Experimental setup This section describes the vertical linear actuator circuit used in the measurements. The rig was built for the experimental verification of the dynamic calculations performed by Bálint in his PhD thesis, [20]. §2.2.1 describes the elements of the test rig, §2.2.2 presents the data acquisition process and finally §2.2.3 gives the material properties of the hydraulic fluid.
2.2.1 Elements of the circuit The hydraulic supply unit consists of a gear pump type Balcanar C54X, V g = 24[cm3 /rev.], pn = 175[bar] (’1’ in Fig.2.1) and a pilot operated pressure relief valve ’2’ type Danuvia DB10-2.10/315. The primary drive was a three-phase asynchronous motor, whose revolutionary speed could be varied between 360[rpm] and 1800[rpm]. The filter ’3’ ensures the required quality of the hydraulic fluid, which was especially important when the proportional valve was used for piston controlling (details are given later in this section). The two-position valve ’4’ was merely used for allowing other measurements (denoted by ’9’) with the same hydraulic supply and its position was not changed during the measurements. For piston controlling, two types of valves were used, both denoted by ’5’ in Fig.2.1. For the sudden load change and piston impingement experiments presented in §2.6, a 3/4 Danuvia 4WE 10G 4.2/G24N valve was used. This valve type does not allow a continuous setting of the flow rate; it simply switches between the three possible positions. However, for the feedback control experiments presented in §5.4, a proportional valve, type Bosch Rexroth 4WRAE 10 E302X/G24K31/A1V was used. In this valve, the spool body displacement and hence the flow-through area is directly proportional to the input signal (in this case voltage, −10[mV] . . . 10[mV]). Thus with this proportional valve, it is possible to vary the flow rate (and thus the velocity of the piston) continuously. However, it is to be noted here that the flow rate is not a linear function of the input signal; if such characteristic is needed, a servovalve should be used. The hydraulic cylinder ’6’ is a double-acting one manufactured by Újpesti Gépelemgyár (group 43), with cylinderto-rod diameter ratio D/d = 50/32 and a stroke of approx. 280[mm]. It is possible to change the load lifted from 0[kg] 1 However,
if the pipe is highly flexible - e.g. in the case of arterial blood flows - friction may play a negligible role compared to the ’damping’ effect
of the pipe diameter change. An example can be found in [49].
7
2. Experimental set-up and mathematical modelling of the system elements
Figure 2.1: Experimental linear actuator circuit, the pipe numbering is given by the length indices (see text for details).
to 210[kg] with the help of metallic cylinders ’7’, whose weight is approx. 40[kg] 2. Two types of pipeline were used for connecting the elements; one of them being metal pipe with outer diameter 15[mm] and wall thickness 2[mm] and the other one was a hydraulic hose (HIDOL), whose inner diameter is 11[mm].
2.2.2 Data acquisition Three pressure transducers and a position transducer were used for data acquisition. An inductive position transducer type Hottinger-Baldwin W100K recorded the displacement of the piston. The type of the pressure transducers was HottingerBaldwin P11. The locations of the pressure measuring taps are shown in Fig.2.1; the two chamber pressures of the hydraulic cylinder (measured where the pipes are connected to the cylinder) and the pressure immediately after the pressure relief valve were recorded. The signals - after amplification by HB MGC digital amplifier - were collected with National Instruments DAQ card. Although little change was experienced, every transducer was re-calibrated before each measurement. The data acquisition process was managed with the help of the LabView software environment and a National Instruments DAQ card. The sampling frequency of the card was 100[kS/s] (kiloSamples/sec), its resolution was 12 bits and the analogue input/output was 0 − 10[V]. As only one A/D converter is built into such cards, multiplexing is used between the channels. In the case of four input channels, the highest available sampling frequency was 50[kHz]. This also introduces a delay of 1/200[kHz] = 5[µs] between the channels, resulting in a total delay of 20[µs] delay between the first and the last sampled channel. In the sudden load change and piston impingement experiments presented in §2.6 the sampling frequency was 1[kHz]. It is assumed from now on that the delay of 20[µs] is negligible compared to the 2 Of
course, the weight of every cylinder was accurately measured but is not listed here.
8
2.2. Experimental setup
sampling interval of 1000[µs]. The resolution of the card was 12 bits, which means that 2 12 = 4096 values were ordered to the measuring interval. This results in an error of approximately 0.0244% due to the quantisation. LabView is a high-level programming language, which provides a quick and easy way of setting up the measurement parameters (sampling frequency, channels, calibration, etc.) and saving the data for later processing and visualisation. In the case of proportional valve control, the signal output was also produced with the appropriate built-in LabView function. It should be noted here that no software filtering was used and all the pressure and displacement measurements presented in this work are the ’original’ data. The vibration measurements presented in §4.5 were managed with the PULSE multichannel dynamic measurement system - which was kindly provided by the Department of Applied Mechanics of the Budapest University of Technology and Economics. The system is capable of recording acceleration and performing Fast Fourier Transform ’on the fly’, with a frequency range up to 10[kHz] and a resolution of 4[Hz]. The author acknowledges the help of Prof. G. Stépán and Mr. T. Gáspár for allowing to use the measuring system and for helping perform the measurements.
2.2.3 Material properties This section describes the identification process of the physical properties of the hydraulic oil used in the experiments. Some elements of the experimental rig are more than 20 years old and the hydraulic fluid inside was several times attenuated (when cleaning), mixed or replaced. Also, aging and temperature change may significantly influence the properties of the fluid. Thus, although fluid parameters are given in [20], it was necessary to determine the density, viscosity and bulk modulus anew.
Density Density is the mass of a unit volume of a substance. It is frequently used to identify an oil, or oil fractions, in the measurement of kinematic viscosity (absolute viscosity divided by density) and is also present in the equation for Reynolds number. During the density measurements, the amount of fluid in a measuring cylinder was systematically increased while its volume and mass were noted. The number of measurements was N = 13. The absolute error bound of mass measurement was ∆m = ±0.05[g] and assuming Student’s distribution with 95% confidence, we have σ m = 0.5[g]/2.228 = 0.2244[g], which gives a relative error between 1.095% and 0.117%. Similarly, we have ∆V = ±0.5[ml] and σV = 0.2244[ml], with relative error of 0.94% − 0.103%. Hence we conclude that both quantities are loaded with
an error of approximately the same order of magnitude and this implies that the least-square technique cannot be applied. Instead, we apply Wald’s technique [33] p.120. Our aim is to find the linear relationship between the mass and the volume of the fluid given by the relationship m(V) = ρV. Given the measured values V i and mi , i = 1 . . . N, we split them into two groups: 1 . . . r and r + 1 . . . N (we also assume that the data set in order following the magnitude of V i ). We compute the centre of mass for both sets, denoted by V 1 , m1 and V 2 , m2 . The slope of the fitted line is then given by ρ=
m2 − m 1 V2 − V1
=
1 PN i=r+1 mi N−r 1 PN i=r+1 Vi N−r
−
−
1 r 1 r
Pr
i=1
Pr
i=1
mi Vi
.
(2.1)
Given the standard deviations σm and σV , the standard deviation of the density σρ can be calculated with the square error spreading rule [33]. The results are "
kg ρ = 876.56 3 m
#
"
# kg and σρ = 0.597 3 . m
The dependence of the density on temperature and pressure is neglected in this work. 9
(2.2)
2. Experimental set-up and mathematical modelling of the system elements
Kinematic and dynamic viscosity Viscosity is the property of a fluid that causes it to resist flow, mechanically being the ratio of shear stress to strain rate. Viscosity may be visualised as a result of physical interaction of molecules when subjected to flow. It is absolutely necessary for hydrodynamic lubrication, and a suitable value is required for many other purposes. If the viscosity is too low, leakage flows increase; if the viscosity is too large, component efficiencies decrease because of additional power loss in fluid friction. The viscosity of liquids decreases significantly with temperature increase and increases, but to a much lesser degree, with increased pressure. The latter case (i.e. the dependence of the viscosity on the pressure) is not considered here. The viscosity of the oil was determined with the help of a Hoeppler falling sphere viscosimeter, type BH2/16478. The diameter of the inner tube was 15.945mm, the measured distance was 100[mm] and the diameter of the chosen sphere was 15.56[mm], resulting in a measuring time greater then 12 seconds for every measurement. The dynamic viscosity is calculated in accordance with the formula µ = t ρball − ρ f luid K,
(2.3)
where the constant K for the chosen ball was K = 0.12837. The density of the ball was ρ ball = 8140kg/m3, the density of the oil was ρoil = 876.56kg/m3, resulting in µ[mPas] = 0.932[mPa] × t[s] and ν[mm2 /s] =
µ = t[s] × 1.06[mm2/s2 ]. ρ
(2.4)
The actual measurement results are given in Fig.2.2, where every point represents the (arithmetic) mean of three measurements. The dependence of viscosity on temperature may be approximated by an equation of the form (see [60] p.11 eq.(2-9)) ν = ν0 e−λ(T −T 0 ) ,
(2.5)
where ν is the kinematic viscosity of the oil, λ is a constant depending on the fluid and T 0 is the reference temperature. By taking the natural logarithm of both sides, one gains a linear relationship between the logarithm of viscosity and the temperature and may use the technique of Wald (see [33]) to determine the unknown parameters. The result is ν0 = 92.79[mm2/s] and λ = 0.0401[C −1] and the curve is plotted in Fig.2.2 with red dashed line. As can be seen, linear approximation does not yield a satisfactory result in terms of accuracy. To improve the quality of approximation, an exponential function of second order was employed, given by 2
ν = ν0 e−λ1 (T −T 0 )+λ2 (T −T 0 ) .
(2.6)
The problem is then how to find ν0 , λ1 and λ2 if the reference temperature is chosen to be T 0 = 20oC. Although it is theoretically solicitous to employ the least square technique in this case 3 , it gives an excellent approximation. The actual values were found to be # mm2 , ν0 = 117.95 s "
λ1 = 0.0608[C −1] and λ2 = 3.45 × 10−4 [C −2 ],
(2.7)
and the resulting curve is plotted in Fig.2.2 with solid green line.
Bulk modulus and sonic velocity The concept of fluid bulk modulus is similar to elasticity in the case of the solid mechanics and it represents the compressibility of the fluid. The definition of the bulk modulus B is that the relative volume change of a fluid exposed 3 The
least square technique assumes that the the abscissa coordinates (in this case, the temperature) are not loaded by error. However, the standard
deviation of the temperature values was σT = 0.5o C.
10
2.2. Experimental setup
90
80
70
ν [mm2/s]
60
50
40
30
20
10 25
30
35
40
45
50
55
60
65
70
75
T [Celsius degree]
Figure 2.2: Kinematic viscosity of the oil as a function of the temperature. Blue points are measured values, red dashed line is (2.5) with ν0 = 92.79[mm2/s] and λ = 0.0401[C −1] and green solid line is equation (2.6) with ν0 = 117.95[mm2/s], λ1 = 0.0608[C −1] and λ2 = 3.45 × 10−4 [C −2 ]. In both cases, T 0 = 20[.oC].
to a pressure change is the same as the relative pressure change compared to the bulk modulus, i.e. ∆V/V = ∆p/B. If the pipe in which the fluid is located is elastic and also responds to the pressure change with some diameter change, this capacity-like effect should also be taken into account. In such cases, the effective bulk modulus of a fluid inside a thin-walled flexible pipe is given by 1 1 d = + , B B f luid δE pipe
(2.8)
where B f luid is the bulk modulus of the fluid, d is the diameter of the pipe, δ is the pipe wall thickness and E pipe is the elasticity modulus of the pipes material. Based on [20], the elasticity modulus of the B35 steel pipe is E B35 = 2.1×1010[Pa] and the elasticity modulus of the rubber hose is E r = 1.09 × 109 [Pa]. The bulk modulus of the hydraulic oil was not known but determined based on the measurements. In [20], the pressure wave in a long steel pipe was measured. The p sonic velocity was found to be a = 1137[m/s], which by taking into account that a = B/ρ and that the reduced bulk modulus was B = 1.135 × 109 [Pa] gives
B f luid =
BδE pipe = 1.615 × 109 [Pa]. δE pipe − Bd
11
(2.9)
2. Experimental set-up and mathematical modelling of the system elements
2.3 Hydraulic cylinder The mechanical model of the hydraulic cylinder is shown in Fig.2.3. The flow rates entering into chamber ’1’ and ’2’ are Q1 and Q2 , respectively. The half stroke length is denoted by L0 , the cylinder diameter is denoted by D and the diameter PSfrag replacements of the rod by d. The overall load force is denoted by F l(oad) . L0
L0
1
2
d
x Fload
D
Ql Q1
Q2
Figure 2.3: Mechanical model of the hydraulic cylinder.
Applying Newton’s second law on the piston we obtain m x¨ = A1 p1 − A2 p2 − F f + Fl ,
(2.10)
with chamber pressures p1,2 , cross-sectional areas A1,2 , F f internal forces (notably friction) and F l external (load) force. Continuity considerations for the fluid in chambers 1 and 2 lead to Q1 Q2
p1 − p2 V01 + A1 x dp1 + and Rh B dt p1 − p2 V02 − A2 x dp2 + . = −A2 x˙ − Rh B dt =
A1 x˙ +
(2.11) (2.12)
with B bulk modulus and V0,i = L0 Ai initial chamber volumes(i = 1, 2) and Rh is the hydraulic resistance, which is used for leakage flow Ql estimation (its actual value is taken from [20]). The above equations express the simple fact that the pressure change in the upper and lower chambers is due to (a) the entering (leaving) volume flow rates Q 1,2 and (b) the varying volume A x˙ and (c) the leakage flow from chamber ’1’ to ’2’. The behaviour of the friction force acting on the sealings on the cylinder and piston rod is highly complex, see e.g. [14]. A series of measurements was performed to determine the friction force. With the help of a proportional valve, different piston velocities were set and the pressures in the upper and lower chambers were recorded. The process was repeated with several loads m, resulting in several pressure differences. It was assumed that after a short initial transient during the rest of the motion (lifting and sinking), the velocity was constant. Thus, if the pressures in the upper and lower chamber are known and load is also given, the friction force can be simply calculated, by setting the acceleration x¨ to zero in (2.10). Fig.2.4 presents the dependence of the friction force on the piston velocity; the different measurements corresponding to different loads are depicted with different colours. The asymmetry of the curves for positive and negative velocities can be attributed to the asymmetry of the sealing element. Fig.2.5 depicts the dependence of the friction force on the pressure difference. As can be seen, the measurements are more stable for positive velocity range (positive friction force). Also, it is clearly seen that the friction force depends linearly on the pressure differences. This is in good agreement with [14], where the friction force was calculated as F f = g(z)pw bD π.
(2.13)
Here g(z) is the friction factor (originally, it was denoted by µ in [14], but we change this to avoid the confusion with viscosity). D denotes the mean diameter of the sealing, b is the width of the sealing element and p w is the sealing 12
1000
1000
800
800
600
600
400
400
200
200
F [N]
0
−200
−200 m = 165 [kg] m = 213 [kg] m = 243 [kg] m = 299 [kg] m = 334 [kg]
−400 −600 −800 −1000 −0.03
0
f
f
F [N]
2.3. Hydraulic cylinder
−0.02
−0.01
0
0.01
0.02
m = 165 [kg] m = 213 [kg] m = 243 [kg] m = 299 [kg] m = 334 [kg]
−400 −600 −800 −1000 10
0.03
15
20
v [m/s]
25
∆ p [bar]
30
35
40
Figure 2.4: Measured friction force as a function of the
Figure 2.5: Measured friction force as a function of the
piston velocity.
pressure difference.
m = 165 [kg] m = 213 [kg] m = 243 [kg] m = 299 [kg] m = 334 [kg] estimation
0.3
0.2
g(z) [−]
0.1
0
−0.1
−0.2
−0.3 −0.5
−0.4
−0.3
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
z × 107 [−]
Figure 2.6: Dimensionless friction factor as a function of the dimensionless velocity.
pressure, which is approximated by the chamber pressure difference. The friction force acting on the piston rod can be calculated by simply changing D to d in (2.13). The parameter z is a dimensionless velocity, defined as
z=
µv , pw b
(2.14)
where µ is the dynamic viscosity of the oil, v is the velocity of the piston, p w is the sealing pressure and b is the width of the sealing. The importance of parameter z is clearly seen in Fig.2.6; independently of the sealing pressure, all the measured curves reasonably coincide with a reasonable accuracy. The friction factor was approximated by a function of the form g(z) = α |z|β sign(z), where the parameters α and beta are determined with the least square technique and were
found to be α = 3.819 × 10−6 and β = −0.5373. However, for small z values, the friction factor is limited by g max = 0.25. Table 2.1 gives the physical properties of the cylinder. 13
2. Experimental set-up and mathematical modelling of the system elements
parameter
value
unit
explanation
D
50
mm
cylinder diameter
d
32
mm
piston rod diameter
2L
270
mm
total stroke of the piston
3 bar/ ms
leakage resistance
mm
width of the sealing element (estimated value)
Rh
1.04 × 10
b
10
7
3
dead volume at chamber 1 (estimated value) dead volume at chamber 2 (estimated value)
V01
19.6
cm
V02
11.6
cm3
Table 2.1: Physical parameters of the hydraulic cylinder.
14
2.4. Pilot-operated pressure relief valve
2.4
Pilot-operated pressure relief valve
The pressure relief valve is an essential part of every hydraulic system. Its function is to limit the maximum pressure in the system. According to the buildup of the valve, there are two main types: direct operated and pilot operated pressure relief valve. The direct operated valve operates with pre-loaded spool or poppet valve, but its use is limited by its poor pressure override characteristics and its unpredictable behaviour in the case of pressure levels close to the opening pressure. To improve this characteristics and to allow the exhaust of small amount of fluid in a stable way, a pilot stage is introduced, as shown in Fig.2.7. In such valves, the pilot poppet valve and the main valve are coupled through the fluid. This makes the dynamic analysis highly nontrivial (see e.g. [60], [91] or [25]). Not only it is difficult to mathematically describe the complex fluid-structure interactions (e.g. the flow forces on a poppet valve), but also the parameter identification (e.g. orifice and damping coefficients) is a challenging task.
Figure 2.7: Rexroth DBW pilot operated pressure relief
Figure 2.8: Simplified sketch of the valve for mathemati-
valve.
cal modelling.
The operation of the pilot operated pressure relief valve can be summarised as follows. Referring to Fig.2.8, the main stage of the valve consists of a spool valve and this main orifice is held closed by a soft seating spring. If the pilot valve is closed, the pressure inside the pressure relief valve is balanced. Should the force due to the inlet pressure p in exceed the force of the hard spring of the poppet valve, the pilot stage opens permitting flow to the return line and hence - due to the orifices - the pressure above the main spool drops and opens the main orifice. Flow is then bypassed to the return line in such a manner that the input pressure is kept at a value set by the pilot spring. To be able to simulate the dynamic behaviour of the relief valve, the motion of the main spool and the pilot stage poppet valve is to be computed. This suggests the need of studying the relationship between the flow rates and pressures in different segments of the valve and the forces on the valve bodies simultaneously. In the case of the main stage spool valve, due to the simple geometry and the chamfered-type valve seat, the flow rate and the forces are easily calculated with the help of standard relationships given in the literature. Contrary to this, there are several issues which make the study of the pilot poppet valve challenging; • the valve seat and the valve body forms an orifice with non-standard geometry and the corresponding flow rate pressure drop relationship is to be determined • there is a fixed orifice (Orifice 2 in Fig.2.8) before the poppet valve whose influence is to be studied and • contrary to the main stage spool valve, where the fluid force acts mainly on the front of the valve body where 15
2. Experimental set-up and mathematical modelling of the system elements
the upstream pressure distribution can be considered as uniform, in the case of the poppet valve, the pressure distribution along the cone changes significantly. For a realistic force estimation, an estimation of this pressure distribution along the cone should be known. The above issues suggest that the full 3D velocity and pressure field is needed for determining the fluid forces on the poppet valve, for example with the help of a CFD code. Keeping in mind that the final aim is to establish a relief valve model for system simulation, this approach cannot be fully employed due to its large computational demand. Instead, analytic calculations will be performed giving a realistic structure of relationships and the parameters are then identified based on the results of CFD computations. During the last decade, several models have been proposed for pressure relief valves with experimental verification, e.g. [72], [16], [76]. The main modelling headlines (flow rate equations, pressure drop calculation, etc.) are given and validated in these articles. However, the author is not aware of any proper description of damping on a poppet valve, which has a fundamental effect on the dynamic performance and stability of the valve. Mostly this damping force is a model parameter, although it is not directly controllable during the manufacturing. In [91], Zung et al. propose a damping model based solely on blueprint dimensions and achieves an excellent agreement with the measurement results. However, due to the geometrical differences, his approach cannot be directly employed in our case. This section is organised as follows. The mathematical model of the relief valve consists of the Newtonian equations of motion for the main stage spool body and pilot stage poppet valve body. These equations are given in §2.4.3. One of the force components in these equations is the fluid force acting on the valve bodies, which is determined in §2.4.2. However, to calculate e.g. the pressure force on a valve body, the pressure distribution is needed and also the flow rate through the valve should be known. These relationships are determined in §2.4.1 with the help of analytical and numerical techniques. Finally, the static performance curve of the relief valve is computed and the result is compared to that one in the catalogue in §2.4.4.
2.4.1 Pressures and flow rates Referring again to Fig.2.8, we assume the fluid to be incompressible everywhere inside the valve. The pressure losses are concentrated into the orifices as the cross-sectional areas of Channel 1 and 2 are large enough to neglect the pressure loss inside4 . Because of the large diameters, we also assume uniform pressure distribution in Channel 1 and 2. The laminar or turbulent nature of the flow in an orifice influences the flow rate - pressure drop relationship fundamentally. If the flow is laminar, the dependence of the flow rate on the pressure difference is linear while if the flow p is turbulent, we have Q ∝ ∆p. Orifices 1, 2 and 3 are with short, sharp-edged seats, suggesting turbulent flow. Similarly,
due to the sharp edges and the complex geometry, the flow through the main stage (spool valve) is also turbulent, (see e.g. [25]). These volume flow rates are calculated by using the orifice equation: Q1 = Q2 = Q3 = Q pl =
Cd,1 A1 Cd,2 A2 Cd,3 A3
q √ 2 ∧ ρ
q √ 2 ∧ ρ
q √ 2 ∧
Cd,pl A pl (x pl )
ρ
p1 − p 2 p3 − p 2
q √ 2 ∧
Qout = Cd,out Aout (x sp ) 4 e.g.
pin − p1
ρ
p2 − pout
q √ 2 ∧ ρ
pin − pout
de f.
√ ∧
de f.
√ ∧
de f.
√ ∧
= R1
= R2 = R3
de f.
= R pl
de f.
= Rout
pin − p1 ,
(2.15)
p1 − p 2 ,
(2.16)
p3 − p 2 ,
(2.17)
√ ∧
p2 − pout
√ ∧
and
pin − pout .
for Channel 1, with L = 60[mm], D = 4[mm], λ = 0.02[−], ρ = 877[kg/m 3 ] and Q = 10[dm3 /min], the pressure loss is ∆p = λ DL
0.0145[bar]
16
(2.18) (2.19) ρ Q2 2 A2
=
2.4. Pilot-operated pressure relief valve
Cd,i represents the discharge coefficient, Ai is the flow-through area and
√ ∧
x represents the signed square root function5.
Orifice 1 and 3 are identical restriction elements, thus their parameters are the same. The flow-through areas are fixed for Orifice 1-3, while for the pilot valve and the main spool, we have (see Fig.2.16 and Fig.2.17) x pl and A pl = πx pl sin α pl d2,in − sin 2α pl 2 x sp Aout = πx sp sin α sp d sp − sin 2α sp . 2
(2.20) (2.21)
The discharge coefficient C d includes the contraction coefficient (the ratio of the upstream area and the area of the vena contracta) and also the velocity coefficient (modelling the viscous losses in the orifice jet). In general, the discharge coefficient is difficult to compute but solutions have been made for orifices with regular geometry and it was shown that for sharp-edged orifices with turbulent flow the theoretical value of C d = π/(π + 2) ≈ 0.611 can be used, see [60] p.42. For other cases, C d can be obtained by steady-state measurements as described in [75]. Orifice 1 and 3 are with regular sharpedged geometry, thus the theoretical value will be used; C d,1 = 0.61 for fully developed turbulent flow. The geometry of the main stage is also regular and the same value will be used for turbulent flow. It is not trivial how to estimate the discharge coefficient for Orifice 2 as the inflow side is a regular converging channel and the outflow is a sudden diameter step. Also, the calculation of the flow through the pilot valve body and the seat is not straightforward due to the varying gap width. To identify the discharge coefficients and to compute the flow forces on the pilot valve body a series of steady-state CFD computations of the whole pilot stage was performed with the commercial CFD code CFX 5.7. Valve displacements between 0.1[mm] and 0.5[mm] were set with pressure differences between 1[bar] and 100[bar]. As shown in Fig.2.9 and Fig.2.10, the computational domain started at Channel 1 and ended at Channel 2, in both cases 10[mm] after entrance in order to be able to prescribe stable uniform pressure field. The pressure at the inlet of Channel 1 was then systematically increased to 1[bar], 2[bar], 5[bar], 10[bar] and so on up to 100[bar] for every valve displacement. The symmetry of the geometry was exploited; only one half of the valve was built up and also the geometry was simplified: as under steadystate conditions the main spool is at rest, Orifice 3 was omitted in these computations. The SST model was used for turbulence modelling [8]. The mesh size varied from 0.5 million elements to 1.2 million elements (smaller valve openings require finer resolution around the valve seat). For a given pair of input pressure p 1 and output pressure pout , the pressure before the pilot valve p2 , the mass flow rate m ˙ and the force components on the pilot poppet valve were extracted from the results. Plotting the flow rate versus p1 − p2 gives the Q(∆p) relationship of fixed Orifice 2, while plotting p2 − pout vs. the flow rate gives flow-through relationship of the pilot valve. A typical pressure distribution and velocity vector plot is shown in Fig.2.11. We observe a low-pressure region underneath the seat as described by earlier experimental studies, e.g. [75], [82] or [12]. This phenomenon is also confirmed by experiments, e.g. dissolved air in the oil was found in the measurements of Stone [75]. The actual configuration differs significantly from the usual geometry in the corresponding literature as there is a nozzle close to the poppet valve whose jet has a considerable impact on the flow field around the valve body. Also, a large vortex is present around the cone (see Fig.2.11). These two effects result in a strongly varying pressure distribution on the surface of the cone, contrary to the usual configuration in the literature, where uniform upstream pressure distribution is assumed. Figures 2.12 and 2.13 depict the pressure and velocity distribution for 40[bar] set pressure along the polyline shown in Fig.2.11, respectively. In both figures, dashed vertical red lines represent the entrance and exit of the nozzle. The black dashed vertical line depicts the position of the valve seat. The pressure and velocity distribution do not differ qualitatively for different valve displacements. The pressure drops quickly in the nozzle and between the nozzle exit and the nose of the cone it is approximately constant (as it is well-known from the theory of jets, see e.g. [70] p.171). As the ∧ 5 √ x
= sign(x)
√ |x|
17
2. Experimental set-up and mathematical modelling of the system elements
jet hits the nose of the cone, there is an increase in the pressure due to the deceleration of the fluid.
Valve body
Orifice 2
Channel 2 Outflow Channel 1
Valve seat
Inflow
Figure 2.9: The full CFD geometry.
Figure 2.10: CFD geometry, a zoom close to the nozzle and nose of the pilot valve body.
Negative pressure region Polyline
Figure 2.11: Pressure distribution and velocity vector plot for x pl = 0.3mm and pin = 40bar. Note that the pressure interval between −2 . . . 3[bar] is exaggerated and red regions correspond to pressures between 3 . . . 40[bar].
18
2.4. Pilot-operated pressure relief valve
45
100
0.2mm 0.3mm 0.4mm 0.5mm
40 35
90
70
30
60
v [m/s]
25
p [bar]
0.2mm 0.3mm 0.4mm 0.5mm
80
20 15
50 40 30
10
20
5
10
0
0
−5 50
−10 50
52
54
56
58
60
62
64
66
52
54
x [mm]
56
58
60
62
64
66
x [mm]
Figure 2.12: Pressure distribution along the cone body
Figure 2.13: Velocity distribution close to the cone
with sharp-edges seat for several valve body displace-
body for several valve body displacements.
ments. Discharge coefficients During an unsteady process the flow presumably passes in and out of laminar and turbulent regions. Thus, we need an orifice model which describes both regimes simultaneously. Such a model is described in [60] for short orifices, where the discharge coefficient can be identified as (p.42) Cd = Cd =
1/2 −1/2 1.5 + 13.74 13.74 Re 64 −1/2 2.28 + Re
for
Re > 50
and
(2.22)
for Re < 50,
where it was assumed that the orifice is a short tube, i.e. L/D ≈ 1. The Reynolds number is defined as Re =
vDh (Q/A0 ) Dh = , ν ν
(2.23)
where A0 is the flow-through area and Dh is the hydraulic diameter. For a circular orifice with diameter d, the hydraulic diameter is Dh = d. Another model for the discharge coefficient is described in [15] by Borutzky, where the equation Cd = cturb √
√
Re √ Re + Ret
(2.24)
is proposed, with parameter values cturb = 0.61 and Ret = 9.33 for sharp-edged orifices. The two models (2.22) and (2.24) give nearby values for the practical range of the Reynolds number. Both models provide a linear relation through the orifice for small velocities while for turbulent flows, they match the conventional square root characteristics. Also, the transition from the laminar to the turbulent regime is smooth in both cases. The price for this is that the flow rate cannot be computed explicitely in (2.15) - (2.19) and iteration is needed at every time step. In what follows, the latter relation (2.24) is employed as (a) it is simple, (b) there is no need to switch back and forth between two functions while iterating and (c) the parameters employed (cturb and Ret ) have physical meaning. The dependence of the discharge coefficient on the Reynolds number is also important from the viewpoint of dynamic numerical simulation. Assume that during an unsteady process, the pressure difference changes its sign in one of the orifice equations, say in (2.15). Existence and uniqueness problems arise around p in − p1 = 0 due to the infinite slope of the square root function around zero. However, if the discharge coefficient is assumed to depend on the Reynolds number rather than being constant, the slope around zero flow rate is well-defined and no problems arise when integrating the dynamical equations. 19
2. Experimental set-up and mathematical modelling of the system elements
Although the primary aim of the CFD computations was to study the flow forces on the poppet valve, additional valuable information could be gained concerning to the behaviour of fixed Orifice 2. Fig.2.14 depicts the discharge coefficient for Orifice 2 as a function of the Reynolds number. It is clearly seen that the discharge coefficient varies significantly and that the Borutzky formula (2.24) gives a satisfactory approximation. The parameters c turb and Ret were calculated by nonlinear least-square technique, i.e. those parameter values were calculated, at which the sum of the square of the errors between the CFD data points and (2.24) is minimal. The effect of constant and Re-dependent discharge coefficient is demonstrated in Fig.2.15, where the red curve indicates (2.16) with C d = 0.78 (mean value of the CFD runs) and Borutzky’s formula. It is evident that the dependence of the discharge coefficient on the Reynolds number should be taken into account in order to be able to compute the flow rate through the orifice. 2.5
0.8 0.7
2
0.6
Q [l/min]
Cd[−]
0.5 0.4 0.3
x=0.1mm x=0.2mm x=0.3mm x=0.4mm x=0.5mm Borutzky
0.2 0.1 0 0
500
1.5
1
x=0.1mm x=0.2mm x=0.3mm x=0.4mm x=0.5mm Cd mean Cd Borutzky
0.5
1000
0 0
1500
5
10
15
20
25
30
35
40
p1−p2 [bar]
Re[−]
Figure 2.14: Discharge coefficient of Orifice 2 as
Figure 2.15: Volume flow rate vs. pressure difference
a function of the Reynolds number and the fitted
for Orifice 2. Crosses denote CFD data points, solid
Borutzky approximation (2.24), with cturb = 0.86 and
red line represents (2.16) with an average C d value
Ret = 30.
while solid blue line denotes (2.16) with (2.24). The discharge coefficient for the cone-seated orifice is harder
pout
to identify as it depends on both the pressure difference and the gap width. Defining a Reynolds number with the help of the gap width
Fpl
H
H = x pl sin α (see Fig.2.16) and plotting the discharge coefficient as
x pl
a function of it results in different curves for every valve displace-
β
αpl
vout
ment rather than one unique curve. Thus, we adopt the approach of [39]. Define the loss coefficient ζ as ζ=
vin p2
p2 − pout , ρ 2 2v
where v =
4Q . 2 π d2,in
(2.25)
By plotting log ζ vs. the dimensionless valve displacement log(x pl /d2,in )
d2,in
and using the least square technique, we obtain
Figure 2.16: Geometry for calculating loss coeffi-
ζ=
cient and the control volume for fluid forces cal-
0.369 2.3 x pl /d2,in
for 0.01 ≤
x pl ≤ 0.1. d2,in
(2.26)
Note that based on experiments, [39] gives
culation. ζ = 0.6 +
0.15 (x pl /d2,in )2
for 0.125 ≤ x pl /d2,in ≤ 0.4.
(2.27)
Unfortunately, the cone angle is not given in [39], thus no direct comparison between the measured and calculated values is possible. Table 2.2 compares CFD results to (2.26). 20
2.4. Pilot-operated pressure relief valve
x pl [mm]
0.2
0.3
0.4
0.5
CFX
841.7
316.1
137.8
111.4
(2.26)
786.1
309.4
159.6
95.5
Table 2.2: Comparison of loss coefficient computed with CFX and fitted formula (2.26)
Having defined the relationships between volume flow rates and pressures, we turn now to the actual computation issues. We have four unknown flow rates (Q1..3 and Q pl ) and three unknown pressures p1..3 . Beside the three orifice equations (2.15)-(2.17) and the pilot valve restriction equation (3.11), we have the continuity equations due to the incompressibility Q1 = Q2 and Q2 = Q pl + Q3 , which means that one more equation is needed. This is simply Q 3 = A sp v sp , (where v sp is the velocity of the spool body, obtained from the equation of motion) which is again a continuity equation for the chamber above the main spool body. The equation to be solved is Q1 = R pl
q
pin − pout − σQ21 + A sp v sp ,
where σ =
R1 2 + R 2 2 , R1 2 R2 2
(2.28)
which can be solved e.g. by function iteration. Once Q1 is known, the rest of the unknowns can be easily computed. Numerical experiments show that for reasonably good initial conditions, 3 or 4 iterations steps are sufficient for function iteration to achieve 0.1% accuracy in Q1 even when considering the discharge coefficients to be Re-dependent.
2.4.2 Fluid forces on the valve bodies The fluid forces are sometimes called Bernoulli forces [7] as they arise from the dynamics of the fluid flow. Steady state flow forces are due to (a) the unbalanced pressure distribution around the body, (b) the wall shear distribution on the surface of the body and (c) the momentum forces, i.e. the angle of the Vena Contracta as the fluid enters and leaves the valve. Unsteady flow forces are due to the acceleration of the fluid in the chamber. In what follows, these forces are calculated in the case of the main spool valve and the pilot poppet valve. In the case of the main spool valve, due to the simple geometry, the formulae available in the literature can be directly applied. In the case of the pilot stage, detailed computations are carried out and the results of the obtained formulae are compared to the CFD results. The steady flow force acting on the main spool can be obtained with the help of the momentum equation [51] as shown in Fig. 2.17. One obtains F f,sp =
ρv2in Ain
−
ρv2out Aout
cos θ + pin Ain − pout Aout cos θ,
(2.29)
where θ is the jet angle, and Aout is given by (2.21). For sharp-edged valve and housing pairs (typically in the case of servovalves), the jet angle
pout
Ff,sp
vout θ
xsp seat L 0 inlet
spool body α sp
pin
L vin
strongly depends on the spool displacement. Due to the chamfered valve seat, we assume θ ≈ α sp , which is the usual assumption in the correspond-
Figure 2.17: Geometry and control volume
ing literature (see [60]). The pressure integral and the viscous forces along
for the calculation of the fluid forces on the
the chamfered part of the valve body are neglected due to the shortness of
main spool.
the seat. Neglecting the flow rate through the pilot stage: Q pl Qout →
Qin ≈ Qout , taking into account that Qin = vin Ain , Qout = vout Aout and employing (2.19), (2.29) turns into ! cos θ 1 − F f,sp = pin Ain − pout Aout + ρQ2out Ain Aout " !# !# " Aout Aout Aout 2 2 = pin Ain 1 + 2Cd,out − cos α sp − pout Aout 1 + 2Cd,out − cos α sp . Ain Ain Ain 21
(2.30)
2. Experimental set-up and mathematical modelling of the system elements
In the corresponding literature it is often assumed that the momentum force can be neglected. Based on the terms in square brackets in the last row of (2.30) it is easy to decide whether this holds in our case or not. By exploiting (2.21), we have Aout /Ain ≈ 4 x sp /d sp sin α sp , which varies6 in the range of 0 . . . 0.435. It is easy to calculate that as x sp varies from 0mm to 4mm, the factor of pin Ain varies from 1 to 0.8485 and that of pout Aout from 0.095 to 0.652. Hence we conclude that the momentum forces in the case of the main spool have to be taken into account. The transient flow force can be obtained by evaluating the unsteady integral part of momentum principle (equation (7.9) in [51]). Uniform mean axial pressure distribution is assumed inside the control volume and only the inlet part of the control volume is considered with cross-section Ain and length L, see Fig. 2.17. This gives s ! Z 2 p dQ ∂ d∆p = Cd,out ρ L π d sp sin α sp (ρv)dV ≈ ρL , Ft,sp = ∆p x˙ sp + x sp ∂t V dt ρ dt
(2.31)
therefore we conclude that the transient force is directly proportional to the valve velocity x˙ sp and the pressure changes d(∆p)/dt. Note that mechanically thinking, the velocity term represents damping force. Following Merritt, we neglect the pressure term in what follows as there is little evidence to indicate that it contributes substantially to valve dynamics, see [60] p.104. The quantity L is the axial length between the mounting manifold and the valve front, i.e. L = L 0 + x sp and is called the damping length and it is a strictly positive quantity. Let us now calculate the fluid force acting on the pilot poppet valve. The main difficulty is that the flow field around the conical body is too complicated to assume some simple pressure distribution and then integrate it, recall Figures 2.12 and 2.13. Thus, we give only a rough estimation here and emphasise that this issue needs further study but it does not belong to the scope of the thesis. Applying the steady momentum equation, we have Z Z F f,pl = (vρ) (vdA) + p dA − S = F f,pl,mom. + F f,pl,pres. − F f,pl,visc. , A
(2.32)
A
where the bold letters denote vector quantities and the forces denote only the axial component of the total (vectorial) force. The control volume is depicted in Fig. 2.16. The easiest bit is to calculate the momentum forces: ! Z cos α pl 1 , − F f,pl,mom = (vρ) (vdA) = ρQ2 Ain Aout A
(2.33)
where Ain and Aout are the cross-sectional areas of inlet and outlet, the latter one is given by (2.20). The pressure force can be obtained by integrating the pressure distribution over the control volume’s surface. As only the axis-wise component is needed, we have F f,pl,pres. =
Z
Ain
p dA − cosα
Z
Aout
p dA = Ain p2 − Aout pout cos α pl ,
(2.34)
where α pl is the half-cone angle of the pilot valve body. The above relationship assumes constant pressure distribution along Ain , which is far from reality. However, let us not forget that first the volume flow rate through the whole pilot stage is computed by virtue of (2.28), where only well-defined input and output pressures and discharge coefficients are needed. Pressure p2 is calculated from the volume flow rate and thus gives a good approximation of the average pressure around the cone body. The pressure and the viscous forces were also extracted from the CFD results and are depicted on Fig 2.18 and Fig 2.19. For small pressure differences, their values are in the same order of magnitude but for larger pressure differences pressure forces become dominant. In Fig 2.18 also (2.34) is given and we conclude that it gives only limited results and the error may be up to 400% for larger valve openings. A possible explanation of this phenomena is that (2.34) 6x
sp,max
= 4mm, d sp = 26mm, α = 45o
22
2.4. Pilot-operated pressure relief valve
25
20
1.8
0.5mm 0.4mm 0.3mm 0.2mm
1.6 1.4
0.5mm 0.4mm 0.3mm 0.2mm
1.2
F [N]
F [N]
15 1 0.8
10 0.6 0.4
5
0.2
0 0
20
40
pin [bar]
60
80
0 0
100
20
40
p [bar]
60
80
100
in
Figure 2.18: Pressure force on the pilot valve. Points
Figure 2.19: Viscous force on the poppet valve, results
denote CFD computations, solid lines represent (2.34).
of CFD computations.
assumes constant pressure distribution around the cone while in the reality, we have a pressure peak close to the nose of the valve body (due to the stagnation point), followed by a slight depression (the fluid accelerates again), see Fig 2.12. As the valve opens and the mean pressure drops, these deviations become larger in relative terms and this is not taken into account in (2.34). An analytical description of the viscous force requires that the velocity field around the cone body is known. However, such a formula seems to be cumbersome to obtain, thus we neglect the viscous force and accept that this introduces an error of approx. 10% compared to the pressure force.
2.4.3 Equations of motion Main valve Applying Newton’s 2nd law to the main spool valve, we have m sp x¨ sp = F f,sp + Ft,sp − p3 A sp − s sp x0,sp + x sp − K sp x˙ sp ,
(2.35)
where x sp is the displacement of the spool, m sp is the mass of the spool, F f,sp represents the fluid forces defined by (2.30), Ft,sp is the transient flow force given by (2.31). The pressure force on the upper part of the spool is A sp p3 , K sp is the viscous damping coefficient (yet to be defined), s sp is the spool spring constant and x0,sp is the precompression of the spring. Table 2.3 gives the parameters of the main spool. The spool mass, diameter and length were directly measured. Spring stiffness was calculated by 20 measured displacement and force data pairs, with the help of the least square technique. The radial clearance between the spool and the housing was measured with the help of a micrometer 10 times. The damping on the body is important for quick and stable operation. If the damping is too low, oscillation die slowly; insufficient damping in nonlinear cases may lead to self-excited oscillations. However, if the damping is too large, the valve opens and closes slowly and may not be suitable for overpressure prevention. It is useful to calculate the eigenfrequency of the undamped system f and the critical damping K crit , which separates periodic and aperiodic motions (see [58], p.61) ω=
r
# " s sp rad , = 646.6 m sp s
f =
ω = 102.9[Hz] and 2π
Kcrit. = 2m sp ω = 84.4
Ns m
.
(2.36)
In our case, damping forces arise in three different ways: (a) due to the oil film between the spool surface and the housing, (b) due to the oil volume between the mounting manifold and the spool, because as the spool moves, this amount of oil 23
2. Experimental set-up and mathematical modelling of the system elements
parameter
value
unit
explanation 3
ρ
878
kg/m
ν
40
mm2 /s
kinematic viscosity of the oil
µ
0.03512
kg/ms
dynamic viscosity of the oil
m sp
0.0653
kg
main spool mass
d sp
22
mm
main spool diameter
−4
2
density of the oil
A sp
5.31 × 10
x sp,max
4
mm
maximal spool displacement
s sp
27.3
kN/m
spool spring stiffness
L sp
40
mm
total length of the spool in contact with the bore
α sp
45
degrees
seat and spool chamfer angle
Cr
0.1
mm
clearance between the spool and the bore
m
main spool cross-sectional area
Table 2.3: Physical parameters of the main spool.
has to be accelerated and decelerated and finally (c) due to Orifice 3. The viscous damping coefficient due to the clearance between the valve body and the housing can be estimated in the following way. Assuming a perfectly centred spool in a bore (see [60] or (3.24)) the damping force on the main spool reads Fd,sp =
πd sp L sp µ de f. x˙ sp = K sp x˙ sp = 0.971[N s/m] x˙ sp, Cr
(2.37)
where d sp is the diameter of the spool, L sp is the total spool length in contact with the bore and C r is the clearance between the spool and the bore. Recall that the transient flow force was also proportional to the velocity (see (2.31) and the corresponding text) Ft,sp = Cd,out ρLπd sp sin α
s
2 ∆p x˙ sp = 32.74[N s/m] x˙ sp, ρ
(2.38)
where L = 0.02[m], ∆p = 40[bar] and an average value of C d = 0.8[−] was assumed. Additional damping force arises on the main spool due to Orifice 3. To get a feeling of the order of magnitude of this force, let us calculate the flow rate through Orifice 3: QOri f ice 3 = A sp x˙ sp = Cd AOr.3
s
2 (p3 − p2 ), ρ
(2.39)
from which we conclude that the reaction force on the spool is F p3 = A sp p3 = A sp p2 +
A3sp ρ 2 Cd2 A2Or.3
x˙2sp = A sp p2 + 471 × 103 [N s2 /m2 ] x˙2sp ,
(2.40)
where the only new parameter is the flow-through area of Orifice 3. As its diameter is 0.6[mm], we have A Or.3 = 0.282[mm2]. The resulting damping term is extremely large and it is clear that even if the above estimation is inaccurate, Orifice 3 gives the dominant damping force. However, another interesting feature is that this damping is proportional to the square of the spool velocity, which means that it does not influence the stability of the equilibrium of the spool. Thus, although the damping through Orifice 3 gives the dominant part of the damping force, in theoretical (stability) calculations, the smaller linear damping terms given be (2.37) and (2.38) cannot be omitted.
24
2.4. Pilot-operated pressure relief valve
Poppet valve The equation of motion for the pilot cone is m pl x¨ pl = F f,pl − Fd,pl − s pl x0,pl + x pl ,
(2.41)
where F f,pl stands for the fluid forces given by (2.32), F d,pl denotes the damping force (yet to be defined), s pl is the pilot spring constant and x0,pl is its precompression. Table 2.4 gives the parameters of the pilot valve. The pilot mass and its geometry was directly measured. Spring stiffness was calculated by 20 measured displacement and force data points, with the help of the least square technique. parameter
value
unit
explanation
m pl
10.6
g
pilot equivalent mass, see (2.42)
r2,in
2.8
mm
inlet housing diameter
x pl,max
6
mm
maximal pilot displacement
s pl
811
kN/m
pilot spring stiffness
α pl
21
degree
half cone angle
β
45
degree
seat chamfer angle
D spring
13
mm
spring outer diameter
d spring
4
mm
spring steel wire diameter
n
10
−
number of spring turns
Table 2.4: Physical parameters of the pilot valve, see Fig. 2.16.
Note that the mass of the cone body mcone = 4.1[g] is comparable to the mass of the spring, which is m spring = 19.6[g]. In such cases, the mass of the spring should also be taken into account when calculating the eigenfrequency as follows. Let us calculate the overall kinetic energy of a system consisting of a mass m m moving with velocity vm (t) attached to a spring with specific mass m s [kg/m] and total length of L0 , whose other end is fixed. We assume a linear velocity distribution along the spring, i.e. v s (x, t) = vm (t) Lx0 ; v s (0, t) = 0 and v s (L0 , t) = vm (t). The total kinetic energy is E=
1 mm v2m + 2
Z
x=L0 x=0
1 2 1 1 v2 v s (x, t)m s dx = mm v2m + m s m2 2 2 2 L0
Z
x=L0
x2 dx =
x=0
m spring 1 2 , vm mm + 2 3
(2.42)
hence we conclude that the energetically equivalent mass can be calculated as m pl = mcone + m spring /3 = 10.6[g]. Contrary to the main spool, the pilot cone body is poorly damped. One way damping occurs is that the spring has relatively large dimensions and as the valve body moves, so does the spring and this gives rise to a drag force. Let D spring denote the main diameter of the spring, d spring denote the diameter of the steel wire of the spring and n stand for the number of turns. Let us assume that the velocity distribution is linear along the axis of the spring and let us also assume that the total surface of the spring perpendicular to the flow (n D spring πd spring ) is distributed continuously along the length. The drag force on the spring reads F=
Z
L 0
D spring πd spring n ρ dx, cD (Re) v2 (x) 2 L spring
(2.43)
where cD = 6/Re = 6ν/d v(x) is the drag coefficient on a cylinder [88], p.183. Assuming the again linear velocity distribution along the spring, i.e. v(x) = v pl x/L spring , we have F=
Z
L spring 0
D spring πd spring n 6ν ρ 3 v(x)2 dx = µD spring πn v pl = 0.0215[N s/m] v pl. d spring v(x) 2 L spring 2 25
(2.44)
2. Experimental set-up and mathematical modelling of the system elements
Again, it is useful to calculate the eigenfrequency of the undamped system f and the critical damping K crit , which gives ω=
r
# " s pl rad , = 8746 m pl s
f =
ω = 1392[Hz] and 2π
Kcrit. = 2mω = 185.4
Ns m
.
(2.45)
We conclude that the pilot valve is poorly damped and oscillations after e.g. an opening transient will die slowly. Although the flow rate through the pilot stage is small, due to the nonlinearities in the systems, the pilot stage may excite the main stage resulting in large volume flow rate oscillations. The stability of the pilot stage is studied later on in this work.
2.4.4 Static performance curve The evaluate the steady-state performance of the model, (2.35) and (2.41) were solved simultaneously under steady-state conditions, i.e. x˙ = 0 and x¨ = 0. For a given overall pressure drop p in − pout , the computation was performed as follows. First, by assuming constant discharge coefficients and 0.1[mm] valve displacements, with the help of (2.28), the pilot stage flow rate Q1 was calculated. As under steady-state conditions Q3 = A sp v sp = 0, we have Q1 = Q2 = Q pl . Then, the pressures p1..3 were determined by virtue of the orifice equations (2.15)-(2.18) and then the steady-state fluid forces were calculated as described in §2.4.2. Once the total forces are known, the valve displacements x pl,sp can be re-calculated by means of (2.35) and (2.41) (under steady-state assumptions: x pl,sp = 0 and v pl,sp = x˙ pl,sp = 0), which results in a balanced spring force. With the new valve displacements, new discharge coefficients and flow rates (Q 1 , Q2 and Q pl ) are determined and the process is restarted. Numerical experiments show that after a maximum number of five iterations 7 , both flow rates and valve displacement changes were below 0.1%. Fig.2.20 depicts the numerically calculated performance curves, which can be compared to those ones given in the catalogue, see Fig.2.21. The parameter of the curves is the pilot spring precompression, i.e. the preset pressure level. As seen, the main features of the performance are described with a reasonable accuracy; the sudden opening characteristics (small Q values) and the slight increase in the losses as the preset pressure is increased. However, due to the poor inner loss estimation, the calculated curves are too sharp in the small flow rate region and also the slight backward bending in the case of large flow rates is not captured (bottommost curve on the right-hand side figure, large flow rate). The dynamic performance of the model is presented in §2.6. 350
300
[bar]
out
200
p −p
150
in
in
p −p
out
[bar]
250
100
50
0 0
50
100
150
200
Q [l/min]
Q [l/min]
Figure 2.20: Computed static performance curve of the
Figure 2.21: Static performance curve of the relief valve,
relief valve.
taken from the catalogue.
7 However,
if the valve displacements and flow rates corresponding to one set of parameters (say, pressure difference) are already known, it takes
only one ore two more iterations to calculate the steady-state parameters corresponding to a slightly varied set of parameters.
26
2.5. Transmission line modelling
2.5 Transmission line modelling In the analysis of fluid control systems, the behaviour of the transmission lines connecting the components is essential to the understanding and prediction of the system’s dynamic characteristics. Even in steady-state analysis, it is rarely sufficient to consider the pipes (notably the fluid inside) as perfectly rigid, massless and lossless elements. The importance of dynamic pipe flow is confirmed by the fact that the first symposium on water hammer was held in 1933 by the American Society of Mechanical Engineers and that the earliest attempts to describe the observed water hammer effects are due to Joukowsky (1889) and Allievi (1903) by means of the one-dimensional wave equation. Nowadays, a large spectrum of models and tools are available, starting from the analytical solutions of the 1D wave equations to the most up-to-date CFD codes. It is clear that a ’too simple’ model is not sufficiently accurate for engineering applications and that it is not (yet) possible to incorporate a full 3D CFD code when analysing a whole system with several elements. Thus, the first step is to summarise the demands we make for modelling the unsteady pipe flow in a hydraulic system. • 1D model. The supreme aim of the model is to predict the longitudinal velocity and pressure distribution along the pipe, thus, the only important velocity component is the axial one. • Laminar flow. In hydraulics, the flow velocity is typically around a few [m/s]-s, which implies laminar flow. Whereas the criterion for the occurrence of turbulence in steady flow is simply the magnitude of the Reynolds number, in unsteady flows neither the criterion for, nor the manner how turbulence occurs is well understood. A note on the paper from Scarton and Rouleau [67] suggests an empirical relationship for the onset of turbulence in a compressible fluid undergoing steady-periodic laminar oscillation, experimental validated by Karam and Franke [42]. • Nonlinear convective terms are negligible. Again, this is a consequence of the low velocity assumed a priori and the relatively high wave celerity. • Rigid pipe. If needed, effects due to the elasticity of the conduit can be taken into account by using the effective bulk modulus of the fluid and the pipe [34] instead of using the isothermal sonic velocity in the fluid.
2.5.1 Governing equations Consider the flow of a compressible fluid in a circular duct with constant cross-sectional area A and inner radius R. Thus, axis-symmetry and the lack of radial velocity components can be assumed. On the other hand, only such small perturbations are to be treated, which allow us to apply an acoustic (first-order) approximation. This way, the continuity equation becomes ∂ρ ∂v + ρ0 = 0, ∂t ∂x
(2.46)
where x is the axial coordinate, t stands for time, ρ0 is the reference density and v(x, r, t) is the only, axial velocity component along the x-axis, while r denotes the radial coordinate. Pressure p(x, t) and density ρ(x, t) are assumed to be uniform across the pipe, the viscosity ν is assumed to be constant (independent from the other state variables). For a detailed explanation, see [34] or [51]. The equation of motion becomes ! 1 ∂p 1 ∂ ∂v ∂v r . =− +ν ∂t ρ0 ∂x r ∂r ∂r
(2.47)
The terms in equation (2.47) can either be absolute velocities or perturbations of a mean value. No-slip boundary condition is applied, i.e. v(x, R, t) = 0 and symmetry:
∂v ∂r (x, 0, t)
= 0. For a slightly compressible fluid in an elastic pipe, the equation 27
2. Experimental set-up and mathematical modelling of the system elements
of state can be written as B dp = a2s = , dρ ρ0
where
1 1 D = + . B f luid δE pipe B˜
(2.48)
In the above equation, B˜ is the modified bulk modulus, δ is the thickness of the pipe wall, D is its diameter and E pipe is the elasticity modulus of the pipe material and B f luid is the bulk modulus of the fluid. If the flow is assumed to be barotropic (i.e., the density depends only on the pressure), the continuity equation can be written as ∂p ∂v + ρ0 a s 2 = 0. ∂t ∂x
(2.49)
The viscosity term on the right-hand side of (2.47) represents the losses in the flow. For brevity, let us denote this term by S. Finally, we end up with a system of linear partial differential equation in the form 0 ∂ p 0 a∂ x p , = + ∂t ρav a∂ x 0 ρav ρaS
(2.50)
where ∂ x represents the differentiation with respect to x. Also, for the sake of simplicity, a denotes the sonic velocity e.g. a s in (2.49) and ρ stands for the constant fluid density, i.e. ρ0 in (2.49). The above form will be useful later on.
2.5.2 Friction modelling Traditionally, steady or quasi-steady friction terms are incorporated in hydraulic transient softwares but experimental results show that the steady-state friction model fails to predict both the attenuation and phase of the wall shear stress (and thus the pressure and velocity histories) in the case rapid hydraulic transients. These discrepancies are due to the deviations of the velocity profile; even if the flow remains laminar under the transient process, its slope at the walls may oscillate in a wide range. This has led to the development of unsteady friction models. One of the most earliest models is developed by Daily et al. [24] in which the friction factor depends on the instantaneous mean velocity and instantaneous local acceleration. Brunone in [18] deduced an improved version, in which the convective acceleration is also taken into account and gives a relatively good match between the measured and computed values. However, the main disadvantage of the Brunone model is that the key parameter k (Brunone friction coefficient) has to determined experimentally, i.e. by trial and error. Zielke in [90] derived a model for frequency dependent friction for unsteady laminar flow. The wall shear is related to the instantaneous mean velocity and to the weighted sum of the past velocity values. The advantage of this approach is that no empirical coefficients are included, however, the numerical evaluation of the friction term is extremely expensive from a computational point of view. Several attempts were made to improve the cost of Zielke’s model and to develop weights for turbulent flow (e.g. [84] or [85]) but its main drawback is still its high computational cost. Beside the one-dimensional (1D) models, which in different ways indirectly incorporate the instantaneous cross-sectional velocity profile, several 2D models have been proposed including e.g. [89] and [83]. The interested reader should consult [13], where an excellent overview is given. After the literature survey it is concluded that the two competitive model for 1D unsteady friction modelling are the Brunone-Vítkovsky model [13] and Zielke’s model [90]. In [13] numerical and experimental comparison is given and it is shown that - although the difference is not significant - the Brunone-Vítkovsky model gives closer values to the measured results. Unfortunately, this excellent agreement is due to the possibility of ’software tuning’ by Brunone’s coefficient k. Also, in [19] the authors report a numerically unstable solution when using an explicit finite-difference form of the Brunone model together with the method of characteristics. In this section, we present the quasi-steady model, the Brunone model and the Zielke model. In the case of the Brunone model, such a form is presented, which can be handled by the numerical techniques presented in §2.5.3 in a convenient and stable way. It is also shown how Zielke’s model can be implemented in a numerically inexpensive way. 28
2.5. Transmission line modelling
Quasi-steady friction modelling λ Traditionally, the steady laminar friction is modelled as S s = − 2D v |v|, where λ is either a constant number (0.02 for crude
computations) or it is computed from λ = ∂ ∂t The Brunone-Vítkovsky model
64 Re .
Inserting this into (2.50) gives
p 0 + ρav a∂ x
a∂ x β
p = 0 , 0 ρav
where β =
32ν . D2
(2.51)
The original Brunone model relates the unsteady friction part Su to the instantaneous local acceleration ∂v/∂t and instantaneous convective acceleration ∂v/∂x. A modified version proposed by Vítkovsky in [13] should be used for the general case λu =
! kD ∂v ∂v + a sign(v) , v|v| ∂t ∂x
in which k is Brunone’s friction coefficient defined as √ 0.0476 C∗ where C ∗ = k= 2 7.41 Re− log(14.3 Re−0.05 )
(2.52)
for laminar flow and
(2.53)
for turbulent flow.
For later numerical purposes, it is convenient to express (2.52) as " ! # ∂(ρav) ∂v ∂(ρav) de f. k ∂(ρav) ∂(ρav) λu = −α + a sign(v) sign −γ . ρaSu = − ρav|v| = − 2D 2 ∂t ∂x ∂x ∂t ∂x
(2.54)
The two sign terms in (2.54) make the relationship cumbersome to handle. If they were known, coefficient γ would also be known a priori. The key step and the main simplification is that the coefficient γ is computed using the values from the previous time step. The benefit of this approximation is that in this case, (2.50) remains linear and takes the following form ∂ ∂t The Zielke model
0 p + (1 + α) ρav a∂ x
0 a∂ x p = . 0 β + γ∂ x ρav
(2.55)
Zielke’s friction model was developed for laminar flow conditions. By taking the Laplace transform of the momentum and continuity equation, the velocity distribution can be computed for arbitrary boundary conditions (for details, see [90] or [38]). In the Laplacian domain the wall shear τˆ 0 turns out to be τˆ 0 (s) =
ρR s vˆ (s), q J1 i νs R − 2
where ρ denotes the density of the fluid, R is the inner radius of the pipe, i = viscosity of the fluid and J1 (z) =
z JJ01 (z) (z)
(2.56) √
−1 is the complex unity, ν is the kinematic
denotes the modified Bessel function (J0 (z) and J1 (z) stand for the Bessel function
of the zeroth and first order). To calculate the friction factor, we need the inverse Laplace transformation of (2.56). The condition of the existence of such a function requires that τˆ 0 (s) tends to zero while the real part of s tends to infinity (see e.g. [45]). However, this condition holds only if we associate the s to vˆ (s), which physically means that not the mean velocity vˆ (s), but its derivative with respect to time s vˆ (s) determines the friction. The inverse is then formally given by
τ0 (t) =
Zt
φ(t − u)
∂v (u)du, ∂t
ˆ = where φ(s)
0
29
ρR . q J1 i νs R − 2
(2.57)
2. Experimental set-up and mathematical modelling of the system elements
Thus we see that the wall shear is the weighted integral of the past acceleration values
∂v ∂t (u)
with weight function φ(t − u),
and the integration is performed from the beginning of the computation u = 0 to the present time step u = t. Calculating ˆ means the evaluation of the following integral: the inverse Laplacian function of φ(s) σ+iω Z
1 φ(t) = lim ω→∞ 2πi
φ(s)e st ds.
(2.58)
σ−iω
ˆ It holds for φ(s), that F(s) → 0 ha |s| → ∞ and the function is analytical over the complex plane, except the origin and q s the series of poles along the real axis z j = i νj R, for which we have J1 z j = 2. The location of these poles is easy to
compute (numerically) if one notices that the distance between them tends to π. The first ten is given in Table 2.5. z1
5.1356
z6
21.117
z2
8.4172
z7
24.270
z3
11.620
z8
27.421
z4
14.796
z9
30.569
z5
17.960
z10
33.717
Table 2.5: The first ten poles of function
1 . J1 (z j )−2
By virtue of the Jordan-lemma, instead of evaluating directly the integral (2.58), one can use the Residual Theorem (the residuals are given by the −1th coefficient of the corresponding Laurent-series). The origin has to be handled separately, as it is a multiple pole. The details are given in [90] or [38], the result is φ (τ) = where T =
ν t R2
∞
4νρ 2νρ X −z2j τ 4νρ 2νρ + e = + W(T ), R R j=1 R R
(2.59)
denotes the dimensionless time (in [90] τ is used instead of T but this notation might be confused with the
wall shear τ0 ). The wall shear is then given by τ0 (t) =
Zt
4νρ 2νρ ∂v v(t) + φ(t − u) (u)du = ∂t R R
Zt
W(t − u)
∂v (u)du. ∂t
(2.60)
0
0
Note the nice result that if we compute the friction coefficient, we have Zt 8τ0 64 1 ∂v 1 + λ= 2 = W(t − u) (u)du , ρv Re 2v ∂t
(2.61)
0
which means that the well-known quasi-steady friction term is regained. Let us turn now to the implementation issues. First, let us assume that only the first n poles are taken into account, i.e. z j , j = 1 . . . n. In this case, by virtue of (2.59), the second ’unsteady’ part in (2.60) can be written as I(t) =
Zt 0
n
X ∂v W(t − u) (u)du = ∂t j=1
Zt
e−z j
2
(t−u) ∂v
∂t
0
n de f. X
(u)du =
I j (t).
(2.62)
j=1
Let us assume now that in the current step of the computation, the quantities I j (t) are already known and I j (t + ∆t) is to be computed in the next integration step. The integral can be evaluated in two parts: from 0 to t and the from t to t + ∆t:
I j (t + ∆t) =
Zt 0
e
t+∆t Z 2 ∂v e−z j (t+∆t−u) (u)du. (u)du + ∂t ∂t
−z j 2 (t+∆t−u) ∂v
t
30
(2.63)
2.5. Transmission line modelling
In the second term, if ∆t is ’small enough’ (compared to t), one can use the integral mean value theorem and replace the integral by the value at the middle of the interval u = t +
∆t 2
multiplied with the length of the interval:
t+∆t ! Z ∆t ∂v ∂v ∆t 2 2 ∆t 2 e−z j (t+∆t−u) (u)du ≈ e−z j (t+∆t−(t+ 2 )) t+ ∆t ≈ e−z j 2 (v (t + ∆t) − v(t)) . ∂t ∂t 2
(2.64)
t
The first term in (2.63) is easier to evaluate by exploiting that e−z j I j (t + ∆t) = e−z j
2
∆t
I j (t) + e−z j
2
(t+∆t−u)
2 ∆t 2
= e−z j
2
∆t −z j 2 (t−u)
e
. The result is
[v (t + ∆t) − v(t)] .
(2.65)
In the actual calculation, the I j values are stored separately at every grid point and updated after each time step. This means that the wall shear has a delay of one time step, which can be usually acceptable (especially if the temporal resolution is fine). The benefit is that updating the I j values is quick as only the algebraic manipulation defined by (2.65) has to be performed. As the I j values are updated independently from the continuity and momentum equation, (2.50) turns now into ∂ ∂t
p 0 + ρav a∂ x
a∂ x β
p = 0 , ρav δ
where β =
32ν D2
n
and δ =
16ν X I j, D2 j=1
(2.66)
where δ is some constant recalculated after each time step with the help of (2.65) and (2.62). However, if this technique is not acceptable because e.g. a more accurate evaluation of I(t) is needed, the approach of Kagawa et al. [41] can be used, where by differentiating (2.60), a new set of ordinary differential equations arise for
dτ0 dt ,
which then has to be integrated
simultaneously with the continuity and the momentum equation.
2.5.3 Numerical techniques for transmission line modelling In what follows, several numerical techniques are given to solve the set of partial differential equations (2.50). First of all, the well-known method of characteristics (MOC) is briefly described and it is pointed out that there is a rigid connection between the spatial and temporal discretisation, which is a severe limitation when trying to connect the pipe to a variable step size ODE solver. To detach the temporal and spatial resolution, the method of lines technique is used, i.e. first, discretising the equations in space and then integrate the resulting system of ordinary differential equations. Several spatial schemes and time marching techniques are compared (in pairs) in terms of accuracy and computational effort with the help of two analytical test cases. Method of characteristics Let µ be an unknown multiplicator, L1 and L2 stand for the continuity and momentum equation given by (2.50), then ! ∂ µ ∂ 2 (2.67) ρa v + p = µS(x, t). L1 + µL2 = (p + µv) + ∂t ∂x ρ If µ = ±ρa, the above equation can be rewritten along the two characteristic lines C + : dx/dt = a and C− : dx/dt = −a as ordinary differential equations (for details, see Halász et al. [34]): C+ : C− :
d (p + ρav) = S (x, t) and dt d (p − ρav) = − S (x, t). dt
(2.68) (2.69)
Referring to Fig.2.22, we introduce the notations p L(e f t) = pnj−1 , pC(entre) = pn+1 and pR(ight) = pnj+1 and similar ones j λ for velocity v. Note that the left-hand sides of (2.69) can be integrated explicitely, while the friction term S(x, t) = − 2D v|v|
31
2. Experimental set-up and mathematical modelling of the system elements
C tn+1 C−
C+ tn
x1
x2
L x j−1 ∆x
∆t MOC ∆t MOC /2
R xj
∆x
x j+1
xN
xN+1
Figure 2.22: Computational grid for the method of characteristics.
has to be integrated approximately: C+ : C− :
λ∆t vC |vL | and 2D λ∆t pC − pR − ρa (vC − vR )=ρa vC |vR |. 2D
pC − pL + ρa (vC − vL )= − ρa
(2.70) (2.71)
After rearranging, we end up with pL − pR + ρa (vL + vR ) i h vC = ∆t (|vL | + |vR |) ρa 2 + λ2D
and
pC = pL + ρavL − ρavC
! λ ∆t 1+ |vL | . 2D
(2.72)
For the first and last grid point, we have
λ∆t n n v |v | and 2D 2 2 λ∆t n n = {p + ρav}nN − ρa v |v |. 2D N N
{p − ρav}n+1 = {p − ρav}n2 + ρa 1
(2.73)
{p + ρav}n+1 N+1
(2.74)
Once pressure, velocity or some algebraic relationship between them is given, the values at the next time level can be computed. The strength of MOC is the accurate and quick integration of (2.50) but the price for this is the constant and rigid timestep. Now let us assume that he pipe is connected to a hydraulic cylinder. The flow rate from the pipe inside the cylinder’s chamber is a parameter of the cylinder’s equation of motion, while the chamber pressure is a boundary value for the pipe solver. However, these values are needed at the same time level, which is cumbersome to achieve as the MOC marches with a fixed ∆t while the ODE integrator is likely to demand an another and much smaller time step. Numerical experiments show that the pressure relief valve model described in §2.4 requires a time step of ∆t = 2 × 10 −4 [s] for a stable integration (using a fourth-order explicit Runge-Kutta scheme). Assuming the sonic velocity to be a = 1000[m/s] and L = 2[m] pipe length, the required number of grid points is N = L/a/∆t = 100. An obvious solution is to interpolate between the grid points based on the actual timestep, as shown with red in Fig.2.22 in the case of ∆t desired = ∆t MOC . The consequence of such an approach is depicted in Fig.2.23, where a few periods of a standard water hammer problem were calculated with unit (black lines) and half step size (red lines). Not only the accuracy of the scheme is destroyed but also artifical damping is introduced as the maximum pressure and velocity values decrease continuously. Also, the interpolation itself increased the computational effort by a factor of approx. 5. Method of lines In this section, several spatial and temporal discretisation schemes are presented for solving the continuity and momentum equation together with the friction models given by (2.51) or (2.55). Zielke’s friction model is not given directly in a similar form but (2.50) has to be used in such cases, with an appropriate function recalculating the source term S after each time step by virtue of (2.65). In what follows, only the general description of the procedure is given not showing the steps in detail. 32
2.5. Transmission line modelling
11.5 11
pu [bar]
10.5 10 9.5 ∆t ∆t/2
9 8.5
0
0.01
0.02
0.03
0.04
0.05
0.06
vd [m/s]
0.2 0.15 0.1 0.05 ∆t ∆t/2
0 −0.05
0
0.01
0.02
0.03
0.04
0.05
0.06
t [s]
Figure 2.23: Test case for MOC, unit step size (black lines) and interpolation with half step size (red lines). Parameters values used are as follows. Pipe: L = 10[m], a = 1000[m/s], ∆x = 1[m]. Initial state: p(x, 0) = 10[bar], v(x, 0) = 0[m/s]. Boundary conditions: upstream: vu = 0.1[m/s], downstream: pd = 10[bar].
Let us start with (yet) another form of (2.50) which is useful in the following calculations. For brevity, without restricting the general case, let us neglect the the source term S in (2.50) and rewrite the system with the help of the new
variable y = (p, v)T as
∂y(x, t) ∂y(x, t) +A = B, ∂t ∂x
(2.75)
where, formally speaking, A is a matrix and B is a vector of linear operators. For example, if the Brunone-Vítkovky model is employed, we have (see (2.55)) −1 1 0 0 A = 0 1 + α a∂ x
a∂ x β + γ∂ x
0 and B = . 0
(2.76)
(Let us emphasise again that this is only a general formal notation to explain the basic idea.) The first step is to choose a spatial discretisation scheme, which approximates the spatial derivative at a grid point with the linear combination of some neighbouring grid points, thus it can be represented by a differentiation matrix D. After the spatial discretisation is performed, the boundary conditions are taken into account, which results in (a) completely removing rows of the differentiation matrix and y (corresponding to the equations for the boundary values) and (b) removing columns of the differentiation matrix into a separate vector b. The semi-discrete form of (2.75) can be written as dy(t) = D y(t) + b, ∂t
(2.77)
where y(t) includes the nodal time histories of pressure and velocity. For example, if upstream velocity and downstream pressure are prescribed, y(t) takes the following form: y(t) = p1 (t), p2 (t), p3 (t), . . . , pN−1 (t), pN (t), v2 (t), v3 (t), . . . , vN−1 (t), vN (t), vN+1 (t) . 33
(2.78)
2. Experimental set-up and mathematical modelling of the system elements
(2.77) is a system of ordinary differential equations, which can be integrated in t to give the pressure and velocity distribution in the pipe. The ultimate goal of this section is to choose such spatial and temporal discretisation schemes which provide a quick and stable way of solving (2.50), plus the unsteady friction models. However, there is a fascinating diversity of differentiation and integration schemes, from the simple two-point Euler scheme to the meshless finite element methods or from the classic Runge-Kutta schemes to the modified extended backward differentiation schemes. The following guidelines help to restrict our attention to a few of them. First of all, it should be noted that when applied to a problem, spatial and temporal schemes are not ’working’ independently and the spatial algorithm plays an important role in determining the temporal efficiency. An illustrive example is given in [32] where it is shown that • the error decreases if the order of the spatial approximation is increased, motivating high-order spatial formulations but • the cost of the computations increase with increasing the order of accuracy. A global optimum is reached for a finite value of the order of the temporal scheme pt.s.. [32] suggest 4 ≤ pt.s. ≤ 6 for most practical problems. • The optimal order for spatial and temporal operators is similar and • spatial discretisation error accumulates linearly in time. Let us turn now to the desired features of the spatial scheme. Consider (2.77) with N + 1 grid points on an equidistant grid and assume that the boundary conditions were already set. For example, in the simplest case, prescribed pressure and/or velocity values at the boundaries result in the two nodes excluded on the boundary from y and put into the source vector b. The equilibrium is given by y0 = −D−1 y and if we introduce a new variable z via the linear transformation y = T (z + y0 ), where the columns of T consist of the eigenvectors of D, the modal equation reads ∂z = Λ z, ∂t
j = 1 . . . 2N,
(2.79)
where Λ is a diagonal matrix containing the eigenvalues of D (if D is regular). System (2.79) is a system of decoupled linear ordinary differential equations and if the solution of the original system (2.77) remains bounded (as expected in real engineering applications), the same must hold for (2.79). System (2.79) is stable if and only if T OL, the step was too large, the results are rejected and the step is re-taken with a reduced stepsize, otherwise it is accepted. This step size reducing is repeated until the step is accepted. Then, if err T OL, the stepsize was inefficiently small and it is increased by a maximum factor of two for the next step. We now turn to a special aspect of system simulation concerning to the boundary conditions. As all integration schemes are implicit, they require boundary conditions from the "future", i.e. the value of g(t + ∆t) (for ESDIRK methods, even at intermediate stages). In reality, these values are not available, because these boundary values are supplied by other system elements (e.g. other pipe, reservoir, check valve, piston, etc.) whose state is also not known beyond the actual time level. One way to overcome these difficulties is to store a few old values and extrapolate or, as in our tests, simply consider the boundary values to be constant up to the next time step.
2.5.4 Comparison and results The aim of the benchmark tests presented here is to decide which combination of spatial and temporal schemes is the most efficient and accurate one. If both the steady and unsteady friction terms are neglected, the problem reduces to the solution of the wave equation and we are in the enviable situation of knowing the analytical solution, i.e. f (x − at) for any initial condition f(x,0). Two test functions were chosen, a smooth one f 1 (x, t) = sin (2π/L(x − at)) and another one with a
steep front f2 (x, t) = 1/2 tanh(10(x − at)). The integration lasted for ten complete periods and the Euclidean L 2 norm of the error is used for visualisation of the behaviour. CDS3
3
CDS3
10 −1
10
err
L2
CPU
2
−2
10
1
10
0
10
10
2
CDS5
10
3
2
CDS5
10
10 −1
10
err
L2
CPU
2
−2
10
1
10
0
10
10
2
FEM
10
3
2
FEM
10
10 −1
10
err
L2
CPU
2
−2
10
1
10
0
10
10
2
2
10
10
N
N
Figure 2.25: Smooth test case, f1 (x, t) = sin (2π/L(x − at)). Black circles denote BEu, purple stars stand for BDF2, the blue plus signs and the red crosses denote ESDIRK3 and ESDIRK4, respectively. CPU time is given in seconds. It has to be noted first that the Chebyshev collocation technique forced the ODE integrator to take extremely small steps, which were not acceptable from a practical point of view. We speculate that the reason is the strongly varying grid size, which gives rise to stiffness. The rest of the tests are depicted in Fig.2.25 and Fig.2.26, where black circles denote 39
2. Experimental set-up and mathematical modelling of the system elements
0
CDS3
2
10
CDS3
10
CPU
err
L2
1
−1
10
−2
0
10
−1
10
0
10
10
2
CDS5
10
2
10
2
CDS5
10
10
CPU
err
L2
1
−1
10
−2
0
10
−1
10
0
10
10
2
FEM
10
2
10
2
FEM
10
10
CPU
err
L2
1
−1
10
−2
10
0
10
−1
10
10
2
2
10
10
N
N
Figure 2.26: Steep test case, f2 (x, t) = 1/2 tanh(10(x − at)). Black circles denote BEu, purple stars stand for BDF2, the blue plus signs and the red crosses denote ESDIRK3 and ESDIRK4, respectively. CPU time is given in seconds.
BEu, purple stars stand for BDF2, the blue plus signs and the red crosses denote ESDIRK3 and ESDIRK4, respectively. It is clearly shown that the use of high-order integration schemes is beneficial, especially as the number of grid points increases. As to the spatial schemes, the finite element formulation was found to be the most accurate one. However, for the smooth test function and large number of grid points, the CDS5 scheme is competitive with FEM both in terms of accuracy and CPU time. The CPU time for all of the schemes increases with approximately the same rate and are in the same order of magnitude, for the steep case the extra effort needed by the higher-order schemes is five times the effort of the lower ones. Thus we conclude that it is advantageous to use higher-order temporal and spatial schemes together. Fig.2.27 depicts the tests of unsteady friction models. All the results were computed with FEM technique and ESDIRK4 scheme. The parameters for the unsteady test case were chosen to be L = 40[m], d = 20[mm], a = 1319[m/s], ρ = 1000[kg/m3] and ν = 1.46 × 10−6 [m2 /s]. Initially, the velocity of the fluid was v(x, 0) = 0.1[m/s] with p(x, 0) = 10[bar] pressure. Beginning at t = 0[s], the inlet velocity was linearly decreased to zero in 0.01[s] closing time while the pressure remained constant at the other end of the pipe. All computations were performed with N=50 grid points and the temporal tolerance was set to T OL = 10−3 . In Fig.2.27, black line corresponds to the frictionless case, red line depicts the quasi-steady case and the dashed blue line denotes Zielke’s model. With an appropriate Brunone friction coefficient value k it could be set that, the two unsteady friction models were within line thickness the same, thus the Brunone-Vítkovsky model is not shown. As seen in Fig.2.27, the unsteady friction model generates a significantly larger damping than the quasi-steady one. Also, a slight phase shift appears, which is also reported in the literature. Compared to the frictionless case, the extra computational work for quasi-steady and unsteady friction modelling was below 1% in terms of CPU time. It has to be noted here, that a weakness of these techniques is that excitations with large gradients (e.g. fast valve closing, water hammer problems) produce serious numerical oscillations. It should be emphasises again that the experimental verification of the unsteady friction models is beyond the scope of this work. This is mainly because 40
2.5. Transmission line modelling
Figure 2.27: Comparison of the friction models. Black line corresponds to the frictionless case, red line depicts the quasi-steady case and the dashed blue line denotes Zielke’s model. recently in [13], Bergant et al. conducted several experiments to verify the Brunone-Vítkovsky and Zielke models and observed an excellent agreement with theory.
2.5.5 Handling pipe junctions The last issue corresponding to the transmission line modelling is the handling of the pipe junctions. Independently of the actual solution techniques (MOC or FEM with ESDIRK4), the following techniques was used. Let us assume that N pipes meet at a single point. In this case, we know that the pressure must be the same and that the (signed) sum of the mass flow rates has to be zero. For the sake of simplicity, we assume that the direction of every pipe is such that the outgoing volume flow rate is positive. In this case, one can use the characteristic lines to gain information from the inner part of the pipe: p − ρavi = βi
i = 1...N
(2.99)
and, because of the continuity, we know that N X
Ai vi = Qin .
(2.100)
i=1
This is a linear system of equations for unknowns vi , i = 1 . . . N and p. The solution is v1 =
Qin +
PN
Ai i=2 ρai (βi P N Ai a1 i=1 ai
− β1 )
.
Once one velocity (v1 ) is known, it is easy to calculate the pressure and the other velocities. 41
(2.101)
2. Experimental set-up and mathematical modelling of the system elements
2.6 Numerical simulation of hydraulic transients in the experimental rig In this section, the transient response of the experimental rig presented in §2.2.1 is studied both experimentally and numerically. The transient scenarios are sudden load change and piston impingement against the cylinder top, both of which were studied in the case of several pump revolutionary numbers, i.e. supply volume flow rates. Fig.2.28 presents the measured pressure histories and piston displacement for several pump revolution numbers. Referring to Fig.2.1 in §2.2.1, in the four upper panel, blue lines correspond to p 1 , pressure at the input of the relief valve, green line represents p2 , which is the lower chamber pressure and the red line is the upper chamber pressure p 3 . The bottommost panel depicts the displacement of the piston for the last case (n = 1230[rpm]). 100 p [bar]
n=502 [1/min] 50 0
0
0.5
1
1.5
2
2.5
3
3.5
4
100 p [bar]
n=800 [1/min] 50 0
4
4.5
5
5.5
6
6.5
7
7.5
8
8.5
3
3.5
4
4.5
100 p [bar]
n=1035 [1/min] 50 0
0
0.5
1
1.5
2
2.5
100 p [bar]
n=1230 [1/min] 50
x [mm]
0 0.5 300
1
1.5
2
2.5
3
3.5
4
3
3.5
4
n=1230 [1/min]
200 100 0 0.5
1
1.5
2
2.5 t [s]
Figure 2.28: Measured pressure histories and piston displacement for several supply flow rates (i.e. pump revolution number). Note that the time scales are different and that the piston displacement at the bottommost panel corresponds to the pressure histories shown in the fourth panel. Referring to Fig.2.1 p 1 (relief valve input pressure) is denoted by blue line, lower chamber pressure p2 is the green line and upper chamber pressure p3 is the red curve. See text for further details.
Let us briefly describe the measurements. Initially, the piston is at rest as the cylinder control valve (the threeposition valve ’5’ in Fig.2.1) is in the middle ’hold’ position. As the gear pump is running with a prescribed constant angular velocity, the flow returns to the tank through pipes 5 and 6 and p 1 increases slightly due to the inner friction losses. Note that the noise of the gear pump (pressure pulsation) occurs only at p 1 as p2 and p3 are mechanically separated from the gear pump. The mass held by the piston is 120[kg]. As the cylinder control valve ’5’ is switched to the ’up’ position, 42
2.6. Numerical simulation of hydraulic transients in the experimental rig
the piston moves upwards, pipe 1 and pipe 2 are connected and pipe 3 is connected to the return tank. The pressure difference p1 − p2 is due to the pipe friction and the losses in other elements, e.g. the filter or the control valve. As the piston reaches the position x = 138[mm], additional load appeared, resulting in a total load of m = 120 + 183 = 303[kg]. The upward motion continues until the cylinder control valve is switched back to the ’hold’ position. In the middle range the piston is again at rest with the increased load, the flow is discharged back to tank. The piston is then restarted by switching its control valve ’5’ into the ’up’ position. The control valve is held in the ’up’ position even when the piston hits the cylinder top. After the impingement, the relief valve opens and the flow is bypassed to the tank while the pressure is kept at the preset value. If the pressure in the upper chamber vanishes (which is actually the case) and the friction forces inside the cylinder are neglected, the required pressure to lift the load is p2,req. =
mg 302[kg] × 9.81[m/s2] = 25.6[bar]. = π 2 2 −6 Alower 4 (50 − 32 ) × 10
(2.102)
As the relief valve was set to 67[bar], there was no danger of load releasing. Note that although the piston was at rest in the middle range (say, t = 2[s] in the case of n = 502[1/min]), the pressure difference p 2 − p3 deviates in the four measurements, ranging between 27[bar] and 32[bar]. This difference is much smaller at the end of the resting period (say, t = 3[s] in the case of n = 502[1/min]), which is due to the pressure equalisation through leakage flow. However, a deviation of approx. 3[bar] remains, which is probably due to the indefinite friction force at zero velocity (see Fig.2.4 in §2.3). Let us turn now to the issues of numerical simulation. The first question was whether to implement the models in a system simulation software package or to start developing an in-house code. On the one hand, commercial ’boxed’ softwares (e.g. AmeSim, ITI-Sim, Modelica) are shipped with built-in development tools and thus one implements only the element itself, which is then easily fit into the general software framework. However, if the element does not ’fit’ into this framework, it may be cumbersome to handle. For example, during his fourth-year project, the author was given the task to implement a distributed parameter pipe model in AmeSim. As AmeSim is developed mainly for simulating lumped parameter systems, the implementation process slow and cumbersome. On the other hand, in-house codes have to be built up from the very beginning with significant extra work but give high flexibility to the developer, which is especially important for a scientific work. Based on these considerations (and taking into account the financial aspects) an in-house code was developed. The mathematical models presented in this chapter were implemented in object-oriented Matlab environment. Although a Matlab code runs slower than a compiled code (C++ or Fortran) by a factor of 5 to 10 (according to the author’s experience), one can significantly reduce the coding effort by using Matlab’s built-in high level mathematical, graphical and data I/O functions. Moreover, taking into account the speed of developments of the hardware side, the advantage of compiled codes in speed is equivalent to, say, two or three years of hardware (notably CPU clock) development. The model consisted of six distributed-parameter pipes (those ones, whose length is depicted in Fig.2.1), the piston model presented in §2.3 and the relief valve model presented in §2.4. The only issue not given yet is the modelling of the gear pump. Based on [46], the mean volume flow rate of a gear pump can be calculated as ! π2 2 2 cos α = 0.0248[l] × n[1/min], qm = B 2πn m z + 1 − 12
(2.103)
where B = 63[mm] is the width of the gear, m = 2.6[mm] (estimated value) stands for the modul, z = 9 is the number of teeth, α = 20o is the pressure angle (estimated value) and n is the revolution number of the gear. The gear pump was disassembled and the geometry was measured, which resulted in the numerical value given in (2.103). Note that the catalogue gives Q[l/min] = 0.022[l] × n[1/min] and if the volume flow rate is recalculated from the velocity of the piston, 43
2. Experimental set-up and mathematical modelling of the system elements
the flow rate at n = 1230[1/min] is Q = 26.7[l/min]. Thus, the theoretical mean flow rate gives a reasonable accuracy and the estimated volumetric efficiency is 87.5%, which is realistic. However, to take into account the flow pulsation, i.e. the excitation of the gear pump, the instantaneous flow rate should be known. Prof. L. Kovács in his lecture notes [46] gives Q pump (t) = B 2πn R2 − r2 − x(t)2 ,
(2.104)
where R is the outer radius of the gear, r is the inner radius (in contact) and x(t) is the distance of the contact point from the main point (the two gears are identical). It can be shown that the velocity of the contact point is constant, given by v x = ωra (for details, see [46] or [86] p.30). Inserting this into (2.104) gives the flow rate as a function of time, which can be (numerically) expanded into Fourier series and implemented in the code as a velocity boundary condition for the pipe connected to the gear pump. Figure series 2.29 and 2.30 present the measured and calculated pressure and displacement histories for n = 1230[rpm] pump revolution number.
p1 [bar]
100
50
p2 [bar]
0 0.8 100
p3 [bar]
(c)
(d)
(e)
0.9
1
1.1
1.2
1.3
0.9
1
1.1
1.2
1.3
0.9
1
1.1
1.2
1.3
0.9
1
1.1 t [s]
1.2
1.3
50
0 0.8 200
x [mm]
(b)
50
0 0.8 100
PSfrag replacements
(a)
100
0 0.8
Figure 2.29: Measured and calculated pressure and position histories for sudden load change. The period of time while the pressure relief is opened is marked by the red arrows. The beginnings and ends of periods (a)-(e) are not well-defined thus the corresponding arrows are only illustrative.
Fig.2.29 presents the transients due to piston start and sudden load change. Green lines denote measured data and the blue lines the calculated ones. The very first period denoted by (a) is merely used for calculating the steady-state flow rate and pressure p1 (p2 and p3 are known a priori from force balance). In the second period denoted by (b), the piston is at rest and only the initial mass of 120[kg] loads the hydraulic piston. At the moment of starting approx. t = 0.976[s], the cylinder control valve connects the low-pressure pipe 2 at the gear pump and relief valve to the high-pressure pipe 3 at the lower chamber pressure (see Fig.2.1). It should be noted that this valve opening is assumed to be infinitely fast, i.e. full opening in one numerical time step. However, it is clearly seen in the diagram of p 1 that in reality this opening 44
2.6. Numerical simulation of hydraulic transients in the experimental rig
is slower. At the end of period (c), the starting transients die and apart from the difference of approx. 4[bars] in p 2 , the measured and calculated pressure values and the displacement match. At approx. t = 1.08[s], additional masses of 183[kg] bump against the piston and the load suddenly increases. At the moment of contact, not only the mass increases in the left-hand side of the Newtonian equation of motion of the piston (2.10) (m x¨ = . . . ), but the additional mass also has to be accelerated and this represents an additional force demand. This d’Alembert-type force was estimated from the time needed for the additional masses to reach the velocity of the moving piston. Based on the experiments, it was assumed that after the contact, ∆t = 0.003[s] was needed to reach the lifting velocity v = 0.38[m/s] of the previous period. Due to this large force, the pressure before the piston (p1 and p2 ) increases and opens the relief valve, the period of time while the relief valve opens is denoted by red arrows in the figure. After the relief valve opens, the pressure drops quickly and thus the valve closes again. During the rest of period (d) the (increased) mass moves with the same velocity as before the bump. After the transients die in the second part of period (d), we observe again a reasonably good agreement between the measured and calculated pressure histories. Finally, the control valve is switched back to the middle ’hold’ position and the the piston stops. The transient pressure histories are captured poorly in this region (beginning of period (e)), which is due to the insufficient modelling of the control valve closure and fluid column separation in pipes 5 and 6.
p1 [bar]
100
50
0 2.4 100
(a)
(b)
(c)
(d)
2.5
2.6
2.7
2.8
2.9
3
3.1
3.2
2.5
2.6
2.7
2.8
2.9
3
3.1
3.2
2.5
2.6
2.7
2.8
2.9
3
3.1
3.2
2.5
2.6
2.7
2.8 t [s]
2.9
3
3.1
3.2
p2 [bar]
80 60 40
20
3
p [bar]
2.4 30
10
0 2.4 300
x [mm]
PSfrag replacements
200
100 2.4
Figure 2.30: Measured and calculated pressure and position histories in the case of piston impingement against the cylinder top. The period of time while the pressure relief is opened is marked by the red arrows (see text for details). The beginnings and ends of periods (a)-(d) are not well-defined thus the corresponding arrows are only illustrative.
In Fig.2.30 the transients due to piston impingement against the cylinder stop are presented. Green lines denote measured data and blue lines stand for the calculated values. The very first period denoted by (a) is merely used for calculating the steady-state flow rate and pressure p1 (p2 and p3 are known a priori from force balance). In the beginning period (denoted by (b) in the Figure), the piston is at rest, the area weighted difference of p 2 and p3 equals the gravity 45
2. Experimental set-up and mathematical modelling of the system elements
force acting on the mass. The fluid delivered by the gear pump returns to the tank via the cylinder control valve (element ’5’ in Fig.2.1) and the pressure difference p1 − p0 equals the overall friction losses (pipeline plus elements e.g. filter, valves). The loss coefficients are estimated based on the catalogue data. At t ≈ 2.6[s], the cylinder control valve is instantly switched to the ’up’ position. There is a slight depression in p 2 as pipe 2 is connected to pipe 1, followed by a pressure peak of approx. 70[bar]. At the same moment, pipe 4 is connected to pipe 5 and after a depression wave, the fluid returns to the tank via the cylinder control valve. The losses in this return section (pipe 3,4,5 and control valve 5) are slightly overestimated, see (c) in p3 . In interval (c), the piston moves upwards with constant velocity (v = 0.38[m/s], defined by the gear pump flow rate and lower piston area). The pressure force due to the lower chamber pressure p 2 equals the sum of the upper chamber pressure force, the friction force plus the gravity. Finally, the difference between p 1 and p2 is simply the overall pressure loss due to the elements between the gear pump and the cylinder (pipes, filter and valves). At t ≈ 2.94[s], the piston hits the cylinder stop and suddenly stops. This induces a large pressure peak in p 2 and a depression in p3 , which drops below the ambient pressure, presumably the fluid column separates as there is no wave reflection. After the impingement, the cylinder control valve remains in the same ’up’ position thus the pressure relief valve opens and the fluid is bypassed to the tank. This interval is denoted by (d). In this interval, p 1 equals the set pressure of the relief valve. Theoretically, in such cases the fluid between the cylinder and the pressure relief valve is at rest, thus p1 should equal p2 . However, measurements show that there is a slight difference between these pressures, which is probably due to the leakage flow at the piston sealings, resulting in a flow rate from the gear pump towards the piston. Finally, it should be noted that the pressure transients are captured with a limited accuracy; the peaks are usually overestimated while the attenuation and the reflections are poorly predicted. However, it should be noted that up to the author’s experiences, these pressure histories were uncertain even in the measurements; when the test case was repeated, only a qualitative agreement was found with the previous ones and notably during the transient periods, the actual pressure histories deviated significantly.
2.7 Summary and conclusions Numerical simulation allows the engineers to answer ’what if ?’-type questions already in the phase of design, thus shortening the time and price of development. Also, extreme operating conditions can be studied in a quick, easy and safe way. To this end, realistic mathematical models of the elements are needed, which can be especially challenging to formulate for transient purposes. In the case of the hydraulic system studied in this work, the transmission line modelling and the dynamic operation of the pilot operated pressure relief valve were found to be the two issues where the models described in the literature needed further development. The following studies were performed and conclusions were drawn. 1. A mathematical model has been presented, which is capable of capturing the dynamic behaviour of the pilotoperated pressure relief valve. The model consists of two Newtonian equations of motion for the main stage and the pilot stage. The main challenge was to give an estimation of the fluid forces on the pilot valve body, which is a conical body with sharp, chamfered-type seat. It was found that the pressure distribution around the valve body is hard to describe analytically and thus we experienced a considerable discrepancy between the estimated fluid forces and those ones calculated with CFD. However, the inaccuracies in the pilot stage model do not lead necessarily to the malfunction of the whole model: (a) The purpose of the pilot stage is simply to open the main stage, which, due to its regular geometry, is relatively easy to describe mathematically. The opening pressure of the pilot stage can be accurately set by the pilot 46
2.7. Summary and conclusions
spring precompression parameter and once the main stage opens, the overall behaviour of the relief valve (notably pressure difference - flow rate characteristics) is dominated by the main stage. (b) Moreover, the main stage is largely overdamped (by Orifice 3, see Fig.2.8). Thus, even if the dynamic behaviour of the pilot stage is not captured accurately, the dynamics of the main stage is dominated by this strong damping effect. However, it was also found that the pilot stage is poorly damped and this may lead to self-excited oscillations, sometimes referred to as valve chatter. An evident question is how to increase the damping on the pilot stage in an easy and inexpensive way. A possible way is given in §3 and a theoretical analysis of valve chatter is given in §4. 2. Transmission line modelling is essential in capturing the pressure waves (e.g. highest pressure peaks) in the case of transient events. However, the traditional way of friction modelling and numerical solution of the fluid mechanical equations may not yield a satisfactory result. (a) It can be hard to achieve compatibility between the large and fixed temporal stepsize of the method of characteristics (MOC) and the small and varying stepsize of ordinary differential equation (ODE) solvers (solving e.g. the ODEs of the relief valve or the piston). In §2.5.3 it was shown how a finite element method together with a high-order implicit Runge-Kutta scheme can be employed for accurately calculating the transient in a transmission line. The technique requires little extra computational effort (compared to the MOC) and adjusts its timestep automatically to the accuracy demands of the user. (b) It was also shown in §2.5.2 how the two most widespread unsteady friction models available in the literature can be implemented in a quick and stable way. Originally, Zielke’s model requires the bookkeeping of several nodal velocity values from the past and it is computationally expensive to evaluate the friction term. By employing some simple properties of the exponential function, only a few extra quantities have to be stored and also the evaluation of the friction term becomes quick. In the case of the Brunone-Vítkovsky model the main difficulties are the numerical instability reported in the literature and the efficient implementation. With the help of one minor approximation (the spatial derivatives in (2.52) are computed on the previous time level) it can be fully incorporated into the variable stepsize FEM-ESDIRK framework and the implicit time integration guarantees the stability of the scheme. 3. Finally, in §2.6 the transient behaviour of a simple linear actuator circuit was numerically computed and the results were compared to the measurements. In general, it can be concluded that piston displacements can be accurately predicted while pressure transients are hard to capture. (a) The proper description of valve openings and closings is important. Fast valve opening induces larger pressure peaks which disappear slowly. (b) Negative relative pressure, fluid column separation and cavitation need further study as neglecting them results in unrealistic pressure wave reflections. (c) The pressure losses in minor elements (elbows, T-junctions, orifices) need also further studies; the loss coefficients (ζ) ’tuned ’ for one measurement may not give a satisfactory agreement for another test case. (d) In the case of sudden load change, the force needed for accelerating the additional load (d’Alambert-type force) was estimated by means of trial-and-error based on the measurements. It would be important to describe this force theoretically for general modelling purposes.
47
2. Experimental set-up and mathematical modelling of the system elements
Issues in which the numerical modelling was successful: (a) It was shown that the models described in this section are suitable for transient simulation. (b) The dynamic model of the relief valve was found to be suitable for dynamic simulation (in both test cases, the realistic dynamic performance of the relief valve model was essential). (c) Although the attenuation of the transmission line model is still smaller in the numerical simulations than in the measurements, the unsteady friction models give much more realistic results. (Simulations with quasi-steady models were also performed but are not presented here.)
48
Chapter 3
Analysis of the flow around a poppet valve in the case of long conical seat 3.1 Introduction Long, flat poppet valve seats are desirable from several points of view. They are dynamically more stable, reliable from a leakage point of view, they help minimise permanent marks on the poppet valve. An overview of the corresponding literature can be found in [12]. Hereby we summarise the main results of the literature corresponding to poppet valve design and analysis. Stone [75] in 1960 studied experimentally the discharge coefficients and steady flow forces on a poppet valve. He assumed linear pressure drop along the cone seat and used the equation to determine the theoretical force acting on the valve body. He found a low-pressure region underneath the poppet valve and he realized that in this region there was a release of dissolved air from the oil. He also found that if the valve exhausts to a chamber filled with fluid, it tended to oscillate at a frequency of 1300[Hz]. He deduced that discharge coefficients needed to be better understood and effects of viscosity, cavitation, fluid rotation needed to be further studied. Takenaka [80] in 1964 studied the performance of several hydraulic control valves, one of them being the pilot operated pressure relief valve. Based on his experiments, he gave formulae for the discharge coefficient for several ranges of the Reynolds number. He also gave an empirical formula for the thrust coefficient based on one particular valve configuration. He established dynamic equations for unsteady behaviour predictions, but unfortunately no experimental validation was given. Urata [82] developed a high-level mathematical approach and studied poppet valves both with and without laps under laminar and turbulent flow conditions. His results were compared to experiments. The basic idea was to assume that the velocity distribution in the gap follows boundary layer theory. With these velocity profiles, integrating the NavierStokes and continuity equations gave the pressure distribution and thus the theoretical force. He found that when the valves have long seats, the pressure distribution on it has to be taken into account in order to have an accurate estimation of force. He found slightly negative pressure at the inlet of the seat and thus deduced that cavitation is likely to occur. Bergada and Watton [12] presented a detailed analysis of poppet valves with long conical seats and developed a set of equations describing the flow rate, the pressure distribution and the thrust on the poppet valve. The novelty of their work is the assumption of laminar flow around the valve body and the explicit relationships obtained for flow rate, pressure distribution and fluid forces. Their new set of equations was compared to CFD and experimental results and an 49
3. Analysis of the flow around a poppet valve in the case of long conical seat
excellent agreement was found. The authors conclude that the pressure relief valve design barely changed for more than thirty years despite the instabilities and other problems related to the body and seat shape.
3.1.1 Outline of this chapter As shown in the literature overview, little theoretical work was done about the flow around poppet valves with long, conical seat. In particular it is not known how to estimate the forces on the valve body under dynamic conditions. The study of Bergada and Watton [12] assumed that the valve body is at rest. In this section an attempt is made to take into account the effect of moving valve body, resulting in formulae for calculating the fluid forces under dynamic conditions. The main goal of this chapter is to calculate the fluid forces on the valve body, i.e. the force due to the pressure and the wall shear distribution. In §3.2, the basic pressure-flow rate relationships are derived, resulting in a closed-formed pressure distribution along the conical gap with moving poppet valve. In §3.3, this exact pressure distribution is integrated along the generatrix of the cone and thus the axial pressure force is obtained. As the resulting formulae are rather lengthy, simplification possibilities are considered. In §3.4, the damping force due to the wall shear is studied, which results in a damping force. From the viewpoint of the dynamic behaviour and stability, damping forces are important. Three different approaches for calculating the viscous force are presented and their results are compared. The analytical solutions are compared to CFD results, similar to those described in §2.4. The actual geometry studied here is exactly the same as that one in §2.4, apart from the valve seat, which is 3.2[mm] long with the same angle as the cone. This allows us a direct comparison of the two valve types in §4.
3.2 The laminar flow through a pilot valve with long conical seat Usually the flow through various valves is calculated with the help of a discharge coefficient C d , in the form of Q ∝ p Cd × area × ∆p, which assumes the flow to be turbulent. This is actually true for spool valves and for poppet valves with chamfered-type seats. However, if both the valve body and the seat is cone-typed and the gap between them is small
enough, there is a narrow channel between the housing and the body of the valve with a considerable length, in which laminar flow develops. Hence, the concept of discharge coefficient cannot be employed in this case. Deciding whether the flow in the gap is laminar or turbulent is not as unambiguous as in the case of e.g. pipe flows, but in [82], Urata has found that the transition is approximately at Re ≈ 2500, where Re = Q/πr0 ν with Q being the volume flow rate, r0 denoting the radius of the inlet hole and ν standing for the kinematic viscosity. Valves with laminar flow have several advantages [12]: • improved stability, • much more predictable flow characteristics (pressure and velocity distribution, etc.) and • linear pressure - flow rate relationship, preferable from a control point of view. In poppet valves with the standard chamfered seat turbulent flow develops, contrary to valves where the length of the seat is much larger than the valve opening gap. In what follows, we adopt the results in [12], where Bergada and Watton derive analytical relationships for poppet valves with conical seat assuming laminar velocity profile in the gap. These results are valid mainly for valves with relatively long chamfered seats and large cone angles 2α ≈ 60−120 degrees. As in our case, the cone angle is 2α = 42 degrees and the gap length is also not as long as in [12], we verify the analytical results with numerical CFD computations. The actual geometry with the notations is given in Fig.3.1. 50
3.2. The laminar flow through a pilot valve with long conical seat
PSfrag replacements L Figure 3.1: Geometry of the poppet valve. α = 21o , L = 3.21[mm], r2,in = 2.8[mm], r1,in = r2,in − H cos α, r2,out =
r2,in + L sin α, r1,out = r2,out − H cos α and H = x pl sin α.
In [12], Bergada and Watton assume a Poiseuille velocity profile in the gap between the valve body at rest and the seat. Their calculation is extended for the case of moving valve body. The velocity profile v(h, y) = −
1 dp h h (H − h) + v pl µ dy 2 H
(3.1)
satisfies the boundary conditions v(0, y) = 0 and v(H, y) = v pl . Integrating (3.1) through the gap gives the flow rate Q=
Z
H 0
v(h, y) 2π (r2 − h cos α) dh =
dp(y) f (y) + g(y), dy
(3.2)
where f (y) and g(y) are defined later on. Exploiting the incompressibility and continuity (i.e. Q = const), the above equation can be separated and integrated to give the pressure distribution as p(y) = pin − ∆p
F (y) G(L) + F (y) − G(y) F (L) F (L)
(3.3)
where p(0) = pin , p(L) = pout and ∆p = pin − pout is the pressure difference. After performing the calculations, we obtain H π H3 π r2,in + y sin α − cos α µ 6 µ 2 ! 2 g(y) = v pl Hπ r2,in + y sin α − H cos α 3 Z y 1 6µ F (y) = d˜y = 3 z(y) f (˜y) H π sin α Z0 y µ g(˜y) d˜y = −v pl 2 (6y + z(y)H cot α) G(y) = y) H 0 f (˜ r2,in − H2 cos α z(y) = ln r2,in − H2 cos α + y sin α f (y) = −
(3.4) (3.5) (3.6) (3.7) (3.8)
Fig.3.2 depicts the velocity distribution for several seat lengths and valve velocities. As seen, for small velocities or short seat lengths, the distribution is close to linear and the cone angle does not change this significantly. However, for long seat lengths, the pressure distribution becomes highly nonlinear. The velocity of the valve body v pl changes the pressure distribution but to a much lower degree. A secondary effect is that as the cone angle increases, the deviation of the pressure distribution from linear increases, which is due to the strong change in the flow-through area and thus the velocity in the gap. Obviously, (3.2) omits several issues, e.g. the effect of the fluid jet leaving Orifice 2, the effect of the contracting streamlines at the entrance of the conical channel or simply the deviation from the laminar velocity profile. (3.2) gives the basic relationship between the parameters, which can then be ’tuned’ to give more accurate results. The results of CFD 51
3. Analysis of the flow around a poppet valve in the case of long conical seat
L=5[mm], vpl=0[m/s]
L=5[mm], vpl=5[m/s]
50
50
o
p [bar]
o
20 o 40 o 60 o 80
40 30
o
20 o 40 o 60 o 80
40 30
30
20
20
10
10
10
0
2
4
0
6
0
L=100[mm], v =0[m/s]
2
4
0
6
o
20 o 40 o 60 80o
40 30
30
20
20
10
10
10
0
100
0
50
y [mm]
20 o 40 o 60 80o
40
20
50
6
50 o
20 o 40 o 60 80o
0
4
pl
50 o
0
2
L=100[mm], v =−5[m/s]
pl
50
30
0
L=100[mm], v =5[m/s]
pl
40
20 o 40 o 60 o 80
40
20
0
p [bar]
L=5[mm], vpl=−5[m/s]
50
100
0
0
y [mm]
50
100
y [mm]
Figure 3.2: Pressure distribution along the conical valve seat, described by (3.3) for several valve lengths, valve velocities and cone angles.
simulations are used to refine (3.2). The velocity-dependent part of (3.2) is neglected as in the CFD computations, the valve body was at rest. Neglecting the material properties and the geometric parameters (which are prescribed values) and keeping in mind that laminar flow was assumed (which means that the dependence on the pressure difference has to be linear) we conclude that only the exponent of the valve displacement is a free parameter to be adjusted to the CFD results. As a first step, we set v pl = 0 and expand (3.2) into Taylor series around x pl = 0, and introduce a new, dimensionless parameter e = x pl /r2,in . This gives Q≈
3 r2,in sin4 α π (p2 − pout ) e3 + O x4pl . L sin α µ 6 ln 1 +
(3.9)
r2,in
Next, we want to find a new exponent for e, denoted by β, which results in a better approximation of the CFD results. The new function Ξ(e) contains only that part of (3.9) which depends on e, i.e. Q
de f.
Ξ(e) =
π µ
(p2 −
4 r3 pout ) sin6 α 2,inL sin α ln 1+ r
= eβ .
(3.10)
2,in
The unknown parameter β can be identified by computing log Ξ and log e and then employing linear regression, which gives β = 2.9. Finally we conclude that if the pilot valve is at rest, the flow rate can be computed as ! 3 π r2,in x pl 2.9 de f. sin4 α Q pl x pl , p2 , pout = (p2 − pout ) = R pl (p2 − pout ) µ 6 ln 1 + L sin α r2,in r2,in
52
(3.11)
3.2. The laminar flow through a pilot valve with long conical seat
2.5
−2
10
x=0.1mm x=0.2mm x=0.3mm x=0.4mm x=0.5mm
2
−3
Ξ [−]
Q [l/min]
10
1.5
1
−4
10
0.5
0.025
0.05
0.075
0.1
0.15
0 0
0.2
5
10
15
20
25
30
35
40
p2−pout [bar]
e [−]
Figure 3.3: Function Ξ(e) (defined by (3.10)) as a func-
Figure 3.4: Volume flow rate through the pilot valve
tion of dimensionless pilot valve displacement e; CFD
vs. pressure difference for several valve displacements.
data points are denoted by crosses, the blue line repre-
CFD data points are denoted by crosses, the lines rep-
3
sents the theoretical result e and the dashed green line
resent (3.11).
β
represents the tuned relationship e , with β = 2.8991. Fig.3.3 shows the actual values of log Ξ(log e) computed from CFD results. Note that for every valve displacements e, ten independent points are depicted corresponding to the ten different pressure steps. It is clearly shown that the values are independent of the pressure. However, for x pl = 0.2mm (e = 0.0714) we have a slight deviation in the data points (green crosses), which is probably due to the insufficient error tolerance set in these CFD runs. In Fig.3.4, the volume flow rate vs. the pressure difference is depicted, for the five valve displacements. Note that the flow is clearly laminar. Again we see a deviation from the linear trend in the case of x pl = 0.2mm and also, for the x pl = 0.4mm case, although the trend is linear, there is an error of up to approx. 25%. For the other three displacements, we observe a reasonable agreement between the theory and CFD results. Fig.3.5 depicts a typical pressure distribution and velocity vector plot with the same parameter values as Fig.2.11, i.e. x pl = 0.3mm and pin = 40bar. By comparing the two figures we immediately see that the long valve seat generates higher restriction and thus the pressure distribution around the cone is much more stable, smaller deviation occurs. The difference between short chamfered-type seats and long conical ones is maybe even more striking if we compare the pressure and velocity distribution along the cone body depicted in Figures 2.12 and 2.13 to Figures 3.6 and 3.7, respectively. However, note that the slightly negative pressure at the outlet of the gap still exists and thus this region is a potential source for cavitation.
53
3. Analysis of the flow around a poppet valve in the case of long conical seat
Negative pressure region
Polyline
Figure 3.5: Pressure distribution and velocity vector plot for x pl = 0.3mm and pin = 40bar. Note that the pressure interval between −2 . . . 3[bar] is exaggerated and red regions correspond to pressures between 3 . . . 40[bar].
60
100
0.1mm 0.2mm 0.3mm 0.4mm 0.5mm
50
0.1mm 0.2mm 0.3mm 0.4mm 0.5mm
90 80 70
40
v [m/s]
p [bar]
60
30
20
50 40 30 20
10
10 0
0 50
52
54
56
58
60
62
64
−10
66
x [mm]
50
52
54
56
58
60
62
64
66
x [mm]
Figure 3.6: Pressure distribution along the cone body
Figure 3.7: Velocity distribution close to the cone body
with sharp-edged seat for several valve body displace-
for several valve body displacements.
ments.
54
3.3. Pressure force
3.3 Pressure force Next the axial force on the valve due to the pressure distribution is calculated. We assume uniform pressure distribution pin up to the radius r1,in and similarly, constant pout is assumed beyond r1,out . In between, we have (3.3). The force due to the pressure distribution is given by 2
F p =Fin + F seat − Fout = r1,in πpin + 2π sin α
Z
L 0
p(y) r1,in + y sin α dy − r1,out 2 πpout .
(3.12)
The integral term in (3.12) gives Z L ∆p L p(y) r1,in + y sin α dy = pin L r1,in + sin α − r1,in A + B sin α 2 z(L) 0 µ 6L + Hz(L) cot α r1,in A + B sin α − v pl 2 z(L) H µ µ + v pl 2 r1,in 3L2 + AH cot α + v pl 2 sin α 2L3 + HB cot α H H
(3.13) (3.14)
where
A= B=
Z
Z
L
z(y)dy = L [1 + z(L)] + z(L)w,
(3.15)
0 L
yz(y)dy = 0
o 1n −2Lw + L2 [1 + 2z(L)]) + 2w2 z(L) 4
and w =
r2,in − H2 cos α sin α
(3.16)
and z(y) is defined by (3.8). However, in the case of small angles, linear pressure distribution can be assumed. This assumption leads to a less lengthy expression; 2
F p = r1,in πpin − r1,out
2
# L2 pin + pout + sin α (pin + 2pout ) . πpout + 2π sin α r1,in L 2 6 "
(3.17)
x=0.1mm 150
140
100
100
F [N]
x=0.1mm x=0.2mm x=0.3mm x=0.4mm x=0.5mm
120
F in F out F seat F sum
50
0 80
F [N]
−50
0
5
10
15
20
25
30
35
40
60
x=0.5mm 10 40
F [N]
5 20
0
−20
F in F out F seat F sum
0
−5
0
5
10
15
20
∆ p [bar]
25
30
35
−10
40
0
5
10
15
20
∆ p [bar]
25
30
35
40
Figure 3.8: Pressure force on the pilot cone body as a
Figure 3.9: Pressure force components on the pilot cone
function of pressure difference for several valve displace-
body for (top) x pl = 0.1mm and (bottom) x pl = 0.5mm.
ments. Crosses denote CFD results, solid lines are ac-
Blue line is the force due to the input pressure, red line
cording to (3.12) and dashed lines represent (3.17)
is the force along the seat and green line is the output pressure force. Light blue line represents the sum of the force components.
Fig.3.8 and Fig.3.9 shows the full solution described by (3.12) (solid line), the linear pressure distribution given by (3.17) (dashed line) and the results of the CFD computations for several valve positions (crosses). It is clearly seen that for 55
3. Analysis of the flow around a poppet valve in the case of long conical seat
the actual geometrical configuration (seat length and cone angle) there is only a slight difference between the full integral of the pressure distribution and the force calculated by means of the linear pressure drop assumption. We conclude that when analysing poppet valves with long, conical seat, it is essential to take into account the pressure distribution beneath the valve seat in order to have a realistic estimation of the pressure force.
3.4 Viscous force In this section the aim is to give an estimation of the viscous shear force on the body due to the motion to decide whether it can be neglected or not. Consider the narrow gap between the cone body and the house. Once the velocity distribution in this gap is given, one can calculate the wall shear by integrating it over the surface to obtain the friction force. As the velocity of the cone is a no-slip boundary value for the fluid velocity distribution, it is supposed to have a significant effect on the velocity profile, hence - indirectly - the velocity of the valve body influences the force on the cone. This suggests that a part of the fluid friction force appears in the equation of motion of the valve body as a damping force; F f riction, f luid = f (v pl ). Direct calculation of the wall shear in the case of a moving valve body requires the need for handling moving grid in CFD computation. This is beyond the scope of this work. Instead, we present three analytical approaches for calculating the wall shear and compare them to each other and also at the resting valve body case to steady-state CFD computations.
3.4.1 Estimation 1. - Equivalent cylinder The simplest approach is to replace the cone with a cylinder of radius r 1 and length l, moving in a circular pipe of radius r2 with velocity v pl . The width of the gap between the cylinder and the pipe wall is H = r 2 − r1 . By assuming linear pressure distribution along the gap, the Navier-Stokes equations for the fluid in the gap simplify to ! ∆p 1 ∂ ∂v ∂p =− =µ r , ∂y L r ∂r ∂r
(3.18)
where y is the coordinate along the gap and the pressure drop between the upstream and downstream is ∆p = p out − pin . The general solution of (3.18) is given by v=−
∆pr2 + C1 ln r + C2 . 4lµ
(3.19)
The boundary conditions are v(r1 ) = v pl and v(r2 ) = 0, thus we see that v pl H(H + 2r1 ) ∆p − , 4L µσ σ i ln(r1 + H) ∆p h 2 r1 ln(r1 + H) − (r1 + H)2 ln r1 + v pl , C2 = − 4L µσ σ C1 =
(3.20) (3.21)
where σ = ln rr12 . The wall shear is
τw = µ Thus, the force due to the wall shear is
! ∆p H + 2r1 µ dv = r1 − H − v pl . dh r=r1 2L 2r1 σ r1 σ
H r1 + F = ∆p −πr12 + π 4 σ
Note that for small gaps H ≈ 0, we have
H 2
L − 2πv pl µ . σ
L π F = ∆p r1 (H − 3r1 ) − 2πv pl µr1 . 4 H 56
(3.22)
(3.23)
(3.24)
3.4. Viscous force
This simple model suggests that the viscous force consists of a pressure-dependent part and a velocity-dependent part. When varying the gap width H, these two components behave in the opposite way; as the gap width decreases, the pressure-dependent part of the viscous force vanishes while the velocity-dependent part increases rapidly. This latter phenomenon represents a strong damping force for small valve openings.
3.4.2 Estimation 2. - Poiseuille flow between conical walls Next, we consider the laminar flow through a conical channel of height H, length L, half-cone angle α and inlet cone radius r1,in . The coordinate along the cone is denoted by y and the one across the gap is H, see Fig.3.1. The assumed laminar velocity profile is
!−1 µ dp(y) 1 dp(y) h v(y, h) = − H − h − 2v pl , µ dy 2 H dy
(3.25)
2
∂ v which satisfies the 1D Navier-Stokes equation 0 = − d p(y) dy + µ ∂h2 and the boundary conditions v(y, 0) = 0 and v(y, H) = v pl .
The flow rate is given by Z
H
dp(y) f (y), dy h=0 where (due to the assumption of incompressible flow) Q is a constant. Integrating (3.26) gives Z y Z y d˜y g(˜y) p(y) = pin + QF (y) − G(y), where F (y) = and G(y) = d˜y. f (˜ y ) y) 0 0 f (˜ Q=
v(y, h) 2π (r2 − h cos α) dh = g(y) +
(3.26)
(3.27)
Taking into account that p(L) = pout , we have Q = (G(L) − ∆p)/F (L) and (3.27) turns into p(y) = pin − ∆p
F (y) G(L) + F (y) − G(y). F (L) F (L)
(3.28)
Thus the wall shear is
H dp(y) ∂v(y, h) µ =− − v pl . (3.29) ∂h h=H 2 dy H The first term represents the wall shear due to the pressure gradient while the second part corresponds to the force due to τ(y) = µ
the velocity of the cone body. Z Fv =
L 0
τ(y) 2π r1,in + y sin α dy
! Z L Q g(y) 2πµL L − r1,in + y sin α dy =v pl r1,in + sin α + Hπ H 2 f (y) f (y) 0 ! ! H(L) µπ G(L) = − ∆pHπ r1,in + sin α + v pl 2L rm,in + H(L) − I(L) . F (L) H F (L)
Here f = and
" # 2 g = v pl Hπ r2,in + y sin α − H cos α 3
H3 π H cos α − 2 r2,in + y sin α 12 µ
rm,in z = ln rm,in + y sin α Z y d˜y 6 µ 1 = 3 z F (y) = f (˜ y ) π sin α H 0 Z y g(˜y) G(y) = d˜y = − (6y + Hz cot α) y) 0 f (˜ Z y 6µ y˜ zrm,in + y sin α d˜y = − H(y) = 2 3 f (˜y) H sin απ Z0 y g(˜y) cot α y˜ I(y) = d˜y = − 6y2 + H 2 z cot2 α − 2H r2,in z + y sin α f (˜y) sin α 0 57
(3.30)
3. Analysis of the flow around a poppet valve in the case of long conical seat
After a lengthy calculation, (3.30) simplifies to Fv,exact = ∆p Hπ
! !# " H L L L H L cos α − sin α − v pl πµ cos α + sin α . 8 r1,in + sin α + 6 2 z(L) H 2 2 z(L)
(3.31)
However, if we assume linear pressure distribution along the gap, we have L L L Fv,linear = ∆pHπ r1,in + sin α − v pl πµ 2 r1,in + sin α . 2 H 2
(3.32)
The difference between the exact viscous force defined by (3.31) and that one given by the linear pressure distribution (3.32) is governed by the two factors multiplying the pressure term ∆p Hπ and the velocity term 2v pl πµL/H. The ratio of the two corresponding factors are depicted in Fig.3.10 and Fig.3.11. To be specific, the quantities Kp = Kv pl =
H 2
L z(L)
cos α − r1,in +
4 r1,in +
sin α
L 2 sin α L 2 sin α
and +3
r1,in +
L 2
H 2
(3.33)
cos α +
L z(L)
sin α
sin α
(3.34)
are plotted, respectively. The parameter values are r1,in = 2.8[mm] and L = 3.21[mm], the cone angle α and the valve displacement x are free parameters. Recall that as H = x sin α, for a given cone angle, the x axis is directly proportional to the gap width. The figures reveal that for our case (α = 21o ), the assumption of linear pressure distribution causes an error of less then 10%, but for moderate valve openings x < 0.5[mm], we have less then 3% of error. Again we conclude that the viscous force consists of two parts behaving counter-wise as the gap width decreases, similar to section §3.4.1. 90
90
0.98
80
0.98
8
6
70
1.02
α [degrees]
1
6
1.0
02
1.
8
1.0
40
1.1
1.0 4
0.9
8
50
1.1
1.1
1.1
1.08
1.08
1.08
50
40 1.06
1.06
1.06
1.04
1.04
1.04
30
30
1.0 6
1.0
1
8
20
1.02
20
1.06
10
0.1
0.2
0.3
0.4
0.5
0.6
10
1.04
1.02 0
1.02
1.02
1.02
1.04
0
1.12
1.12
1.12 60
1.04
60
α [degrees]
80
1
1
0.9
0.9
70
0.96
0.96
0.96
1.02 0.7
0.8
0.9
0 1
x [mm]
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x [mm]
Figure 3.10: Viscous force pressure term error due to lin-
Figure 3.11: Viscous force velocity term error due to lin-
ear pressure drop assumption, K p defined by (3.33).
ear pressure drop assumption, Kv pl defined by (3.34).
3.4.3 Estimation 3. - 2D flow in a diverging channel The previous two models assumed the lack of radial velocity component in the gap between the cone body and the housing. In this section, we take into account that the area of the flow-through channel is changing along the axial direction and allow the presence of radial velocity component. Note that the flow-through area between the cone and the housing is a ring-shaped area, with increasing diameter (and with constant width) along the symmetry axis. Hence, the first step when analysing this flow is to transform the Navier-Stokes equations into conical cylindrical coordinates. The equations are 58
3.4. Viscous force
given in standard textbooks of fluid mechanics (e.g. [70]) or mathematics (e.g. [45]) and are rather lengthy and elaborate. However, a simpler approach can be exploited, which does not require this transformation, yet it is capable of capturing the essential features of the flow (i.e. increasing flow-through area and moving boundary). On the right-hand side of Fig.3.12, this alternative geometry is depicted, that is the flow between two parallel plates, in the segment of 2ϕ angle. (Axis z is the normal of the walls.) L
z,H
H
a
r1
Ain
r0
PSfrag replacements
2j
x
a
Ain
H r,x
s
Figure 3.12: (Left) The cone-ring-shaped flow-through area between the pilot poppet valve and the housing. (Right) The divergent channel used for wall shear estimation.
The next question is how to connect the parameters (a, s, ϕ) to the original parameters (L, r 0 , r1 ). For example, by requiring that the inlet and outlet areas are the same, we have ((r0 + H)2 − r02 )π = 2ϕaH and ((r1 + H)2 − r12 )π = 2ϕ(a + s).
Furthermore, forcing the cone mantel and the circle section to be the same, we have (r 0 + r1 )πL = ϕ((a + s)2 − a2 ). The geometry satisfying these conditions is ϕ=π
r1 (r1 + H) − r0 (r0 + H) , L(r0 + r1 )
a=
L (2r0 + H)(r0 + r1 ) , 2 r1 (r1 + H) − r0 (r0 + H)
and
s=L
r0 + r 1 , r0 + H + r 1
(3.35)
which, by assuming that the gap width H is negligible compared to the radial dimensions, turns into a = r 0 / sin α, ϕ = π sin α and s = L. We are now in the position to apply the steady Navier-Stokes equations given in cylindrical coordinates, neglecting the circumferential velocity component, together with the continuity equation: ! " # ∂ 2 vr ∂vr 1 ∂p ∂ 1 ∂ ∂vr (rvr ) + 2 + vz =− +ν vr ∂r ∂z ρ ∂r ∂r r ∂r ∂z " 2 # 2 ∂vz 1 ∂vz ∂ vz ∂vz ∂vz 1 ∂p vr + vz =− +ν + 2 + ∂r ∂z ρ ∂z ∂r2 r ∂r ∂z ∂vz 1 ∂ (rvr ) + . 0= r ∂r ∂z
(3.36) (3.37) (3.38)
By introducing the dimensionless variables r ξ˜ = , L
η=
z , H
u(ξ, η) = Vvr (r, z) and v˜ (ξ, η) = δVvz (r, z),
where δ = H/L, the above equations take the following form: ! " h ∂p ∂u 1 2 ∂u =− 2 δ δ u +v + ∂η ρV ∂r Re ∂ξ˜ ! " 1 3 ∂˜v h ∂p ∂˜v + δ2 u + v˜ =− 2 δ ∂η ρV ∂z Re ∂ξ˜ ∂ ˜ ∂ ˜ ξ˜v , 0= ξu + ∂η ∂ξ˜
# ! ∂u 1 ∂u u ∂2 u + + − ∂η2 ∂ξ˜2 ξ˜ ∂ξ˜ ξ˜2 ! # 2 1 ∂˜v ∂v ˜ 2∂ v + + δ ∂η2 ∂ξ˜2 ξ˜ ∂ξ˜
(3.39)
(3.40) (3.41) (3.42)
where Re = VH/ν. Suppose now that we consider the flow in a narrow gap, i.e. δ 1. Thus, we can neglect the terms
of order δ3 and δ2 , but, to have a two-dimensional flow, we have to keep the terms of order δ and unity. After introducing 59
3. Analysis of the flow around a poppet valve in the case of long conical seat
˜ = R/H and v = δ˜v = vz /V, we arrive at the final dimensionless form of the momentum and another rescaling, i.e. ξ = ξ/δ continuity equations, typically encountered in boundary layer theory, e.g. [70]: u
where g =
H ∂p ρV 2 ∂r
∂u ∂u 1 ∂2 u +v = −g + ∂ξ ∂η Re ∂η2 ∂p 0= ∂z ∂ ∂ (ξu) + (ξv) , 0= ∂ξ ∂η
(3.43) (3.44) (3.45)
and Re = VH/ν, where the reference velocity is yet to be defined. Note that (3.44) implies that
the pressure is constant along the z axis, as it is well known from the theory of boundary layers [11] or [70]. If the reference velocity is chosen to be the axial component of the pilot velocity (i.e. V = v pl sin θ), the boundary conditions are vr (r, 0) = V, vr (r, H) = 0, vz (r, 0) = V tan θ and vz (r, H) = 0, which is equivalent to u(ξ, 0) = 1, u(ξ, 1) = 0, v(ξ, 0) = tan θ and v(ξ, 1) = 0. The continuity equation (3.45) is automatically satisfied if the velocity components are defined as u = ξ −1 ∂η Ψ and v = −ξ−1 ∂ξ Ψ. Furthermore, we assume the stream function to be in the form Ψ(ξ, η) = ξF(η). Inserting this into (3.43), we arrive at the boundary value problem 1 000 1 F − FF 00 = −g + ξ Re
(3.46)
with F 0 (0) = 1,
F 0 (1) = 0,
F(0) = tan θ
and
F(1) = 0.
(3.47)
Note that the ordinary differential equation (3.46) is of third order, while we have four boundary conditions. This implies that one parameter in (3.46) has to be reserved as a free parameter. Physically, this means that (a) one cannot set the pressure gradient g and the poppet valve velocity at once, or (b) if the problem is to calculate the flow while both the pressure gradient and the valve velocity are given, one cannot assume a solution in the form Ψ(ξ, η) = ξF(η). Once (3.46) with the appropriate boundary conditions is solved, the wall shear can be evaluated as ! ! V 00 ∂vr ∂vz 1 de f. V + Fξ (0) − 2 Fξ (0) = µ τ˜ ξ . τw,r = µ = µ ∂z ∂r z=0 h h ξ
(3.48)
Note that the radial coordinate ξ is a parameter in (3.46), which means that by solving the boundary value problem, one is able to compute the wall shear at a given ξ location. The actual force due to the wall shear distribution is given by Z a+s Z ξ1 de f. Fdrag,pl = 2ϕrτw,r dr = 2ϕµVh ξτ˜ ξ dξ = 2ϕµVh WS(Re, ξ), (3.49) ξ0
a
where WS(Re, ξ) represents the value of the integral of the wall shear (in dimensionless coordinates), which is to be evaluated by solving the boundary value problem (3.46) on a grid ξ 0 , . . . , ξ1 and then numerically integrate the results. Summing up, the presented technique is capable of calculating the wall shear and the viscous flow with moderate effort. A simple boundary value problem has to be solved, whose solution is the velocity profile at a given location along the generatrix of the cone. Once this solution is found, the wall shear can be simply extracted. The viscous force can be computed by solving the boundary value problem for several coordinates, store the results and then numerically integrate it. Once the boundary value solver is written, this goes very quickly if the solution is initiated from the previous profile. (The author used a grid of 20 points and the complete calculation was performed within 5 seconds on a 1.4 MHz PC in Matlab environment.) However, a significant drawback is that the pressure difference and the body velocity cannot be prescribed at once; if one is chosen to be a user-prescribed parameter, the other one becomes a parameter of the problem which also has to be computed while solving the boundary value problem. The reason for this is the special form of the stream function Ψ(ξ, η) = ξF(η). Avoiding this assumption makes it possible to handle both the pressure difference and 60
3.4. Viscous force
the surface velocity as free parameters - the price is the solution of a partial differential equation rather than an ordinary one. 11
15
0.179
10 5
0.1
9
0.15
10
8
5
u0 [m/s]
6 5
10
0.1
h [mm]
15
0.1
1
7
4 5
1
3
0.05
2
0.1
1 0 −4
−2
0
2
4
6
8
0.1
10
0.2
1
0.3
0.4
0.1 0.5 0.6
0.7
0.8
0.1 0.9
1
1
x [mm]
v [m/s]
Figure 3.13: Velocity profiles computed with 2D di-
Figure 3.14: Difference between Poiseuille profile and
verging flow method (solid line) and Poiseuille profile
2D diverging channel profile, defined by (3.50), given in
o
(dashed line). Parameters: v pl = 10[m/s], α = 21 ,
percents.
r1,in = 2.81[mm], L = 3.21[mm] and x pl = 0.5[mm]
3.4.4 Comparison of the models Let us give a brief comparison of the three wall shear models. First, note that despite its simplicity, the ’equivalent cylinder’ approach captures the essential behaviour of the viscous force, namely that it consists of two parts; the pressuredependent part and the velocity-dependent part. The pressure part (for small gaps) increases proportionally to the gap width, while the velocity-dependent part is inversely proportional to the gap width. The reason for the second phenomena is that if the body moves with a constant velocity, the smaller the gap width is, the larger the slope of the velocity profile at the moving surface becomes. Thus, as the second component of the force is velocity-dependent, from the mechanical point of view, this is a damping force with highly nonlinear characteristics. Also note that comparing (3.23) and (3.31) or (3.32), we conclude that the two free parameters of the cylinder (r 1 cylinder radius and H gap width) can be determined in such a way, that the viscous force on the cylinder and on a given cone equals. Thus, one can call the adjusted cylinder an ’equivalent cylinder’. Let us compare the second case (Poiseuille velocity profile) and the third case (2D flow in a diverging channel). The error of the velocity profiles may be defined as ev = max
! v2D − vP , v pl
(3.50)
where v2D is the 2D solution, vP is the Poiseuille solution and v pl is the cone body velocity. Fig.3.13 depicts two actual velocity profiles and Fig.3.14 gives the error (difference) defined by (3.50) over a range of valve displacements and velocities in percents. We observe a reasonable agreement, which increases if either the surface velocity or the gap width (body position) is increased. For the practical range of these parameters u 0 < 2[m/s] and x < 0.5[mm], the difference is below 1%. Finally, Fig.3.15 presents a direct comparison between CFD results and equation (3.32) in the case of cone body at rest. Interestingly, the theory fails to predict even the trend of the force in the case of the smallest gap (x = 0.1mm). We 61
3. Analysis of the flow around a poppet valve in the case of long conical seat
4 x=0.1mm x=0.2mm x=0.3mm x=0.4mm x=0.5mm
3.5 3
v
F [N]
2.5 2 1.5 1 0.5 0 0
5
10
15
20
∆ p [bar]
25
30
35
40
Figure 3.15: Viscous force on the cone body at rest. Crosses denote CFX runs and solid lines represent (3.32) with α = 21o , r1,in = 2.8mm and l = 3.21mm.
speculate that this can be caused by inaccurate parameter estimation, as if the gap width is fixed (rather then h = x sin α), the trend fits well. For the other four displacements, the author’s opinion is that the results are within reasonable accuracy especially if we take into account that in the CFD runs, maximum 10 cells could be generated across the gap and the wall shear depends on the derivative of the velocity profile at the wall - which might be highly inaccurate for such small number of grid points.
3.5 Summary The main aim of this chapter was to study the flow rate, the pressure drop and the fluid forces in the case of a poppet valve with long, conical seat. The work was a direct extension of the one performed by Bergada and Watton in [12], with the difference that the cone body was moving in this study. The following results were found and conclusion were drawn. 1. Analytical and explicit relationships were developed for calculating the pressure drop and the flow rate through a moving conical poppet valve with long seat. The results - for the particular geometry - were verified with CFD computations and a reasonable agreement was experienced. 2. It was concluded that for the actual valve geometry the linear pressure drop along the gap gives a reasonably good estimation. In general, the range of parameters were this holds has to be verified with the help of the analytical formulae explicitly given in §3.2. 3. Closed-formed expressions were given for the pressure force on the valve body. The results were compared to CFD calculations and a reasonably good agreement was found. 4. By using three models of different complexity, the viscous damping force due to the wall shear distribution was also computed for a general poppet valve. It was concluded that the viscous force consists of a pressure-dependent part and a velocity-dependent part. As the valve body closes and the gap width H decreases, these two components behave counter-wise; the pressure-dependent part of the viscous force vanishes while the velocity-dependent part increases rapidly. This latter phenomenon represents a strong damping force for small valve openings. The numerical results were compared to CFD calculations in the case of valve body at rest and a limited agreement was 62
3.5. Summary
found. It is not clear at this point whether the poor resolution of the CFD grid, the assumptions behind the analytical formulae or these two effects together cause the significant difference for small valve openings.
63
Chapter 4
An analysis of valve chatter 4.1 Introduction The valve chatter phenomenon is charaterized by violent oscillation of the safety valve body, which causes extensive damage to both the valve body and the seat, contributes to fatigue failure of pressure vessel and piping components and produces considerable noise. In a review article [63] in 1991, Petherick and Birk stated that their results of the literature review suggest that the pressure relief valves (PRV) currently installed are based on designs more than 30 years old. Moreover, the lack of understanding of dynamic characteristics of valves has slowed the development of improved PRV dynamic computer models. Maybe one of the first papers studying valve chatter is that one of Kasai [43], who mainly focused on the effect of inlet and outlet piping and viscosity of the oil. Sheffel in [71] conducted several tests corresponding to the influence of the geometry of the closing element and presented some graphs. Roughly speaking, two methods of dealing with the problem were found in literature. The first applies linear stability analysis to determine whether the valve is dynamically stable, e.g. [59] or [26]. The second considers the actual response of the valve piston determined from the numerical solution of the dynamic equations in the time domain, e.g. [87], [62], [16] or [17]. Recently, Hayashi et al. in [35] presented a detailed nonlinear analysis of a direct-operated relief valve together with an upstream chamber and an upstream pipe connecting the chamber and the relief valve. The flow in the pipe was modelled by the unsteady Bernoulli equation, the capacitive characteristic of the chamber gave another differential equation for pressure and the equation of motion for the valve body resulted in a system of four first order ordinary differential equations. Linear stability analysis was analytically studied and then nonlinear analysis was numerically performed by direct simulation and it was shown how period-doubling bifurcations lead to chaos. Unfortunately, little effort was made to interpret their interesting result to engineering guidelines. Another recent and interesting paper is presented by Eyres et al. in [29], where the nonlinear dynamics of a hydraulic damper has been studied thoroughly with the help of novel continuation techniques for discontinuous systems. These two examples show the increasing usage of techniques of nonlinear dynamics in hydraulic valve analysis.
64
4.1. Introduction
4.1.1 Outline of this chapter Referring to Fig.4.2, the aim of this chapter is to compare the dynamic characteristics of the same conical poppet valve with chamfered-type short seat (which results in turbulent discharge flow) and long, conical seat (which allows laminar flow in the gap). The system studied in this chapter is depicted in Fig.4.1. It consists of an upstream pressure vessel, into which a constant rate of liquid volume flow Qin enters. As the pressure in the volume increases, it reaches the set pressure of the relief valve, which opens and discharges the fluid. The pipe in between is assumed to be short enough to be neglected. The main interest is on the linear stability of the valve, notably the required damping for stable operation, the effect of spring precompression and set pressure or inlet volume flow rate. p0
x0 x
x0
Qout x
x0 p0 Qout
p Qin
1
p0 x Qout
B, V p1
Figure 4.1: Pressurised system with relief
p1
Figure 4.2: Geometry of the poppet valves.
valve. In this section, we are interested in the ’core’ mechanisms leading the valve chatter and our aim is to give clear and easy-to interpret results. This requires largely simplified models including only those effects that lead to instability. The starting point is the modelling equations presented in the previous chapters §2.4 and §3. To make the analysis easier and to reduce the parameter space, our first step is to rescale the equations to obtain a form including as few parameters as possible. This way, the physical meaning of the parameters is lost but simplicity is gained. These dimensionless equations are then analysed and the results are first interpreted in term of the dimensionless parameters, as these are universal and not parameter-specific. Later, the results are also given in terms of physical parameters. Finally, two small numerical studies are presented demonstrating the effect of spring precompression and the case when the valve oscillate so heavily that it hits the valve seat (impact dynamics). The first two parts of this chapter §4.2 and §4.3 are analogous. In §4.2.1 and §4.3.1, the equations of motion are developed and suitable dimensionless variables are defined. In §4.2.2 and §4.3.2, the linear stability analysis is performed and stability boundaries are computed, still by means of the dimensionless parameters. Sections §4.2.3 and §4.3.3 are rather technical; the aim of the lengthy calculations is to determine whether the way of loosing stability is ’soft’ or ’hard’. ’Soft’ stability loss means that as the parameter value (say, flow rate) is increased beyond a critical value, oscillations appear ’slowly’ and the amplitude of the oscillation is proportional to the square root of the distance from the critical parameter value (left-hand side in Fig.4.3). Contrary, ’hard’ stability loss is the case when immediately beyond the critical parameter value a large-amplitude oscillation is encountered. Also, in the latter case hysteresis is present, which means that the critical parameter value depends on the direction of its variation, see the right-hand side of Fig.4.3. If hard stability loss is encountered in the system, the critical parameter value obtained by linear analysis cannot be trusted since hysteresis is present. In the final part of this section direct comparison of the two seat types is given in terms of stability in §4.4.1. Also two short numerical studies are presented on the effect of spring precompression in §4.4.2 and impact dynamics in §4.4.3. The aim of these short studies is to highlight those effects which are not studied in detail here but might have a significant influence. Finally, in §4.5, the experimental results of the vibration diagnostics of the pilot-operated pressure relief valve 65
4. An analysis of valve chatter
amplitude
amplitude
stable equilibrium
stable equilibrium
stable periodic motion
stable periodic motion unstable equilibrium
critical value
unstable periodic motion
parameter
critical value
unstable equilibrium
parameter
Figure 4.3: (Left) ’Soft’ stability loss (supercritical Hopf bifurcation). (Right) ’Hard’ stability loss (subcritical Hopf bifurcation.
described in §2.4 are presented. This part is intended to show that the pilot stage of the relief valve is in unstable operation and its frequency can be found in spectrum.
4.2 Poppet valve with chamfered-type seat 4.2.1 Governing equations Consider the system depicted in Fig.4.1 together with the sharp-seated relief valve in the left in Fig.4.2. The equation of motion for the valve body is m x¨ + K x˙ + s (x + x0 ) − A (p1 − p0 ) = 0,
(4.1)
where K is the viscous damping coefficient (not depicted in the figures). The physical parameters are given in Table 2.4. The pressure dynamics in the chamber is given by B p˙ 1 = V
s 2 Qin − Cd Aout , ) (p − p 1 0 ρ
where
H Aout = 2πH r2,in − cos α , 2
(4.2)
with H = x sin α. (For details, see §2.4.) Dimensionless variables are defined as follows. The relative displacement of the valve and the relative pressure in the chamber is denoted by y 1 = x/xre f and y3 = (p − p0 )/pre f . Time is re-scaled as τ = ωre f t. Using these variables, equations (4.1) and (4.2) take the following form y001 = − y03 =
A pre f K 0 s (y1 + e) + y3 y1 − 2 mωre f mωre f mω2re f xre f
Qre f B √ qin − y1 (1 − σ y1 ) y3 . pre f ωre f V
and
(4.3) (4.4)
The dimensionless flow rate is qin = Qin /Qre f and the dimensionless precompression of the spring is e = x0 /xre f . The actual values of the reference quantities are normally chosen in such a way that the equations become as simple as possible. It is indeed possible to get rid of three parameters with an appropriate choice of ω re f , xre f and pre f . However, the resulting reference values are far from the physical reality (e.g. x re f = 552[m] and pre f = 1.8 × 108 [bar]). Thus, we choose such reference values which result in dimensionless values close to unity and we accept the presence of one additional parameter: 66
4.2. Poppet valve with chamfered-type seat
ωre f =
r
s = m
s
811 × 103 [N/m] = 8746[s−1] 10.6 × 10−3 [kg]
f =
ωre f = 1392[Hz] 2π
xre f = 0.23[mm] pre f =
mω2re f xre f 2 π r2,in
(4.5) (4.6)
= 75.73[bar]
Qre f = Cd 2πr2,in xre f sin α
s
(4.7)
2 pre f = 9.14[l/min] ρ
(4.8)
Qre f B = 0.0152[−] pre f ωre f V xre f sin α cos α = 6.244 × 10−4 [−] σ= 2r2,in β=
(4.9) (4.10)
The mass of the cone body mcone = 4.1[g] is comparable to the mass of the spring, which is m spring = 19.6[g]. In such cases, the mass of the spring should also be taken into account when calculating the eigenfrequency in (4.5) as given by (2.42) in §2.4.3. The result is m = mcone + m spring /3 = 10.6[g]. The spring stiffness was measured and was found to be s = 811[kN/m]. The reduced bulk modulus of the fluid is B = 1.55 × 10 9 [Pa], the volume of the pipe is
V = 2.356 × 10−4 [m3 ]. The geometrical parameters of the actual valve are α = 21o and r2,in = 2.8[mm]. Note that in the physical system, the highest flow rate occurring in the pilot stage is approx. 1[l/min], thus the interesting range of q in is qin ≈ 0 . . . 0.1. Due to the definition of xre f in (4.6), the precompression of the spring is around e ≈ 0 . . . 1, with e = 1 in the case of direct comparison with measurements. The damping on the valve body is highly nontrivial to estimate. Based on (2.44) in §2.4.3 the estimated damping coefficient is K = 0.0215[N s/m]. This gives k=
K 0.0215[N s/m] = = 2.32 × 10−4 , mωre f 10.6[g] × 8746[1/s]
(4.11)
which is an extremely low damping value; note that referring to (2.45) 185.4[N s/m] would be needed to obtain overdamped (i.e. aperiodic) motion.
4.2.2 Linear analysis with no spring precompression The pilot valve body displacement is typically below 1[mm], hence y 1 < 4, thus in (4.4) we have σy1 < 2.5 × 10−3 , which can be neglected if compared to unity. Consider the first order form of system (4.3) and (4.4) with σ = 0, given by
y01 =y2 y02 = − ky2 − y1 − e + y3 . √ y03 =β qin − y1 y3
(4.12)
The equilibrium is given by y2,e = 0, y3,e = y1,e + e and y21,e y1,e + e = q2in . If e , 0, an analytical solution seems to
be hard to achieve. In order to allow analytical calculations, we consider the case of e = 0. The e , 0 case is studied numerically later on in §4.4.2.
2/3 Consider (4.12) with no spring pre-compression, i.e. e = 0. The equilibrium is given by y e = (qin , 0, q2/3 in ). Our
first aim is to study the local stability of the equilibrium. Expanding (4.12) into Taylor series around y e , dropping the nonlinear terms and computing the characteristic polynomial of the linear coefficient matrix gives D(λ) = λ3 + λ2 (k + z) + λ (kz + 1) + 3z = 0, 67
with z =
β 1/3 q . 2 in
(4.13)
4. An analysis of valve chatter
By inserting λ = iω into (4.13) it can be shown that self-excited oscillation appears in system (4.12) below a critical value of the damping coefficient k, given by1 kcrit
# " q 2 1 2 2 2 1 + z + 8z . = − 1+z + 2z
(4.14)
The frequency of the appearing oscillation is
ω2 = zkcrit + 1 =
" # q 1 1 + z2 2 + 8z2 . 1 − z2 + 2
(4.15)
ω
kcrit 0.7
1.6
stable
0.6
1.5
PSfrag replacements
0.5 1.4
0.4 0.3
unstable
1.3
PSfrag replacements
0.2
1.2 1.1
0.1 1
2
3
4
5
z
1
Figure 4.4: Critical damping, defined by (4.14).
2
3
4
5
z
Figure 4.5: Frequency of the appearing cycle, defined by (4.15).
Fig.4.4 and Fig.4.5 depict the critical damping coefficient and the frequency of the appearing cycle, respectively. Note that the physically relevant range of parameters is 0 ≤ qin ≤ 0.1, i.e. 0 ≤ z ≤ 0.0035. However, a larger interval is √ plotted for the sake of completeness. The maximum of the damping is located at z = 1, k max = −1 + 3 = 0.732. The √ frequency of the cycle is unity (i.e. the eigenfrequency of the spring-mass system) at z = 0 and tends to 3 as z → ∞. We also check whether the velocity of the eigenvalue is non-zero when crossing the imaginary axis, i.e. at the critical point we have
d D(λ,z) (ω4 − 3)(ω2 − 3) d( 1.323 we have a ’hard stability loss’ (or subcritical bifurcation) and an unstable limit cycle co-exists with a stable equilibrium, see the right panel of Fig.4.3. At the stability change z = 1.323, the vanishing Lyapunov coefficient gives rise to a curve corresponding to fold of periodic orbits, which is a global bifurcation and we shall not study here. The actual formula of a is not given here as it is rather lengthy and does not provide a deeper understanding of the phenomena.
70
4.3. Poppet valve with long conical seat
4.3 Poppet valve with long conical seat In this section an analysis is presented, which is analogous to that one given in the previous section. The only difference is that conrary to the sharp chamfered-type valve seat, we consider now the case of long, conical seat. The train of thought is exactly the same as in the previous chapter. However, the recurring equations are given again for the sake of completeness.
4.3.1 Governing equations Consider the system depicted in Fig.4.1 together with the long, cone-seated relief valve in the right of Fig.4.2. The equation of motion for the valve body is m x¨ + K x˙ + s (x + x0 ) − A (p1 − p0 ) = 0,
(4.30)
and the pressure dynamics in the chamber is given by p˙ 1 =
B Qin − R (p1 − p0 ) x3 V
(4.31)
Dimensionless variables are defined as follows. The relative displacement of the valve and the relative pressure in the chamber is denoted by y1 = x/xre f and y3 = (p − p0 )/pre f . Time is re-scaled as τ = ωre f t. Using these variables, equations (4.30) and (4.31) take the following form y001 = − y03 =
A pre f K 0 s (y1 + e) + y − y3 mωre f 1 mω2re f mω2re f xre f
Rx3re f B qin − y3 y31 ωre f V
(4.32)
(4.33)
The dimensionless flow rate is qin = Qin /Qre f and the dimensionless precompression of the spring is e = x0 /xre f . The actual values of the reference quantities are chosen in a similar way as in the previous section: s r ω s 811 × 103 [N/m] −1 ωre f = = = 8746[s ] f = = 1392[Hz] m 2π 10.6 × 10−3 [kg] π sin4 α 1 = 0.7144[ms/kg] µ 6 ln 1 + L sin α r2,in r 3 ωV = 0.232[mm] = BR mω2 xre f = = 75.73[bar] A
R=
xre f pre f
Qre f = R x3re f pre f = 3.95[l/min] β=
(4.34) (4.35)
(4.36) (4.37) (4.38)
Qre f B = 0.00657[−] pre f ωre f V
(4.39)
The reduced bulk modulus of the fluid is B = 1.55 × 109 [Pa], the volume of the pipe is V = 2.356 × 10−4 [m3 ]. The
geometrical parameters of the actual valve are α = 21o , l = 3.21[mm] and r2,in = 2.8[mm]. Note that in the physical system,
the highest flow rate occurring in the pilot stage is approx. 1[l/min], thus the interesting range of q in is qin ≈ 0 . . . 0.25. Similarly, the precompression of the spring is around x0 ≈ 0 . . . 0.3[mm], thus we have e ≈ 0 . . . 1.3.
71
4. An analysis of valve chatter
4.3.2 Linear analysis with no spring precompression Consider the system
The equilibrium is given by y2,e = 0, y3,e
y01 =y2
(4.40) y02 = − ky2 − y1 − e + y3 . y03 =β qin − y31 y3 = y1,e + e and y31,e y1,e + e = qin . If e , 0, an analytical solution seems to be
hard to achieve. However, by using implicit differentiation it is easy to show that for positive e and q in values only one solution for y1,e exists for all qin and e > 0. To allow analytical calculations, we consider the case of e = 0. The e , 0 case is studied numerically in §4.4.2. √ √ Consider (4.40) with no spring pre-compression, i.e. e = 0. The equilibrium is given by y e = ( 4 qin , 0, 4 qin ). Expanding (4.40) into Taylor series around ye , dropping the nonlinear terms and computing the characteristic polynomial of the linear coefficient matrix gives D(λ) = λ3 + λ2 (k + z) + λ(kz + 1) + 4z = 0,
3/4 with z = β qin .
(4.41)
By inserting λ = iω into (4.41) it can be shown that self-excited oscillation appears in system (4.40) below a critical value of the damping coefficient k, given by kcrit
" # q 1 2 2 2 4 = − 1+z + 1 + z + 12z . 2z
(4.42)
The frequency of the appearing oscillation is
" # q 1 2 2 2 4 1−z + ω = zkcrit + 1 = 1 + z + 12z . 2 2
(4.43)
Note that if k = kcrit , we have λ1,2 = ±iω and λ3 = −(kcrit + z). We conclude that unless z , 0, the third eigendirection is stable.
4.3.3 Two-dimensional centre manifold First, let us check whether the velocity of the eigenvalue is non-zero when crossing the imaginary axis, i.e. at the critical point we have d D(λ,z) (ω2 − 4)(ω4 − 4) d( 0 , H(x, x˙) < 0
(5.3)
√ |z|, p s denotes the supply pressure, pr is the relief pressure, c is the valve orifice coefficient
(including the density, valve area gradient and the discharge coefficient) and u is the control signal. Depending on the type of the control valve, the control signal u can have several forms. In this study, we consider two particular forms of u, given by u = −sign H(x, x˙) u = −H(x, x˙)
for two-stage flow control valve (relay-like behaviour) and
(5.4)
for proportional valve or servo-valve (proportional behaviour),
(5.5)
where H(x, x˙) is the indicator function, which is chosen to be a simple PD controller: H(x, x˙) = S (x − xd ) + D x˙,
(5.6)
with xd desired position and S stands for the proportional gain (as p stands for pressure). We can simplify the above equations significantly by introducing dimensionless variables and parameters. Let x be replaced by x˜ = x/L, where L represents the ’initial length’ of the chambers (i.e. V0 = AL) and we also define p˜ dimensionless pressure as p˜ = √ (p − p0)/(p s − p0 ). The dimensionless flow rate is represented by q = Q/Qmax , where Qmax = c p s − pr . A compressible fluid with bulk modulus B in a chamber with length L and cross-sectional area A forms a spring with stiffness AB/L. √ The eigenfrequency of this fluid spring together with the mass m yields ω = AB/Lm. After rescaling the time with this frequency, inserting the above defined dimensionless variables, dropping the tilde and introducing the general system variables y = ( x˜, x˙˜, p˜ 1 , p˜ 2 )T , we arrive at the system y01 = y2 y03 = with
y02 = α(y3 − y4 ) − ky2
1 (γ q1 − y2 ) α(1 + y1 )
if u G(y3 ) q1 = u G(1 − y3 ) if
y04 =
h(y) > 0 h(y) < 0
1 (γ q2 + y2 ) , α(1 − y1 )
−u G(1 − y4 ) if and q2 = −u G(y4 ) if
h(y) > 0
(5.7)
,
(5.8)
h(y) < 0
where ...0 represents the differentiation with respect to the dimensionless time τ = ω t and the dimensionless switching function is given by
h(y) = sy1 + dy2
(5.9)
(assuming xdesired = 0). The dimensionless physical parameters are defined by α=
p s − pr , B
γ=
Qmax , V0 ω
84
k=
K , mLω2
(5.10)
5.3. Relay control
par.
value
unit
par.
value
unit
par.
value
unit
m
200
kg
ρ
870
kg/m2
α
m2
B
1.125 × 109
−
A
1.32 × 10−3
6.13 × 10−3
Pa
ω
136
rad/s
0.4
m
ν
3
ps
70
bar
pr
1
bar
L V0 Qmax
5.28 × 10 30
−4
m
l/min
70 × 10
−6
2
m /s
2.78 × 10
γ
−3
−
Table 5.1: Physical and dimensionless parameters.
and the dimensionless control parameters are s = S L and d = DLω. Typical physical parameters are summarised in Table5.1. Note that although these parameters correspond to a given system, the aim of the table is to give an estimation for the order of magnitude of the dimensionless parameters.
5.3 Relay control Consider system (5.7) and (5.8) together with the relay feedback PD controller u = −sign h(y1 , y2 ) = −sign (sy1 + dy2 ) .
(5.11)
The physical explanation of the relay-like behaviour is that a two-way directional control valve is employed to control the piston and it is assumed that the switching time is significantly less then the characteristic time scale of the system. Mathematically speaking, we are interested in the dynamic behaviour of the discontinuous system of nonlinear differential equations, which can be written in a general form as F1 (y) if h(y) > 0 0 y = . F2 (y) if h(y) < 0 PSfrag replacements
(5.12)
y2
y(0)
F1
h(y )=
y(T 1 + T 2 )
0
y1 F2 y(T 1 )
Figure 5.2: Switching dynamics and boundary conditions for continuation. Referring to Fig.5.2, the so-called switching hyperplane (see [54]) defined by S := {y : h(y) = 0} divides the phase space into two subregions with different continuous dynamics F 1 (y) and F2 (y), respectively. In some cases there exists ˆ which is simultaneously attracting from both sides and the solution trajectory is constrained a part of S denoted by S, ˆ until it reaches its boundary or the equilibrium. A popular control technique is to tune the control to evolve within S,
parameters to generate such sliding motion. However, it can be shown that sliding mode is not possible in this case; consider the projection of the two vector fields onto the switching plane: h∇h, F 1 i and h∇h, F2 i. S is attracting from the 85
5. Stability of controlled hydraulic systems
’positive side’ (h > 0) if this projection is negative, while the switching plane is attracting from the ’negative side’ if the corresponding projection is positive. However, it turns out that the value of this projection is the same for both sides h∇h, F1 i = h∇h, F2 i = sy2 + dα(y3 − y4 ) − d k y2 .
(5.13)
As the same quantity cannot be positive and negative at the same time, sliding mode control is not possible (for any friction model) unless other state variables are included in the control law, e.g. pressure or acceleration as in [61]. Another important issue is that the pressure-to-bulk modulus ratio α (defined by (5.10)) is a small parameter; (α < 0.01), which means that the solutions evolve on two different time scales. The ’fast’ subsystem consists of the two equations for y03 and y04 (pressure in the chambers) as the magnitude of their derivatives is of order α −1 and the ’slow’ subsystem is given by the equations for position and velocity (y 1 and y2 ). In the next section we have a closer look at the slow-fast nature of the system.
5.3.1 Application of singular perturbation theory for relay-like control If the fluid becomes more and more stiff, parameter α (defined by (5.10)) decreases and finally vanishes. If α = 0, the ordinary differential equations for the chamber pressures (y03 = . . . and y04 = . . . in (5.7)) become algebraic equations and the dimension of the ODE2 system reduces from four to two. In this brief section, with the help of geometric singular perturbation theory (for details, see [30]) we give an explanation why the fluid compressibility has a fundamental influence on the dynamics of the system. We consider system (5.7) and (5.8) in the standard form of singular perturbation theory (for being coherent with the literature [30] and [69], we forget our previous notations). The ’slow’ form of the equations is u˙ = F (u, v, α) . α v˙ = G(u, v, α)
(5.14)
As F and G are of O(1), it can be seen that v ∈ 2 are the ’fast’ variables corresponding to the pressure variables and u ∈ 2 are the slow ones; position and velocity. If we switch to the fast time scale t = τ/α, we obtain the equivalent fast system
0 u v0
= α F (u, v, α) = G(u, v, α)
.
(5.15)
By taking the α → 0 limit, the slow system (5.14) turns into a lower-dimensional ODE of the slow variables (the reduced system) evolving on a surface called (critical manifold) given by G(u, v˜ , 0) = 0 or explicitely v˜ = m 0 (u). Also, this critical manifold consist entirely of the equilibria of the limiting layer problem, i.e. (5.15) in the case of α = 0. Also note that the this layer problem is the flow of the fast variables parametrised by the slow variables. The situation is depicted in Fig.5.3. System (5.14) and (5.15) can be decomposed into two complementary systems with lower dimensions; the flow of the slow variables on the critical invariant manifold and the complementary flow of the fast variables. If the flow of the fast variables is stable, the trajectories tend to the critical manifold and then they are forced to evolve on it. Thus, we first have a look at the stability of the complementary flow of the fast variables, given by ∂G δv0 = δv. ∂v ˜v=m0 (u)
In our case the critical manifold consists of two parts, on both sides of the switching manifold, given by 2 y 2 y2 2 if h(y) > 0 if h(y) > 0 1 − γu γu y 2 y 2 and y4 = . y3 = 2 1− 2 if h(y) < 0 if h(y) < 0 γu
2 Ordinary
γu
Differential Equation
86
(5.16)
(5.17)
5.3. Relay control
fast dynamics
v’=G(u,v)
old anif m g chin swit
critical manifold v=M(u)
y3 - y4
y2
slow dynamics
y1
u=F(u,M(u))
Figure 5.3: Slow and fast manifold. Starting at an arbitrary initial state (close enough to the critical manifold), in the initial state the fast dynamics is dominant and the trajectory soon lands on the critical manifold. As the critical manifold is invariant, the trajectory cannot leave it and thus after the initial transient state, the restricted dynamics on the critical manifold is sufficient to study.
Along the critical manifold, we have !2 y 2 , y3 + y4 = 1 and y3 − y4 = −sign(h) 1 − 2 γ
(5.18)
and the eigenvalues of the layer problem turn out to be f λ1,2 =−
γ2 ≤ 0, 2α (1 ± y1 ) |y2 |
(5.19)
thus we conclude that the complementary (fast flow) always tends to the critical manifold (the superscript ’f’ denotes ’fast’). Thus, it is sufficient to study the restricted dynamics, which is the differential equation for y 1 (position) and y2 (velocity) with (5.18). Without giving here the actual calculations, it can be shown that on both sides of the switching manifold, the linear part of the system close to the origin (y1 , y2 ) = (0, 0) is equivalent and the eigenvalues are λ1s = 0 and λ2s = −k (the superscript ’s’ denotes ’slow’). Thus we conclude that the incompressible (restricted) system is linearly neutral (it is on the border of stability) and the stability of the equilibrium is governed by the pressure dynamics - thus the compressibility. When analysing the stability of such a hydraulic system, the fluid compressibility cannot be neglected.
5.3.2 Asymptotically smooth analysis Having shown that the pressure dynamics has to be included in the analysis, our aim is to predict the onset of oscillations. If the critical parameter values are known, the controller parameters s and d can be chosen easily. We consider again system (5.7) and (5.8) together with the relay feedback PD controller (5.11). Unfortunately, standard linear stability analysis cannot be performed as the system is discontinuous (not only F 1 and F2 are not differentiable at the origin, but they also do not equal). One way to overcome this problem is to ’smoothen’ the discontinuity with an appropriate function equipped with a parameter ε such that when taking the ε → 0 limit, we regain the original discrete system (see e.g. [53]). Linear stability analysis of the smooth system is then performed to obtain the stability boundary as a function of the smoothing parameter. Then, by taking the ε → 0 limit one obtains the critical parameter values for the discontinuous 87
5. Stability of controlled hydraulic systems
system. The smoothing functions σL(e f t) and σR(ight) are σR,L (y, ε) = 1/2 1 ± tanh (h(y)/ε) and this way, the system can be re-written as a fourth-order continuous nonlinear system y01
= y2
y02
= α(y3 − y4 ) − ky2 i p 1 1 h √ = γ −σR y3 + σL 1 − y3 − y2 α 1 + y1 i 1 1 h p √ = γ σR 1 − y 4 − σ L y 4 + y 2 . α 1 − y1
y03 y04
(5.20)
The equilibrium is given by ye = (0, 0, 1/2, 1/2)T (thus, we can drop the signum and absolute function in G). Note that by smoothing the discontinuity, we have introduced an equilibrium artifically; there exists no equilibrium for the original system (5.7) together with relay control. The linear coefficient matrix (Jacobian) of (5.20) evaluated at the equilibrium y e yields
A y=ye
0 0 = sγ − √2αε sγ √
2αε
whose characteristic polynomial is given by D(λ) =
1
−k
0 α
1 + √dγ2ε dγ 1 √ 1 + α 2ε
− α1 P4
i=0
− √γ2α 0
0 −α , 0 − √γ2α
(5.21)
ai λi with a4 = 1. The equilibrium is locally attracting (stable)
if every eigenvalue of the Jacobian is with negative real part. However, as a parameter (in this case, d) is varied, the eigenvalues of the Jacobian also move on the complex plane and at a critical parameter value d H , a real eigenvalue may cross the imaginary axis, which leads to repelling (unstable) behaviour. If a pair of complex conjugate eigenvalues cross the imaginary axis, oscillation (periodic orbit) is born in the system, which is called Hopf bifurcation or dynamic instability. Once the characteristic polynomial is known, the critical parameter value can be computed by assuming purely imaginary solutions and substituting iω into λ, which leads to the condition a 21 − a1 a2 a3 + a0 a23 = 0. Performing the actual computations, the critical parameter value and the frequency of the appearing oscillation is given by √
√ 2 2sα + kγ dH = √ −ε 2αγ 2kα + γ 2sα
2
and ω =
√ 2γ 2sα , √ ε 2kα + γ
√
(5.22)
where subscript H refers to Hopf3 . The results allow the qualitative description of the influence of the system parameters on the stability of the control. As we are interested in the properties of the original discontinuous system (ε → 0), we focus on the first term in dH . If d > dH , the system is stable and the control law is able to drive the mass into the desired position. However, there exists a critical value of the differential gain d H , below4 which the system is unstable and oscillations occur. Increasing of γ (raising Qmax ) stabilises the system, while increasing α - large supply pressure or highly compressible fluid - destabilises the system. Large controller gain s also worsens the performance, but it can be shown that if s < d k, the control is unconditionally stable for all positive α, k, d and γ values.
5.3.3 Piecewise continuation Once the linear analysis is done and the behaviour is understood, the nonlinear effects should be studied. (It is often forgotten that nonlinear terms may completely change the results obtained by linear analysis.) Apart from a few extremely simple or special cases, nonlinear analysis requires some numerical solution. This section shows how to find the periodic 3 The
parameter value d exists above and beyond d H but the oscillation and thus the frequency appears only below or above d H . Thus, defining a
’critical frequency’ ωH would be meaningless and confusing. 4 The bifurcation is subcritical, see the numerical results in the next section.
88
5.3. Relay control
orbits of (5.7) and (5.8) numerically in a convenient and efficient way. If an initial value solver (e.g. some Runge-Kutta scheme) would be applied, a large number of steps would be required to converge to the periodic orbit as the solution is usually started from an initial condition chosen by chance. Contrary to this, if the problem is handled as a boundary value problem (BVP), only the periodic orbit will be computed. In addition, when integrating with an initial value solver, one has to include an event handling algorithm to locate the points where the trajectory hits the switching surface (smoothing the discontinuity results in stiff equations, which are expensive to solve). This is completely avoided in the case of the BVP formulation. The continuation software AUTO [27] is a general-purpose program to study ordinary differential equations qualitatively. In this case, AUTO was initiated with the data of a periodic orbit obtained by an initial value solver. Then, one can vary the parameters systematically and compute the periodic orbits, the period of the solutions and locate special points (e.g. turning points, period doubling, etc...). Moreover, the boundary value formulation is insensitive to the stability of the orbits, so unstable periodic orbits can also be easily found. In general, a limit cycle satisfies the boundary conditions h(y(0)) = 0, h(y(T 1)) = 0, h(y(T 2)) = 0 and y(0) = y(T 1 + T 2 ) where T 1 and T 2 represent the time intervals spent at the ‘right’ and ’left’ half-plane, respectively (see Fig.5.2). The first three conditions state that the trajectory hits the switching surface at T = 0, T 1 , T 1 + T 2 while the fourth condition ensures that the orbit is closed. It was numerically verified that there is a limit cycle present at α = 0.01, k = 0.5, γ = 0.5, s = 10 and d = 0.1 and that the orbit is symmetric (T = T 1 = T 2 ). To remove the unknown period T from the boundary conditions, another rescaling of the time is performed: ζ = τ/T . This way, the period of the limit cycle appears as an additional parameter in the right-hand side of the equations and it is a parameter computed for every limit cycle. The results of the continuation process are given in Fig.5.4. In the numerical computations, the differential gain d −4
0.8
4
x 10
max(y1)
T
0.6
0.4
2
0.2
0
0
0.1
0.2
0.3
0
0.4
0
0.1
d
0.2
0.3
0.4
0.3
0.4
d
−3
PSfrag replacements
1.1
2.5
1
2
0.9
1.5
0.8
1
0.7
0.5
0.6
0
max(y1 )
x 10
max(y3)
max(y2)
3
0
0.1
0.2
0.3
0.5
0.4
d
0
0.1
0.2
d
Figure 5.4: Piecewise continuation, computed with AUTO (thick solid line) and the analytically calculated critical value at d ≈ 0.278 (dashed vertical line). Parameter values: α = 0.01, k = 0.5, γ = 0.5 and s = 10.
was chosen as a free parameter and it was found that the oscillation shrinks as d is increased and vanishes at d H ≈ 0.278,
just as predicted by (5.22), when ε → 0. Contrary to the analogue Hopf bifurcation in smooth systems, the period of the 89
5. Stability of controlled hydraulic systems
appearing limit cycle is zero, which is again predicted by (5.22), when ε → 0. Note that after the limit cycle appears, the amplitude of the pressure oscillations (max(y3 )) grows rapidly while the oscillations in the position remain small. This experience is consistent with the result of the previous chapter, namely that the instability appears in the pressure dynamics (’fast’ variables y3 and y4 ). It was numerically verified that the smoothing technique described in the previous chapter gives correct results and (5.22) gives correct values for the critical value of the differential gain d H .
5.4 Proportional control Let us now turn to the case of proportional control, where the control signal u is directly but nonlinearly proportional to p the flow rate: Q ∝ u ∆p. In the experimental rig, such a proportional control valve is mounted (type Bosch Rexroth 4WRAE). The train of thought is essentially same as in the previous section; smooth the discontinuity, perform linear
stability analysis and verify the results with numerical computations.
5.4.1 Theoretical analysis of analogue control Consider system (5.7) and (5.8) together with the proportional PD controller u = −h(y1 , y2 ) = −(sy1 + dy2 ),
(5.23)
which corresponds to the behaviour of a proportional hydraulic valve. The system of differential equations has a set of de f.
equilibrium points; y1,e = 0, y2,e = 0 and y3,e = y4,e = pe . Physically this means that the control valve is closed, no flow is entering or leaving the piston and the chamber pressures are arbitrary as long as they equal. In this case, the vector field is not continuously differentiable, but the two functions F 1 and F2 match on the sliding surface (along h(y) = 0). As a first step, we employ the same smoothing technique as described in the previous chapter. In this case, the linearisation is performed around the parametric equilibrium ye = [0 , 0 , pe , pe ]T , which leads to the stability condition p √ pe + 1 − p e s 2 d ≥ − . p k γ 1 + 2 √ pe 1 − p e
(5.24)
The above expression reveals that if d is large enough, the set of equilibria is stable and the control is able to drive the system into the desired position. When decreasing d, the first point losing the stability is given by y e,0 = [0 , 0 , 12 , 12 ] and √ this happens at the critical parameter value dcrit = s/k − 2/γ. The frequency of the bifurcating limit cycle is given by √ ω2 = 2γs/k. Keeping in mind that the smoothing method may miss special properties of discontinuous systems, we follow the work of Leine et al [54]. Assume that y˜ e is a fixed point of the dynamical system (5.12), it is located at the switching hypersurface S, e.g. h(˜ye ) = 0 and that the vector field is continuous along S, i.e. F 1 (˜ye ) = F2 (˜ye ). We can compute two Jacobian matrices on both sides: J1 (˜ye ) and J2 (˜ye ) (we do not emphasise here the dependence of J on the parameters). ˜ which is set-valued at y˜ e and is the convex hull of J1 (˜ye ) and J2 (˜ye ): The idea is to define a ’generalised Jacobian’ J, J˜ = {(1 − η)J1 (˜ye ) + ηJ2 (˜ye ), ∀η | 0 ≤ η ≤ 1} .
(5.25)
(5.25) gives the set of values, which the generalised Jacobian can attain on the switching surface. Having defined the set-valued Jacobian, we can compute its set-valued eigenvalues and by varying η from 0 to 1, a unique ’path’ of each eigenvalue is given ’during’ the jump at the discontinuity. ’Discontinuous bifurcation’ is then defined as the case when 90
5.4. Proportional control
−3
2.21
4
2.2
max(y )
3
T
1
2.19 2.18
2
1
2.17 2.16
x 10
1
1.2
1.4
1.6
1.8
0
2
1
1.2
1.4
s 0.01
2
1.6
1.8
2
0.8
max(y )
0.006
3
2
max(y )
1.8
0.9
0.008
0.004
0.7
0.6
0.002 0
1.6
s
1
1.2
1.4
1.6
1.8
0.5
2
1
1.2
s
1.4
s
Figure 5.5: Piecewise continuation, computed with AUTO (solid line represents stable, dash-dot line represents unstable periodic orbits) and the analytically calculated critical value at s ≈ 1.464 (dashed vertical line). Parameter values: α = 0.01, k = 0.5, γ = 0.5 and d = 0.1.
the set-valued eigenvalues contain a value on the imaginary axis for some η. After performing the actual computations, we arrive at the two critical parameter values: √ s 4 2 d1 = − k 3 γ
√ 2 s and d2 = − . k γ
(5.26)
For d1 ≤ d ≤ d2 a pair of complex eigenvalues cross the imaginary axes as η is varied, while for d < d 1 , the eigenvalues are located on the right-half complex plane for 0 ≤ η ≤ 1. Thus, we conclude that for a stable fixed point, we must have √ d > s/k − 2/γ, which is the same result as obtained by the smoothing method. An interesting result is that again, by setting dk > s, one can ensure unconditional stability similar to the relay case.
5.4.2 Numerical results for analogue control To verify the theoretical results numerically, the periodic solution of system (5.7), (5.8) with (5.23) was computed by AUTO, as described in the previous section. The results are given in Fig.5.5. Note that the continuation was performed by varying s (instead of d as in the relay case) due to numerical reasons. By setting α = 0.01, k = 0.5, γ = 0.5 and d = 0.1, the critical parameter value, the frequency and the half-period of the appearing limit cycle is s √ √ T π 2 2γscrit scrit = k d + ≈ 1.439 and = ≈ 2.184. ≈ 1.464, ω = γ k 2 ω
We observe an excellent agreement with the numerical results given in Fig.5.5. The continuation also reveals that although
one would expect so, the phenomenon is not analogous to the Hopf bifurcation in continuous systems as an unstable limit cycle was found below the critical parameter value. Although this unstable limit cycle cannot destroy the stability of the fixed point(s) (for d > d1 ), its presence is somewhat unexpected in the light of ’smooth bifurcation theory’. 91
5. Stability of controlled hydraulic systems
y
PSfrag replacements t u
ti
ti+2
ti+1
t ts
ts
ts
Figure 5.6: Process of digital data acquisition and control.
5.4.3 Effect of digital control The stability analysis presented in the previous section was interesting from a mathematical point of view but little was known whether the same boundaries or at least tendencies can be found in a real experimental rig. Unfortunately, the author was unable to reproduce any of the theory with the help of the experimental rig. If the differential gain was switched off or was kept at a relatively small value, the positioning was stable. As the proportional gain was increased, the precision improved but after a critical value the system began to oscillate. In such cases, theory (5.24) suggests to add some differential gain to stabilise the system. However, when increasing d, the situation got even worse; the oscillation became wilder. Moreover, in some cases a proportional gain value, with which the system was stable with no differential gain became unstable when switching on the differential gain. Another interesting issue was that the frequency of the oscillation was always about twice the sampling frequency. These phenomena are typical for stability problems in digital control, see [78]. So far it was assumed that the control is ’analogue’, i.e. the PC performing the data acquisition, computing the control signal values and sending the signals to the proportional valve is ’fast enough’ not to effect the stability. However, reality is different and the real process is depicted in Fig.5.6. With a given sampling period t s a ’screenshot’ of the physical reality is taken (say) at time ti (e.g. the actual position is sampled) and based on this sample, the control law is evaluated. The control signal itself is constant (zeroth-order holder) and is valid from t i to ti+1 . Thus, there is a delay of about t s as the control signal valid at ti+1 was computed on the bases of a signal at ti . For the ease of calculations, this model neglects the time needed for control law evaluation, which increases the dimension of the problem and makes it hard to develop closed-form solutions. Our starting point is (5.7) and (5.8) with proportional control (5.23). By linearising around the equilibrium y e,0 = 1 2
[0 , 0 , , 21 ] and introducing the new variables z1 = y1 , z2 = y2 , z3 = y3 − y4 and z4 = y3 + y4 (after a lengthy but straightforward calculation), we arrive at 0 z1 0 z = 0 2 √ 2γ − α s z3
1 −
√
−k
2γ α
d−
2 α
0 α 0
z1 z 2 z3
and z04 = 0.
(5.27)
Assume that the dimensionless sampling period is τ s . The control is a zeroth order holder, which means that in the time 92
5.4. Proportional control
interval τ j ≤ τ ≤ τ j+1 (τ j+1 = τ j + τ s ), the control law is constant and its value is given by h i hτ = − s z1 τ j + d z 2 τ j
for τ j ≤ τ ≤ τ j+1 = τ j + τ s .
(5.28)
(5.27) is a linear system of ordinary differential equations. After introducing the new time scale T = τ/τ s and some lengthy calculations, the solution of (5.27) - taking into account the initial conditions z 1 (τ j ) = z10 , z2 (τ j ) = z20 and z3 (τ j ) = z30 - is z1 (T ) = τ s BT +
A1 ω1 T A2 ω2 T + +C e e ω1 ω2
z2 (T ) = B + A1 eω1 T + A2 eω2 T z3 (T ) =
!
(5.29) (5.30) !
A1 ω1 T A2 ω2 T 2τ s BT + e + e +C − α ω1 ω2
√
!
k 2γ , hτ τ s T + α 2
(5.31)
where ω1,2 =
√ τs −k ± k2 − 8 , 2
hτ = − (s z10 + d z20 ) √ ω1 A1 = 4τ s 4τ s z20 + 2αω2 z30 + 2hτ γ(2τ + kω2 ) ω1 − ω 2 √ ω2 A2 = 4τ s 4τ s z20 + 2αω1 z30 + 2hτ γ(2τ + kω1 ) ω2 − ω 1 γ B = √ hτ 2 √ z10 2γ α C= + z20 − khτ . τ2 2τ s 4τ s The solution can be given in a compact form: zn+1 = T zn + b,
(5.32)
which is an iteration. Successful control means zn → zd(esired) as n → ∞. The stability of such an iteration depends on the eigenvalues of the matrix T , which are called often characteristic multipliers and are usually denoted by µ. The characteristic multipliers are the roots of the polynomial equation det (T − Iµ) = p3 µ3 + p2 µ2 + p1 µ + p0 = 0,
(5.33)
where I stands for the identity matrix (with appropriate dimensions). If all the characteristic multipliers lie within the complex unit circle, the iteration is stable. There are three basic ways of loosing stability [50]: 1. µ ∈ and µ = 1; if one the characteristic multipliers is exactly unity (real number), we have a fold bifurcation and two new equilibria are born (or vanish). 2. µ ∈ and µ = −1; a cycle of period two is born (every second iteration of the map gives the same value), called flip bifurcation. 3. µ ∈ and |µ| = 1, i.e. two complex conjugate characteristic multipliers lie on the complex unit circle. This case is called Neimark-Sacker bifurcation and corresponds to the birth of a limit cycle. The above conditions (i.e. the number of eigenvalues lying inside the complex unit circle) can be verified with the help of the argument principle, see [45] for theoretical background or [79] for a similar application. In what follows, we verify whether these bifurcations can occur in our system for special parameter values. However, although the calculations are 93
5. Stability of controlled hydraulic systems
performed without any restriction on the parameters, to avoid large formulae, only those ones corresponding to the special case of k = 0 are given here explicitly. The coefficients pi in (5.33) for k = 0 are p3 = 1,
(5.34) n o √ √ √ √ 1 −2 + 2γ (d + sτ s ) − 4 + 2dγ cos( 2τ s ) − sγ sin( 2τ s ) , p2 = (5.35) 2 √ √ √ √ √ 1 p1 = −2 + e−ı 2τs sγ ı − 2τ − eı2 2τs ı + 2τ s + 4 cos( 2τ s ) and (5.36) 2 i o √ h √ √ 1n −2 + 2γ sτ s + d cos( 2τ s ) − 1 − sγ sin( 2τ s ) , (5.37) p0 = 2 √ where ı stands for the complex unity, ı = −1. To locate the multipliers on the unit circle, we do not use the argument principle but due to the low order of the problem, apply the Moebius transformation µ=
λ+1 , λ−1
(5.38)
which maps the complex unit circle onto the left-half plane. Table 5.2 gives some of the coordinates of the mapping. µ
ı
−ı
λ
−ı
ı
− 21 ı
1 5 (−3
1 2ı
+ 4ı)
1 5 (−3 −
4ı)
1
−1
1 2
−∞
0
−3
− 21
− 31
Table 5.2: Original and image coordinates of the Moebius transformation (5.38).
Inserting (5.38) into the characteristic equation (5.33), we have q3 λ3 + q2 λ2 + q1 λ + q0 = 0,
(5.39)
where q3 = p 0 + p 1 + p 2 + p 3 , q1 = 3p0 − p1 − p2 + 3p3
and
q2 = 3p0 + p1 − p2 − 3p3 ,
(5.40)
q 0 = p0 − p1 + p2 − p3 .
(5.41)
Fold bifurcation occurs if µ = 1, which is equivalent to λ → −∞. Having a look at (5.39) we see that if q 3 → 0 satisfies this condition, which gives s f old = 0 for k = 0. This condition is meaningless as the proportional gain must be positive to have a meaningful feedback control. The next case is the flip bifurcation, which occurs if µ = −1 or equivalently λ = 0. Assuming k = 0, this gives 4 cos √τs2 1 s f lip = √ √ 2γ τ s cos √τs2 − 2 sin The above equation has zero slope if τ s =
τs √ 2
.
(5.42)
√ 2πN (N is some integer), where the actual value of the function is 2/(Nπγ).
Thus we conclude that the maximum value of s f lip decreases as τ s increases, thus the maximal allowable proportional gain increases as the sampling period increases. The same expression (5.42) tends to infinity if τ s satisfies the following transcendental equation: τs τs √ = tan √ 2 2
→
τs =
√
2 × [0, 4.49341, 7.72525, 10.9041, 14.0662, . . .],
which means that the stability map is qualitatively same with a period of approximately π in the τ s direction. 94
(5.43)
5.4. Proportional control
The last and most interesting case of loosing stability is the Neimark-Sacker bifurcation (which is analogous to the Hopf case for continuous systems). Inserting λ = ıω into (5.39), we have q 0 q3 = q1 q2 , which gives √ 4d 2 + dγ sin2 √τs 2 . sNS = √ √ 2 + 2dγ sin( 2τ s ) − 2dγτ s The roots of the above equation are τ s =
equation
(5.44)
√ 2πN, with some integer N. The locations of ’vertical slopes’ satisfy the √
1
2τ s
1+
√
2 dγ
√ = sin( 2τ s ).
(5.45)
Fig.5.7 depicts the critical s values for α = 0.00467, γ = 0.23, d = 1 and k = 0. Red curves correspond to flip bifurcation µ = −1 given by (5.42) and blue curves represent Neimark-Sacker bifurcation, given by (5.44). It is very interesting that the effect of the proportional gain is counter-wise in the two cases. If the red curve is crossed upside down, one real multiplier crosses the complex unit circle at −1 while if the blue curve is crossed upwards, two complex multipliers cross the unit circle. Hence, the stable ’islands’ are in between the two curves, as depicted in Fig.5.7. Another interesting result is that if the control is assumed to be continuous, the condition for unconditional stability was that kd > s, which defines the minimally sufficient ’overall’ damping in the system for a given proportional gain. If no damping is present in the system, the control cannot be successful for any proportional gain s. However, the digital stability map shows that there are stable regions even for k = 0 and s = 0, thus the delay may have a stabilising effect in a narrow parameter range of s, τ s . This conclusion is somewhat surprising but similar phenomena was observed in the case of digitally controlled robots, see [77]. s 5 PSfrag replacements 4 5 3 stable 2
1
τs 5
10
15
20
Figure 5.7: Stability boundaries of digital control for α = 0.00467, γ = 0.23 and d = 1. Red curves correspond to flip bifurcation (birth of an oscillation with 2τ s period), while blue curves correspond to Neimark-Sacker bifurcation. (Compare this figure to the numerically computed one in Fig.5.17.)
Next the general stability maps for k , 0 and d , 0 are presented. These maps are computed numerically by simply evaluating the elements of the linear coefficient matrix T in (5.32). Although the explicit equations for the 95
5. Stability of controlled hydraulic systems
stability boundary curves for s f and sNS were computed, these formulae are so lengthy that it was much quicker and easier to simply calculate the characteristic multipliers numerically. Numerical experiments show that when varying k and d simultaneously, rather ’exotic’ stability charts may arise, thus the effect of viscous damping k and differential gain d are separately presented in Figure series Fig.5.8-5.13 and Fig.5.14-5.19, respectively. All these computations were performed with α = 0.0043 and γ = 0.23. As expected, the presence of damping in the system k , 0 stabilises the control, compare Fig.5.8 and Fig.5.9. Note that k does not denote the general nonlinear damping but the slope of the damping force at the given desired position. This means that for a slightly damped system k → 0 and for the Stribeck-type friction force we have k → ∞. Thus, parameter k was varied between 0 and 1000. √
Recall that for the continuous case (i.e. τ s = 0), we found that the condition for stable operation was s < k(d + 2/γ), which gives s < 6.149k with the actual parameter values. Note that the upper bound of the stable range in
Fig.5.8-5.13 matches exactly with this value. (This is not surprising but it is a nice verification of the correctness of the results.) Another interesting feature is that as k is increased, for moderate values of k (k = 1, Fig.5.11) the appearing instability is flip bifurcation, which means that the frequency of the oscillation is the half of the sampling frequency, thus these oscillations are extremely dangerous to the system. For large k values, say k = 10, 1000, Fig.5.12 and 5.13, the frequency of the oscillation born via the Niemark-Sacker bifurcation is typically low and is close to the eigenfrequency of the system. Corresponding to the effect of the differential gain, our first observation is that surprisingly 5, increasing the differential gain destabilises the system. Moreover, the main stricture of the stability boundary is around small proportional gains, where the system is expected to be more stable. For example, let us suppose that the our system works with τ s = 3 and s = 0.5. Having a look at Fig.5.14, we see that the control is stable but suppose that the overshoot is high and we want to decrease it, so a differential gain (’damping’) is added, d = 0.5. The result is that the formerly stable control becomes unstable and cannot fulfil its task, see Fig.5.16 with τ s = 3 and s = 0.5 as before.
5 Physically,
the differential gain can be represented as an artificial damping.
96
5.4. Proportional control
5
5 unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4
4
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
0
2
4
6
8
10
τ
12
14
16
18
unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
s
s
4.5
0
20
0
2
4
Figure 5.8: k = 0, d = 0
4
4
3.5
3.5
3
3
2.5
2.5
s
s
τ
12
14
16
18
20
2
2
1.5
1.5
1
1
0.5
0.5
0
2
4
6
8
10
τ
12
14
16
18
unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
0
20
0
2
4
Figure 5.9: k = 0.1, d = 0
6
8
10
τ
12
14
16
18
20
Figure 5.12: k = 10, d = 0
5
5 unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
4
4
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
2
4
6
8
10
τ
12
14
16
18
unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
s
s
10
5 unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
0
8
Figure 5.11: k = 1, d = 0
5
0
6
0
20
Figure 5.10: k = 0.5, d = 0
0
2
4
6
8
10
τ
12
14
Figure 5.13: k = 1000, d = 0
97
16
18
20
5. Stability of controlled hydraulic systems
5
5 unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4
4
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
0
2
4
6
8
10
τ
12
14
16
18
unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
s
s
4.5
0
20
0
2
4
Figure 5.14: k = 0, d = 0
4
4
3.5
3.5
3
3
2.5
2.5
s
s
τ
12
14
16
18
20
2
2
1.5
1.5
1
1
0.5
0.5
0
2
4
6
8
10
τ
12
14
16
18
unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
0
20
0
2
4
Figure 5.15: k = 0, d = 0.1
6
8
10
τ
12
14
16
18
20
Figure 5.18: k = 0, d = 10
5
5 unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
4
4
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
2
4
6
8
10
τ
12
14
16
18
unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
s
s
10
5 unstable − 0 stable µ unstable − 1 stable µ unstable − 2 stable µ stable − 3 stable µ
4.5
0
8
Figure 5.17: k = 0, d = 1
5
0
6
0
20
Figure 5.16: k = 0, d = 0.5
0
2
4
6
8
10
τ
12
14
Figure 5.19: k = 0, d = 1000
98
16
18
20
5.4. Proportional control
The effect of sampling delay was also studied experimentally on the experimental rig described in §2.2. Fig.5.20 and Fig.5.21 presents the calculated and experimentally determined stability chart, i.e. the maximal allowable proportional gain as a function of the sampling frequency. It should be emphasised from the very beginning that this is only a qualitative comparison as the theoretical model described in §5.2 omits several issues of the real system such as gravity, non-equal chamber pressures, non-constant pressure source, etc. However, the similarity of the two figures is striking. It should be noted that the position of the piston was given in terms of percents of the total stroke and the output signal was measured in volts, thus the dimension of the proportional gain was [V/%]. The sampling frequency was systematically increased from 1Hz to 50Hz and the stability boundary was then searched by increasing the proportional gain. The onset of instability was rather uncertain, thus for every sampling interval t s two S values are given; the blue line represent the ’always stable’ values and the red line depicts the ’always unstable’ values of the proportional gain. In between, the behaviour depends on secondary issues such as the initial and the desired position or the direction of the motion (upwards or downwards). The frequency of the appearing oscillation was found to be the half of the sampling frequency in the middle frequency range. For small frequencies (around t s = 0.5−1[s]), once the system began to oscillate, quickly the piston collided against the cylinder end stop and (due to safety reasons) the control was turned off. For high sampling frequencies, (above 20Hz), the frequency of the oscillation was very hard to identify and also the motion was
80
80
70
70
60
60
50
50
S [V/%]
S [V/%]
often not periodic.
40
40
30
30
20
20
10
10
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 0
0.1
0.2
ts [s]
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
t [s] s
Figure 5.20: Numerically computed critical proportional
Figure 5.21: Measured stability boundary, the stable re-
gain as a function of the sampling interval. The stable
gion is beneath the blue line. Above the red line the sys-
area is filled with red.
tem always oscillates, between the red and the blue line the behaviour is unpredictable (sometimes stable, sometimes unstable).
99
5. Stability of controlled hydraulic systems
5.5 Summary In this section the mathematical analysis of the PD control of a simple hydraulic positioning system was presented. In the first part, the control valve was an on/off valve, while in the second part a proportional valve was used for piston controlling. The reason for studying the on/off valve was that recently, parallel systems of such simple, cheap and robust valves are tried to use to replace expensive and sensitive servovalves. The problem with the analysis of such valves is that the resulting mathematical model is discontinuous and thus standard linear stability analysis cannot be performed. The first part of this chapter summarises the mathematical techniques that can be used for designing such systems and predicting the possible instabilities. The second part deals with the PD control via a proportional valve. Stability boundaries are analytically computed and based on the experiments, it is shown that digital sampling and digital control may lead to additional stability problems, whose solution is sometimes opposite to one’s engineering sense. The main results of this chapter are: 1. With the help of geometric singular perturbation theory, it was explained how the fluid compressibility leads to qualitatively different behaviour from the incompressible case in a simple hydraulic positioning system with relay feedback control. The analysis also reveals that at the onset of the instability only the ’fast’ pressure variables are affected, thus their dynamics should not be omitted when designing the control. 2. An analytical ’smoothing’ technique was applied to cope with the switching nature of the relay-like control. With the help of the technique, linear stability analysis can be performed and closed-formed results can be obtained for the critical parameter values. The correctness of the results was verified numerically with the software package AUTO. Also, it was shown how such a software can be used for quickly and efficiently find stability boundaries in such or more complex (realistic) systems. It was shown that the system is unconditionally stable if the controller satisfies s < dk, where s is the proportional gain, d is the differential gain and k is the linear viscous damping in the system. 3. The stability of the PD control via a proportional valve was studied, analytical formulae were given for the critical parameter values. Interestingly, similar to the previous case, unconditional stability is obtained if the controller satisfies s < dk, where s, d and k are as before. 4. The effect of digital sampling and digital control was studied and it was shown that (a) if correctly chosen, the delay can stabilise an a priori unstable control system, (b) as the sampling interval increases, the maximal allowable proportional gain (and thus the precision of the positioning) decreases and (c) increase in the differential gain easily destabilises the system. Experimental verification was also presented and a satisfactory qualitative agreement was found.
100
Chapter 6
Papers presenting the results of this thesis [1] Cs. H˝os and L. Kullmann. Sliding mode PD control of a simple hydraulic positioning systems (in Hungarian). In 12th International Conference in Mechanical Engineering (OGÉT 2004), pages 133–137, Csíksomlyó, April 2004. Hungarian Technical Scientific Society of Transylvania. [2] Cs. H˝os and L. Kullmann. Stability of the relay control of a simple hydraulic system. In Proc. of the Fourth Conference on Mechanical Engineering (Gépészet’04), pages 357–361, 2004. [3] Cs. H˝os and L. Kullmann. Theoretical and numerical study on the PD controlling of a simple hydraulic system. In C.R. Burrows, K.A. Edge, and D.N. Johnston, editors, Proc. of the Bath Workshop on Power Transmission and Motion Control, pages 211–221. Professional Engineering Publishing Limited, 2004. [4] Cs. H˝os and L. Kullmann. Variable time step methods for 1D hydraulic transmission line modelling. In Proc. of the 22nd IAHR Symposium on Hydraulic Machinery and Systems, pages B07–3, Stockholm, Sweden, June 2004. [5] Cs. H˝os and L. Kullmann. Dynamic modelling of a pilot-operated pressure relief valve. In C.R. Burrows, K.A. Edge, and D.N. Johnston, editors, Proc. of the Bath Workshop on Power Transmission and Motion Control. John Wiley and Sons, Ltd, 2005. accepted. [6] I. Vaik and Cs. H˝os. Computation of 1D unsteady pipe flow and the evaluation of the wall shear stress. In Proc. of the Fourth Conference on Mechanical Engineering (Gépészet’04), pages 332–336, 2004.
101
Chapter 7
Bibliography [7] R.T. Anderson and P.Y. Li. Mathematical modelling of a two spool flow control servoalve using a pressure control pilot. In ASME Symposium on Modelling and Control Electrohydraulic Systems, Orlando, FL., November 2000. [8] Ansys Inc. CFX 5.7 Users’ Guide, 2004. [9] K.J. Åström and T. Hagglund. The future of PID control. Control Engineering Practice, 9:1163–1175, 2001. [10] W. Backe.
What will be the future of fluid power.
Fluid Power WWW Virtual Library,
http://fluid.power.net/fpn/wwwvl/docs/rdpapers.php3, 2005. [11] G.K. Batchelor. An Introduction to Fluid Dynamics. Cambridge University Press, 2000. [12] J.M. Bergada and J. Watton. A direct solution for flowrate and force along a cone-seated poppet valve for laminar flow conditions. Proc. Instn Mech. Engrs, Journal of Systems and Control Engineering, 218:197–210, 2004. [13] A. Bergant, A.R. Simpson, and J. Vítkovsky. Developments in unsteady pipe flow friction modelling. Journal of Hydraulic Research, 39(3):249–255, 2001. [14] S. Bisztray-Balku. The effect of piston sealing friction and casing imperfections in hydraulic elements (in Hungarian). PhD thesis, Budapest University of Technology and Economics, 1984. [15] W. Borutzky, B. Barnard, and J. Thoma. An orifice model for laminar and turbulent conditions. Simulation Practice and Theory, 10:141–152, 2002. [16] K.K. Botros, G.H. Dunn, and J.A. Hrycyk. Riser-relief valve dynamic interactions. Journal of Fluids Engineering, 119:671–679, September 1997. [17] K.K. Botros, G.H. Dunn, and J.A. Hrycyk. Riser-relief valve dynamic interactions (extension to a previous model). Transactions of the ASME Journal of Pressure Vessel Technology, 120:207–212, May 1998. [18] B. Brunone, U.M. Golia, and M. Greco. Some remarks on the momentum equation for fast transients. Int. Meeting on Hydraulic Transients with Column separation, 9th Round Table, pages 140–148, 1991. [19] M.B. Bughazem and A. Anderson. Problems with simple models for damping in unsteady flow. In Int. Conf. on Pressure Surges and Fluid Transients, 1996. [20] L. Bálint. The effect of load change on the operation of a hydraulic cylinder (in Hungarian). PhD thesis, Budapest University of Technology and Economics, 1986. 102
BIBLIOGRAPHY
[21] C. Canudas de Wit, H. Olsson, K.J. Åström, and P. Lischinsky. A new model for control of systems with friction. IEEE Transactions on Automatic Control, 40(3):419–425, March 1995. [22] J. Chen, W.E. Dixon, J.R. Wagner, and D.M. Dawson. Exponential tracking control of a hydraulic proportional directional valve and cylinder via integrator backstepping. In Proceedings of IMECE’02. ASME International Mechanical Engineering Congress and Expo, November 2002. [23] M.H. Chiang, C.C. Chen, and T.N. Tsou. Large stroke and high precision pneumatic-piezoelectric hybrid positioning control using adaptive discrete variable structure control. Mechatronics, 15:523–545, 2005. [24] W.L. Daily, W.L. Hankey, R.W. Olive, and J.M. Jordaan. Resistance coefficients for accelerated and decelerated flows through smooth tubes and orifices. Trans. of ASME, 78:1071–1077, July 1956. [25] K. Dasgupta and R. Karmakar. Dynamic analysis of pilot operated pressure relief valve. Simulation Modelling Practice and Theory, 10:35–49, 2002. [26] W. D’Netto and D.S. Weaver. Divergence and limit-cycle oscillations in valves operating at small openings. Journal of fluids and structures, 1:3–18, 1987. [27] E.J. Doedel, T.F. Fairgrieve, B. Sandstede, A.R. Champneys, Y.A. Kuznetsov, and X. Wang. Auto97: Continuation and bifurcation software for ordinary differential equations. http://indy.cs.concordia.ca/auto, 1999. [28] W.S. Don and A. Solomonoff. Accuracy enhancement for higher derivatives using chebysheb collocation and a mapping technique. SIAM J. Scientific Computing, 18(4):1040–1055, July 1997. [29] R.D. Eyres, P.T. Piiroinen, A.R. Champneys, and N.A.J. Lieven. Grazing bifurcations and chaos in the dynamics of a hydraulic damper with relief valves. BCANM Preprint 2004.29. [30] N. Fenichel. Geometric singular perturbation-theory for ordinary differential-equations. Journal of differential equations, 31(1):53–98, 1979. [31] A.S. Garett and J.E. Bobrow. Experiments and simulations on the nonlinear control of a hydraulic servosystem. unpublished report, 1999. [32] E. Hairer, S.P. Nørsett, and G. Wanner. Solving Ordinary Differential Equations. Springer-Verlag, 1991. [33] G. Halász and A. Huba. Technical Measurements (in Hungarian). M˝uegyetem Kiadó, 2003. [34] G. Halász, G. Kristóf, and L. Kullmann. Fluid flow in pipe systems (in Hungarian). M˝uegyetem Kiadó, 2002. [35] S. Hayashi, T. Hayase, and T. Kurahashi. Chaos in a hydraulic control valve. Journal of fluids and structures, 11:693–716, 1997. [36] R.H.A. Hensen, M.J.G van de Molengraft, and M. Steinbuch. Friction induced hunting limit cycles: A comparison between the lugre and switching function model. Automatica, 39:2131–2137, 2003. [37] Cs. H˝os, A.R. Champneys, and L. Kullmann. Bifurcation analysis of surge and rotating stall in the Moore-Greitzer compression system. IMA Journal of Applied Mathematics, 68:1–24, 2003. [38] Cs. H˝os, L. Kullmann, and M. Lohász. On the friction factor in unsteady, laminar flow. In Third Conference on Mechanical Engineering (Gépészet 2002), volume 1, pages 352–356, 2002. 103
Bibliography
[39] I.E. Idelchic. Handbook of hydraulic resistance (in Russian). Economic and Energy Publishing, Moscow, 1960. [40] M. Jelali and A. Kroll. Hydraulic servo-systems - Modelling, identification and control. Springer, 2003. [41] T. Kagawa, I.-Y. Lee, A. Kitagawa, and T. Takenaka. High speed and accurate computing method for frequencydependent friction in laminar pipe flow for characteristics method. Trans. Jap. Soc. Mech. Engrs B, 49:2638–2644, 1993. [42] J.T. Karam and M.E. Franke. The frequency response of pneumatic lines. Trans. ASME, Journal of Basic Engineering, 89:371–378, 1967. [43] K. Kasai. On the stability of a poppet valve with elastic support. Bulletin of the JSME, 12(53):1094–1098, 1969. [44] T. Knohl and H. Unbehauen. Adaptive position control of electrohydraulic servo systems using ann. Mechatronics, 10:127–143, 2000. [45] G.A. Korn and T.M. Korn. Handbook of Mathematics for Engineers (in Hungarian). M˝uszaki Könyvkiadó, 1975. [46] L. Kovács. Advanced hydraulic machinery. Lecture notes on the Technical University of Budapest, 2000. [47] G.G. Kremer. Enhanced robust stability analysis of large hydraulic control systems via a bifurcation-based procedure. Journal of the Franklin Institute, 338:781–809, 2001. [48] L. Kullmann. Stability analysis of a hydraulic servomechanism (in hungarian). Master’s thesis, ELTE TTK, 1971. supervisor: M. Farkas. [49] L. Kullmann, S. Till, and F. Heged˝us. Comparison of vain material models for arterial blood flow simulations. In BUDAMED 2005 Conference for Biomedical and Clinical Engineers, Budapest, october 2005. [50] Y.A. Kuznetsov. Elements of Applied Bifurcation Theory. Springer, 1998. [51] T. Lajos. Basics of fluid mechanics (in Hungarian). M˝uegyetemi Kiadó, 2004. [52] R.P. Lambeck. Hydraulic Pumps and Motors. Marcel Dekker, Inc., 1983. [53] R.I. Leine, D.H. Van Campen, A. De Kraker, and L. Van Den Steen. Stick-slip vibrations induced by alternate friction models. Nonlinear Dynamics, 16:41–54, 1998. [54] R.I. Leine, D.H. Van Campen, and B.L. Van der Vrande. Bifurcations in nonlinear discontinuous systems. Nonlinear Dynamics, 23:105–164, 2000. [55] T.Y. Lin, Y.C. Pan, and C. Hsieh. Precision-limit positioning of direct drive systems with the existence of friction. Control Engineering Practice, 11:233–244, 2003. [56] Y. Liu and H. Handroos. Sliding mode control for a class of hydraulic position servo. Mechatronics, 9:111–123, 1999. [57] Y. Liu and H. Handroos. Sliding mode control of hydraulic position servo. Facta Universitatis, Mechanical Engineering, 1(9):1217–1230, 2002. [58] Gy. Ludvig. Dynamics of Machines (in Hungarian). M˝uszaki Könyvkiadó, 1983. [59] G. MacLeod. Safety valve dynamic stability: An analysis of chatter. Transactions of the ASME Journal of Pressure Vessel Technology, 107:172–177, May 1985. 104
BIBLIOGRAPHY
[60] H.E. Merritt. Hydraulic Control Systems. John Wiley and Sons, 1967. [61] M. Mihajolv, V. Nikoliˇc, and D. Antiˇc. Position control of an electro-hydraulic servo system using sliding mode control enhanced by fuzzy PI controller. Facta Universitatis, Mechanical Engineering, 1(9):1217–1230, 2002. [62] A.H. Nayfeh and H. Bouguerra. Non-linear response of a fluid valve. International Journal of Nonlinear Mechanics, 25:433–449, 1990. [63] P.M. Petherick and A.M. Birk. State-of-the-art review of pressure relief valve design, testing and modelling. Journal of pressure vessel technology-transactions of the ASME, 1991. [64] W.H. Press, S.A. Teukolsky, W.T. Wetterling, and B.P. Flannery. Numerical Recipes in C++ - The Art of Scientific Computing. Cambridge University Press, 2000. [65] B. Manhartsgruber R. Scheidl and G. Riha. Modern hydrostatic drive systems in the light of dynamical systems theory. In 4th Conference on Dynamical Systems - Theory and Applications. Lodz - Bronislawow, Poland, 8-9 December 1997. [66] K. Sato, K. Nakamoto, and A. Shimokohbe. Practical control of precision positioning mechanism with friction. Precision Engineering, 28:426–434, 2004. [67] H.A. Scarton and W.T Rouleau. Axisymmetric waves in compressible newtonian liquids contained in rigid tubes: steady-periodic mode shapes and dispersion by the method of eigenvalleys. Journal of Fluid Mechanics, 8(3):595– 621, 1973. [68] R. Scheidl, gestions
on
M. Garstenauer, the
future
B. Manhartsgruber,
development of
and G. Riha.
oil-hydraulics.
Fluid
Some comments and sugPower
WWW
Virtual
Library,
http://fluid.power.net/fpn/wwwvl/docs/rdpapers.php3, 2005. [69] R. Scheidl and B. Manhartsgruber. On the dynamic behaviour of servo-electric drives. Nonlinear Dynamics, 17:2147–268, 1998. [70] H. Schlichting. Boundary Layer Theory. McRaw-Hill Book Company, 1968. [71] G. Sheffel. Dynamisches verhalten eines direktgesteuerten kegelsitzventils unter dem einfluss der geometrie des schliessenelements. Ölhydraulik und Pneumatik, 22(5):280–448, 1978. [72] Y.C. Shin. Static and dynamic characteristics of a two stage pilot relief valve. ASME Journal of Dynamic Systems, Measurement, Control, 113(6):280–288, 1991. [73] J.J. Shu. A finite element model and electronic analogue of pipeline pressure transients with frequency-dependent friction. Trans. of the ASME, Journal of Fluid Engineering, 125:194–199, 2003. [74] G. Stepan and Gy. Haller. Micro-chaos in digital control. Journal of Nonlinear Science, 6(5):415–448, 1996. [75] J.A. Stone. Discharge coefficients and steady state flow forces for hydraulic poppet valves. ASME Journal of Basic Engineering, 82:144–152, 1960. [76] S.T.Tsai, A. Akers, and S.J. Lin. Modelling and dynamic evaluation of a two-stage two-spool servovalve. ASME Journal of Dynamic Systems, Measurement, Control, 113(12):709–713, 1991. [77] G. Stépán. Retarded Dynamical Systems: Stability and Characteristic Functions. John Wiley & Sons, 1989. 105
Bibliography
[78] G. Stépán. Nonlinear vibrations. Lecture Notes, 2000. Budapest University of Technology and Economics. [79] R. Szalai. Nonlinear dynamics of high-speed milling. Master’s thesis, Budapest University of Technology and Economics and University of Bristol, 2002. [80] T. Takenaka. Performances of oil hydraulic control valves. Bulletin of the JSME, 7(27):566–576, 1964. [81] L.N. Trefethen. Spectral Methods in Matlab. SIAM, 2000. [82] E. Urata. Thrust of a poppet valve. Bulletin of the JSME, 12(53):1099–1109, 1969. [83] I. Vaik and Cs. H˝os. Computation of 1d unsteady pipe flow and evaluation of the wall shear stress. In Fourth Conf. on Mechanical Engineering (Gépészet 2004), volume 1, pages 332–336, 2004. [84] A.E. Vardy and J.M.B. Brown. Transient, turbulent, smooth pipe flow. Journal of hydraulic Research, IAHR, 33(4):435–456, 1995. [85] A.E. Vardy and J.M.B. Brown. On turbulent, unsteady, smooth-pipe flow. In Int. Conf. on Pressure Surges and Fluid Transients, pages 289–311, Harrogate, England, 1996. BHR Group. [86] I. Vörös. Machine Elements III (in Hungarian). Tankönyvkiadó, 1956. [87] J. Watton. The stability and response of a two stage pressure rate controllable relief valve. Journal of fluid control, pages 55–60, 1990. [88] F.M. White. Viscous Fluid Flow. McGraw-Hill, 1991. [89] D.J. Wood and J.E. Funk. A boundary-layer theory for transient viscous losses in turbulent flow. Journal of basic Engineering, ASME, 92(4):865–873, 1970. [90] W. Zielke. Frequency-dependent friction in transient pipe flow. Journal of Basic Engineering, ASME, 90(1):109– 115, 1968. [91] P.S. Zung and M.H. Perng. Nonlinear dynamic model of a two-stage pressure relief valve for designers. ASME Journal of Dynamic Systems, Measurement, Control, 124:62–66, 2002.
106
