MISKOLCI EGYETEM GÉPÉSZMÉRNÖKI KAR ÁRAMLÁS- ÉS HŐTECHNIKAI GÉPEK TANSZÉKE

MISKOLCI EGYETEM GÉPÉSZMÉRNÖKI KAR ÁRAMLÁS- ÉS HŐTECHNIKAI GÉPEK TANSZÉKE

KÉTDIMENZIÓS NYĺRÓÁRAMLÁSOK SZÁMĺTÁSA A TURBULENS ÖRVÉNYDIFFÚZIÓ DIFFERENCIÁLEGYENLETÉNEK MEGOLDÁSÁVAL Ph.D. ÉRTEKEZÉS

Készítette: KÖNÖZSY LÁSZLÓ okl. gépészmérnök SÁLYI ISTVÁN GÉPÉSZETI TUDOMÁNYOK DOKTORI ISKOLA, GÉPÉSZETI ALAPTUDOMÁNYOK SZAKTERÜLET, TRANSZPORTFOLYAMATOK ÉS GÉPEIK TÉMACSOPORT

A doktori iskola vezetője: DR. PÁCZELT ISTVÁN Az MTA rendes tagja Témavezető: DR. CZIBERE TIBOR Az MTA rendes tagja

MISKOLC-EGYETEMVÁROS, 2004.

TARTALOM 1. BEVEZETÉS .............................................................................................................................................. 2 1.1. AZ IRODALOM ÁTTEKINTÉSE ............................................................................................................. 3 1.2. A FELADAT MEGFOGALMAZÁSA ..................................................................................................... 11 2. A NUMERIKUS ELJÁRÁS ELMÉLETI HÁTTERE............................................................................. 12 2.1. A TURBULENS MOZGÁSOK ALAPEGYENLETEI.............................................................................. 12 2.2. AZ ALKALMAZOTT TURBULENCIA-MODELL ................................................................................. 15 3. A NUMERIKUS ELJÁRÁS...................................................................................................................... 18 3.1. A FELADAT PARCIÁLIS DIFFERENCIÁLEGYENLET-RENDSZERE................................................ 18 3.2. A DISZKRETIZÁCIÓ ÉS A NUMERIKUS MEGOLDÁS ....................................................................... 24 4. A NUMERIKUS ELJÁRÁSSAL SZÁMÍTOTT EREDMÉNYEK ......................................................... 33 4.1. TURBULENS ÁRAMLÁS KÖRKERESZTMETSZETŰ CSŐBEN.......................................................... 33 4.2. TURBULENS ÁRAMLÁS KOAXIÁLIS FORGÁSFELÜLETEK ÁLTAL HATÁROLT FORGÁSSZIMMETRIKUS TÉRBEN...................................................................................................... 33 4.3. A SÍKÁRAMLÁS SPECIÁLIS ESETE.................................................................................................... 33 4.4. A SZÁMÍTÁSI EREDMÉNYEK ÉRTÉKELÉSE .................................................................................... 54 5. KÉTDIMENZIÓS MÁSODRENDBEN FOLYTONOS ORTOGONÁLIS GÖRBEVONALÚ KOORDINÁTAHÁLÓZAT LÉTREHOZÁSA ........................................................ 55 5.1. A NUMERIKUS ELJÁRÁS ALGORITMUSA......................................................................................... 55 5.2. A NYÍRI-FÉLE SIMÍTÓ ELJÁRÁS ALKALMAZÁSA........................................................................... 55 5.3. AZ ELSŐ ÉS MÁSODIK DERIVÁLTAK SZÁMÍTÁSA A DISZKRÉT PONTOKBAN ......................... 57 5.4. A HERMITE-POLINOM ALKALMAZÁSA............................................................................................ 58 5.5. KÉT KONTÚR KÖZÉ BEÍRHATÓ KÖR KÖZÉPPONTJÁNAK SZÁMÍTÁSA ...................................... 59 5.6. AZ ORTOGONÁLIS TRAJEKTÓRIÁK MEGHATÁROZÁSA ............................................................... 61 5.7. A SZÁMÍTÁSI TAPASZTALATOK ÖSSZEFOGLALÁSA .................................................................... 62 6. ÖSSZEFOGLALÁS .................................................................................................................................. 63 7. ÚJ TUDOMÁNYOS EREDMÉNYEK..................................................................................................... 64 FONTOSABB JELÖLÉSEK........................................................................................................................ 65 FÜGGELÉK ................................................................................................................................................. 67 IRODALOMJEGYZÉK ............................................................................................................................... 76 AZ ÉRTEKEZÉS TÉMÁJÁBAN MEGJELENT TUDOMÁNYOS PUBLIKÁCIÓK....................................... 82 AZ ÉRTEKEZÉS TÉMÁJÁVAL KAPCSOLATOS KÉZIRATOK ÉS DOKUMENTÁCIÓK ......................... 84

1

1. BEVEZETÉS A turbulens áramlásokkal kapcsolatos elméleti- és kísérleti kutatások a Miskolci Egyetem Áramlás- és Hőtechnikai Gépek Tanszékén több mint egy évtizede folynak. Az elméleti kutatások alapját a Dr. Czibere Tibor által kidolgozott új sztochasztikus turbulencia-modell jelenti, amely a Kármán-féle kétdimenziós hasonlósági hipotézis háromdimenziós esetre történő kiterjesztése. Az értekezés tárgya numerikus módszer kidolgozása két egymást nem metsző koaxiális forgásfelület által határolt forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens áramlás sebességterének meghatározására. A feladatot a középsebességtérben érvényes turbulens örvénydiffúzió új, Czibere-féle parciális differenciálegyenletének segítségével oldjuk meg, ami egy új differenciálegyenlet numerikus megoldását jelenti. A kidolgozott numerikus eljárást koaxiális forgásfelületek által határolt térben kialakuló perdületmentes turbulens áramlások esetén teszteltük. A számítások hitelességét az analitikus megoldásokkal és a mérési eredményekkel való összehasonlítás igazolja. A parciális differenciálegyenletek numerikus megoldása során diszkrét pontokból álló geometriai adatrendszer előállítása szükséges. A Dr. Nyíri András simító eljárását felhasználva egy új kétdimenziós másodrendben folytonos ortogonális görbevonalú koordinátahálózat generáló eljárást dolgoztunk ki. A számítások eredményeképpen a meridiánmetszet diszkrét pontjainak koordinátái és az ívelemek hossza mellett a pontokbeli első és második deriváltak által a koordinátavonalak görbületei is numerikusan ismertek lesznek. Az értekezés 1. fejezetében a vonatkozó irodalom áttekintése után a feladat megfogalmazása következik. A 2. fejezet a megoldás során alkalmazott turbulencia-modell elméleti hátterét ismerteti. A 3. fejezet az áramlási feladat alapegyenleteit és a megoldásra kidolgozott numerikus eljárást tárgyalja. A 4. fejezet a számítási eredményeket tartalmazza. Az 5. fejezet a kétdimenziós másodrendben folytonos ortogonális görbevonalú koordinátahálózat generáló eljárást ismerteti. A 6. fejezet az értekezés összegzése. A 7. fejezet az új tudományos eredményeket foglalja magában. A Függelék az értekezés főbb pontjainak kiegészítését tartalmazza.

2

1.1. AZ IRODALOM ÁTTEKINTÉSE A műszaki gyakorlatban legtöbbször turbulens áramlásokkal találkozunk, amelyet elsőként Osborne Reynolds angol fizikus-mérnök vizsgált, és kísérletekkel kimutatta, hogy egy kritikus áramlási sebesség felett minden áramlás turbulenssé válik, amelyben a viszkózus folyadék részecskéi egymással állandóan keverednek. A folyadékot egy átlátszó csőben megfestette, és a kritikus áramlási sebesség felett a folyadékrészek elkeveredtek a vizsgálat kiindulópontjául szolgáló festetlen folyadékkal. A kísérletek során kiderült, hogy minél nagyobb a csőátmérő, annál kisebb az a kritikus sebességérték, amelynél fellép a turbulencia jelensége, és az áramlási tér bármely pontjában a sebesség- és a nyomásértékek időben gyorsan, nagy frekvenciával változnak. Reynolds (1895) [1] a kísérletei alapján feltételezte, hogy a) a turbulens mozgásra is érvényes a Navier-Stokes-féle mozgásegyenlet; b) a különböző méretű csövekben áramló folyadékok turbulens áramlás esetén is mechanikailag hasonlóan mozognak. J. Boussinesq (1877) [2] a newtoni folyadék nyírófeszültségének analógiája alapján bevezette az örvényviszkozitás fogalmát, miszerint a turbulens nyírófeszültség a ν t örvényviszkozitás és a turbulens középsebesség deformációsebességének szorzata. Ludwig Prandtl (1904) [3] kísérleti eredményekből kiindulva kidolgozta a határrétegelméletet, miszerint kis viszkozitású folyadékok áramlása esetén a) a közeg belső súrlódását csak a szilárd falak közelében, az ún. határrétegben kell figyelembe venni; b) a határrétegen kívül az áramlási tér bármely pontjában a súrlódási erők elhanyagolhatók, ezért az áramlás potenciálosnak tekinthető. A határrétegen belül az áramlás hasonló a viszkózus folyadék csőben való áramlásához, ezért megkülönböztethetünk lamináris és turbulens határréteget, amelyek között az átmenet egy adott szakaszon belül változó helyen megy végbe. Prandtl a helytől függő keveredési úthossz fogalmából kiindulva meghatározta a körkeresztmetszetű csőben kialakuló turbulens áramlás sebességeloszlását, és feltételezése szerint a keveredési úthossz a faltól mért távolsággal egyenesen arányos. Prandtl feltevéseit a síkfal- és a körkeresztmetszetű cső falának közelében kialakuló turbulens áramlás esetén a kísérletek jól igazolták. Kármán Tódor (1930) [4] kísérleti tapasztalatok alapján felállította a kétdimenziós hasonlósági hipotézisét, miszerint a) a falközeli viszkózus alaprétegtől eltekintve a turbulens sebességeloszlás független a közeg viszkozitásától; b) teljesen kifejlődött turbulens áramlás esetén a tér különböző pontjaiban kialakuló áramképek mechanikai hasonlóságot mutatnak. Kármán szerint a keveredési úthossz az áramlás irányába eső sebességkomponens normális irányú első deriváltjának abszolút értékével egyenesen és a normális irányú második deriváltjának abszolút értékével fordítottan arányos. Az elvégzett kísérletek szerint a körkeresztmetszetű csőben kialakuló turbulens áramlás sebességeloszlásának Kármán-féle megoldása a cső egész keresztmetszetében a mérésekkel jól megegyező értékeket mutat. G.I. Taylor (1935) [5] a turbulens sebességingadozást valószínűségi változónak tekintette, amely a pillanatnyi sebesség és az időbeli átlagsebesség különbsége. Taylor feltételezte, hogy a turbulencia jelensége homogén és izotrop. A turbulens áramlási folyamatok korrelációs tenzorokkal történő leírása esetén a megoldandó differenciálegyenletek száma minden esetben kevesebb, mint az ismeretlenek száma, ezért a problémát a Reynolds-féle turbulens látszólagos feszültségtenzor meghatározása jelenti, vagyis a turbulens áramlási folyamatok matematikai modellezése szükségessé válik. A. Kolmogorov (1941) [6] bevezette a k ún. turbulens kinetikus energia fogalmát, és felfogása szerint ennek transzportja minden turbulens áramlásban jelen van a nagyobb örvényektől a kisebb örvények felé haladva. A turbulens áramlásban az egymással érintkező folyadékrészecskék között folyamatos impulzuscsere valósul meg és a kialakuló örvények a folya3

dék energiájának jelentős részét felemésztik. Kolmogorov számára indokolt volt a hossz-, az idő- és a sebességléptékek bevezetése kapcsán az ε = − dk/dt turbulens kinetikus energiadisszipáció mennyiségének definiálása is, ami azt a sebességet fejezi ki, amellyel a nagyobb örvények a kisebb örvények irányában leadják az energiájukat [7]. A Reynolds-féle látszólagos turbulens feszültségek Smagorinsky-féle (1963) [8] közelítése csak kétdimenziós turbulens áramlások vizsgálata során vezetett eredményre kicsiny osztásközű diszkrét halózat alkalmazása esetén. B. Launder és D. Spalding (1972) [9] kidolgozták a napjainkban széleskörben elterjedt k-ε turbulencia-modellt, amelyben a turbulens áramlást leíró mozgásegyenletek kiegészülnek a k turbulens kinetikus energia és az ε turbulens kinetikus energiadisszipáció transzportját reprezentáló differenciálegyenletekkel. A k-ε modell a Boussinesq-féle örvényviszkozitási hipotézisen alapuló két-egyenlet modell, mert az örvényviszkozitást a k és az ε mennyiségek segítségével definiálják. A k-ε modellnek számtalan változata született, amelyek alkalmazásaira vonatkozóan bőséges irodalom áll rendelkezésre [10], [11]. M. Griebel, T. Dornseifer és T. Neunhoeffer (1998) [12] német szerzők a műszaki gyakorlatban elterjedt turbulencia modelleket az 1. ábra szerint csoportosították.

Turbulens Áramlás Szimulációja

Direkt Numerikus Szimuláció(DNS)

Turbulencia Modellek

Finom Struktúra Modellek

Algebrai Modellek (Null-Egyenlet Modellek)

Statisztikai Modellek

Örvényviszkozitási Modellek

Feszültség Transzport Modellek

Egy-Egyenlet Modellek

Két-Egyenlet Modellek (k-epszilon)

Standard k-epszilon

Kis Reynolds-Szám Modellek

Két-Réteg Modellek

k-epszilon RNG Módszerek

1. ábra A műszaki gyakorlatban elterjedt turbulencia-modellek csoportosítása

4

A turbulencia-kutatás a mérnöki alapkutatások egyike, ezért a turbulencia jelenségével kapcsolatos irodalom meglehetősen bőséges. Hans Breuer vizsgálata [13] szerint 1986-ban a laboratóriumok és a kutatóintézetek természettudományos eredményeit több mint 10000 folyóirat publikálta és a közlemények számának éves növekedési üteme 4-6% volt. Csermely Péter, Gergely Pál, Koltay Tibor, Tóth János a "Kutatás és Közlés a Természettudományokban" [14] című átfogó munkája napjaink szakirodalom-kutatásával foglalkozik. Az 1999-ben végzett kimutatások szerint a természettudományokban évente átlagosan 40-50 ezer tudományos folyóirat és 200000 szakkönyv jelenik meg. Itt fontos megjegyeznünk, hogy a XX. századi számítástechnika fejlődése az információ robbanásszerű terjedését tette lehetővé, mert létrejött a világhálózat, az InterNet. Ezért napjaink egy-egy tudományterületének irodalomkutatása nemcsak a megjelenő szakkönyvek, folyóiratok és konferencia-kiadványok áttekintését jelenti, hanem az InterNeten található elektronikus folyóiratok, a kutatott tudományterülettel foglalkozó web-oldalak és a letölthető szakcikkek beható ismeretét feltételezi. A következőkben nem célunk a turbulencia irodalmának teljes körű áttekintése, hanem kizárólag azzal a szűk területtel foglalkozunk, amely az értekezés kérdésfelvetéseit érinti. A direkt numerikus szimuláció (DNS) összenyomhatatlan közeget feltételezve a kontinuitási egyenlet és a háromdimenziós Navier-Stokes-féle skalár mozgásegyenletek időfüggő megoldását szolgáltatja a turbulens áramlásban előforduló örvényméretek figyelembevételével. A direkt numerikus szimuláció során a Kolmogorov által bevezetett hossz-, idő- és sebességléptékeket alkalmazzák. A mozgás minden jellemző léptékének megoldása explicit módon történik. A k turbulens kinetikus energia Kolmogorov-féle k−5/3 törvénye [7] alapján a turbulens áramlás legkisebb örvényeinek mérete - amelyekben a kinetikus energia belső energiává történő disszipációja történik - a folyadék υ kinematikai viszkozitásának háromnegyedik hatványával arányos. A számítások során a Kolmogorov-féle hálózat csomópontjainak számát a kinematikai viszkozitás által definiált hossz- és időléptékek határozzák meg. A direkt numerikus szimuláció alkalmazásával kapcsolatos numerikus tapasztalatok [10], [12] azt mutatják, hogy a Reynoldsszám növelésével a hálózati elemek és a szükséges időlépések száma együttesen növekszik, ezáltal a számítógépes futtatás időigényessé válik. A módszer előnyei, hogy a) a mozgás pontos egyenleteinek numerikus megoldását szolgáltatja, ezért elvben helyesen tükrözi a turbulens áramlási problémákat; b) egyszerűbb áramlási feladatok esetén napjainkban kezdi átvenni a laboratóriumi mérések szerepét, mert a turbulens áramlások jellemzőivel kapcsolatban olyan információkat is szolgáltat, amelyek a mai mérési technikákkal hozzáférhetetlenek. A módszer hátránya, hogy korunk számítógépeinek teljesítménye mellett is nehézséget jelent a korlátozott memóriakapacitás, mert a számítások és a kapott eredmények feldolgozása nagy tárterületet igényelnek. A direkt numerikus szimuláció az egyetlen eddig született eljárás, amely mellőzi a klasszikus értelemben vett modellalkotást [15], [16], [17], [48]. A large eddy szimuláció (LES) mint számítási eljárás a finom struktúra modellek csoportjába sorolható. A számítási módszer egy szűrési eljárás segítségével a turbulens áramlásban fellépő kis és nagy örvénystrukturákat különválasztja, és a numerikus vizsgálatok középpontjában azok a nagyobb örvények állnak, amelyek a turbulens tömeg-, impulzus- és hő transzportjáért felelősek. A large eddy szimuláció - ugyanúgy mint a direkt numerikus szimuláció - a turbulens áramlási feladatok háromdimenziós és időfüggő megoldását szolgáltatja. A módszer előnyei, hogy a) ugyanazon turbulens áramlási probléma esetén kevesebb számítási műveletet igényel, mint a direkt numerikus szimuláció, ezért nagyobb Reynolds-szám tartományban is kielégítő számítási eredményeket szolgáltat, amelyek az elvégzett mérésekkel jó egyezést mutatnak; b) nagyobb időlépések is alkalmazhatók, mint amekkorák a direkt numerikus szimuláció esetén megengedhetők. A módszer hátrányai, hogy a) a sebességtér nagyobb méretű örvényeinek meghatározása során szűrési eljárást kell alkalmazni, amelynek következtében a skalár mozgásegyenletekben speciális Reynolds-féle látszólagos feszültségek jelennek meg, ezért a módszer alapvető problémája ezen feszültségek modellezése; b) mivel a fal közelében minden 5

esetben kis örvények találhatók, ezért az itt kapott számítási eredményeket döntően befolyásolja a Reynolds-féle feszültségek számítására alkalmazott turbulencia-modell [11], [12], [48]. Az egyre nagyobb teljesítményű számítógépek megjelenése kedvez a turbulens áramlási folyamatok direkt numerikus szimuláció és large eddy szimuláció útján történő megközelítésének, azonban ezek a módszerek a hétköznapi gyakorlat számára túlméretezettek [48]. A turbulencia-modellek következő csoportját a feszültség transzport modellek képezik. Osborne Reynolds nyomán a turbulens áramlást leíró alapegyenletek időátlagát a turbulens ingadozások periódusidejéhez mérten kellően nagy időintervallumban képezzük. Az időátlagolás következtében a mozgásegyenletekben az ismeretlen turbulens ingadozási komponensek szorzatainak időátlagai jelennek meg, amelyek a Reynolds-féle látszólagos feszültségtenzor elemeit alkotják. A fellépő új ismeretlenek száma minden esetben meghaladja a megoldandó differenciálegyenletek számát, ezért ezeket az ismeretlen turbulens ingadozási komponenseket modellezni kell, vagyis a megoldandó egyenletrendszert zárttá kell tenni. A modellezés során az ismeretlen mennyiségeket tapasztalati formulákkal helyettesítik vagy ezen ismeretlen turbulens ingadozásokra nézve újabb differenciálegyenleteket vezetnek be. A feszültség transzport modellek a Reynolds átlagolt Navier-Stokes (RANS) egyenletek megoldásán alapulnak, és a Reynolds-féle látszólagos feszültségek minél pontosabb meghatározására törekszenek. A Reynolds feszültség egyenlet modell a feszültség transzport modellek közé tartozik, mert a Reynolds-féle mozgásegyenletekben fellépő ismeretlen látszólagos turbulens feszültségeket transzportegyenletek bevezetésével modellezik. A Reynolds-féle feszültségtenzor szimmetriájából adódóan háromdimenziós turbulens áramlások esetén hat, kétdimenziós esetben három transzportegyenletet kell megoldani a kontinuitási és a mozgásegyenleteken túl. A megoldandó egyenletek mellett szükséges legalább egy járulékos egyenlet bevezetése is, amely a legtöbb esetben az ε turbulens kinetikus energiadisszipációjára vonatkozik, hogy az adott turbulens áramlási feladatban szereplő egyik lépték meghatározható legyen. A modellegyenletekben ugyan nem szerepel a k turbulens kinetikus energia, de az esetek döntő többségében mégis szükséges kiszámítani, hogy az adott feladat peremfeltételrendszere zárt legyen. Az additív egyenletek konstansait kísérleti úton vagy direkt numerikus szimuláció segítségével határozzák meg, illetve finomítják [11], [12], [48]. A Reynolds feszültség egyenlet modell előnye, hogy jól alkalmazható egyszerűbb és összetettebb turbulens áramlási feladatok esetén, valamint az átlagjellemzők kiszámítása mellett a Reynolds-féle feszültségek meghatározását is lehetővé teszi. A Reynolds feszültség egyenlet modell hátránya, hogy a megjelenő járulékos parciális differenciálegyenletek (transzportegyenletek) miatt rendkívül számításigényes. A turbulens áramlások ezúton történő modellezése J.C. Rotta (1951) [18] nyomán ismeretes. Az 1970-es években e kutatási területen K. Hanjalić, B.E. Launder, G.J. Reece, W. Rodi, V. Schumann, M. Gibson [19], [20], [21], [22], [23] értek el kiemelkedő eredményeket. Az 1980as évektől kezdődően a modell jelentősébb kutatói M. Leschziner (1989) [24], B. Launder (1990) [25], C. Meneveau (1991) [26], M. Ferge (1992) [27], A.V. Johansson, M. Hallbäck, S. Jakirlic, K. Hanjalic és W. Rodi (1994) [28], [29], [30]. A Reynolds feszültség egyenlet modell továbbfejlesztett változatai a nyomás-feszültség modellek (pressure-stress models) és a diffúziós modellek (diffusion models), amelyek napjaink kutatási tárgyát képezik [12]. Az örvényviszkozitási modellek a Boussinesq-féle örvényviszkozitási hipotézisen alapulnak és a Reynolds átlagolt Navier-Stokes (RANS) egyenlet modellek közé sorolhatók. J.O. Hinze (1959) [31], H. Tennekes, J. Lumley (1972) [32] és H. Schlichting (1982) [3] szerint a Boussinesq-féle örvényviszkozitási hipotézisen nyugvó turbulencia-modellek háromdimenziós áramlások esetén nem szolgáltatnak kielégítő eredményeket, mert a pontosabb számítások érdekében a Reynolds-féle feszültségek izotrop approximációja helyett a Reynolds-féle feszültségek tenzorikus, anizotrop leírása szükséges. A Reynolds-féle feszültségek algebrai modelljei (ARSM) a Boussinesq-féle örvényviszkozitási hipotézisen alapuló turbulencia-modellek hiányosságainak kiküszöbölésére törekszenek, és a Boussinesq-féle örvényviszkozitást nem al6

kalmazzák. A Reynolds-féle látszólagos turbulens feszültségtenzort nem-lineáris algebrai approximációkkal helyettesítik, további differenciálegyenletek bevezetése nélkül. A Reynolds-féle feszültségek algebrai modellezésében J. Lumley (1970) [33], W. Rodi (1976) [34], D.C. Wilcox, M.W. Rubesin (1980) [35], C.G. Speziale (1987) [36], és az 1990-es években T.-H. Shih, J. Zhu, J. L. Lumley (1995) [37] értek el kiemelkedő eredményeket. Az örvényviszkozitási modellek csoportján belül algebrai (null-egyenlet) modelleket, egy-egyenlet és két-egyenlet modelleket különböztethetünk meg a kontinuitási és a mozgásegyenleteken kívül bevezetett egyenletek száma szerint. A Boussinesq-féle örvényviszkozitási hipotézisen alapuló turbulencia-modellekben a turbulens áramlást leíró mozgásegyenletek hasonló alakot öltenek, mint a numerikusan jól kezelhető lamináris áramlás esetén [48]. Az algebrai (null-egyenlet) modellek a legegyszerűbb turbulencia-modellek, mert a domináns turbulens nyírófeszültséget skalár örvényviszkozitási tényezővel definiálják. Előnyük, hogy könnyen algoritmizálhatók és jó leírást adnak egyszerűbb turbulens áramlási feladatok esetén. Hátrányuk, hogy csak a folyadék átlagjellemzőit- és a domináns turbulens nyírófeszültségeket veszik figyelembe, és leváló áramlások esetén nehezen alkalmazhatók. Ilyen típusú modellek a L. Prandtl-féle (1925) [38] keveredési úthossz modell, a Baldwin-Lomax (1978) [39] és a Cebeci-Cousteix-féle (1999) [40] algebrai turbulencia-modellek. Az egy-egyenlet modellek a k turbulens kinetikus energiára vagy az ε turbulens kinetikus energiadisszipáció sebességére vonatkozóan vezetnek be egy-egy transzportegyenletet, és ezen jellemzők segítségével definiálják az örvényviszkozitást. Az egy-egyenlet modellek többségükben a k turbulens kinetikus energiára tartalmaznak transzportegyenletet. Az első ilyen turbulencia-modell L. Prandtl (1945) [41] nevéhez fűződik. P.R. Spalart és S.R. Allmaras (1992) [42] az örvényviszkozitásra vonatkozóan írtak fel új transzportegyenletet. A két-egyenlet modellek a k turbulens kinetikus energiára és az ε turbulens kinetikus energiadisszipáció sebességére vonatkozóan egyaránt tartalmaznak egy-egy transzportegyenletet, és ezek megoldása révén határozzák meg az örvényviszkozitást. Az első ilyen turbulenciamodell A. Kolmogorovtól (1942) [7] származik. A műszaki gyakorlatban a legelterjedtebb a B. Launder- és a D. Spalding-féle (1972) [9] k-ε turbulencia-modell. A standard k-ε modell, a kis Reynolds-szám modellek, a két-réteg modellek (two-layer models), és a k-ε renormalizációs módszerek a két-egyenlet modellek továbbfejlesztett változatai. Napjainkban a kutatók érdeklődésének középpontjában egyrészt a nem-lináris k-ε modell kidolgozása áll [12], [100]. Az egy-egyenlet és két-egyenlet modellek előnye, hogy egyszerűségükből adódóan a műszaki gyakorlatban jól alkalmazhatók, hatékonyan algoritmizálhatók. Hátrányuk, hogy alkalmazzák a Boussinesq-féle örvényviszkozitási hipotézist, és ebből következően olyan turbulens áramlások esetén, ahol a deformációsebesség tenzorában hirtelen változások lépnek fel, így görbült felületek közötti áramlások esetén, szekunder áramlások jelenlétében, a háromdimenziós áramlások döntő többségében és a leváló határrétegáramlások számítása esetén pontatlan megoldást eredményeznek [9], [11], [12], [48]. Az előzőek összefoglalásaképpen elmondható, hogy a XIX. század óta született turbulencia-modellek csupán egy-egy meghatározott áramkép esetén érvényesek. Czibere Tibor (1991-2001) [43], [44], [45] a múlt század első felében megalkotott klasszikus elméletekből kiindulva a Kármán-féle kétdimenziós hasonlósági hipotézist háromdimenziós esetre terjesztette ki, mert a turbulencia jelensége mindig háromdimenziós. Korábban a turbulens feszültségeket a viszkózus feszültségekhez hasonlóan a deformációsebességgel arányos mennyiségként próbálták meghatározni. A Czibere-féle új sztochasztikus turbulenciamodell szerint a Reynolds-féle látszólagos turbulens feszültségtenzor nem a deformációsebesség tenzorával van kapcsolatban, mert a turbulencia forrása a középsebességtér örvényességében keresendő. Czibere állandó sűrűségű folyadék turbulens áramlása esetén a Friedmannféle vektorvonal megmaradási tétel alapján új örvénytételeket vezetett le, miszerint a) az ör7

vényvonalak megmaradási sajátosságot sem a középsebességtérben, sem a turbulens ingadozás sebességterében nem mutatnak, mindkét sebességtérben az örvények szétszóródnak, a környezetbe diffundálnak; b) az örvénydiffúziót egyrészt a folyadékviszkozitás, másrészt a turbulencia okozza, következésképpen az örvénydiffúzió turbulens áramlásban eltűnő viszkozitás mellett is fellép [43], [44], [45]. A turbulens mozgás elméletében a fő problémát a látszólagos turbulens feszültségtenzor és a turbulens középsebesség kapcsolatának meghatározása jelentette. A turbulens ingadozás sebességterében felírt turbulens örvénydiffúzió parciális differenciálegyenlete alapján a turbulens ingadozás sebességtere és a középsebesség örvényessége között a kapcsolat egyértelműen kimutatható. Ebből az következik, hogy a sebességingadozás által meghatározott látszólagos turbulens feszültségtenzor a középsebességtér örvényességével és nem a deformációsebességgel van kapcsolatban. Kalmár László a "Numerikus módszer radiális szivattyúkban kialakuló súrlódásos áramlás meghatározására" (1997) [46] című Ph.D. értekezésében radiális átömlésű szivatytyúkban kialakuló súrlódásos áramlás jellemzőinek (sebesség- és nyomáseloszlás, valamint a súrlódási veszteségek eloszlásának) meghatározására alkalmas numerikus eljárást dolgozott ki. A Czibere-féle új sztochasztikus turbulencia-modell alapján a radiális átömlésű szivattyúk járókerekének lapátozott terében kialakuló turbulens sebességprofil, és a lapátcsatorna átáramlási keresztmetszetére vonatkozó átlagos nyomás és a súrlódási veszteség számítására alkalmas összefüggéseket vezetett le. Janiga Gábor a "Kétdimenziós turbulens nyíróáramlások számítása sík, valamint enyhén görbült falakkal határolt csatornákban" (2002) [48] című Ph.D. értekezésében elsőként dolgozott ki numerikus eljárást a turbulens áramlások Reynolds-féle mozgásegyenletének megoldására a Czibere-féle új sztochasztikus turbulencia-modell alapján. Az értekezés új numerikus módszert ismertet, amely egyenes alkotójú és görbült falak által határolt sík- és forgásszimmetrikus perdületmentes turbulens áramlások sebességterének számítására alkalmas. A turbulens mozgást leíró kontinuitási és skalár mozgásegyenletek diszkretizációja véges térfogatok módszerével történt és az áramlástani feladat lineáris algebrai egyenletrendszer megoldására vezett. Janiga a) a lamináris áramlási problémát leíró Navier-Stokes-féle mozgásegyenletek megoldása során alkalmazott SIMPLE nyomáskorrekciós algoritmust terjesztette ki az új turbulencia-modellen alapuló Reynolds-féle mozgásegyenlet megoldására; b) síklapok közötti turbulens áramlás esetén a mozgásegyenletnek egy új aszimptotikus megoldását is levezette; c) az aszimptotikus és numerikus megoldás során a sztochasztikus turbulencia-modell alkalmasan választott l léptékfüggvényét egy negyedfokú polinommal közelítette, és ennek felhasználásával a λ ellenállástényezőre új implicit formulát kapott. Az értekezésben bemutatott számítási eredmények a mérésekkel és a napjainkban sokszor a laboratóriumi méréseket helyettesítő direkt numerikus szimulációval összevetve nagyon jó egyezést mutattak. A műszaki gyakorlatban kétdimenziós lamináris áramlások sebességterének számítása esetén sokszor a Navier-Stokes-féle mozgásegyenletekből rotációképzéssel adódó örvénydiffúzió parciális differenciálegyenletét alkalmazzák és a kontinuitási egyenlet alapján bevezetik az áramfüggvényt. A lamináris áramlási feladatok ezúton történő megoldása során a megoldandó differenciálegyenletek rendszáma ugyan megnő, de az ismeretlen függvények száma eggyel csökken, ezért a numerikus eljárás zárt egyenletrendszerének felállítása kevesebb gondot jelent. A nyomás mennyisége a Navier-Stokes-féle mozgásegyenletekből a rotációképzés során eltűnik és az ismeretlen sebességkomponensek az áramfüggvényből differenciálással származtathatók. Az 1930-as években a lamináris örvénydiffúzió parciális differenciálegyenletének peremfeltételére vonatkozóan A. Thom (1933) [49] javasolt egy elsőrendű, a fal irányításától független összefüggést. A Thom-féle elsőrendű peremfeltétel alkalmazása elterjedt, mert a számítási tapasztalatok azt mutatják, hogy az alkalmazott numerikus eljárást stabilizálja. Az 1960-as évek második felében indult fejlődésnek a lamináris áramlások sebességterének az örvénydiffúzió parciális differenciálegyenletének segítségével történő numerikus megoldása. Ezt különböző 8

lamináris áramlási feladatok esetén K. Aziz, J. Hellums (1967) [50], G. Hirasaki, J. Hellums (1970) [51], A.J. Chorin (1973) [52], G. Mallinson, G. de Vahl Davis (1973, 1977) [53], [55], P. Roache (1976) [54], S. Richardson, A. Cornish (1977) [56], S. Dennis, D. Ingham, R. Cook (1979) [57], U. Ghia, K. Ghia, C. Shin (1982) [58], T. Gatski, C. Grosch, M. Rose (1982) [59], A. Wong, J. Reizes (1984) [60], S. Dennis (1985) [61], T. Gatski, C. Grosch (1985) [62], T. Gatski, C. Grosch, M. Rose (1989) [63], C.A.J. Fletcher (1991) [64], Mo-Hong Chou, Weicheng Huang (1996) [65], M.T. Nair, T.K. Sengupta (1996) [66] ismertté vált szerzők alkalmazták. Az értekezés a lamináris áramlásokra alkalmazott, az irodalomban ismert módszerek kiegészítéseként, a turbulens örvénydiffúzió új, Czibere-féle differenciálegyenletének megoldásán alapuló numerikus eljárást ismertet két egymást nem metsző koaxiális forgásfelület által határolt forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens áramlási feladatok számításán keresztül. Fontos megjegyeznünk, hogy az értekezés tárgyát képező differenciálegyenletet az angol nyelvű irodalomban legtöbbször örvénytranszport egyenletnek (vorticity transport equation) nevezik, amíg a német nyelvű irodalomban az örvénydiffúzió differenciálegyenlete (Wirbeldiffusion Differentialgleichung)1 az elterjedtebb elnevezés. A továbbiakban a német nyelvű irodalomban használt elnevezéshez igazodunk. A lamináris örvénydiffúzió parciális differenciálegyenletének diszkretizált alakja az irodalomban jól ismert váltakozó irányok módszerével megoldható. A váltakozó irányok módszerét2 elsőként D.W. Peaceman és H.H. Rachford (1955) [67] alkalmazta parabolikus- és elliptikus típusú parciális differenciálegyenletek megoldása esetén. Az ADI-módszer alapgondolata, hogy egy közelítő faktorizáció segítségével a többdimenziós feladatot egydimenziós feladatok szorzatára vezeti vissza. A módszer előnyei, hogy a) könnyen programozható, mert a kapott implicit együtthatómátrixok tridiagonális sávszerkezetűek; b) számítástechnikai szempontból memóriatakarékos. A módszer hátrányai, hogy a) nehezen alkalmazható olyan parciális differenciálegyenletek megoldása esetén, amelyek együtthatóiban nagy ugrások vagy erős gradiensek vannak; b) nincs olyan matematikai eljárás, amelynek segítségével előre megállapítható lenne, hogy a módszer hatékony lesz-e vagy sem. A váltakozó irányok módszerét J. Douglas és J.E. Gunn (1964) [68] általánosították. Az 1960-as évek óta a módszernek számtalan továbbfejlesztett változata született: a G.I. Marcsuk-féle (1976) [69] "szeleteléses módszer" (splitting method) vagy az N.N. Yanenko-féle (1979) [70] "frakcionális lépés módszere" (fractional step method). E módszer részletes ismertetése A.R. Mitchell (1969) [71], D.M. Young (1979) [72], A.R. Mitchell, D.F. Griffiths (1980) [73], C.A.J. Fletcher (1991) [74], C. Hirsch (1992) [75] és Stoyan Gisbert (1997) [76] neves szerők könyveiben megtalálható. A váltakozó irányok módszerét a középsebességtérben érvényes turbulens örvénydiffúzió új, Czibere-féle parciális differenciálegyenletének numerikus megoldása során körkeresztmetszetű csőben illetve nyugvó párhuzamos falak közötti résben kialakuló turbulens áramlások számítása esetén korábban már alkalmaztuk [P.9], [P.10], [P.11], [P.12], [KD.2]. Az értekezés tárgyát képező turbulens áramlási feladat numerikus megoldása során célszerű ortogonális görbevonalú koordinátahálózatot alkalmazni. Nyíri András (1990) [82] a vízgépek lapátozott terére vonatkozó háromdimenziós potenciálos áramlási feladat megoldása során egy ortogonális görbevonalú koordinátarendszer előállítására alkalmas eljárást dolgozott ki. A Nyíri-féle módszer kétdimenziós esetben különösen jól alkalmazható, amikor az áramlástechnikai gép meridiánmetszetét határoló peremgörbék analitikusan adottak. A műszaki gyakorlatban előforduló feladatok esetén a vizsgált áramlástechnikai gép meridiáncsatornájának peremgörbéi diszkrét pontjaik által könnyen megadhatók. A diszkrét pontokban vett függvényértékeket viszont véletlenszerű hibák terhelik, mert az áramlástechnikai gép lapátozott terének 1

N.J. Kotschin, I.A. Kibel, N.W. Rose: Theoretische Hydromechanik, Akademie-Verlag, Berlin, Band II., pp. 271-274. 1955. 2 Az angol nyelvű irodalomban Alternating Direction Implicit method, röviden "ADI-method". 9

geometriai adatait a szerkesztési rajzról pontatlanul lehet leolvasni. A nehézség kiküszöbölése érdekében célszerű a diszkrét ponthalmazon simítást végezni, a pontokban az első és második deriváltakat kiszámítani, és a másodrendű folytonosságot előírni. Az értekezés 5. fejezete a Nyíri-féle [81], [83] simító eljárást felhasználva egy új kétdimenziós másodrendben folytonos ortogonális görbevonalú koordinátahálózat generáló eljárást ismertet, amely diszkrét perempontokkal adott meridiáncsatornák cellákra bontására alkalmas. A hálógenerálás önálló szakterületet képez és számtalan eljárás ismeretes [78], [79], [80], [84], viszont az alkalmazott módszerek a másodrendű folytonosságtól eltekintenek. A Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszékén a Magdeburgi Egyetem együttműködésével 1990-től turbulens áramlásokkal kapcsolatos kísérleti kutatások is folynak. A mérések módszereivel és eredményeivel több publikáció [87-99] is foglalkozik, és e tématerületet Szabó Szilárd (2001) [47] habilitációs értekezésében tárgyalja. Az értekezésben vizsgált turbulens áramlási feladat numerikus módszerrel történő megoldása során kapott eredményeket a 4. fejezetben J. Nikuradse (1932) [85] és J. Laufer (1954) [86] méréseivel, valamint a Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke és a Magdeburgi Egyetem együttműködésével elvégzett mérésekkel [97] hasonlítjuk össze.

10

1.2. A FELADAT MEGFOGALMAZÁSA Az értekezés célja numerikus módszer kidolgozása két egymást nem metsző koaxiális forgásfelület által határolt forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens áramlás sebességterének számítására. A feladatot a középsebességtérben érvényes turbulens örvénydiffúzió új, Czibere-féle parciális differenciálegyenletének segítségével oldjuk meg. A kontinuitási egyenlet alapján bevezetjük az áramfüggvényt, és az örvényvektor zérustól különböző kerületi irányú komponensének egyenletéből az áramfüggvényre egy elliptikus típusú inhomogén parciális differenciálegyenletet nyerünk. Ezzel az áramlási feladat numerikus megoldására alkalmas parciális differenciálegyenlet-rendszer adódik. Az örvénydiffúzió egyenletében szereplő domináns turbulens nyírófeszültség a sztochasztikus turbulencia-modell algebrai egyenletével számítható. A feladat ezúton történő megoldását az indokolja, hogy a) a turbulens örvénydiffúzió parciális differenciálegyenletében a nyomás mint ismeretlen mennyiség nem szerepel; b) a szóban forgó differenciálegyenlet új differenciálegyenlet. Homogén, összenyomhatatlan folyadék teljesen kialakult turbulens áramlására korlátozzuk vizsgálatainkat. A turbulens áramlási feladat alapegyenleteit ortogonális görbevonalú koordinátarendszerben írjuk fel és a számítások során ívhosszkoordinátákat használunk. A numerikus megoldás során a cellákra bontott meridiáncsatorna ívelemeit és a koordinátavonalak görbületeit számszerűen ismernünk kell. Ezért kézenfekvő célkitűzés egy másodrendben folytonos kétdimenziós ortogonális görbevonalú koordinátahálózat generáló eljárás kidolgozása is. A kidolgozott numerikus eljárásokkal kapott eredményeket analitikus megoldásokkal és az irodalomból vett mérésekkel hasonlítjuk össze.

11

2. A NUMERIKUS ELJÁRÁS ELMÉLETI HÁTTERE A következőkben áttekintjük az összenyomhatatlan folyékony kontinuumok izotermikus turbulens áramlását leíró differenciálegyenleteket. A 3. fejezetben kidolgozott numerikus eljárás a középsebességtérben érvényes turbulens örvénydiffúzió új, Czibere-féle parciális differenciálegyenletének megoldásán alapszik, ezért ismertetjük a számítások során használt új háromdimenziós sztochasztikus turbulencia-modell elméleti alapjait. A vizsgált folyadékáramlás skalár differenciálegyenleteit ortogonális görbevonalú koordinátarendszerben írjuk fel. 2.1. A TURBULENS MOZGÁSOK ALAPEGYENLETEI A folyékony kontinuumok turbulens mozgása Euler-féle szemléletmódban két – helytől és időtől függő – sebességtér szuperpoziciójaként írható le. Az egyik a folyékony kontinuum ingadozások nélküli, de időben változó fő mozgását jelenti, a másik az időben gyorsan és sztochasztikusan változó sebességingadozás, amely a turbulencia jelenségét fejezi ki. Ezért a pillanatnyi mozgásállapot minden egyes jellemzőjét szokásos módon a középérték és az ingadozás összegeként képezzük. A középértéken egy adott mozgásjellemzőnek a turbulens mozgás periódusidejéhez mérten kellően nagy időintervallumban vett integrál-középértékét értjük. A vT (r,t) pillanatnyi sebességtérben az összenyomhatatlan folyékony kontinuum mozgása esetén érvényes a tömegmegmaradást kifejező ∇ ⋅ vT = 0

(2.1)

kontinuitási egyenlet. A Stokes-féle molekuláris viszkozitási törvény a viszkózus folyadék pillanatnyi turbulens mozgására is érvényes, ezért a pillanatnyi turbulens feszültségtenzor: FT = − pT I + σ T = − pT I + η (v T ! ∇ + ∇ ! v T ),

(2.2)

ahol pT a pillanatnyi nyomásérték, σ T = η (v T ! ∇ + ∇ ! v T ) a pillanatnyi turbulens súrlódási feszültségtenzor, amely az FT pillanatnyi feszültségtenzor deviátora, I az egységtenzor, és η a dinamikai viszkozitási tényező. A viszkózus folyadék turbulens áramlása esetén érvényes a mozgásjellemzők pillanatértékeivel értelmezett általános mozgásegyenlet:

∂ vT 1 + ( v T ⋅ ∇) v T = g + Div FT , ∂t ρ ahol a g vektor a térfogati erőket jelenti. A konzervatív erőtér esetén - amely potenciálos és stacionárius ( g = −∇U és ∂ U ∂ t = 0 ) - a viszkózus folyadék turbulens mozgása esetén érvényes a pillanatértékekkel felírt

∂ vT 1 + (v T ⋅ ∇ )v T = −grad U + Div (− pT I + σ T ) ∂t ρ

(2.3)

12

Navier-Stokes-féle mozgásegyenlet, ahol U az erőtér potenciálja. Az Ω T = rot v T jelölést bevezetve és a (2.3) mozgásegyenlet minden tagjának rotációját véve a pillanatértékekkel értelmezett turbulens sebességtérben érvényes

∂ ΩT + (v T ⋅ ∇ )Ω T − (Ω T ⋅ ∇ )v T = υ∆ Ω T ∂t

(2.4)

Helmholtz-Thomson-féle örvénytételt kapjuk, ahol υ = η ρ a folyadék kinematikai viszkozitási tényezője. A folyékony kontinuum izotermikus turbulens mozgását a pillanatértékekkel értelmezett vT (r,t) turbulens sebességtérben a (2.3)-(2.4) egyenletek írják le, amelyekből idő szerinti átlagolással a v(r,t) időbeni középértékekkel értelmezett sebességtérben érvényes differenciálegyenletek levezethetők. A v(r,t) időbeni középértékekkel értelmezett sebességtérben az összenyomhatatlan folyékony kontinuum izotermikus mozgása esetén érvényes a tömegmegmaradást kifejező ∇⋅v = 0

(2.5)

kontinuitási egyenlet. A v(r,t) sebességtérben turbulens mozgást végző folyadékelem Reynolds-féle mozgásegyenletét szokásos módon a (2.3) mozgásegyenlet időbeli átlagolásával nyerjük: ∂ v 1 + (v ⋅ ∇ )v = −grad U + Div − pI + σ − ρ v ′ ! v ′ , (2.6) ∂t ρ

[

(

)]

ahol a felülvonás időbeni középértéket jelent. A σ súrlódási feszültségtenzorra érvényes a Stokes-féle összefüggés: σ = η ( v ! ∇ + ∇ ! v) .

(2.7)

Az (2.6) egyenlet jobb oldalán a szögletes zárójelben álló utolsó tag a Reynolds-féle látszólagos turbulens feszültségtenzor: FR = − ρ ( v′ ! v′) , amely a viszkózus feszültségtenzorhoz hasonlóan felírható az FR = − p R I + σ R

(2.8)

alakban is, ahol p R a turbulens sebességingadozással együttjáró impulzuscsere okozta nyomástöbblet, amely az FR turbulens feszültségtenzor első skalárinvariánsának harmada: pR =

ρ v1′v1′ + v 2′ v 2′ + v 3′ v 3′ , 3

(

)

és σ R az FR feszültségtenzor deviátora:

(

 13 2v1′v1′ − v 2′ v 2′ − v 3′ v 3′  v ′2 v1′ σR = −ρ   v 3′ v1′ 

) 1 3

v1′v 2′ 2v ′2 v ′2 − v 3′ v 3′ − v1′v1′ v 3′ v ′2

(

) 1 3

v1′v 3′ v 2′ v 3′ 2v 3′ v 3′ − v1′v1′ − v ′2 v ′2

(

  .  

)

13

Az előzőekből következően a (2.6) középsebességtérben érvényes Reynolds-féle mozgásegyenlet a σ viszkózus feszültségtenzornak és a σ R turbulens feszültségtenzor deviátorának segítségével a következő alakban írható fel:

ρ

∂ v + ρ ( v ⋅ ∇) v = − ρ ∇Π + Div (σ + σ R ) , ∂t

(2.9)

ahol a Π = U + ( p + p R ) ρ egyenlettel értelmezett ún. teljes potenciál, amely az erőtér potenciáljának és a nyomáspotenciálnak az összege [45]. A turbulens középértékekkel értelmezett v(r,t) sebességtérben kialakuló Ω = rot v örvénytérre felírható a turbulens örvénydiffúzió differenciálegyenlete, amely a rot v × dr = 0 egyenlettel meghatározott örvényvonalak megmaradási sajátosságaira vonatkozik, és a Helmholtz-Thomson-féle (2.4) örvénytétel időbeni átlagolásával adódik:

(

)

∂Ω + (v ⋅ ∇ )Ω − (Ω ⋅ ∇ )v = υ∆ Ω + ∇ × v ′ × Ω ′ , ∂t

(2.10)

'

ahol Ω = rot v' a sebességingadozás örvényterét jelenti. A (2.10) örvénytétel szerint a v(r,t) turbulens középsebességtérnek a rot v × dr = 0 egyenlettel értelmezett örvényvonalai a υ → 0 határesetben sem maradnak meg, és az örvénydiffúzió mértékét a közeg viszkozitása, valamint döntő mértékben a turbulencia okozta sebességingadozás határozza meg [44], [45].

14

2.2. AZ ALKALMAZOTT TURBULENCIA-MODELL A számításaink során a Czibere-féle sztochasztikus turbulencia-modellt alkalmazzuk. A turbulencia-modell az áramlás középsebességteréhez egy q1′ , q 2′ , q 3′ ún. természetes koordinátarendszert rendel (2.1. ábra), amelynek bázisvektorait a v és rot v vektorok definiálják. A második koordináta-irányt a v × rot v, a harmadik koordináta-irányt a rot v vektor negatívja határozza meg. Az első koordináta-irány az előző kettő vektoriális szorzataként adódik, vagyis e1′ = e′2 × e′3 . Az alkalmazott turbulencia-modell szerint a természetes koordinátarendszerben az örvényerősség és a domináns turbulens nyírófeszültség között egyértelmű kapcsolat adható meg.

2.1. ábra Az ortogonális görbevonalú természetes koordinátarendszer a háromdimenziós határrétegben Az alkalmazott új sztochasztikus turbulencia-modell szerint a Reynolds-féle látszólagos turbulens feszültségtenzor az FR = ρκ 2 l 2 HΩ 2

(2.11)

összefüggés alapján értelmezhető, ahol ρ a folyadék sűrűsége, κ = 0,40704 a Kármánkonstans, l a turbulencia alkalmasan választott hosszúság dimenziójú léptékfüggvénye és H a turbulencia ún. hasonlósági tenzora, amelynek elemei konstansok: α 1 H =  1 β  µ ϑ

µ ϑ  , γ 

(2.12)

15

és az Ω örvényerősség az Ω = rot v vektor nagyságát határozza meg, amely a természetes koordinátarendszerben a következőképpen írható:

Ω=

1 ∂ (v1′ H 1′ ) , H 1′ H 2′ ∂ q 2′

(2.13)

amely az Ω örvényvekor q3′ koordináta-irányú Ω 3 ′ komponensének negatívja [44]: Ω = Ω 3′e′3 = −Ω e′3 . Az Ω előjele az 2.1. ábra szerinti természetes koordinátarendszer választása esetén pozitív. A természetes koordinátarendszer ilyen megválasztásával a q′2 koordináta-irány a faltól kifelé mutat és a sebesség nagysága a faltól távolodva növekszik. A következőkben jelölje az alkalmazott turbulencia-modell (2.11) látszólagos turbulens feszültségtenzor H hasonlósági tenzorát nem tartalmazó részét:

Θ (q1′ , q 2′ , q3′ , t ) = ρκ 2 l 2 Ω 2 ,

(2.14)

amely a turbulens feszültségtenzor első sorának illetve első oszlopának második eleme, és ez a nyíróáramlásban az ún. domináns turbulens nyírófeszültség. Az FR látszólagos turbulens feszültségtenzor és annak σ R deviátora a q1′ , q 2′ , q 3′ természetes koordinátarendszerben a FR = Θ (q1′ , q2′ , q3′ , t )H

σ R = Θ (q1′ , q ′2 , q 3′ , t ) H ∗

és

(2.15)

alakban írható, ahol H* a H hasonlósági tenzor deviátora: α ∗ ∗ H =  1  µ

1 µ   13 (2α − β − γ ) β ∗ ϑ  =  1 ϑ γ ∗   µ

1 3

1 (2 β − γ − α ) ϑ

µ  . ϑ  1  γ α β − − ( 2 ) 3

(2.16)

A számítások során a sebességvektorok által definiált természetes koordinátarendszer használata meglehetősen nehézkes, mert a koordinátarendszer irányai a sebességeknek megfelelően pontról pontra változnak, és egy adott pont esetén időben is változhatnak. Ezért a számítások során célszerű a sebességektől független q1, q2, q3 ortogonális görbevonalú koordinátarendszert alkalmazni. Ekkor az FR turbulens feszültségtenzort és annak σ R deviátorát az új koordinátarendszerben transzformált alakban kell felírnunk: FR = Θ (q1 , q 2 , q 3 , t ) G

és

σ R = Θ (q1 , q2 , q3 , t ) G ∗ ,

(2.17)

ahol a G illetve a G ∗ tenzor a természetes koordinátarendszerben értelmezett H hasonlósági tenzor illetve H ∗ deviátorának a q1 , q 2 , q 3 rendszerbe való transzformáltját jelenti: G = E ⋅ H ⋅ ET

és

G ∗ = E ⋅ H ∗ ⋅ ET . 16

Az E tenzor skalár elemei a v és az Ω = rot v vektorok skalár komponenseivel számíthatók, vagyis a bázisvektorok fizikai koordinátái a következők:  v1   − λ Ω1  ; Ω  1 − λ2  v Ω 2  1  v2 E21 = − λ ; Ω  1 − λ2  v Ω 3  1  v3 E31 = λ − ; Ω  1 − λ2  v E11 =

1

E12 =

v 2 Ω 3 − v3 Ω 2 ; v Ω 1− λ

E13 = −

Ω1 ; Ω

E 22 =

v3 Ω 1 − v1 Ω 3 ; v Ω 1− λ

E 23 = −

Ω2 ; Ω

E32 =

v1Ω 2 − v 2 Ω 1 ; v Ω 1− λ

E33 = −

Ω3 , Ω

1

2

1

2

1

2

ahol λ = (vΩ )/ ( v Ω ) . A q1′ , q 2′ , q 3′ természetes koordinátarendszer bázisvektorainak Eij ( i , j = 1,2,3 ) fizikai koordinátáit a két koordinátarendszer bázisvektorai egyértelműen meghatározzák:  e1 ⋅ e1′  E =  e 2 ⋅ e1′  e ⋅ e′  3 1

e1 ⋅ e′2 e 2 ⋅ e′2 e 3 ⋅ e′2

e1 ⋅ e′3   e 2 ⋅ e′3  e 3 ⋅ e′3 

és

 e1 ⋅ e1′  E =  e1 ⋅ e′2  e ⋅ e′  1 3 T

e 2 ⋅ e1′ e 2 ⋅ e ′2 e 2 ⋅ e ′3

e 3 ⋅ e1′   e 3 ⋅ e′2  . e 3 ⋅ e ′3 

A σ viszkózus feszültségtenzor és a σ R turbulens feszültségtenzor deviátora előbbi kifejezéseivel a (2.9) mozgásegyenlet az alkalmazott új háromdimenziós sztochasztikus turbulencia-modellnek megfelelően a következő alakban írható:

ρ

∂ v + ρ ( v ⋅ ∇) v = − ρ ∇Π + η∆ v + Div (Θ G ∗ ), ∂t

(2.18)

amelyből rotációképzéssel

∂Ω 1 + (v ⋅ ∇ )Ω − (Ω ⋅ ∇ ) v = υ ∆ Ω + ∇ × Div (Θ G ∗ ) ∂t ρ

(2.19)

a (2.10) turbulens örvénydiffúzió új, Czibere-féle differenciálegyenletének a középsebességtérben érvényes másik alakjához jutunk, amely a későbbi diszkretizáció tárgyát jelenti. A (2.18) és (2.19) differenciálegyenlet két különböző lehetőséget nyújt a turbulens áramlási feladatok megoldására. A numerikus megoldás során az utóbbi esetet választottuk, mert a (2.19) differenciálegyenletben a Π teljes potenciál mint ismeretlen mennyiség nem szerepel. Az alkalmazott turbulencia-modell szerint a (2.11) látszólagos turbulens feszültségtenzor nem a deformációsebesség tenzorával, hanem a középsebességtér örvényességével áll összefüggésben. Czibere Tibor összenyomhatatlan folyadékra a Friedmann-féle vektorvonal megmaradási tétel [44] alapján kimutatta, hogy a) az örvényvonalak megmaradási sajátosságot sem a középsebességtérben, sem a turbulens ingadozás sebességterében nem mutatnak, mindkét sebességtérben az örvények szétszóródnak, a környezetbe diffundálnak; b) a (2.19) egyenlet alapján az örvénydiffúziót egyrészt a viszkozitás, másrészt a turbulencia okozza, következésképpen az örvénydiffúzió turbulens áramlásban eltűnő viszkozitás mellett is fellép [44]. 17

3. A NUMERIKUS ELJÁRÁS 3.1. A FELADAT PARCIÁLIS DIFFERENCIÁLEGYENLET-RENDSZERE A következőkben két egymást nem metsző koaxiális forgásfelület által határolt forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens áramlás sebességterének numerikus módszerrel történő meghatározásával foglalkozunk. A vizsgált áramlás parciális differenciálegyenlet-rendszerét először koordinátás alakban a választott q1, q2, q3 ortogonális görbevonalú koordinátarendszerben írjuk fel, és ezután a numerikus eljárással történő számítások megkönnyítésének érdekében a meridiánsíkon bevezetjük az ívhosszkoordinátákat. Az ívhosszkoordinátákkal felírt alapegyenletek diszkretizációját a véges differenciák módszerének segítségével végezzük. A levezetéseket a Függelék tartalmazza. Az időben állandósult turbulens áramlás a forgástengelyre illeszkedő meridiánsík két helykoordinátájának függvényeként írható le. Az áramlást koaxiális forgásfelületek alkotta falak határolják, ezért a számítások során alkalmazott q1, q2, q3 ortogonális görbevonalú koordinátarendszert úgy választjuk meg, hogy a q1, q2 koordinátavonalak a forgásszimmetrikus tér meridiánsíkjában legyenek és az áramlást határoló szilárd falfelület a q2 = konst. koordinátafelület legyen. A 3.1. ábrán az áramlási feladat alkalmasan választott ortogonális görbevonalú számítási koordinátarendszere látható, amelyen az áramlást határoló szilárd fal meridiángörbéje a q1 irányú koordinátavonallal egybeesik és a q3 irányú koordinátavonalak körök.

3.1. ábra A számítási forgásszimmetrikus ortogonális görbevonalú koordinátarendszer forgásfelület által határolt áramlási térben A forgásszimmetrikus térben kialakuló perdületmentes stacionárius áramlás sebességtere a meridiánsík q1, q2 koordinátáival a következőképpen írható le: v = v1 (q1 , q2 )e1 + v2 (q1 , q2 )e 2 , és a hozzátartozó örvényvektor kerületi irányú komponense:

Ω=

1  ∂ (v 2 H 2 ) ∂ (v1 H 1 )  − . H 1 H 2  ∂ q1 ∂ q2 

(3.1) 18

Az összenyomhatatlan (ρ = állandó) folyadék áramlása esetén érvényes  ∂ (v1 H 2 H 3 ) ∂ (v 2 H 1 H 3 )  1 + =0 H 1 H 2 H 3  ∂ q1 ∂ q2 

(3.2)

kontinuitási egyenlet alapján a meridiánsíkon bevezetjük az áramfüggvényt, amellyel a q1 és q 2 irányú sebességkomponensek: v1 =

1 ∂Ψ ; H 3 H 2 ∂ q2

v2 = −

1 ∂Ψ ; H 3 H 1 ∂ q1

alakban írhatók, és ezeket behelyettesítve a (3.1) összefüggésbe, a Ψ (q1 , q 2 ) áramfüggvényre vonatkozó másodrendű parciális differenciálegyenlet adódik (Függelék, F.I.):  1 ∂ H2 ∂ 2Ψ ∂ 2Ψ 1 ∂ H 3  ∂Ψ  + +  − + 2 2 (H 1∂ q1 ) (H 2 ∂ q 2 )  H 2 H 1 ∂ q1 H 3 H 1 ∂ q1  H 1∂ q1  1 ∂ H1 1 ∂ H 3  ∂Ψ  +  − = −Ω H 3 . H H q H H q H q ∂ ∂ ∂ 2 3 2 2  2 2  1 2

(3.3)

A koaxiális forgásfelületek által határolt forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens ármalás sebességterét a középsebességtérben érvényes turbulens örvénydiffúzó (2.19) differenciálegyenletének segítségével számítjuk, amelynek skalár differenciálegyenletei forgásszimmetrikus áramlás esetén a következő alakban írhatók: ∂ 1 H 3 H 2 ∂ q2

{[Div (Θ G )] H }= 0 , 3

(3.4)

{[ (

)] H }= 0 ,

(3.5)

∗

3

∂ 1 Div Θ G ∗ H 3 H 1 ∂ q1

3

3

1  ∂ (v1Ω H 2 ) ∂ (v 2 Ω H 1 ) + = H 1 H 2  ∂ q1 ∂ q2  υ  ∂  H 2 ∂ (Ω H 3 ) ∂ =   + H 1 H 2  ∂ q1  H 3 H 1 ∂ q1  ∂ q 2 +

(

{[ (

∂ 1 Div Θ G ∗ ρ H 1 H 2 ∂ q1

)]

2

}

H2 −

 H 1 ∂ (Ω H 3 )    +  H 3 H 2 ∂ q2 

∂ ∂ q2

{[Div (Θ G )] H } , ∗

1

1

(3.6)

)

ahol Div Θ G ∗ vektor skalár komponensei:

19

[Div (Θ G )] = ∂H(Θ∂Gq ) + ∂H(Θ∂Gq ) + ∗ 11

∗

1

1

G − G +Θ   H 2 H1 ∗ 11

∗ 22

∗ 12

1

2

2

 2 ∂ H1 ∂ H2 G − G ∂ H3 1 ∂ H 3   , + + G12∗  + H 3 H 1 ∂ q1 ∂ q1  H 1 H 2 ∂ q 2 H 3 H 2 ∂ q 2  ∗ 11

∗ 33

[Div (Θ G )] = ∂H(Θ∂Gq ) + ∂H(Θ∂Gq ) + ∗ 12

∗

∗ 22

2

1

1

2

2

G − G ∂ H 2 G − G ∂ H3  2 ∂ H1 1 ∂ H 3    . +Θ  + + G12∗  + H H ∂ q H H ∂ q H H ∂ q H H ∂ q  1 2 2 1 1 3 1 1 2 3 2 2   ∗ 22

∗ 22

∗ 11

∗ 33

A G* deviátor elemeinek kiszámításához a q1′ , q 2′ , q 3′ természetes és a q1 , q 2 , q 3 számítási koordinátarendszerek között transzformációt kell végeznünk. Kétdimenziós feladat feltételezésével az E transzformációs tenzor elemei a következőképpen egyszerűsödnek: E11 = E21 =

v1 v12 + v22 v2

;

E12 =

;

E22 = −

v12 + v22

E 31 = 0 ;

v2 v12 + v22 v1 v12 + v22

E 32 = 0 ;

Ω ; Ω Ω ; Ω

E13 = 0 ; E 23 = 0 ; E33 = −

Ω . Ω

A G ∗ = E ⋅ H ∗ ⋅ E T deviátor elemei a v1 és v 2 sebességkomponensekkel kifejezhetők: G11∗ = α ∗

v12 v 22 vv Ω + β +2 21 2 2 , ∗ 2 2 2 2 v1 + v 2 v1 + v 2 v1 + v 2 Ω

v 22 v12 vv Ω G = α∗ 2 + β∗ 2 −2 21 2 2 , 2 2 v1 + v 2 v1 + v 2 v1 + v 2 Ω ∗ 22

∗ =γ ∗, G33

v1v 2 v 22 Ω  v12  G = G = (α ∗ − β ∗ ) 2 − − v1 + v 22 Ω  v12 + v 22 v12 + v 22 v1 v2 Ω ∗ G13∗ = G31 = −µ −ϑ , 2 2 Ω 2 v1 + v2 v1 + v22 ∗ 12

∗ 21

∗ G23 = G32∗ = − µ

v2 v12 + v 22

  , 

v1 Ω +ϑ , Ω v12 + v 22

ahol α * , β * , γ * , µ és ϑ az alkalmazott sztochasztikus turbulencia-modell alkalmasan választott konstansai. A számítások során célszerű a µ illetve ϑ állandókat zérus értékűnek megválasztani, mert ebből következik, hogy a Div Θ G ∗ vektor harmadik skalár komponense zérus lesz, tehát [Div Θ G ∗ ]3 = 0 esetén a (3.4)-(3.5) transzportegyenletek eltűnnek.

(

)

(

)

20

A q1′ , q 2′ , q 3′ természetes koordinátarendszerben a (2.14) algebrai egyenlettel értelmezett Θ domináns turbulens nyírófeszültséget szintén transzformálni kell a q1 , q 2 , q 3 számítási koordinátarendszerbe. A domináns turbulens nyírófeszültség skalár mennyiség, amely a meridiánsík q1 , q 2 koordinátáival kifejezve:

Θ (q1 , q 2 ) = ρ κ 2 l 2 Ω 2 .

(3.7)

A turbulencia l léptékfüggvényét a meridiánmetszet minden egyes q2 irányú trajektóriája mentén meg kell határoznunk, amelynek értéke a falakon zérus, mert azokon a turbulencia jelensége nem lép fel. Ha az áramlás irányára merőleges ξ ívhosszkoordinátát a meridiánmetszet középső áramfelületétől, vagyis a beírható körök középpontjától számítjuk, akkor a számítási tapasztalatok alapján célszerű a léptékfüggvény negyedfokú polinommal történő közelítése: l (ξ ) =

4S h

 4S − 1  ξ  2    (ξ − ξ A )(ξ B − ξ ) , 1 − S  h   

(3.8)

ahol a h = ξ B − ξ A a meridiáncsatorna ívhosszban mérhető szélességét jelenti, ezért ξ B a q 2 irányú koordinátavonalon vett trajektóriának és a külső peremnek a metszéspontja, ugyanígy ξ A a belső peremmel való metszéspontot jelöli (3.2. ábra). A ξ ívhosszkoordinátát a meridiáncsatorna középső áramfelületétől mérjük, és ennek következtében a csatorna két falához a ξ A = − h / 2 és ξ B = h / 2 ívhosszkoordináták tartoznak.

3.2. ábra A turbulencia léptékfüggvényének meghatározása egy trajektória mentén A Θ (q1 , q 2 ) domináns turbulens nyírófeszültség függvény zérushelye nem feltétlenül esik a középső áramfelületre, viszont jó közelítéssel annak környezetében található. Az S paraméter a számítási tapasztalataink szerint 0,25 < S < 2 értékhatárok között változhat. Czibere Tibor szerint [45] az alsó határérték körkeresztmetszetű csőáramlás esetén, a felső határérték forgó hengerek közötti Couette-áramlás esetén a kísérletekkel jól egyező eredményekre vezetett. 21

Koordinátafelületek Koordinátafelületek koordinátavonalai q2 , q3 q1 = konst.

A koordinátafelületek göbületei

1 1 ∂ H2 1 1 = = ; R12 H2 H1 ∂ q1 R13 H 3 H1 q3 , q1 1 1 ∂ H1 1 1 = = ; R21 H1 H2 ∂ q2 R23 H 3 H2 q1 ,q2 1 1 ∂ H1 1 1 = = ; R31 H1 H3 ∂ q3 R32 H2 H3 3.1. táblázat A koordinátafelületek görbületeinek összefüggései

q2 = konst. q3 = konst.

∂ H3 ∂ q1 ∂ H3 ∂ q2 ∂ H2 ∂ q3

A numerikus eljárás alkalmazásának illetve a számítások megkönnyítésének érdekében a meridiánsíkon bevezetjük a H1dq1 = ds1 és a H 2 dq 2 = ds2 ívhosszkoordinátákat, és a H 3 dq 3 = rdϕ ívhosszkoordinátát a meridiánsíkra merőlegesen. A 3.3. ábrán a forgásszimmetrikus számítási koordinátarendszer koordinátafelületeinek görbületi viszonyai láthatók, mert az áramlást határoló forgásfelület egyben koordinátafelület is.

G2

PG1 = R12 PG2 = R21 PO1 = R13 PO2 = R23

q1 P

P q3

r O1

q2

α

O2

G1 3.3. ábra A forgásszimmetrikus tér ortogonális görbevonalú koordinátafelületeinek görbületi viszonyai Az Ω (s1 , s 2 ) függvény ívhosszkoordinátákkal felírva (Függelék, F.II.):

Ω (s1 , s2 ) =

(

∂ v2 ∂ v1 v2 v − + − 1 , ∂ s1 ∂ s 2 R12 R21

(3.9)

)

és a Div Θ G ∗ vektor skalár komponensei: 22

[Div (Θ G )] = ∂ (Θ∂ sG ) + ∂ (Θ∂ sG ) + Θ  G R− G ∗ 11

∗

∗ 12

∗ 11

1

[Div (Θ G )] ∗

(

1

)

(

2

)



∗ 22

+

12

 2 G11∗ − G33∗ 1   , + G12∗  + R13  R21 R23 

∗ ∗ ∗  G 22  2 ∂ Θ G12∗ ∂ Θ G 22 − G11∗ G 22 − G33∗ 1    . + +Θ  + + G12∗  + 2= R23 ∂ s1 ∂ s2  R12 R13    R21

Ezeket felhasználva a forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens áramlás sebességterének numerikus eljárással történő meghatározására alkalmas diszkretizálandó differenciálegyenlet-rendszer a (3.3) és (3.6) egyenletek alapján (Függelék, F.III.):

1  ∂Ψ  1 1  ∂Ψ ∂ 2Ψ ∂ 2Ψ  1   + +  − +  − = −Ω r , 2 2 ∂ s1 ∂ s 2  R12 R13  ∂ s1  R21 R23  ∂ s 2

(3.10)

 ∂ 2 Ω ∂ 2 Ω  1  1  ∂Ω  Ω − υ  2 +  +  + + 2 ∂ s 2  R12 R13  ∂ s1  ∂ s1   1 1  ∂Ω Ω  ∂ 2 r ∂ 2 r  1 1  ∂r  1 1  ∂ r       + + +  2 + 2 +  − +  −  = R R s r R R s R R ∂ ∂ s s ∂ ∂  21  12  21 23  2 2 13  1 23  ∂ s 2    1  ∗ ∗ ∗  G 22  2 ∂ Θ G22 − G11∗ G 22 − G33∗ 1 1  ∂ Θ G12∗ 1       + = + + Θ + + G12∗  +   ρ R12  ∂ s1 ∂ s2 R R R R 21 23 12 13     

v1

v ∂Ω ∂Ω  ∂ v1 ∂ v 2 v1 + v2 +  + + + 2 ∂ s1 ∂ s 2  ∂ s1 ∂ s 2 R12 R21

(

+

(

)

(

)

)

∗ ∗   G 22 − G11∗ G 22 − G33∗ 1    ∗  2     − Θ + + G +   12  R23   R21  R12 R13    ∗  G11∗ − G22  2 G11∗ − G33∗ 1       − +Θ  + + G12∗  + R R R R  21 12 13 23     

∗ ∂ 2 Θ G12∗ ∂ 2 Θ G 22 ∂ + + 2 ∂ s1 ∂ s 2 ∂ s1 ∂ s1

1 − R21 −

) (

(

(

)

(

 ∂ Θ G11∗ ∂ Θ G12∗ +  ∂ s2  ∂ s1

)

(

)

)

∗   G11∗ − G 22  2 G11∗ − G33∗ 1 + + G12∗  + Θ   R13   R12  R21 R23

∂ 2 Θ G12∗ ∂ 2 Θ G11∗ ∂ − − 2 ∂ s 2 ∂ s1 ∂ s2 ∂ s2

    ,  

(3.11)

ahol a v1 és v 2 sebességkomponensek az áramfüggvény segítségével felírva: v1 =

1 ∂Ψ ; r ∂ s2

v2= −

1 ∂Ψ . r ∂ s1

A G ∗ = E ⋅ H ∗ ⋅ E T deviátor elemei az alkalmazott turbulencia-modell µ illetve ϑ konstansainak zérus értékű választása mellett a v1 és v 2 sebességkomponensekből számíthatók: G11∗ = α ∗

v12 v 22 vv Ω β + +2 21 2 2 , ∗ 2 2 2 2 v1 + v 2 v1 + v 2 v1 + v 2 Ω

v 22 v12 vv Ω G = α∗ 2 + β∗ 2 −2 21 2 2 , 2 2 v1 + v 2 v1 + v 2 v1 + v 2 Ω ∗ 22

∗ G33 = γ∗,

23

v1v 2 v 22 Ω  v12  G = G = (α ∗ − β ∗ ) 2 − − v1 + v 22 Ω  v12 + v 22 v12 + v 22 ∗ G13∗ = G31∗ = G 23 = G32∗ = 0 . ∗ 12

∗ 21

  ,  (3.12)

A (3.10) differenciálegyenlet megoldásával a Ψ áramfüggvényből a sebességkomponensek differenciálással származtathatók. Az Ω függvény a (3.9) összefüggés alapján és a turbulens örvénydiffúzió (3.11) differenciálegyenletében fellépő G* deviátor elemei a v1 és v2 sebességkomponensekkel meghatározhatók. A Θ domináns turbulens nyírófeszültség a turbulenciamodell (3.7) algebrai egyenletével számítható. Következésképpen a (3.10) és (3.11) differenciálegyenletek a Ψ és Ω ismeretlen függvények meghatározásának szempontjából egy zárt parciális differenciálegyenlet-rendszert alkotnak. A turbulens örvénydiffúzió (3.11) differenciálegyenletéből az is látható, hogy a turbulencia jelensége eltűnő viszkozitás mellett is fellép. Ez akkor jelentős, ha az áramlás teljes egészében turbulenssé válik, és ekkor a molekuláris viszkozitás hatása a turbulens impulzuscseréhez képest elhanyagolható. 3.2. A DISZKRETIZÁCIÓ ÉS A NUMERIKUS MEGOLDÁS A (3.10) és (3.11) differenciálegyenletek diszkretizációját ortogonális görbevonalú koordinátahálózaton végezzük és a tömörebb írásmód érdekében az (i,j) indexrácson bevezetjük a 3.4. ábra szerinti jelöléseket. A koordinátahálózat egy vizsgált P csomópontja körül elhelyezkedő pontokat az angol nyelvű irodalomból átvett - N, S, W, E (észak, dél, nyugat, kelet) égtájaknak megfelelő jelölésekkel látjuk el [101]. A számításokhoz szükséges geometriai adatrendszert az 5. fejezetben kidolgozott kétdimenziós másodrendben folytonos hálógeneráló eljárás segítségével állítjuk elő, amely szolgáltatja a számítási pontok koordinátái mellett a koordinátavonalak görbületeit is. A későbbiek során felhasználjuk, hogy a (3.10) és (3.11) differenciálegyenletekben megjelenő görbületi sugarak numerikusan ismertek.

3.4. ábra A P csomópont körül elhelyezkedő pontok jelölése az (i,j) indexrácsban

24

A Ψ (s1 , s 2 ) eloszlást a (3.10) differenciálegyenlet diszkretizációjából nyerjük, amelynek első lépése, hogy a bevezetett jelöléseket alkalmazva a differenciálegyenletben szereplő első és másodrendű parciális deriválatakat véges differencia formulákkal közelítjük [101]:

Ψ E (s P − sW ) −Ψ P (s E − sW ) + Ψ W (s E − s P ) Ψ N (s P − s S ) −Ψ P (s N − s S ) + Ψ S (s N − s P ) + + 1 1 (s P − sW ) ⋅ (s E − sW ) ⋅ (s E − s P ) (s P − s S ) ⋅ (s N − s S ) ⋅ (s N − s P ) 2 2  1 1  Ψ E −ΨW  1 1  Ψ N −Ψ S   +  − +  − = −Ω P rP . R R s − s R R E W  12  21 13  P 23  P s N − s S

(3.13)

A numerikus megoldás során a Ψ (s1 , s 2 ) áramfüggvény az ortogonális görbevonalú koordinátahálózat egy P belső és annak négy szomszédos pontjában ismeretlen, ezért meg kell határoznunk az alábbi algebrai egyenlet együtthatóit és jobboldali forrástagját: ASΨ S + AWΨ W + APΨ P + AEΨ E + ANΨ N = QP ,

(3.14)

ahol a Ψ S ,Ψ W ,Ψ P ,Ψ E ,Ψ N diszkrét pontokban vett értékeknek az AS , AW , AP , AE , AN együtthatóit és a Q P jobboldali forrástagot a (3.13) rendezése után kapjuk: 2 1  1 1    , − − (s P − s S )⋅ (s N − s S ) s N − s S  R21 R23  P 2 1  1 1    , AW = − − (s P − sW )⋅ (s E − sW ) s E − sW  R12 R13  P

AS =

  1 1 AP = −2 + ,  (s P − sW )⋅ (s E − s P ) (s P − s S )⋅ (s N − s P ) 2 1  1 1    , AE = + − (s E − sW )⋅ (s E − s P ) s E − sW  R12 R13  P AN =

2 1  1 1    , + − (s N − s S )⋅ (s N − s P ) s N − s S  R21 R23  P Q P = −Ω P rP .

Az áramlási feladat peremfeltételei mérések alapján meghatározhatók, amelyek zárttá teszik a (3.10) differenciálegyenlet közelítő megoldását. A vizsgált áramlás nyugvó koaxiális forgásfelületek által határolt térben megy végbe, ezért a folyadék sebessége a falon megegyezik a fal sebességével, vagyis a v1 és v2 sebességkomponensek eltűnnek. Ebből következik, hogy a falon Ψ = állandó és a normális irányú deriváltja eltűnik (∂Ψ/∂n = 0). A belső peremen az áramfüggvény értéke zérus és a külső peremen az átáramló q térfogatárammal egyezik meg. A belépő és kilépő oldalakon a meridiánsebesség-eloszlás mérések alapján meghatározható, vagy a belépő oldalon a meridiánsebesség-eloszlást a mérések alapján előírjuk és a kilépő oldalon a sebességkomponensek normális irányú deriváltjait zérusnak tekintjük, ami azt jelenti, hogy a kilépésnél az áramfüggvényre vonatkozóan másodrendű Neumann-típusú peremfeltételt is megadhatunk. Az előzőekből következik, hogy a számítási hálózat belső pontjaiban érvényes (3.14) algebrai egyenlet a tartomány peremén a peremfeltételeknek megfelelően módosul. 25

A (3.14) algebrai egyenletet a számítási hálózat minden egyes csomópontjában felírjuk, tekintettel a tartomány peremén érvényes peremfeltételekre, és így a (3.10) paricális differenciálegyenlet numerikus megoldását lineáris algebrai egyenletrendszer megoldására vezetjük vissza. A koordinátahálózat diszkrét pontjaiban a Ψ (s1 , s 2 ) áramfüggvény közelítő értékeit a 3.5. ábra szerinti sávos együtthatómátrixú lineáris egyenletrendszer megoldása során nyerjük.

3.5. ábra A Ψ (s1 , s 2 ) eloszlás esetén a megoldandó ötsávos egyenletrendszer általános alakja A Ψ (s1 , s 2 ) eloszlás ismeretében a v1 és v 2 sebességkomponensek a tartomány egy adott belső P pontjában véges differencia formulák alkalmazásával számíthatók:

(v1 )P

=

1 Ψ N −Ψ S , rP (s N − s S )

és

(v 2 )P

=−

1 Ψ E −Ψ W , rP (s E − sW )

(3.15)

valamint a pontosabb közelítés érdekében célszerű a peremeken másod- vagy annál magasabb rendű hibataggal rendelkező haladó illetve retrográd differenciákkal számolni. A v1 és v 2 sebességkomponensek ismeretében és a (3.9) összfüggés felhasználásával a tartomány belsejében az Ω (s1 , s 2 ) eloszlás a következő

Ωp =

(v2 )E − (v2 )W (v1 )N − (v1 )S s E − sW

−

sN − sS

v v  +  2 − 1   R12 R21  P

(3.16) 26

centrális differencia-formula alapján számítható. Az Ω (s1 , s 2 ) függvény peremfeltételeit a Ψ (s1 , s 2 ) eloszlásból származtathatjuk. A kilépésnél az áramlástani feladatnak megfelelően Neumann-típusú peremfeltételt is előírhatunk. A G* deviátor elemeit a számítási hálózat diszkrét pontjaiban a v1 és v 2 sebességkomponensek és az Ω (s1 , s 2 ) eloszlás ismeretében az alábbiak szerint határozzuk meg:

(G )

∗ 11 P

(G )

∗ 22 P

 v2   v2   vv  = α ∗  2 1 2  + β ∗  2 2 2  + 2 2 1 2 2   v1 + v 2  P  v1 + v 2  P  v1 + v 2  P  v2   v2   vv  = α ∗  2 2 2  + β ∗  2 1 2  − 2 2 1 2 2   v1 + v 2  P  v1 + v 2  P  v1 + v 2  P = γ ∗,

(G ) (G ) = (G ) ∗ 33 P

∗ 12 P

∗ 21 P

 vv  Ω = (α ∗ − β ∗ ) 2 1 2 2  − P ΩP  v1 + v 2  P ∗ ∗ = G 23 P = G 32 P = 0 ,

(G ) = (G ) ( ) ( ) ∗ 13 P

∗ 31 P

ΩP , ΩP ΩP , ΩP

 v12 v 22  2 − 2 2 2  v1 + v 2 v1 + v 2

  , P

ahol az alkalmazott turbulencia-modell H hasonlósági tenzorának α = −3,76 , β = −1,65 , γ = −2,705 főátló menti elemeihez a H* deviátor α * = −1,055 , β * = 1,055 és γ * = 0 főátló menti elemei tartoznak. A trajektóriák mentén a turbulencia léptékfüggvényének (3.8) ismeretében S = 1 paraméterérték választása mellett a

Θ P = ρ κ 2 l P2 Ω P2

(3.17)

domináns turbulens nyírófeszültség diszkrét pontbeli értékeit a turbulencia-modell algebrai összefüggése alapján számítjuk. Az l léptékfüggvény (3.8) szerinti megválasztásával a Θ (s1 , s 2 ) domináns turbulens nyírófeszültségre Dirichlet peremfeltételek adódnak. A továbbiakban a turbulens örvénydiffúzió (3.11) differenciálegyenletét diszkretizáljuk, és újabb iterációt végzünk az Ω (s1 , s 2 ) eloszlás meghatározására. A diszkretizáció első lépésében a tömörebb írásmód miatt bevezetjük az alábbi jelöléseket: a=

∂ v1 ∂ v 2 v v + + 1 + 2 , ∂ s1 ∂ s 2 R12 R21 1 1 b= + , R12 R13 1 1 c= + , R21 R23

1  ∂ 2r ∂ 2r  1 1  ∂r  1 1  ∂r    d =  2 + 2 +  − +  − , r  ∂ s1 ∂ s 2  R12 R13  ∂ s1  R21 R23  ∂ s 2  valamint a (3.11) differenciálegyenlet jobboldali turbulens tagját jelölje:

27

(

) (

)

(

)

(

)

∗ ∗  1  ∂ Θ G12∗  ∂ 2 Θ G12∗ ∂ Θ G 22 ∂ 2 Θ G 22 ∂e e + + − + + +    2 ∂ s2 ∂ s1 ∂ s 2 ∂ s1 ∂ s1  R12  ∂ s1   ∂ 2 Θ G12∗ ∂ Θ G12∗ ∂ 2 Θ G11∗ ∂ f  1  ∂ Θ G11∗ f − − − + + − ,   2 R21  ∂ s1 ∂ s2 ∂ s 2 ∂ s1 ∂ s 2  ∂ s2 

T=

1 ρ

(

)

(

)

(

)

(

)

ahol ∗  G ∗ − G11∗ G 22  2 − G33∗ 1   , e = Θ  22 + + G12∗  + R23  R21  R12 R13  ∗ ∗  G11∗ − G22  2 G11∗ − G33 1   , + + G12∗  + f =Θ  R R R R  21 12 13 23  

és ezek alapján a diszkretizálandó (3.11) differenciálegyenlet a v1

 ∂ 2Ω ∂ 2Ω  ∂Ω ∂Ω ∂Ω ∂Ω + v2 + aΩ − υ  2 + +b +c + dΩ  = T 2 ∂ s1 ∂ s2 ∂ s1 ∂ s2 ∂ s2  ∂ s1 

(3.18)

tömör alakban írható. Ez az egyenlet mind lamináris, mind turbulens áramlás számítása esetén érvényes. Ha az áramlás lamináris, akkor l = 0 és vele együtt Θ = 0, tehát T = 0, vagyis az előző egyenlet jobboldala eltűnik. Ha az áramlás teljesen egészében turbulenssé válik, akkor a molekuláris viszkozitás hatása a turbulens impulzuscseréhez képest elhanyagolható, ezért az előző egyenlet eltűnő viszkozitás mellett is alkalmas turbulens áramlás számítására. A (3.18) átrendezése után az ismeretlen Ω (s1 , s 2 ) mennyiségre vonatkozóan a −υ

∂ 2Ω ∂ 2Ω ∂Ω ∂Ω − υ + (v1 − υ b ) + (v 2 − υ c ) + (a − υ d )Ω = T 2 2 ∂ s1 ∂ s2 ∂ s1 ∂ s2

(3.19)

differenciálegyenlet adódik, ahol az első és másodrendű parciális differenciálhányadosokat a szokásos módon véges differencia formulákkal közelítjük és az alábbi

Ω E (s P − sW ) − Ω P (s E − sW ) + Ω W (s E − s P ) − 1 ( ) ( ) ( ) s − s ⋅ s − s ⋅ s − s P W E W E P 2 Ω N (s P − s S ) − Ω P (s N − s S ) + Ω S (s N − s P ) −υ + 1 ( )( )( ) 2 sP − sS ⋅ s N − s S ⋅ s N − s P Ω − ΩW Ω − ΩS + (v1 − υ b )P E + (v 2 − υ c )P N + (a − υ d )P Ω P = TP s E − sW s N − sS −υ

(3.20)

differenciaegyenletet nyerjük. Az Ω (s1 , s 2 ) eloszlás a számítási hálózat egy adott P belső és annak négy szomszédos pontjában ismeretlen, ezért meg kell határoznunk a következő algebrai egyenlet együtthatóit és a jobboldali turbulens forrástagját: AS* Ω S + AW* Ω W + AP* Ω P + AE* Ω E + AN* Ω N = TP ,

(3.21) 28

ahol az Ω S , Ω W , Ω P , Ω E , Ω N diszkrét pontokban vett értékeknek az AS* , AW* , AP* , AE* , AN* együtthatóit a (3.20) differenciaegyenlet rendezése után nyerjük:

(v −υ c )P 2υ − 2 , (s P − s S )⋅ (s N − s S ) s N − s S (v −υ b )P 2υ AW* = − − 1 , (s P − sW )⋅ (s E − sW ) s E − sW

AS* = −

  1 1 AP* = 2υ  +  + (a − υ d )P ,  (s P − sW )⋅ (s E − s P ) (s P − s S )⋅ (s N − s P ) (v −υ b )P 2υ AE* = − + 1 , (s E − sW )⋅ (s E − s P ) s E − sW (v −υ c )P 2υ AN* = − + 2 , (s N − s S )⋅ (s N − s P ) s N − s S ahol  1 1   , +  R12 R13  P

 1 1   , +  R21 R23  P

(v1 − υ b )P = (v1 ) p − υ 

(v1 )E − (v1 )W

(v 2 −υ c )P = (v2 )P − υ 

(v2 )N − (v2 )S

 v v  +  1 + 2  − s E − sW s N − sS  R12 R21  P υ  r (s − sW ) − rP (s E − sW ) + rW (s E − s P ) rN (s P − s S ) − rP (s N − s S ) + rS (s N − s P ) −  E 1P + + 1 rP  s − s ⋅ s − s ⋅ s − s s − s ⋅ s − s ⋅ s − s ( ) ( ) ( ) ( ) ( ) ( ) P W E W E P P S N S N P 2 2

(a −υ d )P

=

+

 1 1  rE − rW  1 1  rN − rS    +  − +  − , R R s − s R R s − s W S   12  21 13  P E 23  P N  és a (3.20) differenciaegyenlet jobboldali TP turbulens forrástagja: TP =

1 ρ

 1   R12

(Θ G ) (s + ∗ 12 E

P

[(Θ G ) +

∗ 22 NE

(

  P

(

) (

 Θ G12∗ E − Θ G12∗ ⋅ s E − sW 

(

)

W

+

)

(Θ G ) − (Θ G ) ∗ 22 N

s N − sS

(

− sW ) − Θ G12∗ P (s E − sW ) + Θ G12∗ 1 ( )( )( ) 2 s P − sW ⋅ s E − sW ⋅ s E − s P

(

∗ − Θ G 22

) ]⋅ (s SE

(s NW

) (

∗ 22 S

[(

) (s W

)

E

(

− sP )

∗ ∗ − s SW ) + Θ G22 NW − Θ G 22 − s SW )⋅ (s E − sW )⋅ (s NE − s SE ) NW

)

(

) (

 1   Θ G11∗ E − Θ G11∗ W Θ G12∗ N − Θ G12∗  ⋅  −  + R s − s s N − sS E W  21  P 

)

S

 + eP  + 

+

) ]⋅ (s SW

NE

− s SE ) e E − eW + − s E − sW

 + fP  −  29

(Θ G ) (s − ∗ 12 N

[(Θ G ) −

∗ 11 NE

P

(

)

(

− s S ) − Θ G12∗ P (s N − s S ) + Θ G12∗ 1 ( )( )( ) 2 s P − sS ⋅ s N − sS ⋅ s N − s P

(

− Θ G11∗

) ]⋅ (s SE

(s NW

[(

)

) (s S

N

(

− sP )

− s SW ) + Θ G11∗ NW − Θ G11∗ − s SW )⋅ (s E − sW )⋅ (s NE − s SE ) NW

−

) ]⋅ (s SW

NE

− s SE )

−

fN − fS   , (3.21) s N − sS 

ahol a számítási hálózat P belső pontjában ∗ ∗  G ∗ − G11∗ G 22  2 − G33 1 e P = Θ P  22 + + G12∗  + R23  R12 R13  R21

  ,  P

∗  G11∗ − G 22  2 G11∗ − G33∗ 1   . fP =ΘP  + + G12∗  + R13 R R  R12  21 23  P

A (3.21) algebrai egyenletet a számítási hálózat minden egyes csomópontjában felírjuk, tekintettel a tartomány peremén érvényes peremfeltételekre, és így a turbulens örvénydiffúzió (3.11) paricális differenciálegyenletének numerikus megoldását szintén lineáris algebrai egyenletrendszer megoldására vezetjük vissza. Az Ω (s1 , s 2 ) függvény közelítő értékeit a 3.6. ábra szerinti sávos együtthatómátrixú lineáris egyenletrendszer megoldása során nyerjük.

Az Ω (s1 , s 2 )

3.6. ábra eloszlás esetén a megoldandó ötsávos egyenletrendszer általános alakja 30

A diszkretizáció során kapott ötsávos ritka együtthatómátrixú lineáris egyenletrendszereket a Nyíri-féle sávos együtthatómátrixú lineáris egyenletrendszerek megoldására kidolgozott új módszer [83] segítségével oldottuk meg. A Nyíri-féle módszer alkalmazását az indokolta, hogy az eljárás kevesebb műveletet igényel, mint a sávos Gauss-elimináció, és pontosság szempontjából a vizsgált hidrodinamikai feladat esetén hatékonynak bizonyult. Az iteráció Ω r = 0 választásával kezdődik, ezért először a koaxiális forgásfelületek által határolt térben kialakuló potenciálos áramlást számítjuk. Ebben az esetben a (3.10) differenciálegyenlet jobb oldala eltűnik, így a számítási koordinátarendszerben érvényes forgásszimmetrikus Laplace-egyenletet kell megoldanunk. A véges differenciák módszerén alapuló numerikus eljárás algoritmusa röviden a következő lépésekben foglalható össze: • A (3.10) differenciálegyenlet numerikus megoldásával egy-egy iterációs lépésen belül a Ψ (s1 , s 2 ) eloszlást határozzuk meg. • A Ψ (s1 , s 2 ) eloszlás ismeretében a v1 és v 2 sebességkomponensek számíthatók, amelyek* * ből az Ω (s1 , s 2 ) eloszlás és a G* deviátor G11* , G12* , G22 , G33 elemei meghatározhatók. • Az Ω (s1 , s 2 ) eloszlás és a turbulencia alkalmasan választott l léptékfüggvényének ismeretében egy-egy iterációs lépésen belül - a (3.7) algebrai egyenlet alapján - a Θ (s1 , s 2 ) domináns turbulens nyírófeszültség eloszlás számítható. Az l léptékfüggvény (3.8) szerinti megválasztásával a Θ (s1 , s 2 ) domináns turbulens nyírófeszültség peremfeltételei adódnak. • Az Ω (s1 , s 2 ) eloszlás meghatározására újabb iterációt végzünk, ezért a turbulens örvénydiffúzió (3.11) differenciálegyenletének diszkretizációja során előálló lineáris algebrai egyenletrendszert megoldjuk. A kapott Ω (s1 , s 2 ) eloszlást r-el megszorozzuk, amely a következő iterációs lépésben a (3.10) differenciálegyenlet diszkretizációjából származó lineáris algebrai egyenletrendszer jobboldali oszlopvektorát szolgáltatja. • A Ψ (s1 , s 2 ) és az Ω (s1 , s 2 ) függvények rekurzív számítását addig ismételjük, amíg az ismeretlen mennyiségek két iterációs lépésben vett különbségének sorösszeg normája egy előírt ε küszöbértéknél kisebb érték lesz:  n  max  ∑ Ψ i ,(mj ) −Ψ i ,(mj−1)  < ε , i  j =1 

és

 n  max  ∑ Ω i(,mj ) − Ω i(,mj−1)  < ε , i  j =1 

ahol n az ismeretlenek száma és az m felsőindex az iterációs lépés számát jelöli. A numerikus tapasztalataink azt mutatják, hogy ε ≤ 10 −6 vagy ennél kisebb nagyságrend előírása javasolt. A kidolgozott numerikus eljárás a Czibere-féle új sztochasztikus turbulencia-modell alkalmazásán alapszik. A domináns turbulens nyírófeszültségek figyelembevételével a (3.10)(3.11) parciális differenciálegyenlet-rendszer numerikus megoldása során a véges differenciák módszerét alkalmaztuk. Ezt az indokolta, hogy a középsebességtérben érvényes turbulens örvénydiffúzió (3.11) Czibere-féle parciális differenciálegyenlete új differenciálegyenlet, ezért kézenfekvő volt a legegyszerűbb numerikus eljárás alkalmazása. A kidolgozott számítási módszer rendkívül egyszerű és hatékony, mert kevés számítógépi memóriaigénye mellett igen gyors, hiszen egy-két perc futási idő alatt képes koaxiális forgásfelületek által határolt térben a turbulens áramlás sebességterének számítására. A számítások során nagyon pontosan kell ismernünk a számítási koordinátahalózat csomópontjaiban a görbületi sugarakat, ezért nélkülözhetetlen az 5. fejezetben ismertetett kétdimenziós másodrendben folytonos hálógeneráló eljárás alkalmazása. Ügyelnünk kell a hálósűrűség és a turbulencia léptékfüggvényének megválasztására is. A kidolgozott számítási eljárás további előnye, hogy turbulens örvénydiffúzió (3.11) differenciálegyenletében a nyomás mennyisége nem szerepel az ismeretlen fizikai mennyiségek 31

között, ezért a számítások során nem kell nyomáskorrekciós algoritmust alkalmazni, mint a szóban forgó turbulencia-modell Reynolds-féle mozgásegyenletének numerikus megoldása során, ami további számítási időt venne igénybe. A numerikus eljárás nemcsak görbe alkotójú forgásfelületek közötti térben kialakuló perdületmentes stacionárius turbulens áramlások számítására alkalmas, hanem egyenes alkotójú forgásszimmetrikus térben kialakuló turbulens csőáramlás esetén is használható. Ha viszont a (3.3) és (3.6) ortogonális görbevonalú koordinátarendszerben felírt parciális differenciálegyenleteket a H3 = 1 választása mellett írjuk fel, akkor a kétdimenziós áramlás esetén érvényes differenciálegyenlet-rendszert kapjuk. Az áramfüggvényre vonatkozó (3.3) differenciálegyenlet forgásszimmetrikus és H3 = 1 mellett kétdimenziós potenciáláramlások számítására is alkalmas. A számításainkat saját készítésű C nyelvű számítógépes program segítségével végeztük. A kapott eredményeket viszkózus folyadék turbulens áramlása esetén analitikus megoldásokkal és mérésekkel hasonlítottuk össze, amelyeket a 4. fejezetben ismertetünk.

32

4. A NUMERIKUS ELJÁRÁSSAL SZÁMÍTOTT EREDMÉNYEK 4.1. TURBULENS ÁRAMLÁS KÖRKERESZTMETSZETŰ CSŐBEN A numerikus eljárást először hosszú egyenes körkeresztmetszetű csőben teljesen kialakult turbulens áramlás esetén teszteltük. A kapott eredményeket a Czibere-féle analitikus megoldással [45], valamint J. Nikuradse [85] és J. Laufer [86] méréseivel hasonlítottuk össze. A numerikus megoldás tárgyát képező (3.11) turbulens örvénydiffúzió parciális differenciálegyenlete szerint a turbulencia jelensége eltűnő viszkozitás mellett is fellép, ezért a számításokat ebben az esetben is elvégeztük. A 4.1.1.-4.1.8. ábrák a v~ átlagsebesség és a D csőátmérő által definiált Re = v~D / υ Reynolds-számok mellett, a vm sebességmaximummal dimenziótlanított sebességprofilokat ábrázolnak. A körkeresztmetszetű csőben kialakuló turbulens sebességprofil szimmetrikus, ezért a szimmetriatengely és a fal közötti sebességeloszlást vázoltuk. 4.2. TURBULENS ÁRAMLÁS KOAXIÁLIS FORGÁSFELÜLETEK ÁLTAL HATÁROLT FORGÁSSZIMMETRIKUS TÉRBEN A numerikus eljárással kapott eredményeket a FLUENT v6.1 programrendszerrel számított eredményekkel hasonlítottuk össze - a standard k-ε turbulencia-modell alkalmazásával - azonos koordinátahálózatokon és azonos peremfeltételek mellett. A számítások során a Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke és a Magdeburgi Egyetem együttműködésével elvégzett LDV (Laser Doppler Velocimetry) lézeres méréstechnika által szolgáltatott mérési eredményeket [97] is figyelembe vettük. A belépésnél a 4.2.5. ábrán látható sebességeloszlást feltételeztük Q = 68,8 [m3/h] térfogatáram mellett. A numerikus eljárással történő számításokat a 4.2.21.-4.2.24. ábrák szerinti különböző sűrűségű másodrendben folytonos koordinátahálózatokon végeztük. A 4.2.6.-4.2.20. ábrák a forgásszimmetrikus áramlási tér egy-egy jellemző keresztmetszetében axiális irányú sebességprofilokat ábrázolnak. A 4.2.1.4.2.2. ábrákon a meridiánsík pontjaihoz rendelt színek a tengelyirányú sebességkomponensek nagyságát ábrázolják, így a kialakuló áramlási stukturák jól láthatók. 4.3. A SÍKÁRAMLÁS SPECIÁLIS ESETE A értekezés tárgyát a síkbeli turbulens nyíróáramlások számítása ugyan nem képezi, viszont a 3. fejezetben kidolgozott numerikus eljárás síkáramlások sebességterének meghatározására is alkalmas. Ha a (3.3) és (3.6) ortogonális görbevonalú koordinátarendszerben felírt alapegyenleteket a H3 = 1 választása mellett írjuk fel, akkor a síkáramlás esetén érvényes, diszkretizálandó differenciálegyenlet-rendszert kapjuk. A párhuzamos síkfalak között kialakuló turbulens áramlás számítását [P.12] a v~ átlagsebesség és a DH hidraulikai átmérő által definiált Re = v~DH / υ Reynolds-számok esetén a Janiga-féle aszimptotikus megoldással [48] hasonlítottuk össze. A kialakuló sebességprofil szimmetrikus, ezért a szimmetriatengely és a fal közötti sebességeloszlást vázoltuk, ahol b a síkcsatorna szélességének felét jelöli.

33

v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (viszkózus áramlás)

0,4

Mérési eredmények J. Nikuradse

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.1. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Nikuradse [85] mérésével (Re = 4000) v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (eltűnő viszkozitás mellett)

0,4

Mérési eredmények J. Nikuradse

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.2. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Nikuradse [85] mérésével (Re = 4000)

34

v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (viszkózus áramlás)

0,4

Mérési eredmények J. Nikuradse

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.3. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Nikuradse [85] mérésével (Re = 23300) v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (eltűnő viszkozitás mellett)

0,4

Mérési eredmények J. Nikuradse

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.4. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Nikuradse [85] mérésével (Re = 23300)

35

v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (viszkózus áramlás)

0,4

Mérési eredmények J. Laufer

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.5. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Laufer [86] mérésével (Re = 40000) v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (eltűnő viszkozitás mellett)

0,4

Mérési eredmények J. Laufer

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.6. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Laufer [86] mérésével (Re = 40000)

36

v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (viszkózus áramlás)

0,4

Mérési eredmények J. Nikuradse

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.7. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Nikuradse [85] mérésével (Re = 105000) v vm

1 0,8 Analitikus megoldás Czibere T.

0,6

Numerikus megoldás (eltűnő viszkozitás mellett)

0,4

Mérési eredmények J. Nikuradse

0,2 0 0

0,2

0,4

0,6

0,8

1

r R 4.1.8. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Czibere-féle analitikus megoldással [45] és J. Nikuradse [85] mérésével (Re = 105000)

37

4.2.1. ábra A FLUENT v6.1 programrendszerrel számított turbulens sebességeloszlás tengelyirányú komponensei Q = 68,8 [m3/h] belépő térfogatáram mellett

4.2.2. ábra A numerikus eljárással számított turbulens sebességeloszlás tengelyirányú komponensei Q = 68,8 [m3/h] belépő térfogatáram mellett

38

4.2.3. ábra A mérési pontok helyei [97]

4.2.4. ábra A tengelyirányú sebességkomponens mért eloszlása [97]

39

3

2,4

1,8

1,2

0,6

0 0

0,01

0,02

0,03

0,04

0,05

4.2.5. ábra Az összehasonlító számítások során felvett belépő sebességeloszlás a Czibere-féle analitikus megoldás alapján Q = 68,8 [m3/h] térfogatáram mellett

2,8

2,24

1,68 Standard k-epszilon modell 1,12

Numerikus megoldás

0,56

0 0,003

0,0124

0,0218

0,0312

0,0406

0,05

4.2.6. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,019422 [m] helyen kezdődő trajektória mentén

40

2,8

2,24

1,68 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

1,12

0,56

0 0,007

0,0156

0,0242

0,0328

0,0414

0,05

4.2.7. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,027925 [m] helyen kezdődő trajektória mentén

2,4

1,92

1,44 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,96

0,48

0 0,016

0,0232

0,0304

0,0376

0,0448

0,052

4.2.8. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,038319 [m] helyen kezdődő trajektória mentén

41

2

1,6

1,2 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,8

0,4

0 0,03

0,035

0,04

0,045

0,05

0,055

4.2.9. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,047831 [m] helyen kezdődő trajektória mentén

1,8

1,44

1,08 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,72

0,36

0 0,034

0,0386

0,0432

0,0478

0,0524

0,057

4.2.10. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,050381 [m] helyen kezdődő trajektória mentén

42

3,2

2,56

1,92

1,28 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,64

0 0,072

0,0756

0,0792

0,0828

0,0864

0,09

4.2.11. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,072825 [m] helyen kezdődő trajektória mentén

3

2,4

1,8

1,2

Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,6

0 0,072

0,0756

0,0792

0,0828

0,0864

0,09

4.2.12. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,0898 [m] helyen kezdődő trajektória mentén

43

2,4

1,92

1,44 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,96

0,48

0 0,069

0,0728

0,0766

0,0804

0,0842

0,088

4.2.13. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,096753 [m] helyen kezdődő trajektória mentén

2

1,6

1,2 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,8

0,4

0 0,065

0,0686

0,0722

0,0758

0,0794

0,083

4.2.14. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,102407 [m] helyen kezdődő trajektória mentén

44

1,6

1,28

0,96 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,64

0,32

0 0,055

0,059

0,063

0,067

0,071

0,075

4.2.15. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,109727 [m] helyen kezdődő trajektória mentén

1,9

1,52

1,14

0,76 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,38

0 0,04

0,0452

0,0504

0,0556

0,0608

0,066

4.2.16. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,120531 [m] helyen kezdődő trajektória mentén

45

2,4

1,92

1,44

0,96 Standard k-epszilon modell 0,48

0 0,027

Numerikus megoldás

0,034

0,041

0,048

0,055

0,062

4.2.17. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,135175 [m] helyen kezdődő trajektória mentén

2,4

1,92

1,44

0,96 Standard k-epszilon modell Numerikus megoldás Mérési adatok [97]

0,48

0 0,022

0,03

0,038

0,046

0,054

0,062

4.2.18. ábra Sebességeloszlás a 4.2.23. ábra szerinti z = 0,145654 [m] helyen kezdődő trajektória mentén

46

3,2

2,56

1,92 1617 csomópont 4100 csomópont 5740 csomópont 5940 csomópont

1,28

0,64

0 0,072

0,0756

0,0792

0,0828

0,0864

0,09

4.2.19. ábra A hálósűrítés hatása a számított sebességeloszlások esetén a 4.2.23. ábra szerinti z = 0,072825 [m] helyen kezdődő trajektória mentén

3

2,4

1,8 1617 csomópont 4100 csomópont 5740 csomópont 5940 csomópont

1,2

0,6

0 0,072

0,0756

0,0792

0,0828

0,0864

0,09

4.2.20. ábra A hálósűrítés hatása a számított sebességeloszlások esetén a 4.2.23. ábra szerinti z = 0,0898 [m] helyen kezdődő trajektória mentén

47

4.2.21. ábra A számítások során használt 1617 csomópontból álló koordinátahálózat

4.2.22. ábra A számítások során használt 4100 csomópontból álló koordinátahálózat 48

4.2.23. ábra A számítások során használt 5740 csomópontból álló koordinátahálózat

4.2.24. ábra A számítások során használt 5940 csomópontból álló koordinátahálózat 49

1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (viszkózus áramlás)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.1. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 23300) 1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (eltűnő viszkozitás mellett)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.2. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 23300)

50

1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (viszkózus áramlás)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.3. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 27000) 1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (eltűnő viszkozitás mellett)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.4. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 27000)

51

1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (viszkózus áramlás)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.5. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 43800) 1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (eltűnő viszkozitás mellett)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.6. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 43800)

52

1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (viszkózus áramlás)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.7. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 105000) 1 0,8 0,6

Analitikus megoldás Janiga G.

0,4

Numerikus megoldás (eltűnő viszkozitás mellett)

0,2 0 0

0,2

0,4

0,6

0,8

1

4.3.8. ábra A numerikus eljárással számított dimenziótlan sebességprofil összehasonlítása a Janiga által kidolgozott aszimptotikus megoldással síkfalak közötti áramlás esetén [48] (Re = 105000)

53

4.4. A SZÁMÍTÁSI EREDMÉNYEK ÉRTÉKELÉSE A körkeresztmetszetű csőben kialakuló turbulens áramlás esetén a numerikus eljárással számított eredmények a Czibere-féle analitikus megoldással [45], valamint J. Nikuradse [85] és J. Laufer [86] méréseivel jó egyezést mutatnak. A számított eredmények és a mért értékek 10 −2 ÷ 10 −4 nagyságrendben megegyeznek [KD.2]. A vizsgált meridiánmetszet mind egy szivattyú, mind egy turbina meridiáncsatornájának modellje lehet. A Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke és a Magdeburgi Egyetem együttműködésével elvégzett mérések során az LDV (Laser Doppler Velocimetry) lézeres méréstechnika a koaxiális forgásfelületek által határolt mérőtérben csak korlátozottan volt alkalmazható és speciális korrekciók voltak szükségesek [47], [97]. A mérési nehézségekből adódik, hogy a belépő meridiánsebesség-eloszlás pontosan nem ismert, ezért volt szükséges a numerikus eljárással történő számítások során a belépő sebességprofilt felvenni a mérésnél ismert belépő térfogatáram mellett. A kapott eredmények szerint a mért és számított sebességterek hasonló áramlási struktúrát mutatnak. A 4.2.12. ábrán a számítási eredmények mellett azokat a mérési adatokat tűntettünk fel, amelyek esetén a legkevesebb korrekciót kellett elvégezni. Az elsődleges célunk az volt, hogy a Czibere-féle sztochasztikus turbulencia-modellen alapuló számítási eljárást a FLUENT v6.1 programrendszer által ismert standard k-ε modell számítási eredményeivel hasonlítsuk össze. Az eredmények a műszaki gyakorlatban a jó számszerű egyezések miatt elfogadhatók. A síkfalak között kialakuló turbulens áramlás számítása [P.12] szintén a kidolgozott numerikus eljárás hatékonyságát igazolja. A kapott eredményeket a Janiga Gábor által kidolgozott síkáramlásra vonatkozó aszimptotikus megoldással hasonlítottuk össze, amelyek korábban direkt numerikus szimuláció eredményeivel voltak összehasonlítva [48]. A számítási eredmények és az aszimptotikus megoldás összehasonlítása ugyanúgy, mint a körkeresztmetszetű csőben végzett számítások esetén is 10 −2 ÷ 10 −4 nagyságrendben egyezést mutatnak.

54

5. KÉTDIMENZIÓS MÁSODRENDBEN FOLYTONOS ORTOGONÁLIS GÖRBEVONALÚ KOORDINÁTAHÁLÓZAT LÉTREHOZÁSA A következőkben egy kétdimenziós másodrendben folytonos ortogonális görbevonalú koordinátahálózat generáló eljárást ismertetünk. A meridiánmetszetet határoló két peremgörbét az (r,z) síkon diszkrét pontjaikkal adjuk meg és a Nyíri-féle simító eljárást [81] alkalmazzuk, mert a diszkrét pontokban megadott függvény értékeit véletlenszerű hibák terhelik. A pontokbeli első és második deriváltakat a negyedfokú Newton-féle interpolációs polinom segítségével számítjuk és a kapott pontokra a másodrendű csatlakozási feltételt kielégítő ötödfokú Hermite-polinomot illesztünk, előírva ezzel a másodrendű folytonosságot. A tengelyirányú koordinátavonalak diszkrét pontjait a beírható körök középpontjainak számításával határozzuk meg. A kapott pontokat simítjuk, a pontokbeli első és második deriváltakat a negyedfokú Newton-féle interpolációs polinom alkalmazásával számítjuk, és a pontokra ismét Hermitepolinomot illesztünk. Ezt az eljárást addig folytatjuk, amíg kellő mennyiségű tengelyirányú koordinátavonalat kapunk. Az iránymezőből az ortogonális trajektóriák számíthatók. 5.1. A NUMERIKUS ELJÁRÁS ALGORITMUSA 1. A meridiánmetszetet határoló két peremgörbét diszkrét pontjaikkal adjuk meg. 2. A Nyíri-féle simító eljárás alkalmazása a diszkrét pontokban. 3. A diszkrét pontokban az első és második deriváltak számítása a negyedfokú Newton-féle interpolációs polinom alkalmazásával. 4. A hat feltételt kielégítő másodrendben csatlakozó ötödfokú Hermite-polinom illesztése a diszkrét pontokra. 5. A beírható körök középpontjainak számítása. 6. A 2.-3.-4. pontokban leírtak alkalmazása a beírható körök középpontjaira. 7. Az 5.-6. pontokban foglaltak ismétlése, amíg kellő mennyiségű koordinátavonalat kapunk. 8. Az iránymezőből az ortogonális trajektóriák meghatározása a középvonaltól kezdődően. 9. Az ívelemek és koordinátavonalak görbületeinek számítása.

5.2. A NYÍRI-FÉLE SIMÍTÓ ELJÁRÁS ALKALMAZÁSA Az n darab nem egyenközű z 0 , z1 ,  , z n−1 (zi < z i +1 ) pontokban adottak az η0 , η1 ,  , η n−1 simítandó függvényértékek, és keressük azokat az r0 , r1 ,, rn −1 simított függvényértékeket, amelyek minimalizálják (n > 2m ) esetén a következő S összeget: S=

n − (m −1)

∑ (δ k =0

m

rk

) + ∑ ε (r 2

n −1

k

k =0

k

− ηk ) , 2

(5.1)

55

∂S = 0 szükséges feltétel alapján ker∂ rk

és ezt az összeget minimalizáló rk értékrendszert a essük, ahol

δ m rk =

δ m−1rk +1 − δ m−1rk z k + m − zk

(5.2)

a k-adik m-ed rendű osztott differencia és megállapodás szerint δ 0 rk ≡ rk , valamint ε k > 0 a simítás mértékét meghatározó simítási paraméter. A minimumot eredményező rk (k = 0,  , n − 1) értékeknek ki kell elégíteniük az alábbi lineáris algebrai egyenletrendszert:

∑ a( ) r k

λ

k +λ

λ

= ε kη k ,

(5.3)

ahol a megoldandó sávos együtthatómátrixú lineáris egyenletrendszer együtthatóit az −1

ak(λ )

m  m  = ∑ ∏ (z k − z k −i + j )⋅ ∏ (z k +λ − z(k +λ )−i + j ) + δ k ,λ ε k   i −imin j = 0 j =0 j ≠i  j ≠i  imax

(5.4)

összefüggés alapján számítjuk, amely a λ indexű sáv k-adik elemét jelenti, és a Kroneckerszimbólum a következőképpen változhat:

δ k ,λ = 0 ,

ha

λ ≠ 0;

δ k ,λ = 1 ,

ha

λ =0;

valamint ak(λ ) = 0 , ha (k + λ < 0 ) vagy (k + λ > n − 1). A sávos együtthatómátrixú lineáris egyenletrendszert a Nyíri-féle direkt módszerrel oldjuk meg [83], és úgy építjük fel, hogy az (5.3) összefüggésben szereplő indexek az alábbi határok között érvényesek: (0 ≤ k < n ) , (− m ≤ λ ≤ m) , (0 ≤ i ≤ m) , (0 ≤ k − i < n ) , [0 ≤ (k − i ) + m < n ], [0 ≤ j = λ + i ≤ n ]. Az (5.1)-(5.3)-(5.4) összefüggésekben szereplő ε k simítási paraméter értéke pontról pontra változhat vagy állandó érték is lehet, amelyet célszerű úgy megválasztani, hogy a n

K = ∑ (ri − ηi )

2

(5.5)

i =0

simított és a simítandó pontok különbségének négyzetösszege kis érték legyen. A Nyíri-féle simító eljárást és a lineáris egyenletrendszerek megoldására vonatkozó új direkt módszert a [81], [82], [83] munkák tartalmazzák.

56

5.3. AZ ELSŐ ÉS MÁSODIK DERIVÁLTAK SZÁMÍTÁSA A DISZKRÉT PONTOKBAN A diszkrét pontokbeli első és második deriváltakat a negyedfokú Newton-féle interpolációs polinom alkalmazásával számítjuk, amelynek együtthatói az (5.2) összefüggés szerinti osztott differenciák. Az első öt ponttól az utolsó öt pontig a következő

[

{

]}

Pk (z ) = rk + (z − z k )⋅ δ rk + (z − z k +1 ) ⋅ δ 2 rk + (z − z k + 2 )⋅ δ 3 rk + δ 4 rk ⋅ (z − z k +3 )

(5.6)

haladó negyedfokú Newton-polinom érvényes, és a tömörebb írásmód érdekében bevezetjük a Pk (z ) = p1 (z )

(5.7)

jelölést, ahol (0 ≤ k < n − 5) esetén p1 (z ) = rk + p2 (z ) ⋅ (z − z k );

p2 (z ) = δ rk + p3 (z )⋅ (z − z k +1 ) ;

p3 (z ) = δ 2 rk + p4 (z )⋅ (z − z k + 2 );

p 4 (z ) = δ 3 rk + δ 4 rk ⋅ (z − z k +3 ) .

A diszkrét pontokban számított első deriváltakat az (5.7) kifejezésből Pk′ (z ) = p1′ (z )

(5.8)

differenciálással származtatjuk, ahol p1′ (z ) = p2 (z ) + p′2 (z ) ⋅ (z − z k ) ; p3′ (z ) = p4 (z ) + p′4 (z ) ⋅ (z − z k + 2 ) ;

p2′ (z ) = p3 (z ) + p3′ (z ) ⋅ (z − z k +1 ) p4′ (z ) = δ 4 rk ;

és a diszkrét pontokban számított második deriváltakat az (5.8) kifejezésből Pk′′(z ) = p1′′(z )

(5.9)

szintén differenciálással származtatjuk, ahol p1′′(z ) = 2 p′2 (z ) + p′2′ (z )⋅ (z − z k ); p3′′ (z ) = 2 p ′4 (z ) + p ′4′ (z )⋅ (z − z k + 2 ) ;

p′2′ (z ) = 2 p3′ (z ) + p3′′(z ) ⋅ (z − z k +1 ); p′4′(z ) = 0 .

Az utolsó öt pontban az első és második deriváltak számítása esetén az (5.6) negyedfokú Newton-féle interpolációs polinom következő

{

[

]}

Qk (z ) = rk + (z − z k )⋅ δ rk −1 + (z − z k −1 )⋅ δ 2 rk − 2 + (z − z k − 2 )⋅ δ 3 rk −3 + δ 4 rk −4 ⋅ (z − z k −3 )

(5.10)

retrográd alakjából indulunk ki, és a tömörebb írásmód érdekében bevezetjük a jelölést, ahol (n − 5 ≤ k ≤ n ) esetén

Q k (z ) = q1 (z )

(5.11)

57

q1 (z ) = rk + q2 (z ) ⋅ (z − z k );

q2 (z ) = δ rk −1 + q3 (z )⋅ (z − z k −1 );

q3 (z ) = δ 2 rk −2 + q4 (z )⋅ (z − z k −2 ) ;

q4 (z ) = δ 3 rk −3 + δ 4 rk −4 ⋅ (z − z k −3 ) .

Az utolsó öt diszkrét pontban számított első deriváltakat az (5.11) kifejezésből Qk′ (z ) = q1′ (z )

(5.12)

differenciálással származtatjuk, ahol q1′ (z ) = q2 (z ) + q′2 (z )⋅ (z − z k ) ;

q ′2 (z ) = q3 (z ) + q 3′ (z ) ⋅ (z − z k −1 ) ;

q3′ (z ) = q4 (z ) + q′4 (z ) ⋅ (z − z k − 2 );

q4′ (z ) = δ 4rk − 4 ;

és az utolsó öt diszkrét pontban számított második deriváltakat az (5.12) kifejezésből Qk′′ (z ) = q1′′(z )

(5.13)

szintén differenciálással származtatjuk, ahol q1′′(z ) = 2q′2 (z ) + q′2′ (z ) ⋅ (z − z k ) ; q3′′(z ) = 2q′4 (z ) + q′4′ (z )⋅ (z − z k −2 );

q2′′ (z ) = 2q3′ (z ) + q3′′(z ) ⋅ (z − z k −1 ) ; q′4′ (z ) = 0 .

A diszkrét pontokban vett első és második deriváltak ezúton történő meghatározása a későbbi számítások szempontjából előnyös, mert analitikus függvénnyel végzett ellenőrző számításaink szerint az első és második deriváltak 10−9÷10−10 nagyságrendben megegyeznek. Az ellenőrző számításokat a Függelék F.IV. pontja tartalmazza. 5.4. A HERMITE-POLINOM ALKALMAZÁSA A simított diszkrét pontokban ismerjük az első és második deriváltakat, ezért ezekre a hat feltételt kielégítő másodrendben csatlakozó ötödfokú Hermite-polinomot illesztünk. Az (rk , z k ) és az (rk +1 , zk +1 ) pontokon áthaladó χ (z ) ötödfokú Hermite-polinom a χ ′(zk ) = rk′ és a χ ′(z k +1 ) = rk′+1 , valamint a χ ′′(z k ) = rk′′ és a χ ′′(z k +1 ) = rk′′+1 feltételeknek eleget tesz:

χ (z ) = rk + δ rk (z − z k ) + δ 2 rk (z − z k )2 + δ 3 rk (z − z k )3 + δ 4 rk (z − z k )3 (z − z k +1 ) +

+ δ 5 rk (z − z k )3 (z − z k +1 )2 ,

(5.14)

amelyből differenciálással a következő

[

]

χ ′(z ) = δ rk + 2δ 2 rk (z − z k ) + 3δ 3 rk (z − z k )2 + δ 4 rk 3(z − z k )2 (z − z k +1 ) + (z − z k )3 +

[

]

+ δ 5 rk 3(z − z k ) (z − z k +1 ) + 2(z − z k ) (z − z k +1 ) 2

2

3

(5.15)

összefüggés adódik, amelyből ismét differenciálással a

58

[

]

χ ′′(z ) = 2δ 2 rk + 6δ 3 rk (z − z k ) + δ 4 rk 6(z − z k )(z − z k +1 ) + 6(z − z k +1 )2 +

[

+ δ rk 6(z − z k )(z − z k +1 ) + 12(z − z k ) (z − z k +1 ) + 2(z − z k ) 5

2

2

3

]

(5.16)

ötödfokú Hermite-polinom második differenciálhányadosát kapjuk. Az (5.14)-(5.16) összefüggésekben szereplő osztott differenciák a Hermite-polinom feltételrendszeréből adódnak:

δ rk = rk′ , 1 δ 2 rk = rk′′ , 2 − r r rk′ rk′′ δ 3 rk = k +1 k 3 − − , 2 (z k +1 − zk ) (z k +1 − zk ) 2(zk +1 − zk ) r −r rk′′ 2r ′ + r ′ δ 4 rk = −3 k +1 k 4 + k k +1 3 − , (zk +1 − zk ) (zk +1 − zk ) 2(zk +1 − zk )2 r −r r′ + r′ r ′′ − r ′′ δ 5 rk = 6 k +1 k 5 − 3 k k +1 4 + k +1 k 3 . (z k +1 − z k ) (z k +1 − z k ) 2(z k +1 − z k ) Az előzőek alapján egy tetszőleges z pontban a görbület az alábbi g (z ) =

χ ′′(z )

[1 + χ ′(z ) ]

2 3/ 2

(5.17)

formula szerint számítható. A koordinátahálózat generálása során, a másodrendű folytonosságot a hat feltételt kielégítő ötödfokú Hermite-polinom alkalmazásával biztosítottuk. 5.5. KÉT KONTÚR KÖZÉ BEÍRHATÓ KÖR KÖZÉPPONTJÁNAK SZÁMÍTÁSA Az 5.1. ábra alapján, ötödfokú Hermite-polinomok által adottak a χ1 (z ) és a χ 2 (z ) egymást nem metsző görbék, és keressük a χ1 (z1 ) és χ 2 (z 2 ) értinési pontokon áthaladó kör (r0 , z0 ) középpontját. Az iteratív számítások során a χ1 (z ) vezérgörbe χ1 (z1 ) pontját és ezáltal a z1 értéket rögzítettnek tekintjük, és az adott χ 2 (z ) görbén keressük χ 2 (z 2 ) érintési pontot, tehát a z 2 értéket a kezdeti értékének ismeretén túl ismeretlenként kezeljük.

59

5.1. ábra A két kontúr közé beírható kör geometriai adatai Az r0 , z 0 , z 2 ismeretlenek esetén a χ1 (z ) és a χ 2 (z ) görbék közé beírható kör érintési pontjaiban vett érintők iránytangensére vonatkozóan a következő z1 − z 0 , r0 − χ1 (z1 ) z0 − z 2 tgα 2 = χ 2′ (z 2 ) = , χ 2 (z 2 ) − r0 tgα1 = χ1′ (z1 ) =

(5.18) (5.19)

2 2 2 2 egyenletek írhatók fel. Az R 2 = (z1 − z 0 ) + [χ1 (z1 ) − r0 ] és az R 2 = (z 2 − z 0 ) + [χ 2 (z 2 ) − r0 ] beírható kör esetén érvényes kör egyenletei alapján az alábbi

(z1 − z0 )2 + [χ1 (z1 ) − r0 ]2 = (z2 − z0 )2 + [χ 2 (z 2 ) − r0 ]2

(5.20)

egyenlőség áll fenn. Az (5.18)-(5.19)-(5.20) egyenleteket zérusra rendezve a következő f1 (r0 , z 0 , z 2 ) = χ1′ (z1 ) ⋅ [r0 − χ1 (z1 )] + z 0 − z1 = 0 , f 2 (r0 , z 0 , z 2 ) = χ 2′ (z 2 ) ⋅ [χ 2 (z 2 ) − r0 ] − z 0 + z 2 = 0 ,

f 3 (r0 , z 0 , z 2 ) = z 22 − z12 + 2 z 0 (z1 − z 2 ) + 2r0 ⋅ [χ1 (z1 ) − χ 2 (z 2 )] − χ12 (z1 ) + χ 22 (z 2 ) = 0 , 60

egyenletek írhatók fel, és a megoldás közelítő értékeit a Newton-Raphson módszer segítségével keressük az alábbi lineáris egyenletrendszer alakjában: * * J (i ) ⋅ t (i ) = −f (i ) (r0 , z 0 , z 2 ) ,

(5.21)

ahol

J (i )

 ∂ f 1 (i )   ∂ r0  ∂ f (i ) = 2  ∂ r0 (i )  ∂ f3   ∂ r0

∂ f1 ( ) ∂ z0 i ∂ f2 ( ) ∂ z0 i ∂ f3 ( ) ∂ z0 i

∂ f1 ( )   ∂ z2  i ∂ f2 ( )   ∂ z2  i ∂ f3 ( )   ∂ z2  i

az i-edik iterációs lépésben számított Jacobi-mátrix. Az (5.21) lineáris egyenletrendszert a teljes főelemkiválasztással történő Gauss-elimináció segítségével oldjuk meg. Az r0 , z 0 , z 2 ismeretleneket az i+1-edik iterációs lépésben az alábbi r0(i +1) = t1(i ) + r0(i ) ;

z 0(i +1) = t 2(i ) + z 0(i ) ;

z 2(i +1) = t 3(i ) + z 2(i ) ;

formulák szerint számítjuk, és a Newton-Raphson iterációt addig folytatjuk, amíg az i+1-edik és az i-edik iterációs lépésben számított ismeretlen értékek különbségének abszolút értéke kisebb lesz egy adott ε küszöbértéknél: r0(i +1) − r0(i ) < ε ;

z 0(i +1) − z 0(i ) < ε ;

z 2(i +1) − z 2(i ) < ε .

Az (r,z) meridiánsíkon a tengelyirányú görbevonalú koordinátavonalakat a beírható körök középpontjainak számításával határozzuk meg, amelyeket simítunk, kiszámítjuk a pontbeli első és második deriváltakat, és a pontokra ötödfokú Hermite-polinomot illesztünk. Az eljárást addig folytatjuk, amíg kellő mennyiségű tengelyirányú koordinátavonalat kapunk. 5.6. AZ ORTOGONÁLIS TRAJEKTÓRIÁK MEGHATÁROZÁSA Az (r,z) meridiánsíkon a tengelyirányú koordinátavonalak mentén ismert az iránymező, mert minden pontban ismert az első és második derivált értéke.

5.2. ábra Az ortogonális trajektóriák meghatározása 61

Az 5.2. ábra alapján az első lépésben elindulunk a kiválasztott k-adik koordinátavonalra rT′ (z ) = −

1 rk′ (z )

(5.22)

merőleges iránnyal a következő k+1-edik koordinátavonalig, és a kapott távolságot megfelezzük. A felezési pontból a k+1-edik koordinátavonalra merőleges iránnyal haladunk tovább a k+1-edik koordinátavonal metszéspontjáig. A kapott pontot trajektória pontnak tekintjük. Ezt az eljárást a meridiánmetszetet határoló két peremgörbe között elhelyezkedő összes koordinátavonal esetén elvégezzük. A számított trajektória pontokat simítjuk, a pontokbeli első és második deriváltakat kiszámítjuk, és a kapott pontokra ötödfokú Hermite-polinomot illesztünk. Az ortogonális trajektóriák számítását célszerű a középvonaltól indítani. 5.7. A SZÁMÍTÁSI TAPASZTALATOK ÖSSZEFOGLALÁSA A Nyíri-féle simító eljárás fontos eleme a kidolgozott numerikus eljárásnak, mert a másodrendű folytonosság előírásával a pontokbeli első és második deriváltak ismeretét nagy pontossággal megköveteljük. A számítási tapasztalataink azt mutatják, hogy nem analitikusan adott függvény diszkrét pontjaiban számított első és második deriváltak simító eljárás alkalmazása nélkül pontatlanok, és ebből adódóan a kapott pontokra másodrendben illeszkedő ötödfokú Hermite-polinom oszcillációt mutat. A Nyíri-féle simító eljárás alkalmazásával az előbb említett oszcilláció elkerülhető, az első és második deriváltak 10−9 nagyságrendben pontosak leszek, és ezáltal a hálógenerálás során az öröklődő numerikus hibák mértékét csökkentjük. A kidolgozott kétdimenziós másodrendben folytonons hálógeneráló eljárás geometriai adatrendszer létrehozására alkalmas, amely később felhasználható egy adott áramlástani probléma diszkretizált egyenleteinek megoldása során. A számítások végeredményeképpen a cellákra bontott meridiánmetszet minden pontjában az első és második deriváltak által a koordinátavonalak görbületei ismertek lesznek. Az eljárás előnye, hogy a határoló peremgörbéket kellő mennyiségű diszkrét ponttal kell megadni. A módszer számítástechnikai szempontból könnyen algoritmizálható, hatékonyan bővíthető, ezért a mérnöki gyakorlatban jól alkalmazható. A 3. fejezetben tárgyalt turbulens áramlási feladat numerikus megoldásához szükséges geometriai adatrendszert az ismertetett kétdimenziós másodrendben folytonos hálógeneráló eljárással hoztuk létre. A feladatot saját készítésű C nyelvű programcsomag segítségével oldottuk meg. Az elkészített forráskód a GAMBIT szoftver számára is képes adatrendszert szolgáltatni, így a kidolgozott eljárás alapján létrehozott másodrendben folytonos koordinátahálózaton a FLUENT programrendszerrel áramlási feladatok megoldása is lehetővé vált. A Függelék V. pontjában a hálógenerálásra vonatkozó példák találhatók.

62

6. ÖSSZEFOGLALÁS Az értekezés a Czibere-féle új sztochasztikus turbulencia-modell [44] alapján kidolgozott numerikus eljárást ismertet, amely két egymást nem metsző koaxiális forgásfelület által határolt forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens áramlások sebességterének számítására alkalmas. A feladatot a középsebességtérben érvényes turbulens örvénydiffúzió új, Czibere-féle parciális differenciálegyenletének segítségével oldottuk meg. A kontinuitási egyenlet alapján bevezetjük az áramfüggvényt és az örvényvektor zérustól különböző kerületi irányú komponense egyenletéből az áramfüggvényre egy elliptikus típusú inhomogén parciális differenciálegyenletet nyerünk, és ezzel az áramlási feladat numerikus megoldására alkalmas parciális differenciálegyenlet-rendszer adódik. Az örvénydiffúzió differenciálegyenletében fellépő domináns turbulens nyírófeszültséget a sztochasztikus turbulencia-modell algebrai egyenletével számítottuk. A numerikus eljárást koaxiális forgásfelületek által határolt forgásszimmetrikus térben kialakuló perdületmentes turbulens áramlások esetén teszteltük. A kapott eredményeket mérésekkel is összehasonlítottuk. Az áramlási feladat alapegyenleteit ortogonális görbevonalú koordinátarendszerben írtuk fel és a meridiánsíkon a véges differenciák módszerén alapuló számítások megkönnyítésének érdekében ívhosszkoordinátákat használtunk. Az ívhosszkoordinátákkal felírt parciális differenciálegyenlet-rendszerben az áramlási tér koordinátafelületeinek görbületei szerepelnek, ezért a görbületek számítására alkalmas kétdimenziós másodrendben folytonos ortogonális görbevonalú koordinátahálózat generáló eljárást dolgoztunk ki, amely a számításokhoz szükséges geometriai adatrendszert előállítja. A számítások eredményeképpen a számítási hálózat diszkrét pontjainak koordinátái és az ívelemek hossza mellett a pontokbeli első és második deriváltak által a koordinátavonalak görbületei is numerikusan ismertek. Az alkalmazott új sztochasztikus turbulencia-modell alapján kidolgozott numerikus eljárás nemcsak forgásszimmetrikus, hanem kétdimenziós áramlások számítására is alkalmas. Ha az ortogonális görbevonalú koordinátarendszerben felírt alapegyenletekben a harmadik Lamé-féle metrikus együtthatót egynek választjuk, akkor a kétdimenziós áramlási feladatok esetén érvényes parciális differenciálegyenlet-rendszerhez jutunk. A párhuzamos síkfalak között kialakuló viszkózus folyadék turbulens áramlása esetén végzett számításaink eredményeit [P.12] a Janiga által kidolgozott aszimptotikus megoldással [48] hasonlítottuk össze. Az ismertetett numerikus eljárás általánosan alkalmazható forgásszimmetrikus térben kialakuló perdületmentes stacionárius áramlási feladatok megoldásának esetén. A kapott számítási eredmények más számítási módszerekkel és mérési eredményekkel jó egyezést mutatnak. A jövőben ezért célszerűnek mutatkozik a kétdimenziós másodrendben folytonos ortogonális görbevonalú koordinátahálózat generáló eljárás háromdimenziós esetre történő kiterjesztése és a sztochasztikus turbulencia-modellt háromdimenziós áramlások számítására alkalmazni. Az értekezés elkészítése több tudományterület együttes alkalmazását követelte meg, ezért köszönet illeti a Sályi István Gépészeti Tudományok Doktori Iskola, Gépészeti Tudományok tématerületén oktatókat, akik előadásaikkal hozzájárultak e munka elkészítéséhez. Külön köszönet illeti Dr. Czibere Tibor akadémikust, aki a kutatómunka irányítását és a folyamatos konzultációt biztosította, Dr. Nyíri András egyetemi tanárt, aki a hálógeneráló eljárás kidolgozásában segített, Dr. habil. Szabó Szilárd tanszékvezető egyetemi tanárt, aki a kutatáshoz szükséges színvonalas munkakörülményeket teremtett. Köszönet illeti az Áramlás- és Hőtechnikai Gépek Tanszéke valamennyi munkatársát, akik gondolataikkal és hasznos szakmai tanácsaikkal az értekezés elkészítését segítették. Végül külön köszönet illeti a szüleimet, akik minden tekintetben stabil családi háttérként voltak jelen e disszertáció létrejöttének folyamatában. 63

7. ÚJ TUDOMÁNYOS EREDMÉNYEK I. A turbulens örvénydiffúzió Czibere-féle parciális differenciálegyenletének [44] megoldásával numerikus eljárás kidolgozása két egymást nem metsző koaxiális forgásfelület által határolt forgásszimmetrikus térben kialakuló perdületmentes stacionárius turbulens nyíróáramlás sebességterének számítására. Az eljárás a lamináris áramlásra alkalmazott, az irodalomból ismert módszer kiterjesztése turbulens áramlás esetére. A módszer hatékonyságát más módszerekkel és a mérési eredményekkel való összehasonlítás mutatja. I.1. Az örvénydiffúzió differenciálegyenletét a kontinuitási egyenlettel kiegészítve nyert parciális differenciálegyenlet-rendszer megoldását - az új turbulens tag figyelembevételével - véges differenciák módszerével diszkretizált két lineáris algebrai egyenletrendszer numerikus megoldására vezettük vissza. A turbulens tagban szereplő domináns turbulens nyírófeszültséget a turbulencia-modell algebrai egyenletével számítottuk. A numerikus megoldás végrehajtására C nyelvű forráskód készült. I.2. A kidolgozott numerikus eljárás alkalmas: a) forgásszimmetrikus és perdületmentes lamináris áramlások sebességterének számítására; b) lamináris és turbulens kétdimenziós áramlások sebességterének meghatározására; c) forgásszimmetrikus illetve speciális esetben kétdimenziós potenciálos áramlások sebességterének számítására. II. Tekintettel arra, hogy ortogonális görbevonalú koordinátarendszert használunk és a diszkretizált lineáris algebrai egyenletrendszerekben fellépnek a koordinátavonalak görbületei, szükségessé vált egy alkalmas hálógeneráló eljárás kidolgozása, amely a hálózati pontok létrehozása mellett ezeket a görbületeket is szolgáltatja. II.1. Új numerikus eljárás kidolgozása kétdimenziós másodrendben folytonos ortogonális görbevonalú koordinátahálózat létrehozására iteráció alkalmazásával.

64

FONTOSABB JELÖLÉSEK AS , AW , AP , AE , AN a, b, c, d, e, f E F FR FT F' G G* g H H* h H1, H2, H3 I k l p pT p' q1 , q 2 , q3 q1′ , q 2′ , q 3′ r r R R12, R21, R13, R23 Re s1 s2 t v vT v' vm v~ α , β , γ ,ϑ , µ α* , β* ,γ * ε η κ υ υT ρ σ

a megoldandó lineáris algebrai egyenletrendszerek együtthatói a tömörebb írásmód érdekében bevezetett jelölések a transzformáció tenzora feszültségtenzor Reynolds-féle látszólagos turbulens feszültségtenzor pillanatnyi turbulens feszültségtenzor feszültségtenzor az ingadozások sebességterében a H tenzor transzformáltja a számítási koordinátarendszerbe a H* tenzor transzformáltja a számítási koordinátarendszerbe a térfogati erők összessége a turbulencia Czibere-féle hasonlósági tenzora a H hasonlósági tenzor deviátora a meridiáncsatorna szélessége egy adott keresztmetszetben a Lamé-féle metrikus együtthatók egységtenzor turbulens kinetikus energia az alkalmazott turbulencia-modell léptékfüggvénye a turbulens középnyomás pillanatnyi turbulens nyomás turbulens nyomásingadozás koordinátairányok a számítási koordinátarendszerben koordinátairányok a természetes koordinátarendszerben helyvektor a vizsgált pont forgástengelytől mért távolsága a körkeresztmetszetű cső sugara a forgásszimmetrikus tér koordinátafelületeinek görbületi sugarai Reynolds-szám az áramlás irányába eső ívhosszkoordináta az áramlás irányára merőleges ívhosszkoordináta idő turbulens középsebesség pillanatnyi turbulens sebesség a turbulens sebesség ingadozása sebességmaximum a csatorna keresztmetszetére vonatkoztatott átlagsebesség a turbulencia H hasonlósági tenzorának elemei a H* tenzor főátló menti elemei disszipációsebesség a folyadék dinamikai viszkozitása Kármán-konstans a folyadék kinematikai viszkozitása a Boussinesq-féle örvényviszkozitás a folyadék sűrűsége súrlódási feszültségtenzor 65

σR ξ Θ Π Ω Ψ

Reynolds-féle látszólagos turbulens feszültségtenzor deviátora trajektória menti paraméter a léptékfüggvény esetén domináns turbulens nyírófeszültség teljes potenciál örvényerősség, amely a rot v vektor nagyságát fejezi ki a meridiánsíkon bevezetett áramfüggvény

Indexek i j SW S SE W P E NW N NE

vízszintes irányú rácsindex függőleges irányú rácsindex South-West (Dél-nyugat) South (Dél) South-East (Dél-kelet) West (Nyugat) a számítási hálózat vizsgált csomópontja East (Kelet) North-West (Észak-nyugat) North (Észak) North-East (Észak-kelet)

66

FÜGGELÉK F.I. A (3.3) differenciálegyenlet levezetése. A (3.2) összenyomhatatlan közeg esetén érvényes kontinuitási egyenlet alapján a meridiánsíkon bevezetjük az áramfüggvényt, amellyel a sebességkomponensek: v1 =

1 ∂Ψ , H 3 H 2 ∂ q2

és

v2 = −

1 ∂Ψ . H 3 H 1 ∂ q1

Ezeket behelyettesítve a (3.1) örvényerősség nagyságát kifejező egyenletbe a következő 1 H1H 2

 ∂  H 2 ∂Ψ  ∂   +   ∂ q1  H 3 H 1 ∂ q1  ∂ q 2

 H 1 ∂Ψ    = −Ω ,  H 3 H 2 ∂ q 2 

differenciálegyenletet nyerjük, és a differenciálások elvégzése után a  ∂ H 2 1 ∂Ψ H 2 ∂  1 ∂Ψ  1 ∂ H 3 H 2 ∂Ψ   − 2 + +  q H H q H q H q q H q ∂ ∂ ∂ ∂ ∂ ∂ H 1 3 1  1 1  1 1 1  1 3 1 3 H ∂  1 ∂Ψ  1 ∂ H 3 H 1 ∂Ψ  ∂ H 1 1 ∂Ψ  − + + 1  = −Ω ∂ q 2 H 3 H 2 ∂ q 2 H 3 ∂ q 2  H 2 ∂ q 2  H 32 ∂ q 2 H 2 ∂ q 2  1 H1H 2

alakban adódik. A szögletes zárójelet felbontva ∂  1 ∂Ψ  1 ∂ H 2 1 ∂Ψ 1 1 ∂ H3 1   − 2 + H 2 H 1 ∂ q1 H 3 H 1 ∂ q1 H 3 H 1 ∂ q1  H 1 ∂ q1  H 3 H 1 ∂ q1 H 1 ∂  1 ∂Ψ  1 ∂ H 1 1 ∂Ψ 1 1 ∂ H3   − 2 + + H1 H 2 ∂ q 2 H 3 H 2 ∂ q2 H 3 H 2 ∂ q2  H 2 ∂ q2  H 3 H 2 ∂ q2

∂Ψ + ∂ q1 1 ∂Ψ = −Ω , H 2 ∂ q2

ezután az egyenlet mindkét oldalát H 3 -mal megszorozva ∂ 2Ψ

(H 1 ∂ q1 )

2

+

+

1 ∂ H 2 1 ∂Ψ 1 ∂ H 3 1 ∂Ψ − + H 2 H 1 ∂ q1 H 1 ∂ q1 H 3 H 1 ∂ q1 H 1 ∂ q1

∂ 2Ψ 1 ∂ H 1 1 ∂Ψ 1 ∂ H 3 1 ∂Ψ + − = −Ω H 3 , 2 (H 2 ∂ q2 ) H1 H 2 ∂ q2 H 2 ∂ q2 H 3 H 2 ∂ q2 H 2 ∂ q2

és az összevonások után a (3.3) differenciálegyenletet nyerjük:  1 ∂ H2 ∂ 2Ψ ∂ 2Ψ 1 ∂ H 3  ∂Ψ  + +  − + 2 2 (H 1∂ q1 ) (H 2 ∂ q 2 )  H 2 H 1 ∂ q1 H 3 H 1 ∂ q1  H 1∂ q1  1 ∂ H1 1 ∂ H 3  ∂Ψ  +  − = −Ω H 3 .  H 1 H 2 ∂ q2 H 3 H 2 ∂ q2  H 2 ∂ q 2 67

F.II. A (3.1) örvényerősség nagyságát kifejező összefüggés kifejtése tagonként. 1  ∂ (v 2 H 2 ) ∂ (v1 H 1 )  − =  ∂ q 2  H 1 H 2  ∂ q1 ∂v ∂ H2 ∂v ∂ H1  1   H 2 2 + v 2 = = − H 1 1 − v1 H1H 2  ∂ q1 ∂ q1 ∂ q2 ∂ q 2  v2 ∂ H 2 v1 ∂ H 1 1 ∂ v2 1 ∂ v1 = − + − . H 1 ∂ q1 H 2 ∂ q 2 H 2 H 1 ∂ q1 H 1 H 2 ∂ q 2

Ω=

A 3.1. táblázat és a 3.3. ábra alapján az előző összefüggés ívhosszkoordinátákkal felírt alakja:

Ω=

∂ v 2 ∂ v1 v 2 v − + − 1 . ∂ s1 ∂ s 2 R12 R 21

F.III. A (3.6) differenciálegyenlet konvektív és viszkózus tagjának kifejtése. A (3.6) differenciálegyenlet konvektív tagjának kifejtése: 1  ∂ (v1Ω H 2 ) ∂ (v 2 Ω H 1 ) + = H 1 H 2  ∂ q1 ∂ q2  ∂ v1 ∂ H2 ∂ v2 ∂ H1  ∂Ω ∂Ω 1   Ω H 2 = = + v1 H 2 + v1 Ω + Ω H1 + v2 H 1 + v2 Ω H1H 2  ∂ q1 ∂ q1 ∂ q1 ∂ q2 ∂ q2 ∂ q 2  1 ∂ v1 1 ∂Ω 1 ∂ H2 =Ω + v1 + v1 Ω + H 1 ∂ q1 H 1 ∂ q1 H 2 H 1 ∂ q1 1 ∂ v2 1 ∂Ω 1 ∂ H1 +Ω + v2 + v2 Ω = H 2 ∂ q2 H 2 ∂ q2 H1H 2 ∂ q2 = v1

v1 ∂ H 2 v2 ∂ H 1  1 ∂Ω 1 ∂ Ω  1 ∂ v1 1 ∂ v2 Ω . + v2 +  + + + H 1 ∂ q1 H 2 ∂ q 2  H 1 ∂ q1 H 2 ∂ q 2 H 2 H 1 ∂ q1 H 1 H 2 ∂ q 2  A (3.6) differenciálegyenlet viszkózus tagjának kifejtése:

 ∂  H 2 ∂ (Ω H 3 ) ∂  H 1 ∂ (Ω H 3 )    +   =  ∂ q1  H 3 H 1 ∂ q1  ∂ q 2  H 3 H 2 ∂ q 2   H Ω ∂ H 3  υ  ∂  H 2 ∂ Ω H 2Ω ∂ H 3  ∂  H1 ∂ Ω   +   = = + + 1  H 1 H 2  ∂ q1  H 1 ∂ q1 H 3 H 1 ∂ q1  ∂ q 2  H 2 ∂ q 2 H 3 H 2 ∂ q 2   υ  ∂ H 2 1 ∂ Ω H 2 ∂ 2Ω ∂ H 2 Ω ∂ H 3 ∂ Ω H 2 ∂ H 3 = + + + +  H1 H 2  ∂ q1 H1 ∂ q1 H1 ∂ q12 ∂ q1 H 3 H1 ∂ q1 ∂ q1 H 3 H1 ∂ q1

υ H1H 2

+

H 2Ω ∂  1 ∂ H 3  H 3 ∂ q1  H 1 ∂ q1

 H 2Ω 1 ∂ 2 H 3 ∂ H1 Ω ∂ H 3 ∂ Ω H1 ∂ H 3  − + + + H 1 H 32 ∂ q12 ∂ q2 H 3 H 2 ∂ q2 ∂ q2 H 3 H 2 ∂ q2  68

 1 ∂ H 3  H 1Ω 1 ∂ 2 H 3    − = 2 2   H 2 ∂ q2  H 2 H 3 ∂ q2   ∂ 2 Ω  1 ∂ H2 ∂ 2Ω 1 ∂ H3  1 ∂ Ω  +  + + =υ  + 2 2  (H 1∂ q1 ) (H 2 ∂ q2 )  H 2 H1 ∂ q1 H 3 H1 ∂ q1  H1 ∂ q1

+

H 1Ω ∂ H 3 ∂ q2

 1 ∂ H1 ∂2H3 1 ∂ H3  1 ∂ Ω Ω  ∂2H3  +  + + + +  2 (H 2 ∂ q 2 )2  H 1 H 2 ∂ q 2 H 3 H 2 ∂ q 2  H 2 ∂ q 2 H 3  (H 1∂ q1 )  1 ∂ H2 1 ∂ H 3  1 ∂ H 3  1 ∂ H1 1 ∂ H 3  1 ∂ H 3     +  − +  −  .  H 2 H 1 ∂ q1 H 3 H 1 ∂ q1  H 1 ∂ q1  H 1 H 2 ∂ q 2 H 3 H 2 ∂ q 2  H 2 ∂ q 2   A (3.11) turbulens örvénydiffúzió parciális differenciálegyenlete ívhosszkoordinátákkal az előbb kifejtett konvektív és viszkózus tagok alapján van felírva.

69

F.IV. Numerikus differenciálás a Newton-féle interpolációs polinom alkalmazásával, öt ponton átmenő negyedfokú parabola segítségével. A pontok sorszáma 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40.

x koordináták

y = f (x ) = 5 x 4 + 6 x + 5 5,0000000000 5,2290914810 5,4919608151 5,7797207517 6,0876894179 6,4117132968 6,7475717833 7,0908442746 7,4369631695 7,7814737551 8,1202366266 8,4497572653 8,7672646963 9,0707349893 9,3587923546 9,6305328592 9,8854442031 10,1231878112 10,3436056049 10,5466395223 10,7011057318 10,8290614120 10,9693659838 11,1222495918 11,2877307630 11,4655569695 11,6551335616 11,8554868900 12,0653857049 12,2836781145 12,5098733766 12,7450211309 12,9926983694 13,2597384458 13,5561486603 13,8944175582 14,2886175764 14,7547427468 15,3129139837 16,0001300138

0,0000000000 0,0381801427 0,0819558734 0,1297175123 0,1803989875 0,2328363623 0,2857091235 0,3376434561 0,3873923212 0,4340109812 0,4769252639 0,5159195138 0,5510424628 0,5825091764 0,6106163468 0,6356811679 0,6580132567 0,6778888888 0,6955533574 0,7112184774 0,7227696574 0,7321101106 0,7421243537 0,7527766523 0,7640151256 0,7757703825 0,7879541537 0,8004607643 0,8131782359 0,8260101563 0,8389073268 0,8519095162 0,8651844541 0,8790443152 0,8939141249 0,9102689762 0,9285649635 0,9492308235 0,9727311152 1,0000050005 F1. táblázat Analitikus ellenőrzőfüggvény által felvett diszkrét pontok

70

y′ Abszolút hiba y ′ = 20 x 3 + 6 Analitikus számítás Numerikus számítás 6,0000000000 6,0000000000 0,0000000000 6,0011131217 6,0011131217 0,0000000000 6,0110095671 6,0110095671 0,0000000000 6,0436541795 6,0436541795 0,0000000000 6,1174173522 6,1174173522 0,0000000000 6,2524540906 6,2524540906 0,0000000000 6,4664470193 6,4664470193 0,0000000000 6,7698480371 6,7698480371 0,0000000000 7,1627410884 7,1627410884 0,0000000000 7,6350541856 7,6350541856 0,0000000000 8,1696065420 8,1696065420 0,0000000000 8,7464763245 8,7464763245 0,0000000000 9,3464565845 9,3464565845 0,0000000000 9,9531046320 9,9531046320 0,0000000000 10,5533944671 10,5533944671 0,0000000000 11,1374550199 11,1374550199 0,0000000000 11,6981506274 11,6981506274 0,0000000000 12,2302509797 12,2302509797 0,0000000000 12,7300973996 12,7300973996 0,0000000000 13,1951373513 13,1951373513 0,0000000000 13,5514392360 13,5514392360 0,0000000000 13,8480038867 13,8480038867 0,0000000000 14,1744783287 14,1744783287 0,0000000000 14,5315593842 14,5315593842 0,0000000000 14,9194046156 14,9194046156 0,0000000000 15,3374777657 15,3374777657 0,0000000000 15,7843694603 15,7843694603 0,0000000000 16,2577035417 16,2577035417 0,0000000000 16,7544259940 16,7544259940 0,0000000000 17,2716152891 17,2716152891 0,0000000000 17,8078807356 17,8078807356 0,0000000000 18,3654636254 18,3654636254 0,0000000000 18,9525750561 18,9525750561 0,0000000000 19,5850832680 19,5850832680 0,0000000000 20,2862220073 20,2862220073 0,0000000000 21,0847883021 21,0847883021 0,0000000000 22,0127850674 22,0127850674 0,0000000000 23,1058828065 23,1058828065 0,0000000000 24,4080768982 24,4080768982 0,0000000000 26,0003000315 26,0003000315 0,0000000000 F2. táblázat Az első deriváltak diszkrét pontokbeli abszolút hibája

A pontok sorszáma 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40.

71

y ′′ Abszolút hiba y ′′ = 60x 2 Analitikus számítás Numerikus számítás 0,0000000000 0,0000000000 0,0000000000 0,0874633978 0,0874633978 0,0000000000 0,4030059111 0,4030059111 0,0000000000 1,0095979798 1,0095979798 0,0000000000 1,9526276815 1,9526276815 0,0000000000 3,2527662965 3,2527662965 0,0000000000 4,8977821951 4,8977821951 0,0000000000 6,8401862068 6,8401862068 0,0000000000 9,0043686315 9,0043686315 0,0000000000 11,3019319081 11,3019319081 0,0000000000 13,6474624408 13,6474624408 0,0000000000 15,9703766832 15,9703766832 0,0000000000 18,2188677485 18,2188677485 0,0000000000 20,3590164354 20,3590164354 0,0000000000 22,3711393788 22,3711393788 0,0000000000 24,2454328334 24,2454328333 0,0000000000 25,9788867596 25,9788867596 0,0000000000 27,5720007335 27,5720007336 0,0000000000 29,0276683794 29,0276683793 0,0000000001 30,3499033557 30,3499033558 0,0000000001 31,3437586595 31,3437586595 0,0000000000 32,1591128426 32,1591128425 0,0000000001 33,0449133813 33,0449133814 0,0000000001 34,0003612949 34,0003612948 0,0000000001 35,0231467287 35,0231467288 0,0000000001 36,1091811819 36,1091811819 0,0000000000 37,2523049000 37,2523048999 0,0000000001 38,4442461110 38,4442461111 0,0000000001 39,6755306005 39,6755306004 0,0000000001 40,9375666986 40,9375666987 0,0000000000 42,2259301775 42,2259301774 0,0000000001 43,5449894275 43,5449894277 0,0000000001 44,9126483770 44,9126483769 0,0000000001 46,3631344851 46,3631344851 0,0000000000 47,9449477617 47,9449477618 0,0000000000 49,7153765419 49,7153765421 0,0000000001 51,7339734864 51,7339734863 0,0000000000 54,0623493769 54,0623493769 0,0000000000 56,7723493487 56,7723493487 0,0000000000 60,0006000615 60,0006000615 0,0000000000 F3. táblázat A második deriváltak diszkrét pontokbeli abszolút hibája

A pontok sorszáma 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40.

72

F1. ábra Az y = f (x ) = 5 x + 6 x + 5 ellenőrzőfüggvény grafikonja 4

A diszkrét pontokbeli numerikus differenciálásnál az ellenőrző számításaink során azért az y = f (x ) = 5 x 4 + 6 x + 5 ellenőrzőfüggvényt választottuk, mert a dolgozatban vizsgált meridiánmetszetet határoló peremgörbéket az F1. ábrán vázolt ellenőrzőfüggvényhez hasonló meredekségű görbeívek alkotják. Az F3. táblázatban foglalt számítások szerint az alkalmazott negyedfokú Newton-féle interpolációs polinom céljainknak megfelel, mert a diszkrét pontokban számított második deriváltak abszolút hibája 10 −9 ÷ 10 −10 nagyságrendű.

73

F.V. Példák a hálógenerálás esetén.

F2. ábra 1617 csomópontból álló koordinátahálózat

F3. ábra 2541 csomópontból álló koordinátahálózat 74

F4. ábra Hálógenerálás Venturi-csatorna esetén

F5. ábra Nagy osztásközű hálózat szárnyszelvény esetén 75

IRODALOMJEGYZÉK [1]

[2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

[13] [14] [15] [16]

[17]

O. Reynolds: On the dynamical theory of incompressible viscous fluids and the determination of the criterion, Philosophical Transactions of the Royal Society of London, Series A, Vol. 186, 123. 1895. J. Boussinesq: Théorie de l'Ecoulement tourbillant, Mem. Présentés par Divers Savants Acad. Sci. Inst. Fr., Vol. 23, pp. 46-50. 1877. H. Schlichting: Boundary Layer Theory, McGraw-Hill, New York, 1982. Theodore von Kármán: Mechanische Ähnlichkeit und Turbulenz, Proc. 3rd Congress of Applied Mechanics, pp. 86-93. 1930. G.I. Taylor: Statistical theory of turbulence, Proceedings of the Royal Society of London, Series A, Vol. 151, 421. 1935. A. Kolmogorov: Local structure of turbulence in incompressible viscous fluid for very large Reynolds number, Doklady Akademija Nauk, SSSR, Vol. 30, pp. 299-303. 1941. A. Kolmogorov: The equations of turbulent motion in an incompressible fluid, Izv. Sci. USSR Phys., 6, pp. 56-58. 1942. J. Smagorinsky: General circulation model of the atmosphere, Mon. Weather Rev., Vol. 91, pp. 99-164. 1963. B. Launder, D. Spalding: Mathematical Models of Turbulence, Academic Press, New York, 1972. B. Mohammadi, O. Pironneau: Analysis of the K-Epsilon Turbulence Model, Chichester, UK: Wiley, 1993. H.K. Versteeg, W. Malalasekera: An introduction to computational fluid dynamics, The finite volume method, Longman Scientific & Technical, London, pp. 41-84. 1995. M. Griebel, T. Dornseifer, T. Neunhoeffer: Numerical Simulation in Fluid Dynamics, A Practical Introduction, Society for Industrial and Applied Mathematics, Philadelphia, pp. 153-171. 1998. Hans Breuer: SH atlasz Fizika, Springer Hungarica, Bp., p. 11, 1993. Csermely Péter, Gergely Pál, Koltay Tibor, Tóth János: Kutatás és Közlés a Természettudományokban, Osiris Kiadó, Bp., 1999. T.-H. Shih, J. L. Lumley: A critical comparison of second order closures with direct numerical simulations of homogeneous turbulence, AIAA Journal 31(4):663-70. 1993. J. Eggels, F. Unger, M. Weiss, J. Westerweel, R. Adrian, R. Friedrich: Fully developed turbulent pipe flow: A comparison between direct numerical simulation and experiments, J. Fluid Mech., Vol. 268, pp. 175-209. 1994. H. A. Carlson, G. Berkooz, J. L. Lumley: Direct numerical simulation of flow in a channel with complex, time-dependent wall geometries: A pseudospectral method, Journal of Computational Physics 121:155-75. 1995. 76

[18] J.C. Rotta: Statistische Theorie nichthomogener Turbulenz, Zeitschrift für Physik, Vol. 129, pp. 547-572. 1951. [19] K. Hanjalic, B.E. Launder: A Reynolds stress model of turbulence and its application to thin shear flows, J. Fluid Mech., Vol. 52, pp. 609-638. 1972. [20] B.E. Launder, G.J. Reece, W. Rodi: Progress in the development of a Reynolds-stress turbulence closure, J. Fluid Mech., Vol. 68, pp. 537-566. 1975. [21] K. Hanjalic, B.E. Launder: Contribution towards a Reynolds-stress closure for lowReynolds number turbulence, J. Fluid Mech., Vol. 74, pp. 593-610. 1976. [22] V. Schumann: Realizabilty of Reynolds stress turbulence models, Phys. Fluids, Vol. 20, pp. 721-725. 1977. [23] M. Gibson, B. Launder: Ground effects on pressure fluctuations in the atmospheric boundary layer, J. Fluid Mech., Vol. 86, pp. 491-511. 1978. [24] M. Leschziner: Second moment closure for complex flows, In Proc. Int. Forum on Math. Mod. of Processes in Energy Systems, Sarajevo, 1989. [25] B. Launder: Phenomenological modelling: Present and future, In J. Lumley (Ed.), Whither Turbulence: Turbulence at the Crossroads, Springer-Verlag, Berlin, 1990. [26] C. Meneveau: Analysis of turbulence in the orthonormal wavelet representation, J. Fluid Mech., Vol. 232, pp. 469-520. 1991. [27] M. Farge: Wavelet transforms and their applications to turbulence, Ann. Rev. Fluid Mech., Vol. 24, pp. 395-457. 1992. [28] A.V. Johansson, M. Hallbäck: Modelling of rapid pressure strain in Reynolds-stress closures, J. Fluid Mech., Vol. 269, pp. 143-168. 1994. [29] S. Jakirlic, K. Hanjalic: On the performance of the second moment high and low Reynolds number closures in reattaching flows, In Int. Symposium on Turbulence, Heat and Mass Transfer, Lisboa, 1994. [30] W. Rodi: Current trends in turbulence modelling, Computational Aeronautical Fluid Dynamics, Clarendon Press, Oxford, pp. 225-245. 1994. [31] J.O. Hinze: Turbulence, McGraw-Hill, New York, 1959. [32] H. Tennekes, J. Lumley: A First Course in Turbulence, MIT Press, Cambridge, 1972. [33] J. Lumley: Toward a turbulent constitutive equation, J. Fluid Mech., Vol. 41, pp. 413434. 1970. [34] W. Rodi: A new algebraic relation for calculating Reynolds stresses, Z. Angew. Math. Mech., Vol. 56, pp. 219-221. 1976. [35] D.C. Wilcox, M.W. Rubesin: Progress in turbulence modelling for complex flow fields including effects of compressibility, NASA TP-1517, 1980. [36] C.G. Speziale: On nonlinear k-l and k-ε models of turbulence, J. Fluid Mech., Vol. 178, pp. 459-475. 1987. [37] T.-H. Shih, J. Zhu, J. L. Lumley: A new Reynolds Stress algebraic equation model, Computer Methods in Applied Mechanics and Engineering 125:287-302. 1995. 77

[38] L. Prandtl: Über die ausgebildete Turbulenz, Z. angew. Math. Mech., Vol. 5, pp. 136139. 1925. [39] B.S. Baldwin, H. Lomax: Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows, AIAA Paper, pp. 78-257, 1978. [40] T. Cebeci, J. Cousteix: Modeling and Computation of Boundary-Layer Flows, SpringerVerlag, Berlin, 1999. [41] L. Prandtl: Über ein neues Formelsystem für die ausgebildete Turbulenz, Nachr. Akad. Wiss. Göttingen, Math-Phys. Kl., pp. 6-19. 1945. [42] P.R. Spalart, S.R. Allmaras: A one-equation turbulence model for aerodynamic flows, AIAA Paper, Reno NV, pp. 92-439. 1992. [43] T. Czibere: Mechanische Ähnlichkeit der turbulenten Flüssigkeitsströmungen, Publ. Univ. of Miskolc. Series C, Mechanical Engineering, Vol. 50, pp. 3-12. 1999. [44] T. Czibere: Three Dimensional Stochastic Model of Turbulence, Journal of Computational and Applied Mechanics, Vol. 2., No. 1., pp. 7-20. 2001. [45] Czibere Tibor: Turbulencia Kutatások 1991-2001, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, Kézirat, 1991-2001. [46] Kalmár László: Numerikus módszer radiális szivattyúkban kialakuló súrlódásos áramlás meghatározására, Ph.D. értekezés, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, 1997. [47] Szabó Szilárd: Modellalkotási példák az áramlás- és hőtechnikai tudományterületen, Habilitációs értekezés, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, pp. 40-69. 2001. [48] Janiga Gábor: Kétdimenziós turbulens nyíróáramlások számítása sík, valamint enyhén görbült falakkal határolt csatornákban, Ph.D. értekezés, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, 2002. [49] A. Thom: The flow past circular cylinders at low speeds, Proc. Roy. Soc. London Ser. A, 141, pp. 651-666. 1933. [50] K. Aziz, J. Hellums: Numerical solution of the three-dimensional equations of motion for laminar natural convection, Phys. Fluids, 10, pp. 314-324. 1967. [51] G. Hirasaki, J. Hellums: Boundary conditions on the vector and scalar potentials in viscous three-dimensional hydrodynamics, Quart. Appl. Math., 28, pp. 293. 1970. [52] A.J. Chorin: Numerical study of slightly viscous flow, J. Fluid Mech., Vol. 57, pp. 785796. 1973. [53] G. Mallinson, G. de Vahl Davis: The method of the false transient for the solution of coupled elliptic equations, J. Comput. Phys., 12, pp. 435-461. 1973. [54] P. Roache: Computational Fluid Dynamics, Albuquerque: Hermosa, 1976. [55] G. Mallinson, G. de Vahl Davis: Three-dimensional natural convection in a box: A numerical study, J. Fluid Mech., 83, pp. 1-31. 1977. 78

[56] S. Richardson, A. Cornish: Solution of three-dimensional incompressible flow problems, J. Fluid Mech., 82, pp. 309-340. 1977. [57] S. Dennis, D. Ingham, R. Cook: Finite-difference methods for calculating steady incompressible flows in three dimensions, J. Comput. Phys., 33, pp. 325-339. 1979. [58] U. Ghia, K. Ghia, C. Shin: High-Re solutions for incompressible flow using the NavierStokes equations and a multigrid method, J. Comput. Phys., 48, pp. 387-411. 1982. [59] T. Gatski, C. Grosch, M. Rose: A numerical study of the two-dimensional Navier-Stokes equations in vorticity-velocity variables, J. Comput. Phys., 48, pp. 1-22. 1982. [60] A. Wong, J. Reizes: An effective vorticity-vector potential formulation for the numerical solution of three-dimensional duct flow problems, J. Comput. Phys., 55, pp. 98-114. 1984. [61] S. Dennis: Compact explicit finite-difference approximation to Navier-Stokes equations, In Soubbaramayer & J. Boujot (Eds.), 9th Internat. Conf. Numer. Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 128, Berlin, Heidelberg: Springer-Verlag, pp. 23-36. 1985. [62] T. Gatski, C. Grosch: Numerical study of the two- and three-dimensional unsteady Navier-Stokes equations in velocity-vorticity variables using compact difference schemes, In Soubbaramayer & J. Boujot (Eds.), 9th Internat. Conf. Numer. Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 128, Berlin, Heidelberg: SpringerVerlag, pp. 235-239. 1985. [63] T. Gatski, C. Grosch, M. Rose: The numerical solution of the Navier-Stokes equations for three-dimensional unsteady, incompressible flows by compact schemes, J. Comput. Phys., 82, pp. 298-329. 1989. [64] C.A.J. Fletcher: Computational Techniques for Fluid Dynamics, Vol. II. (Specific Techniques for Different Flow Categories), Springer-Verlag Berlin Heidelberg, pp. 373396. 1991. [65] Mo-Hong Chou, Weicheng Huang: Numerical Study of High-Reynolds-Number Flow Past a Bluff Object, International Journal for Numerical Methods in Fluids, Vol. 23, pp. 711-732. 1996. [66] M.T. Nair, T.K. Sengupta: Onset of Asymmetry: Flow Past Circular and Elliptic Cylinders, International Journal for Numerical Methods in Fluids, Vol. 23, pp. 13271345. 1996. [67] D.W. Peaceman, H.H. Rachford: The numerical solution of parabolic and elliptic differential equations, J. Soc. Ind. Appl. Math., Vol. 3, pp. 28-41. 1955. [68] J. Douglas, J.E. Gunn: "A general formulation of alternating direction methods: parabolic and hyperbolic problems", Numerische Mathematik, Vol. 6, pp. 428-453. 1964. [69] G.I. Marcsuk: A gépi matematika numerikus módszerei (Parciális differenciálegyenletek), Műszaki Könyvkiadó, Budapest, pp. 147-178. 1976. [70] N.N. Yanenko: The Method of Fractional Steps, Springer-Verlag, New York, 1979. 79

[71] A.R. Mitchell: Computational Methods in Partial Differential Equations, John Wiley, New York, 1969. [72] D.M. Young: Nagy lineáris rendszerek iterációs megoldása, Műszaki Könyvkiadó, Budapest, pp 404-447. 1979. [73] A.R. Mitchell, D.F. Griffiths: The Finite Difference Method in Partial Differential Equations, Wiley-Interscience, New York, 1980. [74] C.A.J. Fletcher: Computational Techniques for Fluid Dynamics, Vol. I. (Fundamental and General Techniques), Springer-Verlag, Berlin, pp. 251-256. 1991. [75] C. Hirsch: Numerical Computation of Internal and External Flows, Vol. I. (Fundamentals of Numerical Discretization), John Wiley & Sons, New York, pp. 437-482. 1992. [76] Stoyan Gisbert, Takó Galina: Numerikus módszerek III., ELTE TypoTEX, Budapest, pp 273-282. 1997. [77] Y. Zhuang, X.-H. Sun: A high order ADI method separable generalized Helmholtz equations, Advances in Engineering Software, Vol. 31, pp. 585-591. 2000. [78] E. Whittaker: The Calculus of Observations, Blackie & Son Ltd., London & Glasgow, 1954. [79] Carl de Boor: A Practical Guide to Splines, Applied Mathematical Sciences, Vol 27, Springer-Verlag, New York, pp. 235-239. 1978. [80] Joe F. Thompson, Z.U.A. Warsi, C. Wayne Mastin: Numerical Grid Generation (Foundations and Applications), Elsevier Science Publishing Co., Inc., New York, 1985. [81] Nyíri András: Empírikus függvények simítása, Alkalmazott Matematikai Lapok, Bp. 6, 1980. [82] Nyíri András: Vízgépek lapátozott terére vonatkozó háromdimenziós potenciálos áramlási feladat megoldása, Doktori értekezés, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, pp. 43-59, pp. 74-107. 1990. [83] A. Nyíri: Surface Fitting and a New Direct Method for Solving Block Band Linear System, Computers and Mathematics with Applications, Vol. 38, pp. 161-173. 1999. [84] C.A.J. Fletcher: Computational Techniques for Fluid Dynamics, Vol. II. (Specific Techniques for Different Flow Categories), Springer-Verlag, Berlin, pp. 81-127. 1991. [85] J. Nikuradse: Gesetzmäßigkeiten der turbulenten Strömung in glatten Rohren, Forschungsheft 356, VDI-Verlag GMBH, Berlin NW7/1932, L. 2-35. [86] J. Laufer: The Structure of Turbulence in Fully Developed Pipe Flow, NACA Report 1147, 1954. [87] Sz. Szabó, A. Juhász, A. Farkas, K. Baumel: Berendezés koaxiális csatornában kialakuló sebességprofil mérésére, GÉP, 4-5. szám, pp. 96-97. 1998. [88] T. Czibere, G. Janiga, K. Dombi: Turbulent flow through circular annuli, microCAD'1999, International Computer Science Conference, Miskolc, Section M: 4348, 1999. 80

[89] T. Czibere, G. Janiga, K. Dombi: Comparing the computed velocity distributions based on two different algebraic turbulence models with experimental results, 2nd International Conferene of Ph.D. Students, Section Proceedings: 29-38, Miskolc, 1999. [90] T. Czibere, G. Janiga, Sz. Szabó, H.J. Kecke, R. Praetor: Vergleich von mit unterschiedlichen Turbulenz-modellen berechneten Geschwindigkeitsprofilen, microCAD'1999, International Computer Science Conference, Miskolc, Section M: 7174, 1999. [91] T. Czibere, Sz. Szabó, A. Juhász, A. Farkas, G. Janiga: Measuring Turbulent Pipe Flow, microCAD'1999, International Computer Science Conference, Miskolc, Section M: 7579, 1999. [92] Sz. Szabó, A. Juhász, A. Farkas, G. Janiga: Determination the Velocity Distribution with High Accuracy, 11th Conference on Fluid and Heat Machinery and Equipment, Budapest, CD-ROM: 1-3, 1999. [93] T. Czibere, G. Janiga: Computation of Laminar-Turbulent Transition in Pipe Flow, microCAD'2000, International Computer Science Conference, Miskolc, Section M: 4550, 2000. [94] Sz. Szabó, H.J. Kecke: Untersuchung der Messungsmöglichkeiten der Geschwindigkeitsverteilung innerhalb eines Koaxialen Kanals mittels LDV-Methoden, microCAD'2000, International Computer Science Conference, Miskolc, Section M: 145-150, 2000. [95] T. Czibere, Sz. Szabó, A. Farkas, R. Praetor, G. Janiga: Experimental and Theoretical Analysis of a Coaxial Cylindrical Channel Flow, microCAD'2001, International Computer Science Conference, Miskolc, Section N: 25-32, 2001. [96] G. Janiga, Sz. Szabó, A. Farkas, H.J. Kecke, R. Praetor, B. Wunderlich: Untersuchung der turbulenten ebenen Strömung in einem gekrümmten Kanal, microCAD'2001, International Computer Science Conference, Miskolc, Section N: 33-40, 2001. [97] Sz. Szabó, H. J. Kecke: Experimentelle Bestimmung der Geschwindigkeitsverteilung in einem strömungsmaschinen-typischen Kanal mittels Laser-Doppler-Velocimetrie (LDV), Technisches Messen, Vol. 68, pp. 131-139. 2001. [98] G. Janiga, Sz. Szabó, A. Farkas, Á. Szabó: Development of a measuring device and probes for measuring turbulent velocity profiles, Proceedings of Third Conference on Mechanical Engineering, Budapest, Springer Hungarica: 1-5, 2002. [99] G. Janiga: Computation of turbulent flow in an S-shaped channel, Conference on Numerical Methods and Computational Mechanics, 2002. [100] D.C. Wilcox: Turbulence modeling for CFD, DCW Industries, Inc., La Canada, California, 2000. [101] J.H. Ferziger, M. Perić: Computational Methods for Fluid Dynamics, Springer-Verlag, Berlin, pp. 39-66. 2002. 81

AZ ÉRTEKEZÉS TÉMÁJÁBAN MEGJELENT TUDOMÁNYOS PUBLIKÁCIÓK [P.1] Determination of the Trajectories on the Meridional Section of an Impeller, microCAD'2000, International Computer Science Conference, Miskolc, pp. 101- 105. 2000. Társszerző: Nyíri András. [P.2] Ortogonális trajektóriák meghatározása egy járókerék meridián-metszetében, Doktoranduszok Fóruma, Szekciókiadvány, Miskolc, pp. 58-63. 2000. [P.3] Determination of the Trajectories on the Meridional Section of a Mixed-Flow Hydraulic Machine, microCAD'2001, International Computer Science Conference, Miskolc, pp. 59-63. 2001. [P.4] A Computational Method for Trajectories on a Meridional Section the Boundaries of Which are Given by Discrete Sets of Points, 3rd International Scientific Conference, Herl'any, Acta Mechanica Slovaca, 3/2001, pp. 127-132. 2001. Társszerző: Nyíri András. [P.5] Mesh Generation for the Meridional Section of a Mixed-Flow Hydraulic Machine with Given Discrete Points on its Boundary, 3rd International Scientific Conference of Ph.D. Students, Miskolc, pp. 259-264. 2001. [P.6] Ortogonális görbevonalú koordinátahálózat létrehozása tetszőleges perempontokkal adott meridiánsíkcsatornák esetén, Doktoranduszok Fóruma, Szekciókiadvány, Gépészmérnöki Kar, Miskolc, pp. 110-117. 2001. [P.7] Ortogonális görbevonalú koordináta-hálózat létrehozása, GÉP LIV. évfolyam, pp. 3839. 2003/1. [P.8] A New Mesh Generation Method for Meridional Section of Turbomachinery, International Conference of Water Service Science, Brno-Úbislav, pp. 49-56. 2003. [P.9] Fully-Developed Turbulent Pipe Flow Based on the Vorticity Transport, International Conference of Water Service Science, Brno-Úbislav, pp. 57-65. 2003.

82

[P.10] Comparison of the Analytical and Numerical Solution of Fully-Developed Turbulent Pipe Flow, 4th International Scientific Conference of Ph.D. Students, Miskolc, pp. 147-152. 2003. [P.11] Numerical Computation of Fully-Developed Turbulent Pipe Flow Based on the Streamfunction-Vorticity Formulation, The 12th International Conference on Fluid Flow Technologies, Conference on Modelling Fluid Flow (CMFF'03), Budapest, Hungary, pp. 638-645. 2003. [P.12] Numerical Computation of Fully-Developed Turbulent Channel Flow Based on the Streamfunction-Vorticity Formulation, 5th International Scientific Conference, Herl'any, Acta Mechanica Slovaca (ISSN 1335-2393), pp. 449-454. 3/2003.

83

AZ ÉRTEKEZÉS TÉMÁJÁVAL KAPCSOLATOS KÉZIRATOK ÉS DOKUMENTÁCIÓK [KD.1] Körkeresztmetszetű csőhöz csatlakozó koaxiális forgásfelületek által határolt tér geometriai adatrendszerének előállítása a Czibere-féle analitikus megoldással és a kidolgozott numerikus eljárással, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, Kézirat, 2003. [KD.2] Turbulens áramlás számítása nyugvó párhuzamos falak közötti résben és hosszúegyenes körkeresztmetszetű csőben a turbulens örvénydiffúzió Czibere-féle parciális differenciálegyenletének felhasználásával, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, Kézirat, 2003. [KD.3] Forgásszimmetrikus potenciáláramlás körkeresztmetszetű csőhöz csatlakozó koaxiális forgásfelületek által határolt térben, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, Kézirat, 2003. [KD.4] Forgásszimmetrikus és perdületmentes turbulens áramlás sebességterének számítása numerikus módszerrel körkeresztmetszetű csőhöz csatlakozó koaxiális forgásfelületek által határolt térben a Czibere-féle új sztochasztikus turbulencia-modell alapján, Miskolci Egyetem, Áramlás- és Hőtechnikai Gépek Tanszéke, Miskolc-Egyetemváros, Kézirat, 2003. [KD.5] Report of Researches, DAAD-MÖB Report, Universität Siegen, Institut für Fluidund Thermodynamik, 2003.

84