A RÁCS-BOLTZMANN MÓDSZER ALKALMAZÁSA EGY- ÉS KÉTFÁZISÚ ÁRAMLÁSI PROBLÉMÁK MODELLEZÉSÉRE PhD tézisfüzet MAYER GUSZTÁV MTA, KFKI, Atomenergia Kutatóintézet
Témavezet˝o: DR. HÁZI GÁBOR, tudományos f˝omunkatárs, MTA, KFKI, AEKI Tanszéki konzulens: DR. ASZÓDI ATTILA, igazgató, BME, NTI
BUDAPEST 2009
1
A kutatások el˝ozménye A Three Mile Island-i atomer˝om˝u 2. blokkjának (TMI-2) 1979-ben történt balesete rámutatott a nukleáris er˝om˝uvekben zajló kétfázisú áramlások megismerésének jelent˝oségére. Ebben az id˝oszakban kezdték el fejleszteni a szimulátorokban és a biztonsági analízisekben használt kétfázisú, egydimenziós rendszerkódokat. 2000-ben felmerült az igény, hogy az AEKI is rendelkezzen egy saját fejlesztés˝u, valós idej˝u, nagylépték˝u, a paksi szimulátorba illeszthet˝o rendszerkóddal. Ennek a kódnak (RETINA) [10, 11, 12, 13] a fejlesztése során, melyben magam is részt vettem, számos problémával szembesültünk. A tesztek rávilágítottak arra, hogy a kétfázisú kódokban használt korrelációk pontossága sok esetben megkérd˝ojelezhet˝o. A legtöbb probléma a két fázist elválasztó határfelületen zajló tömeg, impulzus és energia transzfer modellezéseinél adódik [10]. A RETINA kód kifejlesztése után figyelmem az egy- és kétfázisú termohidraulikai folyamatok numerikus leírásának részletesebb megismerésére irányult. Napjainkban, az intenzíven fejl˝od˝o számítógépek korszakában, egyre nagyobb szerephez jutnak a nagy számításigény˝u folyadékáramlásokat modellez˝o programrendszerek. Ezeket a számításos folyadékdinamikát modellez˝o programkódokat az angol terminológia után a szaknyelvben általánosan CFD (Computational Fluid Dynamics) kódoknak nevezik.
Vizsgálati módszerek Az elmúlt évtizedben egy új numerikus módszer látott napvilágot a folyadékáramlások szimulációjának területén, nevezetesen a rács-Boltzmann módszer, vagy ahogy az angolszász irodalomban olvasható, a „Lattice Boltzmann Method” - LBM. A módszer a Boltzmann egyenletb˝ol közvetlenül származtatható és alkalmas a Navier-Stokes egyenletek numerikus megoldására. Egyik el˝onye, hogy viszonylag könnyen kiterjeszthet˝o kétfázisú áramlási problémák vizsgálatára, így segítségével a két fázist elválasztó határfelületen zajló folyamatok részletesen vizsgálhatók. További el˝onye még, hogy könnyen programozható és alkalmas komplex geometriák vizsgálatára is [Treeck és társai, 2007]. Az el˝obb vázolt kedvez˝o tulajdonságai miatt választottam dolgozatom központi témájának a rács-Boltzmann módszert. Kutatásom során, használatával olyan kérdésekre kerestem a választ, amelyek részint a módszer jellegzetességeihez, részint a nukleáris iparban zajló alapkutatásokhoz kapcsolódnak.
Célkituzések ˝ Turbulens folyadékdinamikai szimulációkat alapvet˝oen három f˝o szinten lehet megvalósítani. Ezen szintek a számítás részletességében és pontosságában térhetnek el jelent˝osebb mértékben egymástól. Mindegyik közös tulajdonsága, hogy a Navier-Stokes egyenleteket oldják meg, egy arra alkalmas numerikus módszerrel. Az els˝o ilyen szintet direkt numerikus szimulációnak nevezik. Lényege, hogy közvetlenül oldja meg a Navier-Stokes egyenleteket, bárminem˝u egyszer˝usítés nélkül, mégpedig úgy, hogy az áramlást a legkisebb örvények méretéig felbontja. Ez a módszer adja a legpontosabb megoldást, ugyanakkor ez igényli a legnagyobb számítási kapacitást is. Emiatt segítségével csak alacsony Reynolds számú turbulens áramlások vizsgálhatók. A második megközelítésmódot nagyörvény szimulációnak nevezik. Jellemz˝oje, hogy még mindig igen jól felbontja az áramlást, és a felbontott skálák alatt - ott ahol a folyamatok univerzálisan viselkednek - általában egy egyszer˝u turbulencia modellt alkalmaz. A harmadik és egyben legelterjedtebben használt (gyakorlati szempontból legjelent˝osebb) szint - turbulens áramlások modellezésére - az ún. Reynolds Átlagolt Navier-Stokes (RANS) szimuláció. El˝onyei közé tartozik, hogy segítségével igen komplex és bonyolult geometriák is modellezhet˝ok a ma rendelkezésre álló er˝oforrásokkal. Hátránya viszont az, hogy - a pontosság biztosítása érdekében - nagy figyelmet kell fordítani a turbulenciamodellek, illetve ezen modellek paramétereinek megválasztására. 2
A rács-Boltzmann módszerrel már mind a három - az el˝obb említett - szinten végeztek sikeres szimulációkat. Munkám alapvet˝o célkit˝uzése egy olyan két- és háromdimenziós rács-Boltzmann alapú kód kifejlesztése volt, amely a szóban forgó szintek közül alkalmas direkt numerikus és nagyörvény szimuláció megvalósítására. A nukleáris iparban nagy jelent˝osége van a kétfázisú áramlások modellezésének, így további célként szerepelt egy olyan kód kidolgozása is, amely alkalmas kétfázisú folyadékok direkt numerikus szimulációjára is. Az elkészített kódokat m˝uködésük helyességének szempontjából ellen˝orizni kell. Ennek megfelel˝oen végtelen kiterjedés˝u (periodikusan csatolt) sík lapok között gyorsítottam a folyadékot, melyben a gyorsítóer˝o vektora párhuzamos volt a falakkal. Ezt a tesztfeladatot azért alkalmazzák széleskör˝uen, mert egyszer˝u analitikus megoldása van lamináris esetben, ezáltal a kód m˝uköd˝oképessége könnyen ellen˝orizhet˝o. Bár a rács-Boltzmann módszerrel el˝ottem többen is végeztek sikerrel direkt numerikus szimulációt (a módszer alkalmazhatósága erre a szintre bizonyított volt), azonban a kód m˝uköd˝oképességének demonstrálásához, célul t˝uztem ki végtelen síklapok közötti szimulációk elvégzését turbulens esetre is. További érv a síklapok közötti direkt numerikus szimulációk elvégzésére, hogy ez az egyik legjobban dokumentált tesztfeladat az irodalomban. A VVER-440 típusú f˝ut˝oelemköteg belsejében zajló termohidraulikai folyamatok mind a mai napig intenzíven kutatott témája a nukleáris iparnak. A valós rendszer geometriája igen bonyolult és a benne zajló áramlás Reynolds száma is 200000 feletti, így egy teljes köteg finom skálás modellezésére nincs mód. Lehet˝oség van azonban egy VVER-440 típusú geometriához hasonló végtelen kiterjedés˝u pálcakötegszakasz vizsgálatára az üzemi Reynolds számnál jóval alacsonyabb Reynolds számok mellett. A kifejlesztett kóddal megvizsgáltam a fent említett geometriában az áramlás f˝obb tulajdonságait. El˝oször lamináris esetre ellen˝oriztem a kódot. A következ˝okben a pálcakötegszakasz direkt numerikus szimulációját végeztem el alacsony - de már turbulens - Reynolds szám tartományban. A f˝o kérdés az volt, hogy visszaadja-e az áramlás alapvet˝o tulajdonságait (másodlagos áramlás, áramlási pulzáció) a szimuláció. Ezt követ˝oen magasabb (12000 és 24000) Reynolds számok mellett végeztem nagyörvény-szimulációs vizsgálatokat. Ezen Reynolds számok mellett mérési adatok is rendelkezésemre álltak, így azok jó támpontot nyújtottak az eredmények ellen˝orzésére. Turbulens esetben nem csak a sebességek várhatóértékei, hanem a fluktuáló komponensek szorzatának várhatóértékei is fontos jellemz˝ok. Ezek a Reynolds-féle feszültségtenzor elemei, melyek vizsgálatát szintén elvégeztem. A kétfázisú numerikus vizsgálatok esetében az egyik legfontosabb tesztfeladat annak igazolása, hogy visszakapható-e a helyes fázisegyensúlyi görbe a szimulációk során. Ennek vizsgálatához - metastabil állapotból kiindulva - kezdetben szükséges egy kis s˝ur˝uségperturbáció, mely hatására fázisszeparáció jön létre. A rendszerben állandósult állapotban kialakulnak maximális és minimális s˝ur˝uség˝u pontok, melyek a fázisegyensúlyi görbén két pontot adnak adott h˝omérséklet mellett. A rács-Boltzmann módszerekben a legáltalánosabban használt kétfázisú módszert Shan és Chen [Shan és Chen, 1993] publikálta 1993-ban. Azonban Qian és Chen [Qian és Chen, 1997], valamint Imre és Házi [Imre és Házi, 2002] egy, illetve két dimenzióban megfigyelték, hogy a kialakult fázisegyensúlyi görbe még periodikus peremfeltétel használata esetében is függ az alkalmazott felbontástól. Megállapították, hogy a Shan-Chen modell csak adott felbontás felett ad pontos fázisegyensúlyi görbét. Felmerült tehát annak igénye, hogy a kifejlesztett kétfázisú, háromdimenziós programrendszerrel megvizsgáljam ezt a problémát három dimenzióban. Ezen vizsgálataimban arra kerestem a választ, hogy vajon jelentkezik-e a véges felbontás hatása a háromdimenziós rács-Boltzmann szimulációk során, és ha igen, akkor mekkora minimális felbontást kell választani ennek elkerüléséhez. A kutatók már az els˝o kétfázisú rács-Boltzmann modellek kapcsán megfigyelték, hogy a kialakult folyadék-g˝oz határfelület relatíve szélesre adódik. Például egy Shan-Chen modellel szimulált 20 rácsegység átmér˝oj˝u buborék esetében stabilitási okokból nemigen lehet három vagy négy rácstávolságnál keskenyebb határfelületet modellezni. A gyakorlatban azonban egy makroszkopikus méret˝u buborék esetében a folyadék-g˝oz határfelület - a kritikus pont kis környezetét kivéve - a nanométeres tartományba 3
esik. Felmerült tehát az igény annak a kérdésnek a megválaszolására, hogy a kialakult határfelület profilja megegyezik-e az elméleti, illetve kísérleti adatokkal, és ha igen, akkor a határfelület vastagságát mérföldk˝onek tekintve mekkora valós fizikai méretnek felel meg egy rács a valóságban.
Új tudományos eredmények tézisei 1. Els˝oként végeztem el egy végtelen kiterjedés˝u háromszög elrendezés˝u pálcaköteg szubcsatorna szakaszának direkt numerikus szimulációját a D3Q19-es [1, 2] és a D3Q27-es rácsmodellel [2, 3, 4], 3680-as Reynolds szám mellett az LBM3D1P, saját fejlesztés˝u, egyfázisú, háromdimenziós, rács-Boltzmann alapú programkóddal. a.
A program helyes m˝uködésének ellen˝orzéséhez analitikus megoldásokkal vetettem össze a rács-Boltzmann módszer eredményeit lamináris esetben. A számítások során mind a D3Q19-es, mind a D3Q27-es rácsmodellekkel kapott eredményeket megvizsgáltam [2, 4]. Mindkét rácsmodell esetén azt tapasztaltam, hogy az axiális sebesség relatív hibája a fal mellett a legnagyobb, és a faltól távolodva jelent˝osen csökken. A maximális hiba nagysága 2% alatti.
b.
A leellen˝orzött programmal direkt numerikus szimulációt végeztem 3680-as Reynolds számú turbulens áramlás esetén. Az id˝oben átlagolt laterális sebességek másodlagos áramlást mutattak [2, 3, 4].
2. Els˝oként végeztem el egy végtelen kiterjedés˝u háromszög elrendezés˝u pálcaköteg szubcsatorna szakaszának nagyörvény szimulációját a D3Q19-es [1, 2] és a D3Q27-es rácsmodellel [1, 2, 3, 4, 5, 6], 12000-es és 24000-es Reynolds számok mellett az LBM3D1P programrendszerrel. a.
A keresztirányú sebességek várhatóértékei a D3Q27 rácsmodell esetében másodlagos áramlási cellákat mutattak a szubcsatorna egy-egy harminc fokos szegmensében, a mérésekkel megegyez˝oen [2, 3, 4, 5, 6].
b.
A D3Q19-es modellel végzett szimulációk esetében az axiális sebességek várhatóértéke a - D3Q27-es modellel ellentétben - min˝oségileg helytelen eredményt szolgáltatott [2, 3]. A szimuláció nem adta vissza a rendszer geometriája által el˝oírt 60 fokos szimmetriát. Ez mind a nagyörvény, mind a direkt numerikus szimuláció esetében megfigyelhet˝o volt, így a nagyörvény szimulációban alkalmazott turbulenciamodellnek a rossz várhatóértékben játszott szerepe kizárható.
c.
A Reynolds feszültségek keresztirányú (UV) komponense - a falak sz˝uk környezetét kivéve - jól követi a Trupp és Azad [Trupp és Azad, 1975] által publikált mérési eredményeket [5, 6]. Megállapítható, hogy a felbontás növelésével a pontosság növekszik. A f˝oirányú normál feszültségek (RMS - root mean square) a falaktól távolabb szintén jó egyezést mutatnak a mérésekkel, azonban a fal melletti tartományban - a falfüggvény, illetve a fal melletti felbontás hiányában - jelent˝os eltérés figyelhet˝o meg.
3. Stacioner egyensúlyi állapotban, három dimenzióban, egy x-y síkban kialakuló, z tengelyre mer˝oleges, végtelen kiterjedés˝u folyadékfilm folyadék-g˝oz határfelületét modelleztem a rács-Boltzmann módszerrel. A vizsgálatokhoz az LBM3D2P, háromdimenziós, kétfázisú modellezésre is alkalmas programkódot alkalmaztam, amely a Shan-Chen modellt használja [7, 8]. a.
A numerikus számítások azt mutatták, hogy a rács-Boltzmann módszerrel kapott folyadékg˝oz határfelület profilja igen jó egyezést mutat az elméleti tangens hiperbolikus profillal és a kísérleti eredményekkel. 4
b.
Mindez lehet˝ové tette, hogy a határfelület vastagságát felhasználva meghatározzam egy cella méretét a rács-Boltzmann módszerben. Összevetettem az általam kapott eredményeket argonra végzett molekuláris dinamikai szimulációkkal. Azt találtam, hogy egy rácsBoltzmann cella mérete hozzávet˝olegesen 0,36 nm.
c.
Abban az esetben, amikor a határfelület vastagságát figyelembe kell venni, a Shan-Chen rács-Boltzmann módszerrel csak a kritikus pont közvetlen közelében lehet makroszkopikus méret˝u határfelületet modellezni.
4. Els˝oként vizsgáltam meg három dimenzióban a véges rácsfelbontás hatását a kialakult fázisegyensúlyi görbére a rács-Boltzmann módszeren belül egy kocka alakú, szemközti oldalain periodikusan csatolt térrészben, [9] az LBM3D2P - kétfázisú programkóddal. A szimulációkhoz a Shan és Chen által javasolt módszert alkalmaztam. a.
A szimulációs eredmények három dimenzióban is igazolták - az egy és a két dimenzióban már megfigyelt - véges felbontás hatását a kialakult fázisegyensúlyi görbére.
b.
Azt találtam, hogy a kialakult fázisegyensúlyi görbe már 40x40x40-es rácsfelbontás esetében is jól követi az állapotegyenletnek megfelel˝o értékeket.
c.
A háromdimenziós szimulációk során kapott fázisegyensúlyi görbe igen jól közelíti a kétdimenziós vizsgálatok fázisegyensúlyi görbéjét, tehát megállapítható, hogy abban az esetben, amikor csak a fázisegyensúlyi görbét szeretnénk meghatározni, a nagy számítási igény˝u, háromdimenziós számítások helyett elegend˝o a jóval kevesebb számításokat igényl˝o kétdimenziós számítások elvégzése is.
A tézispontokhoz kapcsolódó tudományos közlemények • [1] P. Kávrán, G. Mayer, G. Házi, Nuclear Energy for New Europe 2003, Portoroz, Slovenia, September, 2003 • [2] G. Mayer, G. Házi, Direct numerical simulation of longitudinal flow along triangular array of rods using the lattice Boltzmann method, The 14th International Conference on Discrete Simulation of Fluid Dynamics in Complex Systems, Kyoto, Japan, August 22-26, 2005 • [3] G. Házi, G. Mayer, Flow in rod bundles, Nuclear Energy for New Europe 2005, Bled, Slovenia, September 5-8. • [4] G. Mayer, G. Házi, Direct numerical and large eddy simulation of longitudinal flow along triangular array of rods using the Lattice Boltzmann method, Mathematics and Computers in Simulation, Vol.72, 173-178, 2006 • [5] G. Mayer, J. Páles, G. Házi, Large eddy simulation of subchannels using the lattice Boltzmann method, Annals of Nuclear Energy, 34, 140-149, 2007 • [6] G. Mayer, Large eddy simulation of a fuel rod subchannel, ICONE15, 15th International Conference on Nuclear Engineering, April 22-26, 2007, Nagoya, Japan • [7] G. Mayer, G. Házi, J. Páles, A. R. Imre, B. Fischer and T. Kraska, On the system size of lattice Boltzmann simulations, International Journal of Modern Physics C., Vol.15, No. 8, 2004 • [8] G. Mayer, G. Házi, A. R. Imre, B. Fischer, J. Páles and T. Kraska, Lattice Boltzmann and molecular dynamics simulation: on the size of mesoscopic systems, 3rd International Workshop "Global Phase Diagrams", Sept. 14-19, Odessa, Ukraine, 2003 5
• [9] G. Mayer, G. Házi, A. R. Imre, T. Kraska and L. V. Yelash, Lattice Boltzmann simulation of vapour-liquid equilibrium on 3D finite lattice, International Journal of Modern Physics C., Vol. 15, No. 3, 2004
További tudományos közlemények • [10] G. Házi, G. Mayer, I. Farkas, P. Makovi and A. A. El-Kafas, Simulation of a small loss of coolant accident by using RETINA V1.0D code, Annals of Nuclear Energy, Vol.28, 1583-1594, 2001 • [11] G. Házi, G. Mayer and I. Farkas, RETINA V1.0D Kétfázisú áramlásokat modellez˝o programrendszer, Magyar Energetika 1, 2001 • [12] I. Farkas, G. Házi, G. Mayer, A. Keresztúri, Gy. Hegyi and I. Panka, First experience with a six-loop nodalization of a VVER-440 using a new coupled neutronic-thermohydraulics system KIKO3D-RETINA V1.1D, Annals of Nuclear Energy, 29, 2235-2242, 2002 • [13] G. Házi, I. Farkas and G. Mayer, RETINA: A new simulator code for two-phase flow modeling, Int. Conf. on Supercomputing in Nuclear Applications SNA2000, Toranamon-Pastoral, Tokyo, Sept 4-7, 2000 • [14] E. Végh, G. Házi, J. S. Jánosy, G. Mayer, L. Seregi, B. K. Szabó, Graphical Simulator of the ETTR-1 Research Reactor, 15th European Simulation Multiconference, Jun. 6-9, Prague, Czezh Rep., 2001 • [15] G. Házi, I. Farkas, G. Mayer, Kétfázisú áramlások modellezésének fizikai és numerikus kérdései, KFKI AEKI Tanulmány, OMFB Részjelentés, 1999 ALK-00093/98/1, 1999 • [16] G. Házi, I. Farkas, G. Mayer, P. Makovi, RETINA: Two-phase Thermohydraulics for NPP Simulators (1999), AEKI Progress Report, pp. 23 • [17] G. Házi, I. Farkas, G. Mayer, RETINA: Termohidraulikai és szoftver rendszerterv, ALK 00093/98/2, 1999 • [18] G. Házi, G. Mayer, RETINA: Benchmark feladatok el˝ozetes kidolgozása, OMFB részjelentés, ALK-00093/98/3, 1999 • [19] G. Házi, I. Farkas, G. Mayer, RETINA: A fizikai és numerikus tesztek eredményei, OMFB Zárójelentés 2000, ALK-00093/98 • [20] A. R. Imre, G. Mayer, G. Házi, F. Römer, R. Rozas, T. Kraska, Determination of the spinodals of pure fluids from interfacial properties: molecular dynamics an Lattice Boltzmann Simulations, European Symposium on Applied Thermodynamics, Cannes, May/June 2008 • [21] G. Mayer, Coherent Structures in rod bundle flows: Lattice Boltzmann simulation, The Fourth International Conference for Mesoscopic Methods in Engineering and Science, July 16-20, 2007, Munich, Germany • [22] A. R. Imre, G. Mayer, G. Házi, R. Rozas, T. Kraska, Estimation of the liquid-vapor spinodal from interfacial properties obtained from molecular dynamics and lattice Boltzmann simulations, The Journal of Chemical Physics, 128, 114708, 2008 • [23] T. Kraska, L. V. Yelash, A. R. Imre, G. Házi, G. Mayer, Fluid-fluid phase transitions in pure binary liquids, DFG-OTKA Workshop on Promoting Hungarian-German Cooperations in Physics, May 2003, Budapest 6
• [24] A. R. Imre, G. Mayer, G. Házi, R. Rozas, T. Kraska, Estimation of spinodals in pure fluids from interfacial properties obtained from equilibrium molecular dynamics and lattice Boltzmann simulations, DPG-Frühjahrstagung, Berlin, Februar 2008 • [25] A. R. Imre, G. Mayer, G. Házi, R. Rozas, T. Kraska, Abschätzung der Spinodalen von fluiden Reinstoffen aus Grenzflächneigenshaften: Molekulardynamische und Lattice Boltzmann Simulationen, Bunsentagung, Saarbrücken, May 2008 • [26] G. Házi, A. R. Imre, G. Mayer and I. Farkas, Lattice-Boltzmann methods for two-phase flow modeling, Annals of Nuclear Energy, 29, 1421-1453, 2002 • [27] G. Mayer, G. Házi and I. Farkas, Determination of pressure loss coefficient of a VVER440 steam generator with lattice Boltzmann simulation, Bridging the time-scale gap, Sept 10-13, Konstanz, Germany, 2001 • [28] G. Mayer, Falsúrlódásból származó nyomásesés vizsgálata egy VVER-440 típusú g˝ozfejleszt˝o köpenyterén keresztül, háromdimenziós rács-Boltzmann módszer segítségével, Nukleáris Technikai Szimpózium 2002, Budapest, 2002. október 3-4.
Hivatkozások [Imre és Házi, 2002] Imre, A. R. és Házi, G. (2002). The effect of finite lattice-size in lattice Boltzmann model. Int. J. Mod. Phys C, 13(5):649–657. [Qian és Chen, 1997] Qian, Y. H. és Chen, S. (1997). Finite size effect in lattice-BGK models. Int. J. Mod. Phys. C., 8(4):763–771. [Shan és Chen, 1993] Shan, X. és Chen, H. (1993). Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E, 47(3):1815. [Treeck és társai, 2007] Treeck, C., Hillion, F., Cordier, A., Pfaffinger, M., és Rank, E. (2007). Flow and thermal sensation analisys in an aeroplane passenger cabin using lattice Boltzmann based approach. International Conference for Mesoscopic Methods in Engineering and Science. [Trupp és Azad, 1975] Trupp, A. C. és Azad, R. S. (1975). The structure of turbulent flow in triangular array rod bundles. Nucl. Eng. Des., 32:47.
7
