Multidiszciplináris tudományok, 3. kötet. (2013) 1. sz. pp. 97-106.
OPTIMALIZÁLT LÉPÉSKÖZŰ NEWTON-RAPHSON ALGORITMUS EHD FELADAT MEGOLDÁSÁHOZ Szávai Szabolcs
egyetemi adjunktus, Miskolci Egyetem, Anyagszerkezettani és Anyagtechnológiai Intézet, Szerkezetintegritási Intézeti Tanszék Cím: 3515 Miskolc, Miskolc-Egyetemváros, e-mail:
[email protected] Összefoglalás Csúszva-gördülő elemek, mint például a fogaskerekek, bütykös mechanizmusok és csapágyak gyakorta igen nagy terhelésnek, nagy sebességeknek és csúszásnak vannak kitéve, mikor nem csak a kenőanyagban kialakuló kontaktnyomás, hanem a felületi deformáció és a viszkozitás nyomásfüggése is kérdés. Az kontaktfelületek tribológiai viszonyainak elemzésére Dowson [2]megalkotta az általánosított Reynolds egyenletet. Azonban annak ellenére, hogy számos módszer került kifejlesztésre az EHD probléma vizsgálatára, az erősen nemlineáris feladat megoldása továbbra is kihívást jelent. A numerikus problémák kezelésére optimalizált lépésközű Newton-Raphson algoritmus került alkalmazásra a nem-lineáris egyenletrendszer megoldásához. A javasolt eljárás csökkenti a lépések számát és stabilizálja az megoldáskeresést. Kulcsszavak: optimalizált lépésköz, EHD feladat, tribilógia, Reynolds egyenlet Abstract Rolling-sliding machines such as gears, cams and followers, and bearings, which are often subjected to high loads, high speeds and high slip conditions when not only the pressure distribution in the lubricant is a question but the surface deformation, and variation of the viscosity due to pressure. To analyse the tribological conditions of the contact surfaces, the generalized Reynolds equation was developed by Dowson [2]. Although several methods have been already developed for solving EHD problems, the solution of the highly nonlinear problem is still quite challenging. Handling the numerical problems during the solution quadratic step optimized Newton-Raphson method has been implemented for solving the nonlinear equation system. The proposed technique reduces the step number during the iteration and the solution is stabilized. Keywords: quadratic step optimization, EHD, tribology, Reynolds equation
1. Bevezetés Az 1. ábra mutatja a folyadéksúrlódás állapotában lévő foltszerűen érintkező felületpárok általánosított esetét. A testek egymáshoz viszonyított relatív elmozdulásának következtében a testek közti rést kenőanyag tölti ki, mely mozgásának hatására hidrodinamikai nyomás alakul ki. A kenőanyag mozgását a felületek egymáshoz viszonyított relatív mozgásának hatására a kenőanyagban fellépő csúsztatófeszültség váltja ki. Az érintkező testek kinematikai állapota és egy adott résgeometria mellett, az érintkező felületeken fellépő nyomás-
Szávai, Sz. megoszlás képes egyensúlyt tartani a felületeket összeszorító erővel, megakadályozva a test-test kapcsolatot. Adott esetben a felületeket terhelő nyomáseloszlás, illetve a kenőanyagban fellépő csúsztatófeszültségek hatására képződő hődisszipáció okozta, lokális vagy globális hőmérsékletnövekedés akkorává válhat, hogy a felületek figyelmen kívül nem hagyható deformációját okozhatja, továbbá kihathat a kenőanyag anyagjellemzőire. Látható tehát, hogy amennyiben termo-elasztohidrodinamikai kenési viszonyok közt kívánjuk modellezni az érintkezés során kialakuló körülményeket, egyszerre kapcsoltan kell megoldanunk hidrodinamikai, termodinamikai és szilárdságtani problémát, melyek már önmagukban is, de a különböző kontinuumok anyagjellemzőinek állapotfüggése miatt is erősen nemlineáris rendszert alkotnak. 2. test
2. felület (S2)
z
u2
y h2
Kontakt zóna (Ac)
h1
x
u1 1. felület (S1)
Kontakt zóna határa (c)
1. test
1. ábra. Érintkezési feladat folyadékkenés esetén A brit O. Reynolds 1886-ban közzétett megoldásával elérte, hogy adott résgeometriánál (adott h(x,y)), a térbeli sebességmező és nyomáseloszlás helyett kenéselméleti feladatok megoldásához, elegendő egy résvastagság mentén vett átlagnyomást meghatározni ( p( x, y) ), mely alkalmas a résben lezajló áramlástani jelenség leírására. 1961-ben Dowson és Higginson megalkotta newtoni kenőanyagokra az általánosított Reynolds egyenletet, melyben figyelembe vették a résmenti hőmérsékletkülönbség okozta viszkozitás- és sűrűségváltozást. Az általánosított nem-newtoni Reynolds egyenletet Najji, Bou-Said és Berthe munkáját követve Wolff és Kubo alkotta meg. A kenési problémák leírása során kezelni kell továbbá a folyadékfilm keletkezésének és megszűnésének problémáját is, mely a korábbi egyenletek kavitációs zónára való kiterjesztését teszi szükségessé. A megoldás alapját képező kiterjesztett egyenletek részletezett formában [1]-ben találhatók meg.
98
Optimalizált lépésközű Newton-Raphson algoritmus EHD feladat megoldásához Az [1]-ben található egyenlet a Reynolds által bevezetett, illetve a turbulens és az inercia tagokon kívül, egyéb elhanyagolásokat nem tartalmaz. Ennek formális alakja a következő: R p, , , h, r , t , , eq , u1, u2 ..... xy xy xy p 0
(1)
ahol ( , h, r, t , , eq ....), ( , h, r, t , , eq , u1 , u2 ....) és ( , , h, r, t,W1 ,W2 ....) töb bek között a sűrűséget, h résméretet, az r helyvektort, t időt, viszkozitást és u i Wi felületi sebességeket, nem Newtoni kenőanyagoknál eq egyenértékű csúsztatófeszültséget tartalmazó belső függvények, melyek a Reynolds egyenlet tagjait adják vissza. A Reynolds egyenletet néhány egyszerű esettől eltekintve nem lehet zárt alakban megoldani, nem is szólva annak általánosított alakjáról. Így szükségessé váltak a numerikus módszerekre épülő megoldások, melyek során az váltókat vagy diszkrét pontokban keressük vagy valamilyen approximációs függvény segítségével közelítjük. Ennek következtében az eredeti egyenleteink helyett, változónkként a megoldás keresésére kijelölt pontok vagy az approximáció ismeretlen paramétereinek számával egyező számú egyenletet kapunk. Így a megoldás a kapott egyenletrendszer megoldásával közelítjük.
2. Optimalizált lépésközű Newton-Raphson algoritmus A vizsgált tárgykörben a megoldás nehézségét különösen a Reynolds egyenlet erősen nemlineáris tulajdonsága okozza. Az erősen nemlineáris egyenletek megoldásai az esetek túlnyomó többségében a Newton vagy gradiens módszerre alapulnak Ebből adódóan a diszkretizált Reynolds egyenletet linearizálni kell, hogy a nyomáseloszlás és a hozzá tartozó résalak meghatározható legyen. Tekintsük az egyenletet mint: R R P, , , h, r , t , , eq ,, , u1 , u2 ..... (2)
A megoldás keresése során a (2) változóinak olyan értékét keressük, melyre a R maradványfüggvény értéke 0. Az egyenlet változói nem függetlenek egymástól, hanem azokat különböző egyenletek, mint pl. az anyagegyenlet kapcsolja össze. Ugyanakkor az esetek jelentős részénél a változók közti összefüggéseket nem lehet explicit alakban felírni. Ennek következtében a numerikus megoldás során a változóknak kezdőértéket kell adni, majd a feladatot egy kiválasztott változóra (esetünkben P) nézve meg kell oldani, míg a többi változó (eq, h, …) rögzítve marad. Ezt követően a kiválasztott változó új értékéhez (az új Phez) kell meghatározni a többi változó értékét, melyek a következő iterációs lépés új kiinduló értékeként fognak szerepelni. A Reynolds egyenlet változóinak áttekintése során megállapíthatjuk, hogy a sűrűség, a viszkozitás, a kitöltési tényező és a lineárisan rugalmas anyagmodell mellett a felületek deformációból származó résméret esetén a nyomással és a hőmérséklettel való kapcsolat explicit alakban felírható, így a P-re vett linearizálás során ezek deriváltjait is figyelembe tudjuk venni, mely nagymértékben gyorsíthatja a megoldást és párhuzamosan kerül meghatározásra a nyomáseloszlás és az ehhez tartozó résméret. i
R R P, ( p(P)), h(P), ( p(P)), ( p(P)), i1 R
(3)
99
Szávai, Sz. Ennek linearizált formája egy tetszőleges P=Pj pontban: i
i R i R i R i R i R P Ο 0 (4) R P j ,i 1 R Np Np NhD p Np P p p h p PP j
Ezen egyenlet Pj megoldásnak sorozatával közelítjük az iR=0 maradvány iP*
megoldását a következő szerint:
P j 1 P j P j
[0..1]
(5)
Mivel az EHD feladatok esetén a kiinduló egyenletrendszer többszörösen nemlineáris, igen nehéz olyan kiinduló állapotot találni, mely esetén elkerülhető a megoldás oszcillációja. Ezért egy célszerűen megválasztott α csillapításra van szükségünk, amelyik a leggyorsabb konvergenciát eredményezi. Ennek meghatározására viszont csak általános megfontolások állnak rendelkezésre, melyek adott feladathoz való megfelelőssége nem biztosított. Ezért a csillapítás optimális mértékét úgy érdemes meghatározni, hogy az a maradvány értékének a lehető legnagyobb mértékű csökkenését eredményezze. Természetesen ezt csak egy eljárással lehet biztosítani. Az csillapítás optimális értékének meghatározásához a maradvány négyzet minimuj j mát keresem. Az csillapítás meghatározásakor P és P vektorok állandóak, így az j j R(P +P ) maradványvektor egyváltozós függvény, melyet továbbiakban R()val jelölök. Annak érdekében, hogy ne legyen szükség a minimumkeresés során a deriváltak időigényes meghatározására, a maradvány értékét másodfokú approximációval közelítem a következő módon: Legyen R0 az =0 -hoz tartozó maradvány Kezdő értékként legyen 1=0,6 melyhez tartozó maradvány R1 Az 2 érték meghatározásánál rendelkezésre áll két pontban ( 0=0 és 1 pontban) a maradvány R0 és R1 értéke és természetesen azok (R0)2= R0R0 és (R1)2= R1R1 négyzetössze2 ge. Valamint a kiinduló pontban a maradvány deriváltjait is ismerjük, így a (R())
maradványnégyzet függvény meredeksége is ismert, melyet jelöljük „m”-mel. Az előbbi adatokkal (R())2 másodfokú közelítése alapján: 2
m 12 1 2 R 1 2 R 0 2 m 1
(6)
Amennyiben az m értéke valamilyen oknál fogva nem állna rendelkezésre, az [0..1] tartományon a következő súlyozott átlaggal a következőképpen is felvehető 2:
2
100
1 R 0 R 0 0 R1 R1
R 0 R 0 R1 R1
(7)
Optimalizált lépésközű Newton-Raphson algoritmus EHD feladat megoldásához
R2() R0 2
R2
1
2
R1 2
m
2
1
2. ábra. Maradványnégyzet első közelítése
R2() R0 2
R1 2 1
Ri-12 m Ri 2 i
i-1
1
3. ábra. Első lépéssorral kapott három pont Ehhez a csillapításhoz tartozó maradvány érték legyen R2. Így három -R értékpár lehetőséget ad az R() maradványérték R2()=RR négyzetének másodfokú közelítéséhez, mely mellett továbbra is rendelkezésre áll a kiinduló pontban az (R())2 maradványnégyzet függvény meredeksége. A parabolikus közelítéshez azonban az egyik adat felesleges Abban az esetben, ha 2=1± fennáll, ahol egy megválasztott elegendően kicsi hibahatár, az 2 az optimális lépéscsillapításnak tekinthető. Abban az esetben, ha 2<1 és (R0)2 <(R2)2 fennáll az 3 értékét az [0;(R0)2] és [2;(R2)2] pontok valamint az [0;(R0)2] pontbeli m érték alapján határozzuk meg. Tesszük ezt mindaddig, míg i<i-1 és (Ri)2 <(R0)2 nem teljesül (3. ábra).
101
Szávai, Sz. Legyen ekkor:
1=i-1, (R1)2=(Ri-1)2 és 2=i, (R2)2=(Ri)2. A fenti eljárás során mindenképpen előáll egy olyan [0;(R0)2], [1;(R1)2], [2;(R2)2] ponthármas, ahol (R1)2 és (R2)2 értékek közül legalább egy kisebb, mint (R0)2. Ezután a soron következő i értékét már két célszerűen megválasztott, a korábbi pontok közül a két legkisebb, [k;(Rk)2],[l;(R l)2] és az [i-1;(R
i-1)
2
] pontra fektetett
2
c j
j
j 0
alakú másodfokú interpolációval előállított parabola minimumhelye, vagy tengellyel vett metszéspontja adja (4. ábra -6. ábra). Belátható, hogy a rendelkezésre álló ponthalmazból ez három egymás melletti pont, az interpoláció k, l és i-1 sorrendjétől függetlenül, 3 különböző esetet jelöl ki, melyeket a 4. ábra, a 5. ábra és a 6. ábra mutat. Speciális eset az, amikor a 3 pontban ugyanakkora a maradványnégyzet: (Rk)2=(R l)2=(R i-1)2. Ekkor legyen i=(i-1+k)/2.
2 A 4. ábra és az 5. ábra, által mutatott esetekben, mikor a d 2 ci i d 2 0 az 3 i 0 értékét a másodfokú interpolációs görbe minimumpontja jelöli ki.
R2() Rk 2 Ri-12
Rl 2
k
l
i
i-1
4. ábra. Három ponton keresztül a maradványnégyzet közelítésének 1. esete
102
Optimalizált lépésközű Newton-Raphson algoritmus EHD feladat megoldásához
R2() Rk2
Rl2 Ri-12
k
l
i-1
i
5. ábra. Három ponton keresztül a maradványnégyzet közelítésének 2. esete
R2() Rk2
Rl2 Ri-12
k
l
i-1
i
6. ábra. Három ponton keresztül a maradványnégyzet közelítésének 3. esete
R2() Rm2 Rk2
Rl2 Ri-12
k
l
i-1
mi i’
7. ábra. Három ponton keresztül a maradványnégyzet közelítésének 4. esete 103
Szávai, Sz. Ha a 6. ábra által mutatott eset áll fenn, vagy a három pont egy egyenesbe esik, azaz 2 d 2 ci i d 2 0 , az i értékét a görbe tengellyel vett metszéspontja adja, ahol i 0
2
c i
i
0 . Mivel a görbének két metszéspontja van, az lesz az i, amelyik nagyobb,
i 0
mint 0 és közelebb áll az i-1-hez. Az 5. ábra. és 6. ábra által mutatott esetben a meghatározott i kívül esik a három pont által felvett tartományon. Ekkor előfordulhat, hogy a tartomány határa és az új i érték közé esik egy korábban felvett m érték, ahol már ismert az (Rm)2 értéke. Ebben az esetben, a meghatározott i helyett az i=m –t veszünk, ahogy azt a 7. ábra is mutatja. A fenti eljárással viszonylag gyorsan megtalálható az az érték, ami mellett a megoldás felé vezető, optimális lépés tehető. Mivel a lépés iránya kötött, az értékét nem kell nagy pontossággal meghatározni, az optimalizált lépés hossz 15%-os pontosságú meghatározása elegendő a kívánt hatás eléréséhez.
3. Egy EHD feladat megoldása Az eljárás hatékonyságának bemutatását a Houpert, L. G. és Hamrock, B. J. [3] által, 1986ban közölt cikkben található példán keresztül mutatom be. A megoldás során a kidolgozott optimalizált lépésközű Newton-Rapshon került alkalmazásra. Az algoritmus hatékonyságának szemléltetésére két kiragadott lépésben látható a hibanégyzet alakulása az optimális lépésköz meghatározása során (8. ábra és a 9. ábra). A számítás elején az algoritmus megakadályozza az iteráció elszállását, míg a későbbiekben hatékonyan gyorsítja a megoldás menetét.
8. ábra. R2-α értékek lépésenként a megoldás későbbi fázisában 104
Optimalizált lépésközű Newton-Raphson algoritmus EHD feladat megoldásához
9. ábra. R2-α értékek lépésenként a megoldás későbbi fázisában A nyomáseloszlás megoldását a 10. ábra mutatja az irodalomban közölt eredményekkel együtt. Összehasonlítva megállapítható, hogy az eredmények jó egyezést mutatnak, különösen a Houpert, L. G. és Hamrock, B. J.[3], 1986-ban közölt eredményével. Hamrock, B. J. és Jacobson, Bo O. Houpert, L. G. és Hamrock, B. J., p-FEM solutions
p (MPa)
2
h (µm10 )
x (mm) 10. ábra. A megoldás során és az irodalomban közölt kialakuló nyomáseloszlás
105
Szávai, Sz.
4. Összefoglalás A mintafeladat megoldása során a teljes Newton-Raphson helyett, annak optimalizált lépésű formáját használom abban az értelemben, hogy a kiinduló paraméter kombinációból a Newton-Raphson megoldás első lépéseként előálló paraméter kombináció irányába elindulva kerestem azt a paraméterkombinációt, mely a legkisebb maradvány értéket szolgáltatja. A lépésköz-optimalizálás a numerikus elaszto-hidrodinamikai feladat Newton-Raphson megoldása során jelentősen csökkenti a N-R lépések számát, illetve a megoldás oszcillációját. A szükséges deriváltak előállításának nagy számításigénye jelentős időt is igényel. Az eljárással minimalizálható az N-R lépések száma, így a deriváltak előállításának gyakorisága is és ezzel együtt a megoldáshoz szükséges idő is. Ezek alapján az optimalizált lépésközű Newton-Rraphson algoritmus kifejezetten hatékony eszköze az erősen nemlineáris feladatok megoldásának.
5. Köszönetnyilvánítás A kutatás a TÁMOP 4.2.4.A/2-11-1-2012-0001 azonosító számú Nemzeti Kiválóság Program – Hazai hallgatói, illetve kutatói személyi támogatást biztosító rendszer kidolgozása és működtetése országos program című kiemelt projekt keretében zajlott. A projekt az Európai Unió támogatásával, az Európai Szociális Alap társfinanszírozásával valósul meg.
6. Irodalom [1] [2] [3]
106
Szávai, Sz.(2009) Efficient p-version FEM Solution for TEHD Problems with New Penalty-Parameter Based Cavitation Model, BALTTRIB’2009, Kaunas, 2009. pp 194-199. Dowson, D., (1962) A Generalised Reynolds Equation for Fluid-Film Lubrication, Int. Journal of Mechanical Science, Pergamon Press. Houpert, L. G. and Hamrock, B. J., (1986) A Fast Approach for Calculating Film Thickness and Pressure in Elastohydrodynamically Lubricated Contacts at High Loads, ASME Journal of Tribology, Vol. 108, No. 3, pp. 411-420.
