SZAKDOLGOZAT. Kondor István

SZAKDOLGOZAT

Kondor István 2008

BUDAPESTI MŰSZAKI ÉS GAZDASÁGTUDOMÁNYI EGYETEM ÁRAMLÁSTAN TANSZÉK

CSATORNÁBAN LEJÁTSZÓDÓ HŐÁTADÁS NAGY ÖRVÉNY  SZIMULÁCIÓJA

Kondor István Konzulens: Lohász Máté Márton

2008

2

­­FELADATKIÍRÁS­­

3

Alulírott Kondor István, a Budapesti Műszaki és Gazdaságtudományi Egyetem hallgatója  kijelentem, hogy ezt a szakdolgozatot meg nem engedett segítség nélkül, saját magam  készítettem, és a szakdolgozatban csak a megadott forrásokat használtam fel. Minden olyan  részt, amelyet szó szerint vagy azonos értelemben, de átfogalmazva más forrásból átvettem,  a forrás megjelölésével egyértelműen megjelöltem. ……………………………… Kondor István

4

TARTALOMJEGYZÉK 1 BEVEZETÉS.....................................................................................................................11 1.1 IPARI HÁTTÉR........................................................................................................11 1.2 ELMÉLETI HÁTTÉR...............................................................................................11 1.3 A SZAKDOLGOZAT CÉLJA..................................................................................12 2 IRODALMI ÁTTEKINTÉS..............................................................................................12 2.1 DIREKT NUMERIKUS SZIMULÁCIÓS EREDMÉNYEK...................................12 2.1.1 A vizsgált mennyiségek.....................................................................................12 2.1.2 Irodalmi előzmények.........................................................................................13 2.2 NAGY ÖRVÉNY SZIMULÁCIÓS EREDMÉNYEK..............................................14 3 MATEMATIKAI HÁTTÉR..............................................................................................14 3.1 A SZÁMÍTÁSI TARTOMÁNY...............................................................................14 3.2 AZ ÁRAMLÁST LEÍRÓ EGYENLETEK...............................................................15 3.3 AZ ÁRAMLÁS PERIODIKUSSÁGA......................................................................16 3.4 A HŐMÉSRÉKLET PERIODIKUSSÁGA..............................................................17 3.5 PEREMFELTÉTELEK.............................................................................................19 4 NUMERIKUS SZIMULÁCIÓ..........................................................................................20 4.1 A NAGY ÖRVÉNY SZIMULÁCIÓ ELVE.............................................................20 4.2 HÁLÓ KÉSZÍTÉS.....................................................................................................21 4.3 NAGY ÖRVÉNY SZIMULÁCIÓ FLUENT PROGRAMMAL..............................22 5 KIÉRTÉKELÉS ÉS VALIDÁLÁS...................................................................................24 5.1 A SZIMULÁCIÓ FUTÁSÁNAK ELLENŐRZÉSE.................................................24 5.2 FELTÉTELES ÁTLAGOLÁS..................................................................................30 5.3 SEBESSÉGMEZŐ VALIDÁLÁSA..........................................................................30 5.4 A HŐMÉRSÉKLET-ELOSZLÁS ÖSSZEHASONLÍTÁSA...................................33 5.4.1 Az UHF peremfeltétel........................................................................................33 5.4.2 A CTD peremfeltétel..........................................................................................36 5.4.3 Korrelációk........................................................................................................37 5.5 KOHERENS STRUKTÚRA.....................................................................................41 6 ÖSSZEFOGLALÁS..........................................................................................................44

5

KÖSZÖNETNYILVÁNÍTÁS Ezúton szeretnék köszönetet mondani konzulensemnek, Lohász Máté Mártonnak,  hathatós segítségéért és hogy mindig rendelkezésemre állt, ha elakadtam. Külön szeretném megköszönni Illyés László Zoltánnak és Mach Évának szíves  segítségüket és végtelen türelmüket.

6

JELÖLÉSJEGYZÉK A

csatorna keresztmetszete

C

Courant­szám

CS

Smagorinsky konstans

cp G k L m˙ Nu p P Pr q qw Q Reb  Reτ Sij t tavg T

állandó nyomáson vett fajhő térbeli átlagolás magfüggvénye turbulens kinetikus energia periódus hossza tömegáram Nusselt­szám nyomás csökkentett nyomás Prandtl szám fajlagos bevezetett hőmennyiség fali hőáram sebességgradiens tenzor második invariánsa átlagos Reynolds­szám turbulens Reynolds­szám alakváltozási tenzor idő átlagolási idő súrlódási hőmérséklet

T 〈 T m 〉 ui, u, v, w u ub

csökkentett hőmérséklet kevert átlaghőmérséklet sebességkomponensek súrlódási sebesség átlagsebesség

xi, x, y, z

koordináta irányok

Görög jelölések:     ij  t T

hőátadási tényező nyomás gradiens hőmérséklet gradiens a csatorna magasságának fele Kronecker delta szűrő szélesség időlépés nagysága hőmérséklet­különbség 7



kinematikai viszkozitás

t

örvényviszkozitás



sűrűség

ij

szűrőméret alatti feszültség tenzor

w

fali csúsztató feszültség



transzformált hőmérséklet

ij

örvénytenzor

Indexek: ( )'

ingadozó komponens

( )+

u­val, T­val és ­vel dimenziótlanított

( )*

­val és T­vel deimenziótlanított

( )rms

root mean square

( )max

pontbeli maximum

〈〉

csatorna keresztmetszetében átlagolt

 ( )_mean statisztikailag átlagolt

8

TARTALMI ÖSSZEFOGLALÓ A szakdolgozatomban ismertetem a direkt numerikus szimuláció és a nagy örvény  szimuláció alapvető elméleti hátterét, a csatornaáramlásban való alkalmazásuk fontosabb  irodalmi eredményeit, valamint bemutatom a periodikus  hőátadás  számításának módját.  Ezt  követően  nagy örvény  szimuláció  segítségével  vizsgálom  végtelen  síklapok  közötti  kialakult   turbulens   áramlásban   a   hőtranszportot   különböző   termikus   peremfeltételek  esetén. A sebesség­ és hőmérsékletmező főbb mennyiségeinek pillanatnyi értékeit valamint  időbeli és térbeli átlagait is közlöm. A   számítástechnika   rohamos   fejlődésének   köszönhetően   lehetővé   vált   a   Navier­ Stokes   egyenlet   numerikus   megoldása,   s   elérhetővé   vált   az   áramlások   időfüggő  modellezése. Napjainkban egyre többen végeznek nagy örvény szimulációt, így a módszer  fejlesztése   és   ellenőrzése   különösen   fontos.   A   dolgozatban   vizsgált   áramlás   turbulens  Reynolds­száma   180,   míg   Prandtl­száma   0,71.   A   számítás   során   kétféle   peremfeltételt  alkalmazok a hőmérsékletre vonatkozóan. Az egyik esetben a fali hőáram állandó, míg a  másikban a két fal hőmérsékletkülönbségét tartom konstans értékén. A kiértékelés során feltételes átlagolást végzek, az alkalmazott feltétel a Q­kritérium.  Ellenőrzöm a sebesség­ és hőmérséklet­profilokat, valamint ezek ingadozásait is. Továbbá  vizsgálom   ezek   korrelációit,   az   eredményeket   pedig   DNS   adatokhoz   viszonyítom.   A  koherens struktúrákat a Q­kritérium alapján tanulmányozom.

9

ABSTRACT In my thesis I review the basic theoretical background of direct numerical simulation  and large eddy simulation, the most important results of their usage and I also present a  method to compute periodic heat transfer. After that I study the turbulent heat transfer in a  fully developed channel flow. I analyzed two different thermal boundary conditions with  the help of large eddy simulation. I report the main characteristics of the velocity­ and  thermal field including instantaneous values and both spatial and temporal averages. Thanks   to   the   rapid   development   of   computer   science   it   became   possible   to  numerically solve the Navier­Stokes equation. Consequently, time­dependent modeling of  flows are now available. Nowadays large eddy simulations are performed more frequently,  therefore its development and validation became crucial. In the present study the turbulent  Reynolds   number   of   the   flow   was   180,   while   the   applied   Prandtl   number   was   0.71.   I  examine the effect of two thermal boundary conditions. The first case is the uniform heat  flux (UHF) over the surfaces, while the second one is the constant temperature difference  (CTD) between the top and bottom walls. I   adopt   a   conditional   averaging   method   in   the   evaluation   process,   in   which   the  applied condition is the Q­criteria. I verify the profiles of the velocity and thermal fields  and also their fluctuations. In addition, I study their correlations and compare the results to  DNS data. I examine the coherent structures from the viewpoint of Q­criteria.

10

1 BEVEZETÉS 1.1 IPARI HÁTTÉR A   mérnöki   gyakorlatban   egyre   fontosabbá   válik,   hogy   a   gyártást,   felépítést  megelőzően már pontos ismeretekkel rendelkezzünk a termék leendő tulajdonságairól. Ez  különösen   lényeges   az   áramlástani   és   hőtani   problémák   esetében,   példának   okáért  gondoljunk egy turbina lapátkerékre, amelyet 1100­1600°C körüli hőmérsékletű gáz forgat.  Ez   magasabb   érték,   mint   a   lapátok   anyagának   olvadáspontja,   ezért   a   megfelelő  üzembiztonság   eléréséhez   elengedhetetlen   megfelelő   hűtésről   gondoskodni.   A   belső  hűtőrendszer   tervezésekor   jól   kihasználható   a   turbulens   áramlások   jobb   hővezetése   a  laminárisokhoz   képest,   ám   ennek   elengedhetetlen   feltétele,   hogy   a   turbulenciát   kellő  pontossággal   modellezzük.   Fontos   tudni   például,   hogy   a   forgás   hogyan   hat   a   hűtő  csatornában kialakuló áramlásra és ezek hatására milyen hőmérséklet­eloszlás jön létre.  Hogy   idáig   eljussunk   először   az   alapjelenségeket   kell   megértenünk,   majd   a   lényeges  tulajdonságaikat a modellbe be kell építenünk. Ennek szellemében a szakdolgozat témája a  turbulens hőátadás számítása kialakult csatornaáramlásban.

1.2 ELMÉLETI HÁTTÉR Az   áramlástan   alapegyenletei   régóta   ismertek,   ám   megoldásuk   néhány   egyszerű  esetet leszámítva gyakorlatilag lehetetlen volt. A helyzet az 1980­as években változott meg,  a számítógépek rohamos fejlődésének köszönhetően. Ekkor vált lehetővé az áramlásokat  leíró Navier­Stokes egyenlet (NS) numerikus megoldása. A csatornában kialakult áramlás jól dokumentált, bőséges irodalommal rendelkezik.  Ezek nagy része a Reynolds átlagolt Navier­Stokes (RANS) egyenleten alapuló turbulencia  modelleket használja. Gyorsaságuk mellett legnagyobb hátrányuk, hogy a Navier­Stokes  egyenletet időben átlagolják, így időfüggő szimulációra csak korlátozottan alkalmazhatók. A   direkt   numerikus   szimuláció   (DNS)   során   a   Navier­Stokes,   a   kontinuitási   és   a  transzport egyenleteket is közvetlenül, numerikusan megoldják. Ez a módszer semmilyen  közelítést  nem   alkalmaz,   ugyanakkor  a  pontos  eredmény  feltétele  a   megfelelően   finom  háló, ami a turbulencia legkisebb léptékének felbontására is alkalmas. Következésképpen  számításigénye hatalmas, a mindennapi mérnöki feladatok megoldása során használata a  jelenlegi számítógépekkel nem gazdaságos. A fenti probléma megoldására dolgozták ki a nagy örvény szimulációt (LES). Mivel a  kialakult   áramlás   tulajdonságait   elsősorban   a   nagyobb   örvények   határozzák   meg,   ezért  elegendő olyan háló készítése, mely ezeket felbontja, az ennél kisebb léptékű örvényeket  pedig valamilyen turbulencia modellel közelíti, ezeket Subgrid Scale (SGS) modelleknek  nevezzük.

11

1.3 A SZAKDOLGOZAT CÉLJA A   szakdolgozat   célja   a   csatornában   kialakult   turbulens   áramlás   áramképének  meghatározása,   valamint   a   hőmérséklet­eloszlás   vizsgálata   különböző   peremfeltételek  mellett,   a   kapott   eredmények   összehasonlítása   az   irodalomban   közöltekkel.   A   feladat  elvégzése során ismertetem a periodikus hőátadás számításának módját, a numerikus háló  készítésének   lépéseit   Gambit   szoftverben,   továbbá   a   nagy   örvény   szimuláció  alapkoncepcióját és alkalmazásának mikéntjét a Fluent rendszerben. A dolgozatban bemutatom a témában megjelent fontosabb irodalmi eredményeket,  kiemelve   a   direkt   numerikus   szimulációval   és   a   nagy   örvény   szimulációval   végzett  vizsgálatokat.   Ezt   egy   áttekintés   követi,   melyben   a   számítások   matematikai   hátterét  közlöm. A szimuláció menetének részletes leírására a negyedik fejezetben kerül sor. Ezután  elvégzem   az   eredmények   kiértékelését   és   validálását,   utóbbi   során   direkt   numerikus  szimulációs adatokat tekintek referenciának.

2 IRODALMI ÁTTEKINTÉS 2.1 DIREKT NUMERIKUS SZIMULÁCIÓS EREDMÉNYEK 2.1.1 A vizsgált mennyiségek A csatornaáramlás mindenkori dimenziós jellemzői függnek a csatorna geometriai  alakjától, az áramló közeg anyagjellemzőitől és az áramlási sebességtől, következésképpen  ezekből   az   adatokból   nem   vonhatunk   le   általános   következtetéseket.   Éppen   ezért   a  probléma vizsgálatakor szükségünk van olyan dimenziótlan mennyiségekre, melyekkel a  különböző egyedi áramlások összehasonlíthatóak. Az áramlási sebesség jellemzésére a Reynolds­számot használjuk, kétféle definíció  szerint. Az átlagos Reynolds­számot a (2.1) egyenlet írja le: Reb=

u⋅ 

(2.1)

ahol u az áramlás (x) irányú (lásd 1. ábra) átlagsebesség a keresztmetszetben,   a csatorna  magasságának   fele,    pedig   a   közeg   kinematikai   viszkozitása.   A   (2.2)   összefüggésből  látható, hogy hasonlóan értelmezhető a turbulens Reynolds­szám (Reτ) is: Re=

u⋅ 

(2.2)

ahol az u súrlódási sebességet a  fali csúsztató feszültség segítségével a (2.3) egyenlet adja  meg:

12

u =



w 

(2.3)

A hőmérséklet transzport jellemzésére a Prandtl­számot használjuk a (2.4) definíció  szerint: Pr=

 

(2.4)

ahol  a hőátadási tényező. A   hőmérsékletek   dimenziótlanítására   a   súrlódási   hőmérsékletet   (T)   alkalmazzuk  (2.5) szerint: T =

qw ⋅c p⋅u

(2.5)

Az eredményeket általában a dimenziótlanított fali koordináta függvényében adják  meg, melyet definíció szerint a (2.6) ír le: 

y =

u y 

(2.6)

A Reynolds­szám a legkisebb és legnagyobb örvények léptékének arányát is jellemzi,  tehát   ennek   növekedésével   a   hálót   finomítani   kell.   A   csatornában   kialakuló   áramlás  biztosan turbulens, ha az átlagos Reynolds­száma nagyobb, mint 1800. A kutatások során  általában   ennél   magasabb   értékekkel   dolgoznak,   jellemzően   2800­3300   között.   Ennek  numerikus számítására már az 1990­es évek elején lehetőség nyílt, nem meglepő hát, hogy  erre az esetre vonatkózóan áll a rendelkezésre a legtöbb adat. A szerzők praktikus okok  miatt általában a turbulens Reynolds­számot adják meg, ez a fentiek szerint általában 180. Látható   tehát,   hogy   a   Reynolds­szám   növelésével   drasztikusan   növekedik   a   háló  szükséges cellaszáma és ezzel együtt a számításigény is. Még nehezebb a helyzet, ha a  Prandtl­számot   növeljük,   mivel   a   hőmérséklet­fluktuációk   legkisebb   mérettartománya  fordítottan   arányos  Pr1/2­vel  (bővebben   lásd   [1]   és   [2]).   Ebből   következik,   hogy   a  hőmérséklet­ és a sebességmező felbontásához szükséges cellaszám hányadosa arányos a  Prandtl­szám 3/2­ik hatványával (bővebben lásd [3]).

2.1.2 Irodalmi előzmények Napjainkban   a   turbulens   hőátadásra   egyre   többen   végeznek   direkt   numerikus  szimulációkat különböző geometriai elrendezésekre. A kialakult turbulens csatornaáramlás  a   leggyakrabban   vizsgált   eset,   köszönhetően   egyszerű   alakjának   és   annak,   hogy  segítségével mégérthetjük a szilárd fal és folyadék közötti hőátadás alapvető természetét. A témában az első kutatás Kim & Moin [4] nevéhez fűződik, három Prandtl­számot  vizsgáltak: Pr=0,1; 0,71 és 2; Re=180 mellett még 1989­ben. Állandó térfogati hőforrást  és konstans fali hőmérsékletet tételeztek fel. Az átlaghőmérséklet, a hőmérséklet változás,  valamint a turbulens hőáram profilokat közölték, de a magasabb rendű momentumokra  13

vonatkozó   mérkegegyenletek   tagjait   nem   részletezték.   Később   valamivel   alacsonyabb  Reynolds­szám   mellett,  Re=150,   Kasagi   et   al.   [5]   és   Kasagi   &   Ohtsubo   [6]   végeztek  számításokat  Pr=0,025  és  0,71­re. Ők prezentálták a   hőmérséklet változás, a turbulens  hőáram és disszipációik mérlegegyenleteinek tagjait is. A falak mentén az átlagos hőáramot  állandónak   tekintették,   de   a   pillanatnyi   értéke   változhat   a   helytől   és   időtől   függően.  Kawamura   et   al.   [7]   vizsgálták   a   Prandtl­szám   hatását,  Pr=0,025  és  5  között   nyolc  lépésben  Re=180  mellett.   Mivel   az   eddig   felsorolt   eredményeket   alacsony   Reynolds­ számon   kapták,   érvényességük   általánosságát   óvatosan   kell   kezelni.   Később   [8]  tovább  emelték   a   Reynolds­számot,   egészen   395­ig.   Mindkét   esetben   állandó   fali   hőáramot  alkalmaztak   fali   peremfeltételnek.   Wikström   &   Johansson   a   falak   közötti   állandó  hőmérsékletkülönbség feltételt vizsgálták  Re=265  esetén. A peremfeltételek hatását Seki  et al. [9] részletesen ellenőrizték  Re=180  és  395  mellett, miközben a Prandtl­szám 0.71  volt. Napjainkban Fukugata et al. [10] folytattak kutatásokat a bőrsúrlódás csökkentése és a  hővezetés   javítása   érdekében,   ennek   során   az   állandó   hőáram   és   állandó  hőmérsékletkülönbség feltételt is alkalmazták.

2.2 NAGY ÖRVÉNY SZIMULÁCIÓS EREDMÉNYEK A nagy örvény szimuláció egy erős eszköz a mérnökök kezében az időfüggő áramlási  problémák megoldásához. Az utóbbi néhány évben gyors fejlődésen ment keresztül és az  iparban   is   egyre   szélesebb   körben   alkalmazzák.   Mivel   a   kicsi   örvényeket   közelítő  eljárásokkal modellezi, ezért fontos kérdés, hogy mennyire pontos a számítás és az adott  problémánál melyik SGS modell használata javasolt. A módszer részletes bemutatására a  negyedik fejezetben kerül sor. Chatelain   et   al.   [11]   a   hőmérséklet­eloszlás   vizsgálatakor   különböző   numerikus  módszereket hasonlított össze. Munkájukban rámutatnak, hogy jobb eredmény érhető el, ha  a hőmérséklet transzport egyenletére regularizáló, a sebességére pedig középpontos sémát  alkalmazunk,   mintha   mindkettőre   az   utóbbit   használnánk.   Montreuil   et   al.   [12]   öt  különböző SGS modellt vizsgált a skalár transzport leírására, Fick­féle és nem­Fick­féle  módszereket egyaránt.

3 MATEMATIKAI HÁTTÉR 3.1 A SZÁMÍTÁSI TARTOMÁNY A   dolgozatban   két   párhuzamos,   végtelen   kiterjedésű,  2  távolságú   síklap   közötti  kialakult áramlást tanulmányozok. A tényleges számítási tartomány ennek egy 3 széles  és 5 hosszú, téglatest alakú része, amely az (1. ábrán) látható. A későbbiekben használt  koordinátarendszer a következőképp definiálható:  x  az áramlás iránya,  y  a fali normális  14

irány, míg z az előző kettőre merőleges irányba mutat.

(1. ábra: Vizsgált geometria, forrás: [3])

A kiválasztott térről feltételezzük, hogy a belépéstől elegendően messze van, tehát az  az   áramképet   már   nem   befolyásolja,   az   áramlás   kialakultnak   tekinthető.   A   vizsgált  geometria oldalarányai úgy lettek meghatározva, hogy minden irányban azonos cellaszámot  véve   a   felbontott   örvénnyesség   anizotrópiájához   illeszkednek,   ezért   a   numerikus  szimuláció   során   a   cellák   hossz­   és   keresztirányú   méretének   (x+  és  z+)   vizsgálatát  mellőzöm, csak az y+ értékének megfelelő beállítását ellenőrzöm.

3.2 AZ ÁRAMLÁST LEÍRÓ EGYENLETEK Az   áramlástanban   a   folyadékok   mozgását   irányító   összefüggés   a   Navier­Stokes  egyenlet (3.1): ∂ui ∂u ∂2 u i 1 ∂p u j i =−   ∂ xi ∂t ∂xj ∂ xj∂ x j

(3.1)

és a kontinuitási egyenlet (3.2): ∂ui =0 ∂ xi

(3.2)

ahol xi a koordináta­irányokat, ui a a folyadékrész adott irányba eső sebességét, p a nyomást  és  t  az   időt   jelöli.   Bár   a   mérnöki   gyakorlatban   nagy   sikereket   ért   el   a   Navier­Stokes  egyenlet,   mégis   a   hét   legfontosabb   matematikai   probléma   között   szerepel,   mivel   a  megoldások  egzisztenciája   és   unicitása,   ami   a   lehetséges   megoldások   keresésének  alapfeltétele,   a   mai   napig   nem   bizonyított.   Ennek   ellenére   a   numerikus   megoldások  pontosan egyeznek a mérésekkel, így használata kézenfekvő. A (3.1) egyenlet implicit módon tartalmaz egy elhanyagolást, nevezetesen a sűrűséget  15

állandónak feltételezi, így a természetes konvekciót figyelmen kívül hagyjuk.

A koordinátákat ­val, a változókat ­vel, u­val és T­val, a súrlódási hőmérséklettel  dimenziótlanítjuk. Így a Navier­Stokes egyenlet a (3.3), a kontinuitási egyenlet pedig   a  (3.4) alakban írható: 2  ∂ p 1 ∂ ui ∂ p u =− ∗  ⋅  ∗⋅ i1 ∂ t∗ ∂ x ∗j ∂ x i Re ∂ x∗2 ∂ x1 j

∂ui

 j

∂ ui

∂ui ∂ x ∗i

=0

(3.3) (3.4)

ahol a (3.3) jobb oldalának utolsó tagja az áramlás irányú átlagos nyomás gradiens.

3.3 AZ ÁRAMLÁS PERIODIKUSSÁGA Az áramlást a kialakultnak tekinthető tartományban periodikusnak tekintjük, amihez  a   szükséges   feltételeket   először   Patankar   et   al.   [13]   vezette   le 1977­ben. Nem szabad azonban megfeledkezni arról, hogy a turbulencia, mint hely­ és  időfüggő jelenség sosem lehet periodikus. Ez a tulajdonsága, bár még nem nyert explicit  bizonyítást [14], feltehetően kaotikus viselkedéséből adódik. Mivel azonban a turbulencia  jellemzőinek statisztikai átlagát vizsgáljuk, így a periodizálás nem okoz súlyos hibát. Tekintsük két kellően nagy síklap közötti áramlásnak egy olyan részét, mely elég  messze van a belépéstől ahhoz, hogy a belépési problémát elhanyagolhassuk. Az így kapott  térben az áramlás kialakultnak tekinthető. Vegyük észre, hogy ez a definíció ekvivalens  azzal   amikor   végtelennek   tekintjük   a   csatorna   oldalait.   A   koordinátarendszert   a   3.1  fejezetben tárgyaltaknak megfelelően vegyük fel. Ekkor nyilvánvaló, hogy az áramlásban a  sebesség   nagysága   független   az  x  koordinátától,   amit   dimenziótlan   alakban   a   (3.5)  összefüggés fejez ki. ∂u =0 ∂ x∗

(3.5)

A (3.5) feltétel szerint az  x  irányú periodikusság végtelen rövid lehet. Az áramlás  modellezésekor azonban figyelmbe kell venni a turbulencia hosszléptékét is, a vizsgálati  tartományt úgy kell felvenni, hogy az minden irányban nagyobb legyen, mint a legnagyobb  örvények adott irányú kiterjedése. A 3.1 fejezetben ennek szellemében jártam el. Lamináris   esetben   ezt   a   v=0   feltétel   egészíti   ki,   de   turbulens   áramlásnál   ez   a  pillanatnyi   sebességekre   nem   igaz,   csak   az   időátlagokra,   ezért   más   peremfeltételt   kell  keresni. Most a vizsgált tartományt osszuk fel az áramlás irányában n darab L hosszúságú  részre.   A   (3.5)   egyenlettel   összhangban   felírható   (3.6),   mely   a   sebességmező  16

periodikusságát fejezi ki és turbulens esetben is helyes. u x , y , z =u xL , y , z =u  x2L , y , z =...=u  xn⋅L , y , z v  x , y , z =v xL , y , z =v x2L , y , z =...=v  xn⋅L , y , z  w  x , y , z=v  xL , y , z=w x 2L , y , z =...=w x n⋅L , y , z 

(3.6)

Habár  u,  v  és  w  nagysága   folyamatosan   változik   az   áramlási   térben,   a   (3.6)  összefüggés   lehetővé   teszi,   hogy   a   sebességkomponenseket   csak   egy  L  hosszúságú  térrészben vizsgáljuk. Továbbá a kialakult áramlás jellemzőit anélkül határozhatjuk meg,  hogy a belépési problémával foglalkoznánk. Ahhoz,   hogy   a   csatornában  x  irányú   áramlás   jöjjön   létre,   abban   az   irányban  nyomásesésre van szükség, ebből következik, hogy nem vonatkozhat rá a sebességeket leíró  periodikusság. Ha a nyomást felrajzoljuk az y koordináta függvényében egy tetszőleges  x  és  x+L  helyen vett keresztmetszetben, azt tapasztaljuk, hogy a két görbe alakja hasonló,  csak el vannak tolva egy konstanssal. Ha az x+2L­hez tartozó görbét is berajzoljuk, látható  lesz, hogy az is ugyanannyival van eltolva lefelé, mint az  x+L  helyhez tartozó az  x­ben  levőhöz képest. Következésképpen a (3.7) egyenlet minden x­re teljesül: p x , y , z − p  xL , y , z = p  xL , y , z− p  x2L , y , z =...

(3.7)

Ez az állandó nyomásesés hozza létre az átlagos tömegáramot a csatornában. A (3.7)  összefüggés alapján definiálható  konstans (3.8) szerint: =

p  x , y , z − p  xL , y , z L

(3.8)

Ezután   Patankar   et   al.   [13]   javaslata   alapján   érdemes   a   nyomást   (3.9)   szerint  szétválasztanai: p x , y , z =P x , y , z − x

(3.9)

Belátható,   hogy   a  x  okozza   a   globális   tömegáramot,  P(x,y,z)  pedig   a   helyi  mozgásokért felelős. Utóbbira a sebességeket leíró periodikusság is jellemző, vagyis (3.10)  felírható: P x , y , z =P xL , y , z =P  x2L , y , z=...

(3.10)

3.4 A HŐMÉSRÉKLET PERIODIKUSSÁGA A hőmérsékletmező leírásakor Patankar et al. [13] és Seki et al. [9] gondolatmenetét  követtem. Hasonlóan ahhoz ahogy azt már a nyomás vizsgálatakor láttuk, ha a hőmérsékletet az  y koordináta függvényében felrajzoljuk különböző x, x+L és x+2L helyeken, kiderül, hogy  a kapott görbék hasonlóak, csak egy konstans értékkel el vannak tolva egymáshoz képest.  Ebből adódik, hogy a (3.7)­hez hasonlóan a hőmérsékletre felírható a (3.11) egyenlet: T  xL , y , z −T  x , y , z=T  x2L , y , z−T  xL , y , z =...

(3.11)

17

Ezt követően definiáljuk ­t (3.12) szerint: =

T  x L , y , z −T  x , y , z  L

(3.12)

Könnyen belátható, hogy  (3.13) szerint is felírható: =

q m⋅c ˙ p⋅L

(3.13)

ahol q a fajlagos bevezetett hőmennyiség,  m˙  a tömegáram, cp pedig a folyadék fajhője.  Ezután a hőmérséklet a nyomáshoz hasonlóan szétválasztható, melyet a (3.14) összefüggés  mutat be: T  x , y , z = xT  x , y , z 

(3.14)

T ­ról belátható, hogy (3.15) szerint periodikus: T  x , y , z =T  xL , y , z=T  x2L , y , z=...

(3.15)

A hőmérsékletre vonatkozóan kétféle peremfeltételt vizsgáltam. Az első esetben a fali  hőáram nagysága volt állandó (Uniform Heat Flux – UHF), a másodikban a falak állandó,  de   eltérő   hőmérsékletűek   (Constant   Temperature   Difference   –   CTD)   voltak.   A   két  különböző feltételnél a számítás menete eltérő. Állandó  fali  hőáram  esetén  a    transzformált  hőmérsékletet  célszerű  a   (3.16)­nak  megfelelően bevezetni:  

=T w −T

(3.16)

ahol Tw a fali hőmérséklet. Ekkor a (3.14) összefüggés dimenziótlanítva a (3.17) alakot ölti: 

∗

∗

∗

T  x , y , z =

d 〈 T m 〉 dx

∗

x ∗−  x∗ , y ∗ , z ∗

(3.17)

ahol  〈 T m 〉  az ún. kevert átlaghőmérséklet, mely (3.18) szerint fejezhető ki: 1

〈 T

 m

∫0 u1 T  dy∗ 〉= 1 ∫0 u1 dy ∗

(3.18)

A jelen elrendezésben az áramlás irányú gradiense (3.19) egyenletté egyszerűsödik: d 〈 T m 〉 dx

∗

=

1 〈 u1 〉

(3.19)

ahol   〈 u1 〉   az átlagsebesség a keresztmetszetben. A fenti átalakítások után az energia  egyenlet a (3.20) formában írható fel:   2  u1 ∂ 1 ∂   ∂ u j = ⋅ ∗2   ∗ ∗ ∂t ∂ x j Re⋅Pr ∂ x j 〈 u 1 〉

(3.20)

18

Az állandó hőmérséklet­különbség esetén a transzformált hőmérskletet (3.21) szerint  célszerű bevezetni: =

T −T 1 T 2−T 1

(3.21)

A pillanatnyi hőmérséklet ebben az esetben is szétválasztható, ekkor a (3.14) egyenlet  dimenziótlan formában a (3.22) alakot ölti:  T ∗  ∗ ∗ ∗ T  x , y , z = y   x , y , z  2 

∗

∗

∗

(3.22)

ahol T a két fal közötti hőmérsklet­különbség. Ezt az energia egyenletbe beírva arra a (3.23) összefüggés adódik: ∂  ∂ 1 ∂2   T   u j = ⋅ − u2 2 ∂ t∗ ∂ x ∗j Re⋅Pr ∂ x ∗2 j

(3.23)

3.5 PEREMFELTÉTELEK A   fenti   egyenletek   megoldásához   szükség   van   megfelelő   számú   peremfeltétel  megadására. A sebességekre vonatkozóan a falak felületén a szokásos tapadási törvényt  alkalmazzuk, amit a (3.24) egyenlettel fejezhetünk ki: ui=0  az y=0 és y=2 helyeken

(3.24)

Amint   azt   az   előző   fejezetben   megjegyeztük   az   egyik   esetben   a   falakon   állandó  hőáramot  feltételezünk.  Ez  azt  jelenti,  hogy  a  felületen  időegység   alatt   átmenő  átlagos  hőmennyiség   konstans,   de   a   pillanatnyi   értéke   változhat   térben   és   időben.   Ennek  megfelelően a hőmérsékletre vonatkozó peremfeltétel (3.25) szerint írható: =0  az y=0 és y=2 helyeken

(3.25)

Ezt a premfeltételt használják a DNS kutatások jelentős részénél. Ebben az esetben  feltételezik, hogy a fal végtelenül jó hővezető és ezért a hőmérséklete is mindig állandó. Én  ettől egy kicsit eltérő esetet vizsgáltam, mert csak a fali hőáramot írtam elő, de az előbbi  egyszerűsítéssel   nem   éltem.   A   falakon   használt   másodfajú   peremfeltételeket   a   (3.26)  egyenlet írja le: ∂ =1  az y=0 helyen, ∂ x2 ∂ =−1  az y=2 helyen ∂ x2

(3.26)

19

A második esetben a falakon állandó a hőmérséklet, itt a hőáram változhat. Ekkor a  peremfeltételek a (3.27)­nak megfelelően alakulnak: =−1  az y=0 helyen,

=1  az y=2 helyen

(3.27)

A   feladat   így   már   lezárt   és   megoldható.   A   következő   fejezetben   a   szoftveres  megoldás lépéseit mutatom be, majd a kapott eredményeket értékelem.

4 NUMERIKUS SZIMULÁCIÓ 4.1 A NAGY ÖRVÉNY SZIMULÁCIÓ ELVE Ahogy   a   bevezetésben   már   megemlítettem,   direkt   numerikus   szimuláció   során   az  áramlás összes mérettartományát felbontják, így ez az eljárás adja az elméletileg elérhető  legpontosabb eredményt. Ugyanakkor éppen a teljes felbontás adja a módszer legnagyobb  gyakorlati   hátrányát   is,   mivel   óriási   számítási   kapacitást   igényel.   Ebből   kifolyólag   az  elérhető számítógépes erőforrások a viszonylag alacsony Reynolds­számú tartományokra  korlátozzák a kutatásokat. Mint a neve is sugallja a nagy örvény szimuláció az áramlásnak  csak a nagyobb léptékű részeit próbálja meg leírni. A DNS­hez hasonlóan itt is a Navier­ Stokes   egynletet   oldják   meg,   de   a   háló   térbeli   felbontása   nem   elégséges   a   legkisebb  örvények modellezésére. Ebből adódóan a LES a valóságos áramlás egy olyan közelítését  adja, melyből egy bizonyos méret alatti rész hiányzik. Ugyanakkor dinamikai szempontból  a kis mérettartományok is fontosak lehetnek, ezért az eredmény korrekcióra szorul. Ezt úgy  oldják meg,  hogy egy új tagot adnak  a mozgásegyenlethez,  amit  SGS  (Subgrid  Scale)  modellnek nevezünk. Legfőbb tulajdonsága, hogy csak a LES által felbontott legkisebb  léptékeket módosítja, pontosítja ezek leírását, miközben a nagyobbakat nem befolyásolja.  A   nehézséget   az   adott   geometriához   és   numerikus   módszerhez   megfelelő   SGS   típus  kiválasztása   jelenti.   Ennek   bonyolultsága   és   terjedelme   meghaladja   jelen   dolgozat  célkitűzéseit. Ahhoz,   hogy   csak   a   nagy   örvényeket   vegyük   figyelembe   a   teljes   sebességmezőt  szűrnünk kell. Így a Navier­Stokes egyenletet csak a felbontott mérettartományra oldjuk  meg, a kisebbeket az SGS modell segítségével közelítjük. A LES­ben általánosan használt  szűrő operátort a (4.1) egyenlet mutatja be:   x , t=∫V G  y−x , x ⋅  y , t d y

(4.1)

ahol  G  a magfüggvény, melynek (4.2)­nek megfelelően kompakt tartójúnak kell lennie,  valamint teljesülnie kell (4.2)­nek:

∫V G r , x d r =1

(4.2)

Az előző összefüggésekben  a szűrő szélességére utal. Lényegében a ­nál nagyobb  örvények nagynak tekintjük és kiszámítjuk, az ennél kisebbeket pedig csak modellezzük.  Ha a  G  függvény homogén és izotróp, akkor a szűrt Navier­Stokes egyenlet (4.3) szerint  20

írható fel: ∂ ui ∂ ui uj ∂2 ui ∂ij 1 ∂ p  =− − 2 −  ∂ xi ∂t ∂xj ∂x j ∂xj

(4.3)

A kontinuitási egyenlet lineáris, ezért az továbbra is a (3.5) alakban írható. A (4.3)  összefüggésben szereplő  ij­t tradícionálisan hálóméret alatti (SGS) feszültségnek nevezik.  Valójában szűrőméret alatti feszültség lenne a helyes megnevezése, mivel ­tól függ, ami  viszont, mint később látni fogjuk, elméletileg független a hálótól. Az SGS feszültség (4.4)  egyenlettel fejezhető ki: ij =ui u j−ui uj

(4.4)

A (4.5) összefüggés jobb oldalának első tagját nem lehet kiszámolni, ezért ij­re az  örvényviszkozitás modellt használjuk, így a (4.5) kifejezés adódik: 1 ij −  kk ij =−2 t Sij 3

(4.5)

ahol   az   örvényviszkozitás,   ezt   kell   valahogyan   közelíteni.   Az   első   és   legelterjedtebb  megoldást Smagorinsky dolgozta ki, melyet a (4.6) összefüggés mutat be: 2 t = C S   ∣S∣

(4.6)

ahol  S  az alakváltozási tenzorból számolható és a (4.7) egyenlettel definiálhatjuk: ∣S∣= Sij Sij

(4.7)

  CS  pedig   a   Smagorinsky   konstans,   mely   különböző   módszerekkel   határozható   meg,  például mérési vagy DNS eredményekből.

4.2 HÁLÓ KÉSZÍTÉS Numerikus   szimuláció   során   a   problémát   modellező   differenciálegyenleteket   nem  analitikusan   oldjuk   meg,   hanem   azokat   valamilyen   módszer   segítségével  differenciaegyenletekkel  helyettesítjük.   Mivel   ezek   az   egyenletek   a   tér   adott  pontjaiban  érvényesek, ezért a megoldás sem lesz folytonos. Az eredmény pontossága nagyban függ  attól, hogy a vizsgált tartományt hogyan diszkretizáljuk. Mivel a csatorna geometriája egyszerű, ezért a háló elkészítése nem jelent akkora  problémát, mint egy átlagos mérnöki feladatnál. Ettől függetlenül persze néhány dologra  különös figyelmet kell fordítani. A hálót a Gambit szoftver segítségével készítettem. Az alapelgondolás az volt, hogy  először   40x40x40   cellára   osztom   a   csatornát,   ami   már   elég   ahhoz,   hogy   a   kiszámított  áramlás turbulens legyen. Ennek segítségével meg tudtam határozni, hogy mekkorának kell  lennie a fal melletti első cellának és az áramlás főbb jellemzőiről is nagyságrendi képet  kaptam. Ezt követően létrehoztam egy 72x72x72 elemet tartalmazó hálót, amit később a  végső számításnál használtam. 21

A háló néhány egyszerű lépésben előállítható: ­   A   csatorna   magasságának   felét   ()   egységnyinek   választva   létrehoztam   egy  téglatestet az (1. ábrának) megfelelően. ­ Tekintve a 3. fejezetben leírtakat, az áramlás periodikus, mind x és z irányban, ezért  szükséges   a   megfelelő   oldalak   összekapcsolása.   Két   párt   hoztam   létre,   az   egyik   az  áramlásra   merőleges   két   lapból   áll,   a   másikat   az   áramvonalakkal   párhuzamos   kisebb  téglalapok alkotják. Erre azért van szükség, mert ahhoz, hogy a periodikus peremfeltételt  pontosan definiálhassuk garantálni kell az egyforma hálót a felületpárokon. ­ Ezek után az élek hálózása következett. Mivel a csatorna méreteit úgy határoztuk  meg a 3.1 fejezetben, hogy az már “tartalmazza” a turbulencia anizotrópiáját, ezért x és z  irányban egyforma hosszúságú részekre osztottam az oldalakat. Y irányban sűríteni kell a  hálót, hogy a fali határréteg is megfelelően fel legyen bontva. Itt a “Successive Ratio”  opciót használtam kétoldali sűrítéssel, melynek mértéke mindkét oldalnál 1,05 volt. Így az  első cella mérete 0,010434452 lett. ­   A   három   élen   lévő   háló   egyértelműen   meghatározza   a   testhálót,   így   ennek  generálása már nem jelent problémát. ­ A lapokra vonatkozó peremfeltételeket is itt kell megadni, beállítottam az előbb  említett   két   periodikusságot,   valamint   a   fennmaradó   két   oldalt   alsó   és   felső   falként  definiáltam.

4.3 NAGY ÖRVÉNY SZIMULÁCIÓ FLUENT PROGRAMMAL A szimulációt Fluenttel végeztem, melyben a lamináris és RANS modellek mellett  lehetőség van nagy örvény számításra is. A különböző hőmérséklet peremfeltételek miatt  két   számításra   van   szükség,   az   egyikben   az   átlagos   fali   hőáram,   a   másikban   a   fali  hőmérséklet állandó. A DNS eredményekkel való összehasonlíthatóság miatt a kiindulási  adatok Pr=0,71 és  Re=180. Időfüggő   szimulációt   végeztem   másodrendű   implicit   formulációval.   A   számítás  három dimenzióban, nyomás alapú implicit módszerrel történt. A nagy örvény számításhoz  a sebességek pontosítására a Smagorinsky­Lilly­féle SGS modellt használtam dinamikus  feszültségek opcióval. A Fluent periodikus peremfeltétel esetén lerögzíti a hőmérsékletet  egy   cellában   azért,   hogy   az   ne   növekedhessen   a   végtelenségig.   Viszont   ennek  következtében   mindkét   peremfeltétel   esetében   hibás   lesz   az   eredmény.  Mivel   a   (3.1)  egyenletben   már   feltételeztük,   hogy   a   közeg   összenyomhatatlan,   vagyis   a   hőmérséklet  változása nem módosítja a tulajdonságait, ebből kifolyólag  a hőmérséklet modellezhető  passzív skalárként, s ez megoldja a fenti problémat is. A hőmérsékletmezőre semmilyen  SGS modellt nem használtam, tehát a fel nem bontott skálákat teljesen elhanyagoltam. Ezután a közeg anyagjellemzőit adtam meg. Mivel a kiértékelést dimenziótlanított  mennyiségekkel   végzem,   ezért   az   egyszerűség   kedvéért   a   független   anyagjellemzőket  egységnyinek választottam. A DNS kutatások során az áramlást a periodikus falak között  22

előírt nyomáseséssel hozzák létre,   így  Re  értékét közvetlenül  be tudják állítani. Én   az  átlagos Reynolds­számot állítottam be a tömegáram meghatározásával. Az áramlás irányára  merőleges   keresztmetszet   nagysága   A=2⋅2=12,566 m2 .   Az   átlagsebességet  kg m ub =1   és   a   sűrűséget   =1 3   egységnyinek   feltételezve   a   tömegáramra   (4.8)  s m adódik: m= ˙ A⋅⋅ub =12,566

kg s

(4.8)

A (4.1) szerinti tömegáramot létrehozó nyomásgradienst a Fluent számítja ki és minden  iterációs lépés után korrigálja.

Az átlagos Reynolds­szám értékét úgy vettem fel, hogy az egyezzen a Kasagi [17]  által   találttal,   mivel   az   áramkép   validálásánál   azt   veszem   majd   referenciának.   Ennek  megfelelően  Reb=3245­öt   felvéve   a   kinematikai   viszkozitás   számítható   (4.9)­nek  megfelelően: =

u b⋅ kg =3,08166⋅10−4 ; Re b m⋅s

(4.9)

Mivel a Prandtl­szám adott, így a hőátadási tényező (4.10) szerint kifejezhető: =

 −4 kg =4,34037⋅10 Pr m⋅s

(4.10)

Fontos a peremfeltételek megfelelő megadása is, a két periodikussal nincs további  teendő, de a falakra elő kell írni első esetben az állandó hőáramot, a másodikban az állandó  hőmérsékletet. Az UHF esetében a falakon keresztül bekerülő hőt el kell vonni, különben a  hőmérséklet a számolás során végig emelkedni fog.  Erre a tanszéki fejlesztésű, C­ben írt  átlagoló   program   megfelelő   forrástag   függvényét   használtam.   Ez   annyi   hőt   von   ki  egyenletesen a térfogatból, mint amennyi a falakon beérkezik, így anélkül küszöbölhető ki  a   túlmelegedés,   hogy   a   hőszállítás   jellemzőit   befolyásolnánk,   ugyanakkor   az  x  irányú  átlaghőmérséklet­növekedést elveszítjük. A   nyomás   diszkretizálására   másodrendű,   a   hőmérsékletre   harmadrendű,   míg   a  sebességre   “kötött”   középpontos   differenciáló   (Bounded   Central   Differencing)   sémát  alkalmaztam   a   lehető   legpontosabb   eredmény   érdekében.   A   sebesség   és   nyomás  összekapcsolására   a   SIMPLE   eljárást   használtam.   Az   alulrelaxációs   faktorok  módosításával elméletileg csökkenthető a számításigény, de a próbák azt mutatták, hogy  nem gyorsul a konvergencia, ezért az alapbeállítással futtattam a számítást. Az inicializálás során a hőmérsékletet nullára állítottam, az x irányú sebességet pedig  1­re. Ezt követően a három sebességkomponenst állítottam be egy ad­hoc függvénnyel.  Ezzel egy instabil állapotot hoztam létre, ami garantálta, hogy turbulens lesz az áramkép.  Minderre azért volt szükség, mert a viszonylag alacsony Reynolds­szám miatt előfordulhat,  hogy   a   lamináris   áramlásként   inicializált   számítás   “nem   csap   át”   turbulensbe,   hanem  megmarad egy instabil lamináris állapotban. Az alkalmazott függvény a (4.11) összefüggés  23

volt: f  x , y , z =0,3⋅sin3⋅x⋅y⋅z

(4.11)

Mindkét   szimuláció   során   figyelemmel   követtem   az  y=1  egyenletű   sík   (azaz   a  középsík)   átlaghőmérsékletét,   a   falak   melletti   áramlás   irányú   nyírófeszültség   átlagos  értékét,   valamint   egy   keresztmetszetben   az  y  és  z  irányú   sebességkomponensek  négyzetösszegének a síkra vett térbeli átlagát. A CTD peremfeltétel esetében vizsgáltam  továbbá az áramlási térben előforduló legnagyobb hőmérsékletet is. Az   időlépés   nagyságát   úgy   kell   meghatározni,   hogy   a   leggyorsabban   mozgó  részecske se haladhasson egy cellánál többet egy lépés alatt. Ellenkező esetben a számítás  pontatlanná   válik.   Az   előzetes   szimulációk   alapján   úgy   találtam,   hogy   az   előforduló  m legnagyobb sebesség nagysága  ∣v∣max =1,35 .  Ennek felhasznlásával az időlépés (4.12)  s szerint meghatározható:  t= ahol 

5  1 ⋅ ≈0,16 s 72 ∣v∣max

(4.12)

5   a cellák x irányú hossza. 72

A szimuláció során először a tranzíciós szakaszt kell kiszámolni, amelynek során az  áramlás az inicializált állapotból eljut a kialakultnak tekinthető fázisba. Annak eldöntése,  hogy ez pontosan mikor következik be igen nehéz, a fentebb említett figyelemmel kísért  mennyiségek  változásából   próbáljuk   meg   meghatározni.   Mindez   azért   lényeges,   mert   a  turbulencia statisztikai jellemzőit úgy tudjuk meghatározni, hogy a vizsgált mennyiségeket  időben átlagoljuk. Amennyiben ezt túl korán kedjük el valószínűleg nagy hibát fog okozni  a tranziens szakasz jelentős változékonysága. Az átlagolást a korábban említett tanszéki  program   végzi.   Ennek   működéséhez   felhasználói   memória­helyeket   (User   Defined  Memory) kell létrehozni Fluentben, ezekben tárolódnak az átlagolt értékek. Ezeket úgy  számolja   a  kód,  hogy  minden   időlépés   végén  lekérdezi   az  adott  mennyiség  pillanatnyi  értékét és azt súlyozva hozzáadja az előző lépésbelihez. Az is fontos, hogy mennyi időn  keresztül végezzük az átlagolást. Elméletileg végtelenül hosszan kellene, általánosságban  elmondható, hogy minél tovább, annál pontosabb lesz az eredmény. Ideális minimumnak  tekinthető Kawamura [16] összefüggése, melyet a (4.13) egyenlet prezentál: t avg=

 u2

(4.13)

Ez az időmennyiség jelen esetben a csatorna kb. 50 átöblítésének felel meg. Kasagi [17]  valamivel hosszabb átlagolási időt javasol (4.14), ami közelítőleg 75 átáramlást jelent: t avg=

5 u

(4.14)

24

5 KIÉRTÉKELÉS ÉS VALIDÁLÁS 5.1 A SZIMULÁCIÓ FUTÁSÁNAK ELLENŐRZÉSE A futtatás során a számítás ellenőrzésére az előző fejezetben említett néhány alpvető  mennyiség alakulását követtem figyelemmel. A (2. ábrán) az UHF, a (3. ábrán) a CTD  peremfeltételhez tartozó eredmények időlépésenkénti alakulása látható.

(2.a ábra: Az y és z irányú sebességkomponensek négyzetösszegének a síkra vett térbeli átlaga)

25

(2.b ábra: Az y=1 egyenletű középsík átlaghőmérséklete)

(2.c ábra: Az x irányú fali nyírófeszültség két falon vett átlaga)

26

(3.a ábra: Az y és z irányú sebességkomponensek négyzetösszegének a síkra vett térbeli átlaga)

(3.b ábra: Az y=1 egyenletű középsík átlaghőmérséklete)

27

(3.c ábra: Az x irányú fali nyírófeszültség két falon vett átlaga)

(3.d ábra: Az áramlási térben előforduló legmagasabb hőmérséklet)

28

A   CTD   peremfeltétel   esetén   az   átlagolás   során   vizsgáltam   az   áramlási   térben  előforduló legmagasabb hőmérsékletet (lásd 3.c ábra), de nem találtam a Chatelain et al.  [11] által jelentett irreális értékhez hasonló jelenséget. Valószínű, hogy ez a numerikus  eljárásukból adódó hiba volt. Először   a   tranzíciós   szakaszt   számítottam   ki,   aminek   hossza   közelítőleg   1000  időlépésnek,   azaz   körülbelül   10   átáramlásnak   adódott.   Ezt   követően   kapcsoltam   be   a  tanszéki   átlagoló   programot,   melyet   további  13000   időlépésen   keresztül   futtattam.   Ez  2080s áramlási időnek és kb. 130 teljes átáramlásnak felel meg. Ellenőrzés   szempontjából   fontos   mennyiségek   továbbá   az  x  irányú   sebesség,   a  csatornán keresztül haladó tömegáram, a fal melletti első cellában y+ értéke és a turbulens  Reynolds­szám.   Az   első   kettőnek   a  t=2240s  időpillanathoz   tartozó,   utóbbi   kettőnek   az  átlagolás során kapott értékét, és ezek célértékét tüntettem fel az (1. táblázat)­ban: (1. táblázat: Ellenőrzött mennyiségek)

x irányú sebesség [m/s] Peremfeltétel

Valós

Cél

UHF

1,000012

1

CTD

0,999959

1

Tömegáram [kg/s] Valós

Cél

y+ az első  cellában Valós

Reτ

Cél

Valós

Cél

12,56599 12,566 0,9257

1

177

180

12,56597 12,566 0,9272

1

178

180

Lényeges   még   a   Courant­szám   cellákban   vett   értékének   eloszlása   is.   Ezzel  ellenőrizhető a számítás időbeli felbontása. Az adott áramlás esetén ennek az átlaga egy  körüli   kell,   hogy   legyen.   Továbbá   törekedni   kell   arra,   hogy   minél   kevesebb   cellában  kapjunk egy fölötti értéket. A (4.a ábrán) látható hisztogram prezentálja a Courant­szám  eloszlását a teljes vizsgált tartományban. Látható, hogy az előbbi feltételeknek megfelelő  mértékben sikerült eleget tenni, az átlagos Courant­szám C=0,8893376, maximális értéke  Cmax=1,664092  és a cellák több, mint 95%­a 1,2 alatt van. A (4.b ábra) mutatja az y+  értékének alakulását a két fal mellett. Maximuma  y+max=1,549662; átlaga  y+=0,9190032;  egy cellában sem kisebb az értéke mint 0,6; tehát a fal melletti réteg felbontása megfelelő.

29

(4.a ábra: A cellák Courant­számának eloszlása a teljes áramlási térre vonatkoztatva)

(4.b ábra: y+ megoszlása a falak melletti első cellarétegekben)

30

5.2 FELTÉTELES ÁTLAGOLÁS A számítás során a tanszéki program segítségével ún. feltételes átlagolást végeztem,  melyről bővebb leírás található Lohász et al. [16] cikkében. A feltétel a Q­kritérium volt,  melynek   segítségével   a   vizsgált   mennyiségeket   négy   különböző   csoportra   bontottam.  Annak függvényében, hogy egy cella az adott pillanatban melyik  Q  osztályba tartozik, a  változókat külön átlagolja. Így minden csoportra külön eredmény adódik, melyeket az adott  osztályba való tartozás idejével súlyozni kell, majd a vizsgált mennyiségeket csoportonként  összegezni. Az így kapott időbeli átlagokat térben vízszintes síkonként újból átlagoltam.  Az   ilyenformán   adódó   eredményeket   használtam   a   sebesség­   és   hőmérséklet­profilok  elkészítéséhez. A Q­kritériumról bővebben az 5.5 fejezetben lesz szó.

5.3 SEBESSÉGMEZŐ VALIDÁLÁSA A   sebességek   vizsgálatakor   Kawamura   [19]   és   Kasagi   [20]   DNS   eredményeihez  hasonlítom a számításomat. A sebességprofilok a hőmérséklet peremfeltételtől függetlenek,  ezért itt nem közlöm őket külön­külön. Az  x  irányú átlagsebességet  u­val dimenziótlanított formában és a fali koordináták  függvényében a (5. ábra) mutatja be:

31

Átlagsebesség 25

u_mean+

20

15

10

LES Kasagi Kawamura

5

0 0

20

40

60

80

100

120

140

160

180

y+ (5.a ábra: Dimenziótlan sebességprofil a fali koordináták függvényében)

Átlagsebesség 25

u_mean+

20

15

10

LES Kasagi Kawamura

5

0 0.18

1.8

18

180

y+

(5.b ábra: Logaritmikus dimenziótlan sebességprofil a fali koordináták függvényében)

32

A   három   sebességkomponens   rms   értékét   a   fali   koordináták   függvényében   a   (6.  ábra) szemlélteti, a dimenziótlanítás itt is  uτ­val történt.

Sebesség-ingadozások rms értéke 3.5

u'+ v'+ w'+ u'+ (Kasagi) v'+ (Kasagi) w'+ (Kasagi) u'+ (Kawamura) v'+ (Kawamura) w'+ (Kawamura)

3

u'+ v'+ w'+

2.5

2

1.5

1

0.5

0 0

20

40

60

80

100

120

140

160

180

y+ (6. ábra: Sebességkomponensek ingadozásának rms értéke  a fali koordináták függvényében)

A   két   előző   ábrából   jól   látható,   hogy   a   profilok   jellege   megfelelő,   azonban  pontosságuk hagy kívánni valót maga után. A maximális relatív hiba megközelíti a 25%­ot.  Ennek oka feltehetően az, hogy ilyen kis cellaszám mellett nagy örvény szimulációval csak  kisebb   tartományt   vizsgálva   kapható   pontosabb   eredmény.  A   kisebb   modellezett   tér  azonban felveti a kérdést, hogy elegendően nagy­e ahhoz, hogy a turbulencia hosszléptékét  magába foglalja. A (7. ábra) egy pillanatnyi hőmérséklet­eloszlást mutat a két falon. Jól  megfigyelhetők a felső lapon lévő, a tartomány hosszának felénél határozottan nagyobb  örvények, melyekből arra következtethetünk, hogy ezek alapján  3­nál rövidebbre nem  vehető a vizsgált térfogat hosszúsága.

33

(7. ábra: Pillanatnyi hőmérséklet­eloszlás a falakon UHF peremfeltétel esetén)

5.4 A HŐMÉRSÉKLET­ELOSZLÁS ÖSSZEHASONLÍTÁSA 5.4.1 Az UHF peremfeltétel A   transzformált   hőmérséklet   eloszlása  Tτ­val   dimenziótlanítva   a   csatorna  keresztmetszetében a (8. ábrán) látható.

34

Hőmérséklet-eloszlás 20 18 16 14

T_mean

12 10

T_mean T_mean (Kasagi) T_mean (Kawamura)

8 6 4 2 0 0

20

40

60

80

100

120

140

160

180

y+

(8.a ábra: Dimenziótlanított hőmérséklet­eloszlás  a fali koordináták függvényében)

Hőmérséklet-eloszlás 20 18 16 14

T_m ean

12 10 8

T_mean T_mean (Kasagi) T_mean (Kawamura)

6 4 2 0 0.18

1.8

18

180

y+

(8.b ábra: Logaritmikus dimenziótlanított hőmérséklet­eloszlás  a fali koordináták függvényében)

35

A transzformált hőmérséklet fluktuációja a (9. ábrán) látható. A DNS eredményektől  való   eltérés  a  különböző   peremfeltételből   adódik.  Jól  látható,   hogy  az   általam  vizsgált  esetben   a   hőmérséklet   ingadozása   a   falon   nem   nulla,   mivel   nem   használtam   azt   az  egyszerűsítést, hogy a fal hővezetése végtelen nagy. Így nagyobb gradiensek alakulhatnak  ki, melyek hatékonyabb hővezetést tesznek lehetővé a fal közelében, ami magyarázza a  görbék eltolt helyzetét.

Hőmérséklet-ingadozások rms értéke 3.5

T' T' (Kasagi) T' (Kawamura)

3

2.5

2

T'

1.5

1

0.5

0 0

20

40

60

80

100

120

140

160

180

y+ (9. ábra: Transzformált hőmérsklet ingadozásának rms értéke Tτ­val dimenziótlanítva a fali koordináták függvényében)

36

5.4.2 A CTD peremfeltétel A transzformált hőmérséklet eloszlása a csatorna keresztmetszetében a (10. ábrán)  látható.

Hőmérséklet-eloszlás 1 0.9 0.8 0.7

T_mean

0.6 0.5 0.4

T_mean T_mean (Seki)

0.3 0.2 0.1 0 0

0.5

1

1.5

2

y

(10. ábra: Transzformált hőmérséklet­eloszlás  a csatorna keresztmetszetében)

37

A transzformált hőmérséklet ingadozását a (11. ábra) szemlélteti.

Hőmérséklet-ingadozás rms értéke 0.09 0.08

T'

0.07

CTD

0.06 0.05 0.04 0.03 0.02 0.01 0 0

20

40

60

80

100

120

140

160

180

200

y+

(11. ábra: Transzformált hőmérsklet ingadozásának rms értéke Tτ­val dimenziótlanítva a fali koordináták függvényében )

5.4.3 Korrelációk A két peremfeltétel kereszt­korrelációit a könnyebb áttekinthetőség kedvéért közös  grafikonban ábrázolom. A (12. ábrán) az x és y irányú sebességek korrelációja látható. A csatorna belsejében  kb.  y =60­nál   nagyobb   értékeknél   a   nyírófeszültség   már   lineárisan   közelíthető,   mely  analitikusan levezethető. Az általam kapott eredmény ennél magasabb, ami arra utal, hogy  mivel   az   SGS   feszültséggel   nincsenek   korrigálva   az   adatok,   ott   annak   negatívnak   kell  lennie. +

38

Sebesség korreláció 0.9

uv+ uv+ (Kasagi) uv+ (Kawamura)

0.8 0.7 0.6

uv+

0.5 0.4 0.3 0.2 0.1 0 0

20

40

60

80

100

120

140

160

180

y+ (12. ábra: x és y irányú sebességkomponensek korrelációja)

Az  x  irányú sebesség és a hőmérséklet korrelációját a (13. ábra), míg az  y  irányú  sebesség és a hőmérséklet korrelációját a (14. ábra) szemlélteti.

Sebesség-hőmérséklet korreláció 1.2 1 0.8

u-T

0.6

CTD UHF

0.4 0.2 0 0

20

40

60

80

100

120

140

160

180

200

y+

39

(13. ábra: u'' korreláció a fali koordináták függvényében)

Sebesség-hőmérséklet korreláció 0.6

0.5

0.4

CTD UHF

v-T

0.3

0.2

0.1

0 0

20

40

60

80

100

120

140

160

180

200

y+ (14. ábra: v'' korreláció a fali koordináták függvényében)

Az   UHF   és   a   CTD   peremfeltételek   között   a   legnagyobb   különbség   a   csatorna  közepénél   mutatkozik.   Állandó   hőáram   esetén   az  v''  kereszt­korrelációs   együttható   a  középsíknál nullához tart, mert a falra merőleges irányú turbulens hőáram inverz módon  szimmetrikus   a   csatorna   középsíkjára.   Az   állandó   hőmérséklet­különbség   esetén  u''  hasonló okok miatt tart nullához. A másik két esetben, az UHF­nél u'', míg CTD­nél v''  közel állandó marad egészen a csatorna közepéig. Ennek oka, hogy a megfelelő irányú  turbulens hőáram is konstans. Mivel a korábban többször említett tanszéki átlagoló az UHF peremfeltétel esetére  lett   kifejlesztve,   ezért   működése   korlátozott   volt   a   CTD   számítása   során.   Ennek  következtében   a   (13.   és   14.   ábra)   elkészítéséhez   lineáris   transzformációkat   kellett  alkalmaznom. A   (15.   ábra)   a   sebességmező   és   a   hőmérséklet   pillanatnyi   eloszlását   mutatja   a  t=2240s  időpontban.   A   két   kép   inverz   hasonlóságából   jól   látható,   hogy   a   konvektív  hővezetés egyértelműen domináns, a konduktív hőtranszport szinte elhanyagolhatóan kicsi.

40

(15.a ábra: Pillanatnyi sebességeloszlás a periodikus falakon)

(15.b ábra: Pillanatnyi hőmérséklet­eloszlás a periodikus falakon)

41

5.5 KOHERENS STRUKTÚRA A   koherens   struktúra   alapötlete,   hogy   a   sebességkomponensek   ingadozó   tagját  felbontják   egy  koherens   mozgásra   és  egy   turbulens  háttérre.  Az  az   elképzelés,  hogy   a  koherens   mozgás   sokkal   könnyebben   leírható,   mint   az   utóbbi,   ugyanakkor   dinamikai  szempontból is meghatározóbb. Jeong & Hussian [16] szerint a koherens struktúrát térben  koherens, időben fejlődő örvénylő mozgás alkotja. A turbulencia koherens struktúráinak vizsgálatakor a  Q­kritériumot használtam az  örvények azonosítására, ahol Q az (5.1) egyenlet szerint a sebességgradiens tenzor második  invariánsa: 1 Q= ij  ji −S ij S ji  2

(5.1)

ahol ij az örvénytenzor, melyet definíció szerint az (5.2) összefüggés ír le: 1 ij = ∂ i u j−∂ j ui  2

 (5.2)

Sij pedig az alakváltozási tenzor, mely az (5.3) szerint fejezhető ki: 1 S ij = ∂ i u j ∂ j u i 2

 (5.3)

Bővebb leírásért lásd [17], valamint [18], utóbbiban más leírási módszerek (például  2) részletes bemutatása is megtalálható. A   szimuláció   során   négy  Q  osztályt   definiáltam,   ezek   határait   ad­hoc   módon   a  következőkben állapítottam meg: Első osztály:

­

Második osztály:

0      <  Q < 0,01

Harmadik osztály:

0,01  <  Q < 0,1

Negyedik osztály:

0,1

< Q < 0

<  Q < 

A (16. ábrán) adott pillanatban kialakuló Q= állandó felületek láthatóak.

42

(16.a ábra: Q iso­felületek; kék: Q=0; zöld: Q=0,01)

(16.a ábra: Q iso­felületek; sárga: Q=0,1; piros: Q=0,2)

43

A   (2.   táblázat)   a  Q  osztályokba   tartozás   százalékos   arányát   mutatja,   a   két  peremfeltétel mellett feltüntettem Lohász Máté [17] eredményeit is. Ez részben megerősíti a  konzulensem által felfedezetteket, miszerint az örvényesség nem a leggyakoribb folyadék  állapot a turbulenciában. Ő úgy találta, hogy a  blokkolt  csatornaáramlásban az örvények  valószínűsége 40%, lásd [17] 143­166 oldalán. Meglepő módon a  bordát  nem tartalmazó  csatornában ugyanazzal a számítási módszerrel ennél valamivel magasabb értéket kaptam,  42%­ra   adódott   az   örvények   aránya   az   áramlásban.   Az   előzőek   ellent   mondanak   az  intuíciónak, miszerint a  borda  felkavarja az áramlást, sok új örvényt létrehozva. Mindez  látszólag arra utalhat, hogy jobb hűtést lehet elérni borda nélkül, de ez nyilvánvalóan nem  igaz. Ez valószínűleg azzal magyarázható, hogy a konzulensem által vizsgált csatornában  az   átlagos   örvények   nagyobb   intenzitásúak,   amik   jobb   hatásfokkal   szállítják   el   a   hőt.  Sajnos   ennek   vizsgálata   nem   megoldható,   mert   a  Q  osztályokra   előírt   határok   nem  konzekvensek, ad­hoc módon lettek meghatározva. Továbbá nem szabad megfeledkezni a  borda felületnövelő hatásáról sem, ami szintén a jobb hőtranszportot segíti. (2. táblázat: Q osztályok százalékos aránya a különböző peremfeltételek esetén)

Vizsgált eset

QI [%]

QII [%]

QIII [%]

QIV [%]

UHF

57,79

9,96

22,04

10,20

CTD

57,83

9,90

21,93

10,35

45° 1

59,7

38,8

1,4

0,09

90° 1

60

37,6

1,9

0,08

1: lásd [17] 147. oldalán

A (3. táblázatban) található a  Q  osztályok átlagos  x  irányú sebessége. Jól látható,  hogy minél nagyobb intenzitású egy örvény, annál lassabban halad az áramlás irányában. (3. táblázat: Q osztályok átlagos x irányú sebessége)

Peremfeltétel

QI

QII

QIII

QIV

UHF

1,0032

1,0254

0,9992

0,9548

CTD

1,0032

1,0253

0,9992

0,9546

44

6 ÖSSZEFOGLALÁS A szakdolgozatom során megismerkedtem a direkt numerikus szimuláció és a nagy  örvény szimuláció alapvető elméleti hátterével, a csatornaáramlásban való alkalmazásuk  fontosabb   irodalmi   eredményeivel,   valamint   elsajátítottam   a   periodikus   hőátadás  számításának módját. Ezt   követően   nagy   örvény   szimuláció   segítségével   vizsgáltam   végtelen   síklapok  közötti   kialakult   turbulens   áramlásban   a   hőtranszportot   különböző   termikus  peremfeltételek esetén. A sebesség­ és hőmérsékletmező főbb mennyiségeinek pillanatnyi  értékeit   valamint   időbeli   és   térbeli   átlagait   is   tanulmányoztam.   A   kapott   eredmények  alátámasztják a nagy örvény szimuláció létjogosultságát, ugyanakkor a kapott nagy relatív  hiba miatt a további eredmények fenntartással kezelendőek. A hiba oka feltehetően az adott  cellaszám mellett vizsgált túl nagy tartomány volt, ezt célszerű volna újabb számításokkal  ellenőrizni. A   kiértékelés   során   feltételes   átlagolást   végeztem   a  Q­kritérium   alapján.   Az  örvényesség vizsgálatakor kimutattam, hogy a turbulens áramlásban nem az örvények a  leggyakoribbak,   valamint   hogy   minél   nagyobb   intenzitású   egy   örvény,   az   áramlás  irányában annál lassabban halad.

45

IRODALOMJEGYZÉK [1] Batchelor, G. K., 1959, “Small­Scale Variation of Convected Quantities Like  Temperature   in   Turbulent   Fluid:   Part   1.   General   Discussion   and   the   Case   of   Small  Conductivity” J. Fluid Mech., Vol. 5, pp. 113­133. [2] Tennekes, H., Lumley, J.L., 1972, “A First Course in Turbulence”  MIT Press,  Cambridge, MA, pp. 96­97.  [3]   Kasagi,   N.,   Iida,   O.,   1999,   “Progress   in   Direct   Numerical   Simulation   of  Turbulent Heat Transfer ” Proceedings of the 5th ASME/JSME Joint Thermal Engineering   Conference , San Diego, California [4] Kim, J., Moin, P., 1989, “Transport of Passive Scalars in a Turbulent Channel  Flow” Andre, J. ­C. (Ed.), Turbulent Shear Flows 6, Springer, Berlin, pp. 85–96. [5] Kasagi, N., Tomita Y. and Kuroda A., 1992, “Direct Numerical Simulation of  Passive Scalar Field in a Turbulent Channel Flow” ASME J. Heat Transfer, 114:598–606. [6] Kasagi  N.,  Ohtsubo Y., 1993, “Direct  Numerical  Simulation  of Low   Prandtl  Number Thermal Field in a Turbulent Channel Flow”  Durst, F. et al. (Eds.), Turbulent   Shear Flows 8, Springer, Berlin, pp. 97–119. [7] Kawamura, H., Ohsaka, K., Abe H. and Yamamoto K., 1998, “DNS of Turbulent  Heat Transfer in Channel Flow with Low to Medium­high Prandtl Number” Int. J. of Heat  and Fluid Flow, 19:482–491.  [8] Kawamura, H., Abe H. and Matsuo Y., 1999, “DNS of Turbulent Heat Transfer  in Channel Flow with Respect to Reynolds and Prandtl Number Effect” Int. J. of Heat and  Fluid Flow, 20:196–207.  [9] Seki, Y., Abe, H., Kawamura H., 2003, “DNS of Turbulent Heat Transfer in a  Channel   Flow   with   Different   Thermal   Boundary   Conditions   ”  The   6th   ASME­JSME  Thermal Engineering Joint Conference [10] Fukagata, K., Iwamoto K. and Kasagi, N., 2005, “Novel Turbulence Control  Strategy   For   Simultaneously   Achieving   Friction   Drag   Reduction   and   Heat   Transfer  Augmentation   ”  Proc.   4th   Int.   Symp.   Turbulence   and   Shear   Flow   Phenomena   Williamsburg, Virginia, pp. 307­312.  [11]   Chatelain,   A.,   Ducros,   F.   and   M´tais,   O.,   2004,   “LES   of   Turbulent   Heat  Transfer:   Proper   Convection   Numerical   Schemes   for   Temperature   Transport   ”  Int.   J.  Numer. Meth. Fluids; 44:1017–1044 [12] Montreuil, E., O. Labbé, O. and Sagaut, P., 2005, “Assessment of non­Fickian  Subgrid­Scale Models for Passive Scalar in a Channel Flow ” Int. J. Numer. Meth. Fluids,  49:75–98 [13] Patankar, S. V., Liu, C. H. and Sparrow, E. M., 1977, “Fully Developed Flow  and   Heat   Transfer   I   Ducts   Having   Streamwise­Periodic   Variations   of   Cross­Sectional  Area” Journal of Heat Transfer, 99:180­186 46

[14] Mathieu, J., Scott, J., 2000, “An Introduction To Turbulent Flow”  Cambridge  University Press, pp. 23­27 [15]   Kim,   S­E,   2004,   “Large   Eddy   Simulation   Using   Unstructered   Meshes   and  Dynamic   Subgrid­Scale   Turbulence   Models”,  American   Institute   of   Aeronautics   and   Astronautics Paper, 2004­2548 [16] Jeong, J., Hussian, F., 1995, “On the Identification of a Vortex”,  Journal of  Fluid Mechanics, vol. 285, pp. 69­94 [17] Lohasz, M., M., Rambaud, P., Benocci, C., 2005, “Conditional Averaging of the  Fully Developed Stationary Ribbed Duct Flow Using Q Criteria ”, Proceedings of the iTi  Conference in Turbulence, pp. 285­288 [18] Lohasz, M., M., 2008, “Large Eddy Simulation of Heat Transfer in Ribbed  Ducts”, Ph.d. Thesis, pp. 40­44

INTERNETES HIVATKOZÁSOK [19] http://murasun.me.noda.tus.ac.jp/turbulence/poi/text/Poi180_4th_C.dat [20] http://www.thtlab.t.u­tokyo.ac.jp/DNS/CH122_PG.WL3

47