PhD. DISSZERTÁCIÓ
Veress Árpád
Budapest
2004
Budapesti Műszaki és Gazdaságtudományi Egyetem Repülőgépek és Hajók Tanszék Budapest, 1111, Sztoczek u. 6 J. ép. 4. em. Tel: (36)-1-463-1922, Fax: (36)-1-463-30-80, E-mail:
[email protected]
Numerikus módszerek és alkalmazások a hő- és áramlástechnikai gépekben lezajló folyamatok modellezésére PhD értekezés
Készítette: Veress Árpád okleveles gépészmérnök Témavezető: Dr. Sánta Imre
Budapest, 2004-01-06 A dolgozat bírálatai és a védéskor készült jegyzőkönyv a BME Közlekedésmérnöki Karának Dékáni Hivatalában a későbbiekben megtekinthető.
ii
Nyilatkozat
Alulírott Veress Árpád kijelentem, hogy ezt a doktori értekezést magam készítettem és abban csak a megadott forrásokat használtam fel. Minden olyan részt, amelyet szó szerint vagy azonos tartalomban, de átfogalmazva más forrásból átvettem, egyértelműen, a forrás megadásával megjelöltem.
Budapest, 2004-01-06
.........……………………… Veress Árpád
Köszönetnyilvánítás Ezúton szeretném köszönetemet kifejezni mindenek előtt szüleimnek, akik elindítottak és támogattak pályafutásomon. Szerencsésnek mondhatom magam, hogy olyan feleségem van aki megteremtette számomra a biztos hátteret a disszertáció elkészítéséhez, önzetlenül, teljes mértékben átvállalva minden, egy család életviteléhez szükséges tevékenységet különösen Bálint (1998) és Ágnes (1999) nevű gyermekeink gondozását. A von Kármán Intézetben eltöltött tanulmányaim során Van den Braembussche professzor úr témavezetésével sikerült új eredményeket elérnem a többfokozatú centrifugálkompresszorok lapátozásának tervezése során. Köszönöm türelmét és a munka tökéletesebb megvalósítása érdekében tett kritikai észrevételeit. Köszönettel tartozom a BME Gépészmérnöki Kar Áramlástan Tanszék kollektívájának, hogy szakmai tanácsaikkal, iránymutatásukkal segítették e munka elkészültét, valamint hogy a FLUENT numerikus áramlástani szoftver használatára gépidőt biztosítottak számomra. Ezúton szeretnék szintén köszönetet mondani a BME Közlekedésmérnöki Kar Repülőgépek és Hajók Tanszék kiváló munkatársainak, tanszékvezetőmnek: Dr. Rohács József professzor úrnak, aki baráti beszélgetések során bíztatott, segítette pályafutásom, témavezetőmnek: Dr. Sánta Imre úrnak, hogy hathatós szakmai segítségnyújtás mellett, a lehetőségekhez képest maximálisan tehermentesített a tanszéki feladatok alól a doktori iskola és a dolgozat elkészítésének ideje alatt. Külön köszönet illeti meg Dr. Gausz Tamás urat szakmai tanácsaiért és útmutatásaiért, ami nagyban hozzájárult e munka elkészüléséhez. Köszönöm, ……...........………………………
Veress Árpád
iii
iv
Összefoglalás Az olyan áramlástani jelenségeket leíró nem lineáris parciális differenciálegyenleteknek, mint például az Euler- vagy a Navier−Stokes-egyenleteknek, komplexitásuk miatt ez idáig nem létezik zárt alakú, általános érvényű megoldása. A számítógépek gyors fejlődésének, valamint napjaink összetett mérnöki tevékenységeivel szemben támasztott egyre magasabb szintű elvárásoknak köszönhetően azonban egyre inkább előtérbe kerülnek az áramlástan numerikus módszerei (angolul CFD (Computaional Fluid Dynamics)), amelyek hathatós segítséget nyújtanak az áramlástani-mérnöki problémák megoldásában. A disszertációban új eljárások kifejlesztésén keresztül kapcsolódhatunk bele a hő- és áramlástechnikai gépek tervezésébe saját fejlesztésű, kereskedelmi és a von Kármán kutatóintézetben (VKI) kidolgozott áramlástani szoftverek segítségével. Az új tervezési módszerek bemutatása egy konkrét berendezés kifejlesztésén keresztül történik, azonban az eljárás során megállapított irányelvek és főbb lépések általános érvényűen alkalmazhatók a hasonló típusú berendezések tervezésében. Az újfajta eljárásrendszerek kidolgozása mellett fontos szempont volt a saját numerikus programok kompaktsága és továbbfejleszthetősége a minél szélesebb áramlástani problémák megoldhatóságának érdekében. Az értekezés alapvetően két fő részre bontható. Az első a gáz halmazállapotú közegek, míg a második a folyadékok áramlásának modellezésével foglalkozik. Mindkét részben a fejezetek felépítése hasonló. A matematikai modellek és a numerikus módszerek részletes ismertetése után az eljárások validálása, érvényesítése következik. A gyakorlati szempontból fontos utolsó alfejezetekben pedig a különféle hő- és áramlástechnikai gépekben lezajlódó áramlástani folyamatok modellezése során − a numerikus analízis vagy a tervezés által nyújtott eredményekből − vonhatunk le mérnöki szempontból fontos következtetéseket a vizsgált tervezési eljárásrendszerre vonatkozóan. Az összenyomható ideális áramlás numerikus modellezésére egy újfajta konvergencia gyorsító karakterisztikus és egy forgó gépek modellezése esetén alkalmazható peremfeltétel számítási eljárás kidolgozásának segítségével került sor. A peremfeltételek adaptálása Jameson-féle mesterséges disszipációs módszer, valamint Roe által közelített Riemann (upwind) módszer alkalmazásával történt. A numerikus módszerek összehasonlítása érdekében érdemesnek tűnt kidolgozni és validálni Van Leer fluxus vektor megosztáson alapuló módszerét is. A forgógépekben lezajlódó folyamatok modellezésre kifejlesztett peremfeltétel számítási eljárás alkalmazhatóságának bemutatása egy transzszonikus áramlásba helyezett 2D-s lapátrács karakterisztikus görbéinek meghatározásán keresztül történt. Modern kereskedelmi (CFX-TASCflow) és a VKI-ban kidolgozott numerikus programok segítségével, egy új eljárásrendszer kifejlesztésének során, a többfokozatú centrifugálkompresszorok összekötő csatorna lapátozásának tervezése, optimalizálása valósítható meg. A számítógépes modellezés alapja egy 3D-s állandó lapátterhelésen alapuló lapáttervezési módszer. A VKI-s Euler alapú inverz tervező eszköz és negatív lapátelhajlítás alkalmazásával jelentős javulás érhető el a tervezési specifikációk tekintetében. A folyadékok és gázok kompaktabb, egy sémában történő modellezésének érdekében a véges térfogat módszerre átdolgozott Chorin pszeudo-kompresszibilitási módszere alkalmazható leginkább az összenyomhatatlannak feltételezett 2D-s áramlás Euleregyenleteinek numerikus megoldására. A program megvalósításához kapcsolódóan egy új abszorpciós szilárd fal és egy statikus nyomás bemeneti peremfeltétel számítási eljárás kifejlesztése fontos részét képezi a számítás gépidő csökkentésének és a környezeti beáramlások modellezésének irányába tett lépéseknek. A numerikus módszer alkalmazásával új tervezési irányelvek kidolgozása vált lehetővé a fordító típusúhoz hasonló folyadék sugárszivattyúk numerikus optimalizálásában. A számítógépes modellezések arra a következtetésre vezettek, miszerint a torokban alkalmazott letörés, az áramlásterelő fül, továbbá egy megfelelően megválasztott geometria-méret kombináció segítségével optimalizálható a hasonló típusú szivattyúk folyadékszállítása. v
vi
ABSTRACT Title of thesis: NUMERICAL METHODS AND APPLICATIONS FOR FLOW CALCULATIONS IN THE TURBO AND FLUID MACHINES Árpád VERESS, MSc, 2004 Nonlinear partial differential equations such as Euler and Navier-Stokes (NS) have no general closed form solution as yet. However, in consideration of the high level evolution of computer technology and expectations arisen from the industry, numerical methods for fluid flow (CFD (Computational Fluid Dynamics)) have come into prominence to solve fluid dynamic problems in engineering sciences. In this work a couple of new design procedures are developed in the field of turbo and fluid machinery using own codes, commercial software and numerical optimization tools developed at VKI (von Kármán Institute). New design methodologies are demonstrated through the development of a certain machine, but the governing principles and design steps are valid for all machines of similar type. Besides the development of the new design procedures, particular care has been taken of the compactness and further improvability (3D, viscous flow, etc.) of own codes. The thesis can be divided into two main parts. The first one deals with compressible and second one with incompressible flow modelling. The structure of both parts is similar. After rigorous mathematical deduction and description of the numerical methods, the next step is the validation. Finally, by the application of own codes and other CFD techniques in the design of a new product, a conclusion has been made concerning the correctness of the new design procedures. In case of compressible flow, a convergence accelerator characteristic type boundary condition (BC) and a special inlet boundary condition for rotating machines application are developed. A 2D cell centered finite volume method with a computationally effective artificial dissipation method, originally developed by Jameson [9], and Roe approximated Riemann solver [7] are reconstructed for the adaptation of the BCs. In addition a 2D Van Leer flux splitting upwind method is also developed for the validation of own codes. The application of the special inlet boundary condition for the rotating machines is tested on a 2D DCA (Double Circular Arc) cascade, which is made by the projection of the outer axissymmetric stream surface of the transonic axial rotor into 2D. With using of modern commercial CFD codes (CFX-Tascflow) and an inverse design program developed at VKI, a new design strategy is established for the deswirl vane design in the return flow channel for multistage centrifugal compressors. As a first step a constant blade loading condition has been made and an analytical procedure is developed to determine the blade camber line. As a result of the use of Euler based VKI inverse design program and introduction of negative lean, significant improvements can be observed in the design specification namely the loss coefficient and pressure recovery factor, which parameters were determined by using CFX-Tascflow NS solver. In order to be able to use the same scheme for compressible and incompressible flow modelling, Chorin’s pseudo-compressibility method is developed for the incompressible flow by means of 2D finite volume discretization based on the Euler equations. Beside the code, a soft solid wall convergence accelerator technique and static pressure special iterative inlet boundary condition are also developed. By using the own code explained in this paragraph, a new design strategy is developed for the optimization of the return type jet pumps. It has been verified by numerical simulations that the chamfered throat, flow driver lug and a special variation of the pump geometry determine a certain jet pump configuration, which belongs to the maximum flow transportation. vii
viii
Tartalomjegyzék Bevezetés ..................................................................................................................................1 Összefoglalás .............................................................................................................. 2 I. Fejezet: Bevezetés .................................................................................................................5 1.1. Áramlások matematikai modellezése ..........................................................................5 1.2. Alapegyenletek ............................................................................................................8 1.3. Numerikus módszerek az áramlástanban ....................................................................8 1.3.1. Véges differenciák módszere..............................................................................9 1.3.2. Véges elemek módszere .....................................................................................9 1.3.3. Véges térfogat módszere ....................................................................................10 II. Fejezet: Numerikus módszerek és alkalmazásuk az összenyomható áramlás modellezésére .........................................................................................................11 2.1. Alapegyenletek – Euler-egyenletek .............................................................................11 2.1.1. Euler-egyenletrendszer konzervatív formája ......................................................12 2.1.2. Euler-egyenletrendszer kvázi-lineáris alakja ......................................................12 2.1.3. Euler-egyenletrendszer karakterisztikus alakja ..................................................14 2.1.4. Navier-Stokes-egyenletrendszer konzervatív formája........................................16 2.2. Véges térfogat módszer ...............................................................................................18 2.2.1. Konvektív tagok térbeli diszkretizációja ............................................................19 2.2.1.1. Mesterséges viszkozitás módszere .........................................................19 2.2.1.2. Áramlásirányú (upwind) módszerek ......................................................21 a. Van Leer fluxus vektor megosztáson alapuló módszere (Flux Vector Splitting Method) ...............................................................................22 b. Roe által közelített Riemann eljárás (Roe’s approximate Riemann solver)..................................................................................................23 2.2.2. Magasabb rendű térbeli diszkretizáció (MUSCL (Monotone Upstream Schemes for Conservation Laws)) .....................................................................26 2.3. Negyed rendű Runge-Kutta módszer az idő iterációra................................................29 2.4. Peremfeltételek ............................................................................................................30 2.4.1. Deconinck módszere a peremfeltételek magadására ..........................................31 2.4.2. Extrapolációs módszer a peremfeltételek magadására .......................................36 2.4.3. Karakterisztikák módszere a peremfeltételek magadására.................................36 2.5. Validáció......................................................................................................................41 2.5.1. Laval-cső ............................................................................................................43 2.5.2. Körprofilt magába foglaló 2D-s csatorna ...........................................................44 2.6. 2D-s transzszonikus rotor karakterisztika meghatározása...........................................49 2.7. Tervezési eljárásrendszer kifejlesztése a többfokozatú centrifugál kompresszorok összekötő-csatorna lapátozásának tervezésére ..........................................................52 2.7.1. Bevezetés ............................................................................................................52 2.7.2. Referencia geometria..........................................................................................54 2.7.3. Állandó terhelésű lapát tervezése (ÁLT (Állandó Lapát Terhelés)) ..................54 2.7.4. Inverz lapáttervezés ............................................................................................57 2.7.5. Tervezési specifikációk ......................................................................................59 2.7.6. Lapátelhajás hatása .............................................................................................60 2.7.7. Értékelés .............................................................................................................62
ix
III. Fejezet: Numerikus módszerek és alkalmazásuk az összenyomhatatlan ideális áramlás modellezésére...................................................................................................... 69 3.1. Pszeudo-kompresszibilitási módszer .......................................................................... 70 3.1.1. Alapegyenletek, véges térfogat módszer, peremfeltételek és validálás ............. 70 3.2. Konvergencia gyorsító abszorpciós szilárd fal peremfeltétel ..................................... 75 3.2.1. Alapegyenletek................................................................................................... 75 3.2.2. A τ és a µ paraméterek meghatározása........................................................... 78 3.2.3. Konzisztencia vizsgálat...................................................................................... 79 3.2.4. Időlépések és a konvergencia lefutás vizsgálata ................................................ 80 3.3. Tervezési irányelvek kidolgozása a fordító típusúhoz hasonló folyadéksugár szivattyúk numerikus optimalizálására ..................................................................... 81 IV. Fejezet: Tézisgyűjtemény .................................................................................................. 85 4.1. Az összenyomható áramlások numerikus modellezése .............................................. 85 4.2. Az összenyomhatatlan áramlások numerikus modellezése ........................................ 87 Irodalomjegyzék....................................................................................................................... 89
x
Ábrajegyzék 1. ábra. Másodpercekénti lebegőpontos műveletek száma a megjelenési idő függvényében 1 2. ábra. Cella központú véges térfogat háló modell ...............................................................18 3. ábra. Információterjedés szub-szonikus belépő perem esetén a karakterisztikák segítségével................................................................................................................38 4. ábra. Információterjedés szub-szonikus kilépő perem esetén a karakteriszti-kák segítségével................................................................................................................39 5/a. ábra. Mach-szám eloszlás a Laval-csőben a karakterisztikus peremfeltételekkel figyelembe vett Roe által közelített másodrendű Riemann-féle módszerrel.............41 5/b. ábra.Karakterisztikus és extrapolációs peremfeltételek összehasonlítása az iteráció szám tekintetében a Laval-csőben történő áramlásra ................................................41 6. ábra. Fent: Lökéshullám mintázat DCA lapátrácsban (Schlieren felvétel [69]). Állandó Mach-szám vonalak DCA lapátrácsban; középen: centrális diszkretizáció Jameson mesterséges disszipációs módszerével, lent: Roe által közelített másodrendű Riemann-féle módszerrel, MinMod limiterrel és a karakterisztikus peremfeltételekkel .....................................................................................................42 7. ábra. 2D Laval-cső geometria a numerikus hálóval ...........................................................43 8. ábra. Jameson mesterséges disszipációjával figyelembe vett centrális diszkretizációs módszer Laval-csőben történő validálása..................................................................43 9. ábra. A karakterisztikus peremfeltételekkel figyelembe vett Roe által közelített másodrendű Riemann-féle módszer Laval-csőben történő validálása ......................44 10. ábra. Másodrendű Van Leer-féle fluxusvektor megosztáson alapuló módszer Lavalcsőben történő validálása...........................................................................................44 11. ábra. Körprofilt magába foglaló 2D-s csatorna a 142x40-s numerikus hálóval.................45 12. ábra. Állandó Mach-szám vonalak a körprofilt magába foglaló 2D-s csatornában. Folytonos vonal: centrális diszkretizáció Jameson-féle mesterséges disszipációs módszerrel, szaggatott vonal: Fluent .........................................................................45 13. ábra. Állandó Mach-szám vonalak a körprofilt magába foglaló 2D-s csatornában a 142x40-s hálón. Folytonos vonal: Roe által közelített másodrendű Riemann-féle módszer a karakterisztikus peremfeltételek alkalmazásával, szaggatott vonal: Fluent .........................................................................................................................46 14. ábra. Állandó Mach-szám vonalak a körprofilt magába foglaló 2D-s csatornában a 142x40-s hálón. Folytonos vonal: Van Leer másodrendű módszere, szaggatott vonal: Fluent ..............................................................................................................46 15. ábra. A cellánkénti Mach-szám normált eltérésének abszolút értéke: fent: centrális diszkretitáció Jameson-féle disszipációval és Fluent, középen: másodrendű Roeféle módszer és Fluent, lent Van Leer másodrendű módszere és Fluent eredményei között .....................................................................................................48 16. ábra. 2D-s kompresszor rotor karakterisztika.....................................................................50 17. ábra. 3D-s kompresszor rotor modell .................................................................................51 18. ábra. A-C többfokozatú kompresszor.................................................................................52 19. ábra. Kompresszor fokozat meridián keresztmetszete .......................................................53 20. ábra. Visszatérő csatorna lapát tervezése Rothstein módszere alapján ..............................54 21. ábra. Az analitikai lapáttervezés ellenőrző felülete ............................................................55 22. ábra. Lapátszög eloszlás a csatorna külső és belső fala mentén.........................................56 23. ábra. Kiindulási lapát Mach-szám eloszlása (+) a belső (bal oldal), illetve külső csatorna faltól 12,5%-ra eső profilkontúron ..............................................................57 24. ábra. Inverz újratervezett lapát előírt (o) és számított (+) Mach-szám eloszlása a belső (bal oldal), illetve a külső csatorna faltól 12,5%-ra eső profilkontúron ....................58 25. ábra. Eredeti (szürke árnyalatú) és az újratervezett (+) profil a belső (bal oldal), illetve a külső csatorna faltól 12,5%-ra eső lapát-szegmensben .............................................58 xi
26. ábra. Eredeti (szürke árnyalatú) és az újratervezett (-) lapátprofil..................................... 59 27. ábra. A Csatorna- (bal oldal) és a lapátörvény kialakulása [77] ........................................ 60 28. ábra. A Csatorna- (bal oldal) és a lapátörvény kialakulása................................................ 61 29. ábra. Módosult nyomás eloszlás negatív lapátelhajlítás esetén ......................................... 61 30. ábra. Lapátelhajlítás hatása a csatorna örvényre (A4. ábra alapján).................................. 62 31. ábra. Az új 3D-s lapáttervezési eljárásrendszer ................................................................. 62 A1. ábra. Sebességvektorok a lapát szívott oldalán elhelyezkedő hálófelületen a referencia geometria esetén........................................................................................................ 64 A2. ábra. Sebességvektorok a lapát szívott oldalán elhelyezkedő hálófelületen ÁLT lapáttervezés esetén........................................................................................................... 65 A3. ábra. Sebességvektorok a lapát szívott oldalán elhelyezkedő hálófelületen az inverz módszerrel újratervezett ÁLT lapát-tervezés esetén................................................. 66 A4. ábra. Sebességvektorok a lapátterjedtség közepén elhelyezkedő hálófelületen az inverz módszerrel újratervezett ÁLT lapáttervezés esetén .................................................. 67 A5. ábra. Állandó statikus nyomásgörbék a lapátterjedtség közepén elhelyezkedő hálófelületen az inverz módszerrel újratervezett ÁLT lapáttervezés esetén............. 68 32. ábra. Torlóponti (baloldal) és statikus nyomás bemeneti peremfeltételek hatása az állandó statikus nyomás görbékre ............................................................................. 72 33. ábra. Statikus és torlóponti bemeneti peremfeltételek hatása az állandó nyomás görbékre .................................................................................................................... 72 34. ábra. Állandó statikus nyomás görbék (A: saját eredmény, B: FLUENT) ........................ 73 35. ábra. Állandó statikus nyomás görbék (folytonos vonal: saját eredmény, szaggatott vonal: FLUENT) ....................................................................................................... 74 36. ábra. A cellánkénti statikus nyomás normált eltérésének abszolút értéke: saját és a Fluent program között............................................................................................... 74 37. ábra. Diszkretizációs séma az abszorpciós fal peremfeltételre .......................................... 76 38. ábra. Rugós csillapító lengőrendszer ................................................................................. 77 39. ábra. Konvergencia lefutások............................................................................................. 80 40. ábra. Peremfeltétel alkalmazások eredményeinek összehasonlítása.................................. 81 41. ábra. 1,9 mm-s fúvóka átmérőjű alap modell fordító sugárszivattyú metszeti (bal oldal) és robbantott ábrája ................................................................................................... 81 42. ábra. Áramvonalak az alap modell fordító sugárszivattyúban. A keverőtér beömlő keresztmetszete A: letörés nélkül, B: letöréssel........................................................ 82 43. ábra. Sugárszivattyú alap geometria hálózása a letöréssel+ és a 0,5-mm hosszú áramlásterelő füllel*.................................................................................................. 82 44. ábra. Sebességeloszlás [m/s] (fent) és áramvonalak (lent) a sugárszivattyúban; A: Alap geometria, B: Optimalizált geometria....................................................................... 83
Táblázatok jegyzéke 1. táblázat. Eltérés a CFD programok Mach-szám eredményei között.................................... 47 2. táblázat. A tervezési folyamat eredményeinek összehasonlítása ......................................... 60 3. táblázat. Lapátelhajlítás hatása............................................................................................. 61 4. táblázat. Eltérés a CFD programok statikus nyomás eredményei között............................. 74 5. táblázat. A különféle geometriai modellek és szállítóképességeik ...................................... 83
xii
Jelölésjegyzék Latin betűk, jelölések: • • • • • • • • • • • • • • • • • • • •
A: A: ai,j: B: bi,j: C: Cp: Cp: Cv: c: c: D: dm: ds: d(2): d(4): E: e: v ex : v ey :
F fluxus Jacobi-féle mátrixa; keresztmetszet [m2]; A mátrix (i,j) eleme; G fluxus Jacobi-féle mátrixa; Bvmátrix (i,j) eleme; H fluxusvektor Jacobi-féle mátrixa; állandó nyomáson vett fajhő [J/(kgK)]; nyomásnövekedési tényező; állandó térfogaton vett fajhő [J/(kgK)]; hangsebesség [m/s]; húrhossz [m]; mesterséges disszipáció; meridionális irányú elemi távolság; lapáthossz irányú elemi távolság; másodrendű mesterséges disszipációs tag; negyedrendű mesterséges disszipációs tag; torlóponti energia [J/kg]; statikus energia [J/kg]; x irányú egységvektor; y irányú egységvektor;
• • • • • • • • • • • • • • • • • •
F: Fv: G: Gv: hv: H: ~ H: k: L: M: M: m& : m: m: Np: Nf: n: v n:
konvektív fluxusvektor (x komponens); diffúzív fluxusvektor (x komponens); konvektív fluxusvektor (y komponens); diffúzív fluxusvektor (y komponens); entalpia [J/kg]; Euler-egyenlet teljes fluxusvektora; numerikus fluxusfüggvény; merevségi tényező; C mátrix bal oldali sajátvektor mátrixa; konzervatívból primitív formába való átalakítás transzformációs mátrixa; Mach-szám; tömegáram [kg/s]; meridionális irányú vektor; meridionális hossz [m], meridionális irány; a számítási tér összes pontjainak száma [db]; az ellenőrző térfogatot határoló oldalak száma [db]; normális irány; kifelé mutató normál egységvektor ( n x , n y ) komponensekkel;
• • • • • •
P: P: p: pto: R: R:
primitív változók vektora; = p ρ , a mesterséges kompresszibilitási módszer paramétere [J/kg]; statikus nyomás [Pa]; torlóponti nyomás [Pa]; sugár [m] (sugár-irány); C mátrix jobb oldali sajátvektor mátrixa; xiii
• • • • • • • • • • • • • • • • • • • • • • •
R: Rc: ℜ: Re: r: S: s: v s: s: T: Tto: t: U: uv: V: V& : V: v: W: W: w: x,y,z: z:
specifikus gázállandó (levegőre: 287,2 [J/(kgK]); görbületi sugár [m]; reziduum; Reynolds-szám; sugár irány [mm]; amplitúdó erősítés mérőszáma; görbe vonal hossza [m]; v n -ra merőleges, vele jobbrendszert alkotó egységvektor; tangenciális irányú vektor; statikus hőmérséklet [K]; torlóponti hőmérséklet [K]; idő [s]; konzervatív változók vektora; x irányú sebességkomponens [m/s]; sebességvektor (u, v, w) komponensekkel [m/s]; térfogatáram [l/h]; abszolút sebesség nagysága [m/s]; y irányú sebességkomponens [m/s]; karakterisztikus változók vektora; relatív sebesség [m/s]; z irányú sebességkomponens [m/s]; Descartes-féle térbeli változók; lapátszám [db];
Görög betűk, jelölések: • • • • • • • • • • • • • • • • • • • • • • xiv
αk : α: α: β: β: β: ∆: ε: ε: Γ : γ: κ: λ: Λ: µ: µ: π: θ: ρ: τ: τ i,j : Ω:
Runge-Kutta együttható; sebességvektor iránya [°]; csillapítási tényező; kompressziós paraméter; áramlás- és lapátszög [°] (sugártól mérve); mesterséges kompresszibilitási tényező; két állapot közötti eltérés; stabilizáló tag hangoló paramétere; szállítóképesség; Ω ellenőrzőfelületet határoló oldalhossz [m]; fajhőviszony; rekonstrukciós paraméter; C Jacobi-féle mátrix sajátértéke; sajátérték mátrix; molekuláris viszkozitás [m2/s]; súlyozó paraméter; nyomásviszony; kerületi irányú lapátprofil-helyzetszög (adott sugártól mérve) [°]; sűrűség (kg/m3); súlyozó-paraméter; viszkózus feszültség [Pa]; ellenőrző felület [m2];
ω: ω:
frekvencia [1/s]; veszteségi tényező;
• • • • • • • • •
: ∗: ∼: ˆ: +: -: á: B: B:
• • • • • • • • • • • •
be: I: i,j: in: ki: l: L: n: n: out: össz: R:
• •
R: S:
• • • • •
s: st: t: táp: to:
statikus állapotjelzők; peremfeltételként adott; primitív alak; Roe-féle paramétervektor-elem; pozitív (sajátérték); negatív (sajátérték); áramlás; peremen érvényes paraméter (Boundary); az adott időpillanathoz tartozó, a peremen érvényes paraméterek (mesterséges kompresszibilitási módszer); belépés; belső számítási térben érvényes paraméter (Interior); térbeli indexek; bemenet (input); kilépés; lapát; bal (Left) oldal a Riemann probléma esetén; v n irányában; időbeli index; kimenet (output); össz-folyadék; előző időlépéshez tartozó, a peremen érvényes paraméterek (mesterséges kompresszibilitási módszer); jobb (Right) oldal a Riemann probléma esetén; a következő időlépéshez tartozó, a peremen érvényes paraméter (mesterséges kompresszibilitási módszer); v s irányában; statikus állapotjelzők; tangenciális; táp-oldali; torlóponti;
• • Indexek:
Rövidítések: • • • • • • • • •
BÉ: belépő él; Ct: állandó (Constant); KÉ: kilépő él; k: kompresszor fokozat; LE, le: belépőél (leading edge); lv: lapátvastagság; nyo: nyomott oldal; szo: szívott oldal; TE, te: kilépőél (trailing edge);
xv
xvi
Bevezetés A mai mérnöki tervezés tárgyai rendkívül összetettek, drágák és a velük szemben támasztott biztonsági és tartóssági követelmények magasak. Az olyan modern technológiák, mint például a járműgyártás-gépgyártás, az űrhajózás, vagy a mikroelektronikai és mechanikai rendszerek gyártása megkövetelik az eddigi mérnöki analitika módszereinek gyökeres változtatását és megfelelő továbbfejlesztését. Újfajta, speciális igények jelennek meg napról napra az ipar minden területén, keresve a választ eddig nem ismert problémákra. A jelenségek megértése és behatóbb vizsgálata olyan Flops (Floating point operations per second) 14 ismert és új matematikai 10 Másodpercenkénti lebegõpontos mûveletek száma modellek felírását, illetve 13 felállítását teszi szükségessé, 10 A világ leggyorsabb szuperszámítógépe amelyekkel a bonyolult terSzuperszámítógép Mikroprocesszor (párhuzamos) mészeti-fizikai rendszerek 12 10 IBM SP Intel x86 sokkal részletesebben sziAlpha INTEL PARAGON CM-5 mulálhatók és elemezhetők. A 1011 modellekben előforduló CM-2 (nemlineáris) parciális diffe1010 renciál-egyenletek (szilárdtestfizika, áramlástan, terCRAY-C90 109 modinamika) leírják az adott NEC SX-3 CRAY X-MP CRAY 2 természeti folyamatokat. A cél CRAY 1 108 CYBER 205 ezen egyenletek minél CRAY Y-MP IBM 3090 pontosabb és átfogóbb 107 megfogalmazása, adaptálása és CDC 7600 megoldása. CDC 6600 106 IBM STRETCH Ehhez nyújtanak hathatós segítséget a numerikus mód105 szerek – elsősorban a zárt IBM 7090 alakú általános megoldás 104 hiánya esetén –, amelyek annál IBM 704 IBM 701 általánosabban és intenzí103 vebben alkalmazhatók, minél UNIVAC Megjelenés éve nagyobb a számítógépek kapacitása. Mivel azonban 1950 1960 1970 1980 1990 2000 2010 ezen a területen a fejlődés 1. ábra. Másodpercekénti lebegőpontos műveletek rendkívül erős (1. ábra), ezért a megoldandó feladatok is egyre száma a megjelenési idő függvényében növekvő komplexitással fogalmazódhatnak meg. A rendelkezésre álló kereskedelmi szoftverek akkor érthetők meg igazán, akkor alkalmazhatók tudatosan és átfogóan, ha a programok szelleme teljes mértékben ismert. Ez, pedig csak akkor lehetséges, ha minden lépésnél tudatában van a felhasználó annak, hogy mit tesz, és hogy ez a tevékenység milyen kihatással van a számítás oldalán [20]. A különféle áramlások modellezése kitüntetett figyelmet érdemel, nem csak azért, mert igen összetett folyamatokról van szó, hanem legfőképpen azért, mivel az élet minden területén találkozunk vele (pl. közlekedés, ipar, csillagászat, gyógyászat, sport). A különféle numerikus áramlástani alkalmazások (pl.: [37-39,24]) napjaink korszerű ipari felhasználása szempontjából elsődleges fontosságúak. Ebben a munkában a két legáltalánosabb áramláson: a gázok és a folyadékok mechanikáján keresztül kapcsolódhatunk bele e folyamatok matematikai modellezésébe. 1
Összefoglalás A disszertáció célja új tervezési irányelvek kidolgozása a hő- és áramlástechnikai gépek területén saját fejlesztésű, kereskedelmi és a VKI (von Kármán Institute) kutatóintézetben kidolgozott áramlástani szoftverek segítségével. Az újfajta tervezési eljárásrendszerek kifejlesztése mellett nagy figyelmet fordítottam az általam megírt numerikus programok kompaktságára és továbbfejleszthetőségére a minél szélesebb körű áramlástani probléma megoldhatósága érdekében. Az első, az ismereteket megalapozó bevezető fejezet után a disszertáció két fő részre tagolódik. Az első rész a gáz halmazállapotú, míg a második a folyadékok áramlásának modellezésével foglalkozik. Mindkét részben a fejezetek felépítése hasonló. A matematikai modellek és a numerikus módszerek részletes ismertetése után az eljárások validálása következik áramlástani iskolapéldákon keresztül bemutatva. A gyakorlati szempontból fontos utolsó alfejezetekben, pedig a különféle hő- és áramlástechnikai gépekben lezajlódó áramlástani folyamatok modellezése során – a numerikus analízis vagy a tervezés által nyújtott eredményekből – vontam le mérnöki szempontból fontos következtetéseket a vizsgált tervezési eljárásrendszer hatékonyságáról. Az összenyomható áramlás numerikus modellezése során kidolgoztam egy újfajta karakterisztikus peremfeltétel számítási eljárást, amelyben a numerikus peremfeltételeket a cella normális irányú Mach-szám nemlineáris egyenletének segítségével számítottam ki (2.4.3. alfejezet). A módszer alkalmazásával kb. 20 %-kal csökkent az iteráció szám és 49,22 %-kal csökkent (a 2.4.3. alfejezetben megadott geometriai és az áramlástani perem, illetve kezdeti feltételek esetén) az állandósult állapot eléréséhez szükséges idő a hagyományos extrapolációs peremfeltételhez képest. Az eljárás adaptálása érdekében számítógép által kezelhető alakra hoztam Roe által közelített Riemann eljárás 2D-re kidolgozott változatát [55], illetve szintén a modern, nagy pontosságú upwind (high resolution upwind methods) módszer megvalósítása érdekében a hasonlóan széles körben elterjedt van Leer fluxus vektor megosztáson alapuló módszert (2.2.1.2. alfejezet) [76]. Szintén az összenyomható áramlás numerikus modellezése során kidolgoztam egy olyan, az axiális kompresszor rotorok matematikai modellezésébe integrálható bemeneti peremfeltétel számítási eljárást, amely figyelembe veszi a rotor fordulatszámát és a futólapátozás mögött kialakult nyomásnövekedést (2.6. alfejezet). A két paraméter határozza meg a belépő tömegáramot, illetve a belépő áramlás irányát. Az egyre inkább elterjedésben lévő transzszonikus fokozatok legjobb numerikus szimulációja érdekében (egyszerű programozhatóság, alacsony gépidő, illetve kis számítógépi kapacitás melletti megfelelő pontosság), valamint az újfajta peremfeltétel alkalmazására egy 2D-s, Euler-alapú, cella központú véges térfogat módszert dolgoztam ki. Az eljárás időben állandósult ideális folyamatok modellezésére alkalmas, amelyben egy negyedrendű Runge-Kutta módszer biztosítja (2.3. alfejezet) – egy tetszőleges kezdeti állapotból kiindulva – az állandósult állapot elérését. A centrális diszkretizáció okozta páros-páratlan oszcilláció elkerülését, illetve a szakadási felület közelében jelentkező pontatlanságok megszüntetését egy beépített numerikus szenzor segítségével Jameson által kidolgozott mesterséges disszipáció biztosítja (2.2.1.1. alfejezet) [35]. A módszer segítségével meghatároztam egy 2D-s DCA ((Double Circular Arc) kettős körívből álló profil [43]) profilú előperdítés nélküli rotor metszet karakterisztikáját. A numerikus analízis során a 2D-s lapátrácsot a lapátozás külső felületén elhelyezkedő áramfelület síkba vetítésével állítottam elő. Modern kereskedelmi és a VKI kutatóintézetben kidolgozott numerikus-tervezési és optimalizálási programok segítségével új tervezési eljárásrendszert fejlesztettem ki a többfokozatú centrifugálkompresszorok összekötő csatorna lapátozásának tervezésére (2.7. alfejezet). A számítógépes modellezés során először egy állandó lapátterhelésen alapuló, lapáttervezési módszert dolgoztam ki (2.7.3. alfejezet). Az analitikus eljárás során zérusértékű cirkuláció feltételezéséből indultam ki, amelyet a szívott és a nyomott oldal, valamint két állandó sugarú kör által közrezárt ellenőrzőfelületre írtam fel. Az állandó lapátterhelés 2
feltételét a lapát két oldalán kialakult sebességkülönbség állandóvá tételével értem el. Ezek után numerikus-integrálással meghatároztam a lapátszöget a csatorna szívott és nyomott oldalán. A lapátvastagság eloszlást a vázvonal első 75%-ban egy ellipszis kis tengelyével párhuzamos metszékeivel modelleztem. A lapát másik felében lineáris eloszlást alkalmaztam, a maximális vastagságtól az előírt kilépő él vastagságáig. Numerikus áramlástani vizsgálatok arra engedtek következtetni, hogy jelentős javulás érhető el a veszteségi és a nyomásnövekedési tényezők tekintetében ( ω , C p ), ha a lapátok belépő élét kiterjesztem a lapát nélküli diffúzor kimenetéig. Ebben az esetben megszűnt a fordítókönyökbe történő belépéskor kialakult leválási buborék. A numerikus analízist minden esetben a CFXTASCflow Navier-Stokes alapú kereskedelmi szoftver segítségével végezetem, amely az áramlástani gépekben lezajlódó folyamatok modellezésének napjaink élenjáró eszköze. A további optimalizáció érdekében egy 3D-s inverz lapáttervezési eljárás segítségével tovább tudtam javítani a tervezési paramétereken (2.7.4. alfejezet). A szoftver a VKI kutatóintézetben kifejlesztett Euler-alapú numerikus tervező-eszköz [18,19]. A program az általam meghatározott Mach-szám eloszlás alapján úgy módosította a lapátprofilt, hogy sikerült megszüntetni a visszatérő ágban a lapát szívott oldalán idáig minden esetben előálló leválást. A további veszteségcsökkentés érdekében megvizsgáltam a lapátelhajlításnak a szekunder áramlás okozta lapát- és csatorna-örvény intenzitás csökkentésére gyakorolt hatását. Ennek eredményeként további javulást sikerült elérnem negatív lapátelhajlítás esetén a tervezési specifikációk tekintetében ( ω , C p ) (2.7.6. alfejezet). Az új tervezési eljárásrendszer kidolgozását egy konkrét berendezés kifejlesztésén keresztül mutatom be, azonban az eljárás során megállapított irányelvek és főbb lépések általános érvényűen alkalmazhatók a hasonló típusú összekötő csatornák lapátozásának tervezésében. A folyadékok és gázok ugyanazon sémában történő modellezése érdekében átdolgoztam Chorin [11] pszeudo-kompresszibilitási módszerét a véges térfogat diszkretizációs eljárásra az összenyomhatatlannak feltételezett áramlás Euler-egyenleteinek numerikus megoldására (3.1. alfejezet). Az egyenleteknek ebben a formában való felírása formailag megegyezik az összenyomható közegek alapegyenleteivel (2.1.1. alfejezet), így a megoldás menetének hasonlósága miatt ugyanaz az eljárásrendszer alkalmazható mindkét esetben, ami általánosabb értelemben használható numerikus módszert jelent a különféle szimulációs alkalmazásokban. Az összenyomható áramlás modellezéséhez hasonlóan (2.2. alfejezet) a diszkretizálást szintén a kétdimenziós cellaközpontú véges térfogat módszer segítségével végeztem el, amelyben az időben beállt folyamat elérésére a kis memóriaigényű negyedrendű Runge-Kutta módszert alkalmaztam. A numerikus séma stabilitását Jameson által bevezetett mesterséges viszkozitás segítségével értem el [35]. A módszer nagy előnye, hogy könnyen kibővíthető a harmadik térbeli koordinátával, illetve kiterjeszthető a diffuzív fluxusok bevezetésével a valóságos áramlás, kettős időléptetés (dual time stepping) alkalmazásával pedig a tranziens folyamatok modellezésére. A numerikus eljáráshoz kapcsolódóan új abszorpciós szilárd fal peremfeltétel számítási eljárást dolgoztam ki (3.2. alfejezet). A módszer alkalmazásával kb. 40 %-kal csökkent az iteráció szám, mialatt az időben beállt folyamat eléréséhez szükséges idő 59,5 %-kal szintén csökkent (a 3.2.4. alfejezetben megadott geometriai és az áramlástani perem, valamint a kezdeti feltételek esetén) a hagyományos szimmetria szilárd fal peremfeltétel alkalmazásához képest. A minél átfogóbb mérnöki-áramlástani problémák modellezésére − szintén az előzőekben ismertetett numerikus eljárásban történő alkalmazásra − kidolgoztam egy új, specifikus peremfeltétel számítási eljárást, amelynek alkalmazása esetén a bemeneti torlóponti nyomás helyett a statikus nyomás megadásával lehet a környezeti beáramlásokat szimulálni (3.1.1. alfejezet). A 3.1. és 3.2. alfejezetekben az ideális összenyomhatatlan áramlás modellezésére kidolgozott numerikus eszköz segítségével új tervezési irányelveket dolgoztam ki a fordító 3
típusúhoz hasonló sugárszivattyúk numerikus optimalizálására (3.3. alfejezet). Számítógépes modellezések során bebizonyítottam, hogy a torokban alkalmazott letörés és az áramlásterelő fül alkalmazása esetén növelni lehet a szivattyú folyadékszállítását. Igazoltam továbbá, hogy a minimális gyártási módosítások mellett egyértelműen megállapítható a szivattyú beömlő keresztmetszete után érvényes olyan geometria kombináció, melynek segítségével maximalizálható a folyadékszállítás. A tervezési eljárásrendszer kidolgozása egyben a VISTEON Hungary Kft. részére végzett fejlesztői munka része [78] − amelynek célja egy új termék előállítása − volt, azonban az eljárás során megállapított irányelvek általános érvényűen alkalmazhatók minden hasonló típusú szivattyú tervezésében.
4
I. Fejezet: Bevezetés A 19. sz. eleje óta ismertek azok a fizikai törvényszerűségek, amelyek a folytonos áramlás dinamikáját írják le. Ezek a törvényszerűségek azon a megfigyelésen alapulnak, amely szerint három alapvető elv valósul meg az áramlás során. Ezek rendre a tömeg-, az impulzus és az energia megmaradás, amelyek parciális differenciálegyenletek (PDE) formájában öltenek matematikai alakot. A Navier-Stokes (NS) impulzus-megmaradási egyenleteket M. Navier és G. Stokes tudósokról nevezték el, akik egymástól függetlenül, elsőként írták fel őket. Az egyenletek teljes mértékben leírják a folytonos áramlástani jelenségeket, függetlenül attól, hogy milyen egyéb külső hatás éri a vizsgált rendszert. A NS-egyenletrendszer1 nemlineáris, másodrendű, hibrid (parabolikus és hiperbolikus) típusú, napjainkig nem létezik általános érvényű zárt alakú megoldása. A számítógépek gyors fejlődésének köszönhetően (1. ábra) azonban, egyre inkább előtérbe kerülnek az áramlástan numerikus módszerei, amelyek segítségével számokkal helyettesíthetők a kiinduló egyenletek állapotváltozói és egy megfelelően választott eljárás segítségével meghatározhatók a keresett áramlástani paraméterek. A különféle numerikus módszerek alkalmazásakor igen fontos, hogy tisztában legyünk az alkalmazott matematikai modell tulajdonságaival, az általa szolgáltatott megoldások helyességével, vagyis minden esetben ismerni kell a vonatkozó fizikai rendszer által elfogadható eredményeket. Az alkalmazott matematikai és fizikai ismereteken kívül a megfelelő programozási nyelv magas szintű alkalmazása jelentősen hozzásegít hatékonyabb numerikus módszerek kifejlesztéséhez, amelyeknél a cél mindig a lehető legpontosabb megoldás a lehető legrövidebb gépidő alatt.
1.1. Áramlások matematikai modellezése A disszertációban a kontinuum-mechanika alapján értelmezett matematikai modelleket használtam, azonban érdemes megemlíteni napjaink egyre inkább előtérbe kerülő olyan áramlás-modellezési problémáit, amelyekben a lokális Knudsen szám 0,1 fölé emelkedhet. Ilyen például a MEMS (Mikro-Electronical and Mechanical Systems) rendszerekben létrejövő, vagy a föld körüli orbitális pályán az űrjárművek körül kialakult hiperszonikus áramlás. Ezekben az esetekben a kontinuum-mechanika alapján felírt NS- ( Kn = 0 ,0 ÷ 0 ,1 ) vagy Euler-egyenletek ( Kn = 0,0 ) elveszítik érvényességüket és a molekuláris kapcsolatokat figyelembevevő Boltzmann-egyenletek, vagy pl. a magasabb rendű Chapman-Enskog kiterjesztéssel előállítható Burnett- vagy szuper−Burnett-egyenletek segítségével modellezhető az áramlás [80]. Visszatérve a kontinuum-mechanika elvén történő közelítésre a különféle áramlások modellezése mind más és más követelményt támaszt a matematikai modellekkel szemben. Lamináris áramlás esetén nincs különösebb nehézség a NS-egyenletek numerikus megoldásakor, mivel a modell nem tartalmaz empirikus kifejezéseket. Sajnos azonban, a lamináris áramlás igen kis területét fedi le a különféle mérnöki-áramlástani problémáknak. A turbulens áramlások esetén a NS-egyenletek közvetlen numerikus integrálása szintén elvégezhető, de fontos a megfelelően kicsi tér- és időlépték. Az úgynevezett DNS (Direct Numerical Simulation) (Orszag, Patterson [51], Rogallo [56], Spalart [67] és mások) módszerben a pillanatnyi változók meghatározásán keresztül kaphatunk kielégítően pontos eredményt a NS-egyenletrendszerre. A megfelelően kicsiny tér és időlépték feletti megoldás
1
A modern irodalomban a NS-, illetve az Euler-egyenletek alatt nem csak az impulzus-megmaradás egyenletei értendőek, hanem beletartoznak a tömeg- és az energia-megmaradás egyenletei is.
5
óriási számítógépi kapacitást igényel, még az olyan egyszerű teszteseteknél is, mint például az egyenes csatornában történő áramlás. Ekkor, a háló csomópontjainak becsült száma megközelítőleg, N DNS ≅ 0,088 Reh9 / 4 , amelyben Re h az áramlás középsebességén és csatorna magasságán alapuló Reynolds-szám. A mai szuperszámítógépek teljesítményének tekintetében, a hálópontok számának előzőekben történő meghatározása miatt csak pár tízezerre korlátozódik az a Reynolds-szám, amellyel még numerikus szempontból gazdaságosan szimulálható az áramlás [79]. A gépidő és a memória csökkentését szem előtt tartva, nehéz feladat a különféle turbulens intenzitással rendelkező áramlás paramétereit minden léptékben pontosan meghatározni, ezért érdemesnek tűnik csak a nagyobb léptékű örvényeket figyelembe venni, amelyek a legtöbb energiát tartalmazzák. Az áramlásban kialakuló kisebb örvények modellezésére pedig, egy megfelelően választott pl. al-háló lépték (Sub-Grid Scale (SGS)) modell alkalmazható. Ezen az elméleten alapszik a nagy örvények szimulációja (Large Eddy Simulation (LES)) módszer (Smagorinsky [66], Germano [25], Zang [81], Shen és Yue [64] és mások). A háló léptékű áramlás (örvények) meghatározása a mozgásegyenletek paramétereinek térbeli szűrésének segítségével, míg az al-háló léptékkel jellemzett áramlástani paraméterek meghatározása Boussinesq [6] közelítésén alapuló turbulens viszkozitás segítségével történik. Az eljárás azon az egyszerű megfigyelésen alapszik, miszerint a nagyobb örvények áramlásfüggők, míg a kisebbek jóval inkább univerzálisak, függetlenek a geometriától. A LES módszer a korszerű CFD (Computational Fluid Dynamics (számítógépes áramlásmodellezés)) egyik legígéretesebben fejlődő módszere, azonban a DNS módszerhez hasonlóan, a nagyobb Reynolds-számok estében nem alkalmazható számítástechnikai szempontból gazdaságosan, így az ipari alkalmazása ( Re jellemző = Ο ( 10 6 − 10 9 ) ) még szintén várat magára [44].
Napjaink gyakorlati-mérnöki alkalmazásaiban leginkább használható, legpontosabb közelítést adó matematikai modellje a Reynolds átlagolt Navier−Stokes-egyenletek (RANS (Reynolds Averaged Navier−Stokes equations)) [23]. A NS-egyenletek áramlástani paraméterei egy jellemző idő alatt, azzal a céllal átlagoltak, hogy megszüntessék a turbulens fluktuációból adódó bizonytalanságot, megőrizve magára a turbulenciára jellemző időléptéken kívüli időfüggő folyamatok leírását. Az így előálló új egyenletrendszer formálisan megegyezik a lamináris áramlás NS-egyenleteivel, azonban néhány új tag megjelenik benne. Ezek a tagok a Reynolds-feszültségek, amelyeket az alapegyenletek nemlineáris tagjai eredményeznek. Mivel nem ismert a kapcsolat az egyenletbeli átlagolt paraméterek, illetve a Reynolds-feszültségek között, ezért a különféle turbulencia modellek egyik fő célja, hogy határozottá tegyék a RANS egyenletet. Az utóbbi harminc évben a turbulencia modellek széles skáláját dolgozták ki. Boussinesq közelítését figyelembe véve (ami szoros kapcsolatot feltételez a lamináris és a turbulens feszültségek között) a Reynolds-feszültségek a turbulens viszkozitás segítségével fejezhetők ki. A turbulens áramlásból adódó viszkozitás a molekuláris viszkozitás matematikai modelljének segítségével állítható elő, az eredő viszkozitás pedig, e kettő összegeként írható fel. A molekuláris viszkozitás magának az áramló közegnek az állandó tulajdonsága, ellenben a turbulens viszkozitással, ami az áramlástani paraméterek függvényeként, térben és időben változhat. A turbulencia modellek célja tehát, hogy függvénykapcsolatot teremtsenek az áramlástani paraméterek és a turbulens viszkozitás között. Prandtl, 1925-ben vezette be az első modellt, amely a keveredési úthossz elmélete szerint biztosítja a kapcsolatot a turbulens nyírófeszültség és a sebesség-gradiens között. Kifinomultabb turbulencia modellek esetén új egyenletek bevezetésével biztosítható a RANS egyenletrendszer egyértelműsége. Ilyen például a turbulens mozgási energiára, illetve ennek a disszipációs rátájára felírt transzport egyenlet, amelyek már egy magasabb rendjét képviselik az egyenletszámmal definiált turbulencia modelleknek. Ennek értelmében megkülönböztethetők az egy egyenletes (pl. Prandtl [52], Baldvin-Barth [2] és SpalartAlmaras [68]), a két egyenletes (pl. k-ω, k-є [79]), illetve a magasabb rendű, mint például a Reynolds-feszültség modellek. Visszatérve Prandtl keveredési úthossz modelljére, ami pl. a 6
Cebeci-Smith [8], illetve a Baldvin-Lomax [3] modellekkel együtt képezi a turbulencia modellek legegyszerűbb kategóriáját az algebrai vagy zéró egyenletes modellek csoportját. Ezek a módszerek nem tartalmaznak explicit hozzáadott egyenletet a RANS egyenletrendszer megoldhatóságának érdekében. Végezetül fontos megemlíteni az idáig egyik csoportba sem tartozó, a turbulencia-modellezés ígéretesen fejlődő ágát, az un. sztochasztikus turbulenciamodelleket [14,36]. Az előzőekből következően tehát, a turbulens áramlás modellezésekor döntő szerepe van azoknak a függvényeknek, amelyek határozottá teszik a RANS egyenletrendszert. Elsődleges fontosságú a kapcsolatteremtő egyenletek fizikai jelentésének megismerése, mivel döntő jelentőségűek a turbulens áramlás modellezésében. Sajnos, napjainkig ez a tudás behatárolt. Nagy mennyiségű empirikus és tapasztalati összefüggés használata jellemzi a turbulens áramlás modellezését, aminek köszönhetően jelentős eltérés is adódhat a különféle modellek eredményei és a valóságos áramlás paraméterei között. Ez a pontatlanság különösen igaz az átváltási szám körüli szimulációk esetében. Napjaink komplex mérnöki problémáinak megoldása felveti az olyan újfajta modellek iránti igényt, amelyek a legszélesebb áramlástani jelenségek leírását képesek megvalósítani [44]. Az áramló közegek dinamikájának következő nagyobb csoportja a leválás nélküli nagyobb Reynolds-számú áramlás. A kísérleti tapasztalatok azt mutatják, hogy a lamináris és a turbulens feszültségek hatása csak a szilárd falhoz közeli rétegben jelentős, míg ezen a vékony rétegen kívül a nem viszkózus, un. konvektív hatások dominálnak. A numerikus szimulációk tekintetében, a viszkózus és a nem viszkózus áramlást érdemes külön számolni, és egy rekurzív formula segítségével összehangolni az eredményeket. A különféle határréteg elméletek területén napjainkig elért magas szintű fizikai és matematikai ismeretek következtében sokféle, pontos számítási eljárást dolgoztak ki [7,63 és mások]. Az áramlás határrétegen kívüli tartományának modellezésében, egészen az 1970-80-as évek végéig, a potenciálos áramlásmodellezés [31,4,60,72,34,48,49,50 és mások] játszotta a főszerepet. Az izentrópikus feltétel következtében az ilyen modellek érvényessége az összenyomható áramlás esetén a transzszonikus sebesség határától megkérdőjeleződik. A tisztán potenciálos áramlás kielégíti ugyan a tömeg- az impulzus- és az energiamegmaradás egyenleteit, de a szakadás következtében előálló ugrás kapcsolatot (pl. Rankine−Hugonoit-egyenlet) nem. A modell gyenge lökéshullámok esetén még használható, mivel a szakadáson keresztül viszonylag kis entrópianövekedés lép fel, vagyis a lökéshullám által okozott irreverzibiltás jó közelítéssel elhanyagolható. A potenciálos modellek érvényességének egy másik, a hangsebesség feletti tartományban fellépő korlátja a nem egyértelmű megoldások megjelenése. A módszerekkel kapcsolatos problémákról bővebb információ Salas [60,61] munkáiban található. A potenciálos áramlásmodelleknek a különféle szakadások (lökéshullámok, kontaktszakadások és az örvényfelületek) megjelenésével fellépő hátrányain az Euler-egyenletek segítségével lehet javítani2. Az Euler-egyenletek a NS-egyenletekből származtathatók a viszkózus és a hővezetési tagok elhanyagolásával, és képezik a nem viszkózus áramlás legmagasabb fokú approximációját. Az Euler-egyenletek elsősorban a határrétegen kívüli és a leválás nélküli nagy Reynolds-számú áramlás modellezésére használhatók. Mivel napjaink ipari-áramlástani folyamataiban lezajló folyamatokra a nagy Reynolds-szám jellemző, ezért ebből a szempontból az Euler-egyenletek tűnnek a legalkalmasabb közelítésnek. A megmaradási törvények matematikai tulajdonságai jelentősen megváltoznak a viszkózus és a hővezetési tagok elhagyásával. A parciális differenciálegyenletek rendje kettőről egyre csökken, így a peremfeltételek száma szintén csökken. A nemlineáris konvektív tagok tekintetében nem történik változás. Az előálló megoldások halmaza egy szélesebb skálán értelmezett, ami egyenes következménye a diffuzív tagok elhagyásának, vagyis a szakadásos (a disztribúció elmélet szerinti gyenge) megoldások megjelenésének. Ennek 2
Az “Euler-egyenletek” kifejezés minden esetben a kiindulási egyenletekre utal.
7
köszönhetően az Euler-egyenletek bizonyos értelemben nehezebben kezelhetők, mint a NSegyenletek [44]. A különféle numerikus módszerek tekintetében a viszkózus és a nem viszkózus áramlásmodellek esetén az elvégzendő műveletek száma ez utóbbi esetben a turbulencia modellek hiánya miatt kisebb, így a számításhoz szükséges gépidő is rövidebb. Másrészt, a legtöbb komplex 3D-s áramlás, a megfelelő szimmetria-feltételek, illetve közelítések figyelembevételével, visszavezethető 2D-s folyamatokra. Az Euler-egyenletek esetén pl. a határréteg és ezáltal a szekunder-áramlások modellezésének hiánya miatt a 2, illetve a 3D-s eredmények közötti eltérés kisebb, mint a viszkózus áramlás-modellezések esetén. Mindezek figyelembevételével, − elsősorban az előtervezési folyamatok lerövidítésének érdekében − a kisebb gépidő-igény jól kompenzálhatja a lehető legnagyobb approximációs fokú 2D-s ideális áramlástani modellekben rejlő közelítést. Az utóbbi negyven évben nagy számú numerikus módszer kifejlesztésére került sor, amelyekben különféle közelítéseket alkalmaztak a nemlineáris megmaradási törvények közelítésére, azonban a fejlődés, főleg a többdimenziós problémákra, még nem mondható lezártnak. Minden módszernek megvan az előnye, és a hátránya is. A modern áramlástan numerikus módszerei között nem található olyan egységes modell, amely egyaránt alkalmazható lehetne minden áramlástani problémára.
1.2. Alapegyenletek Mint ismeretes az Euler-egyenletek integrálása után az egyenletrendszer megoldásai között előfordulhatnak különféle szakadásos megoldások: örvényfelületek (slip lines vagy vortex sheets), kontakt diszkontinuitások (contact discontinuities) és lökéshullámok (shock waves) [28], ezért a disszertációban szereplő, általam kidogozott numerikus módszerek a nem viszkózus áramlás legmagasabb fokú approximációjával, az Euler-egyenletek segítségével modellezik az áramlást. Olyan egységes numerikus modell megalkotását tűztem ki célul, amely egyaránt alkalmas az összenyomható és az összenyomhatatlan áramlások 2D-s modellezésére (2.2.1.1. és 3.1. alfejezetek). A továbbfejleszthetőség szempontjából fontos szempont az alapegyenletekben rejlő olyan viszkózus és 3D-s irányba történő kiterjesztés lehetősége, amellyel pontosabb numerikus szimulációk valósíthatók meg. Az Euler-egyenletek a NS-egyenletekből származtathatók a viszkózus és a hővezetési tagok elhagyásával. Ekkor a parciális differenciálegyenletek elsőrendűvé fajulnak, de a nem lineáris konvektív tagok jelen vannak. Mint ahogyan arról már korábban szintén volt szó, ez a fajta approximáció elsősorban a leválás nélküli nagy Reynolds-számú (1/Re → 0) és konvekciós jellegű problémák esetén ad igen jó közelítést a különféle áramlásokra.
1.3. Numerikus módszerek az áramlástanban A különféle áramlástani jelenségeket leíró nem lineáris parciális differenciálegyenleteknek, mint például az Euler- vagy a Navier−Stokes-egyenleteknek, komplexitásuk miatt, napjainkig nem létezik általános érvényű zárt alakú megoldása, ezért a különféle áramlástani-mérnöki problémák megoldása során hathatós segítséget nyújtanak az áramlástan numerikus módszerei (angolul CFD (Computational Fluid Dynamics)). A különféle eljárások alkalmazásakor, a megfelelő pontosság biztosítása mellett, az alapegyenletekben szereplő áramlástani paraméterek a zárt alakú folytonos megoldás helyett, diszkrét értékekkel helyettesíthetők. Ezért, először a folytonos számítási tartomány véges számú diszkrét felosztását, a hálózást kell elvégezni. A háló pontjainak segítségével kerülnek meghatározásra az áramlástani paraméterek. A differenciál vagy integrál formában rendelkezésre álló alapegyenletek a diszkretizálást követően kerülnek a számítógép által jól 8
kezelhető alakra. A lineáris vagy nem lineáris algebrai egyenletek azonban a kiinduló egyenletek meghatározott fokú approximációját is jelentik. A folytonosról diszkrét alakra történő átalakítás során, általában, egy adott pont körüli Taylor-sorba fejtés maradék tagjának segítségével állapítható meg a közelítés rendje. A hálóméret nullához tartása esetén a numerikus diszkretizáció hibájának is a zérushoz kell tartania, a sebessége, pedig arányos a diszkretizáció approximációjának rendjével. Az eredmények pontosságát a hálózás minősége, homogenitása is jelentős mértékben befolyásolja. A különféle diszkretizációs technikák három főbb csoportba oszthatók: a Véges Differenciák Módszere (VDM), a Véges Elemek Módszere (VEM) és a Véges Térfogat Módszere (VTM).
1.3.1. Véges differenciák módszere A VDM módszer a parciális differenciálegyenletek numerikus megoldásának egyik legrégebbi módszere. Az algebrai egyenletek a Taylor-sorba fejtés segítségével, illetve a numerikus deriválás differenciál operátorra történő közvetlen alkalmazásával állíthatók elő. A pontosság rendje (hány, illetve milyen távoli szomszédos elemeket képes az eljárása figyelembe venni egy pont számítása során, vagyis hányadik kitevőn szerepel a független változó a Taylor-sor maradék tagjában) tetszőleges lehet, de a megoldás jóságát nagymértékben befolyásolja a diszkretizáció során felvett pontok száma is. A pontosság növelésével az algebrai egyenletrendszer megoldása során előálló megnövekedett sávszélességű mátrixok invertálása jelentős akadályba ütközhet a számítógépi oldalon. A véges differenciák módszer eredményének pontossága szintén szoros összefüggésben van a numerikus háló minőségével is. Hirtelen hálóméret változás, illetve az eltorzult alakzatok szintén csökkenthetik a pontosságot [29]. A VDM módszer strukturált, uniform hálózás esetén alkalmazható jó számítógépi hatásfokkal, a számítási tér rendezett, legtöbbször négyszög (2Dben) alakú felosztása mellett, ami egyben megkönnyíti a számítások programozhatóságát. Amennyiben a háló nem homogén, illetve görbe vonalú elemekből áll, akkor célszerű mind a geometria, mind az alapegyenletek tekintetében áttérni a fizikai térből egy olyan számítási térbe, ahol teljesülnek az uniformitás feltételei. Ez szintén jelentős számítógépi kapacitást igényel. A véges differencia módszer valóságos áramlásra vonatkozó alkalmazásával kapcsolatban bővebb információ a [22]-ben található.
1.3.2. Véges elemek módszere A VEM módszer először a szilárdságtanban terjedt el, majd onnan került át az áramlástani folyamatok modellezésébe. A számítási tér cellákra vagy elemekre van felosztva, amelyek összessége alkotja a numerikus hálót. Matematikai szempontból a legelegánsabb módszer, mivel a megoldást függvények segítségével közelíti. A közelítés egy függvénytér választásával kezdődik. A függvénytérnek meg kell választani egy kvázi ortogonális alterét (bázisát) úgy, hogy ezek úgynevezett véges tartójú függvények legyenek, vagyis a tartomány korlátos része felett lehetnek zérustól különböző értékűek. Ezeket a függvényeket definiálni kell, ezért meg kell határozni azt a tartományt, ami felett a bázis függvény nem zérus. Végeredményben ez a tartomány-meghatározás jelenti az áramlási tér diszkretizálását. A háló lehet strukturált vagy nem strukturált. Ez utóbbi nagyobb lehetőséget biztosít a geometria leírására, könnyen lehet a hálózást automatizálni. A háló cellái nem lapolhatják át egymást, és ki kell tölteniük a rendelkezésre álló számítási teret. A legtöbb VEM módszer segítségével felírt anyag-megmaradási problémában kétféle fontosabb megközelítési elv létezik. Az első a variációs módszer, amelyben egy integrál (funkcionál) írható fel az egész számítási térre vonatkozóan. Ez a funkcionál tartalmazza az ismeretlen függvényeket, illetve ezek deriváltjait. A megoldás során meghatározandó az a 9
függvény, amely kielégíti a funkcionál minimumát. A második módszerben differenciálegyenletek írhatók fel egy elemi térrészre, megfelelő peremfeltételekkel, amelyek egyértelműen leírják a tartomány állapotát. A második módszer népszerűbb, de a variációs megközelítés sokkal általánosabb, egységesebb értelemben és egyszerűbb kifejezésekkel (pl. alacsonyabb rendű deriváltakkal) alkalmazható a különféle véges elemes mérnöki alkalmazásokban. Nem minden fizikai problémát lehet azonban a funkcionálanalízis segítségével leírni, így sajnos az áramlástani, pl. a hiperbolikus típusú egyenleteket sem. Mind a variációs, mind a differenciál-egyenletes megközelítési elvekben igen nehéz a funkcionál minimumát adó, illetve a differenciál-egyenletet (a hozzá tartozó peremfeltétellel együtt) kielégítő függvényt találni, ezért a megoldást csak közelíteni lehet. A közelítő módszerek közül a próba függvényeket alkalmazók a legelterjedtebbek. Az áramlástanban, a közelítési elveknek megfelelően kétféle módszer terjedt el; a variációs módszer (pl. a Rayligh-Ritz módszer), és a súlyozott hibák módszere (pl. a Galerkin (SUPG (Streamline Upwind Petrov Galerkin)), a lekisebb négyzetek és a kollokáció módszere) [53]. Ez utóbbi a differenciálegyenletes megközelítésnek felel meg, azonban a legtöbb esetben az alapegyenleteket integrál egyenletre alakítva diszkretizálják. Azokban az estekben, amikor a nem folytonos megoldások jelennek meg az áramlásban (pl. lökéshullámok, kontaktfelületek), kizárólag a megmaradási törvények integrál alakjának van fizikai jelentése, mivel a szakadásos áramlástani paraméterek deriváltjai nincsenek értelmezve. A VEM módszerek fejlődése nem tekinthető lezártnak, sőt, a bennük rejlő matematikai egzaktságnak, eleganciának, illetve a geometriai flexibilitásnak köszönhetően ígéretes kutatási területe napjaink numerikus áramlástanának (lásd [1,12]).
1.3.3. Véges térfogat módszere A véges térfogat módszer 1971-ben McDonald [46] által került bevezetésre, azon a megfigyelésen alapulva, miszerint a megmaradási törvények integrál alakban képesek kezelni a nem folytonos megoldásokat is. A szakadások akkor és csak akkor érvényes megoldásai a kiindulási parciális differenciálegyenleteknek, amennyiben un. gyenge megoldásként értelmezhetők, vagyis az egyenletek integrál alakjának is megoldásai. Ez numerikus szempontból legegyszerűbben az integrál egyenletek közvetlen diszkretizálásával érhető el. Hasonlóan a VEM módszerhez, az áramlási tér cellákra van felosztva. A cellák egyaránt lehetnek strukturáltak vagy nem strukturáltak, ami nagyobb rugalmasságot biztosít a véges differenciák módszerével szemben. A véges térfogatok, a geometriai diszkretizációnak megfelelően, egy vagy több cella szuperpozíciójával állíthatók elő. Lehetnek átlapoltak is, ellentétben a cellákkal, amelyek a konzisztencia biztosítása érdekében nem lapolhatják át egymást. A megmaradási törvények diszkretizálása az ezekre az ellenőrző térfogatokra történő egyenletfelírást jelenti. Az ismeretlen paraméterek a véges térfogatok feletti átlagolt paraméterek formájában öntenek alakot. A hálópontokhoz tartozó ellenőrző területek alak és elhelyezkedés lehetősége nagy szabadságot biztosít a számítási tér függvényreprezentációjának szempontjából. Ez a tulajdonság nem igaz a VDM és a VEM módszerekre, és talán ennek köszönhető a VTM növekvő népszerűsége a különféle mérnöki alkalmazásokban. Az integrál egyenletek közvetlen diszkretizálása a fizikai mennyiségek (tömeg-, impulzus- és energia-) megmaradásának biztosítása mellett szükségtelenné teszi a fizikai térből a számítási térbe történő transzformációt. A véges térfogat módszer magába foglalja a VEM módszer geometriai és a VDM diszkretizálási flexibilitását. Napjainkban, a hangsúly, a VDM módszer elterjedt használatáról egyre inkább áttevődik a VTM használatára, amely azonban nem közelíti meg a VEM módszer matematikai komplexitását. Sajnos a magasabb rendű deriváltak kezelése szintén nem könnyű feladat a VTM módszerrel. Talán ennek is köszönhetően, a VEM módszer szigorú matematikai eleganciája miatt (például az alacsonyabb rendű deriváltak előállítása a magasabb rendűekből a súlyozott hibák módszerével), fejlődésének üteme meghaladja a VDM és a VTM módszerekét [44]. 10
II Fejezet: Numerikus módszerek és alkalmazásuk az összenyomható áramlás modellezésére 2.1. Alapegyenletek – Euler-egyenletek Az Euler-egyenletek esetén eltekintünk a gázok belső súrlódásából, illetve a hővezetéséből adódó hatások figyelembevételétől. A parciális differenciálegyenlet rendszer elsőrendű hiperbolikus típusú, amelyben a nemlinearitás a konvektív tagok formájában van jelen. Az Euler-egyenletek (a tömeg-, az impulzus-, és az energiamegmaradás) felírhatók konzervatív, kvázi lineáris (nem konzervatív), illetve karakterisztikus formában. Konzervatív forma esetén a függő változók a sűrűség a sebesség-komponensek és a torlóponti energia segítségével fejezhetők ki. A kvázi lineáris forma esetén a függő (primitív) változók a sűrűség, a három sebesség komponens, illetve a statikus nyomás. A karakterisztikus formánál a független változók állandóak a hozzájuk tartozó karakterisztikus görbén. Elméleti szempontból mindhárom alak ekvivalens, mivel az egyik forma a másikból matematikai átalakításokkal állítható elő. Numerikus szempontból azonban a különféle formák nem ekvivalensek. A konzisztencia megtartása mellett, a különféle szakadások (lökéshullámok, kontaktfelületek, illetve örvényvonalak) modellezésének képessége az Euler-egyenletek integrálása után előálló egyenletrendszer belső tulajdonsága, ami azonban csak a konzervatív forma segítségével biztosítható. A primitív változók gradiensei ugyanis nincsenek definiálva a különféle szakadási felületeken, ezért a konzervatív alak segítségével írhatók le helyesen azok az áramlástani folyamatok, amelyekben előfordulhatnak a különféle diszkontinuitások. Ezáltal, függetlenül az eredményben megjelenő nem folytonos fizikai (primitív) állapotváltozóktól, az időben beállt folyamat létezése, illetve elérése esetén, valósul meg a tömeg, az impulzus és az energia egyensúly az ellenőrzőtartomány felett. A konzervatív formában felírt Euler-egyenletrendszer integrálását követően a különféle szakadások (lökéshullámok, kontakt- és örvényfelületek) az egyenletrendszer által elfogadható megoldások. Ennek alátámasztására szolgál a konzervatív alakú Euleregyenletrendszer diszkontinuitást tartalmazó, végtelen kis ellenőrzőfelület feletti integrálása után levezethető szakadási kapcsolata (jump relation), ami a jól ismert Rankine−Hugonoitegyenlet formájában ölt alakot. Mint ismeretes, az Euler-egyenletek leírják a szakadások nélküli izentrópikus áramlást, amelyben egy áramvonal mentén az entrópia állandó. Crocco szerint, az időben állandó, külső erők hatásától mentes nem viszkózus áramlásban az egyik áramvonalról a másikra azonban megváltozhat az entrópia. Az elmélet arra vezethető vissza, hogy a lokális sebességre merőleges entrópianövekedés kapcsolatban áll az örvényességgel, vagyis az entrópiaváltozás örvényességet, illetve az örvényesség entrópiaváltozást idéz elő. Ennek egyenes következménye például egy változó intenzitású lökéshullám után a hullámra merőleges sebességcsökkenések különbözősége miatt eltérő entrópianövekedés vagyis örvényesség megjelenése, még akkor is, hogyha a lökéshullám előtt teljesült is az örvénymentesség feltétele [28]. Matematikai szempontból, a konzervatív formában felírt Euler-egyenletrendszernek éppúgy lehet érvényes megoldása az entrópia csökkenéssel járó folyamat, mint az entrópia növekedéssel járó. Az expanziós lökéshullámok esetén, amelyeken keresztül entrópia csökkenés áll elő, szintén megoldása lehet az egyenletrendszernek, mivel a hőátadási tagok hiánya miatt az Euler-egyenletek megfordítható folyamatot modelleznek. Ezért nem lehetséges a különbségtétel a pozitív entrópia ugrással járó (kompressziós hullámok) és a negatív entrópia ugrással járó szakadások (expanziós hullámok) között. Lax [41] azonban megmutatta, hogy az olyan nem fizikai eredmények, mint például az expanziós hullámok elkerülése érdekében, egy pótlólagos disszipációs mechanizmust (más néven entrópiafeltételt (entropy condition)) kell explicit vagy implicit módon beépíteni a numerikus eljárásba. Az 11
entrópiafeltétel fogja garantálni, hogy a nem viszkózus egyenletek numerikus megoldása megfelelően reprezentálja a fizikai áramlást a zérushoz tartó viszkozitás határán. Ez egyenes következménye a termodinamika második főtételének. Az entrópia feltételről bővebb információ a [28,42] irodalomban található [44].
2.1.1. Euler-egyenletrendszer konzervatív formája A modern irodalomban az Euler-egyenletek alatt nem csak az impulzus-megmaradás egyenletei értendőek, hanem beletartoznak a kontinuitás és az energia-megmaradás egyenletei is. Az Euler-egyenletek két dimenzióban, konzervatív változókkal, divergencia formában a következőképpen írhatók fel [30]: ∂U ∂F (U ) ∂G (U ) + + =0, ∂t ∂x ∂y
(2.1)
U = [ρ , ρu , ρv , ρE ] ,
(2.2)
ahol T
[
]
(2.3)
[
]
(2.4)
T
F ( U ) = ρu , ρu 2 + p , ρuv , ρuh to , G( U ) = ρv , ρuv , ρv 2 + p , ρvh to
T
és amely könnyen kiegészíthető a harmadik térbeli koordinátával. A (2.3) és (2.4) egyenletekben szereplő konvektív típusú fluxusokban h to a torlóponti entalpia r2 ( h to = E + p / ρ ), E a torlóponti energia ( E = e + V / 2 ), amelyben a lokális sebesség r r r r nagysága: V = V = u 2 + v 2 ( V = uex + ve y ), valamint e = C v T = (1 / (γ − 1))( p / ρ ) a statikus energia. u és v a vízszintes és függőleges irányú sebességkomponens, ρ a sűrűség, T a statikus hőmérséklet, p a statikus nyomás, γ a fajhőviszony és C v az állandó térfogaton r r r vett fajhő. A H = Fe x + Ge y eredő konvektív fluxusvektor bevezetésével a (2.1) egyenlet összevont formában a következőképpen írható fel: ∂U v v + ∇ ⋅ H (U ) = 0 . ∂t
(2.5)
A (2.5) egyenletrendszer Γ oldalakkal határolt Ω ellenőrző felület feletti integrálását, valamint a Gauss−Osztogradszkij-tétel alkalmazását követően a következő egyenlet írható fel: rr ∂ UdΩ + ∫ HndΓ = 0 , (2.6) ∫∫ Γ ∂t Ω r amelyben n = (n x , n y ) a Γ interfészből kifelé mutató egységvektor. Az (2.6) egyenlet jelentése szerint az Ω tartománybeli U konzervatív változók időegység alatti megváltozása v egyenlő az Ω tartomány Γ felületén belépő, illetve eltávozó eredő H fluxussal.
2.1.2. Euler-egyenletrendszer kvázi-lineáris alakja A numerikus módszerek kidolgozása során a konzervatív forma elsőbbséget élvez ugyan, de a teljesség kedvéért, illetve az egyéb alkalmazásokra való tekintettel (pl. peremfeltételek)
12
érdemes foglalkozni az egyenletek primitív formában történő felírásával is: T P = (ρ , u , v , p ) . A (2.1) egyenlet kvázi-lineáris formában a következőképpen írható fel [16]: r ∂U ∂H v ∂U v v + ⋅ ∇ U = + C ⋅ ∇U = 0 , ∂t ∂U ∂t
(2.7)
r v ∂H ∂F r ∂G r = C= ex + ey . ∂U ∂U ∂U
(2.8)
amelyben:
Ekkor a (2.7) egyenletet a következőképpen lehet felírni: ∂U ∂U ∂U +A +B = 0, ∂t ∂x ∂y
(2.9)
amelyben A és B a Jacobi-mátrixok, azaz a fluxus-vektoroknak a konzervatív változók szerinti deriváltja: A=
∂F ∂F = ai , j = i , ∂U ∂U j
[ ]
0 1 αβ − u 2 (3 − γ )u A= v − uv 2 (− λE + 2αβ )u γE − αβ − βu
(2.10) 0 − βv β , u 0 − βuv γu 0
(2.11)
illetve B=
∂G ∂G = bi , j = i , ∂U ∂U j
[ ]
0 0 1 v u − uv B= 2 αβ − v (3 − γ )v − βu 2 (− λE + 2αβ )v − βuv γE − αβ − βv
(2.12) 0 0 , β γv
(2.13)
r2 amelyben α = V / 2 és β = γ − 1 . Legyen M a konzervatívból primitív változós formába
leképező transzformációs mátrix: M=
∂U ∂P
és M −1 =
∂P , ∂U
(2.14)
vagyis 1 u M = v α
0
ρ 0 ρu
0 1 − u / ρ 0 , illetve M −1 = − v / ρ ρ 0 ρv 1 / β αβ 0 0
0 1/ ρ 0 − uβ
0 0 1/ ρ − vβ
0 0 . 0 β
(2.15)
Ezek után a (2.1) egyenlet primitív változós formában a következőképpen írható fel: 13
∂P ~ ∂P ~ ∂P +A +B = 0, ∂t ∂x ∂y
(2.16)
~ ~ amelyben A = M −1 AM és B = M −1 BM : u ρ ~ ∂F 0 u A= = ∂P 0 0 2 0 ρc
0 0 0 1 / ρ u 0 0 u
v ~ ∂G 0 és B = = ∂P 0 0
0 ρ 0 v 0 0 . 0 v 1/ ρ 0 ρc 2 v
(2.17)
2.1.3. Euler-egyenletrendszer karakterisztikus alakja A (2.16) egyenlet kvázilineáris, hiperbolikus típusú elsőrendű parciális differenciálegyenlet rendszer. Az ilyen rendszerek fontos tulajdonsága a karakterisztikus görbék, felületek, illetve hiperfelületek (egy, kettő és három dimenzióban) létezése, ami szoros kapcsolatban áll a hullámterjedés jelenségével. A hiperbolikus egyenletrendszerek részletes ismertetése nem célja e munkának, ezekről bővebb információ a [13] irodalomban található. Az ilyen típusú egyenletek fontos jellemzője, hogy a Vn 0 ~ ~ ~ Cn = Anx + Bn y = 0 0
ρn x
ρn y
Vn
0
0 ρc 2 n x
Vn ρc 2 n y
0 n x / ρ ny / ρ Vn
(2.18)
mátrixnak léteznek valós sajátértékei és lineárisan független sajátvektorai [16]. A (2.18) r r T egyenletben n = (n x , n y ) tetszőleges irányú egységvektor ( n = 1), c a lokális hangsebesség: c=
γp ρ
(2.19)
r r és Vn a V sebességvektor n irányába eső merőleges vetülete: r r Vn = V ⋅ n = un x + vn y . (2.20) r Jelen esetben n tetszőleges, de elöljáróban érdemes megjegyezni, hogy a véges térfogat r ~ módszer esetén n az ellenőrző felületből kifelé mutató egységvektort fogja jelenteni. A C n mátrix sajátértékei a következő algebrai egyenlet segítségével állíthatók elő: ~ det C n − λI = 0 . (2.21)
A (2.21) egyenletben szereplő művelet elvégzése után a négy sajátérték a következő:
λ(1) = Vn , λ(2 ) = Vn , λ(3 ) = Vn + c
és
λ( 4 ) = Vn − c ,
(2.22)
~ amelyek három térbeli dimenzió esetén a Vn sajátértékkel egészülnek ki. A C n mátrix bal ~ oldali sajátvektorai, mint l (α ) sorvektorok a következőképpen írhatók fel: ~ (α ) ~ ~ l C n = λ(α ) l (α ) , (2.23)
vagy összevonás után mátrixos formában: 14
~~ ~ L C n = Λn L ,
(2.24)
0 0 Vn 0 0 V 0 0 n Λn = 0 0 Vn − c 0 0 Vn + c 0 0
(2.25)
amelyben Λn a sajátérték mátrix:
~ és Ln mátrix pedig a bal oldali sajátvektor mátrix:
1 0 0 ny ~ Ln = 0 n x 0 − n x
0 − nx ny − ny
~ − 1 / c 2 ln (1) ~ 0 ln (2 ) = ~ . 1 / ρc ln (3 ) ~ 1 / ρc ln (4 )
(2.26)
~ A Λn sajátérték mátrix a C n mátrix diagonizálásával állítható elő: 0 0 λ(1) Vn 0 0 0 0 ~ ~ ~−1 0 Vn = Λn = Ln C n Ln = 0 0 Vn − c 0 0 Vn + c 0 0 0 0
0
0
(2 )
0
0
λ( 3 )
0
0
λ
0 0 , 0 λ(4 )
(2.27)
amelyben a jobb oldali sajátvektor mátrix: 0 1 0 n ~ ~ y Rn = Ln−1 = 0 − n x 0 0
ρ / (2c ) ρ / (2c ) nx / 2 ny / 2 ρc / 2
− n x / 2 ~ (1) ~ (2 ) ~ (3) ~ (4 ) = rn , rn , rn , rn − n y / 2 ρc / 2
(
)
(2.28)
~ alakban írható fel. Az 1D-s esethez hasonlóan a karakterisztikus változók a Ln mátrixszal balról történő beszorzása után állnak elő [16]: ~ ~ ∂Wn = Ln ∂P illetve ∂P = Ln−1∂Wn = ∑ ~ rn(α )∂Wn(α ) , (2.29) α
vagyis: 1 ∂ρ − 2 ∂p ∂Wn r cr ∂V s (2 ) ∂ W r n = r 1 ∂Wn = , ( ) 3 ∂ + ∂ n V p ∂Wn ρ c (4 ) r r 1 ∂Wn ∂p − n∂V + ρc (1)
(2.30)
r r r r ahol s = n y e x − n x e y n vektorra merőleges, vele jobbrendszert alkotó egységvektor. (Ekkor a karakterisztikus egyenlet a következőképpen írható fel (1D-ben):
∂W ∂W +Λ = 0 .) ∂t ∂x
15
A karakterisztikus változók természetesen érvényesek a konzervatív formában felírt egyenletrendszerre is:
Λ n = Ln C n L−n1 ,
(2.31)
amelyben Cn = Anx + Bn y , vagyis: 0 αβ n − uV x n Cn = αβ n y − vV ( −γE + 2αβ )Vn
nx
ny
Vn + (1 − β )un x
un y − βvn x
vn x − βun y ( γE − αβ )n x − β uV n
Vn + ( 1 − β )vn y ( γE − αβ )n y − βvVn
0 βn x , βn y γVn
(2.32)
valamint
αβ 1− c2 − un + vn y x ρ αβ Ln = − V n + c ρ αβ Vn + c ρ
uβ c2 ny
vβ c2 − nx
ρ
ρ
uβ nx − c
ny −
uβ − nx + c
− ny +
ρ
ρ
ρ
−
vβ c
ρ
vβ c
β
c2 0 l (1) n l n(2 ) = , β l n(3) ρc l (4 ) n β ρc
(2.33)
illetve 1 u −1 Rn = Ln = v α
ρ
nx
ρn y − ρn x
ρ 2c
ρ
ρ
2c
(u + cn x )
(v + cn ) y
2c ρ c2 α + + cVn ρ ( un y − vn x ) 2c β
2c ρ (u − cn x ) 2c . ρ (v − cn y ) 2c γVn
(2.34)
A karakterisztikus változók, az előzőekhez hasonlóan, a következőképpen állíthatók elő [16]: ∂Wn = Ln ∂U
vagy ∂Wn(α ) = l n(α )∂U .
(2.35)
2.1.4. Navier−Stokes-egyenletrendszer konzervatív formája A teljesség kedvéért érdemes megemlíteni a Navier-Stokes (NS) -egyenleteket, amelyek az alapvető fizikai megmaradási törvényekből származtathatók. Történeti szempontból a NSegyenletek a viszkózus impulzus-megmaradási egyenleteket foglalják magukba, felírva Newton II. törvényét egy végtelenül kicsi rögzített ellenőrzőfelületre, de a modern áramlástani irodalomban azonban a NS-egyenletek a tömeg-, illetve az energia-megmaradás egyenleteit is tartalmazzák. A külső erők, hőbevezetés, anyag diffúzió és véges állapotú kémiai reakciók nélküli, instacioner 2D-s Navier−Stokes-egyenletek konzervatív formában a (2.1) egyenlet kibővítésével a következőképpen írhatók fel [30]:
16
∂U ∂F (U ) ∂G (U ) ∂F v (U ) ∂G v (U ) + + = + , ∂t ∂x ∂y ∂x ∂y
(2.36)
ahol
0 0 τ yx τ xx v G = és Fv = . τ τ xy yy ∂T ∂T uτ yx + vτ yy + κ uτ xx + vτ xy + κ ∂y ∂x
(2.37)
A viszkózus fluxusbeli τ ij feszültségek magukba foglalják egyrészt a normál feszültségeket ( i = j ):
τ xx = 2 µ
∂u 2 v v ∂v 2 v v − µ (∇ ⋅ V ) és τ yy = 2µ − µ (∇ ⋅ V ) , ∂x 3 ∂x 3
(2.38)
valamint a nyírófeszültségeket ( i ≠ j ): ∂v
∂u
τ xy = τ yx = µ + , ∂x ∂y
(2.39)
amelyben µ a molekuláris viszkozitás.
17
2.2.
Véges térfogat módszer
A numerikus eljárások kidolgozása során, az 1.3.3. alfejezetben leírtak miatt, a véges térfogat módszert alkalmaztam az algebrai egyenletek előállítására. Mint ahogyan már az 1.3.3. alfejezetben volt róla szó, a véges térfogatok módszere az alapegyenletek integrálása után megjelenő kifejezések diszkretizálásán alapszik. A numerikus fluxusok a Gausz−Osztogradszkij-tétel alkalmazásával, az integrál egyenletek térbeli operátorainak átalakításával állíthatók elő. A következő lépés a (2.6) egyenletben szereplő H fluxus kontúrintegráljának meghatározása. Nyilvánvaló, hogy a kontúrintegrált a fluxusok felületelemre merőleges komponenseivel egyszerűbb meghatározni, mint a koordinátarendszer irányai szerinti összetevők segítségével, ezért teljes fluxus n irányba eső komponensét a következő alakban írtam fel: ρVn ρuVn + pn x H n = Hn = , ρvVn + pn y to ρ Vn h
(2.40)
Vn ahol ( Vn = Vn = ( ue x + ve y ) ( n x e x + n y e y ) = un x + vn y a normális sebességkomponens. Ezek után (2.6) egyenlet a következőképpen írható fel:
∂ ∫∫ UdΩ + ∫Γ H n dΓ = 0 . ∂t Ω
irányú
(2.41)
A folytonosról a diszkrét alakra történő átalakítás során a véges térfogat feletti U megoldásfüggvény az U j diszkrét ismeretlennel közelítettem: Uj =
1 Ωj
∫∫UdΩ ,
(2.42)
Ω
amely a cella központú megközelítésnek felel meg. Ezért, az ismeretlen U megoldásvektor inkább értelmezhető az Ω ellenőrző felület feletti középértéknek, mintsem egy csomóponti függvényértéknek, amely a cella csomóponti közelítés jellemző paramétere. A Ω j ellenőrző felület N f számú kontúrjára elvégzendő összeggel helyettesítve a (2.41) egyenletbeli második integrált, a (2.1) egyenlet a következő j pont körüli szemi-diszkretizált differenciálegyenletre vezethető vissza: i,j+1 nij,1 Γ ij,4
i-1,j nij,4
i-1/2,j
x
R L i+1/2,j Γ ij,2
i,j Ω
R L i,j-1/2 Γ ij,3
y ey
i,j+1/2 Γ ij,1
ij
L
L R
i+1,j nij,2
R nij,3
i,j-1 2.eábra. Cella központú véges térfogat háló modell x 18
[ H n ] j ,k
az Ω j ellenőrzőfelülethez tartozó, a j és k csomópontok között elhelyezkedő, a Γ j ,k oldalhosszra jellemző nem viszkózus fluxus (2. ábra). A [ H n ] j ,k mennyiség az elkövetkezendőkben, mint numerikus fluxusfüggvény fog szerepelni és amely a legalacsonyabb szintű közelítéssel az ellenőrző felület feletti állandó paraméterek segítségével a következőképpen állítható elő: ~ ~ Hn U L + Hn U R ~ ~ L R . (2.44) H n = H n U ,U = 2 amelyben
(
)
( )
( )
L Az egyenletben szereplő U = U j (L: Left (bal)) és U R = U k (R: Right (jobb)) paraméterek a
Γ j ,k határ két oldalán elhelyezkedő állapotváltozók. A diszkretizációról a 2.2.2. alfejezetben található bővebb információ.
2.2.1.
magasabbrendű
térbeli
Konvektív tagok térbeli diszkretizációja
Az Euler-egyenletek jól közelítik a transzszonikus áramlásban lezajló folyamatokat. Az egyenletrendszer numerikus módszerekkel történő megoldása esetén azonban némi disszipációnak kell fellépni ahhoz, hogy a séma megőrizze numerikus stabilitását, illetve csökkenjenek, vagy megszűnjenek a hibás oszcillációk. Ezt alapvetően kétféle módon lehet biztosítani: mesterséges viszkozitás bevezetésével, valamint az úgynevezett áramlásirányú (upwind) módszerek alkalmazásával. Az első módszer esetén a centrális séma alkalmazása mellett egy explicit hozzáadott magasabbrendű disszipáció, un. mesterséges viszkozitás (pl. Jameson-féle mesterséges disszipáció [35]) segítségével biztosítható a numerikus stabilitás. A módszer szimmetrikus tulajdonságokkal rendelkezik a Jacobi-mátrix sajátértékeinek előjelváltására, vagyis nincs tekintettel a hullámterjedés irányítottságára. Ezért, a hiperbolikus egyenletekre jellemző karakterisztikák menti zavarásterjedés nincs figyelembe véve a numerikus diszkretizáció során. A áramlásirányú (upwind) módszerek azonban túlmutatnak ezen a problémán, mivel alkalmazásuk esetén az információterjedésnek megfelelően vannak szétválasztva az áramlástani paraméterek, így biztosítva a zavarások fizikai terjedésének optimális lehetőségét. Az upwind típusú módszereken belül a fluxusvektor megosztáson (flux vector splitting) (pl. Steger and Warming [70], Van Leer [76]), illetve a fluxuskülönbség megosztáson (flux difference splitting) (Godunov típusú módszerek: közvetlen és Roe [55] vagy Rumsey [58,59] által közelített Riemann problémák) alapuló módszerek a legelterjedtebbek. A számítógépi kapacitást tekintve ezek a módszerek költségesebbek ugyan, mint a mesterséges viszkozitás módszere, de nagyobb Mach-számok, valamint az instcioner folyamatok esetén előálló mozgó lökéshullámok esetén pontosabban, hangolás nélkül képesek modellezni az áramlást. A teljesség, illetve napjaink igen széles: a szubszonikustól a szuperszonikus sebességtartományig nagy pontossággal alkalmazható numerikus módszereinek bemutatása érdekében, a mesterséges viszkozitás módszerének részletes ismertetését követően, a fluxusvektor és a fluxus differencia megosztáson alapuló módszerek is tárgyalásra kerülnek. Az új karakterisztikus peremfeltétel számítási eljárást a Roe által közelített Riemann módszerben alkalmaztam. 2.2.1.1.
Mesterséges viszkozitás módszere
Minden másodrendű centrális diszkretizáció alkalmazásakor (stabil időlépés esetén is) fellép némi instabilitás. A probléma gyökere az a nem fizikai megoldás halmaz, amely ugyan szintén kielégíti a diszkretizált egyenletrendszert, azonban hamis eredményt szolgáltat. A jelenségnek köszönhetően a páros és a páratlan hálópontok közötti – a séma által biztosított – folytonos kapcsolat megszűnik. A tisztán numerikus hiba oka, egy másodrendű diszkretizáció 19
alkalmazása egy elsőrendű differenciáloperátorra. Egy másik jellegzetes probléma az oszcillációs jelenség, amely az olyan nagy gradiensek környezetében jelentkezhet, mint például a lökéshullámok vagy a kontaktfelületek. A diszkretizáció nem képes kezelni a szakadás környezetében előálló nagyfrekvenciás jeleket, mivel 2Δx az a legnagyobb hullámhossz, ami a hálón még megoldható. A hamis eredmény, illetve az oszcillációs jelenség elkerülése érdekében érdemesnek tűnik magasabbrendű tagok hozzáadásával javítani az eljáráson. A disszertációban használt mesterséges viszkozitás módszert eredetileg Jameson [35] fejlesztette ki és napjainkig ez a módszer alkalmazható a legjobb számítógépi hatásfokkal a hangsebességhez közeli áramlások szimulációjára. Ezért ezt az eljárást alkalmaztam a transzszonikus lapátrács numerikus áramlástani modellezésében a 2.6. alfejezetben. A (2.43) egyenletrendszer az explicit stabilizáló taggal a következőképpen írható fel: N ∂ 1 f , Uj =− [ H ] Γ − D ∑ n j , k j ,k ∂t Ω j k =1
(2.45)
amelyben a (2.40) kifejezés alapján a már ismert fluxus: ρVn ρuVn + pnx H n = Hn = , ρvVn + pn y to ρVnh
(2.46)
valamint a D = Dx + D y vagyis Dx = d i + 1 , j − d i − 1 , j , illetve D y = d i , j + 1 − d i , j − 1 2
2
2
2
(2.47)
a disszipációs tag. Ez utóbbi két egyenlet jobb oldalán hasonló kifejezések találhatók, például az ( i , j ) és az ( i + 1, j ) csomópontok között felírható:
d
1 i+ , j 2
=
Ωi , j ∆t
(2) ε 1 (U i +1 , j − U i , j ) − ε ( 41) (Ui +2 , j − 3U i +1 , j + 3U i, j − U i −1, j ) . (2.48) i+ , j i+ , j 2 2
Az ε ( 2 ) és ε ( 4 ) tényezők meghatározása a következő nyomásszenzor segítségével történik: ν i ,j =
pi +1 , j − 2 pi , j + pi −1, j pi +1, j + 2 pi , j + pi −1, j
,
(2.49)
így: ε ( 2)1 = κ ( 2) max(ν i +1 , j ,ν i , j ) i+ , j 2
( 4) (4 ) (2 ) , illetve ε i+ 1 = max 0, κ − ε i + 1 , j , 2
2
(2.50)
1 1 és κ ( 4 ) = a numerikus hangoló paraméterek. A másodrendű 4 256 disszipáció, d ( 2 ) , az ε ( 2 ) tényező miatt csak a szakadás közelében kap értéket, elsőrendűvé alakítva a diszkretizációt, ami a monoton (oszcilláció-mentes) lökéshullám meghatározásának feltétele. Ezzel ellentétben a negyedrendű disszipáció, d ( 4 ) , a szakadáson kívül mindenhol érvényes, így biztosítva a séma másodrendűségét, valamint a nemlineáris instabilitások és a páros páratlan szétválasztódás elkerülését. 2.2.1.2. Áramlásirányú (upwind) módszerek amelyekben κ ( 2 ) =
20
A áramlásirányú módszerek esetén, a mesterséges viszkozitás módszer alfejezetben tárgyalt oszcillációs jelenségek elkerülése érdekében nincs szükség új egyenletbeli tagok hozzáadására, mivel a numerikus fluxusok az áramlás irányának, vagyis a karakterisztikus görbék meredekségének (a fizikai információterjedés irányának) megfelelően vannak szétválasztva. Legyenek az U L = U i , j és U R = U i +1 , j
(2.51)
( i , j ) , valamint a ( i + 1, j ) szomszédos cellákra jellemző konzervatív állapotjelzők, amelyeket az ( i + 1 2 , j ) cellahatár választ el egymástól. A (2.8) egyenlet alapján, az n irányba levetített, a teljes fluxus konzervatív változók szerinti deriváltja, a (2.32) egyenlet egyenletben szereplő Cn mátrix a következőképpen állítható elő: Cn =
∂H n . ∂U
(2.52)
Az áramlási irány szerint differenciált (upwind) sémák numerikus fluxusfüggvényei, Harten és mások [27] által, a következőképpen definiálhatók: ~ (2.53) H n U L ,U R = H n U L , ha λ( α ) > 0 és
(
)
(
~ H n U L ,U R
( ) ) = H (U ) , ha λ( R
n
α)
< 0,
(2.54)
amelyben λ( α ) C n mátrix sajátértéke. A fenti feltétel a hangsebesség feletti áramlás esetén azt a fizikai hatást modellezi, miszerint nincs információterjedés az áramlással ellentétes irányba. Tetszőleges, U L és U R paraméterek segítségével, egy hozzájuk közeli U ∗ referenciaállapottal a fluxusfüggvény Taylor-sorba fejtése a következőképpen írható fel: ~ H n U L ,U R = H n U ∗ + Cn+ U ∗ U L − U ∗ + Cn− U ∗ U R − U ∗ +
(
) + Ο ((U
L
( ) ( )( − U ) , (U − U )(U ∗ 2
L
∗
)
R
( )(
)(
−U ∗ , U R −U ∗
) ).
)
2
(2.55)
Az egyenletben szereplő Cn+ és Cn− mátrixok a jobb és bal oldali sajátérték mátrixok: C n+ = L−1 Λ+ L és
(2.56)
Cn− = L−1Λ− L ,
(2.57)
valamint Λ+n pozitív és Λ−n negatív diagonális mátrixok segítségével állíthatók elő. A diagonális mátrixokban szereplő negatív, illetve a pozitív értékek helyére zérust kell írni. Legyen U ∗ U L és U R aritmetikai középértéke: U L + U R / 2 . Ekkor a (2.55) numerikus fluxusfüggvény a következőképpen írható fel:
(
~ Hn
i+
1 2
=
)
1 U i + U i +1 ( U i +1 − U i ) = H n ( U i ) + H n ( U i +1 ) − Cn 2 2 =
1 H n ( Ui ) + H n ( Ui +1 ) − δHn i + 1 . 2 2
(2.58)
A δH n tagnak ki kell elégítenie a kozisztenciát, vagyis: δH n (U ,U ) = 0 . 21
(2.59)
δH n lényegében egy mesterséges disszipációs tag, amely a centrális diszkretizáció stabilizálására szolgál. Az upwind jellegű sémák többségében, a közöttük lévő különbség a δH n tag különféle, precíz meghatározásában rejlik. Ennek értelmében megkülönböztethetők: a. Fluxusvektor megosztáson alapuló módszer (Flux Vector Splitting Technique) b. Fluxuskülönbség megosztáson alapuló módszer (Flux Different Splitting Technique) (Godunov típusú módszerek a Riemann problémák megoldásával) b.1. Riemann megoldók (Exact Riemann Solvers) b.2. Riemann problémák közelítő megoldásai (Approximate Riemann Solvers) Az a. és a b.2. módszerek részletes ismertetése előtt a Riemann problémát pontosan megoldó Godunov (1959) által kidolgozott eljárásokkal kapcsolatban érdemes megjegyezni, hogy a módszerben szereplő térbeli cellák fölötti állandó értékkel közelített megoldások pontos kiszámítása számítástechnikailag költséges és sémában rejlő approximáció miatt nem adnak pontosabb eredményeket, mint a közelítő Riemann megoldók. a. Van Leer fluxusvektor megosztáson alapuló módszere (Flux Vector Splitting Method) A fluxusvektor megosztáson alapuló módszerekben a numerikus fluxus a cellahatár jobb és a baloldali – a megfelelő terjedési irány figyelembevételével megosztott – fluxusok összegeként írható fel [44]: ~ H n U L ,U R = H n+ U L + H n− U R . (2.60)
(
)
( )
( )
A két legelterjedtebb módszer a Steger and Warming [70] és a Van Leer [76] fluxusvektor megosztás módszere. A Steger and Warming módszer esetén azonban, a szónikus, illetve a torlóponti helyeken, vagyis ott, ahol a Jacobi-mátrix sajátértéke előjelet vált, instabilitás, oszcillációs megoldás állhat elő. Ezekben a pontokban ugyanis nem differenciálható a séma osztott fluxusa. Van Leer módszerében [76] nem jelentkezik ez a probléma, mivel fluxusoknak az áramlás irány szerinti differenciálásuk mellett ki kell elégíteniük néhány megkötést. Ennek alapján az osztott fluxusok és deriváltjaik a Mach-szám folytonos polinom függvényei, amelyek a szuperszonikus áramlás esetén a teljes fluxussá alakulnak. A másik fontos tulajdonság a szimmetria, amely megteremti a kapcsolatot a –1 és a 1 Mach-számok m között. Ennek értelmében a fluxusok a (1 ± M n ) m ≥ 2 polinommal arányosak. A 2D-s Van Leer fluxusok levezetése egyenesen következik az egydimenziósból, valamint a normális irányú Mach-szám bevezetésével ( M n = Vn / c ), a − 1 < M n < 1 sebességtartományban a (2.60) egyenletbeli fluxusok a következőképpen írhatók fel [28]: 1 c u + ( ± 2 − M n ) nx γ ρc 2 ± c . Hn = ± (1 ± M n ) v + ( ± 2 − M n ) ny 4 γ 2 2 c 2 ± 2 M − M 2 + V n n γ + 1 γ − 1 2
(2.61)
Ha a Mach-szám cellahatárra merőleges vetülete nagyobb mint egy, azaz M n ≥ 1 , akkor: H n+ = H n és H n− = 0 ,
(2.62)
H n+ = 0 és H n− = H n .
(2.63)
ha M n ≤ −1 , akkor:
22
Könnyen ellenőrizhető, hogy ha a H n+ és H n− fluxusok a (2.60) egyenletnek megfelelően összeadódnak, akkor a (2.46) egyenletbeli teljes normál fluxus áll elő. Az eljárás gyorsan konvergál és monoton tulajdonságokkal (oszcilláció-mentesség) rendelkezik, azonban az elsőrendű térbeli diszkretizációnak köszönhetően igen pontatlan. A nagy disszipáció miatt elsősorban a lökéshullámok pontos helyzetének meghatározása jelent nehézséget. A nagy pontosság (high resolution) elérése érdekében az approximáció rendjét egyről kettőre kell növelni, valamint a másodrendű sémákra jellemző lökéshullám körüli oszcilláció megszüntetése érdekében un. limitereket kell a sémába beépíteni (2.2.2. alfejezet). b. Roe által közelített Riemann eljárás (Roe’s approximate Riemann solver) Roe módszere az egyik legnépszerűbb közelítése az egzakt Riemann módszernek. Az eljárás a fluxuskülönbségek karakterisztikus mező-dekompozícióján alapszik. A (2.3) és (2.4) nem viszkózus fluxusok a konzervatív változók elsőfokú függvényei, ami azt jelenti, hogy: F ( β U ) = βF ( U ) , ∀ β ∈ R .
(2.64)
Képezzük a (2.64) egyenlet β szerinti deriváltját β = 1 helyen: ∂F ( βU ) ∂β
=F = β =1
∂F ∂( β U ) ⋅ = AU . ∂( β U ) ∂β
(2.65)
A (2.65) egyenlet alapján, valamint felhasználva, hogy a Cn mátrix a H n eredő fluxus konzervatív változók vektora szerinti differenciálásával állítható elő, a H n fluxus formálisan a következőképpen írható fel:
Hn =
∂H n U = CnU = RnΛn LnU , ∂U
(2.66)
( )
amelyben Λn Cn mátrix sajátértékeinek diagonál mátrixa: Λn = Diag λ(nα ) , α = 1, 2 , 3, 4 két dimenzióban: λ(n1) = Vn , λ(n2 ) = Vn , λ(n3 ) = Vn + c , λ(n4 ) = Vn − c
(
(c =
)
γRT a hangsebesség),
valamint Rn és Ln a jobb és bal oldali sajátvektor mátrixok Rn = L−n1 . Az információterjedés irányának megfelelően, a Λn mátrix elemeinek előjele szerint egy pozitív és egy negatív mátrixra bontható fel, mialatt a negatív és pozitív elemek helyére zérus kerül. Ezzel a megosztással az eredő fluxus a következőképpen írható fel:
(
)
(
)
H n = Rn Λ+n + Λ−n LnU = C n+ + C n− U = H n+ + H n− .
(2.67)
Az előzőek értelmében, a cellahatáron megjelenő fluxusfüggvény a bal és a jobb oldali állapotoknak megfelelően alakul: ~ (2.68) H n U L ,U R = H n+ U L + H n− U R ,
(
)
( )
( )
vagy
(
)
( )
( )
( )
( )
( )
( )
( )
UR ~ H n U L ,U R = H n U R − H n+ U R + H n+ U L = H n U R − ∫ Cn+ ( U ) dU . (2.69) UL
Másrészt szintén felírható, hogy
(
)
( )
UR ~ H n U L ,U R = H n U L − H n− U L + H n− U R = H n U L + ∫ C n− (U ) dU . (2.70)
23
UL
A megoldás azon a lineáris hullám dekompozíció elven alapul, amelynek során a jobb és a bal oldalnak létezik egy speciálisan átlagolt állapota (kalappal jelölve), amelyet Roe paramétervektornak neveznek. Az integrálás elvégzését követően az előző két egyenlet segítségével a következő összefüggés adódik:
(
{ ( )
)
( )
(
)}
)(
1 ~ H n U L ,U R = H n U L + H n U R − Cˆ n U L ,U R U R − U L . 2
(
(2.71)
)
A Cˆ n U L ,U R mátrix a következő speciális tulajdonságokkal rendelkezik: 1. Valós értékű sajátvektorai és lineárisan független sajátvektorai vannak. 2. Ha U L = U R = U , akkor Cˆ n U L ,U R = Cˆ n (U ) . 3. H U R − H U L = Cˆ U L ,U R U L − U R . n
( )
n
( )
n
(
(
)(
)
)
Az első feltétel a hiperbolikus típusú egyenletek jellemző tulajdonsága, a második a konzisztencia. A harmadik feltétel biztosítja, hogy nyugvó lökéshullámok esetén, amikor H n U R = H n U L , a (2.71) numerikus fluxusfüggvény a fizikailag releváns
( ) ( ) 1 2 [ H (U ) + H (U ) ] R
L
fluxussá alakuljon. Más szóval ez azt jelenti, hogy a séma képes legyen a diszkontinuitás megoldására egy cellán belül, ha a szakadás hálóirányú [44]. A (2.71) L R egyenletbeli Cˆ n U ,U ∆U vektor a jobb oldali sajátvektorok mátrixának, a n
n
(
)
karakterisztikus változók vektorának és a sajátérték mátrixnak a lineáris kombinációjával állítható elő: 4
ˆ ∆W = R ˆ Λˆ ∆W = ∑ λˆ ( i ) ˆr ( i ) ∆W ( i ) , (2.72) Cˆ n ∆U = Cˆn Rˆn ∆W n = Rˆn Lˆn Cˆn R n n n n n n n n i =1
amelyben Rˆn jobb oldali sajátvektor mátrix rˆn( i ) oszlopvektorok, a ∆Wn vektor pedig a ∆Wn( i ) karakterisztikus változók segítségével állítható elő. Roe megmutatta, hogy ideális gázok esetén Cˆ n mátrix abban az esetben egyezik meg C n Jacobi-mátrixszal ha Cˆ n felírható ρˆ , uˆ , vˆ , és hˆ to változók függvényeként, amelyek a sűrűség négyzetgyökével átlagolt paraméterek formájában, Roe paramétervektoraként állíthatók elő: ρˆ = ρ χ , uˆ = L
h to,R χ + h to ,L uRχ + u L vRχ + vL to ˆ ˆ ,v= , illetve h = , 1+ χ 1+ χ 1+ χ
(2.73)
ahol amihez a számítástechnikailag leggazdaságosabb egy χ = ρR ρL , négyzetgyökszámítás szükséges. A hangsebesség az előbbi paraméterek segítségével szintén kifejezhető:
(
)
ˆ to uˆ 2 + vˆ 2 cˆ = ( γ − 1) h − . 2
(2.74)
Ezek után Roe fluxusfüggvénye a következőképpen írható fel:
(
)
( )
( )
4 1 ~ H n U L ,U R = H n U L + H n U R − ∑ λˆ (ni ) ˆrn( i ) ∆W n( i ) , 2 i =1
(2.75)
amelyben a jobb oldali sajátvektor mátrix:
24
1 uˆ Rˆn = vˆ α
ρˆ 2cˆ
0 ρˆ n y − ρˆ n x ρˆ ( uˆn y − ˆvnx )
ρˆ ( uˆ + cˆn x ) 2cˆ ρˆ ( vˆ + cˆn y ) 2cˆ ρˆ cˆ 2 α + + ˆcVn 2cˆ γ −1
ρˆ 2cˆ
ρˆ ( uˆ − cˆn x ) 2cˆ , (2.76) ρˆ ( vˆ − cˆn y ) 2cˆ ρˆ cˆ 2 α + − ˆcVn 2cˆ γ −1
2 ahol Vˆ = uˆe x + vˆe y , α = Vˆ 2 = uˆ 2 + vˆ 2 2 , és Vˆn = Vˆn = ( uˆe x + vˆe y ) ⋅ ( n x e x + n y e y ) =
(
)
= uˆn x + vˆn y . rˆni oszlop vektorok Rˆn mátrix oszlop elemei ( i = 1, 2, 3, 4 balról jobbra). A karakterisztikus vektor elemei a következőképpen írhatók fel: 1 ∆ρ − 2 ∆p ∆Wn cˆ s ∆V (2 ) ∆W n 1 ∆W n = = + n∆ V + ∆p , ( 3) ∆ W ˆ n ˆ ρ c ∆W ( 4 ) 1 n − n∆ V + ∆ p ρˆ cˆ (1 )
(2.77)
ahol s = n y e x − n x e y n -ral jobbrendszert alkotó, rá merőleges egységvektor, a ∆ -val jelölt paraméterek jelentése: ∆φ = φ R − φ L . Végezetül a Λˆ n sajátérték mátrix kifejtve következőképpen néz ki: (1 ) Vˆn 0 0 0 0 0 0 λˆ n ( 2) 0 Vˆn 0 0 0 0 λˆ n 0 Λˆ n = = (2.78) . ( 3 ) ˆ ˆ ˆ 0 0 λ 0 0 0 V n + c 0 n 0 0 0 ˆ ( 4) ˆ ˆ 0 0 λ 0 V n − c v n A numerikus eljárást azonban célszerű kibővíteni az un. entrópia feltétellel, hogy elkerülhetők legyenek az olyan nem fizikai megoldások, mint például az expanziós hullámok (2.1. alfejezet). Az entrópia korrekció lényegében egy hozzáadott numerikus disszipációt jelent, azokban a kritikus esetekben, amikor a sajátérték előjelet vált (pl. lökéshullám, hangsebesség). A disszertációban alkalmazott módszer Yee (1989) által javasolt eljáráson alapszik, amelyben a sajátértékek λ(n1) - től λ(n4 ) -ig egy un. entrópia-fix függvény segítségével állíthatók elő [32]:
( )
∆ λ(ni ) = ψ λ(ni ) ,
(2.79)
ahol a ψ függvény: ψ (z) = cˆ 2δ valamint
25
z ≥δ, cˆ z z2 2 + δ 2 , ha <δ, cˆ cˆ z , ha
(2.80)
uˆ + vˆ δ = δ 0 + 1 , cˆ
(2.81)
amelyben δ 0 állandó, értéke 0,025. Roe által közelített fluxuskülönbség megosztásán alapuló módszer a van Leer fluxusvektor megosztáson alapuló módszerhez hasonlóan elsőrendű. A séma pontosítása a magasabbrendű térbeli diszkretizáció bevezetésével érhető el (2.2.2. alfejezet).
2.2.2.
Magasabbrendű térbeli diszkretizáció (MUSCL (Monotone Upstream Schemes for Conservation Laws)) Az elsőrendű van Leer és Roe eljárások pontosságának javítása érdekében a magasabbrendű térbeli approximáció alkalmazása elkerülhetetlen. A diszkretizáció pontosságának rendje az U L és U R állapotok megfelelő függvénykapcsolatának figyelembevételével határozható meg a numerikus fluxusok előállítása során. Mint ahogyan már korábban volt róla szó, az első rendű sémák esetében az ellenőrző felület felett állandó értékkel szerepelnek a paraméterek: U iL+1 2 , j = U i , j és U iR+1 2, j = U i +1, j
(2.82)
az ( i + 1 2 , j ) cellahatáron az ( i , j ) , illetve az ( i + 1 , j ) véges térfogat elemek között. A magasabbrendű approximáció érdekében az állandó paraméterek lineárisan változó, másod-, vagy magasabbrendűekkel cserélendők fel. Ennek értelmében felírható U ( x ) megoldás függvény U i pont körüli Taylor-sorba fejtése (folytonos függvények esetén): ∂U U ( x) = U i + ∂x
1 ∂ 2U x i ( x − xi ) + 2 ∂x 2
( x − xi ) 2 + Ο ( ∆x 3 ) ,
xi
(2.83)
továbbá legyen: ∂U ∂x ∂ 2U ∂x 2
xi
xi
=
=
U i +1 − U i −1 δ U + Ο ∆x 2 = i + Ο ∆x 2 és 2∆ x ∆x
(
)
(
)
(2.84)
U i+1 − 2U i + U i−1 δ i2U 2 + Ο ∆ x = + Ο ∆x 2 . 2 ∆x 2 ∆x 2
(
)
(
)
(2.85)
A (2.84) és a (2.85) egyenletek segítségével a (2.83) egyenlet a következőképpen írható fel: δ iU δ i2U U ( x) = U i + ( x − x i ) + 2 ( x − xi ) 2 + Ο ∆x 3 , ∆x 2 ∆x
( )
(2.86)
amely az x ∈ [ xi − (1 2) ∆x , xi + (1 2 ) ∆x ] tartományban érvényes. A cellaközpontú eljárás során az állapotváltozók az ellenőrzőfelület feletti átlagolt paraméterek formájában kerülnek meghatározásra, vagyis: Ui =
1 ∆x
1 xi + ∆ x 2
∫ U ( x ) dx .
(2.87)
1
xi − ∆ x 2
26
A (2.86) egyenlet a (2.87) egyenletbe való behelyettesítése, és az integrálás elvégzése után a következő összefüggés adódik: Ui = Ui +
δ i2U + Ο ∆x 3 . 24
(
)
(2.88)
A (2.88) egyenlet (2.86) egyenletbe való visszahelyettesítése után, pedig felírható, hogy: 2 2 δ iU ( x − xi ) + δ i U2 ( x − xi ) 2 − ∆x + Ο ∆x 3 . ∆x 2 ∆x 12
(
U ( x) = U i +
)
(2.89)
Az egyszerű kezelhetőség érdekében érdemes beszorozni a szögletes zárójelben lévő tagot 3κ paraméterrel. Ennek értelmében, κ = 1 3 esetén a (2.89) egyenlet harmadrendűen pontos és másodfokú közelítést jelent U ( x ) térbeli diszkretizációja tekintetében. Ha azonban κ ≠ 1 3 , de − 1 ≤ κ ≤ 1 , akkor a (2.89) egyenlet lineáris, a levágási hiba nagysága κ értékétől függően különböző lehet, ami az approximáció rendjét a legtöbb esetben másodrendűre csökkenti. A többdimenziós problémák esetén az ismeretlen paraméterek vektorának lineáris vagy másodfokú előállítását minden egyes cellahatár normális irányba el kell végezni. A numerikus fluxusfüggvényre csak a cellahatáron van szükség, ezért az x = x i ± ∆x 2 (uniformis hálózás feltételezésével), illetve az U = U helyettesítéssel a (2.89) egyenlet a következőképpen írható fel:
U
R 1 i+ 2
1 = U i +1 − ( 1 − κ ) ∆ 3 + (1 + κ ) ∆ 1 és i+ i+ 4 2 2
L
U
1 i+ 2
= Ui +
1 ( 1 − κ ) ∆i − 1 + (1 + κ ) ∆i + 1 , 4 2 2
(2.90)
(2.91)
L
R
amelyben az U i + 1 és U i + 1 paraméterek az ( i , j ) és az ( i + 1 , j ) cellák között elhelyezkedő
(i +1 2 , ∆
i+
3 2
2
2
j ) cellahatár jobb és bal oldalán megjelenő értékeket jelentik, továbbá = U i + 2 − U i +1 , ∆ 1 = U i +1 − U i , valamint ∆ 1 = U i − U i −1 . A (2.90) és a (2.91) i+ i− 2
2
kifejezések a magasabbrendű sémák egy általános alakja, amelyet gyakran neveznek a másodrendű sémák κ osztályának. Az egyoldalú másodrendű diszkretizálás a κ = −1 , a Fromm séma κ = 0 , a harmadrendűen pontos (upwind) diszkretizáció κ = 1 3 , a Leonard séma κ = 1 2 , és a három pontos centrális differencia séma κ = 1 esetén áll elő. Sajnos azonban a magasabbrendű upwind eljárások esetén, hasonlóan a centrális sémákhoz, hamis oszcillációk jelennek meg a lökéshullámok közelében, amelyekért a numerikus diszkretizáció egyenértékű egyenletében megjelenő páratlan rendű deriváltak tehetők felelőssé (a párosak okozzák a numerikus disszipációt). A nem oszcilláló megoldások érdekében komplexebb módszereket kellett bevezetni. Két nagy osztálya létezik az ilyen jellegű eljárásoknak; az ENO (Essentially Non Oscillatory), illetve a TVD (Total Variation Diminishing) sémák. Mivel igen pontosan, oszcilláció-mentesen képesek modellezni az áramlást a szakadás környezetében, valamint legalább másodrendűen pontosak azon kívül, ezért gyakran nevezik nagy megoldóképességű módszereknek (high resolution methods). Az ENO sémák esetében a magasabbrendű kiterjesztés miatt nem mindig biztosítható az oszcilláció-mentesség, ezért ebben a munkában a TVD-n alapuló eljárásokat részesítettem előnyben. Az U n paraméter teljes változása (total variation) a következőképpen írható fel:
27
( ) ∑U
TV U n =
∞
i =−∞
A módszer akkor TVD, ha:
(
)
n i +1
− U in .
( )
TV U n +1 ≤ TV U n ,
(2.92)
(2.93)
ami lényegében azt jelenti, hogy nem alakul ki lokális szélsőérték a megoldásban. A monotonitás tulajdonság a TVD-nél erősebb feltétel az oszcilláció-mentesség szempontjából. A TVD, a monotonitás, a pozitivitás tulajdonságokról bővebb információ a [44] irodalomban található. Magasabbrendű upwind TVD sémák A másod, illetve magasabbrendű térbeli approximációk esetén megfelelő korlátok között kell tartani a jobb és bal oldali konzervatív változók gradiensét, hogy ne szűnjön meg a megoldás monotonitása. Ezért a MUSCL módszerek esetén nemlineáris függvényeket, úgynevezett limitereket kell alkalmazni. A legelterjedtebben használt limiterek a Van Albada, a Mulder és a MinMod. Ez utóbbi esetében az extrapolációs eljárás a következőképpen módosul [44]:
U
L i+
U
1 2
R 1 i− 2
= Ui +
1 ( 1 + κ ) ∆i+ 1 + (1 − κ ) ∆i− 1 és 4 2 2
(2.94)
1 (1 − κ ) ∆i+ 1 + (1 + κ ) ∆i − 1 , 4 2 2
(2.95)
=Ui −
ahol:
∆
i+
1 2
= MinMod ∆ 1 , β∆ 1 és ∆ 1 = MinMod ∆ 1 , β ∆ 1 i− i− i+ i+2 2 2 i−2 2
és amelyben a MinMod függvény a következőképpen számítható ki: MinMod ( x , y ) = sgn( x ) max[ 0 , min( x sgn( y ), y sgn( x )) ] . β a kompressziós paraméter: 1 ≤ β ≤
3−κ , amelyet a validációs eljárás során 1 3 értékkel 1− κ
vettem figyelembe. A (2.94), illetve a (2.95) egyenletekben szereplő magasabbrendű interpolációs függvények különböző típusú változókkal kerülhetnek kiszámításra: konzervatív, primitív, illetve karakterisztikus változók formájában. Az 1D-s áramlások tekintetében a karakterisztikus változók egyértelmű előnyt élveznek a primitív és a konzervatív változókkal szemben. Ennek oka, hogy a Riemann invariánsok közül csak az egyik változik meg kismértékben a lökéshullámon, illetve a kontakt-felületen keresztül, ellenben a konzervatív vagy primitív változós esetekben, ahol a paraméterek változása jelentős [44]. 2, illetve 3D-s esetekben szintén a karakterisztikus változók jelentik a legpontosabb közelítést. A számítás gyorsaságának tekintetében azonban a primitív változók előnyt élveznek a konzervatív és a karakterisztikus változókkal szemben, ugyanis ebben az esetben elkerülhető a különböző típusú változók közötti transzformáció. Dadone és Grossmann 1991, valamint Manna 1991 által végzett numerikus vizsgálatok arra engedtek következtetni, hogy a konzervatív változók használata esetén csökkent a konvergencia sebessége, továbbá némely paraméter esetében oszcillációk jelentek meg a szakadások közelében [32]. Ezekből, illetve előzetes 28
összehasonlító elemzésekből következően, a primitív alkalmaztam az interpolációs függvények meghatározásra.
2.3.
változókat
( U = ( ρ ,u , v , p ) )
Negyedrendű Runge−Kutta-módszer az idő iterációra
Az előző fejezetekben a különféle térbeli diszkretizációs módszerek során jutottunk el a megmaradási egyenletek a szemi-diszkretizált alakjáig ((2.43) vagy (2.45) egyenlet):
∂ U j = ℜnj , ∂t
(2.96)
amelyben ℜ j egy tetszőleges j pontbeli n időpillanathoz tartozó reziduumot, vagyis pl. a (2.41) egyenletbeli kontúrintegrál numerikus megfelelőjét jelenti. A közönséges elsőrendű differenciálegyenlet rendszer egyaránt megoldható implicit (pl. Newton linearizáció, relaxáció módszerek) és explicit (pl. haladó Euler, negyedrendű Runge−Kutta-módszerek) integrálás segítségével. Az implicit módszerek nagy előnye a feltétel nélküli stabilitás. Az explicit időintegrálás esetén jelentkező feltételes stabilitás miatt maximálni kell a megengedhető legnagyobb ∆t időlépés értékét. Az időben beállt folyamatok modellezése esetén a végső cél az ℜ reziduum zérus értékre történő beállítása olyan gyorsan amilyen gyorsan lehetséges. Az iterációk során az időben pontos tranziens fizikai jelenségek nem fontosak a számítási eljárás lényegében egy numerikus szállítóeszköznek tekinthető, amely a megoldást az időben beállt folyamat irányába viszi. Ebből a szempontból az implicit módszerek használata a legalkalmasabb, megfelelő korlátozással célszerűen nagyra választva az időlépést. A megfelelő korlátozás a nemlinearitások miatti oszcillációs jelenség megszüntetése érdekében szükséges korlátos időlépést jelent. Az explicit módszerek szintén alkalmasak az időben beállt folyamatok modellezésére a stabilitási feltételnek megfelelő, a maximálisan megengedhető időlépés korlátozásával. Az időben változó, tranziens folyamatok modellezésére az explicit módszerek a legalkalmasabbak. A számítási idő csökkentése érdekében azonban az implicit módszerek is alkalmazhatók a fizikai szempontból megfelelően megválasztott (általában kisebb) időlépés figyelembevételével. A számítógépi kapacitás; a memóriafoglalás, és az iterációs idő, valamint a programozhatóság szempontjából az explicit módszerek előnyt élveznek. Az implicit módszerek esetén a háló csomópontjainak maximális száma, az esetlegesen előálló mátrixműveletek elvégezhetősége miatt korlátozva van, ami a komplex háromdimenziós (esetleg DNS) számításokra jelent korlátot. Elsősorban a kisebb számítógépi kapacitás igény miatt az általam kidolgozott numerikus módszerekben explicit időintegrálást alkalmaztam. Megfelelő megkötések mellett többféle stabil explicit módszer létezik (pl. Haladó Euler, Explicit McCormack prediktor-korrektor, Runge-Kutta-módszerek) az olyan közönséges elsőrendű differenciálegyenletek megoldására, mint például a (2.45) egyenlet: n
N ∂ 1 f = ℜ(U j ,U k ) . Uj =− [ H ] Γ − D ∑ n j , k j , k ∂t Ω j k =1
(2.97)
A lehető legkisebb számítógépi kapacitás, illetve az α k paraméter optimális megválasztása esetén elérhető széles stabilitási intervallum elérésének érdekében a negyedrendű nemlineáris Runge-Kutta-módszert alkalmaztam, amely a következőképpen írható fel: U0 = Un
(
U k = U 0 + α k ∆tℜ U k −1 U n +1 = U m
29
)
k = 1, ..., 4 .
(2.98)
Az eljárás stabilitásvizsgálata nem célja a munkának, mivel erről a szakirodalomban elégséges információ áll rendelkezésre [44]. Általában, az időben negyedrendű pontosság elérése érdekében az α k paraméter a következőképpen állítható elő:
αk =
1 . 4 − k +1
(2.99)
A stabilitási tartomány valós tengely menti legnagyobb kiterjedtsége érdekében az upwind típusú diszkretizációk α k paraméterének meghatározása oly módon történik, hogy maximálja a Courant vagy CFL (Courant-Friedrichs-Levy) számot. Ennek értelmében, az elsőrendű diszkretizációk esetén [44]: α 1 = 0 ,1 , α 2 = 0,26 , α 3 = 0 ,5 és α 4 = 1 , valamint a másodrendű diszkretizációk esetén: α 1 = 0 ,12 , α 2 = 0,25 , α 3 = 0 ,5 és α 4 = 1 .
2.4.
Peremfeltételek
A numerikus eljárások fejlesztése során a peremfeltételek helyes meghatározása elsődleges fontosságú a konvergencia és a pontosság szempontjából, de hatással van a konvergencia sebességére is. A fizikai (elsőfajú) és a numerikus (másodfajú) peremfeltételek száma a karakterisztikák módszerének segítségével kerül meghatározásra. A módszer a hiperbolikus típusú egyenletek azon sajátosságán alapul, miszerint az információterjedés irányának meghatározó szerepe van az áramlás modellezésében. A nem elsőfajú, un. nem fizikai peremfeltételek kiszámítása során a legelterjedtebb eljárások a következők: a karakterisztikák módszere, extrapolációs módszerek, kompatibilitási kapcsolatok módszere, illetve ezek kombinációja. A legnépszerűbbek a különböző rendű extrapolációs módszerek, azonban a karakterisztikák módszerének segítségével a Riemann invariáns állandóvá tételével meghatározott peremen érvényes paraméterek esetén a zavarások intenzitása csökken, ami gyorsabb konvergenciához vezet. A Jacobi-mátrix sajátértékeinek ((2.25) mátrix) a peremen elvégzett lokális egydimenziós analízise során a következő jól ismert eredmények adódnak: 1. Szuperszonikus bemeneti perem: Ha a belépő áramlás cellahatár normál komponensére eső Mach-száma nagyobb, mint egy, akkor mind a négy λ( i ) sajátvektor pozitív. Ebből az következik, hogy mind a négy karakterisztikus görbe befelé mutat, vagyis minden információ kívülről érkezik a számítási tér belseje felé. Négy elsőfajú, vagyis fizikai peremfeltételt kell előírni. 2. Szuperszonikus kimeneti perem: Minden λ( i ) sajátvektor negatív, minden, a peremen érvényes peremfeltételek meghatározásához szükséges információ a számítási tér belsejéből érkezik, ezáltal nincs lehetőség fizikai peremfeltétel megadására. Négy numerikus peremfeltételt kell meghatározni.
3. Szubszonikus bemeneti perem (3. ábra (38. o.)): Ha a belépő áramlás normális irányú Mach-száma kisebb mint egy, de nagyobb mint nulla, akkor három sajátérték pozitív és egy negatív. Három fizikai és egy numerikus peremfeltételt kell előírni. 30
4. Szubszonikus kimeneti perem (4. ábra (39. o.)): Ha a kilépő áramlás normális irányú Mach-száma kisebb mint egy, de nagyobb mint nulla, akkor három sajátérték negatív és egy pozitív. Három numerikus és egy fizikai peremfeltételt kell alkalmazni. 5. Szilárd fal perem: Szilárd fal perem esetén a számítási tér és a perem között nem engedhető meg tömeg-, impulzus- és energiatranszport. Így, a cellahatárra merőleges sebességkomponens értékének nullának kell lennie, aminek értelmében egy pozitív sajátérték adódik. Egy fizikai és három numerikus peremfeltételt kell megadni. Az áramlástani gépekben való alkalmazások miatt csak belső áramlásokkal foglalkoztam, ezért a távoli perem (far-field) feltételt nem specifikáltam a peremfeltétel alfejezetben.
2.4.1. Deconinck módszere a peremfeltételek magadására A módszert eredetileg Chakravarthy [9] fejlesztette ki, azonban Deconinck [17] dolgozta át és terjesztette ki két, illetve háromdimenziós alkalmazásokra. A módszerben a perem és a számítási tér közötti konzisztencia a peremen, illetve a cellahatáron előállított numerikus fluxusfüggvények egyenlővé tételével biztosítható. A peremen érvényes U B vektor ismeretlen paraméterei a fizikai és a numerikus peremfeltételek segítségével felírt [44]: ~ (2.100) H n* U B = H n U B ,U I .
( )
(
)
egyenlet segítségével határozhatók meg. A csillaggal ( * ) jelölt érték a peremfeltétel megkötéseit tartalmazó fluxusfüggvényt jelenti. A felső indexben szereplő B (Boundary: perem), illetve I (Inner: belső) betűk a peremen valamint a számítási tér belső pontjában érvényes, megadott, vagy ismeretlen paramétereket jelentik. Az explicit módszerek esetén a (2.100) egyenlet egy nemlineáris egyenletrendszerre vezet. A módszert a Van Leer fluxus megosztáson alapuló módszerénél alkalmaztam. Szuperszonikus bemeneti perem meghatározása A peremfeltételek alfejezet bevezetésében leírtak alapján, a szuperszonikus bemeneti perem esetén négy fizikai peremfeltételt kell előírni [44]: p1 T to p2 α P= = , p3 p p p to 4
(2.101)
ahol T to a torlóponti hőmérséklet, α = arctan( u / v ) a belépő áramlás iránya, p a statikus nyomás és p to a torlóponti nyomás. A P vektor a következő egyenletek segítségével alakítható át a konzervatív ( U = [ ρ , ρu , ρv , ρE ] T ) változók vektorává: T =T
to
(p
to
/p
)
1−γ γ
, M=
( T to / T −1 ) 2 γ −1
,V=
v = V sin( α ) , ρ = p / ( RT ) és ρE = 31
M , u = V cos( α ) , γ RT p V2 +ρ . γ −1 2
Szubszonikus bemeneti perem meghatározása Szintén a peremfeltételek alfejezet bevezetésében leírtak alapján, a szubszonikus bemeneti perem esetén három fizikai peremfeltételt kell előírni, amelyek rendre a torlóponti hőmérséklet T to , az áramlás iránya α , illetve a torlóponti nyomás p to . Ekkor a peremen érvényes állapotváltozók vektora [44]: to p1 T p2 M P= = , p3 α p to 4 p
(2.102)
amelyek közül a Mach-szám a számítási tér függvényeként állítható elő. A fluxusvektor megosztás (2.60) egyenlete a (2.100) egyenlet figyelembevételével a következőképpen írható fel: ~ H n U I ,U B = H n+ U I + H n− U B = H n* U B = H n+ U B + H n− U B , (2.103)
(
)
( )
( )
vagyis:
( )
( )
( )
( )
( )
H n+ U I = H n+ U B ,
(2.104)
ρc 2 ρc 2 4 (1 + M n ) = 4 ( 1 + M n ) . I B
(2.105)
illetve:
A (2.105) egyenlet bal oldala ismert, mivel a számítási tér állapotváltozóit tartalmazza, ellenben a jobb oldallal, ahol peremen érvényes paramétereket kell meghatározni. Ennek érdekében az energia- és a Poisson-egyenlet segítségével, a (2.105) egyenlet peremen érvényes állapotváltozói felírhatók a Mach-szám függvényeként [44]: 2
a 1 + α n M 1 + κ − 1 M 2 − k = 0 = f ( M ) , 2 αm
(2.106)
ahol: 2
v v α m = 1 + , α n = n x + n y , a = −( γ + 1) / ( 2γ − 2 ) és u u
[
k = ρc ( 1 + M n )
2
] /[p
to
I
λ / RT0
]
B
.
A (2.106) egyenlet változói az M = M B peremen érvényes Mach-szám kivételével ismertek, mivel vagy fizikai peremfeltételként adottak vagy a belső számítási térhez tartozó értékek. A (2.106) egyenlet nemlineáris egyenlet a peremen érvényes Mach-számra, amit az egyenlet könnyen elvégezhető differenciálása miatt célszerű Newton-Raphson módszerrel megoldani. Ha a Mach-szám ismert a sűrűség: ρ , a statikus hőmérséklet: T és nyomás: p analitikus összefüggések segítségével határozhatók meg. Végezetül a konzervatív változók U B vektora szintén könnyen előállítható. Szuperszonikus kilépő perem meghatározása
32
A Jacobi-mátrix minden sajátértéke negatív, ezért nem lehet fizikai peremfeltételt előírni. Négy numerikus (másodfajú) peremfeltétel alkalmazásával extrapolációs technikával (2.4.2. alfejezet) határozhatók meg a peremen érvényes konzervatív áramlástani paraméterek [44]. Szubszonikus kilépő perem meghatározása Mint ismeretes, ebben az esetben egy fizikai és három numerikusan meghatározott peremfeltételt kell megadni. A numerikus irodalomban a statikus kilépő nyomás fizikai peremként való megadása a legelterjedtebb a peremen érvényes primitív változós vektorban [44]: p1 ρ p2 u P = = . p v 3 * p p 4
(2.107)
A peremfeltétel kiszámítása számítástechnikailag könnyebben elvégezhető a fluxusvektorok cellahatár normális és tangenciális irányába eső elforgatott rendszerében, mivel ekkor a tangenciális irányba eső összetevő értéke nulla. Ez egyben a V sebesség vektor e x , e y koordináta rendszerből az n , s rendszerbe történő transzformálását is jelenti. Ekkor a (2.61) egyenletbeli fluxusfüggvény megosztás a következőképpen írható fel: 1 ( γ − 1) V ± 2 c n γ γ ρc 2 ± Hn = ± (1 ± M n ) V s 4 [ (γ − 1)V ± 2c] 2 V 2 n + s 2 2 2 γ −1
(
)
.
(2.108)
A (2.108) egyenlet második komponense az adott p * nyomás, a (2.103) egyenlet második egyenlősége és a (2.40) egyenlet segítségével a következőképpen írható fel: ρc 2c 2 ( γ − 1) Vn + − (1 + M n ) γ I γ 4 ρc ( γ − 1) 2c − (1 − M n ) 2 Vn − = ρVn2 + p* γ B γ 4
[
]
B
.
(2.109)
A (2.109) kifejezés egy nemlineáris egyenletre vezet, amely legegyszerűbben NewtonRaphson módszer segítségével oldható meg az ismeretlen, peremen érvényes M n Machszámra: a1 M n3 + a 2 M n2 + a 3 M n + a 4 = 0 , ahol a1 = γ − 1 , a 2 = 2γ , a 3 = γ + 3 és 2c 2 ( γ − 1) a 4 = 2 − ρc( 1 + M n ) Vn + / p * . γ I γ 33
(2.110)
Ha az M n Mach-szám ismert, a peremen érvényes hangsebesség a (2.108) egyenlet első komponensének, vagyis a (2.105) egyenlet segítségével számítható:
[c] B
[γp ( M = [ ρc( M *
n n
+ 1)
+ 1)
2 2
] ]
B
.
(2.111)
I
Ezek után a normális irányú sebesség és a sűrűség könnyen meghatározható. A (2.108) egyenlet harmadik komponense a (2.105) egyenlet figyelembevételével a következőképpen írható fel:
[Vs ] B = [Vs ] I .
(2.112) Végezetül a n , s térből a Descartes-féle, e x , e y térbe történő inverz transzformáció segítségével az összes keresett, a peremen érvényes áramlástani paraméter meghatározható [44]. Szilárd fal perem meghatározása A peremfeltételek alfejezet bevezetésében leírtak alapján, szilárd fal perem esetén nem engedhető meg tömeg-, impulzus- és energiatranszport a perem és a számítási tér között. Ennek értelmében a keresett numerikus fluxusfüggvény a következőképpen írható fel [44]: 0 pn x * Hn = . pn y 0
(2.113)
Hasonlóan a szubszonikus kilépő perem meghatározásához a (2.108) egyenlet harmadik komponensét, a (2.112) egyenletet érdemes felírni:
[Vs ] B = [Vs ] I .
(2.114)
Az egyszerűség kedvéért a (2.105) egyenlet a számítási térből ismert bal oldalát jelölje H 1 : ρc 2 + 4 ( 1 + M n ) = H1 ( U I ) = H 1 . I
(2.115)
Ekkor a (2.108) egyenlet első komponense a (2.113) és a (2.103) egyenletek figyelembevételével a következőképpen írható fel: ρc 2 4 (1− M n ) = H 1 . B
(2.116)
A (2.108) egyenlet negyedik komponensét szintén a (2.113) és a (2.103) egyenletek segítségével a következőképpen lehet felírni: 2 Vs 2 2 [ ( γ − 1)V n + 2c] ρ c ( 1 + M ) + = n 2 2 γ 2 −1 I
(
)
2 Vs2 2 [ ( γ − 1)Vn − 2c ] = ρc(1 − M n ) + . 2 2 γ 2 −1 B
(
)
(2.117)
34
A (2.114) egyenlet (2.117) egyenletbe való helyettesítését követően kapjuk
{
}
H 1 [ ( γ − 1)Vn + 2c ] I − [ ( γ − 1)Vn − 2c ] B = 0 . 2
2
(2.118)
Mivel H 1 ≠ 0 , ezért a (2.118) egyenletnek 2 megoldása van. Az első
[ ( γ − 1)Vn + 2c ] I = [ ( γ − 1)Vn − 2c ] B ,
(2.119)
amely a (2.108) egyenlet második komponensébe való helyettesítés után, a (2.103) és a (2.113) egyenletek segítségével a következőképpen írható fel: ρc 2c 2 ( γ − 1) Vn + − (1 + M n ) γ I γ 4 ρc ( γ − 1) 2c − (1 − M n ) 2 Vn − = [ p ] B , γ B γ 4
(2.120)
ami zérusértékű fali nyomásra vezet, ezért ez a megoldás fizikai szempontból nem elfogadható. A (2. 118) egyenlet második megoldása
[ ( γ − 1)Vn + 2c] I
= −[ ( γ − 1)Vn − 2c ] B ,
(2.121)
amely γ -val való osztás után
H 2 2c γ − 1 = − Vn H1 γ γ B
(2.122)
ρc ( γ − 1) 2c H 2 = (1 + M n ) 2 Vn + . γ I γ 4
(2.123)
alakban írható fel, amelyben:
A (2.121) egyenlet γ -val való osztása és a (2.120) egyenletbe való helyettesítése után (a (2.115) egyenlet figyelembevételét követően) felírható
[ p] B
= 2H 2 .
(2.124)
A peremen érvényes hangsebesség a (2.122), a sűrűség és a c = γp / ρ segítségével határozható meg c=
H2 γ [ 2 − ( γ − 1) M n ] H 1 és
egyenletek
(2.125)
2 2 2 H1 ρ = [ 2 − ( γ − 1) M n ] . γ H2
(2.126)
A (2.125) és a (2.126) egyenletek beírva a (2.116) egyenletbe, a következő összefüggés írható fel a szilárd fal peremen érvényes normális irányú Mach-számra:
( 1 − γ ) M n3 + 2γM n2 − ( 3 + γ ) M n
=0.
(2.127)
A levezetés során zérusértékű normális irányú sebességkomponens nem volt feltétel, azonban egyszerű behelyettesítés után azonnal megállapítható, hogy csak az M n = 0 gyök az egyenlet megoldása. A normális irányú Mach-szám ismeretében, a (2.125), a (2.126), a (2.124), és a (2.114) egyenletek segítségével minden áramlástani paraméter meghatározható [44].
35
2.4.2. Extrapolációs módszer a peremfeltételek magadására A legelterjedtebb eljárás a másodfajú peremfeltételek meghatározására az extrapolációs módszer, amely lehet nulla, első- vagy magasabbrendű. Gustafsson és Sundstrom (1978) megmutatta, hogy lineáris egyenletek esetén, a másodfajú peremek meghatározására alkalmazott módszer pontosságának legalább egy renddel alacsonyabbnak kell lennie, mint a magának a numerikus eljárásnak, hogy ne csökkenjen a globális számítás pontosságának rendje [28]. Ebből következően a centrális diszkretizáció mesterséges viszkozitással figyelembe vett numerikus eljárás kidolgozásakor a numerikus (másodfajú) peremértékek meghatározására lineáris extrapolációt alkalmaztam. A fizikai (elsőfajú) és numerikus peremfeltételek számának megadása a peremfeltételek alfejezet bevezetésében leírtak szerint történt. Szuperszonikus belépő perem esetén a négy fizikai peremfeltétel a következő volt: torlóponti nyomás, torlóponti hőmérséklet, állásszög és a statikus nyomás. Szubszonikus belépő perem esetén az előző peremfeltétel megadás azzal a különbséggel igaz, hogy a statikus nyomás mint másodfajú perem kerül meghatározásra lineáris extrapoláció segítségével. Szuperszonikus kilépő perem esetén a négy konzervatív változót: U = [ ρ , ρu , ρv , ρE ] T másodfajú peremfeltételként, lineáris extrapoláció segítségével határoztam meg. Szubszonikus kilépő perem esetén a statikus nyomást írtam elő, mint fizikai peremfeltétel, míg a peremen érvényes három konzervatív változót, U = [ ρ , ρu , ρv ] T lineáris extrapoláció segítségével határoztam meg. A szilárd fal perem esetén szimmetria (tükör) elsőfajú peremfeltételt alkalmaztam, vagyis a perem normális irányú sebességkomponensét ellenkező előjellel egyenlővé tettem a számítási térbeliével. Ennek következtében nem alakult ki anyag, impulzus és energia transzport a perem és számítási tér között. A numerikus perem meghatározására szintén lineáris extrapolációt alkalmaztam. A peremfeltételek megadása során figyelmen kívül hagyott paramétereket izentrópikus összefüggések segítségével határoztam meg.
2.4.3. Karakterisztikák módszere a peremfeltételek magadására A karakterisztikák módszere figyelembe veszi a hiperbolikus típusú parciális differenciálegyenletekre jellemző információterjedés kitüntetett irányát, így az alkalmazásával meghatározott peremfeltétel fizikai szempontból a legkorrektebb. Újfajta módszert dolgoztam ki a cella-normális irányába levetített áramlástani paraméterek egydimenziós karakterisztikus peremfeltételének kiszámítására, amelyet a Roe által közelített Riemann megoldóban, a szubszonikus be-, és kilépés esetén alkalmaztam. Minden más esetben extrapolációs módszert használtam a numerikus peremek meghatározására. Szubszonikus bemeneti perem meghatározása
36
A peremfeltételek alfejezet bevezetésében leírtak alapján, a szubszonikus bemeneti perem esetén három fizikai peremfeltételt kell előírni, amelyek rendre a torlóponti hőmérséklet T to , az áramlás iránya α , illetve a torlóponti nyomás p to . A peremen érvényes állapotváltozók vektora: p1 T to p2 M P= = . p 3 αto p4 p
(2.128)
Számítástechnikailag könnyebben meghatározhatók a numerikus peremértékek, ha a fluxusokat felbontjuk a cellahatárral párhuzamos, illetve rá merőleges irányba, ami a sebességvektorok elforgatásának segítségével valósítható meg. Ez lényegében egy lokálisan egydimenziós információterjedésnek felel meg a cellahatárra merőleges irányba, míg a cellahatárral párhuzamos komponens nincs hatással az információterjedésre. A karakterisztikák módszere szerint a Riemann invariánsok állandóak a hozzájuk tartozó karakterisztikus görbék mentén ((2.30) egyenlet alapján): ∂Wn(1) ∂Wn( 2 ) ∂Wn = ( 3) ∂Wn ∂W ( 4 ) n
∂ρ − 1 ∂p = 0 , c2 ∂Vs = ∂ V + 1 ∂ p = 0, n ρc 1 ∂ p = 0, − ∂Vn + ρc
dn = Vn dt dn = Vn dt dn = Vn + c dt dn = Vn − c dt
görbén görbén , görbén görbén
(2.129)
ahol ∂V s = s ∂V = ( n y ;− n x )( ∂u ; ∂v ) . Izentrópikus feltételezéssel élve, a Riemann vagy karakterisztikus változók integrálása után a perem ( B (Boundary)), illetve a számítási tér ( I (Interior)) között a következő összefüggések írhatók fel [28]:
pI p = Bγ , γ ρI ρB
(2.130)
Wn( 2 ) = Vs ,I = Vs ,B ,
(2.131)
Wn(1) =
Wn( 3) = Vn , I +
2c I 2c = Vn ,B + B és γ −1 γ −1
(2.132)
2cI 2c = Vn , B − B . γ −1 γ −1
(2.133)
Wn( 4 ) = Vn , I −
Jelen esetben a (2.133) egyenlet érdemes figyelembe venni, mivel ez tartozik a Jacobi-mátrix negatív λ(n4 ) = Vn − c sajátértékéhez, ami a számítási térből kifelé történő információterjedést reprezentálja (3. ábra). Következő lépésként meghatároztam a Mach-szám és normális komponense közötti összefüggést:
t λ (3) (+ ) V+c t2 n ( 2 ) n λ n(+ ) Vn λ (4)n(− ) Vn-c λ (1) (+ ) V nn
37
v 2 M 2 = M x2 + M 2y = M x2 1 + , u
3. ábra.
(2.134)
Információterjedés szub-szonikus belépő perem esetén a karakterisztikák segítségével
belépõ celahtár
n . UI t1 UB áramlásirány
egyenletek segítségével a sűrűség és a helyi hangsebesség a következőképpen fejezhető ki a Mach-szám függvényeként: p to γ − 1 2 ρ= M 1 + 2 RT to
−
1 γ −1
−
1
2 , illetve c = γRT 1 + γ − 1 M 2 . 2 to
(2.139)
A (2.137) második, a Vn = M n c és a (2.139) egyenletek felhasználásával a (2.133) egyenlet, mint a peremen érvényes ismeretlen Mach-szám függvényeként a következőképpen írható fel: γ −1 − 2 γ −1 2 a( M B− 2 + ) − b1 + MB 2 2 1
−
1 2
−k = f ( M B) = 0 ,
(2.140)
ahol a=
α n γRTBto αm
, b=
2cI 2 γRTBto és k = Vn , I − . γ −1 γ −1
(2.141)
Az áramlási tér belső részéből ismert paraméterek, illetve a megadott elsőfajú peremfeltételek miatt ((2.141) egyenletek), a (2.140) egyenlet a nemlineáris formában megjelenő M B Machszámra megoldható. Az egyenlet differenciálása könnyen elvégezhető, ezért a (2.140) egyenletet Newton-Raphson módszer segítségével oldottam meg, amíg az M B = M Búj feltétel nem teljesült: M =M új B
régi B
(
)
f M Brégi − ' . f M Brégi
(
)
(2.142)
A Mach-szám ismeretében a peremen érvényes sűrűség a (2.139) első, a statikus hőmérséklet a (2.138) első, a statikus nyomás a (2.138) utolsó és a hangsebesség pedig a (2.138) harmadik egyenletének segítségével határozható meg. A fizikai peremértékként megadott belépő áramlás szög ( α ), és a V = Mc ismeretében az u és v sebességkomponensek, végül az U B peremen érvényes konzervatív változók vektora is könnyen előállítható. Szubszonikus kimeneti perem meghatározása A peremfeltételek alfejezet bevezetésében leírtak alapján, a szubszonikus kimeneti perem esetén egy fizikai, a statikus nyomás és három numerikus peremfeltételt kell előírni, illetve maghatározni (4. ábra). Legyen a peremen érvényes primitív változók vektora a következő: p1 ρ p u P = 2 = . p v 3 p4 p
(2.143)
Hasonlóan az alfejezetbeli szubszonikus bemeneti peremhez, először leképeztem a numerikus fluxusfüggvényeket a cellahatárra merőleges irányba. A karakterisztikák módszere alapján a (2.130) egyenlet ismét felírható a
Wn(1 ) =
pI p = Bγ γ ρI ρB
(2.144)
38
egyenlet a Vn karakterisztikán a perem ( B ) és a számítási tér ( I ) között, amiből a peremen előálló sűrűség meghatározható. A statikus hőmérséklet a gáztörvény segítségével számítható ki. A kimeneti cellahatárral párhuzamos sebességkomponens nem játszik szerepet a cellahatárra merőleges irányú információterjedésben, ezért előírtam, hogy a második, szintén
t
λ (3)n ( + ) Vn+c
t2
λ ( 4 ) n( − ) Vn-c
λ ( 2 ) n ( + ) Vn λ
s
UI
.
. n t1 k ilé põ ce lla ny
á lásir
áram
39
(1 )
UB ha t
ár
n(+
) Vn
4. ábra. Információterjedés szubszonikus kilépő perem esetén a karakteriszti-kák segítségével
Vn karakterisztikán érvényes Riemann invariáns egyezzen meg számítási tér és a perem között: Vs ,I = Vs ,B .
(2.145)
cB = γRTB
(2.146)
A peremen érvényes hangsebesség a
összefüggés segítségével határozható meg. Végezetül a Vn + c karakterisztikán állandó harmadik Riemann invariáns Wn( 3) = Vn , I +
2c I 2c = Vn ,B + B γ −1 γ −1
(2.147)
segítségével számítottam ki a peremen érvényes, a cellahatár normális irányú sebességkomponensének nagyságát. A vektor jellegű mennyiségek Descartes rendszerbe történő visszatranszformálásával a (2.143), illetve a (2.2) egyenletbeli primitív és konzervatív változók vektora szintén könnyen előállítható. Az újfajta peremfeltétel számítási eljárás validálását a 2.5. alfejezetben, a Roe által közelített Riemann módszerrel együtt végeztem el. Az új karakterisztikus és az extrapolációs peremfeltételek összehasonlítása érdekében a 2.5.1. alfejezetben ismertetett Laval-cső st 3 geometrián a p = 1140000 Pa , a ρ = 8,09 kg / m , az u = 111,49 m / s , a v = 0,0 m / s és a T st = 490,65o K kezdeti feltételekből kiindulva futtattam le mindkét peremfeltétel számítási eljárásra külön-külön a programot. Az első fajú peremfeltételek mindkét esetben a következők voltak: α in : lineáris eloszlás az alsó és a fölső kontúr iránytangense között, a st p into = 1200000 Pa , és a Tinto = 500,15 o K bemeneti és a pout = 956165,6 Pa kimeneti paraméterek (5/a. ábra).
A disszertációban előforduló numerikus számításokat a Budapesti Műszaki Egyetem EISZK üzemeltetése alá tartozó szuperszámítógépen végeztem el. A Compaq szuperszámítógép egy 4 node-os "serverfarm", amely összesen 16 darab EV5.6 (21164A, 600 Mhz, 8 MByte cache) Alpha CPU-val, 32 GByte memoriával és 0.62 TByte lemezzel rendelkezik. A node-ok Memory Channel-en (Full/100 Mb) keresztül tartják egymással a kapcsolatot. A futtatás eredményeit tekintve a karakterisztikus peremfeltétel alkalmazásának esetén a számítás időszükséglete 50,78 %-a volt az extrapolációs peremfeltétel alkalmazásnak. A konvergencia lefutások során a karakterisztikus peremfeltétel számítási eljárás kevesebb iterációs ciklus alatt érte el a konvergencia-kritérium értékét (5/b. ábra). Az összenyomható áramlások konvergencia kritériumára a sűrűség reziduum normalizált L2 normáját használtam: ∆ρ 1 = log10 ρ Np
2
∆ ρi , ∑ i =1 ρ i Np
(2.148)
amelyben N p a háló pontjainak száma és a ∆ kifejezés két egymást követő idő iteráció közötti − jelen esetben sűrűség − változásra utal. A konvergencia kritérium értékét (-7) minden esetben úgy határoztam meg, hogy a további iterációk ne jelentsenek lényegi változást az áramlástani paraméterek tekintetében.
40
5/a. ábra. Mach-szám eloszlás a Laval-csőben a karakterisztikus peremfeltételekkel figyelembe vett Roe által közelített másodrendű Riemann-féle módszerrel Iteráció szám 0
2500
5000
7500
10000
Sűrűség reziduum L 2 normája
-2 -3
Extrapolációs peremfeltétel
-4
Karakterisztikus peremfeltétel
-5 -6 -7
5/b.
2.5.
ábra. Karakterisztikus és extrapolációs peremfeltételek összehasonlítása az iteráció szám tekintetében a Laval-csőben történő áramlásra
Validáció
A különféle numerikus módszerek és a peremfeltételeik meghatározásához kidolgozott számítási eljárások jóságának megállapítása alapvetően háromféle módszerrel valósítható meg. A legjobb a kísérletekkel történő érvényesítés, azonban a mérőelemek áramlásra gyakorolt hatása, valamint a mérési körülmények előállításának nehézsége vagy a megfelelő mérőberendezések (pl. szuperszonikus szélcsatorna) hiánya miatt igen elterjedt az eredmények analitikai számítással (pl. Laval-cső) vagy egy validációs geometriai modell (pl. körprofilt magába foglaló csatorna) segítségével − kereskedelmi vagy egyéb programok eredményeivel − történő összehasonlítása, validálása. A [69] szakirodalomból rendelkezésemre álló Schlieren felvétel mintájára elkészítettem egy DCA ((Double Circular Arc) kettős körívből álló profil) lapátrács 2D-s geometriáját. Az ismert belépő Mach-szám (M1=1,1) előállítása során, a fényképen látható lökéshullám mintázat igen jó egyezőséget mutat mindkét számítási eljárás (centrális diszkretizáció Jameson-féle mesterséges disszipációs módszerrel (2.2.1.1. alfejezet) és Roe által közelített másodrendű Riemann-féle módszer (2.2.1.2./b. alfejezet) MinMod limiterrel (2.2.2. alfejezet) és a karakterisztikus peremfeltételekkel (2.4.3. alfejezet)) állandó Mach-szám görbéivel (6. ábra). Fontos azonban megjegyezni, hogy az Euler alapú numerikus módszerek nem veszik figyelembe a határréteget, így a valóságban a lapátozás után a határréteg okozta csatorna szűkülés miatt nagyobb sebességre gyorsul fel a közeg.
41
M1=1,1
M1=1,1 1
6. ábra. Fent: Lökéshullám mintázat DCA lapátrácsban (Schlieren felvétel [69]). Állandó Mach-szám vonalak DCA lapátrácsban; középen: centrális diszkretizáció Jameson mesterséges disszipációs módszerével, lent: Roe által közelített másodrendű Riemann-féle módszerrel, MinMod limiterrel és a karakterisztikus peremfeltételekkel
2.5.1.
Laval-cső 42
Ebben a teszt esetben a csatorna A = A( x ) meridionális profilját két fő részre osztottam. Az elsőt, a belépő keresztmetszettől a torokig az
x A( x ) = 1 + 1.51 − 5
2
(2.149)
7. ábra. 2D Laval-cső geometria a numerikus hálóval függvénnyel, míg a másodikat a toroktól a kilépő keresztmetszetig az
A( x ) = 1 + 0.5 1 −
x 5
2
(2.150)
függvénnyel valósítottam meg. A csatorna hossza 10 méter. A 2D numerikus modell esetén a hálópontok száma i irányba 100 és j irányba 15 (7. ábra). A numerikus teszteset futtatása során a következő bemeneti és kimeneti elsőrendű peremfeltételeket írtam elő: α in =lineáris eloszlás az alsó és a fölső kontúr iránytangense to st között, p in = 1200000 Pa , Tinto = 500,15 o K , és p out = 724214,2 Pa . Ennek eredményeként, a csatorna negyedik ötödrészénél egy lökéshullámnak kell megjelennie, amit az analitikai és a 3 numerikus (centrális diszkretizáció Jameson mesterséges disszipációs módszer (2.2.1.1. alfejezet), Roe által közelített másodrendű Riemann módszer (2.2.1.2./b. alfejezet) MinMod limiterrel (2.2.2. alfejezet) és a karakterisztikus peremfeltételekkel (2.4.3. alfejezet), valamint a másod rendű Van Leer fluxusvektor megosztáson alapuló módszere (2.2.1.2./a. alfejezet)
1,7 Mach szám
1,4 1,1
Egzakt eredmény
0,8
Numerikus eredmény
0,5 0,2 0
2
4
6
8
10
Hossz (m)
8. ábra. Jameson mesterséges disszipációjával figyelembe vett centrális diszkretizációs módszer Laval-csőben történő validálása szintén a MinMod limiterrel) számítási eredmény igazol (8-10. ábra). A numerikus módszerek eredményei a Laval-cső középvonalán elhelyezkedő paraméterek ábrázolásán keresztül igen jó 43
1,7
Mach szám
egyezőséget mutat az analitikus módszer eredményivel mindhárom esetben. Érdemes 1,4 hogy a mesterséges disszipációs módszer esetén, a túlzott disszipációnak megjegyezni, köszönhetően a megoldás nem követi a lökéshullámot olyan mértékben, mint a nagy 1,1 Egzakt eredmény pontosságú (high resolution) upwind módszerek esetén. Ezzel magyarázható az MUSCL módszerek alkalmazásának előnye nagyobb Mach-számú áramlásokNumerikus esetén. eredmény 0,8 0,5 0,2 0
2
4
6
8
10
Hossz [m]
9. ábra.
A karakterisztikus peremfeltételekkel figyelembe vett Roe által közelített másodrendű Riemann-féle módszer Laval-csőben történő validálása
2.5.2. Körprofilt magába foglaló 2D-s csatorna A körprofilt magába foglaló 2D-s csatorna alap teszt esetet eredetileg Rizzi és Viviand 1,7
Mach szám
1,4 Egzakt eredmény
1,1
Numerikus eredmény
0,8 0,5 0,2 0
10.
2
4 6 Hossz [m]
8
10
ábra. Másodrendű Van Leer-féle fluxusvektor megosztáson alapuló módszer Laval-csőben történő validálása
[54] dolgozta ki transzszonikus áramlások modellezésére. A négyzet alakú csatorna alsó falának középső részén helyezkedik el a 4,2 % maximális vastagságú körszelet profil (11. ábra). A számítási tér 3 különböző részre van felosztva. Az első részben, x ∈ 3,00, 4 ,35 , a
[
[
]
] közötti vízszintes irányú távolság. Végezetül az utolsó részben, x ∈ [5,65, 8] , a hálózás ritkul
háló a profil belépő éléig sűrűsödik. A középső részben, x ∈ 4,35, 5,65 , állandó pontok a profil kilépő élétől a csatorna kimenetéig. A hálópontok megoszlása a következő:
44 11. ábra. Körprofilt magába foglaló 2D-s csatorna a 142x40-s numerikus hálóval
amely az adott esetén nemlineáris függvénykapcsolatot jelent 2D-s N és Állandó L paraméterekre R 13. ábra. Mach-szám vonalak a körprofilt magába foglaló értékére ∆ξ M ( ∆ξ icsatornában legnagyobb megkívánt értéke) ismeretében. A keresett paraméter Newtona 142x40-s hálón. Folytonos vonal: Roe által közelített Riemann-féle módszer ameg karakterisztikus peremfeltételek Raphson eljárássalmásodrendű vagy fixpontos iterációval határozható (11. ábra). y irányban állandó alkalmazásával, szaggatott vonal: Fluent osztásközöket vettem fel. Az összenyomható, transzszonikus áramlások modellezése esetén a be-, illetve a kimeneti peremfeltételek segítségével 0,85 Mach-számot állítottam elő a csatorna profil nélküli be- és a kilépő részénél. Ennek eredményként a profil hátsó része fölött megjelent egy gyenge lökéshullám, mialatt az áramlás kb. 1,33 Mach-számig expandál. Az áramlás numerikus szimulációi után a 3 numerikus módszer (a centrális diszkretizáció Jameson-féle mesterséges disszipációs módszerrel (2.2.1.1. alfejezet), a Roe által közelített másodrendű Riemann-féle eljárás (2.2.1.2./b. alfejezet), valamint a Van Leer-féle fluxusvektor megosztáson alapuló másodrendű módszer (2.2.1.2./a. alfejezet) szintén a MinMod limiterrel (2.2.2. alfejezet)) számítási eredményei a 12-14. ábrán láthatók. A széleskörűbb összehasonlítás érdekében elvégeztem a modell Fluent kereskedelmi programmal történő analízisét is, amelynek eredményeit szaggatott vonallal ábrázoltam (12-14. ábra). A geometria, a hálóméret, a peremfeltételek és az ideális áramlásra vonatkozó feltételek minden esetben megegyeztek. A minimális eltérés oka a különféle numerikus módszerekben rejlő eltérő közelítési elvből adódhat. A mesterséges disszipációs módszer esetén a lökéshullám megoldása alulmarad a nagy pontosságú upwind módszerekhez képest. A saját és a Fluent programok Mach-szám eredményei közötti eltérés számszerűsítésére a ∆M i , j ∆M i , j = ABS Mi,j összefüggést, vagyis a cellánkénti normált eltérés abszolút értékének kifejezését használtam, amelynek eredményei a 15. ábrán láthatók. Az eltérés minimális, a lökéshullám környezetében tapasztalható a legjelentősebb különbség. A különféle eljárások tekintetében a legnagyobb differencia a centrális módszer, míg a legjobb egyezőség Roe-féle módszer esetén adódott. További összehasonlító elemzéseket végeztem a saját és a Fluent programok eredményei közötti eltérés vizsgálatára a Mach-szám átlagos abszolút normált eltérése: ∆M =
N ∆M i 1 p ABS ∑ N p i =1 Mi
és a Mach-szám normált L2 normája: 2
N ∆M 1 p ∆M i = log 10 ∑ M N p i =1 M i
segítségével. Az 1. táblázatba foglalt eredmények jól tükrözik a validálás során megállapított következtetéseket, miszerint a legnagyobb eltérés a mesterséges disszipációs módszer esetén adódott a Fluent program eredményével történő összehasonlításkor. Az upwind módszerek esetén a Fluent eredményeihez viszonyított eltérés jelentősen kisebb (1 táblázat). Az összenyomhatatlan áramlások modellezésének validálására is felhasználtam a körprofilt magába foglaló 2D-s teszt esetet (3.1.1. alfejezet). Ebben az esetben a következő to bemeneti és kimeneti elsőfajú peremfeltételeket írtam elő: α in = 0° , pin = 100000 Pa és st p out = 98000 Pa .
45
12. ábr
14. ábr
EltérésJameson-FluentRoe-Fluent Van Leer-Fluent 0,008840,00421 0,00364 - 1,648 -2,095 -2,0044
1. táblázat
. Eltérés a CFD programok Mach-szám eredményei között
46
15. ábra. A cellánkénti Mach-szám normált eltérésének abszolút értéke: fent: centrális diszkretitáció Jameson-féle disszipációval és Fluent, középen: másodrendű Roe-féle módszer és Fluent, lent Van Leer másodrendű módszere és Fluent eredményei között
2.6.
2D-s transzszonikus rotor karakterisztika meghatározása
Napjainkban, a transzszonikus axiál-kompresszorokat gyakran alkalmazzák repülőgépeken használt gázturbinás sugárhajtóművekben [62]. Mérnöki szempontból a fő cél a súly, a méret és a költségcsökkentés, mialatt nem változnak vagy javulnak az olyan tervezési paramétereket, mint például a nyomásviszony vagy a teljesítmény felvétel. Többfokozatú kompresszorok esetén ez elsősorban a fokozatszám csökkentéssel valósítható meg, azonban a gépegység nyomásviszonyának állandón tartása, esetleg növelése érdekében a fokozatok tervezési paraméterein kell javítani. A fokozatokban lezajlódó nyomásnövekedés nemcsak a lapátterhelés, hanem a fordulatszám növelésével is javítható, aminek következményeként a relatív belépő sebesség a transz-, vagy a szuperszonikus tartományba kerülhet. A megjelenő gyenge lökéshullámokon keresztüli nyomásnövekedés jelentősen javíthat a futólapátozás torlóponti nyomásviszonyán. Az intenzív lökéshullámok esetén azonban a keletkező veszteség, pl. a lökéshullám mögött hirtelen megnövekedő nyomás esetén előálló leválás miatt nem garantálható az elvárt nyomásnövekedés. Napjainkban, a hő- és áramlástani gépek analízise és tervezése során a 3D-s Euler/NavierStokes-egyenletek numerikus megoldásainak alkalmazása igen elterjedt. Sokféle kereskedelmi (pl. FLUENT, CFXTascflow, STAR-CD, FLOW3D, ADINA-F, COSMIC NASA, ACRi) és saját fejlesztésű numerikus kód áll a felhasználók rendelkezésére. A korszerű, összenyomható áramlást modellező nagy pontosságú (high resolution upwind methods) numerikus eljárások kifejlesztésével és a kompresszorokban lezajlódó folyamatok vizsgálatával többen is foglalkoztak. Daiguji és Yamamoto [15], valamint Matsuo [45], Chakravarthy és Osher [10] nagy pontosságú upwind sémáját Roe által közelített Riemann módszerrel alkalmazták transzszonikus lapátrácsok analízisére. Furukawa és mások [21] olyan nagy pontosságú upwind véges térfogat módszert dolgoztak ki, amely egy harmadrendűen pontos MUSCL (Monotone Upstream Schemes for Conservation Laws) közelítést alkalmaz a Roe által approximált Riemann megoldóban. Yuan és mások (1993) által kidolgozott nagy pontosságú upwind véges differenciák segítségével diszkretizált módszerében egy negyedrendű MUSCL közelítést használ a közelítő Riemann megoldójában. Ezek a numerikus eljárások egy adott peremfeltétel esetén képesek meghatározni az áramlási paramétereket.
47
Az általános numerikus áramlástani szoftverek esetén nehézségekben ütközhet a forgó gépek karakterisztikájának meghatározása. A fő problémát az a kompresszorban lezajlódó jelenség okozza, miszerint a futólapátozás fordulatszáma, illetve a rotor mögött kialakult nyomásnövekedés határozza még a gépegységen keresztüláramló anyagmennyiség nagyságát, illetve a belépő relatív sebesség irányát. Az összenyomható közegek numerikus áramlástani modellezése során kidolgoztam a transzszonikus axiális kompresszor forgó lapátsorának kétdimenziós matematikai modellezésébe integrálható olyan numerikus bemeneti peremfeltétel számítási eljárást, amely figyelembe veszi a forgórész fordulatszámát és a futólapátozás mögött kialakult nyomásnövekedést, ezáltal teremtve meg a kapcsolatot az abszolút és a relatív jellemzők között. Ennek értelmében az általában előírt ( T to , p to , α ), és a másodfajú (numerikus) p bemeneti peremfeltételek közül kettő nem független, vagyis a kialakult nyomásnövekedés és a fordulatszám határozza meg a tömegáramot, illetve a belépő áramlás irányát. A módszer alkalmazásával a kompresszor karakterisztika numerikus módszerekkel történő meghatározása válik könnyebben megvalósíthatóvá. A peremen érvényes statikus nyomást, mint másodfajú peremfeltételt a számítási térből lineáris extrapoláció segítségével határoztam meg. Fizikai peremfeltételként a torlóponti nyomást p to és hőmérsékletet T to írtam elő. A bemeneti peremfeltételek meghatározásának következő lépéseként izentropikus egyenletek segítségével számítottam ki a statikus hőmérsékletet, a Mach-számot és az előperdítés nélküli axiális (abszolút) sebesség nagyságát: p T = T p to
to
γ −1 γ
T to − 12 , T és u = M γRT . M = γ −1
(2.154)
A fordulatszám és a 2D-s lapátprofil radiális helyzete ismert, ezért v kerületi sebesség ismeretében meghatározható az utolsó hiányzó fizikai peremfeltétel, vagyis a bementi relatív sebesség iránya az α . A szilárd fal és a kilépő peremek meghatározására a 2.4.2. alfejezetben leírt extrapolációs peremfeltételt alkalmaztam. A lapátrács számításakor egy egyszerű, ám idáig nem tárgyalt úgynevezett periodikus peremfeltételről is érdemes szót említeni. Mivel a modellezés során két lapát által határolt csatornában végeztem el a számításokat, így a belépő él előtt, tőle kifelé és a kilépő él mögött, tőle kifelé húzódó, egymástól lapátosztásnyira elhelyezkedő kvázi-áramvonalak között biztosítani kell a kapcsolatot. A számítás során, ha a futó koordináták elérik ezt a legszélső helyzetet, akkor a hiányzó peremértékek helyett a lapátosztásnyival eltolt, a számítási térhez tartozó értéket kell megadni. Ennek megfelelően
π kto
7
100%: v=311 [m/s] (kerületi sebesség)
m [ kg / s ] 16. ábra. 2D-s kompresszor rotor karakterisztika lesz biztosítva a végtelen lapátozás feltétele. 48
A disszertáció 2.2.1.1. alfejezetében ismertetett Jamenson-féle mesterséges viszkozitással kiegészített centrális diszkretizációs módszer napjaink, számítástechnikai szempontból, legjobb hatásfokkal használható módszere a transzszonikus áramlások szimulációjában. Ezért, az összenyomható áramlás ugyanazon sémában történő modellezése érdekében (3.1. alfajezet), valamint az előzőekben ismertetett numerikus peremfeltétel újfajta kiszámítása és a 2D-s kompresszor rotor karakterisztika felállítása érdekében ezt a módszert dolgoztam ki. Az eljárás validálását a 2.5. alfejezetben a mesterséges viszkozitás módszerével együtt végeztem el. A 2D-s karakterisztika meghatározására a NACA RM E52C27 [43] transzszonikus egyfokozatú rotort használtam fel. A 2D-s lapátrács előállítását a lapátozás legkülső áramfelületének síkba vetítésével állítottam elő. A 16. ábrán látható jelleggörbe sereg jól demonstrálja a kompresszorokban lezajló áramlástani folyamatokat, a mérésekkel való összehasonlítás a megfelelő adatok hiánya miatt még nem történhetett meg. A határréteg okozta másodlagos áramlások nem modellezhetők Euler alapú eljárásokkal, ezért a számítási idő csökkentése érdekében, a 3D-s megoldó helyett, a következő egyszerűsített 2D-s eljárást dolgoztam a 3D-s kompresszor rotorok karakterisztikájának meghatározására. Először felosztottam a lapátozást tengelyszimmetrikus, a tengellyel párhuzamos áramfelületek segítségével gyűrű alakú szegmensekre. Legyen a belépő A1 (17. ábra) keresztmetszet területarányaival súlyozott közepes sugáron elhelyezkedő áramfelület síkba vetítésével előállított 2D-s lapátrács az analízis alapja. A megoldás végeredményét tekintsük érvényesnek az A1 , állandó keresztmetszettel jelölt áramcsőben. A 17. ábrán látható meridionális metszet áramfelületekkel határolt szegmensei azonban nem állandók, így a megfelelő méretű kilépő keresztmetszet figyelembe vétele miatt, az ehhez a keresztmetszethez tartozó áramlástani paraméterek a következő analitikus iterációval határozhatók meg: u' 2 =
A2 ρ ' 2 ; A1u 2 ρ 2
V ' 22 = u' 22 +v 22 ; V' 22 −V22 T ' 2 = T2 − ; 2C p T' p' 2 = 2 T2
γ
γ −1 p 2 ; p' 2 ρ' 2 = ; RT' 2
ahol a vessző nélküli paraméterek jelentik az A1 keresztmetszethez tartozó 2D-s analízis eredményeit, a vesszős értékek pedig az A2 keresztmetszethez tartozó paramétereket. Az iteráció az első, a kontinuitás egyenlet kielégítéséig tart. A tengelyszimmetrikus felbontás sűrűsége jelentősen befolyásolja az eljárás pontosságát. A 2D-s közelítés számítógépi kapacitásigénye alulmarad a hasonló 3D-s alkalmazásokhoz képest, azért elsősorban az előtervezési feladatokban használható jó hatásfokkal.
A2 49
A1
2.7. Tervezési eljárásrendszer kifejlesztése a többfokozatú centrifugálkompresszorok összekötő-csatorna lapátozásának tervezésére A többfokozatú centrifugálkompresszorok forgórészei a diffúzoron, a fordító könyökön, illetve a visszatérő csatornán keresztül vannak összeköttetésben egymással (18. ábra). A lapátos vagy lapát nélküli diffúzor fő feladata, hogy a lapátos keréktől érkező közeg mozgási energiájának lehető legnagyobb részét statikus nyomásnövekedéssé alakítsa. A lapátozást összekötő-csatorna
18. ábra. A-C többfokozatú kompresszor magába foglaló visszatérő csatorna fő célja, hogy meggátolja az örvényképződést, mialatt a közeget a következő fokozat forgórészének bemenetéhez juttatja. A két komponenst köti össze a 180 fokos fordító könyök (18. ábra).
2.7.1. Bevezetés Az előzőekben említett, forgó mozgást nem végző alkatrészek kompaktsága igen fontos, mivel az áramlási sebességek csökkenése miatt nem csak a súrlódásból származó veszteségek, melyek az áramvonalak hosszával, illetve a sebesség négyzetével arányosak, hanem a gyártási költségek is csökkenthetők. Ez utóbbi igen fontos azoknál a nagynyomású kompresszoroknál (18. ábra), amelyeknél a nagyszilárdságú köpeny költségnövekedése az átmérő harmadik hatványával arányos. Nincs túl sok lehetőség a diffúzor és a fordító könyök (19. ábra) nyomásnövekedési, illetve veszteségi tényezőjének javítására ((2.161) egyenletek (59. o.)). A lapát nélküli diffúzorban lezajlódó nyomásnövekedés a be- és a kimeneti átmérőviszony függvénye és a nyomásnövekedési tényező javítása a befoglaló méretek növekedését vonja maga után. A lapátos diffúzor használatával kompaktabb geometriát érhetünk el, de ennek kedvezőtlen hatása van a működési tartományra. Jelentős veszteség jelentkezhet a fordító könyökben a leválás miatt, ami főként a hirtelen meridián görbület változása miatt alakul ki. Nehéz nyomásnövekedést elérni ebben a szegmensben. A visszatérő csatornában szintén problematikus a nyomásviszony javítása, egyrészt a fordító könyökben kialakult viszonyoknak köszönhető inhomogén belépő áramlás miatt, másrészt a nagy kerületi elcsavarás miatt, amit a lapátok hoznak létre az áramlás radiálistangenciális irányból axiálisba történő visszafordításkor, hogy a következő fokozat számára biztosított legyen a perdületmentes belépés (19. ábra). Thygesen [71] megmutatta, hogyha a lapátokat kiterjesztjük a fordító könyök bemenete és a lapát nélküli diffúzor kimenete közé, nagyobb hatásfok érhető el. Hasonló eredmény jelent meg az [47]-ben. Lapátkiterjesztéssel történő hatásfok növelés esetén ugyanis, a kisebb 52
nyomásveszteség miatt a belépő össznyomás nagyobb hányada fordítódhat nyomásnövekedésre, ami az elvárt nyomásviszony esetén − ugyanakkora hasznos teljesítmény elérésével − csökkenthető kompresszor befoglaló méreteket jelent. A jelentős meridián görbületbe foglalt erősen örvényes áramlásban elhelyezett lapátozás megtervezése nem könnyű feladat. A hagyományos tervezési eljárások nem alkalmasak ilyen körülmények között működő, jó hatásfokú és irányított diffúzorosságú lapátok előállítására. A célom egy olyan tervezési eljárásrendszer kidolgozása a lapátos fordító könyök és fordító könyök lapát belépő él
lapátnélküli diffúzor
visszatérő csatorna
lapát kilépő él
forgórész bemenet
19. ábra. Kompresszor fokozat meridián keresztmetszete visszatérő csatorna kombináció kifejlesztésére, amelyben egyaránt szerepel egy analitikus módszer és egy háromdimenziós inverz tervezési eljárás alkalmazása a lapátgeometria optimalizálására. További új eredménynek tekintendő azon numerikus vizsgálatok eredménye, amelyek a lapátelhajlításon keresztül próbálnak tovább javítani a tervezési specifikációkon. Az új tervezési eljárásrendszer segítségével előállított konfigurációk áramlástani tulajdonságait egy hagyományos (referencia) visszatérő csatornával (20. ábra) [57,65] hasonlítottam össze. A geometriák analizálását a CFX-TASCflow kereskedelmi szoftver segítségével végeztem el a k-ε turbulencia modell és a fal-törvény felhasználásával. CFX-TurboGrid program segítségével készítettem el a háromdimenziós hálót, 128 ponttal az áramlás, 24 ponttal a lapát osztás és 26 ponttal a lapát-terjedtség irányába. Két periodikus perem között helyezkedik el 3D-s lapát egy áramcsőben, amit a csatorna be és kilépő, illetve külső és belső fal pereme határol (A1-A5. ábrák (64.-68. o.)). A tervezés kiindulásaként egy levegő munkaközegű, jellemző geometriájú kompresszorfokozatot modelleztem (19. ábra), majd a VKI-i (von Kármán Institute) CCOD (Centrifugal Compressor Off-Design Ananlysis [73]) egydimenziós számítógépes analízis segítségével határoztam meg a csatorna főbb áramlástani paramétereit (a be- és a kilépésnél), melyeket a további számításokhoz is felhasználtam (2. táblázat, 60. o.). Az áramlási csatorna 19. ábrán látható meridián keresztmetszete nem változik a tervezés során. A csatorna bemeneti keresztmetszetét a lapátos-kerék kilépő szélessége, még a kilépő keresztmetszetét a következő fokozat bemenetének szélessége határozza meg.
53
2.7.2. Referencia geometria A referencia geometria a lapát nélküli diffúzorból, a fordító könyökből és végül a visszatérő csatornából áll. A lapátokat Rothstein [57,65] módszere alapján terveztem meg. A
20. ábra. Visszatérő csatorna lapát tervezése Rothstein módszere alapján lapátok belépő éle a fordító könyök kilépő keresztmetszeténél helyezkedik el, érintő irányban a központi áramvonallal. A kilépő él esetén a lapátszög radiális-axiális irányú, hogy az áramlás iránya megfelelő legyen következő fokozat számára. A szívott és nyomott oldal helyzete, vagyis a lapátvastagság állandó vagy enyhén divergens a központi áramvonal körül (20. ábra). A cél, az áramlás axiális helyzetbe történő fordítása, mialatt az átlagsebesség közel állandó. A lapátok be- és kilépő éle hegyes, a hirtelen keresztmetszet-változás okozta problémák csökkentése érdekében. Csak a belépő él lekerekített. A lapátgeometria azon részei, amelyek leginkább érzékenyek a leválásra jól megfigyelhetők a Navier-Stokes megoldóval végzett analízis során (A1. ábra, 64. o.). Az első leválási zóna közel a fordítókönyök bemenetéhez, a csatorna külső felületén található. Ez elsősorban a meridián sebességcsökkenés miatt alakul ki, amit a fordító könyök jelentős görbülete okoz. Egy másik leválási zóna is megfigyelhető, mégpedig a visszatérő csatorna belső falánál a lapát szívott oldalán a belépő éltől kezdődően. A jelenség egyrészt a meridián sebesség hirtelen csökkenése miatt jöhet létre a fordító könyök kilépő keresztmetszetétől kezdődően, valamint a nem megfelelő állásszög miatt, amit a belső és a külső csatorna fal közötti radiális irányú sebességkülönbség okoz. Az A1. ábrán a 3D-s lapát szívott oldalán a be- és a kilépés, valamint a csatorna külső és belső fala között végigfutó felület ábrázolásakor a sebességvektorok nincsenek a hálózás síkjába vetítve. Ennek köszönhetően a fordító könyök középső részénél hirtelen sebességváltozás látható.
2.7.3. Állandó terhelésű lapát tervezése (ÁLT (Állandó Lapát Terhelés)) A referencia csatornában tapasztalt leválások elkerülésének érdekében, a bevezetésben említett új módszert használtam, ami szerint az örvénymentesítő-lapátot a lapátnélküli diffúzor kimenetéig célszerű kiterjeszteni (19. ábra). A lapát belépő éle elég távol helyezkedik 54
m ds
W nyo
dm
Wszo
δlv
s
R+∆R
R
2πR z
21. ábra. Az analitikai lapáttervezés ellenőrző felülete el a forgórész kilépő keresztmetszetétől, így ez a változás a működési tartományt nem befolyásolja. A kilépő él helyzete, illetve a lapátszám (18) változatlan marad. A komplett 3D-s lapát alakjának, vonalvezetésének, illetve vastagságának meghatározása nem egyértelmű, igen komplex probléma. A lapátterhelés egyaránt függ a lapát geometriájától, a sugár irányú, illetve a meridián irányú görbület változásaitól. A lapát görbületének és vastagságának meghatározására egy analitikai módszert dolgoztam ki, amelynek célja az állandó lapátterhelés (ÁLT) figyelembevétele a tervezés során. Az eljárás végeredménye a lapát vázvonalának közel-optimális kerületi elhelyezkedése és vastagsága a meridián metszet ismeretében mind a külső, mind a belső csatorna falon. Ezek után a 3D-s lapát lineáris interpolációval állítható elő. A tervezési folyamat a β lapátszög-eloszlás (sugártól mérve) meghatározásával kezdődik. A lapát erőt gyakorol ugyan az áramlásra, de jó közelítésnek tekinthető az ideális, örvénymentes áramlás feltételezése: ∇ ×V = 0 .
(2.155)
Legyen nulla a cirkuláció két szomszédos lapát, illetve két állandó sugarú kör által meghatározott ellenőrzőfelület felett (21. ábra). Ekkor, a következő összefüggés írható fel a lapátterhelés és a β áramlási irány között [74].
2π δ lv Wnyo − Wszo = cos β l − z R cos β l
d Wm R tgβ á , dm
(
)
(2.156)
m& dm és Wm = egységnyi lapátterjedtségre. Elhanyagolva a ds ρ 2πR lapátvastagságot ( δ lv = 0 ), továbbá feltéve, hogy β á = β l , a következő összefüggés írható fel az áramlás iránya és a terhelés között:
melyben
cos β l =
2π 1 dβ . Wnyo − Wszo = C t = cos β 2 z cos β dm
(2.157)
Az állandó lapátterhelés (Wnyo-Wszo=Ct) előírásával: KÉ
β = Ct
∫ cos
β dm
(2.158)
BÉ
összefüggésre jutottam. A lapát középvonal irányszög (sugártól mérve) eloszlása a (2.158) egyenlet, belépő és kilépő él közötti integrálása után adódik. A Ct szabad paraméter biztosítja, 55
hogy az áramlás iránya axiális legyen a lapát kilépő élénél. Jó közelítésnek tekintettem azt a feltételezést, miszerint az áramlás érintő irányú azzal a lapát középvonallal, amely a lapátterjedtség közepén helyezkedik el. A belépő, illetve a kilépő él környezetében azonban az állásszög hatás, illetve a lapát-terjedtség irányú, elsősorban a meridionális görbületi hatásokra visszavezethető sebességváltozások miatt az approximáció pontossága kivételt képez. A csatorna külső és belső felületén meghatározandó átlagos lapátszög eloszlás számítása korrekcióra szorul, mint ahogy az a 22. ábrán látható. A módosításra a csatorna külső és belső fal hosszának különbsége, illetve az eltérő lapátterhelés miatt volt szükség. Elliptikus eloszlású korrekciós lapátszög eloszlást használtam a fordítókönyökben, amely kiegyenlítettebbé teszi a lapátterjedtség irányú, meridián sebességváltozásból adódó lapátterhelés különbséget a csatorna külső és belső felülete között. A lapát középvonal β Külső fal mentén Belső fal mentén
m
22. ábra. Lapátszög eloszlás a csatorna külső és belső fala mentén
kerületi irányú θ szög eloszlása a meridián metszet (sugár és tengelyirány), valamint a lapátszög függvényében a következő összefüggéssel számolható:
tan( β l ) dm . R BÉ
KÉ
θ ( m, βl ) =
∫
(2.159)
A θkülső=θbelső feltételnek (nincs lapátelhajlítás) kell teljesülnie a lapát kilépő élénél. Az elliptikus korrekcióban az ellipszis kistengelyének nagyságával lehet beállítani, hogy a kilépő él (θkülső-θbelső)KÉ különbsége egy elfogadható határ alatt maradjon. Az analitikus lapáttervezés utolsó lépése a lapátvastagság meghatározása. A szívott és a nyomott oldal pontjainak távolságát a lapáthossz első 75%-ában szintén egy ellipszis segítségével határoztam meg. A rövidebbik tengely 10%-a a hosszabbiknak, így 7,5% maximális relatív lapátvastagság érhető el mind a külső, mind a belső csatorna falakon. A lapát másik felében lineáris eloszlás alkalmazása tűnt a legcélszerűbbnek a maximális vastagság helyétől az előírt kilépő él vastagságig. Az új analitikus módszert számítógépes program keretében valósítottam meg. A program neve: BLADECONTOUR (31. ábra (62. o.)). A lapát előállítását és a csatornába való elhelyezését követően elvégeztem a 3D-s Navier-Stokes analízist, amelynek egy eredménye az A2. ábrán (65. o.) látható. A csatorna meridián kontúrja torzult, köszönhetően a térben görbült lapátfelülethez közeli háló, illetve a hozzátartozó sebességvektorok lehető legjobb megjelenítésnek. A fordító könyök bemeneténél, a referencia geometriánál tapasztalt leválás 56
megszűnt, de a visszatérő ágban a lapát szívott oldalán a fordító könyök kilépő keresztmetszetétől kezdődően jelen van, ami jól látszik a sebességvektorok irányából. A probléma oka főként a gyors, lokális lapátterhelés változásokra vezethető vissza, amit a csatorna meridián görbületváltozásai okoznak. A jelenség részletes magyarázatára következő alfejezetben kerül sor.
2.7.4. Inverz lapáttervezés A komplex 3D-s áramlás jobb vezetése érdekében inverz lapáttervezési módszer alkalmazásával célszerű tovább tökéletesíteni a lapátot. Ez egy olyan iteratív eljárás, melynek során a kiindulási geometria addig módosul, amíg az elvárt, tervezett sebesség vagy nyomásprofil ki nem alakul a lapát kontúrja mentén. A program 3D-s Euler megoldót használ módosított peremfeltétekkel. Az iterációk során a profilra merőleges sebeség-komponens nagyságának meghatározása az aktuális és az előírt paraméterek segítségével történik. Miután ismert a normális irányú sebesség, az áteresztő modell segítségével a program meghatározza azt az új áramvonalat, amely egyrészt kielégíti az előírt áramlástani paraméter eloszlást, másrészt az új lapátgeometriát is szolgáltatja. A használt eljárást Demeulenaere [18,19] fejlesztette ki. A tervezési folyamat során a felhasználó döntése, hogy a lapátterjedtség mentén hány szegmensben kívánja a Mach-számot előírni mind a szívott, mind a nyomott oldalon. A jelen munkában a csatorna külső, illetve belső felületén elhelyezkedő profilra írtam elő Mach-szám eloszlást, míg a közbenső háló pontokban interpolálással adódott a megkívánt disztribúció. Az inverz tervezési eljárás során jelentkező fő probléma nem csak az olyan optimális Mach-szám eloszlás megtalálása, amely a minimális veszteséget adja a maximális nyomásnövekedés mellett, hanem a tervezés során kialakuló lapátnak ki kell elégítenie a minimális vastagságra, a maximális lapátelhajlásra, illetve a megfelelő be-, illetve kilépési lapátszögre vonatkozó megkötéseket. Jelen esetben mindezek betartása igen nehéz faladat, mivel a lapát rövid, de a vázvonala igen hosszú. Annak ellenére, hogy a relatív lapátmagasság kicsi, olyan áramlástani állapotok alakulhatnak ki a csatornában, amelyek jelentős lapátelhajlást (különbség a θbelső és a θkülső között) eredményezhetnek. A numerikus vizsgálatok eredményei azt mutatták, hogy úgy lehet a legkönnyebben megtalálni egy reális és optimális Mach-szám eloszlást, ha a kiindulási lapát adatait használjuk fel néhány fontosabb áramlástani megfontolás figyelembevételével. A program használatakor nem előírás, de a konvergencia gyorsítása, illetve az említett megkötések betartatása érdekében érdemes olyan
23. ábra. Kiindulási lapát Mach-szám eloszlása (+) a belső (bal oldal), illetve külső csatorna faltól 12,5%-ra eső profilkontúron
57
kiindulási geometriát használni, ami nincs túl messze az elvárttól, olyat például, mint aminek a tervezését az előző alfejezetben ismertettem. A disszertációban használt inverz tervezési eljárás az ÁLT módszerével tervezett 3D-s lapát, valamint az előírt Mach-szám eloszlás implementálásával kezdődik. A megkívánt eloszlás az inverz kód Euler-megoldó része által szolgáltatott eredményekből származtatható (23. ábra, +++ jel). Két kritikus terület figyelhető meg a kiindulási Mach-szám eloszlásban, amelyeknél a tervezés során jelentős javulás érhető el. Az első a fordítókönyök belépésénél található. A csatorna külső falán csökken a meridián sebesség nagysága a konkáv görbületi hatás miatt. A perdület-állandóság miatt az érintő irányú sebesség változatlan, ezért az áramlás tangenciális irányba fordul (növekvő βáramlás), így lokálisan a lapátterhelés is nő (B). Hasonló
24. ábra. Inverz újratervezett lapát előírt (o) és számított (+) Mach-szám eloszlása a belső (bal oldal), illetve a külső csatorna faltól 12,5%-ra eső profilkontúron
keresztmetszetben, de a csatorna belső falán az előző jelenség ellenkezője zajlik le (A). A konvex görbületi hatás miatt a meridián sebesség a fal környezetében megnő és az áramlás axiális-radiális irányba való fordulását eredményezi, minek következtében csökken a lapátterhelés. A második kritikus zóna a fordító könyök kilépésénél található, ahol a görbületi hatások befolyása megszűnik, és az áramlás uniformissá válik a lapátterjedtség mentén. A csatorna belső falán a meridián sebességkomponens hirtelen lecsökken (konvex görbületből belső profil
θ
KÉ
külső profil
KÉ
R
θ
58
BÉ
BÉ
25. ábra. Eredeti (szürke árnyalatú) és az újratervezett (+) profil a belső (bal oldal), illetve a külső csatorna faltól 12,5%-ra eső lapátszegmensben
történő kiáramlás), az áramlás tangenciális irányba fordul, megnövelve lokálisan a lapátterhelést (C). A csatorna külső falán a hirtelen csökkenő lapátterhelés a meridián sebesség komponens növekedésének köszönhető, amint a gáz kiáramlik a konkáv görbületből a sík falra. Az inverz lapát-áttervezési eljárás fő célja tehát, hogy megszüntesse a meridián keresztmetszet okozta lokális sebességváltozásokat, ezáltal szüntetve meg a visszatérő csatornában kialakuló leválásokat. A 23. ábrában körökkel jelölve látható az inverz tervezés számára előírt egyenletes Mach-szám eloszlás. A meridián kontúr nem változik a tervezés során, így nem várható nagy átlagsebesség változás. A lapáthossz mentén az erősen változó helyett, egy közel állandó lapátterhelést írtam elő, különösen ügyelve az KÉ
∫ (W
szo
− Wnyo ) ds
(2.160)
BÉ
integrál állandóságának betartására, elsősorban a kilépő áramlás szögének állandó értéken tartása érdekében. Az inverz megoldó 550 időlépés után érte el a konvergencia kritérium értékét, és az előírt, valamint az aktuális Mach-szám eloszlás teljes egyezőséget mutat, mint ahogy az a 24. ábrán is R.θ látszik. Az eredeti (szürke) és az újra tervezett (+++-el 12,5%-ra a shroud jelölve) lapátprofil összehasonlítása a 25. ábrán belső faltól látható. A legszembetűnőbb változás a fordító könyöknél, illetve, a visszatérő ágban jelentkezik, TE ami a 26. ábrán, az m-Rθ síkban jobban megfigyelhető. A fordítókönyök bemeneténél, a csatorna külső falán megnövekedett a profil KÉ vastagsága. Az így létrejött lokális keresztmetszetcsökkenés kompenzálja a konkáv görbületi hatás okozta meridián sebességcsökkenést és az áramlást LE a csatorna belső fala felé kényszeríti. Másrészt ezzel ellentétesen a konvex görbületi hatás BÉ kompenzálása miatt a csatorna ugyanazon, de a hub a 12,5%-ra belső keresztmetszetében a lapát profil külső faltól elvékonyodott. A csatorna külső részében a fordítókönyök kilépő keresztmetszeténél csökkent a lapátvastagság, ellenben a belső falnál megnövekedett. A jelenség magyarázata ellentétes értelemben teljesen hasonló az előzőekben m tárgyaltakéhoz, kiérve ugyanis a meridián görbület 26. ábra. Eredeti (szürke árnyalatú) változás okozta sebességváltozások hatása alól, és az újratervezett (-) ellenkező lokális változások fejtik ki hatásukat lapátprofil csatorna külső és belső felületén.
2.7.5. Tervezési specifikációk Az olyan áramlástani gépek esetén, amelyekben a nyomásnövekedés elérése a cél, a két legfontosabb paraméter a nyomásnövekedési tényező, illetve a veszteségi tényező:
Cp =
st p out − p inst
p into − p inst
és ω =
to pinto − pout pinto − pinst
.
(2.161)
59
Az összehasonlító számításokat a 2.7.1. alfejezetben már ismertetett CFX-Tascflow Navier-Stokes megoldóval végeztem el. A különféle tervezési esetek főbb áramlástani paraméterei a 2. táblázatban találhatók, melyekben a megoldó által számított paraméterek átlagolt értékek, a lokális tömegáram figyelembevételével kerültek súlyozásra. A számítások
p into [Pa] p inst[Pa] p outto [Pa] p outst [Pa]
NEM KITERJESZTETT 299699,1 159038,5 236665,3 225215,1
299526,6 182298,2 261108,8 258238,5
ÁLT + INVERZ TERVEZÉS 299696,6 174936,4 264954,2 258414,2
ω
0,44813
0,3277
0,27847
Cp m& [kg/s]
0,4705 4,68
0,6478 4,5
0,669106 4,64
TERVEZÉS
ÁLT
2. táblázat. A tervezési folyamat eredményeinek összehasonlítása
során tömegáram állandó értéken tartása volt a cél, de mivel a kimeneten fizikai peremfeltételként a statikus nyomást adtam meg, ezért a tömegáram az össznyomás-veszteség függvényeként alakul. A csatornában kialakult áramlás Mach-száma viszonylag kicsi ( M ≈ 0,28 ), ezért a különböző tervezések kilépő statikus nyomás eltérései nem gyakorolnak jelentős befolyást a Cp és az ω dimenziómentes paraméterekre. A csatornába lépő közeg áramlás iránya, mint peremfeltétel, minden esetben ugyanakkora: 67,6o, aminek következtében tömegáram változás nem lesz hatással az állásszögre. Az eredmények jól tükrözik azon elvárásokat, amelyek a korábbi fejezetekben kerültek tárgyalásra. A statikus nyomásnövekedés jelentősen növekedett állandó lapátterhelés (ÁLT) módszere és a lapátkiterjesztés bevezetésével. A veszteségi tényező változása szintén kedvezően alakult (2. táblázat). Kismértékű, további teljesítménynövekedés érhető el az inverz tervezési eljárás segítségével. Ez egy kicsit meglepő, mivel a Navier-Stokes számítás azt mutatta, hogy a korábbi esetekkel ellentétben, profil szívott oldalán megfigyelhető intenzív leválás megszűnt (A3-A5. ábra (66.-68. o.)).
2.7.6. Lapátelhajás hatása A hosszú és keskeny csatornaméret következtében jelentős szekunder (másodlagos) áramlások [40] alakulnak ki, létrehozva a csatorna- és a lapátörvény szuperpozícióját (27-28,
Rg
nyo
n
külső csatorna fal
v
szo
60
Rg
2
∂p v ~ ∂n R g
belső csatorna fal
n
27. ábra. A Csatorna- (bal oldal) és a lapátörvény kialakulása [77]
külső csatorna fal
nyo
külső csatorna fal
szo
nyo
belső csatorna fal
szo
belső csatorna fal
28. ábra. A Csatorna- (bal oldal) és a lapátörvény kialakulása
30. ábrák). Az előzőekben használt Euler alapú inverz tervezési módszer nem veszi figyelembe a határréteg okozta örvényképződéseket, így a lapátelhajlítás hatásának figyelembevétele csak a Navier-Stokes szoftver segítségével lehetséges. Napjaink intenzíven fejlődő kutatási területei közé tartozik a szekunder áramlások által felemésztett energiának lapátelhajlítással vagy nyilazással való csökkentése [75]. A geometria hosszának, illetve a rövid lapát magasságnak köszönhetően nem valószínű, hogy a lapátnyilazás jelentős befolyást gyakorolna az áramlásra, ezért csak a lapátelhajlítás hatását vizsgáltam. A megvalósítást tekintve, a csatorna külső felületén elhelyezkedő profil ∆θ = 4,0o-os, a szívott lapátoldal felé történő elforgatásával a negatív, míg a nyomott oldal felé történő elforgatással pozitív lapátelhajlítást valósítottam meg. A csatorna belső felületének helyzete változatlan, ezért az új 3D-s lapát az egymáshoz képest eltérő radiális helyzetben elhelyezkedő profilok megfelelő pontjainak összekötésével adódik, ami 27,0o-os maximális lapátelhajlást eredményez a fordító könyök legnagyobb sugárhoz tartozó keresztmetszetében. külső csatorna fal
nyo
szo
belső csatorna fal
29. ábra. Módosult nyomás eloszlás negatív lapátelhajlítás esetén
Megkülönböztetett figyelmet érdemel a negatív lapátelhajlítás. Elvárásaim szerint ekkor a nyomáseloszlás a 29. ábra szerint módosul, így csökkentve a 30. ábra felső részén látható örvény-szuperpozíció intenzitását. A lapátelhajlítás tervezési paraméterekre gyakorolt hatása a 3. táblázatban található. A nyomásnövekedési tényező negatív lapátelhajlítás estén tovább növekedett, ellenben a pozitív esettel, ahol kb. 5% csökkenés volt tapasztalható.
ω
Cp
ÁLT 0,3277 0,6478
ÁLT+INVERZ 0,2785 0,6691
NEGATÍV ELH. 0,2635 0,6817
POZITÍV ELH. 0,2802 0,6380
3. táblázat. Lapátelhajlítás hatása
A 30. ábrán látható sebességvektorok a lapát kilépőél közeli csatorna keresztmetszet (A4. ábra, nyilakkal jelölve (67. o.)) hálófelületére vannak vetítve és jól látszik a szekunder áramlás hatása. A metszet közepén látható a lapát, melytől jobbra is és balra is fél-fél lapátosztás távolságban áramlik a lap síkjából kifelé a közeg. A két szélső kontúron 61
szimmetria peremfeltételek kerültek megadásra, így a teljes áramcsövet megkaphatjuk, ha a jobb oldalt eltoljuk, úgy, hogy a jobb oldali élét a bal oldali szegmens bal oldali élével hozzuk fedésbe. A kiindulási esetben (nincs lapátelhajlítás) az óra mutató járásával ellentétes irányú örvény (30. ábra, fent) alakult ki, melynek intenzitása nemcsak csökkent negatív lapátelhajlítás esetén, hanem lokálisan, a külső csatornafalnál a forgásirány is megváltozott. Pozitív lapátelhajlítás következtében a megnövekedett veszteség az eredeti örvény további erősödésére, illetve a csatorna belső részében egy új, az óramutató járásával ellentétes forgásirányú örvény kialakulására vezethető vissza.
2.7.7. Értékelés A hagyományos többfokozatú cetrifugálkompresszorok összekötő csatornájában a
B L A D E C O N T O U R
3 D
INVERZ MÓDSZER
L A P Á T G E O M E T R I A
L A P Á T E L H A J L Í T Á S
REFERENCIA GEOMETRIA +
-
NS MEGOLDÓ TERVEZÉSI SPECIFIKÁCIÓK 30. ábra. Lapátelhajlítás hatása a csatorna örvényre (A4. ábra alapján)
62
31. ábra. Az új 3D-s lapáttervezési eljárásrendszer
fordítókönyök bemenetének külső felületén, illetve a nem kiterjesztett lapát szívott oldalán, a belépő éltől kezdődően van a legnagyobb valószínűsége intenzív leválási buborék kialakulásának. A jelenséget mindkét esetben elsősorban a meridián görbületi hatások miatt létrejövő hirtelen sebességcsökkenés, illetve a második esetben az ennek következtében előálló állásszög-változás váltja ki. A 31. ábrán látható új eljárásrendszer kidolgozásával egy olyan tervezőeszközt fejlesztettem ki, amely segítségével javítható a veszteségi, illetve a nyomásnövekedési tényező a többfokozatú centrifugálkompresszorok összekötő csatornájában. Az örvénymentesítő-lapát diffúzor kimenetéig történő kiterjesztésével jelentős javulás érhető el a tervezési paraméterek tekintetében. A kiterjesztés mellett szintén fontos a megfelelő diffúzorosság biztosítása a lapáthossz mentén, amelyek következtében lehetővé válik a kompresszor befoglaló méreteinek csökkentése. Az ÁLT és az inverz lapáttervezési eljárások használatával megszűnt a szívott oldali leválás, tovább lehetett növelni a teljesítményt. Megfelelően megválasztott irányú és nagyságú lapátelhajlítás alkalmazásával, csökkentve a szekunder áramlások okozta örvények intenzitását, szintén javulás tapasztalható a nyomásnövekedési, illetve veszteségi tényezők tekintetében. A további optimalizáció érdekében a fordító könyök görbületi sugarának, illetve a visszatérő csatorna diffúzorosságának hatását érdemes megvizsgálni.
63
A Navier−Stokes-megoldó eredményei
A1. ábra. Sebességvektorok a lapát szívott oldalán elhelyezkedő hálófelületen a referencia geometria esetén
64
A2. ábra. Sebességvektorok a lapát szívott oldalán elhelyezkedő hálófelületen ÁLT lapáttervezés esetén
65
A3. ábra. Sebességvektorok a lapát szívott oldalán elhelyezkedő hálófelületen az inverz módszerrel újratervezett ÁLT lapáttervezés esetén
66
67
A4. ábra. Sebességvektorok a lapátterjedtség közepén elhelyezkedő Fig. A5. – Velocity vectors in the mid-span plane for inverse redesign blade hálófelületen az inverz módszerrel újratervezett ÁLT lapáttervezés with CBL esetén
68
A5. ábra. Állandó statikus nyomásgörbék a lapátterjedtség közepén Fig. A6. – Static pressure iso-contours in the mid-span plane for inverse elhelyezkedő redesignhálófelületen blade with CBLaz inverz módszerrel újratervezett ÁLT lapáttervezés esetén
III. Fejezet: Numerikus módszerek és alkalmazásuk az összenyomhatatlan ideális áramlás modellezésére Az összenyomhatatlan közegek áramlását leíró, időfüggő konzervatív Euler- vagy Navier−Stokes-egyenletek numerikus megoldása különös figyelmet igényel, mivel a sűrűség idő szerinti differenciálhányadosa nem szerepel a tömegmegmaradás egyenletében, vagyis valamilyen módon kapcsolatot kell teremteni a sebesség-, és a nyomásmező között. Napjainkban, az összenyomhatatlan áramlások modellezésének két legelterjedtebb módszere a primitív-változók, illetve az örvényesség-áramfüggvény segítségével felírt kiindulási egyenleteken alapszik. Az örvényesség-áramfüggvény módszer alkalmazása esetén az alapegyenletekben szereplő sebességkomponensek az örvényesség, illetve az áramfüggvény bevezetésével helyettesíthetők, amelyekre egy megfelelően választott diszkretizációs eljárás után megoldhatók az egyenletek. Végül a sebességmező könnyen visszaszámolható az áramfüggvény definíciójából. A nyomás nem szerepel ugyan explicit alakban az egyenletekben, de az impulzus-megmaradás egyenletéből levezethető un. Poisson-egyenlet segítségével, amelyben a nyomás már explicit formában szerepel vagy az örvényességáramfüggvény vagy a sebességkomponensek függvényeként, könnyen meghatározható. Az örvényesség-áramfüggvény módszer 3D-s kiterjesztése nem egyértelmű, mivel az áramfüggvény nem létezik a háromdimenziós térbeli áramlásra. Az elkövetkezendőkben három olyan módszerről esik szó, amelyek túlmutatnak ezen a problémán. Az első az örvényesség-potenciál módszer, amelyben a vektor potenciál segítségével általánosíthatók az egyenleteket. A módszer nagyobb számítógépi kapacitást igényel, mint a primitív változós eset. A második módszer a duális-potenciál módszere, amelyben a sebesség egy skalár- és egy vektor-potenciál jellegű mennyiségre van szétválasztva. A megoldás során e potenciálok meghatározása a cél. Az eljárás során a peremfeltételek specifikálása összetett problémát jelent (Morino 1986; Gegg és mások 1989). A harmadik módszer az örvényesség-sebesség módszer, amely a két áramlástani paraméter egyidejű meghatározását jelenti. Az eljárásról bővebb információt Agarwal 1981, Gastski és mások 1982, Gui, Stella 1988 munkáiban találhatunk [33]. Összefoglalásként megállapítható, hogy az örvényesség-áramfüggvény módszer nem tekinthető 3D-s áramlások modellezésének általános eszközének, mivel a numerikus sémák, illetve a peremfeltételek meghatározása igen összetett [33]. A primitív-változós módszer alkalmazása túlmutat az előzőekben tapasztalt hátrányokon, ugyanakkor a térbeli kiterjesztés sem jelen olyan nehézségeket az egyenletek felírása során, mint az örvényesség-áramfüggvény módszer esetén. Az eljárás lényege, hogy az összenyomhatatlan Euler- vagy Navier−Stokes-egyenletek a primitív változóikra (u, v, w, p) oldjuk meg. Az eljárás alkalmazása igen elterjedt mind 2, mind 3D esetén. A módszeren belül két fő kategóriát lehet megkülönböztetni; a nyomás korrekciós módszereket és a kapcsolt módszereket. A nyomás korrekciós módszerek egyaránt alkalmazhatók stacionárius és időfüggő áramlástani egyenletek megoldására. A módszer olyan iteratív eljárásokon alapszik, amelyek kapcsolatot teremtenek a sebesség-, és a nyomásmező között. Egy kezdeti nyomásmező felvétele után az impulzus-megmaradás egyenletének numerikus megoldása során áll elő a sebességmező, amely nem fogja kielégíteni a divergenciamentes tömegmegmaradás egyenletét, ezért korrekcióra szorul. A korrigált sebességprofil szintén hatással van a kiindulási nyomásra, amelyet leggyakrabban az impulzusmegmaradás képletéből levezethető Poisson-egyenlettel lehet újraszámolni, korrigálni (a Poisson-egyenlet teremti meg lényegében a kapcsolatot a sebesség-, és a nyomásmező között). Az iteráció(k) 69
végeredményeként előálló sebesség-, és nyomáseloszlás egyaránt kielégíti a tömeg, és impulzusmegmaradási egyenleteket. A három legelterjedtebb nyomáskorrekciós módszer a MAC (Marker and Cell) módszer (Harlow és Welch 1965), a SIMPLE (Semi Implicit Method for Pressure Linked Equations) típusú módszerek (pl. SIMPLEC, SIMPLER) (Caretto és mások 1972, Patankar 1980) és a PISO (Primitive-Variable Implicit Separator) módszer (Issa 1986) [28]. A kapcsolt módszerek nagy előnye, hogy nincs szükség a sebesség és a nyomás különszámolására, mivel a nyomás is szerepel minden egyenletben. A stacionárius áramlás esetén, az egyenletek felépítése hasonló az összenyomható áramlás egyenleteihez, mert a kontinuitási egyenlete egy időszerinti differenciál-hányadossal, az úgynevezett mesterséges kompresszibilitási taggal egészül ki. Az időben beállt folyamat elérését követően az új tag hatása érvényét veszti, így teljesül a konzisztencia. Az ilyen típusú eljárások szintén használhatók a tranziens folyamatok leírására is, pl. a kettős időlépés módszer alkalmazása esetén. A kapcsolt eljárások egyik legelterjedtebb változata a pszeudo-kompresszibilitási módszer, amelyet Chorin (1967) vezetett be a véges differenciák módszerébe [11]. Az egyenletek formailag és tartalmilag is hasonlítanak az összenyomható áramlás kiinduló egyenleteire, ezért az összenyomható és az összenyomhatatlan áramlások egy sémában történő modellezhetősége érdekében ezt az eljárást használtam fel, a véges térfogat módszerével kombinálva, az új, komplex numerikus séma kidolgozására. A következő alfejezetben kidolgozott numerikus módszer és a hozzátartozó különféle peremfeltétel számítási eljárások az ideálisnak és összenyomhatatlannak feltételezett áramlások időben állandó 2D-s modellezésére alkalmas.
3.1. Pszeudo-kompresszibilitási módszer 3.1.1. Alapegyenletek, véges térfogat módszer, peremfeltételek és validálás A 2D-s összenyomhatatlan Euler-egyenletek dimenziós, konzervatív formában a következőképpen írhatók fel [30]:
∂U ∂F (U ) ∂G (U ) + + =0, ∂t ∂x ∂y
(3.1)
ahol T
T
p p U = [0 ,u , v ] , F ( U ) = u ,u 2 + ,uv és G( U ) = v ,uv ,v 2 + , ρ ρ T
(3.2)
amely a pszeudó-kompresszibilitási módszer alkalmazása esetén a következőként módosul [11]:
[
]
T
[
U = P / β 2 ,u , v , F ( U ) = u ,u 2 + P ,uv
]
T
[
]
T
és G( U ) = v ,uv ,v 2 + P .
(3.3)
A (3.3) egyenletekben β a kompresszibilitási tényező és P = p ρ . A konzisztencia biztosítása érdekében, az időben állandósult állapot elérése esetén a kontinuitási egyenlet első tagja eltűnik, így a kiindulási (3.1) egyenlet áll elő. Numerikus szimulációk arra engedtek következtetni, hogy a β paraméter megválasztása jelentősen befolyásolja a konvergenciát és a pontosságot. Ennek figyelembevételével β = 3 értéket használtam a programok futtatása során. r r r Az összenyomható áramláshoz hasonlóan, a H = Fe x + Ge y eredő konvektív fluxusvektor
bevezetésével a (3.1) egyenlet összevont formában a következőképpen írható fel: 70
∂U v v + ∇ ⋅ H (U ) = 0 . ∂t
(3.4)
Integrálva az (3.4) egyenletrendszert a Γ oldalakkal határolt Ω ellenőrző felület felett, valamint alkalmazva a Gauss-Osztogradszkij tételt a következő egyenlet állítható elő: rr ∂ Ud Ω + H (3.5) ∫Γ ndΓ = 0 , ∂t ∫∫ Ω r amelyben n = (n x , n y ) a Γ interfészből kifelé mutató egységvektor. A (3.5) egyenlet
jelentése szerint a Ω tartománybeli U konzervatív változók időegység alatti megváltozása v egyenlő a Ω tartomány Γ felületén belépő, illetve veltávozó eredő H fluxussal. Az v összenyomható áramláshoz hasonlóan (2.2. alfejezet), a Hn eredő fluxus a normális irányú tagokkal a következőképpen írható fel:
ahol
Vn
Vn β 2 rr H n = Hn = uVn + Pn x , vV + Pn y n rr r r r r ( Vn = Vn = (ue x + ve y ) (n x e x + n y e y ) = un x + vn y
(3.6)
a
normális
irányú
sebességkomponens. A folytonosról diszkrét alakra történő átalakítás során a véges térfogat feletti U megoldásfüggvényt az U j diszkrét ismeretlennel közelítettem. Az Ω j ellenőrző felület N f számú kontúrjára elvégzendő összeggel helyettesítve a (3.5) egyenletbeli második integrált, a következő j pont differenciálegyenletre jutottam:
körüli
szemi-diszkretizált
közönséges
1 Nf ∂ Uj =− ∑ [H n ] j ,k Γ j ,k , Ω j k =1 ∂t
amelyben
[H n ] j ,k
elsőrendű
(3.7)
az Ω j ellenőrzőfelülethez tartozó, a j és k csomópontok között
elhelyezkedő Γ j ,k oldalhosszra leképezett konvektív fluxus. A
[H n ] j ,k
mennyiség a
legalacsonyabb szintű közelítéssel az ellenőrző felület feletti állandó paraméterek segítségével a következőképpen állítható elő: ~ ~ Hn U L + Hn U L ~ ~ L R H n = H n U ,U = . (3.8) 2
(
( )
)
( )
Az egyenletben szereplő U L = U j és U R = U k állapotváltozók a Γ j ,k határ két oldalán elhelyezkedő állapotváltozók. A numerikus stabilitás megőrzése érdekében Jameson-féle mesterséges viszkozitás negyedrendű tagját alkalmaztam [35]: ∂ 1 Uj =− ∂t Ωj
Nf ∑ [H n ] j ,k Γ j ,k − D , k =1
(3.9)
ahol már a 2.2.1.1. alfejezetből ismert D = Dx + D y , vagyis Dx = d
1 i+ , j 2
−d
1 i− , j 2
, illetve D y = d
i , j+
1 2
−d
i , j−
1 2
(3.10)
disszipációs tag. Ez utóbbi két egyenlet jobb oldalán hasonló kifejezések találhatók, például az (i , j ) és az (i + 1, j ) csomópontok között felírható:
71
d
1 i+ , j 2
= −ε (4 1) (U i + 2 , j − 3U i +1, j + 3U i , j − U i −1, j ) . i+ , j 2
(3.11)
Az ε (4 ) hangoló paraméter értéke: 1 128 . A kiindulási egyenletek hiperbolikus típusúak, ezért a fizikai és numerikus peremfeltételek számát a 2.4. alfejezetben leírtaknak megfelelően a karakterisztikák módszere segítségével st , és szimmetria fal perem). A másodfajú peremfeltételek határoztam meg ( pinto , α in , pout értékét a 2.4.2. alfejezetben leírt extrapolációs módszerrel állapítottam meg. A mérnöki alkalmazást tekintve, a 3.3. alfejezetben található sugárszivattyú tervezési irányelveinek kidolgozása során felvetődött egy kérdés: hogyan lehetne a belépő torlóponti helyett a bemeneten ismert statikus nyomást peremfeltételként megadni? A sugárszivattyúk esetén, a fúvókából kiáramló nagysebességű áramlás depressziót hoz létre a keverőtérben, pinto=105 Pa
pinst=105 Pa
32. ábra. Torlóponti (baloldal) és statikus nyomás bemeneti peremfeltételek hatása az állandó statikus nyomás görbékre
aminek következtében a beömlő csatorna bemenetén jelentkező nagyobb környezeti statikus nyomás miatt jön létre a szállítandó közeg keverőtérbe áramlása (41. ábra). Abban az esetben, ha össznyomást adunk meg ezen a belépő peremen, akkor a kialakult áramlási sebesség nagyságától függően a peremen előálló statikus nyomás nem fog megfelelni a valóságban itt megjelenő környezeti nyomásnak. A statikus nyomás bemeneti peremként való megadása érdekében kidolgoztam egy új eljárást. A kérdéses bemenethez közeli keresztmetszetben numerikus integrálással meghatároztam a tömegáramot, majd a belépő keresztmetszet ismeretében egy fiktív belépő sebességet. Homogén, a belépő keresztmetszetre merőleges beáramlást feltételezve, a fiktív sebesség és az általam előírt statikus nyomás segítségével számítottam ki a peremen érvényes össznyomást. Ezek után minden egyéb, a peremen érvényes paramétert a szokásos módon határoztam meg. A dinamikus nyomás ily módon történő visszacsatolásával addig módosult a bemeneti torlóponti nyomás, míg a szorosan összetartozó statikus nyomás- és sebességmező a
pst bemeneti perem pto bemeneti perem
72
33. ábra. Statikus és torlóponti bemeneti peremfeltételek hatása az állandó nyomás görbékre
megkívánt statikus nyomáshoz tartozó értékre állt be (32. ábra). Az időiterációk során a kialakuló áramlás belépő sebességétől függően fog változni az össznyomás, így több időt vesz igénybe az időben beállt folyamat elérése, de a 32. ábrának megfelelően a peremen előálló statikus nyomást sikerült a megkívánt értékre beállítanom. A validálást elvégezve minimális a különbség a hagyományos és az újfajta számítási eljárás eredményei között (33. ábra, a 34. és 35. ábra alapján (a validációs eljárás ismertetése a következő második bekezdésben található)). Az időiterációk során már a 2.3. alfejezetben ismertetett negyedrendű Runge-Kutta módszert alkalmaztam a következő α k paraméterekkel:
α 1 = 1 4 , α 2 = 1 3 , α 3 = 1 2 és α 4 = 1.
(3.12)
A módszer validálását a 2.5.2 alfejezetben leírt teszteset segítségével végeztem el. A saját megoldást a FLUENT kereskedelmi program eredményeivel hasonlítottam össze. A geometria, a hálópontok mérete és elhelyezkedése, valamint a peremfeltételek megadása
A
B
34. ábra. Állandó statikus nyomás görbék (A: saját eredmény, B: FLUENT)
mindkét esetben ugyanaz volt (2.5.2. alfejezet, utolsó bekezdés). A nem viszkózus 2D-s áramlás-modellezés eredményei igen jó egyezőséget mutatnak (34, 35. ábra). A saját és a Fluent programok statikus nyomás-eredményei közötti eltérés számszerűsítésére a ∆pi , j ∆pi , j = ABS pi , j
(3.13)
összefüggést, vagyis a cellánkénti normált eltérés abszolút értékének kifejezését használtam. Az eredmények a 36. ábrán láthatók. Az eltérés minimális, a profil be-, illetve kilépő élénél tapasztalható a legjelentősebb különbség. 73
További összehasonlító elemzéseket végeztem a saját és a Fluent programok eredményei közötti eltérés vizsgálatára a statikus nyomás átlagos abszolút normált eltérése: 1 ∆p = Np
Np
∆p i pi
∑ ABS i =1
(3.14)
35. ábra. Állandó statikus nyomás görbék (folytonos vonal: saját eredmény, szaggatott vonal: FLUENT)
és a statikus nyomás normált L2 normája:
∆p p
= log 10
1 Np
∆pi ∑ i =1 p i Np
2
(3.15)
36. ábra. A cellánkénti statikus nyomás normált eltérésének abszolút értéke: saját és a Fluent program között
segítségével. A 4. táblázatba foglalt eredmények jól tükrözik a korábban már megfogalmazott állítást, miszerint elhanyagolható a különbség a saját és a Fluent program eredményei között. ELTÉRÉS
SAJÁT-FLUENT
∆p ∆p
1,633E-5
p
74
-4,2886
4. táblázat. Eltérés a CFD statikus nyomás között
programok eredményei
3.2. Konvergencia gyorsító abszorpciós szilárd fal peremfeltétel A hiperbolikus típusú Euler-egyenletek fontos tulajdonsága a hullámterjedés jelensége, vagyis a karakterisztikák menti információterjedés. A karakterisztikus görbesereg (állandó együtthatók esetén egyenes-sereg) meredeksége a Jacobi-mátrix sajátértékei segítségével határozható meg, mialatt a függő változók karakterisztikus formában vannak jelen az egyenletekben. Ezek a változók állandóak a hozzájuk tartozó karakterisztikus görbén, vagyis amelyeken az információ terjed (2.1.3. alfejezet). Az un. szimmetria (tükör vagy reflexiós) módszer szilárd fal perem alkalmazása esetén (2.4.2. alfejezet) ezek az információs hullámok vagy zavarások visszaverődnek. A korrekt kitűzésű numerikus problémák esetén akkor áll elő az időben beállt folyamat megoldása, miután az összes zavarás (azon numerikus hibák összessége, amelyeknél az időiterációk során a karakterisztikus változók állandósága nem biztosított a hozzájuk tartozó karakterisztikus görbén) elhagyta a kilépő peremen keresztül a számítási teret. Egy új numerikus peremfeltétel számítási eljárást dolgoztam ki az összenyomható áramlásra kidolgozott [5] alapján, amelynek használatával jelentősen csökkenthető az idő iterációk száma és a számítási idő. Az abszorpciós fal perem-feltétel lényege egy rugócsillapító lengőrendszer (38. ábra), amely az időiterációk során megengedi a falra merőleges nem zérusértékű sebesség komponensek megjelenését, mialatt a zavarások a falon keresztül, visszaverődés nélkül távozhatnak a számítási térből.
3.2.1. Alapegyenletek Az egy dimenziós Euler-egyenletek a (3.3) egyenlethez hasonlóan, pszeudokompresszibilitási alakban és a konzervatív formában a következőképpen írhatók fel:
∂U ∂U = 0, +A ∂x ∂t
(3.16)
0 β 2 P p A= , U = , és P = . ρ u 1 2u
(3.17)
λ(1,2 ) = u ± u 2 + β 2 .
(3.18)
ahol
Az A mátrix sajátértékei:
Az A mátrix sajátvektorai, pedig a következőképpen állíthatók elő: v v l (α ) A = λ(α )l (α ) ,
(3.19)
(l ( ) ,l ( ) )0
(3.20)
1
2
1
β 2
=λ 2u
(l
(1,2 ) (1)
)
,l (2 ) vagy
v l (1) 1 u + c L = v (2 ) = , l 1 u − c
(3.21)
amelyben a c = u 2 + β 2 , az úgynevezett pszeudo hangsebesség. Ekkor a karakterisztikus változók előállítása a (2.35) egyenlethez hasonlóan következőképpen történik: ∂W = L∂U ,
(3.22) 75
illetve
∂W (1) = ∂P + ∂u( u + c ) és ∂W (2 ) = ∂P + ∂u( u − c ) .
(3.23)
∂x = u + c meredekséggel a W (1) ∂t ∂x karakterisztikus változó, míg a befelé jövő karakterisztikus görbén, a =u−c ∂t meredekséggel a W (2 ) karakterisztikus változó állandó. A teljes differenciál fogalma alapján:
A számítási térből kifelé mutató karakterisztikus görbén a
dφ =
∂φ ∂φ dx + dt ∂x ∂t
(3.24)
a (3.23) első egyenletéből kiindulva felírható: ∂x 1. A = u + c karakterisztikus görbén ∂P + ∂u( u + c ) = 0 , vagyis ∂t
∂P ∂u ∂P ∂u + (u + c ) + (u + c ) + (u + c ) = 0 , ∂t ∂t ∂x ∂x
(3.25)
illetve a 37. ábrának megfelelő diszkretizálás elvégzése után PB − P1 u B − u1 Ps − PR u s − u R = 0 . ( ) ( ) ( ) + u + c + u + c + u + c R R R R R R ∆x 2 δt ∆ x 2 δt ( k )
(3.26)
A B pont körüli lineáris interpoláció
PB = Psδ 2 + PRδ 1 , illetve u B = u s δ 2 + u Rδ 1 ( δ 2 + δ 1 = 1 )
(3.27)
alkalmazását követően a (3.25) egyenlet a következő formában írható fel: S kifelé mutató karakterisztika
1
d1 u+c
DX/2
dt(k) B d2 R
37. ábra. Diszkretizációs séma az abszorpciós fal peremfeltételre
PS + (u R + c R )u S = PR + (u R + c R )u R −
δt ( k ) (u R + c R )[{(Psδ 2 + PRδ 1 ) − P1 } + ∆x 2
+ (u R + c R )[{(u s δ 2 + u R δ 1 ) − u1 }] .
(3.28)
Ekkor a következő változókat vezettem be:
δt (k ) (u R + c R ) és a4 = ∆x 2 76
(3.29)
a5 =
(1 − a 4 δ )(P + (u 1
R
R
+ c R )u R ) + a 4(P1 + (u R + c R )u1 ) . 1 + a4 δ 2
(3.30)
Az egyenletben szereplő összes paraméter ismert, mivel egyrészt a belső áramlási térhez tartoznak (1), másrészt a peremen érvényes előző időlépes eredményei (R) (37. ábra). Ezek után a (3.25) egyenlet a következőképpen írható fel:
PS + (u R + a R )u S = a5 .
(3.31)
2. A tömeg nélküli rugós-csillapító lengőrendszer (38. ábra) probléma specifikus differenciálegyenlete a következő formában fogalmazható meg [5]:
∂u u µ ∂P + = , ∂t τ ρc ∂t
(3.32)
ahol c a pszeudo-hangsebesség, u a falra merőleges sebességkomponens, P a (3.3) egyenletben megismert változó ( P = p / ρ ), k = 1 / µτ a merevségi és az α = 1 / µ pedig a csillapítási tényező. A τ és a µ paramétereket a 3.2.2. alfejezetben leírtak alapján úgy választottam meg, hogy minimálisra csökkentsék a visszavert zavarások számát. A 37. ábra alapján elvégeztem a (3.32) egyenlet diszkretizálását: valóságos szilárd fal
a=1/ m
k=1/ mt
numerikus fal perem
38. ábra. Rugós csillapító lengőrendszer us − uR us + uR µ Ps − PR , + = (k ) 2τ ρc R δt ( k ) δt
(3.33)
δt ( k ) δt ( k ) µ 1 + u s − 1 − u R = (Ps − PR ) . 2τ 2τ ρc R
(3.34)
δt ( k ) δt ( k ) µ a1 = 1 + , a 2 = ρc és a3 = 1 − 2τ u R − a 2 PR . 2 τ R
(3.35)
vagy
Legyen:
Az egyenletekben szereplő összes paraméter ismert, mivel egyrészt a belső áramlási térhez tartoznak (1), másrészt az előző időlépes eredményei (R). Ezek után a (3.33) egyenlet a következőképpen írható fel:
uS =
a 2 PS + a3 . a1
(3.36) 77
Végezetül az (3.31) és a (3.36) egyenletek segítségével a perem PS és u S következő (S) időlépésbeli paraméterei könnyen előállíthatók:
PS =
a5 a1 − ( u R + c R ) a3 és a1 + a 2(u R + c R )
(3.37)
a 2 PS + a3 . a1
(3.38)
uS =
3.2.2. A τ és a µ paraméterek meghatározása A τ és a µ paramétereket az egydimenziós akusztikus egyenletek segítségével határozhatók meg [5]: 1 ∂u ∂p + = 0 és ρ0 ∂t ∂x
(3.39)
∂p ∂u + ρ 0c02 = 0. ∂t ∂x
(3.40)
(A ρ 0 és c 0 vonatkoztató paraméterek a dimenziómentesítés következtében jelentek meg az egyenletekben [5]). A (3.32) egyenlet alapján a 1 µ ∂P ∂u u + = ∂t τ ( 1 + µ ) ( 1 + µ ) ρ 0 c0 ∂t egyenletet felírva, továbbá a τ ' = (1 + µ )τ és a µ' =
µ
(1 + µ )
(3.41) paraméterek bevezetése, a
ω P( x ) = p 0 + p sin x kezdeti feltétel megadása, az egyenlet integrálása, valamint a c komplex sebesség valós részének és az u( 0,0 ) = 0 sebességre vonatkozó kezdeti érték figyelembevétele után a következő egyenlet írható fel [5]:
(ωτ ' )µ' p (cos(ωt ) + (ωτ ' ) sin(ωt ) − e −t / τ ' ) . ρ 0 c 0 (1 + (ωτ ' )2 ) valamint a p ( x ,0 ) = P ( x ) kezdeti feltételekkel felírt
u( 0 , t ) = −
A már ismert u = 0 , karakterisztikus egyenlet, vagyis a p + ρ 0 c 0 u = állandó miatt:
p(0 ,t ) + ρ 0 c0 u (0 ,t ) = p(− c 0 t ,0 ) + 0 = P(− c0 t )
(3.42) kifelé menő (3.43)
egyenlet a következő összefüggésre vezet [5]:
p( 0,t ) = p 0 − p sin(ωt ) +
(ωτ ' )µ' p (cos(ωt ) + (ωτ ' ) sin(ωt ) − e −t / τ ' ) . 2 1 + (ωτ ' )
(3.44)
A nyomás és a sebesség azonos ω frekvenciával szinuszosan változik. Legyen S a visszavert és a beérkező hullám amplitúdó-viszonyának mérőszáma:
S = s r si , 78
(3.45)
amelyben: s r (0 ,t ) = p − ρ 0 c0 u = ( p (0 ,t ) − p 0 ) − u (0 ,t ) ,
(3.46)
s i (0 ,t ) = p + ρ 0 c0 u = p .
(3.47)
illetve:
Az időben állandósult megoldás esetén S értéke a következőképpen írható fel [5]:
[2(ωτ ' )µ' ]2 + [2(ωτ ' )2 µ' −1 − (ωτ ' )2 ] 2 1 + (ωτ ' )
2
S =
,
(3.48)
vagy matematikai átalakítások után: 2
S =
1 + (ωτ ) (1 − µ ) . 2 2 1 + (ωτ ) (1 + µ ) 2
2
(3.49)
A (3.49) egyenletben S értéke kizárólag µ -től és ωτ értékétől függ. Magasabb frekvenciák
esetén (ωτ >> 1) :
S =
1− µ . 1+ µ
(3.50)
A fenti egyenletekből következően a hullám-visszaverődés µ = 1,0 esetén minimalizálható. τ értéke egyrészt nem lehet jelentősen kisebb, mint az iteráció időlépése, mert ellenkező esetben problémákat okozhat a levágási hibából adódó pontatlanság. Másrészt, τ értéke nem lehet nagyobb annál az idő-intervallumnál, mialatt a következő hullám a peremhez érkezik. Ezek figyelembevételével és az [5] ajánlásával τ ∆t értékét 100-ra választottam.
3.2.3. Konzisztencia vizsgálat A konzisztencia vizsgálat célja, hogy t→∞ esetén valóban teljesül e az a feltétel, miszerint az időben beállt folyamatnál a szilárd fal peremre merőleges sebességkomponens zérus értékű. A (3.33) egyenlet us − uR us + uR µ Ps − PR + = (k ) 2τ ρc R δt ( k ) δt
(3.33)
alapján, amelyben R = S , azaz az időben beállt folyamat esetén uS + uS = 0; ⇒ u S = 0 2τ
(3.51)
valóban teljesül az a feltétel, hogy a falra merőleges sebességkomponens zérus. Az Euler alapú megoldók másik fontos tulajdonsága, hogy az időben beállt folyamat esetén ugyanakkora a nyomás a szilárd fal peremen, illetve a számítási térben vele szomszédos pontban. Ez a feltétel azonban a (3.28) egyenlet
δt ( k ) (u R + a R )[{(Psδ 2 + PRδ 1 ) − P1 } + PS + (u R + a R )u S = PR + (u R + a R )u R − ∆x 2 + (u R + a R )[{(u s δ 2 + u R δ 1 ) − u1 }]
(3.52)
79
esetén nem teljesül. Az időben beállt folyamat, vagyis a R = S feltétel esetén az előző egyenlet a következőképpen írható fel: 0 = −a4 δ 1 (PS + (u R + a R )u S ) − a4 δ 2 (PS + (u R + a R )u S ) + a4(P1 + (u R + a R )u1 ) ,
(3.53)
vagyis 0 = −a4 (PS + (u R + a R )u S ) + a4(P1 + (u R + a R )u1 ) .
(3.54)
A (3.54) egyenlet esetén Ps = P1 feltétel csak a (3.30) egyenletbeli u1 = u R helyettesítés esetén teljesül, ezért ennek megfelelően alakítottam át az egyenletet.
3.2.4. Időlépések és a konvergencia lefutás vizsgálata A negyedrendű Runge-Kutta módszer alkalmazása miatt az egyenletekben szereplő nem áramlástani paramétereket a következőképpen számítottam ki [5]:
δt ( i ) = ( 1 + α i )
∆t 2
, δ 1,i =
αi 1 1 , δ 2 ,i = , valamint i = 1, 2 , 3, 4 és α i = . 1+ αi 1+ αi imax + 1 − i
A konvergencia gyorsító abszorpciós és a hagyományos szimmetria (2.4.2. alfejezet) szilárd fal peremfeltétel alkalmazások összehasonlítására a körprofilt magába foglaló 2D-s csatorna (2.5.2. alfejezet) tesztesetet használtam. A kezdeti feltételek: v = 0 m s , u = 2 m s,
p st = 98000 Pa , illetve a be- és kilépő peremfeltételek:
α in = 0° ,
st pinto = 100000 Pa és p out = 98000 Pa mindkét esetben ugyanazok az értékek voltak. A
konvergencia kritériumot a (2.148) egyenlethez hasonlóan, de a nyomás reziduum normált L2 normájával határoztam meg:
∆p p
= log 10
1 Np
∆pi ∑ i =1 p i Np
2
(3.55)
Nyomás reziduum L2 normája
és értéke mindkét esetben -7 volt. A programok futtatására a 2.4.3. alfejezetben ismertetett Compaq szuperszámítógépet használtam. A konvergencia gyorsító abszorpciós szilárd fal peremfeltétel alkalmazása esetén a program-futás időszükséglete 40,5 %-a volt a szimmetria fal peremfeltétel alkalmazás időszükségletének. Az iterációszám pedig kb. 43%-kal csökkent szimmetria fal peremfeltétel alkalmazásához képest (39. ábra). A peremfeltétel validálását a 2.5.2. alfejezetben már ismertetett módszerrel, a 3.1.1. Iterációszám 0
10000
20000
30000
-4 -4,5 -5 -5,5 -6 -6,5 -7 -7,5
Abszorpciós szilárd fal modell Szimmetria szilárd fal modell
39. ábra. Konvergencia lefutások
80
40000
alfejezetben leírtakhoz hasonlóan végeztem. A szimmetria és az abszorpciós fal peremfeltételek eredményeinek összehasonlításakor minimális eltérést tapasztaltam az időben beállt folyamat elérését követően (40. ábra).
szimmetria szilárd fal peremfeltétel abszorpciós szilárd fal peremfeltétel
40. ábra. Peremfeltétel alkalmazások eredményeinek összehasonlítása
3.3. Tervezési irányelvek kidolgozása a fordító típusúhoz hasonló folyadéksugár-szivattyúk numerikus optimalizálására A 3.1. alfejezetben tárgyalt numerikus eljárás segítségével irányelveket dolgoztam ki, amelyek figyelembevételével növelhető a folyadékszállítás a fordító típusúhoz hasonló folyadéksugár-szivattyúkban. A 2D-s analízis és direkt numerikus optimalizálás egyben egy ipari-fejlesztői megbízás része volt [78]. Napjainkban, a gépjárművek tüzelőanyag tartályában elterjedten alkalmazzák a fordító folyadéksugár-szivattyút (41. ábra), amelynek feladata a tüzelőanyag eljuttatása a passzív nyeregtankból az aktívba, vagy az aktív nyeregtankból a központi kiürítő tartályba. A fordító beömlő csatorna
keverőtér, diffúzor
fúvóka
torok
szállítandó közeg belépése
tápfolyadék belépés össz-folyadék kilépése
41. ábra. 1,9 mm-s fúvóka átmérőjű alap modell fordító sugárszivattyú metszeti (bal oldal) és robbantott ábrája
elnevezés a tápoldali és a visszaszállítandó közeg vezetékében létrejövő egymáshoz képesti ellentétes irányú folyadékáramlásra utal. A cél minden esetben − a megfelelő üzembiztonság betartásával − a lehető legnagyobb szállítóképesség (a kilépő össz-térfogatáram és a tápoldali belépő térfogatáram hányadosa) elérése, amely az adott üzemi körülmények között a beömlési és a keverőtér optimális geometriai kialakításával érhető el. A számítógépes szimulációkat tekintve, a kétdimenziós numerikus modell és a hálózás elkészítését követően elvégeztem a numerikus modellek futtatási eredményeinek elemzését. A 81
A
B
42. ábra. Áramvonalak az alap modell fordító sugárszivattyúban. A keverőtér beömlő keresztmetszete A: letörés nélkül, B: letöréssel
térfogatáramokat egy, a méréseken alapuló 2D-3D konverziós táblázat segítségével számítottam ki [78]. Az előzetes vizsgálatok során arra a megállapításra jutottam, miszerint a kiindulási, alap geometria jobb oldali 90°-s fordítókönyökének külső ívén leválási buborék alakul ki (42/A ábra). Egy folyadékelemre vonatkoztatott, a nyomásból és a centrifugális erőtérből származó radiális irányú erőegyensúly a következő összefüggésre vezet: V 2 ∂p (3.56) ρ t = , Rc ∂n ahol a Vt a tangenciális irányú sebesség, Rc a görbületi sugár és n a kifelé mutató normális irány. A konkáv görbületi hatás miatt kialakult pozitív nyomásgradiens következtében a külső falhoz közeli rétegekben annyira lecsökkenhet a sebesség, hogy visszaáramlás jöhet létre. Ebben a speciális esetben azonban, a kialakult állapotért a keverőtér kör alakú belépő keresztmetszetének sarkos kialakítása is felelős. A fordítókönyök belépő keresztmetszetéhez érkező, a kontrakció hatásának megszűnése miatti széttartó áramlás okozta sebességcsökkenés + *
43. ábra. Sugárszivattyú alap geometria hálózása a letöréssel+ és a 0,5-mm hosszú áramlásterelő füllel*
hatása váltja ki lényegében a külső fal melletti lokális visszaáramlást. Összefoglalva, a konkáv görbületi hatás és a kontrakciós hatás megszűnésének szuperpozíciója miatt alakul ki a nem kívánt áramlási kép. A torokban alkalmazott letörés segítségével (42/B. ábra) sikerült megszüntetnem a leválást, aminek következtében javult a szállítóképesség (5. táblázat), a görbületi hatások okozta, a külső ív irányába történő áramvonal-ritkulás azonban nem szűnt meg. További numerikus áramlástani vizsgálatok segítségével sikerült bebizonyítanom, hogy a fúvóka kilépő keresztmetszeténél alkalmazott áramlásvezető fül (43. ábra) alkalmazásával sikerült úgy módosítani az áramlást, hogy az pozitív hatással legyen a szállítóképességre (5. 82
A
B
44. ábra. Sebességeloszlás [m/s] (fent) és áramvonalak (lent) a sugárszivattyúban; A: Alap geometria, B: Optimalizált geometria
táblázat). Ebben az esetben azonban a fúvóka és a torok közötti távolság jelentős csökkenése negatív irányban is változtathat a szállítóképességen. A további szállítóképesség növelésének érdekében, a minimális gyártósori változtatások figyelembevételével és az adott üzemi feltételek biztosítása mellett meghatároztam azokat a főbb geometriai méreteket, amelyek megfelelő méretkombinációja estén maximálható a szállítóképesség. Ezek rendre: a beömlő csatorna szélessége, a keverőtér torok átmérője és a keverőtér hossza. Az egyes geometriai jellemzők méretváltozatai a kiindulási konfigurációhoz közel, minden esetben ötféle ekvidisztáns segítségével kerültek meghatározásra. A direkt numerikus optimalizációhoz a [78]-ban kidolgozott programot használtam. A 2D-s analízisnek köszönhetően az 53=125 db számítás, a konvergencia kritériumtól függően, kb. 24 óra alatt lefutott. A folyamat során a legnagyobb szállítóképességű geometria (44/B ábra)
330,69
Alap+letörés+ 0,2 mm terelőnyelv 335,08
Alap+letörés+ 0,5 mm terelőnyelv 338,465
129,8
128,83
127,73
128,45
115,9
2,36
2,57
2,62
2,63
3,18
Modell
Alap
Alap+letörés
V&kiössz [l/h] V& táp [l/h]
307,0
be
ε=
V&kiössz V&betáp
Optimalizált 369,1
5. táblázat. A különféle geometriai modellek és szállítóképességeik
kiválasztásra került (5. táblázat). A 44. ábra sebességeloszlásán jól megfigyelhető, hogy az optimalizált modell esetén a fúvókából kiáramló nagy sebességű közeg mozgási energiájának nagyobb része fordítódott a beömlő csatorna felől érkező közeg gyorsítására, mint a kiindulási, alap geometria esetén. A fordítókönyökben kialakuló áramvonal eloszlás egyenletesebbé vált, így további javulást sikerült elérnem a szállítóképesség tekintetében, valamint tovább csökkent a leválási buborék kialakulásának valószínűsége. A 83
fordítókönyökben kialakult áramlástani viszonyok visszahatnak a torok és a fúvóka közötti térrészben lezajló folyamatokra is. A fordítókönyök belső ívében a konvex görbületi hatás miatt kialakult alacsony nyomásnak köszönhetően nagyobb mennyiségű és sebességű folyadék áramlik az alsó beömlő csatornából a diffúzorba. Ennek következményeként a fúvókát elhagyó áramlás nem vízszintes, hanem enyhén felfelé ívelő, aminek hatására az áramvonalak a fölső diffúzorfal irányába sűrűsödnek. Ez viszont jól kompenzálja a fordítókönyök külső felületén a lokális túlnyomás okozta áramvonal ritkulást. További vizsgálatokat lenne célszerű elvégezni a fúvóka excentricitásának, illetve irányítottságának hatásáról.
84
IV. Fejezet: Tézisgyűjtemény 4.1. Az összenyomható áramlások numerikus modellezése Numerikus eljárásokat dolgoztam ki és validáltam kétdimenziós, nem viszkózus és összenyomható áramlás numerikus modellezésére, amelyeknek kiinduló egyenleteit a konzervatív formában felírt Euler-egyenletek alkotják [30]: ∂ρ ∂ (ρu ) ∂ (ρv ) =0; + + ∂y ∂x ∂t
(
)
∂ (ρu ) ∂ ρu 2 + p ∂ (ρuv ) = 0; + + ∂y ∂x ∂t
(
)
) (
)
∂ (ρv ) ∂ (ρuv ) ∂ ρv 2 + p = 0; + + ∂y ∂x ∂t
(
∂ (ρE ) ∂ ρuhto ∂ ρvhto = 0; + + ∂y ∂x ∂t
1. A lökéshullámok nagypontosságú megoldásának (high resolution method), illetve az 1.1. alpontban kidolgozott újfajta peremfeltétel alkalmazása érdekében a számítógép által kezelhető formára dolgoztam át és validáltam Roe fluxus különbség megosztáson alapuló közelítő Riemann (upwind) módszerét [55] (2.2.1.2/b. alfejezet). (Megj.: a séma másod-, illetve magasabbrendű kiterjesztésekor a lökéshullámok közelében fellépő − a nem monoton tulajdonság okozta − oszcilláció elkerülésének érdekében Mulder-határolót alkalmaztam [44].) 1.1. Új, a karakterisztikák módszerén alapuló be- és kilépő peremfeltétel számítási eljárást dolgoztam ki, melyben a numerikus peremfeltételeket a cella normális irányú Machszám nemlineáris egyenletének segítségével számítottam ki (2.4.3. alfejezet). (Megj.: az áramlástani paramétereket leképeztem a cellából kifelé mutató egységvektor irányába. A szakadások nélküli belépő perem esetén a Riemann invariánsok állandóságának biztosítása érdekében a következő nemlineáris egyenletre jutottam: a( M B− 2 +
γ −1 2
)
−
1 2
− b1 +
γ −1
M B2 2
−
1 2
− k = 0;
A df ( M B ) dM B derivált meghatározása után a Newton−Raphson-eljárás segítségével határoztam meg a peremen érvényes M B Mach-számot. Az eljárás során csökkent a peremről visszavert zavarások intenzitása, amely gyorsabb konvergenciához vezetett.) 2. A folyadékok és a gázok áramlásának ugyanazon sémában történő modellezése és a transzszonikus áramlás egyik legjobb hatásfokú szimulációjának érdekében (egyszerű programozhatóság, kis gépidő, illetve kis számítógépi kapacitás melletti megfelelő pontosság), valamint a 2.1. alpontban leírt újfajta peremfeltétel alkalmazására egy cellaközpontú véges térfogat módszert dolgoztam ki (2.2.1.1. alfejezet). (Megj.: az eljárás időben állandósult ideális folyamatok modellezésére alkalmas, amelyben a negyedrendű Runge-Kutta módszer biztosítja – egy tetszőleges kezdeti állapotból kiindulva – az állandósult állapot elérését. A centrális diszkretizáció okozta oszcilláció elkerülését, illetve a szakadási felület közelében jelentkező pontatlanságok megszüntetését egy beépített numerikus szenzor segítségével Jameson [35] által kidolgozott mesterséges disszipáció biztosítja.) 2.1. Az összenyomható közegek numerikus áramlástani modellezése során új, a transzszonikus axiálkompresszorok forgó lapátsorának 2D-s matematikai 85
modellezésébe integrálható olyan bemeneti peremfeltétel számítási eljárást dolgoztam ki, amely figyelembe veszi a forgórész fordulatszámát és a futólapátozás mögött kialakult nyomásnövekedést, ezáltal teremtve meg a kapcsolatot az abszolút és relatív jellemzők között (2.6. alfejezet). (Megj.: ennek értelmében az általában előírt ( T to , p to , α ), és az ismert ( p ) bemeneti peremfeltételek közül kettő nem független, vagyis a kialakult nyomásnövekedés és a fordulatszám határozza meg a tömegáramot, illetve a belépő áramlás irányát. Az eljárás validálását egy transzszonikus áramlásba helyezett 2D-s DCA ((Double Circular Arc) kettős körívből álló) profilú lapátrács segítségével végeztem el. A módszer alkalmazásaként meghatároztam egy előperdítés nélküli DCA profilú rotor karakterisztikáját. A numerikus analízis során a 2D-s lapátrácsot a lapátozás külső felületén elhelyezkedő áramfelület síkba vetítésével állítottam elő.) 3. Modern kereskedelmi és a VKI-ban kifejlesztett numerikus-tervezési és optimalizálási programok felhasználásával a többfokozatú centrifugálkompresszor két fokozatát összekötő csatorna lapátozás-tervezésének újfajta irányelveit határoztam meg és dolgoztam ki, amelyet egy gyakorlati alkalmazáson keresztül mutattam be (2.7. alfejezet): 3.1. A számítógépes modellezés során egy új, állandó lapátterhelésen alapuló (ÁLT) lapáttervezési módszert fejlesztettem ki (2.7.3. alfejezet). (Megj.: Az analitikus eljárásban zérusértékű cirkuláció feltételezéséből indultam ki, amelyet a szívott és a nyomott oldal, valamint két állandó sugarú kör által közrezárt ellenőrzőfelületre írtam fel. Az állandó lapátterhelés feltételét a lapát két oldalán kialakult sebességkülönbség állandóvá tételével értem el:
δth d 2π Wps −Wss = cosβbl − WmR tg (β fl ) = Ct ; z R cos (βbl ) dm
(
)
Ezek után numerikus-integrálással meghatároztam a lapátszöget a csatorna külső és belső oldalán. A lapátvastagság-eloszlást a vázvonal első 75%-ban egy ellipszis kis tengelyével párhuzamos metszékeivel modelleztem. A lapát másik felében lineáris eloszlást alkalmaztam a maximális vastagságtól az előírt kilépő él vastagságáig.) 3.2. Numerikus áramlástani vizsgálatok arra engedtek következtetni, hogy jelentős javulás érhető el a két tervezést meghatározó paraméter, a veszteségi tényező és a nyomásnövekedési tényező tekintetében:
ω=
to
to
to
st
p in − p out p in − p in
; Cp =
st
st
to
st
p out − p in p in − p in
;
ha a lapátok belépő élét kiterjesztem a lapát nélküli diffúzor kimenetéig. Ebben az esetben megszűnt a fordítókönyökbe történő belépéskor kialakult leválási buborék. (Megj.: a numerikus analízist minden esetben a CFX-TASCflow kereskedelmi szoftver segítségével végeztem, amely az áramlástani gépekben lezajlódó folyamatok modellezésének napjaink élenjáró eszköze.) 3.3. A további optimalizáció érdekében egy 3D-s inverz lapáttervezési eljárás segítségével tovább tudtam javítani a tervezési paramétereket (2.7.4. alfejezet). (Megj.: a szoftver a VKI-ban kifejlesztett numerikus-tervező eszköz [19]. A program az általam előírt Mach-szám eloszlás alapján úgy módosította a lapátprofilt, hogy sikerült megszüntetni a visszatérő ágban, a lapát szívott oldalán idáig minden esetben előálló leválást). 3.4. A további hatásfoknövelés érdekében elvégeztem a pozitív és negatív lapátelhajlítás másodlagos áramlás csökkentésére gyakorolt hatásvizsgálatát. Ennek eredményeként további javulást sikerült elérnem negatív lapátelhajlítás esetén a tervezési paraméterek tekintetében (2.7.6. alfejezet). 86
4.2. Az összenyomhatatlan áramlások numerikus modellezése 4. A folyadékok és a gázok áramlásának ugyanazon eljárásban történő modellezése, valamint Chorin módszerének [11] a Jameson-féle disszipációval [35] figyelembe vett és a véges térfogat modellre kidolgozott változatának megvalósítása érdekében numerikus eljárást dolgoztam ki és validáltam a kétdimenziós nem viszkózus, összenyomhatatlannak feltételezett áramlás numerikus modellezésére (3.1. alfejezet), amelynek kiinduló egyenleteit a konzervatív formában felírt Euler-egyenletek alkotják [30]:
(
∂u ∂v + = 0; ∂x ∂y
)
∂u ∂ u 2 + p ρ ∂ (uv ) + + = 0; ∂t ∂x ∂y ∂v ∂ (uv ) ∂ v 2 + p ρ + + =0; ∂t ∂x ∂y
(
)
(Megj.: az időben beállt folyamat elérésére a kis memóriaigényű negyedrendű Runge-Kutta módszert alkalmaztam. A numerikus séma stabilitását Jameson [35] mesterséges viszkozitás segítségével értem el.) 4.1. Újfajta, alkalmazásfüggő peremfeltételt dolgoztam ki olyan speciális feltételek esetére, amelyekben a bemeneti torlóponti nyomás helyett a statikus nyomás megadásával lehet a környezeti beáramlásokat szimulálni (3.1.1. alfejezet). 4.2. Kidolgoztam egy újfajta abszorpciós fal peremfeltétel számítási eljárást az összenyomhatatlannak feltételezett áramlások numerikus modellezésébe. A módszer alkalmazásával kb. 50 %-kal növelhető a konvergencia sebessége (3.2. alfejezet). (Megj.: az előzőekben tárgyalt új eljárások validálását a FLUENT kereskedelmi szoftver segítségével végeztem el.) 5.A 4. pontban kidolgozott numerikus eszköz segítségével új tervezési irányelveket határoztam meg és dolgoztam ki a fordító típusúhoz hasonló sugárszivattyúk numerikus optimalizálására. A tervezési eljárásrendszer kidolgozása egyben egy ipari-fejlesztői munka része volt [78], amelynek célja egy új termék előállítása, azonban az eljárás során megállapított irányelvek általános érvényűen alkalmazhatók minden hasonló típusú szivattyú tervezésében. A munka során meghatároztam azokat a fontosabb irányelveket, amelyek figyelembevételével növelhető a folyadékszállítás a hasonló típusú sugárszivattyúkban (3.3. alfejezet): 5.1. Az előzetes numerikus analízis során megállapítottam, hogy a folyadéksugár visszafordításakor a könyökben előálló leválási buborék a torok beáramlási keresztmetszetében alkalmazott letörés segítségével megszüntethető. 5.2. Bebizonyítottam továbbá, hogy a fúvóka külső felületének konfúzoros kialakítása minden esetben növelte a folyadékszállítást. 5.3. A direkt numerikus optimalizálás eredményeként egyértelműen meghatározható egy olyan szerkezeti kialakítás, amely a torok beáramlási keresztmetszet átmérőjének, a diffúzor hosszának, illetve a beömlőcsatorna méretének olyan kombinációja esetén áll elő, amikor maximális a szállítóképesség [78].
87
88
Irodalomjegyzék [1] ABEELE V. D., SNYDER D. O. and DEGREZ G.: A Combined Sectral/Finete Element Method for the Direct and Large Eddy Simulation of Turbulent Flows in Complex, TwoDimensional Geometries. Conference on Modelling Fluid Flow CMFF’03, Conference Proceedings, Volume II., ISBN 9634207790, pp. 783-790, 2003. [2] BALDVIN B. S. and BARTH T. J.: A One Equation Turbulence Transport Model for High Reynolds Number Wall-Bounded Flows. NASA TM-102847 and AIAA paper 910610, 1991. [3] BALDWIN B. S. and LOMAX H.: Thin-Layer Approximation and Algebraic Model for Separated Turbulent Flows. AIAA paper, 1978. [4] BALLHAUS W. F. AND BAILEY F. R.: Numerical Calculations of Transonic Flow about Sweept Wings. AIAA paper 72-677, 1972. [5] BODY K. M.: Numerical Wave Propagation and Steady State Solution, Ph.D. thesis, Aerospace Engineering and Scientific Computing, University of Michigan, 1992. [6] BOUSSINESQ J.: Théori de l’Écoulement Tourbillant. Mem. Présantés par Divers Savants Acad. Sci. Inst. Fr.., Vol. 23, pp. 46-50, 1877. [7] CEBECI T. and BRADSHAW P.: Momentum Transfer in Boundary layers. Washington DC: Hemisphere, 1984. [8] CEBECI T. and SMITH A. M. O.: Analysis of Turbulent Boundary Layers. Ser. In Appl. Math. And Mech., Vol XV, Academic Press, Orlando, Fl. 1974. [9] CHAKRAVARTHY S. R.: Euler Equations, Implicit Schemes and Boundary Conditions. AIAA Journal, 21:699-706, 1983. [10] CHAKRAVARTY S. R. and OSHER S.: A New Class of High Accuracy TVD Scheme for Hyperbolic Conservation Laws. AIAA paper No. 85-0363, 1985. [11] CHORIN A. J.: A Numerical Method for Solving Incompressible Viscous Flow Problems. Journal of Computational Physics, 2, 12-26. 1967. [12] CORSINI A., RISPOLI F. and SANTORIELLO A.: A New Stabilized Finite Element Method for Advection-Diffusion-Reaction Equations Using Quadratic Elements. Conference on Modelling Fluid Flow CMFF’03, Conference Proceedings, Volume II., ISBN 9634207790, pp. 791-799, 2003. [13] COURANT R. and HILBERT D.: Methods of Mathematical Physics. Volume I-II, WileyInterscience, New York, 1962. [14] CZIBERE T.: Three Dimensional Stochastic Model of Turbulence. Journal of Computational and Applied Mechanics, Vol. 2 pp. 7-20, 2001. [15] DAIGUJI H. and ZAMAMOTO S.: Numerical Methods for Transonic Cascade Flow Problems. A Collection of Technical Papers: International Symposium on Computational Fluid Dynamics-Nagoya, pp. 548-553, 1989. [16] DECONINCK H.: Upwind Methods and Multidimensional Splitting for the Euler Equations. VKI Lecture Series: Computational Fluid Dynamics, 1991-01, 1991. [17] DECONINCK H. and STRUYS R.: A Consistent Boundary Condition for Cell Centered Upwind Finite Volume Euler Solvers. Numerical Methods for Fluid Dynamic III, Clarendon Press, Oxford, 1988. [18] DEMEULENAERE A.: Conception et development d'une methode inverse pour la gmeration d'aubes de turbomachines, Ph.D Thesis VKI, 1997. [19] DEMEULENAERE A., and VAN DEN BRAEMBUSSCHE R.: Three-Dimensional Inverse Method for Turbomachinery Blading Design, ASME Journal of Turbomachinery, Vol. 120, 1998. [20] ERDŐDY I.: Véges elemes módszer az áramlástanban és a mechanikában. BME jegyzet, megjelenés alatt, 2002.
89
[21] FURUKAWA M. YAMASAKI M. and INOUE, M.: A Zonal Approach for NavierStokes Computations of Compressible Cascade Flow Fields Using a TVD Finite Volume Method. ASME Journal of Turbomachinery, Vol. 113, No. 4, pp. 573-582, 1991. [22] GAUSZ T.: Az örvénytranszport egyenlet közelítő megoldása a véges differenciák módszerével. Oktatási segédanyag, BME, Repülőgépek és Hajók Tanszék, Budapest, 2003. [23] GAUSZ T.: Hő- és áramlástan II, Egyetemi jegyzet, BME, Repülőgépek és Hajók Tanszék, Budapest, 2003. [24] GAUSZ T.: Aerodynamical and Dynamical Investigation of Helicopter Rotors. ICAS Congress, Harrogate, United Kingdom, ICA0182, 2000. [25] GERMANO M., PIOMELLI U., MOIN P. and CABOT W. H.: A Dynamic Subgrid-Scale Eddy Viscosity Model. Phys. Fluids, A3 (7), pp. 1760-1765. 1991. [26] HAFEZ M. M.: Progress in Finite Element Techniques for Transonic Flows. AIAA 6th Computational Fluid Dynamics Conference, pages 243-250, 1983. [27] HARTEN A. LAX P. D. and VAN LEER B.: On Upstream Differencing and Godunov Type Schemes for Hyperbolic Conservation Laws. SIAM rewiev, 25:35-61, 1983. [28] HIRSCH C.: Numerical Computation of Internal and External Flows. John and Wiley and Sons Volume I-II. ISBN 0471917621, 1988. [29] HOFFMANN J. D.: Relationship Between the Truncation Errors of the Centered Finite Difference Approximation on uniform and Non-Uniform Meshes. Journal of Computational Physics, 46:469-477, 1982. [30] HOFFMANN K. A., CHIANG S. T. L. SIDDIQUI M. S. and PAPADAKIS M.: Fundamental Equations of Fluid Mechanics. Engineering Education System, ISBN 09623731-9-2, 1996. [31] HOLST T. L.: An Implicit Algorithm for the Conservative Transonic Full Potential Equation Using an Arbitrary Mesh. AIAA Journal, Volume 17:1038-1045, 1979. [32] HOUTMANN E. M.: Numerical Simulation of Three-Dimensional Compressible Flows. Faculty of Aerospace Engineering, Delft University of Technology HSL TM 96381, 1996. [33] ISMERETLEN, Numerical Methods for the Incompressible Navier-Stokes Equations. Lecture Notes at the National Central University, Taiwan. Letöltve a w3new.ncu.edu.tw /~junwu/me681/chap6.doc internet címről, 2002. [34] JAMESON A.: Iterative Solutions of Transonic Flows Over Airfoils and Wings Including Flows at Mach 1. Comm. Pure and Applied Mathematics., 27:283-309, 1974. [35] JAMESON A. and SCHMIDT W.: Recent Developments in Finite Volume TimeDependent Techniques for Two and Three Dimensional Transonic Flow. VKI Lecture Series: Computational Fluid Dynamics, 1982. [36] JANIGA G.: 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. PhD értekezés Miskolci Egyetem Gépészmérnöki Kar, 2002. [37] KRISTÓF G., LAJOS T., LOHÁSZ M., and RÉGERT T.: Application of Numerical Simulation in Comfort Analyses and Fire Case Studies of Budapest Arena. KlímaváltozásEnergiatudatosság-Energiahatékonyság, Győr, 2003. [38] KRISTÓF G., VARGA L., és KAPÁS N.: A numerikus áramlástan néhány energetikai alkalmazása. Magyar energetika, 2003/2 [39] KRISTÓF G., PÖSZMET I., RÉGERT T., MEZŐSI B., and DÁVIDHÁZI N.: Numerical simulation of a side channel pump. The 12th International Conference on Fluid Flow Technologies (CMFF'03), Budapest, Hungary, 2003. [40] LAJOS T.: Az Áramlástan alapjai I-II. Egyetemi jegyzet, ISBN 9634205860ö, BME Áramlástan Tsz, 1999. [41] LAX P. D.: Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves. SIAM Publication, Philadelphia, 1973.
90
[42] LEVELQUE R. J.: Numerical Methods for Conservation Laws. Lectures at ETH Zurich, Birkhauser Verlag, Lecture Note Series, 1989. [43] LEWIS, G. W. JR.: Experimental Investigation of Axial-Flow Compressor Inlet Stage Operating at Transonic Relative Inlet Mach Numbers, II-Blade Coordinate Data. NACA Research Memorandum, RM E52C27, National Advisory Committee for Aeronautics Washington, 1952. [44] MANNA M.: A Three Dimensional High Resolution Compressible Flow Solver. PhD Thesis in Catholic University of Louvain, 1992. [45] MATSUO Y.: Computations of Three Dimensional Viscous Flows in Turbomachinery Cascades. AIAA Paper No. 91-2237, 1991. [46] MCDONALD P. W.: The Computation of Transonic Flow through Two Dimensional Gas Turbine Cascade. ASME paper 71-GT-89, 1971. [47] MENG S. Y. and JACKSON E. D.: The Continuous Diffusion Crossover System Design. ASME FED-vol3. Return Passages of Multistage Turbomachinery. [48] NYÍRI A.: Potential Flow around the Blades of Hydraulic Machines. JSME Symposion, Tokyo, Vol. 1. pp. 61-70, Vol. 4. pp. 44-45, 1972. [49] NYÍRI A.: Three Dimensional Potential Flow around the Blades of a Hydraulic Machine, 6th Symposion, Magdeburg, pp. 162-167. 1989 [50] NYÍRI A.: Vízgépek lapátozott terére vonatkozó háromdimenziós potenciálos áramlási feladat megoldása. Akadémiai doktori értekezés, pp. 1-115, Miskolc, 1990. [51] ORSZAG S. A. and PATTERSON G.S.: Numerical Simulation of Three-Dimensional Homogeneous Isotropic Turbulence. Phys. Rev. Lett., 28:76-79, 1972. [52] PRANDTL L.: Über die ausgebildete Turbulenz. ZAMM, Vol. 5, pp. 136-139, 1925. [53] RAO S. S.: The Finite Element Method in Engineering. Pergamon Press ISBN 0-08033419-9 Second Edition, 1988. [54] RIZZI A. W. and VIVIAND H.: Numerical Methods for the Computations of Transonic Flows with Shock Waves. Notes on Numerical Fluid Mechanics, 3, Vieweg, Braunschweig, 1981. [55] ROE P. L.: Approximate Riemann Solvers, Parameter Vectors and Difference Schemes. Journal of Computational Physics, 43:357-372, 1981. [56] ROGALLO R. S.: Numerical Experiments in Homogeneous Turbulence. NASA TM 81315, 1981. [57] ROTHSTEIN E.: Experimentelle und Theoretische Untersuchung der Strömungsforgänge in Ruekfuerkanalen von Radialverdichterstufen. Thesis, Aachen, Germany, 1984 [58] RUMSEY C., VAN LEER B. and ROE P. L.: A Grid Independent Approximate Riemann Solver with Applications to the Euler and Navier-Stokes Equations. In AIAA paper 910239, 1991. [59] RUMSEY C., VAN LEER B. and ROE P. L.: A Multidimensional Flux Function with Applications to the Euler and Navier-Stokes Equations. Journal of Computational Physics 105, 306-323, 1992. [60] SALAS M. D. and GUMBERT C. R.: Breakdown of the Conservative potential Equation. AIAA 23rd Aerospace Science Meeting, AIAA paper 85-0367, 1985. [61] SALAS M. D., JAMESON A. and MELNIKE R. E.: A Comparative Study of the non Uniqueness Problem of the Potential Equation. AIAA 6th Computational Fluid Dynamics Conference, pages 48-60, AIAA paper 83-1888, 1983. [62] SÁNTA I.: Transzszonikus kompresszor fokozatok. Oktatási segédanyag, BME, Repülőgépek és Hajók Tanszék, Budapest, 2002. [63] SCHLICHTING H.: Boundary Layer Theory. Series in Mechanical Engineering McGraw-Hill, 4th English Edition, 1979. [64] SHEN L. and YUE D. K. P.: Large-Eddy Simulation of Free Surface Turbulence. J. Fluid Mech., 440, pp. 75-116, 2001.
91
[65] SIMON H. and ROTHSTEIN E.: On the development of return passages of multistage centrifugal compressors. ASME FED-vol3. Return Passages of Multistage Turbomachinery. [66] SMAGORINSKY J.: General Circulation Experiments With the Primitive Equations. Mon. Weather Rev., 93, pp. 99-164, 1963. [67] SPALART P. R.: Direct Numerical Simulation of a Turbulent Boundary layer up to Rθ = 1410 . J. Fluid Mech., 187:61-98, 1988. [68] SPALART, P. R. and ALLMARAS S. R.: A One Equation Turbulence Model for Aerodynamic Flows. AIAA paper 92-439, Reno, NV, 1992. [69] STARKEN H.: Untersuchung der Strömung in eben Überschallverzögerungsgittern. DLRForschungsbericht, 1971. [70] STEGER J. L. and WARMING R. F.: Flux Vector Splitting of the Inviscid Gas Dynamic Equations with Application to Finite Difference Methods. Journal of Computational Physics, 40:263-293, 1981. [71] THYGESEN R.: Optimization of Return Channel Blades for Radial Compressors. ,VKI PR 2000-21, 2000. [72] TINOCO E. N. and CHEN A. W.: Transonic CFD Application to Engine/Airframe Integration. AIAA 22nd Aerospace Science Meeting, AIAA paper 84-0381, 1984. [73] VAN DEN BRAEMBUSSCHE R. A.: The Centrifugal Compressor Off-Design Analysis Program “CCOD”. VKI Course Notes 1995-06/B, 1999. [74] VAN DEN BRAEMBUSSCHE R. A.: Design and Optimisation of Centrifugal Compressor, VKI Course Notes 141, 1990. [75] VAN DEN BRAEMBUSSCHE R. A.: Szóbeli konzultáció, VKI Diploma Course, 2001. [76] VAN LEER B.: Flux Vector Splitting for the Euler Equations. 8th International Conference on Numerical Methods in Fluid Dynamics, Berlin Springer Verlag, 1982. [77] VERESS Á.: Inverse Design on Return Channel for Multistage Radial Compressor. VKI DC Project Report 2001-27, 2001. [78] VERESS Á. (témavezető), BÁNYAI T., BERKE P. and BITVAI I.: Numerical Simulation and Semi-Optimisation on 1.9 Orifice Fuel Jet-Pump. BME, DAS-VISTEON Project Report, 2003. [79] WILCOX D. C.: Turbulence Modelling for CFD. DCW Industries ISBN 0963605151, 1998. [80] YUN K. Y. and AGARWAL R. K.: Navier – Stokes and Burnett Simulations of Flows in Microchannels. AIAA Paper 2001-3045, AIAA 31st Fluid Dynamics Conference and Exhibit, Anaheim, CA, 2001. [81] ZANG Y., STREET R. L., and KOSEFF J. R.: A Dynamic Mixed Subgrid-Scale Model and Its Application to Turbulent Recirculating Flows. Phys. Fluids A5 (12), pp. 31863196, 1993.
92
