Foglalkoztatáspolitikai és Munkaügyi Minisztérium Humánerőforrás-fejlesztés Operatív Program
Dr. Kalmár László – Dr. Baranyi László – Dr. Könözsy László
Hő- és áramlástani feladatok numerikus modellezése
Készült
„A felsőoktatás szerkezeti és tartalmi fejlesztése” CAD/CAM/FEM kompetencia kurzusok projekt keretében
Miskolci Egyetem
1 NAVIER-STOKES EGYENLET MEGOLDÁSI LEHETŐSÉGEI
A Navier-Stokes-féle mozgásegyenlet numerikus megoldása lamináris, kis Reynolds számú áramlások vizsgálata esetén ugyan gyakran nehézségeket okoz, de még bonyolult áramlási tartományok esetén is megoldható. A gyakorlatban előforduló áramlások többsége viszont turbulens, amely áramlások számítása sokkal nagyobb kihívást jelölt a numerikus áramlástannal foglalkozóknak, mint a lamináris áramlásoké. Nem véletlenül tette 1994-ben a következő kijelentést Peter Bradshaw, a turbulencia kutatás egyik kiemelkedő alakja: „Turbulence is the invention of the devil on the 7th day of creation”, amely kissé szabad fordításban magyarul így hangzik, hogy „A turbulenciát az ördög találta ki a teremtés 7. napján”. Először mondjunk néhány szót a turbulens áramlások általános jellemzőiről. - Az ilyen áramlásra a folyadékrészecskék állandó keveredése, véletlenszerű mozgása a jellemző. Ilyen szempontból hasonlít a Brown-féle hőmozgáshoz. - A turbulens áramlás szigorúan véve mindig instacionárius és háromdimenziós. - A turbulens áramlásban általában mind az örvények mérete, mind az ingadozásra jellemző frekvencia nagyon széles skálán mozog. Egy nagy méretű és Reynolds számú szélcsatornában például az áramlás szerkezetére jellemző frekvenciára és ingadozási hullámhosszra a következőt találták: f turb = 1 ÷ 10000 Hz, illetve L=0,1 ÷ 4000 mm. Ilyen széles spektrumú áramlást nagyon nehéz vizsgálni úgy, hogy mind a térbeli és időbeli felbontás elég finom legyen. - A sebesség- és nyomásmérő műszerek általában átlagértéket mérnek. Kivételt képeznek e kijelentés alól a hődrótos anemométerek, amelyek ingadozást is képesek mérni. A turbulens áramlások vizsgálatára számos módszert dolgoztak ki. Ezek közül mi itt a teljesség igénye nélkül csak három fontos módszert kívánunk röviden ismertetni. Ezek az időátlagolt Navier-Stokes egyenlet, a Nagy Örvények Szimulációja és a Direkt Numerikus Szimuláció.
1.1 IDŐÁTLAGOLT NAVIER-STOKES EGYENLET (RANS) Időátlagolt Navier-Stokes egyenletre, amelyet Reynolds egyenletnek is neveznek, a nemzetközi szakirodalomban leginkább a RANS elnevezés terjedt el, amelyet a “Reynolds Avaraged Navier-Stokes” angol kifejezés kezdőbetűiből képeztek. A mérnöki gyakorlatban jelenleg ez a legáltalánosabban használt módszer turbulens áramlások vizsgálatára. Természetesen a magasabb szinteket képviselő Nagy Örvények Szimulációja és a Direkt Numerikus Szimuláció használata is kezd mind általánosabbá válni. A praktizáló mérnökök számára általában az időben átlagolt értékek a fontosabbak; számukra az ingadozás kevésbé érdekes. Mint ismeretes, az 1. fejezetben származtatott Navier-Stokes egyenlet nemcsak lamináris áramlásokra, hanem a turbulens áramlás pillanatnyi sebességére és nyomására is érvényes. Jelöljük a pillanatnyi sebesség vektorát és a pillanatnyi nyomást az alábbi módon: v T ( x , y , z ,t );
pT ( x, y,z,t ) .
(6.1)
Átlagoljuk most ezeket a pillanatnyi értékeket egy olyan T időtartamra, amely jóval hosszabb, mint a turbulens áramlás ingadozására jellemző időállandó, de T sokkal kisebb, mint az instacionárius áramlásra jellemző időállandó. Az átlagolás után kapjuk a v = v ( x, y,z,t ) átlagsebességet és p = p ( x, y,z,t ) átlagos nyomást: T
v=
1 v T dt , T ∫0
T
p=
1 pT dt . T ∫0
(6.2)
[2] szerint ezt az időt célszerű T ≈ 5 s nagyságúra választani. Ezek alapján a (6.1) egyenlettel megadott pillanatnyi sebesség és a pillanatnyi nyomást felbontható egy időbeli átlag és egy ingadozás összegére [4]:
vT ( x, y,z,t ) = v ( x, y,z ) + v ′ ( x, y,z,t ) ,
(6.3)
pT ( x, y,z,t ) = p ( x, y,z ) + p ′ ( x, y,z,t ) ,
(6.4)
ahol v ′ ( x, y,z,t ) , illetve p ′ ( x, y,z,t ) a turbulens sebességingadozás ill. nyomásingadozás. Írjuk most fel a különböző sebességek közötti összefüggést skaláregyenletek segítségével is: vTx = vx + u ′, ,
(6.5)
vTy = v y + v′, ,
(6.6)
vTz = vz + w′,
(6.7)
A v ' és a p ′ ingadozások tehát a (6.3) és (6.4) egyenletek alapján az alábbi módon számíthatók: v ′ ( x, y,z,t ) = vT ( x, y,z,t ) − v ( x, y,z,t ) ,
(6.8)
p ′ ( x, y,z,t ) = pT ( x, y,z,t ) − p ( x, y,z,t ) .
(6.9)
Írjuk fel most az 1. fejezetben származtatott, összenyomhatatlan közegre vonatkozó (1.30) Navier-Stokes egyenletet a pillanatnyi vT sebességvektor és pT nyomás figyelembe vételével: ∂ vT 1 + ( v T ⋅ ∇ ) v T = f − ∇pT + ν∆ v T . ∂t ρ
(6.10)
Az 1. fejezetben bemutatott (1.6) kontinuitási egyenletet írjuk most fel összenyomhatatlan közegre
div vT =
∂ vxT ∂ v yT ∂ vzT + + =0. ∂x ∂y ∂z
(6.11)
Mivel az átlagértékekre is hasonlóan teljesül a kontinuitási egyenlet, azaz a div v =
∂vx ∂v y ∂vz + + =0 ∂x ∂y ∂z
(6.12)
összefüggés felírható, így a (6.3) egyenlet alapján a sebességingadozásokra is hasonló alakban kapjuk a kontinuitási egyenletet: div v ′ =
∂ u ′ ∂ v ′ ∂ w′ + + =0. ∂x ∂y ∂z
(6.13)
Helyettesítsük most be az átlagértékek és ingadozások összegeként előállított pillanatnyi értékek (6.3) és (6.4) alakjait a (6.10) Navier-Stokes egyenletbe: ∂ v ∂ v′ 1 1 + + ⎡⎣( v + v ′ ) ⋅ ∇ ⎤⎦ ( v + v ′ ) = f − ∇p − ∇p ′ + ν∆ v + ν∆ v ′ . ∂t ∂t ρ ρ
A fenti egyenletben a konvektív tagokra vonatkozólag elvégezve a kijelölt műveleteket a következőt kapjuk:
∂ v ∂ v′ + + ( v ⋅∇ ) v + ( v ⋅ ∇ ) v ′ + ( v ′ ⋅ ∇ ) v + ( v ′ ⋅∇ ) v ′ = ∂t ∂t 1 1 = f − ∇p − ∇p ′ + ν∆ v + ν∆ v ′.
ρ
(6.14)
ρ
A (6.14) egyenletben a következő tagokat részletesen is kiírjuk:
( v ⋅∇ ) v′ = vx
∂ v′ ∂ v′ ∂ v′ + vy + vz , ∂x ∂y ∂z
( v′ ⋅ ∇ ) v = u ′
∂v ∂v ∂v , + v′ + w′ ∂x ∂y ∂z
( v′ ⋅ ∇ ) v′ = u ′
∂ v′ ∂ v′ ∂ v′ + v′ + w′ . ∂x ∂y ∂z
Ezek után átlagolni fogjuk a (6.14) egyenletet, de előtte nézzük meg az átlagolás szabályait két tetszőleges f és h függvényekre: f +h = f +h, fh = f h, cf = c f , ahol c egy állandó érték.
Az átlagolást az itt említett szabályoknak és fejezet elején leírtaknak megfelelően végezzük el, figyelembe véve, hogy az ingadozások időátlaga 0 lesz. Néhány számítási részlet: ∂v ∂v = , ∂t ∂t ∂v ′ = 0, ∂t
( v ⋅∇ ) v = ( v ⋅ ∇ ) v , v′ = 0 , ( v ⋅ ∇ ) v′ = ( v ⋅ ∇ ) { =0
⎛
⎞
⎝ =0
⎠
v ′ ⋅∇ ⎟ v = 0 , ( v′ ⋅ ∇ ) v = ⎜ { 1
ρ
∇ ( p + p′) =
1
ρ
∇p ,
ν∆ ( v + v ′ ) = ν∆ v . Ezen időátlagolt tagokat visszaírjuk az időátlagolt (6.14) egyenletbe, akkor a következő egyenletet kapjuk:
∂v 1 + ( v ⋅ ∇ ) v = f − ∇p + ν ∆ v − ( v ′ ⋅ ∇ ) v ′ . ρ ∂t
(6.15)
A −( v′ ⋅ ∇ ) v ′ tagot az elemi vektoranalízis felhasználásával a következő alakba írhatjuk át:
−( v ′ ⋅∇ ) v ′ =
1
ρ
Div σ T ,
(6.16)
ahol σΤ az úgynevezett Reynolds-féle turbulens feszültségtenzor, amely a következő alakú:
⎡ u ′2 u ′v ′ u ′w′⎤ ⎢ ⎥ σ T = − ρ ( v ′ o v ′ ) = − ρ ⎢ v ′u ′ v ′2 v′w′ ⎥ . ⎢ ⎥ ⎢ w′u ′ w′v ′ w′2 ⎥ ⎣ ⎦
(6.17)
Mivel a szorzatokban a tényezők sorrendje nem befolyásolja az elemek értékét; u ′v′ = v ′u ′ , u ′w′ = w′u ′ , stb, így a (6.17) egyenletre tekintve látható, hogy a tenzor szimmetrikus, mivel a főátlójára szimmetrikusan elhelyezkedő elemek megegyeznek egymással. Így tehát a tenzor 6 független elemet tartalmaz. A Reynolds-féle turbulens feszültségtenzort szokás látszólagos feszültségtenzornak is nevezni, mivel elemei nem valóságos feszültségek, hanem a turbulens impulzuscseréből származó erőhatásoknak a folyadék felületére lokalizálva, [4]. A (6.16) összefüggés segítségével a (6.15) egyenletből a következőt kapjuk: ∂v 1 1 + ( v ⋅ ∇ ) v = f − ∇p + ν ∆ v + Div σ T , ∂t ρ ρ
(6.18)
amely nem más, mint a Reynolds-féle mozgásegyenlet. A (6.16) és (6.18) egyenletekben a Div a tenzor divergenciáját jelenti, amelyet úgy számítunk ki, hogy a tenzort jobbról skalárisan megszorozzuk a nabla operátorral: ⎡ u ′2 u ′v ′ u ′w′ ⎤ ⎢ ⎥ ⎡ ∂ ∂x ⎤ 1 1 2 ⎢ Div σT = − ρ ( v ′ o v ′ ) ⋅∇ = −( v ′ o v ′ ) ⋅∇ = − v′u ′ v′ v′w′ ⎥ ⋅ ⎢⎢∂ ∂y ⎥⎥ ⎢ ⎥ ρ ρ ⎢ w′u ′ w′v′ w′2 ⎥ ⎢⎣ ∂ ∂z ⎥⎦ . ⎣ ⎦ A tenzor divergenciáját kiszámítva egy vektort kapunk. Ezt behelyettesítve a (6.18) Reynolds egyenletbe, felírhatjuk Reynolds egyenlet koordináta-egyenleteit: ∂vx ∂v ∂v ∂v + vx x + v y x + vz x = ∂t ∂x ∂y ∂z 1 ∂p ∂ ∂ ∂ = fx − + ν∆vx − u ′2 − u ′v ′ − u ′w′ , ∂x ∂y ∂z ρ ∂x
( ) ( ) ( )
∂v y
∂v y
∂vz ∂v ∂v ∂v + vx z + v y z + vz z = ∂t ∂x ∂y ∂z 1 ∂p ∂ ∂ ∂ u ′w′ − v ′w′ − w′2 . = fz − + ν∆vz − ρ ∂z ∂x ∂y ∂z
(6.21)
∂y
+ vz
∂v y
(6.20)
∂x
+ vy
∂v y
= ∂z ∂ ∂ ∂ 1 ∂p v′u ′ − v ′2 − v′w′ , = fy − + ν∆v y − ρ ∂y ∂x ∂y ∂z
∂t
+ vx
(6.19)
( )
( )
( ) ( ) ( )
( )
A fenti három koordináta-egyenlethez hozzávéve az átlagsebességekre felírt (6.12 ) kontinuitási egyenletet, összesen 4 egyenletünk van. Számoljuk most össze az ismeretleneket: háromdimenziós esetben a v átlagsebesség 3 komponense, a p nyomás, és a turbulens feszültségtenzor 6 független eleme, ez összesen 10 ismeretlent jelent a 4 egyenlettel szemben. Kétdimenziós esetben: 2 sebességkomponens, a nyomás, és a turbulens feszültségvektor 3 független eleme, azaz 6 ismeretlen áll szemben a 3 egyenlettel szemben. Ezekből nyilvánvaló, hogy ez az egyenletrendszer önmagában nem alkalmas a probléma megoldására. A probléma áthidalására turbulencia modelleket használnak. Számos közelítés létezik. Ezek egyike az örvényviszkozitási modell, amely még Boussinesq (1877) nevéhez fűződik. Ennek lényege, hogy a turbulens feszültséget a nyírófeszültséghez hasonló alakban tételezték fel [10]:
u ′v ′ = −ν T
∂U , ∂y
(6.22)
ahol ν T a turbulens viszkozitási tényező, y a falra merőleges koordináta. Kezdetben ν T -t állandónak tekintették. Ez a közelítés csak nagyon kevés probléma megoldásához használható. Körülbelül 50 évvel Boussinesq elmélete után 1925-ben Ludwig Prandtl csatornaáramlásra publikálta a keveredési úthosszra vonatkozó elméletét. Ezeket a módszereket ma nulla egyenletes módszereknek nevezik. Egy újabb lépcsőt jelentett a turbulencia kutatásában az egy egyenletes turbulencia modellek. Prandtl volt az első, aki 1945-ben a turbulens kinetikus energiára vonatkozó egyenletet használta. Az egyik legnépszerűbb ilyen egy egyenletes modell Spalart és Almaras nevéhez fűződik [11]. Egy másik jellegzetes közelítést jelent a két egyenltes turbulencia modell. Ilyen típusú modellek igen nagy számban állnak rendelkezésre, itt csak néhányat tudunk megemlíteni. Közülük talán a K- ε modellt használják leggyakrabban, ahol a K turbulens kinetikus energia transzportjára vonatkozó egyenlet mellett a ε disszipációjára vonatkozó egyenletet is megoldják. Egy másik gyakran használt modell a K- ω modell. Itt ω úgy fogható fel, mint egy nagy örvényre vonatkozó inverz időlépték. Ezt a módszert nagyon gyakran határréteg-áramlásoknál használják. A tapasztalat szerint a leválásos áramlások többségénél ez a módszer jobb eredményt ad, mint a K- ε . A Reynolds feszültség transzportján alapuló modell bizonyos esetekben (pl. leválásos áramlások, nem körkeresztmetszetű csatornában fellépő szekunder áramlások, görbült felületek menti áramlások, stb.) jobban leírja az áramlási jelenség fizikáját, mint a két egyenletes modellek. Az algebrai Reynolds feszültség modellek az örvény viszkozitási és a Reynolds feszültség transzportján alapuló modellek között jelent egy átmenetet. Ebben az esetben a Boussinesq hipotézist kicserélik a Reynolds feszültség anizotróp tenzora transzportjára vonatkozó algebrai közelítésével.
1.2 NAGY ÖRVÉNYEK SZIMULÁCIÓJA (LES) A nagy örvények szimulációjának népszerű angol elnevezése “Large Eddy Simulation”, ennek kezdőbetűiből ered az LES név, amellyel leggyakrabban erre az eljárásra hivatkoznak. Az LES lényege, hogy - a 6.3 fejezetben röviden ismertetésre kerülő Direkt Numerikus Szimulációval szemben - nem kívánják közvetlenül számítani a turbulens spektrumban lévő nagyon széles teljes skálájú áramlást. A módszer csak a nagy méretű áramlási struktúrákat, nagy örvényeket számítja közvetlenül a Navier-Stokes egyenletekből. Ezek azok az örvények, amelyek alapvetően felelősek az impulzus és hő transzportjáért. A kis örvények főleg csak a turbulens energia disszipációt határozzák meg. Az LES módszer egy térbeli szűrő alkalmazásával kiszűri ezeket az örvényeket a Navier-Stokes egyenletből, és hatásukat egy általános érvényű örvény viszkozitási turbulencia modellel veszi figyelembe, amely nem tartalmaz feladattól függő empirikus állandókat. Ez az LES-nél leggyakrabban használt turbulencia modell Smagorinsky nevéhez fűződik. A módszernél először a diszkretizálást kell elvégezni általában véges differenciák vagy véges térfogatok módszerével. Azután egy térbeli szűrő segítségével kiszűrik a jellemző számítási háló ∆ x rácsméreténél kisebb léptékű szerkezeteket (subgrid-scales) az áramlásból. Hatásukat általában a Smagorinsky modell segítségével veszik figyelembe. A nagy örvényeket a Navier-Stokes egyenletből közvetlenül számítják. A módszer gépidőigényét tekintve a RANS és a DNS (Direkt Numerikus Szimuláció) között van, még mindig igen nagy CPU idő igénnyel. Egyesek az LES-t a “szegény ember DNS”-e-nek titulálják. LES módszer különösen előnyös amikor koherens örvények és struktúrák dinamikáját szeretnénk megérteni. A mai jobb CFD kereskedelmi szoftverek LES számításokra is alkalmasak.
1.3 DIREKT NUMERIKUS SZIMULÁCIÓ (DNS) A Direkt Numerikus Szimuláció (DNS) jelenti jelenleg a legmagasabb számítási szintet. A módszernél a NavierStokes egyenletet az ingadozó sebesség és nyomásértékekre vonatkozóan közvetlen módon oldják meg. Ebben az esetbe egyáltalán nem használnak turbulencia modellt. A feladatot nagyon nagy mértékben megnehezíti az a tény, hogy a nagy méretű térbeli áramlásban az örvényeknek mind a mérete, mind a frekvenciája igen széles határok között mozog. Például egy nagy méretű szélcsatornában azt találták, hogy a jelenlévő örvények mérete néhány tized milliméter és több méter között, frekvenciája pedig 1 és 10000 Hz között változott! Természetesen nagyon finom
térbeli háló és nagyon kis időlépcső szükséges még magas rendű numerikus módszer alkalmazása esetén is, hogy ezt az igen széles térbeli és időbeli spektrumot le tudja írni. Ezért ezt a módszert főleg kis Reynolds számú, egyszerű geometriájú csatornaáramlások vizsgálatára alkalmazzák. Természetesen a számítások igen nagy gépidőt igényelnek, így azok nem képzelhetők el párhuzamos számítások (parallel computing) alkalmazása nélkül, amikor is egyszerre gyakran több ezer összekapcsolt számítógép dolgozik ugyanazon probléma megoldásán. Példáként szeretnénk megemlíteni, hogy ma egy párhuzamos áramlásba helyezett körhenger körüli viszonylag kis Reynolds számú áramlás (Re=4000) alapos vizsgálatához mintegy félmillió óra CPU gépidő szükséges (természetesen ez egyetlen gépre vonatkoztatott CPU idő). A DNS alkalmazásával esetenként olyan áramlásokat is ki tudnak számítani, amelyeket nehéz lenne műszerrel megmérni. A módszerrel olyan nagy pontosság érhető el, hogy egyes eredményei felhasználhatók mérőműszerek kalibrálására is. Ugyanakkor alkalmazásával nagyon hasznos információt nyerhetünk a turbulencia modellek tökéletesítéséhez is. A számítógép folytonos forradalmi fejlődése mellett a módszer szerepe a numerikus áramlástanban egyre nagyobb lesz.
