82 | Pomurska obzorja 1/ 2014/ 2 Tehnika Matej Zadravec 1 , Uroš Jeke 2 , Matjaž Hriberšek 1 Numerična simulacija uparjanja kapljevitega filma 1. Uvod V mnogih inženirskih problemih predvsem pri procesih klimatizacije in prezračevanja in marsikaterem procesnem postrojenju se srečujemo s procesom uparjanja in kondenzacije. V vseh teh inženirskih problemih imamo večinoma prisoten vlažen zrak, ki pa vemo vsebuje pri nekem stanju temperature in tlaka določeno količino vlage. Tak vlažen zrak je sposoben sprejeti še dodatno količino vlage vse do točke nasičenja vlažnega zraka, pri kateri se začne vlaga izločati iz vlažnega zraka v obliki kapljevite faze in sicer v obliki kapljic. Točka nasičenja je odvisna od temperature in tlaka pri kateri se nahaja vlažen zrak in zato zasledimo plast kondenzirane vlage na hladnejših površinah, kjer je na teh površinah njihova temperatura nižja od temperature rosišča pri določenem tlaku. Vsa ta stanja je zelo dobro možno opisati s pomočjo analitičnih in empiričnih enačb za enostavne primere, ko imamo znane hitrosti in ostale fizikalne parametre sistema, pri čemer pa je uporaba teh analitično-empiričnih nastavkov v bolj geometrijsko in fizikalno kompleksnejših geometrijah nemogoča zaradi različnih vplivov na sam proces uparjanja in kondenzacije. Dandanes obstaja precej računalniških paketov, ki so sposobni reševati kompleksne inženirske problema s področja dinamike tekočin. Ampak vseh fizikalnih pojavov seveda nimajo vključenih in zajetih, kar pa je seveda možno v same programske pakete vgraditi preko zunanjih računalniških podprogramov, ki jih lahko uporabnik sprogramira in implementira v obstoječi računski paket. V našem primeru reševanja prej omenjenih fizikalnih pojavov smo se odločili za začetni korak implementacije procesa uparjanja preko različnih empiričnih nastavkov izračuna stanj in točke nasičenja vlažnega zraka. Kot testni primer je bil izbran preprost primer izhlapevanja vode iz plasti kapljevine v kanalu, skozi katerega teče zrak. Primer je dobro obdelan eksprimentalno kot tudi numerično v članku Talukdarja [1]. Predstavljen prispevek se navezuje na na magistrsko delo Jekeja [2]. 2. Definicija problema Vključeni teoretični model uparjanja je preizkušen na modelu kratkega 3D pravokotnega kanala prikazanega na sliki 1, katerega geometrija je bila povzeta po izhodiščnem delu [1]. Simulacija je bila izvedena za vrednost Reynoldsovega števila Re=844. Slika 1. – Geometrija numeričnega modela [2]. POVZETEK V večini tehniških in naravnih sistemov je prisoten vlažen zrak. Vlažen zrak lahko pri določenem tlaku in temperaturi sprejme določeno količino vlage, pri čemer pride do točke nasičenja. Kondenzacija na stenah trdnih površin velikokrat ovira normalno delovanje v tehničnih napravah in je nezaželena. Nasproten pojav kondenzacije je uparjanje. Predstavljen je splošen matematično- fizikalni model uparjanja kapljevitega filma in njegova implementacija v obstoječi programski paket za računalniško dinamiko tekočin. Implementacija numeričnega modela uparjanja je bila preverjena na primeru simulacije pravokotnega 3D kanala in primerjana z izvedenim eksperimentom uparjanja kapljevitega filma. Zrak, ki teče nad vodno površino uparja vodo, pri čemer se vodna površina ohlaja, dokler se ne vzpostavi ravnotežje med prevodom toplote iz zraka in latentne toplote, ki se porabi za uparjanje. Spremljali smo kombiniran prenos toplote in snovi. Ugotovljeno je bilo, da je vpliv naravne konvekcije v primerjavi s konvektivnim tokom zanemarljiv in da je bila implementacija modela uparjanja uspešna. Ključne besede: uparjane, računalniška dinamika tekočin, naravna konvekcija, mešana konvekcija. 1 Fakulteta za strojništvo, Univerza v Mariboru 2 Hella Saturnus Slovenija d.o.o. Matej ZADRAVEC, Uroš JEKE, Matjaž HRIBERŠEK: NUMERIČNA SIMULACIJA… Pomurska obzorja 1/ 2014/ 2 | 83 Zaradi zahtevnosti obravnavane problematike bomo modelirali enofazni dvokomponentni tok, pri čemer se bo druga komponenta zaradi majhnih masnih deležev obravnavala ločeno od zraka in sicer s posebno enačbo ohranitve vodne pare v obliki koncentracije vodne pare. Sistem vodilnih enačb, ki jih moramo rešiti, da dobimo podatke o lokalnem koncentracijskem polju vodne pare kot tudi o difuzijskih tokovih na robu območja, je razširjeni sistem Navier-stokes enačb (en.1-6): Upoštevali bomo zgolj časovno ustaljeno stanje, seveda pa je model možno razširiti za simulacijo časovno spremenljivih razmer. Vpliv prisotnosti vodne pare v zraku je upoštevan z dodatnimi vzgonskimi silami v enačbi ohranitve gibalne količine (en.3), temperaturni razteznostni koeficient 𝛽 𝑇 in snovni razteznostni koeficient 𝛽 𝐶 podamo z sledečima enačbama (en.7 in en.8) kjer sta 𝑀 𝑧 in 𝑀 𝑣 molski masi zraka in vode. Pri termičnem razteznostnem koeficientu upoštevamo temperaturo, pri kateri bodo upoštevane lastnosti zraka v numeričnem modelu, in sicer 25°𝐶 , kar pomeni, da je bila vrednost gostote 𝜌 0 = 1,185 𝑘𝑔 / 𝑚 3 . Za sistem zrak-voda smo uporabili vrednost snovne difuzivnosti 𝐷 𝑣 = 2,82 ∙ 10 −5 𝑚 2 /𝑠 [3]. Snovne lastnosti zraka so bile konstantne. Uporabljene lastnosti zraka so bile pri temperaturi 𝑇 = 25 °𝐶 in referenčnem tlaku 𝑝 = 1 𝑏𝑎𝑟 . Vpliv vlage na gibanje zraka bo upoštevan s pomočjo dodatne konvektivne hitrosti na meji med vodno površino in zrakom, ki je posledica procesa izhlapevanja vode. Predpostavili bomo, da je glavni mehanizem prenosa snovi v tanki mejni plasti tik ob vodni površini difuzija, zato lahko upoštevamo model enostranske difuzije – Stefanov tok, [4]. Uparjanje bomo predpisali kot tok z neko hitrostjo glede na normalo na površino. Upoštevan bo tudi prenos toplote med zrakom in kapljevinskim filmom. Na vodni gladini bomo spremljali snovni uparjalni tok, ter lokalne in povprečne vrednosti Nusselt-ovega ter Sherwoodovega brezdimenzijskega števila (en.9 in en.10): 𝑁𝑢 𝑧 = 𝛼 ∙𝐷 ℎ 𝜆 = − 𝜕𝑇 𝜕 𝑦 | 𝑦 =−0,01025 𝑇 𝑤 − 𝑇 𝑚 𝐷 ℎ , (9) 𝑆 ℎ 𝑧 = 𝐾 ∙𝐷 ℎ 𝐷 = − 𝜕𝐶 𝜕𝑦 | 𝑦 =−0,01025 𝐶 𝑤 − 𝐶 𝑚 𝐷 ℎ . (10) kjer se indeks w nanaša na temperaturo oz. koncentracijo na vodni površini in indeks m na povprečno temperaturo oz. koncentracijo po celotnem volumnu kanala. Položaj y=-0,01025 je izbran tako, da sovpada s položajem površine vodne gladine. Na površini vodne gladine je za izračun parcialnega tlaka nasičenja vodne pare uporabljena Antoine-jeva enačba (en.11), ki velja za tlak 𝑝 = 1 𝑏𝑎𝑟 in temperaturno območje 274 𝐾 < 𝑇 < 373 𝐾 : 𝑝 𝑣 ,𝑠 = 𝑒𝑥𝑝 (11,96481− 3984,923 𝑇 −39,724 ) [𝑏𝑎𝑟 ], (11) Koncentracijo vodne pare v zraku (en.12) naposled izračunamo s pomočjo definicije masnega deleža (en.13) in definicije vlažnosti zraka (en.14) 𝐶 = 𝑚 𝑣 𝑉 = 𝜉 𝑣 𝜌 , (12) 𝜉 𝑣 = 𝑚 𝑣 𝑚 = 𝑚 𝑣 𝑚 𝑣 +𝑚 𝑧 = 𝑥 1+𝑥 , (13) 𝑥 = 0,622∙ 𝑝 𝑣 ,𝑠 𝑝 −𝑝 𝑣 ,𝑠 = 𝑥 𝑠 (𝑝 , 𝑇 ). (14) 3. Numerični model Numerični model kanala smo diskretizirali s strukturirano heksaedersko mrežo, ki je imela 105 .896 vozlišč kar nam je omogočalo dovolj natančne rezultate ob sprejemljivih računskih časih. Največji izziv pri implementaciji modela uparjanja je bilo direktno predpisovanje robnih pogojev odvoda spremenljivke, bodisi po komponenti kartezijevega koordinatnega sistema bodisi po normali na površino. Rešitev upoštevanja takih robnih pogojev na vodni površini smo našli v implementaciji podprograma napisanega v programskem jeziku Fortran 90, ki se izračuna izven osnovnega računskega modula reševanja sistemov enačb in se potem podatki preko vmesnika prenašajo iz osnovnega računskega algoritma v ta podprogram. Osnovni koncept je prikazan na sliki (2). 𝜕𝑢 𝜕𝑥 + 𝜕𝑣 𝜕𝑦 + 𝜕𝑤 𝜕𝑧 = 0, (1) [𝑢 𝜕𝑢 𝜕𝑥 + 𝑣 𝜕𝑢 𝜕𝑦 + 𝑤 𝜕𝑢 𝜕𝑧 ] = = − 1 𝜌 0 𝜕𝑝 𝜕𝑥 + 𝜈 0 [ 𝜕 2 𝑢 𝜕 𝑥 2 + 𝜕 2 𝑢 𝜕 𝑦 2 + 𝜕 2 𝑢 𝜕 𝑧 2 ], (2) [𝑢 𝜕𝑣 𝜕𝑥 + 𝑣 𝜕𝑣 𝜕𝑦 + 𝑤 𝜕𝑣 𝜕𝑧 ] = − 1 𝜌 0 𝜕𝑝 𝜕𝑦 + 𝜈 0 [ 𝜕 2 𝑣 𝜕 𝑥 2 + 𝜕 2 𝑣 𝜕 𝑦 2 + 𝜕 2 𝑣 𝜕 𝑧 2 ] + + 𝛽 𝑇 𝑔 (𝑇 − 𝑇 0 ) + 𝛽 𝐶 𝑔 (𝐶 − 𝐶 0 ), (3) [𝑢 𝜕𝑤 𝜕𝑥 + 𝑣 𝜕𝑤 𝜕𝑦 + 𝑤 𝜕𝑤 𝜕𝑧 ] = − 1 𝜌 0 𝜕𝑝 𝜕𝑧 + 𝜈 0 [ 𝜕 2 𝑤 𝜕 𝑥 2 + 𝜕 2 𝑤 𝜕 𝑦 2 + 𝜕 2 𝑤 𝜕 𝑧 2 ], (4) [𝑢 𝜕𝑇 𝜕𝑥 + 𝑣 𝜕𝑇 𝜕𝑦 + 𝑤 𝜕𝑇 𝜕𝑧 ] = 𝜆 0 [ 𝜕 2 𝑇 𝜕 𝑥 2 + 𝜕 2 𝑇 𝜕 𝑦 2 + 𝜕 2 𝑇 𝜕 𝑧 2 ], (5) [𝑢 𝜕𝐶 𝜕𝑥 + 𝑣 𝜕𝐶 𝜕𝑦 + 𝑤 𝜕𝐶 𝜕𝑧 ] = 𝐷 𝑣 [ 𝜕 2 𝐶 𝜕 𝑥 2 + 𝜕 2 𝐶 𝜕 𝑦 2 + 𝜕 2 𝐶 𝜕 𝑧 2 ] + 𝑟 . (6) 𝛽 𝑇 = − 1 𝜌 𝜕𝑝 𝜕𝑇 = 1 𝑇 0 = 1 298,15 𝐾 = 0,00335 𝐾 −1 , (7) 𝛽 𝐶 = − 1 𝜌 𝜕𝑝 𝜕𝐶 = 1 𝜌 0 [ 𝑀 𝑣 𝑀 𝑧 − 1] = = 1∙𝑚 3 1,185 𝑘𝑔 [ 28,960 𝑘𝑔 /𝑘𝑚𝑜𝑙 18,016 𝑘𝑔 /𝑘𝑚𝑜𝑙 − 1] = 0,513 𝑚 3 /𝑘𝑔 , (8) Matej ZADRAVEC, Uroš JEKE, Matjaž HRIBERŠEK: NUMERIČNA SIMULACIJA… 84 | Pomurska obzorja 1/ 2014/ 2 Slika 2. Numerični algoritem reševanja s pomočjo zunanje rutine v programski paket Ansys CFX 13.0 Robni pogoji ki smo jih predpisali na omenjenem modelu se lahko delijo na štiri osnovne sklope:  Za vtok zraka je bil predpisan analitični hitrostni profil zraka v pravokotnem kanalu [6]. Povprečna hitrost 𝑤 𝑎𝑣 je definirana preko Reynoldsovega števila na vstopnem robu, kjer smo upoštevali 𝑅𝑒 = 844, kar nam okarakterizira laminaren tok vlažnega zraka [3]. Temperatura na vstopu je bila izbrana, glede na eksperimentalne in numerične rezultate v izhodiščnem delu [1] in znaša 𝑇 0 = 22,4 °𝐶 . Na vstopnem robu predpišemo še vrednost koncentracije vodne pare. Ta se izračunava kot t.i. dodatna spremenljivka. Vrednost koncentracije vlage v vlažnem zraku je odvisna od temperature na vstopu 𝑇 0 in relativne vlažnosti zraka 𝜑 0 = 53,1 % [1].  Na stranskih stenah in zgornji steni kanala smo predpisali brez zdrsni robni pogoj, kar pomeni da so vrednosti hitrosti zraka enake nič. Na teh stenah je bil predpisan adiabatni robni pogoj. Preko sten ni nobenega toka dodatne spremenljivke za vlago in je odvod koncentracije po normali enak nič.  Vodna površina je definirana kot stena na kateri imamo izvor toplote in dodatno hitrost zraka v domeno s hitrostjo, ki je enaka hitrosti uparjanja Temperatura na vodni površini bo nižja kot pa je temperatura toka zraka. Namesto konstantnih vrednosti temperature in koncentracij na vodni površini se bosta temperatura in koncentracija izračunavali na podlagi robnega pogoja za energijo, ki pravi, da je latentna toplota uparjanja enaka prenosu toplote med zrakom in vodo. Uparjanje vlage je v tem modelu upoštevano kot tok s hitrostjo 𝑣 𝑤 . Hitrost je glede na normalo na površino enaka 𝑣 𝑤 = ( 𝐷 𝑣 1−𝜉 𝑣 ( 𝜕 𝜉 𝑣 𝜕𝑛 ) ) 𝑠 . Pri pojavu uparjanja in kondenzacije imamo zaradi spremembe agregatnega stanja intenziven prenos toplote. Ta toplotni tok vpliva na ravnotežno stanje na vodni gladini. Toplotni tok, ki priteka na vodno površino je vsota prevoda toplote v tekočini tik ob steni in latentne toplote, ki se pri uparjanju porabi, pri kondenzaciji pa sprošča. Predpostavimo, da imamo na vodni površini adiabatno steno, torej je vrednost toplotnega toka enaka 𝑄 ̇ 𝑣 𝐴 = 0. Ob tej predpostavki je na našem robnem pogoju količina toplote potrebna za fazno spremembo pri uparjanju vode enaka senzibilnemu prenosu toplote iz zraka na vodno površino. Robni pogoj za zrak na vodni površini lahko zapišemo kot −𝜆 𝜕𝑇 𝜕𝑛 = 𝑟 0 ∙ 𝐷 𝑣 ∙ 𝜕 𝐶 𝜕𝑛 . Vrednost koncentracije vodne pare tik ob vodni površini temelji na temperaturi vode 𝑇 𝑤 . Temperatura se izračunava iz enačbe −𝜆 𝜕𝑇 𝜕𝑛 = 𝑟 0 ∙ 𝐷 𝑣 ∙ 𝜕 𝐶 𝜕𝑛 in predpostavlja razmere popolnoma nasičenega zraka. Vrednost koncentracije vlage izračunamo iz enačbe 𝐶 = 0,622∙ 𝑝 𝑣 ,𝑠 𝑝 −𝑝 𝑣 ,𝑠 ∙ 𝜌 0 , kjer dobimo vrednost 𝑝 𝑣 ,𝑠 iz Antoine-ove enačbe (en.11). Na spodnji površini ne smemo pozabiti, da je širina vodne površine ožja kot je širina kanala. Na preostankih spodnjega roba, ki ni del vodne površine, predpostavljamo adiabaten robni pogoj, kjer so stene nepropustne za toplotni tok in snovni tok vodne pare.  Na izstopnem robu nimamo nobenim sprememb gradientov, ki bi bili gonilna sila in predpišemo 𝜕𝑢 𝜕𝑥 = 𝜕𝑣 𝜕𝑦 = 𝜕𝑤 𝜕𝑧 = 𝜕𝑇 𝜕𝑧 = 𝜕𝐶 𝜕𝑧 = 0. 4. Rezultati Na vstopnem robu imamo predpisani hitrostni profil kar lahko vidimo na sliki 3 in razvidno je tudi kako se tok zraka v kanalu še razvija in se v 0,6 𝑚 dolžine popolnoma še ne razvije. Slika 3. Hitrost v prerezni ravnini YZ pri x=0,00 m in detajl hitrosti na vstopnem robu. Matej ZADRAVEC, Uroš JEKE, Matjaž HRIBERŠEK: NUMERIČNA SIMULACIJA… Pomurska obzorja 1/ 2014/ 2 | 85 Na vodni površini smo upoštevali mehanizem enostranske difuzije in s tem gibanje vodne pare (uparjanje vode) opisali kot tok s hitrostjo podano na sliki 4 (levo). Vrednosti hitrosti zaradi uparjanja vidimo da so zelo nizke v primerjavi s konvektivno hitrostjo, ki smo jo predpisali toku zraka na vstopu. Za potrditev rezultatov smo izrisali še diagram povprečnih vrednosti hitrosti na vodni površini v prečni smeri na tok (smer x) po dolžini kanala v primeru z upoštevanjem difuzijske hitrosti na medfazni meji in brez upoštevanja. (slika 4 - desno). Če opazujemo krivuljo difuzijske hitrosti ima ta eksponentno obliko. Proti koncu kanala se vrednost približuje konstantni vrednosti, vendar se še ne ustali. To lahko povežemo z razvojem mejne plasti, ki se sodeč po difuzijski hitrosti popolnoma še ne razvije. Vrednosti difuzijske hitrosti so v primerjavi z ostalimi hitrostmi v kanalu nekaj redov manjše, zato je pričakovan zanemarljiv vpliv. Na sliki 5 so prikazane konture temperatur in koncentracijske porazdelitve (slika 6). Iz slik je vidno, da se termična in koncentracijska mejna plast še razvijata. Iz literature je moč zaslediti, da se termična mejna plast popolnoma razvije pri toku z 𝐷 ℎ = 0,0384 𝑚 , 𝑅𝑒 = 844 in 𝑃𝑟 = 0,71 šele po več kot 1 𝑚 . Iz dobljenih rezultatov lahko ugotovimo, da se kljub nižji vrednosti karakterističnega števila koncentracijski profil še ne razvije popolnoma po dolžini kanala 𝐿 = 0,6 𝑚 . Rezultati Nusseltovega in Sherwoodovega števila na sliki 7 so podani za primer brez upoštevanja difuzijske hitrosti na vodni površini in z upoštevanjem difuzijske hitrosti ter so primerjani z vrednostmi, ki jih je v svoji simulaciji pridobil Talukdar [1]. Vrednosti so povprečne vrednosti prečno na smer toka (v 𝑥 smeri) na ravnini 𝑋𝑍 , kjer je 𝑦 = −0,01025 𝑚 . Iz teh rezultatov lahko ugotovimo, da upoštevanje difuzijske hitrosti na vodni površini ne vpliva na končni rezultat. Dobljene vrednosti na vstopnem delu kanala odstopajo od Talukdarjevih rezultatov [1], ko se pomikamo vzdolž kanala (po 𝑧 smeri) se rezultati približujejo in vse bolj ujemajo. Vzrok v delnem razhajanju je v primerjalno redki računski mreži, katere vpliv je največji prav na področju najvišjih snovnih tokov, to je na vstopu v kanal. Slika 4. Hitrosti na vodni površini. Slika 5. Izoterme na ravnini YZ pri x=0,00 m. Slika 6. Konture porazdelitve koncentracij vodne pare na ravnini YZ pri x=0,00 m. Matej ZADRAVEC, Uroš JEKE, Matjaž HRIBERŠEK: NUMERIČNA SIMULACIJA… 86 | Pomurska obzorja 1/ 2014/ 2 Slika 7. Razporeditev lokalnih vrednosti Nusseltovega in Sherwoodovega števila vzdolž kanala. 5. Zaključek Izdelani računalniški podprogram, ki je bil uporabljen v programskem paketu Ansys-CFX [5], je primeren za modeliranje procesa uparjanja na vodnih površinah. Rezultati izvedenih numeričnih simulacij za izbrani testni primer dokazujejo, da se rezultati dobro ujemajo z referenčnimi rezultati [1]. Pri tem je potrebno poudariti, da je za uporabo modela na gostejših računskih mrežah, kjer je vpliv izmenjave toplote tik ob vodni površini lokalno zelo velik, potrebno podrobneje preučiti način vključitve te vrste prenosa toplote v računski algoritem Ansys-CFX, predvsem s ciljem, da zagotovimo monotono konvergenco celotnega sklopa enačb. Literatura 1. Talukdar Prabal, Iskra Conrad R., Simonson Carey J. Combined heat and mass transfer for laminar flow of moist air in a 3D rectangular duct: CFD simulation and validation with experimental data. International Journal of Heat and Mass Transfer, 2008 , 51 , 3091 – 3102. 2. U. Jeke. Modeliranje kondenzacije in uparjanja vlage na trdnih površinah z računalniško dinamiko tekočin. Magistrsko delo študijskega programa 2. stopnje Strojništvo. Fakulteta za strojništvo, Univerza v Mariboru. Avgust 2012. 3. Škerget Leopold. Prenosni pojavi: Prenos gibalne količine, toplote in snovi. Osnutek učbenika., 2012. 4. L. Škerget, M. Hriberšek, J. Ravnik, M. Ramšak. Modeliranje kondenzacije in uparjanja vlage na trdnih površinah svetil z računalniško dinamiko tekočin. Poročilo, Univerza v Mariboru, Fakulteta za strojništvo, Inštitut za energetsko, procesno in okoljsko inženirstvo, 2010. 5. Ansys inc. Ansys CFX 13.0, 2010. 6. Shah R.K., London A. Chapter VII; Rectangular Ducts. V: Shah R.K., London A. Advances in Heat Transfer: L. Laminar Flow Forced Convection in Ducts. New York: Academic Press, 1978. ,000 2,000 4,000 6,000 8,000 10,000 12,000 14,000 16,000 18,000 20,000 00,001 00,010 00,100 Nuz, Shz z/(Re Pr Dh) NuZ Talukdar [13] ShZ Talukdar [13] NuZ brez dodane hitrosti ShZ brez dodane hitrosti NuZ dodana difuzijska hitrost