Univerzitet u Beogradu Masinski fakultet Мilada L. Pezo -l!Jt~"J;i'l~~J.r.в :899 ~ '4 О!;))'). ~с v v NUMERICКA SIМULACIJA КRIZE КLJUCANJA U ISP ARIV АСКIМ CEVIМA doktorska disertacija Beograd, 2011. Komisija za pregled i odbranu: Mentor: Prof. dr V1adimir Stevanovic, Masinski fakultet и Beogradu Clanovi komisije: Prof. dr Branislav Savic, Masinski fakuitet и Beogradu Prof. dr Dragan Tucakovic, Ma$inski fakultet и Beogradu Prof. dr Мilos Banjac, Ma$inski fakultet и Beogradu Dr Мilan Rajkovic, ~" ... ' . ' Institut za nuklearne nauke VINCA, Laboratorija za termotehniku i energetiku Datum odbrane: Predgovor Ovaj rad је nastao posle visegodisnjeg istraiivanja и Laboratoriji za termotehniku i energetiku, Instituta za nukleame nauke VINCA i na Ma5inskom fakultetu Univerziteta и Beogradu pod rukovodstvom Prof. dr Vladimira Stevanovica. Ovom prilikom se zahvaljujem Prof. dr Vladimiru Stevanovicu koji је prihvatio da rukovodi izradom ovog rada, kao i za korisne primedbe, pomoc i sugestije prilikom izrade i uoЫicavanja ovog rada. Beograd, 2011. Мilada L. Pezo Sadrzaj Abstakt i Abstract ii Nomenklatura iii 1. Uvod 1 2. Pregled dosadasnjih rezultata iz oЫasti krize kljucanja 5 3. Standardni modeli krize razmene toplote pri bazenskom kljucanju 24 3.1. NestaЬilnost parnog sloja 24 3.2. Nestabilnost pamog mlaza 25 3.3. Racunanje kriticnog toplotnog fluksa 26 3.4. Dvofazno strujanje u horizontalnom isparivacu 27 4. Razvoj modela krize kljucanja 29 4.1. Opis modela 30 4.1.1. Parametri mikro nivoa 31 4.1.2. Jednacine odrzanja 35 4 .1.3. Konstitutivne korelacije 3 7 4.1.4. Prelaz toplote izmedu zida i dvofazne mesavine 39 4.1.5. Uticaj zagrejaca na kritican toplotni fluks 40 4.1.6. Granicni uslovi 41 4.2. Strujanje и vertikalnoj cevi eetvrtastog poprecnog preseka 42 4.2.1. Opis modela 42 4.3. Strujanje u vertikalnoj cevi okruglog poprecnog preseka 43 4.3.1. Opis modela 43 4.3.2. Jednacine odrZanja mase, kolicine kretanja i energije za slucaj polarno-cilindricnih koordinata 43 5. Numericko resenje sistema diferencijalnih jednacina 46 5.1. Polarno-cilindrisni koordinatni sistem 54 5.1.1. Koeficijenti u diskretizovanim jednacinama 54 5.2. Primenjene numericke mreze 56 6. Prikaz i anal.iza dobljenih rezultata simulacije krize razmene topJote 58 6.1. Bazensko kljucanje 59 6.1.1. Uticaj obradenosti zagrejacke povrsine na pojavu kriticnog toplotnog fluksa 73 6.1.2. Uticaj kontaktnog ugla 77 6.1.3. Racunanje koeficijenta prelaza toplote koriS6enjem razliCitih empirijskih izraza i poredenje sa sopstvenim modelom 79 6.2. Strujanje u vertik:alnoj cevi kvadratnog poprecnog preseka 83 6.3. Strujanje u cilindricnoj vertik:alnoj cevi 87 7. Zakljucak 91 Literatura 93 Biografija 96 Prilog 97 М. Pezo, V. Stevanovic, Numerical prediction of critical heat flux in pool boiling with the two-fluid model, /nternational Journal of Heat and Mass Transfer, 54 (2011) 3296-3303 NUМERICКA SIМULACIJA КRIZE КLJUCANJA U ISP AR1V АСКIМ CEVIМA -Abstrakt- Predmet ovog rada је numericka simulacija i analiza trodimenzionalnog dvofaznog strujanja i istraiivanje mehanizama krize kljucanja и isparivackim cevima. Кriza kljucanja је nepovoljna pojava. Nagli porast temperature zida cevi izaziva termomebanicka naprezanja, koja mogu dovesti do pojave prskotina i pucanja cevi pod pritiskom. Razvijen је matematicki model dvofaznog strujanja u isparivackoj cevi, koji se sastoji od jednacina odrZanja mase, kolicine kretanja i energije za оЬе faze. Sastoji se od sistema parcijalnih diferencijalnih jednacina koje su recavane za definisane granicne uslove. Bilo је neophodno uvesti i modeliranje i numericku simulaciju na mikro nivou, jer је Ыlо potrebno ispitati i proces pojave i rasta mehura, kao i pona.Sanje dvofazne mesavine na zagrejackoj povrsini na makroskopskom nivou. Rezultati numerickiЬ simulacija su uporedeni sa raspolozivim eksperimentalnim rezultatima. Юjucne reci : kljucanje, kritican toplotni fluks, kriza kljucanja, numericka simulacija Naucna oЫast : Ma.Sinstvo Uza naucna oЫast: Termoenergetika 1 NUМERICAL SIМULATION OF BOILING CRISES МECВANISM IN EV APORA TION PIPES -Abstract- The subject of this thesis is numerical simulation and analyses of three dimensional two-phase flow and mechanism of boiling crises in evaporation pipes. Boiling crises is characterized Ьу а dried out heat surface and сао bring physical destruction of the heater. А developed mathematical model of two phase flow in evaporation pipe consists of prescribed mass, momentum and energy conservation equations for liquid and vapor two-phase flow. It consists of set of partial differential equations which were solved for specific boundary conditions. Modeling of the micro scale level was necessary to take into account processes of the ЬuЬЫе rise and behavior of the two-phase mixture at the heating surface. Results of the numerical simulation are compared with similar availaЫe results of the experiments from the literature. Кеу \.Vords : boiling, critical heat flux, boiling crises, numerical simulation Scientific discipline: Mechanical engineering Scientific subdiscipline: Thermal power engineering 11 Nomenklatura а ь Ср Со Csf D Dь F g h h h12 Ја k l Lc п р Pr q qh qь t т и х temperaturska provodnost sirina zone nukleacije specificna toplota koeficijent medufaznog trenja koeficijent и Rohsenow-oj korelaciji precnik. precnik. mehura pri odvaj anju sila и jedinici zapremine gravitaciono ubrzanje entalpija koeficij ent prelaza toplote latentna toplota isparavanja Jakob-ov broj termicka provodnost karak:teristicna dimenzija za bazensko kljucanj e duiina vodene kapilare gustina centra nukleacije pritisak Prandtl-ov broj toplotni fluks toplotni fluks (zapreminski udeo) toplotni izvor па zidu, koji se javlja zbog rasta mehura vreme temperatura brzina koordinata Grcka slova: а г 0 р µ \) (Ј 't zapreminski udeo faze udeo promene faze kontak:tni ugao kvasenja gustina dinamicka viskoznost kinematska viskoznost povrsinski napon relaksaciono vreme promene faze lndcksi: с kondenzac~a е isparavanje g para k indeks faze (k=l ,2) l tecnost L uzgon Nomenklatura m m N/m3 m/s2 J/kg W/m2K J/kg W/mК. m m -2 m Ра W/m2 W/m3 W/m3 s к rn/s m kg/mЗ-s deg kg/m3 Pas m2s-1 N/m s 111 Nomenklatura р deo sat saturacija vм virtuala masa w zid 1 voda 2 para 21 medufazni zasicena tecnost " zasicena para lV 1. Uvod 1. Uvod Predmet istraiivanja и okviru doktorske disertacije је proucavanje mehanizama krize kljucanja. Do krize kljucanja dolazi pri visokim toplotnim fluksevima i intenzivnoj generaciji mehurova. Pri tim uslovima se foпnira sloj dvofazne mesavine sa visokim zapreminskih udelom pare koji sprecava prodiranje tecne faze do zagrejane povrsine. Zagrejacka povrsina se zasu5uje, а s obzirom daje para znatno losiji prenosnik toplote od tecnosti, znaeajno se smanjuje prelaienje toplote sa zagrejacke povrsine na fluid, Sto је praceno naglim porastom temperature zida. Кriza kljucanja је nepovoljna pojava. Nagli porast temperature zida cevi izaziva termomehanicka naprezanja, koja mogu dovesti do pojave naprslina i pucanja cevi pod pritiskom. Takode, pri izrazito visokim porastima temperature moze doCi do topljenja materijala i pregorevanja zida cevi. Pojava krize kljucanja se и literaturi naziva i odstupanje od mehurastog kljucanja ("Departure from Nucleate Boiling"-DNB) ili kriza razmene toplote prve vrste. Toplotni fluks pri kome dolazi do krize kljucanja se naziva kriticni toplotni fluks ("Critical Heat Flux"-CНF). Predvidanje uslova pri kojima moze doCi do krize razmene toplote је narocito znacajno za sigurnost razlicitih vrsta generatora pare. Кriza kljucanja и isparivackim cevima parnog kotla и lozistu dovodi do oStecenja cevi i curenja radnog fluida iz cevnog sistema. Кriza kljucanja и jezgru nuklearnog reaktora dovodi do topUenja nukleamog goriva i sirenja radiokativnosti u okviru nukleamog sistema za proizvodnju pare. Takode, kriza kljueanja је i ogranicenje povecanju koeficijenta prelaienja toplote kod povr5ina koje se hlade kljueanjem tecnosti, odnosno и toplotnim aparatima kod kojih se razmena toplote intenzivira kljucanjem prijemnika toplote. Na ovaj nacin kriza kljueanja је ogranicavajuCi faktor smanjenju povrsina za razmenu toplote, odnosno smanjenju ukupnih dimenzija toplotnih aparata. Potreba za predvidanjem uslova pri kojima dolazi do krize kljucanja uslovila је sprovodenje velikog broja eksperimenata za razlicite strujne uslove i geometrije strujnih kanala. Na osnovu merenja razvijene su brojne empirijske korelacije koje omogucuju odredivanje kriticnog toplotnog fluksa и zavisnosti od razlicitih strujnih i geometrijskih parametara i termofizickih karakteristika fluida. Broj ovih korelacija prelazi nekoliko stotina, sto ukazuje na nepostojanje opsteg pristupa modeliranju pojave koji Ы Ьiо primenjiv na uslove od znacaja za tehnicku praksu. S obzirom na ogranieeni opseg primene razvijenih korelacija i greske и odredivanju kriticnog toplotnog fluksa koje se mogu javiti pri njihovoj primeni, za odredivanje kriticnog toplotnog fluksa se koriste i eksperimentalne baze podataka, takozvane tablice za ocitavanje kriticnih toplotnih flukseva ("look-up" tabele). Razvijene su prvo pedesetih godina proslog veka od strane Akademije nauka и Ыvsem Sovjetskom Savezu, а zatim su prihva6ene, redigovane i dopunjene i od strane naucnih institucija SAD-a i drugih zemalja. lako је ovaj proЫem izucavan и prethodnom periodu, mnogi uticajni faktori i kornpleksnost strujnih i termickih procesa jos uvek nisu dovoljno ispitani. Granicni uslovi, koji se odnose па zid isparivackog strujnog kanala i polje dvofazne mesavine uticu na to daje defmisanje, resavanje, analiza i predvidanje fenomena krize kljueanja izuzetno slozeno. Savremena istraiivanja krize kljucanja, primenom modeme merne tehnike, pokazuju da na procese kljucanja i pojavu pregrevanja uticu Ьitno pojave па mikro nivou па povrsini grejanja, povezane sa pona5anjem dvofazne mesavine na~~ ...... .,. e /Jlt411 V'Y' ".t, ! ~\ ~1 / ) . .. ;.. "" .ЈЈ ... " ":-'" f "'•;Ј. 1..w\"~' 1. Uvod makro nivou. Na kompleksnost procesa na mikro nivou bitno utice stepen obradenosti zagrejacke povтsine, sposobnost kva5enja kljucale tecnosti i prisutne hemijske necistoce. Jedan od rezultata izucavanja pojave na mikro nivou је i saznanje da sa povecanjem hrapavosti zagrevane povr5ine pove6ava se k:riticni toplotni fluks, jer se sa povecanjem hrapavosti pove6ava gustina klijalista mehurova. U okviru teze је razvijen matematicki model kljuCale dvofazne mesavine. Matematicki model se sastoji od jednacina odrZanja mase, kolicine kretanja i energije pare i tecnosti. Medufazno dejstvo је uzeto и obzir primenom odgovarajucih konstitutivnih korelacija. Model је resen primenom odgovarajuceg algoritma za resavanje modela dvofaznog strujanja. Da Ьi se па osnovu ovog modela mogle ispitati i pojave koje se odnose na odredivanje kriticnih vrednosti toplotnog fluksa, pojave zasusenja cevi, kao ј pojave pregrevanja zagrejackih povrsina potrebno је matematickom modelu dodati i jednacine koje se odnose na mikro nivo, tj. па nivo mehura, koji се omoguCiti predvidanje pojave i rasta mehurova. Parametarski је analiziran uticaj obradenosti zida zagrejacke povrsine, toplotnog fluksa, geometrije isparivackih cevi i procesa razmene na medufaznim povrsinama na pojavu krize kljueanja. Na osnovu ove analize је moguce odredivanje kriticnih vrednosti toplotnog fluksa za razlicite uslove kljucanja, kao i predvidanje pojave mehanizma pregrevanja isparivackih cevi. Cilj ovog istraZivanja је visestruk: (а) razvoj fizickog i matematickog modela i odgovarajuceg numerickog postupaka resavanja za spregnutu simulaciju konvektivnog kljucanja i termickih procesa u zagrejackom zidu; (Ь) sprovodenje numerickih simulacija dvofaznog strujanja tecne i рате faze preko zagrejackih povrsina iz1ozenih visokim vrednostima toplotnog fluksa, kao i proracun prostomog nestacionarnog provodenja toplote 1 temperaturskog polja u zidu zagrejackog kanala; (с) razvoj algoritma za odredivanje pojave krize kljucanja na osnovu sracunatog nestacionarnog temperaturskog polja и zagrejackom zidu i strujno-termickih parametara dvofaznog toka u simuliranom strujnom kanalu. Кrajnji cilj је razvoj metode za pouzdano predvidanje uslova pojave krize kljucanja i vrednosti kriticnog toplotnog fluksa. S obzirom da се razvijeni modeli obuhvatiti Ьitne mehanisticke procese koji odreduju pojavu krize kljucanja na mikro i makro nivou (mikro nivo se odnosi na povrsinu zagrejackog zida i granicni sloj dvofazne mesavine, а makro nivo na globalne parametre dvofaznog strujanja и isparivackom kanalu), razvijeni numericki pristup је primenljiv и opstem slucaju za razlicite geometrijske, termicke i strujne uslove (u opsegu masenih flukseva dvofaznog strujanja od nekoliko stotina do nekoliko hiljada kglm2s i masenili udela pare koji odgovaraju mehurastom i penastom dvofaznom toku), kao i razlicite termofizicke karakteristike tecne i рате faze. Postavljeni cilj је od izuzetnog znacaja za sigumost i pouzdanost rada razlicitih vrsta generatora pare. 2 1. Uvod Predlozeni rnodel dvofaznog strujanja se sastoji od jednaCina konzervacije mase, kolicine kretanja i energije za svak:u od faza. U slucaju prostornog trodimenzionog strujanja formira se sistern od osam parcijalnih diferencijalnih jednacina eliptickog tipa (ро jedna skalarna jednacina konzervacije mase i energije i tri projekcije vektorske jednacine konzervacije kolicine kretanja za svaku fazu). S obzirom na nelinearnost proЫema, prisustvo konvektivnih i difuzionih clanova и jednacinama konzervacije i slozenost konstitutivnih korelacija, analiticko ispitivanje stabilnosti resenja је neostvarivo, tako da se polazi od hipoteze da za fizicki realno definisane pocetne i granicne uslove razvijeni sistem jednacina konvergira ka jedinstvenom partikularnorn resenju. Razmena Ьilansnih velicina na razdelnim povrsinama tecne i parne faze se odreduje prirnenom odgovarajucih polu-empirijskib konstitutivnih korelacija. Uslovi generacije pare na mikro nivou na zagrejackom zidu su obuhvaceni pomocu gustine mesta nukleacije mehurova i vremena rasta mehura na zidu. Gustina mesta nukleacije је vafan parametar u procesu modeliranja. Ima znacajan uticaj na dinamiku kljucanja i uslove pod kojima dolazi do krize kljucanja. Na osnovu dosadasnjih istraiivanja moze se zakljuciti da је gustina oukleacije odredena hrapavoscu povrsine na kojoj dolazi do kljucanja. U modelu se preko gustine nukleacije moze uzeti и obzir hrapavost povrsine. Mesta nukleacije mehurova se zadaju pomocu probaЬilistickog modela, tako da се ona Ьiti promenljiva u prostoru, a1i се ukupan broj mehurova na odredenoj povrsini odgovarati gustini nukleacije mehurova. Drugi parametar koji је od znaeaja za lokalne uslove generacije mehurova је vreme rasta mehura na zidu. Veci broj parametara utice na vreme rasta mehura, kao sto su toplotni fluks, ugao kva5enja itd. Njihov uticaj se parametarski analizira odgovarajucim modelom. Za numericko resavanje sistema bilansnih jednacina koriscen је metod konacnih zapremina. Smenom jednacina konzervacije kolicine kretanja и skalarnom oЫiku и jednacinu konzervacije mase, doЬija se jednacina za korekciju pritiska u dvofaznom toku, koja se resava primenom metode SIМPLE (Semi-Irnplicit Method for Pressure- Linked Equations). Numericke simulacije dvofaznog strujanja sa kljucanjem и isparivackoj cevi su sprovedene i za eksperimentalne uslove raspolozive u literaturi. Rezultati numerickih simulacija su uporedeni sa raspolozivim eksperimentalnim rezultatima. U poglavljima 2 i 3 је dat pregled dosada5njih rezultata iz oЫasti eksperimentalnih istraZivanja pojave kriticnog toplotnog fluksa i razvijenih pristupa modeliranja kriticnih pojava и procesu isparavanja. U poglavlju 4 је dat detaljan opis modela dvofaznog strujanja и slucaju bazenskog kljucanja, za slueaj strujanja u vertikalnoj posudi cetvrtastog poprecnog preseka i za slucaj vertikalnog strujanja и cilindricnoj cevi. Ovde su definisani i granicni uslovi za svaki slucaj pojedinacno. Nurnericki postupak resavanja slozenog modela strujanja, spregnutog sa karakteristikama na mikro nivou је dat и poglavlju 5. Rezultati numericke simulacije i analiza strujno-termickih parametara su dati и poglavlju 6. Analiziran је uticaj hrapavosti povrsine na generaciju mehurova i na 3 1. Uvod pojavu kriticnog toplotnog fluksa. Takode је proucavan i uticaj razlicite geometrije i razlicitih granicnih usiova na proces modeliranja. Pretpostavljano је da se kriticni toplotni fluks javlja u trenutku kada dode do naglog temperaturskog skoka u zagrejackom zidu. Znacajan uticaj na kriticne pojave ima i vrednost toplotnog fluksa kome је izlozen zid, tako da је analiziran i nj egov uticaj. Uradena је i numericka ana1iza koriscenja optimalnog izraza za koeficijent prelaza toplote, kao i poredenje rezultata doЬijenih simulacijom i drugih pred1ozenih korelacija. U zakljucku, poglavlje 7 је dat saiet prikaz ostvarenih rezultata i njihov znacaj sa stanovista istraiivanja efikasnosti, sigurnosti i pouzdanosti rada u razlicitim postrojenjima. Predlozeru su i planovi i smernice za dalje istra2ivanje и ovoj oЫasti. 4 2. Pregled dosada5njih rezultata iz oЫasti krize kljucanja 2. Pregled dosadasnjih rezultata iz oЫasti krize kljucanja Nukiyama [1] је jos 1934 godine odredio maksimalnu i miлimalnu vrednost toplotnog fluksa koji se prenese sa metalne zagrejacke povrsine na kljucalu vodu pri atmosferskom pritisku. On је za eksperimente koristio posudu cetvrtastog poprecnog preseka koja је Ьila napunjena vodom, slika 2. 1. Vodu је zagrevao pomocu Zica ( od platine, nikla ili legure nikla i hroma) koje su Ьile uronjene и tecnost. DoЬio је krivu zavisnosti toplotnog fluksa q od razlike temperatura izmedu zagrejacke povrsine i vode ЛТ, slika 2.2. Slika 2. 1 Aparatura koris6ena za Nukiyam-in eksperiment Koriscenaje relacija: q=h·ЛT (2.1) Nukiyama је pokazao da se koeficijent prelaza toplote ne menja monotono sa povecanjem ЛТ. Kada је isparavanj e Ыаgо, uzbrkanost parnih mehurova ima veCi efekat na prenos toplote, tako da је moguce povecanje h i q istovremeno sa povecanjem ЛТ Ukoliko је generacija pare intenzivna onda је veci deo metalne povr5ine prekriven mehurovima, tako da vise nema direktnog kontakta izmedu zagrejacke povrsine i vode. Tada se javlja negativan efekat i postoji samo prenos topote izmedu metalne povrsine i pare. Tako na krivi ЛТ-h, slika 2.2 vrednosti na ordinati se povecavaju sa povecanjem ЛТ do kriticne tacke, zatim se smanjuju sa daljim povecanjem ЛТ. Kako је q proizvod h i ЛТ, ono ne Ьi trebalo da se smanjuje sa smanjenjem h. Diferenciranjem gornje jednacine i izjednacavanjem sa nulom se doЬija maksimalna vrednost toplotnog fluksa q max ( na slici је obelezeno sa poloZзjem Ь). Sa daljim povecanjem ЛТ se smanjuje q i dostize minimalnu vrednost u tacki с. Deo krive Ьс је izuzetno nestabilan i tesko ga је ostvariti u praksi. Pokazano је da maksimalna vrednost toplotnog :tluksa ne zavisi od materijala od koga је napravljena zagrejacka zica. Vrednosti za q max mogu da budu pril icno visoke u odnosu па ЛТ, kojije (20-30)0С. q max je oko (125-210) W/cm2• 5 2. Pregled dosada5njih rezultata iz oЫasti krize kljucanja д! Slika 2.2 Nukiyam-ina kriva Ма i Рал [2] su numericki istraiivali uticaj faktora zagrejacke povrsine na kljucanja za visoke toplotne flukseve, koji se karakterisu postojanjem mikrosloja. Definisali su dva regiona za numericku simulaciju: strujanje vodeno termo-kapilarnim pojavama u tecnom sloju i kondukciju toplote и cvrstom zidu. Rezultati numericke simulacije, slika 2.3 pokazali su da su strujanje и mikrosloju i isparavanje na granici parna-tecna faza veoma efikasan mehanizam za obja5njenje visokog koeficijenta prelaza toplote kljueanja Ыizu vrednosti kriticnog toplotnog fluksa. Pokazano је da veoma tanak zid ili zid sa slabim provodenjem toplote ima znacajan uticaj na strujanje и tecnom sloju i rasporedu temperatura и zagrejackom zidu. Celije и kojima se javlja vrtlog su izazvane и makrosloju i energija se moie efikasnije predati sa zagrejackog zida preko vrtlo:Znib 6elija i isparavanja na granici tecna-pama faza. Ove celij e vrtloga i isparavanje na granici tecna-parna faza dovode do veoma efikasnog mehanizma prenosa toplote i visokog koeficijent prenosa toplote pri mehurastom kljucanju Ыizu kriticnog toplotnog fluksa. Numericki rezultati pokazuju da veoma tanak zid i/ili zid sa slaЫm provodenjem predstavlja veliki temperaturski otpor za prenos toplote и bocnom pravcu. Као posledica se javlja jos veca neuniformnost raspodele temperature и zidu. Тај efekat se guЬi kada debUina zida i/ili kada је termicka provodnost zida veca od neke kriticne vrednosti. Posto је kritican toplotni fluks ogranicavajuca okolnost za mehurasto kljucanje pri visokim vrednostima toplotnog fluksa, ocekivan је znacajan uticaj deЬljine zida i termicke provodnosti na kritican toplotni fluks za tanke zidove i/ili za zid sa slaЬim termickim provodenjem. 6 2. Pregled dosadasnjih rezultata iz oЫasti krize kliucanja (а) Tw - 124 .88 (С) (ЬЈ Т(С) • 11!<.8 1\8.0 11\S.2 107.• 101.8 Slika 2.3 Raspodela vektora brzina (а) i raspodela temperaturskog profila (Ь) и zagrejackom zidu i u Ыizini zida Celata i dr. [3, 4] su razvili metod za predvidanje kriticnog toplotnog tluksa za vertikalno strujanje u cevi cilindricnog poprecnog preseka, slika 2.4. Pretpostavili su da postoji rastojanje od zida na kome је temperatura fluida jed.naka temperaturi zasicenja. То rastojanje se naziva pregrejani sloj i predstavlja jedinu oЫast gde је moguce da se formira mehur. Zbog rasta i akumulacije mehura fonnira se veci mehur ili sloj pare. Кritican toplotni fluks se javlja kada parni mehur razЬije pregrejani sloj i dode и dodir sa zagrejackim zidom. parni merur tecna faza • 1 , • - pojava СНF pojavaCHF - ' • • • Slika 2.4 Sematski prikaz podhladenog kljueanja Ыizu uslova za kritican toplotni fluks Pretpostavlja se da је deЫjina parnog sloja jed.naka precniku mehura i da ne zavisi od toplotnog fluksa, vec od fizickih, teпnohidraulickih i geometrijskih parametara. Pregrejani sloj zavisi i od toplotnog fluksa i od fizickih, termohidraulickih i geometrij skih parametara. 2. Pregled dosadasniih rezultata iz oЫasti krize kljucanja Inoue i Lee [5] su eksperimentalno ispitivali uticaj kriticnog toplotnog fluksa pri malim vrednostima pritiska na karakteristike strujanja, slika 2.9. Eksperimentja raden za vertikalnu cev, а kao radni fluid је koriscena voda. Кritican toplotni fluks је meren za razlicite vrednosti zapreminskog udela teene faze za penasti rezim pri atmosferskom pritisku. Cev је izradena od providnog stakla kako Ьi se pojava mogla pratiti i vizuelno. Strujanje је skoro nepromenljivo tokom eksperimenta, kako Ьi se pratio uticaj vrste strujanja na dinamiku pojave kriticnog toplotnog fluksa .. Pojava је ist:raZivana u mehurastom, cepastom, penastom i maglenom toku. Pokazano је da pojava kriticnog toplotnog fluksa zavisi od rezima strujanja. Uticaj brzine strujanja na pojavu kriticnog toplotnog fluksa је izuzetno mali. Pokazano је dobro slaganje eksperimentalnih rezultata sa korelacijama Zubera i Kutateladsea, cije su jednacine koriscene i и ovom radu. U cepastom rezimu strujanja vrednosti za kritican toplotni fluks su minimalne zbog fluktuacija u strujanju. U penastom rezimu strujanja vrednost kriticnog toplotnog fluksa se povecava sa povecanjem zapreminskog udela parne faze i smanjuje usled smanjenja brzine strujanja. Povecanje vrednosti kriticnog toplotnog fluksa se javlja zbog suzbUanja fluktuacija i povecanja zapreminskog udela tecne faze i brzine strujanja. U delu prelaznog rezima izmedu penastog i maglenog vrednost za kriticni toplotni fluks pokazuje nagli skok (pik). Pojava pika ne zavisi od brzine strujanja. U maglenom rezimu strujanja vrednost kriticnog toplotnog fluksa se smanjuje sa poveeanjem zapreminskog udela рате faze. о 0.1 0.1 " . v..:.1 о 1n . ~ . .... А •11 а llЖ о ' '11r:н~"") • Р.,.-О IMP• D"0.0%1 m !'\ • • • • • () • D • о .~ w, tk1>'m't ол 0.6, 0-'S 114 0.4, 03) 07, 06$ о'' 21 1 0 4, о,, 0.1 0.2 0.3 0.4 o.s х. Slika 2.5 Eksperimentalni rezultati-veza izmedu kriticnog toplotnog fluksa i zapreminskog udela parne faze za penasti rezim strujanja (levo) i lokacija na kojoj se pojavljuje kritiean toplotni fluks (desno) Sturgis i Mudawar [6] su razvili teoretski model za predvidanje kriticnog toplotnog fluksa u dugackom kanalu pravougaonog poprecnog preseka koji је izlozen zagrevanju sa jedne strane, slika 2.6. Pojava kriticnog toplotnog fluksa, tj. velikih parnih mehurova na zagrejackom zidu је predstavljena kao periodicna sinusoida cija 8 2. Pregled dosadaSпiih rezultata iz oЫasti krize kljucanja se amplituda i perioda povecavaju u pravcu strujanja. Modeli strujanja omogucavaju predvidanje kriticne periode па razdelnoj povrsinj korisceпjem analize пestabllnosti brzine faze i prosecпe deЫjine mehшa. Pretpostavljeno је da se preпos toplote odvija samo na okvaseпoj povrsini. Utvrdivanje odnosa izmedu dиZine okva5ene povrsine i duzine zasusene povrsine prekrivene parom, koji su doЬijeni па osпovu vizuelizacije strujaпja, predstavljaju glavni doprinos ovog modela. 4 zagrejaclС> еС:>ь.~ %fdJ -- /- 11 z• tecnost , 1 1 ~~ %: Slika 2.6 Idealizovana granicпa povrsina voda-para koja pokazuje polo2.aj tecne i parne faze i zagrejacke segmente Кritican toplotni fluks se javlja kada se kolicina kretanja рате faze normalna na zid naglo poveca ujedinici vremena i nadmasi silu pritiska na razdelnoj povrsini. Pioro i dr. [7, 8, 9] su ispitivali osnovne parameter koji uticu na koeficijent prelaza toplote pri bazenskom kljucanju (kljueanju и velikoj zapremjni vode), kao sto su: toplotni fluks, pritisak saturacije i termofizicke osoЬine radnog fluida. Narocitu pa:Znju su posvetili karakteristikama zagrejacke povrsine. Proucavali su kako na bazensko kljucanje uticu: termofizicke osoblne materijala zagrejackog zida, dimenzije, deЫjina i obradenost povrsine zida, mikrostruktura materijala zida itd. Predlozili su empirijske izraze na osnovu kojih se odreduje koeficijent prelaza toplote, poglavlje 6. Не i dr. [1 О, 11] su predstavili model pomocu koga је kriva kljucanja odredena numericki na osnovu proracuna deЫjine makrosloja. Pokazano је da isparavanje usled rasta mehura najvise utice na ukupan toplotni fluks. Makrosloj koji sadrZi paru okupira sloj uz zid. NajvaZniji deo ovog modela је da se para ne formira samo na delu gde se javljaju centri nukleacije, vee i na dodiru tecnog klina sa parnom fazom. Opis modela је dat na slici 2. 7. 9 2. Pregled dosadasnjih rezultata iz oЫasti k:rize kJiucanja malcrosloJ iwr• r1 ~ ~ 1 1 о 1 malcrosloJ 1 н Slika 2. 7 Sema modela mak:rosloja [ 1О, 11] Theofanous i dr. [12, 13] је izveo eksperimente pomocu kojih је analizirao mehurasto bazensko kljucanje. Ovi eksperimenti su sprovedeni pod izuzetno dobro kontrolisanim uslovima. Koris6ena је brza infra crvena kamera, sa visokom rezolucUom za posmatranje dinamicko-termickih karakteristika zagrejacke povrsine za razlicite toplotne flukseve, pocev od pocetka nukleacije do pojave krize kljucanja. lspitivani su elektricni grejaci ciji su materijali razlicite hrapavosti. Nanomorfologija i hemijske karakteristike zagrejacke povrsine su snimljene elektronskim mikroskopom i spektoskopski Х zracima. DoЬijeni su eksperimentalni rezultati nukleacije i prenosa toplote kljucanjem za visoke vrednosti toplotnog fluksa i pokazana је jasna razlika uticaja grejaca razlicite hrapavosti povrsine na pojavu krize razmene toplote. Jasno su predstavljeni nastanak, razvoj i dinamika procesa zasusenja, koji kasnije dovodi do pregrevanja. U eksperimentu pod imenom ВЕТА је koris6ena horizontalna ravna ploca koja је uniformno zagrevana. Sprovedeni su eksperimenti bazenskog kljucanja и kontrolisanim uslovima. Uslovi eksperimenta omogucavaju slede6e: (а) konstantan toplotni fluks na zagrejackom zidu, (Ь) direktno posmatranje termickog i dinamickog polja оа zagrejackoj povrsini, ukljucuju6i stvaranje i sirenje zasusenja - ovo је postignuto brzim infracrvenim slikanjem nano filma grejaca od dole, (с) direktno posmatranje polja zapreminskog udela parne faze i prodiranjem tecne faze iz zagrejane tecnosti. 10 2. Pregled dosadasnjih rezultata iz oЫasti krize kljucanja Eksperimentalna aparatura је prikazana na slici 2.8. Test posuda је napravljena od optickog Pirex stakla, а zagreva se prolazom elektricne struje kontrolisanog napona kroz grejac deЫjine 140 nm. Eksperimenti su sprovedeni za dve konfiguracije test sekcije. Za konfiguraciji А, slika 2.8, test sekcija је staklena cetvrtasta posuda, zatvorena na dnu sa zagrejackim elementom, 20х40 mm. Grejac ima dimenzije koje su 8 i 16 puta vece od duZine kapilara, tako da test sekcija predstavlja Ьeskonacnu ravnu plocu koja se uniformno zagreva. Posuda test sekcije је napunjena kljucalom vodom do visine od 2,5 mm i nalazi se na atmosferskom pritisku. Ove dimenzije su mnogo vece od minimalnih koje su zahtevali Gogonin i Kutateladze [ 14]. ZЬog toga је sproveden eksperiment sa test sekcijom podeljenom и 8 celija, dimenzija 1 ст, doЬijenih ubacivanjem jajaste strukture napravljene od folije od nerdajuceg celika deЫjine 200~tm . U kon:figuraciji В, slika 2.8, zidovi test sekcije su smanjeni na sirinu od lmm tako da dozvoljavaju direktno posmatranje. zagrejac (20 mm х 40 mm) 140 :i:m '1 t .....Рј ===.====il 7 1 'VI ,.т Gl.w IR velike brzine (l lH Hi ЧЈ µ n1) KONFIGURACIJA А KONFIGURACIJA В kamera (20 kJh. ЈО µ с) ' c.'u-para laser OQOQlnm 'Р ~" __________ ?------------· ' zagrejac Slika 2.8 Sematski prikaz aparature koja је koriscena u ВЕТ А eksperimentu Test posuda (konfiguracija А) је prikazana na slici 2.8, zajedno sa IR kamerom. Koriscena је Kodakova digitalna kamera sa brzinom od 20000 slika и sekundi za rezoluciji 34х 128 pixela. Dve kamere moraju biti sinhronizovane. Specijalna optika је koriseena и eksperimentu. Мnogo Ъоlја prostoma rezolucija је omogucena sa specijalлim dodatnim socivima. Zagrejacki elementi su napravljeni talozenjem metala na staklu. U glavnoj senJI ВЕТ А experimenta, koriscen је titanijumski film od 140, 270, 300, 450, 500 i 1 ООО nrn rastopljen na deЫjine 130 µm od borosilicijumskog stakla. Film se zagreva pri prolasku elektricne struje Gednosmeme ili naizmenicne), i generacija pare је unifonnna pri uslovima kljucanja. Svaki eksperiment pregrevanja dovodi do ostecenja ~ zagrejaca, sto је veoma Ьitno za osiguranje rada sa ovim zagrejacima. Neophodno jer . · 11 2. Pregled dosada8njih rezultata iz oЫasti krize k}jucanja osigurati da se otkazi ne javljaju pod uslovima neophodnim za postizanje pregrevanja i to zahteva znacajan razvoj procedure proizvodnje i izbor odgovarajuceg materijala. Generacija раге na nanofilmu zavisi od lokalnog elektricnog otpora koji varira sa temperaturom na nacin kako је to prikazano na slici 2.9. Za kljucanje od 1 00-150°С, lokalni toplotni fluks varira oko 5%. Ispod 1 50°С, kada dolazi do foпniranja suvih tacaka i temperatura se naglo povecava, gubi se snaga lokalno i ovo dovodi do variranja toplotnog fluksa za oko 30% i porasta temperature do 380°С. Pregrevanje se desava pre nego sto dode do otkaza i znacajnog gubitka snage. Povrsina glatkog mikrofilma dozvoljava radne uslove sa tacno definisanim karakteristikama. Tipicna mikroskopska elektronska slika је prikazana na slici 2.1 О. Primecuje se slicnost и gustini i hrapavosti. Amplitude, su relativno uniformne sa rms vrednostima +-4 nm. Х zraci pokazuju primarni pik na 38,5° i manje pikove za uglove 53°, 70°, 83°. Ovi podaci pokazuju osnovnu kristalnu orijentaciju (О О 2), zajedno sa malom kolicinom slucajne raspodele (1 О 2), (1 О 3), (1 О 4). о з -;:: о 2 о~~-'-~~_.__~~~~~~~~~~ 100 њо 200 ;гsо зоо зsо ~оо Т e111peratшa {oq Slika 2.9 Zavisnost temperature i elektricnog otpora Grejaci su izlozeni starenju pulsnim zagrevanjem и vazdusnom ili parnom okruZenju. Ispitivanje zagrejacke povrsine posle starenja pokazuje pojavu ostrva sa nehomogenom strukturom. Pokazana је losa konduktivnost па ovim ostrvima na materijalu. Ova ostrva su oksidi koji se foпniraju tokom procesa zagrevanja. Vaquila, [ 15] је pokazao da za temperature ispod 200°С, oksidacija se karakterise prisustvom samo jedne faze oksida : Ti02. Za vise vrednosti temperature se javlja Ti20 3• Нrapavije povrsine zahtevaju dиZi period kljucanja, koji dovodi do oksidacije i do necistoce и vodi koja se koristi za eksperimente. 12 2. Pregled dosadasnjih rezultata iz oЫasti krize kljucanja а ь Slika 2.1 О Povrsina nanofilma Kontaktni ugao је meren i za glatke i za hrapave zagrejace na 60° -75°. Ni jedna od ovih karakteristika se ne menja ako se rukuje sa instalacUom do otkaza zagrejaca. Gustina nehornogenosti se povecava sa povecanjem broja ponovljenih ciklusa zagrevanja i njihovog trajanja. Poznato је iz literature, а to pokazuju i ВЕТА eksperimenti, da parna okolina znacajno ubrzava oksidaciju titanijumskog tankog filma i izaziva formiranje ostrva oksida i Ьidroksidne grupe. Cak i pulsno zagrevanje vazduha utice na prisutni nivo vlage. Ponovljeno kljucanje i pulsno zagrevanje dovodi zagrejac и stanje teskog starenja, sa veoma velikom gustinom nehomogene povrsine. Ovi zagrejaci daju vece vrednosti kriticnog toplotnog fluksa. Као dodatak nanofilmu, koristi se 5 cm debeo zagrejac od bakra, koji se snabdeva strujom i smesten је и maloj test sekciji. Tesko starenje ove bak.arne povrsine se doЬija od delova aluminijum oksida na ispoliranoj bakarnoj povrsini. Na osnovu dobro zadatih pocetnih uslova, zagrejaci su tretirani razlicitim metodama do razlicitih stepena hrapavosti. Najefektniji nacin povecanja hrapavosti povrsine је periodicno zagrevanje do (350-400)0С, do 2 s, zatim naglo hladenje, а zatim ostaviti da stoji vremenski period pre nego sto se test sekcija napuni vodom. Eksperiment је izveden sa vodom velike cistoce. Koriscena је IR kamera za monitoring prisustva i odvajanja gasnih mehurova, preko njihovih otisaka па zagrejacu, koji se lako identifikuju kao svetle tacke. Grejac је kompletno kvasen i eksperimeпt је raden za razlicite vrednosti toplotnog fluksa, koji odgovara naponu grejaca. Kada se svaki nivo toplotnog fluksa staЬilizuje odredeni vremenski period, IR kamera se aktivira. Broj centara nukleacije је bitna karakteristika procesa kljucanja za definisanje topologije i mikrohidrodinamike. Znacajni parametri su jos i velicina mehura, dinamika rasta mehura i termicke karakteristike zagrejacke povrsine и odsustvu i izmedu pojave nastanka mehura. Rezultati snimanja doЬijenih ВЕТА eksperimentom su dati па slici 2.1 1 i 2.12 za povrsine razlicite starosti materijala. Ovo su kompletni originalni rezultati, sa kompletnom temperaturskom skalom i celom grejnom povrsinom. Ove temperature su predstavljane sivom skalom, gde svetlije Ьоје odgovaraju visim temperaturama ( oko 13 2. Pregled dosada5njih rezultata iz oЫasti krize kljucanja 150°С) na slici, а tamnije Ьоје odgovaraju niZim temperaturama (oko 100°С). Na pokretnim slikama se mogla videti dinamika rasta mehura, od nukleacije, preko rasta mehura i odvajanja i nestajanja. Brojenje ovih tamnih delova је radeno automatski uz pomoc odgovarajuceg software, koji је razvijen prilikom ovog istrai.ivanja. (а) (Ь) (с) Slika 2.11 Slika IR termometrije glatkog zagrejaca za tri razliCite vrednosti toplotnog fluksa~ q=406, 536 i 807 kW/m2 (а) (Ь) (с) Slika 2.12 Slika IR termometrije starog zagrejaca za tri razlicite vrednosti toplotnog fluksa, q=348, 1051 i 1517 kW/m2 Na slici 2. 13 је dato kako gustina centra nukleacije zavisi od toplotnog fluksa i primecuju se mnogo vece vrednosti gustine na hrapavim povrsinama nego na glatkim. Razlika и gustini centara nukleacije izmedu hrapavih i glatkih povrsina је oko jedan red velicina. Takode se primeeuje linearna zavisnost toplotnog fluksa, efekata kvaliteta vode i efekata starenja povr5ine. Ostale karakteristike su : neuniformnost povecanja nukleacije sa toplotnim fluksem, kako pregrevanje postaje intenzivnije; cak i za izuzetno velike vrednosti toplotnog fluksa, Ыizu kriticnog toplotnog fluksa, postoji sporadicki distribuirana povrsina koja је pregorevana u odsustvu nukleacije; u svak:om frejmu postoji spektar mehurasto 14 2. Pregled dosadaSnjih rezultata iz oЫasti krize kljucanja bladenje povrsine, ali su znacajniji za sveze zagrejace. Hrapave povrsine zagrejaca se pona5aju mnogo unifoпnnije nego glatke povrsine. Prelaz toplote pri mehurastom bazenskom kljueanju је proucavano i intenzivno mereno u proslosti. Опо se karak:terise zavisnoscu izmedu pregrevanja povтSine i toplotnog fluksa i naziva se kriva kljucanja. Na osnovu rezultata eksperimenta se moze zakljuciti da prosecna temperatura zida za kritican toplotni fluks varira od 22 do 32 К za glatke zagrejace, dok је za hrapave povrsine to variranje od 18 do 25 К. Za ove uslove odgovarajuCi vrh za pregrejace је 60 К, odnosno 40 К za glatke i za grube povrsine respektivno. Pokazana је skoro linearna zavisnost za hrapave povrsine, slika 2.14, dok se za glatke zagrejace, slika 2.13 kriva sastoji od dva dela: и pocetku Ыagi deo sa rastom, а zatim skoro vertikalni skok za visoke vrednosti toplotnog fluksa. 1000 о F1 \1 F4 о (1 F9 800 ~ о '1 \1 (\/- 600 ~ џ\1 о Е ~ о ()." о .:.: - '1" 8 0- 400 о '1 '1џ о 200 ~ џ о \l 'О У' ф о . . 5 10 15 20 25 зо 35 Л Ts (К) Slika 2.13 Кriva kljucanja dobijena na osnovu ВЕТА eksperimenta za glatku povrsinu zagrejaca Na slikama 2.15 i 2.16 su prikazane krive kljucanja dobijene eksperimentalnim putem za hrapavu i za glatku povrsinu zagrejaca. 15 2. Pregled dosada5njih rezultata iz oЫasti krize k.ljucanja Slik:a 2.14 Кriva kljucanja doЬijena na osnovu ВЕТА eksperimenta za hrapavu povrsinu zagrejaca 10 [ill 10 о FI 1 4 о v F4 о "' о F9 l 8 ·Г 8 о о ~ о о о о v • v " = 6 о v ~ е о v ... о " о "' о • о " о ~ v :;;: о о = v = • ovo = 4 d3 " "' v ~ Ov ~ v о • 2 <;pV " vV " " о g Vo 1 2~ ; Ј> i ov = о => о & v о 200 400 600 800 1000 о q (kW/mz) 5 10 15 20 25 зо ~ Ts (К) Slika 2.15 Кriva kljucanja doЬijena eksperimentalno za glatku povrsinu 60 О А11 оо D А.1 о о А1 1 ~ 50 ом [] АЗ - 50 Q м "( Е 40 о з. " о 40 ~ оо ф о ~ " = " ЈО о ... ~ r: ~ о о ф зо о = о ~ о о r: 20 о о = о Е оо r: о ~ 20 о • о о оо~ " ф "' 10 о Оо " = о о .Ё 10 оо 8 ~ о о :, о ;; о о ;, о 200 400 600 800 1!ЈОО 1200 1400 1600 о о q (kW/m2) 5 10 15 20 25 зо \ Ts (К} Slika 2. 16 Кriva kljucanja doЬijena eksperimentalno za hrapavu povrsinu 35 16 2. Pregled dosada.Snjih rezultata iz oЫasti krize kljucanja Kureta i dr. [ 16, 17] su vrsili eksperimente и kanalu kvadratnog poprecnog preseka koji је zagrevan sa jedne strane. Rezultati eksperimenta su prikazani па slici 2. 17. parna faza 1·] 1 i' 1 о 100 struJanje direktno rastoJanJe (mm) о ::-15 poprecna pozicija (mm) (b )10.7m1 (С)2 1 Зms "parna faza 1 Е Е 0.5 8 о Slika 2.17 Trenutna raspodela zapremjnskog udela tecne faze doЬijena eksperimentalno U okviru ovog rada је predlozena nova korelacija za racunanje kriticnog toplotnog fluksa za kanale cetvrtastog poprecnog preseka i cevi malog precnika koristeci podatke za kriticnu vrednost zapreminskog udela tecne faze, bezdimenzionog parametra kriticnog toplotnog fluksa i odnosa zagrejackog oЬima. Tacnost doЬijenih rezultata је do 45 %. Predlozena је eksperimentalna korelacija za odredivanje vrednosti kriticnog toplotnog fluksa: (2.2) Stosic i Stevaпovic [18] su izvl'Sili numericku simulaciju procesa zasusenja na horizontalnoj povrsini bazenskog kljucanja. U ovom radu је modeliran i mikro nivo sa nastankom i rastom mehura i makroskopsko pona.Sanje dvofazne mesavine za veoma velike vrednosti toplotnog fluksa. Mikro nivo na zagrejackoj povrsini је modeliran pomocu parametara Ьitnih za generaciju pare na zagrejackom zidu, kao sto su: gustina centra nukleacije mehura, vreme zadrZavanja mehura па zagrejackom zidu i odredeni nivo slucajnog odablra lokacije centara nastanka mehurova. Ovaj nacin modeliranja је preuzet iz ovog rada i primenjen и doktorskom radu. Toplotni fluks је neuniformno rasporeden ро zagrejackoj povrsini sa pikovima na mestima centara nukleacije. Centri nukleacije su odredeni funkcijom slucajnih brojeva, dok se srednja vrednost centara nukleacije odreduje na osnovu karakteristika materijala i hrapavosti zagrejacke 17 2. Pregled dosadasnjih rezultata iz oЬlasti krize kljucanja povrsine. Makro nivo је modeliran pomocu modela dva fluida. Transportni procesi na razdelnoj povrsini faza su modelirani pomocu odgovarajucih konstitutivnih korelacija. Rezultati su dati za veoma kratak period posle pocetka zagrevanja i generacije mehurova na zagrejackoj povrsini, kao i za kvazi-stacionarne uslove nekoliko sek:undi posle pojave bazenskog kljucanja. Simulacija је radena i za nize vrednosti toplotnog fluksa, ali је fenomen zasusenja povrsine modeliran za velike vrednosti toplotnog fluksa. Uticaj gustine centara nukleacije i vreme zadriavanja mehurova na zidu na dinamiku bazenskog kljucanja је proucavan. Takode је ispitan i uticaj intenziteta toplotnog fluksa na dinamiku bazenskog kljucanja. Pokazano је da smanjenje gustina centara nukleacije i povecanja vremena nastanka mehura na zagrejackoj povrsini dovodi do smanjenja vrednosti kriticnog toplotnog fluksa. Takode, povecanje gustine centara nukleacije dovodi do povecanja vrednosti kriticnog toplotnog fluksa. Ovo znaci da ukoliko Ьi Ьi lo moguce predvideti gustinu centara nukJeacije, na osnovu hrapavosti zagrejacke povrsine, hemUskih i fizickih osoblna i zida i fluida, onda је moguce predvideti odgovarajuce vrednosti kriticnog toplotnog fluksa. 0.9 0.8 o.s 0.8 0.7 0.7 е 0.6 :§: Е 0.8 а o.s ('; r-- : o.s 11 и; ::!: o . .i > о.• 0.3 0.3 02 0.2 0 .1 0.01 0.1 о.о о.о Х = 35 П1111 siri11a (ш} 70 0 .0 7 60 t ,.. 0 .06 •• 50 "' .:. !: 40 ~ 30 20 10 ,.. ..,,, ./ i ~ 11"'" 1~ 411 "'-- 0 .0 5 .s.. 0 .04 ~ ~ о .оз 0 .02 0 .01 о о 0.2 0.4 0.6 0.8 1.0 0.2 0 .4 0 .6 0.8 1 z~p1t111lnskJ udeo ра111е faze Slika 2.18 Uporedna analiza eksperimentalnih rezultata i predlozenog modela za toplotni fluks od 1000 kW/m2 18 2. Pregled dosadasnjih rezultata iz oЫasti krize kljucanja I-liblki i Ishii [19] su koristili gustinu centara nukleacije kao znacajan parameter za odredivanje dodime povrsine и modelu dva tluida. Gustina centara nukleacije је modelirana poznajuci velicinu i ugao kavitacije koji se javljaju na povтSini. Pokazano је da је gustina centara nukleacije funkcija kriticne velicine kavitacije i kontaktnog ugla. Payan-Rodriguez i dr. [20] su razvili metod za predvidanje kriticnog toplotnog fluksa za kljucanje vode и vertikalnim cevima pod standardnim uslovima koji se javljaju и generatorima pare, slika 2. 19. Pretpostavlja se da se kritican toplotni tluks javlja kada tecni sloj na zidu cevi ispari pri velikim vrednostima toplotnog fluksa. Razvijen је model na osnovu teorije sloja zasusenja, uzimajuci и obzir metod faktora oЫika, da Ьi se uzela и obzir i neuniformna raspodela toplotnog fluksa na vertikalnoj cevi, sto је slucaj и generatorima рате. Ovaj model se ne moze primeniti na slucaj eve kada se kritican toplotni tluks javlja pri maHm vrednostima zapreminskog udela рате faze. tecni Щу~; ! podsloj 1 Ui • parni 1 mehur . 1 Slika 2.19 Model sloja zasusenja .... - 118\HJ•I - t)фtrlюti.1 - fe111"er.wr1111 d•rvolJ••• v1•d• oitl ltn11,11a1ur• .... - blti<•• i.pl"'11i flllk~ 4o>•о.оз о.о 0.01 О О.25а О.5а О.75а а Sirina (m) Pocetni nivo оо о о --~---~- Ј_! + 4 5 6 i 7 1 8 Zone 1 9 10 11 1 12 i ' 13 14 15 16 Slika 4.1 Primer bazenskog kljucanja (40х40х40) Pami meburovi se generisu па dnu posude. Za vreme kljucanja nivo tecnosti se menja i predvidanje polo:Zaja nivoa је takode ukljuceno u model. Zapremina iznad nivoa teenosti је ispunjena parom. Dno posude је podeljeno na zone, slika 4.1. Dimenzije zone zavise od hrapavosti zida i odgovaraju gustini centara nukleacije (4х4 zone, kvadratnog oЬlika dimenzija 0.01 m х 0.01 m za gustinu nukleacije 1 1/cm2, 0.003 m х 0.003 m za gustinu nukleacije 11 l/cm2, 0.002 m х 0.002 m za gustinu nukleacije 25 l /cm2, 0.0016 m х 0.0016 m za gustinu nukleacije 39 l /cm2, 0.0014 m х 0.0014 m za gustiпu nukleacije 51 l /cm2). U svakoj zoni se javlja centar nukleacije i svaka zona se sastoji od deset kontrolnih zapremina. Lokacija na kome se generise mehur unutar zone zavisi od funkcije slucajnih brojeva, koja је koriscena da bi se odredilo polozaj mehura. Pretpostavljeno је da se tokom vremena samo jedan mehur generise unutar jedne zone. Pretpostavlja se da se najveci deo toplote prenosi na centar nukleacije, а konvekcija sa zida na jednofazni fluid se zanemaruje. 4.1.1. Parametri mikro nivoa Dimenzije posude su uzete tako da predstave Ьeskonacnu ravnu geometriju. DuZina kapilare vode је uzeta za skaliranje Ьeskonacne geometrije i data је sa: (4.1) 31 4. Razvoj modela krize kljucanja Vrednost za Lc је 2.s·10-3 m za atmosferske uslove. То је nekoliko puta manje od sirine posude koja se simulira и ovom slueaju. Dno је podeljeno na zone. Sirina zone zavisi od gustine centra nukleacije. Odnos izmedu gustine centra nukleacije п i sirine zone Ь је dat sledeeim izrazom: 1 Ь= ј;, (4.2) Za sirinu od О.О 1 m, dobija se da је gustina nukleacije 104 centara /m2. Gustina centara nukleacije zavisi od toplotnog fluksa, hrapavosti zagrejacke povrsine, kontaktnog ugla kva5enja i termo fizickiЬ karakteristika radnog :fluida i zagrejackog zida. Vrednosti za gustine centara nukleacije su izmedu 103 sites/m2 to 107 sites/m2 za atmosferske uslove (18]. Ako је hrapavost povrsine manja, broj centara nukleacije је manji i manjaje gustina centara nukleacije. Na kontaktni ugao kva5enja utice hrapavost povrsine i uslovi povrsine (starost povrsine). Precnik mehura и trenutku odvajanja od zida se racuna na sledeci nacin: (4.3) U tabeli 4.1 su dati parametri mikro nivoa: kontaktni ugao, precnik mehura i gustina centara nukleacije za povтSine razlicite starosti. Нrapavost sveze povrsine је mala, za male vrednosti uduЫjenja; otuda, kljucanje na svezoj povrsini se karak.terise sa malim vrednostima gustine centara nukleacije. Tak:ode, kontak:tni ugao kvasenja se povecava za glatke zagrejace. Uticaj povrsinske hrapavosti је predstavljen preko kontak:tnog ugla fJ. Za grube povrsine, kontaktni ugao od 5 stepeni је pretpostavljen i za atmosferske uslove precnik mehura koji odlazi је 2.6·I04 m. Za visoke vrednosti toplotnog fluksa, i za usJove bliske napustanju mehurastog kljucanja, moze se pretpostaviti da је cela zagrejacka povrsina prekrivena sa mehurovima, koji su dati prostom geometrij skom vezom da је jedan metar kvadratni prekriven sa п mehurova precnika Dь: пп; = 1 (4.4) gde је п gustina centra nukleacije. Za hrapavu povrsinu 2.6·104 m је gustina centra nukleacije l.s·107 m-2 • Za glatke povrsine, kontaktni ugao od 40 stepeni, moze se pretpostaviti Dь=2.08· 10-3m i п=2.3 ·10-5 . Na osnovu prezentovanih centara nukleacije sa zonama datih na slici 4.1 , i za uslove potpune pokrivenosti mehurovima zagrejacke povтSine, sirina zone moZe Ьiti jednak:a preCniku napustajuceg mehura. Sirina zone u milimetarskom opsegu, koja је prikazana za odgovarajuce glatke zagrejace, moze se resiti sa numerickom mrezom za prihvatljiv broj kontrolnih zapremina. Tabela4.l. Parametri mikro nivoa () Dь п gladak zagп~jacki zid 40 2.08·1 о-.) 2.з· 10=> hrapav zagrejacki zid 5 2.6·10 .... 1.s·1 О ' 32 4. Razvoj modela krize kljucanja Drugi parametar koji odreduje dinamiku bazenskog kljucanja је vreme boravka mehura na zagrejackoj povrsini. То је precnik mehura za vreme od rasta mehura do napustanja. Takode se moze predvideti zavisnost izmedu precnika mehura i vremena rasta mehura: (4.5) (4.6) Za kontaktni ugao izmedu 40° i 90 ° је у = 0.1 do у = 0.49, respektivno. Empirijski parametar јЈ = 6. Korelacija za Dь је primenjiva и opstem oЫiku kada se toplota prenosi i rastorn mehurova od zagrejacke povrsine do osnove mehura i od zasicenog tecnog sloja oko mehura. Vreme nastanka mehura је dato и funkciji od ugla kva5enja i Jacob-sovog broja, na osnovu jednacina (4. 5) i (4. 6): (О .0208х 9)2 L~ т =~-,-'--~-;=========~ 4а(уЈа +~r Ја2 + 2/ЈЈа) (4.7) Na osnovu Тhomove korelacije h = 1.9712 (2Р/8687) (Т -Т ) е w sal (4.7а) se doЬija da vrerne rasta i odvajanja mehura na mestu nukleacije zavisi direktno od toplotnog fluksa, Ја broja i kontaktnog ugla. Slika 4.2 pokazuje zavisnost vremena boravka mehura na mestu nukleacije od ugla kvasenja i toplotnog fluksa. Predvidanje gustine centara nukleacije za zagrejace pomocu gore navedenihjednacina (4.2) i (4.4) је и opsegu 107 centara/m2 i mnogo su veCi nego rezultati direktnih merenja, koji su doЬijeni ро prvi put brzim infracrvenim snimanjem uz pomoc nano filma zagrejaca. Rezultati direktnog merenja na hrapavoj povrsini su dati na slici 4.3, i nisu veci od 6·105 centara/m2• Takode, prikazano је povecanje gustine centara nukleacije sa poveeanjem toplotnog fluksa. Ovo poredenje pokazuje da analizirane korelacije nisu pouzdane za predvidanje gustine centara nukleacije na hrapavoj povrsini. 33 ro с: ro .... (Ј) ::i - .r:; ·-ф .5: Е ~ 111 > ·~О с: а. 111 ·~ > о 111 ~ N џ .... ro "О·~ ro Ф N ОЈ Ф ro Е N ф .... > 2.50Е-02 " _ -· 2.ООЕ-02 :\. ' '\ '\ 1 .50Е-02 ' 1 .ООЕ-02 ~ 5.ООЕ-03 . О .ООЕ+ОО 200 1 ' ' ' 1 . 400 4. Razvoj modela krize kljucanja ··-~·--г -т - kontaktлi ugao = 40 s!epeni - - - kontaktлi ugao = 90 stepeni - -;...... 1 -- - - ,__ -- ----- ' 600 800 1000 1200 1400 1600 Toplotлi fluks, kW/m2 Slika 4.2 Zavisnost izmedu vremena rasta mehura i toplotnog fluk.sa za razlicite kontaktne uglove 60 -N Е Q -A1 (..) 50 Q О -АЗ ~ - О -А4 :~ 40 џ оф оо t; ф ~ зо џо - о = 1:1 о о ... ~ 20 о -= о ооо ф (..) t; 10 сов = а ·~ о ::::: о С'! о о 400 800 1200 1600 to1)lot11i л.11ts . KW/m2 Slika 4.3 Eksperimentalna zavisnost gustine centara nuk.leacije i toplotnog fl'uksa za hrapave povr5ine ( Zagrejac Al је pulsno zagrevan na vazduhu, а zagrejaci А3 i А4 su tretirani zagrevanjem i kljueanjem и vodi. ) Za glatke zagrejace su izmerene vrednosti centara nuk.leacije do vrednosti od 105 centara/m2 slika 4.4. Pokazano је da se gore navedene jednacine mogu koristiti za odredivanje gustine centara nukleacije na glatkim povrsinama. Treba naglasiti da vece razlike u gustinama centara nukleacije su pod istim vrednostima kontaktnog ugla kva5enja. Ovo је и suprotnosti sa prethodnim verovanjem da је kontaktni ugao kljucni parametar koji odreduje gustinu centra nuk.leacije. 34 4. Razvoj modela krize kljucanja 10 -N Е v о 8 о ..:: о ф о о ....... \/ о ·~ 6 о /: оО '\! ф ~ о - 'V - !:: 4 OvO о О - F1 ... !': = О;т O -F2 ф v 2 \$5v 'V - F4 !:: - O -F9 -~ 'V') \70 :) о о .200 400 600 800 1000 toplot11i flt1 ks , kW/m2 Slika 4.4 Eksperimentalna zavisnost gustine centara nukleacije i toplotnog f1uksa za glatke povrsine ( Eksperiment F9 је raden sa cistom destilovanom vodom, а Fl, F2 i F4 sa vodom najvece cistoce ) Gustina centra nukleacije и opsegu 105 centara/m2 se moze resiti na domenu sa prihvatljivim brojem kontrolnih zapremina. Za vece vrednosti pritiska, gustina centra nukleacije treba da dostigne mnogo vecu vrednost i rezolucija svakog centra unutar mreze moze biti prakticno postignut. Jedino razumno resenje za visoke vrednosti gustina centara nukleacije је da se primene prosecni centri nukleacije i generacije раге ро citavoj sirini zagrejacke povгSine, ali to ne znaCi da toplotni fluks mora biti konstantan ро sirini zagrejacke povrsine. Za generaciju mehurova na zagrejackoj povrsini se koristi funkcija slucajnih brojeva. Broj mehurova na povrsini zavisi od hrapavosti povrsine, ali mesto na kome se mehur generise је odredeno uz pomoc funkcije slucajnih broj eva. RND funkcija vraca vrednost izmedu О i 1. Vrednost broja odreduje kako RND generise slucajan broj: Tabela 4.2. RND funkcija Ukoliko је Ьтој Rnd generise Manji od nule Isti broj svaki put, koristeci pocetni broj Veci od nule SledeCi slueajan broj u nizu Jednak nuli NajcesCi broj koji se generise Za svaku datu pocetnu vrednost, isti broj iz niza је generisan, jer svaki suksesivan poziv RND funcije koristi prethodni broj kao pocetni za sledeci broj и nizu. 35 4. Razvoj modela krize kliucanja 4.1.2. Jednacine odrianja Modelirana su dva razlicita regiona: dvofazno strujanje u bazenu i konduk:cija и z.agrejackorn zidu. Dvofazno strujanje је modelirano koriseenjern modela dva fluida. ReSavane su jednacine od.ГL.a.Dja mase, koliCine kretanja i energije za оЬе faze posebno, а transportni procesi na razdelnim povrsinama su modelirani pomocu konstitutivnih korelacija. Dvofazno strujanje је modelirano kao polukompresibilno, zanemaren је akusticki efekat. Pretpostavljeno је da nema promene pritiska izmedu tecne i parne faze. PovгSinski napon је zanemaren i isti је pritisak za оЬе faze unutar kontrolne zapremine. Dinamika fluida: Jednacina odrZa.nja mase: (4.8) Jednacina odrZanja kolicine kretanja: Jednacina odrZanja energije: д(akpk~) + д(akpkuk.i4 ) =~[k 84 - Т.и~]+ (-1)k (Г -Г )Т. +(2 - k) · /с дt дх1 дх. k дх. Р* * , е с k qь p,k 1 1 (4.10) gde је k = 1 za tecnost i k = 2 za paru. Pararnetri и, р, h i Т su usrednjene brzine, pritisak, entalpija temperatura, respektivno. Као posledica vremenskog usrednjavanja dobijaju se clanovi: (4.11) sto је tenzor Rejnoldsovih napona i vektor turbulentnog toplotnog fluksa. Ovi turbulentni efekti su zanemareni i pretpostavljeno је laminamo strujanje. Indeks k је jednak 1 za vodu, 2 za paru i 3 za zid. Izvorni clanovi za masu, kolicinu kretanja i energiju su dati и Ьilansnim jednacinama. Intenzitet tranzicije faza, tj. brzina i~paravanja ili kondenzacije su dati sa Ге i Гс. Sila medufaznog trenj a је data sa F12.i· Clanovi F 12,i i F vм,i su uzgonska sila i sila virtualne mase, respektivno. Clan qзк predstavlja zapreminski toplotni fluks sa zida na odgovarajucu fazu. 36 4. Razvoj modela krize kljucanja Da bi se zatvorio sistem jednacina, potrebno је dodati i: (4.12) Energetska jednacina z.a zid: (4.13) gde је qh zapreininski izvor toplote и zidu, elektricru grejac, dok је qь је toplotni ponor и kontrolnoj zapremini na z.agrejackom zidu, koji se javlja zbog rasta mehшa (u isto vreme ova vrednost predstavlja i toplotni izvor и kontrolnoj zapremitll fluida na zidu kad se pojavi rast mehшa) slika 4.5. Pretpostavlja se da se nukleacija mehшa ne javlja, tj topotni izvor qь је jednak nuli, ako zapreminski udeo рате faze na zagrejackoj povrsitll prede 0.9. 4.1.3. Konstitutivne korelacije Da bi se sistem jednacina mogao resiti potrebno је definisati izraze z.a sledeee velicine: izvorlli clan koji potice od medufaznog trenja, uzgonska sila, brzina isparavanja, brzina kondenzacije. Izvorlli clan koji potice od medufaznog trenja se definise na sledeci nacin: (4.14) gde је Cv koeficijent medufaznog trenja, а Dp Је precnik mehшa. Koe:ficijent medufaznog trenjaje z.adat sledecim izrazom: (4.15) Uzgonska sila: ( ~ди1. U -U ,Ј 2,) 1,) дх. 1 ( ~диl . - и -и '1 2.) \,ј дх . Ј (4.16) gde је CL=0.3. Brzina isparavanja i kondenz.acije se racuna na osnovu empirijskog modela koji uzima и obzir relaksaciono vreme i-. Intenzitet isparenja iznosi: 37 Brzina kondenzacije: Za zagrejacku p locu vaii jednacina odrZa.nja Kolicina topote koj a se dovodi ploCi: • 2 • qgr а Ь=Qь 4ь Лх Лу Лz = Qь п ( q: Лх' Лу' Лz') = (4ь Лх Лу Лz) fluid zid 4. Razvoj modela krize k.ljucanja (4.17) (4.18) (4.19) (4.20) (4.21) (4.22) (4 .23) Pretpostavlja se da se nukleacija mehura ne javlja, tj topotni izvor qь је jednak nuli, ako zapreminski udeo parne faze na zagrejackoj povr5ini prede 0.9. -/ !!: \1 qlV Slika 4.5 Promene и zidu usled pojave mehura и dvofaznoj mesavini Ыizu zagrejackog zida 38 4. Razvoj modela krize kljucanja 4.1.4. Prelaz toplote izmed:u zida i dvofazne mesavine Najces6e koris6ena procedura za odredivanje kondukcije izmedu razlicitih faza [33] је da se pretpostavi linearna zavisnost izmedu tacaka Р i Е, slika 4.6: gde se faktor interpolacije definise na slede6i nacin: (ох)". fe= (ox)e neophodno је doЬiti dobro predvidanje toplotnog fluksa na povrsini: (4.24) (4.25) (4.26) koji se koristi и diskretizovanim jednacinama. Dobro odabrana vrednost ke је ona koja dovodi do dobrog toplonog fluksa na povrsini.°C Kontrolna povrsina koja okrufuje tacku Р је popunjena materijalom uniformnog koeficijenta kondukcije kp, а oko tacke Еје materijal koeficijenta kondukcije kE. (4.27) k = ( l - fe + f. Ј- 1 i ako је granica postavyena na sredini izmedu tacaka Р i Е е k p kE k -1 = 0 s(k -1 k -1) ·1· k = 2kPkE е . р + Е 11 е • kp+ks (4.28) Za slucaj zagrejacka povrsina od bakra i dvofazna mesavina voda-para kE = 3.369·10-7 m 2/s. (ох). (ох) ... i 1 1 ~~~~6~~----1-1 ~~---<о~~~~- Р е Е Slika 4.6 Prelaz toplote izmedu zida i dvofazne mesavine Ovde predstavljen model је dat uzimaju6i u obzir efekat generacije pare i dvofaznog strujanja, kao sto је disperzija faze u okviru mesavine. Rezultati su dati za kratak 39 4. Razvoj modela krize kljucanja period vremena od pocetka zagrevanja i generacije pare na zagrejackoj povrsini, kao i z.a k:vazi stacionarne uslove posle nekoliko sekundi od pocetka bazenskog kljucanja. Dopuna z.agrejacke povrsine sa vodom i delimicno kva5enje z.a nize vrednosti toplotnog fluksa ја korisceno, dok је zasusenje posmatrano za visoke vrednosti toplotnog fluksa. Uticaj na gustinu cenatara nukleacije i vreme pojave mehurova na zidu bazenskog kljucanja је istraZivano. Takode, uticaj intenziteta dinamike bazenskog k.ljucanja је analiziran. Numericka sirnulacija pokazuje smanjenje gustine centara nukleacije i povecanje vremena pojavljivanja mehurova na zagrejackoj povrsini (za glatke povrsine) i dovodi do smanjenja vrednosti kriticnog toplotnog fluksa. 4.1.5. Uticaj zagrejaca па pojavu krize razmene toplote Zagrejacka povrsina ima veliki uticaj na gustinu centara nukleacije. Gustina nukleacije na hrapavim povrsinama је veca od 1 О centara/cm2 za toplotni fluks veCi od (400-500) kW/m2 , dok је za slucaj glatkЉ povrsina gustina manja od 10 centara/cm2 z.a toplotni fluks veci od 1 ООО kW/m2. Eksperimentalni rezultati takode pokazuju da hrapavost z.agrejacke povrsine ima uticaj na kriticni toplotni fluks. Ovaj efekat је pokazan na slici 4.7, gde su mereni СНF podaci uporedeni sa podaciroa doЬijenim sa Kutateladze -Zuber korelacijom. Slika 4.7 pokazuje da se kriticni toplotni fluks povecava sa hrapavoscu zagrejacke povrsine, dok se primenom Kutateladze-Zuber korelacije predvida СНF koji odgovara srednje hrapavim povrsinama. Kuttaeladze-Zuber-ova korelacija se zasniva na Kutateladze-ovom broju: (4.29) gde је kriticna brzina parne faze: (4.30) dok је Kutateladze-Zuber-ova jednacina za predvidanje kriticnog toplotnog fluksa data sa: (4.31) Ск је Kutateladze-ova konstanta i iznosi 0.145. Eksperimentalni rezultati su u suprotnosti sa nekim prethodnim modelima koja predvidaju povecanje kriticnog toplotnog fluksa sa smanjenjem gustine centara nukleacije. Moze se zak.ljuciti da је hrapavost zagrejaca povezana sa gustinom cenatar nukleacije. U ovom radu је gustina centara nukleacije prihvacena kao granicni uslov koji pokazuje hrapavost zagrejacke povr8ine i njen uticaj na mik:ro fenomene pojave 40 4. Razvoj modela krize kljucanja kriticnog toplotnog fluksa. UsvajajuCi razliCite vrednosti gustine centara nukleacije u numerickom pristupu, uticaj hrapavosti zagrejaca se uzima u obzir. 1.6 ..------.--------,.- ------. -~ ~ ~.8 ~<) Ео . 6 ~<>о.Јо 0- 0.4 '-----4------------- sre(lвje l1п11нtv glc:нlak Slika 4.7 Odnos izmedu merenih i izracunatih vrednosti kriticnog toplotnog fluksa 4.1.6. Granicni uslovi Na slici 4.8 su dati granicni uslovi za zagrejacki zid, za ostale zidove i na izlazu. Pretpostavlja se da је posuda adijabatski izolovana i da је na izlazu atmosferski pritisak. za adijabatski izolovane zidove ди дv дw --О --О --О дt ' дt , дt z • у о • х р=О " L..---- " 1 1 1 L---- "_ - - - "" "" ,, "" о za fluid 2· 10-зт i z i О.02т za zagrejacki zid о i z i 2·10-Зт z -о дТ -•О дz х •О, х •а у •О. у• а дТ -О дх дТ •О ду Slika 4.8 Granicni uslovi za slucaj bazenskog kljucanja 41 4. Razvoj modela krize kljucanja 4.2. Strujanje u vertikalnoj cevi cetvrtastog poprecnog preseka 4.2.1. Opis modela Cev је cetvrtastog poprecnog preseka, slika 4.9. Zagreva se bocni zid, а ostali zidovi su adijaЬatski izolovani. Granicni uslova na ulazu је profil brzina strujanja tecne faze koji se zadaje pomocu izтaza: v = 0.5 · sin (25·1z" у· 0.001) (4.32) Bocni zid је podeljen na 4х4 zone. U svakoj zoni se stohasticki odreduje polozaj jednog centra nukleacije. Pretpostavljeno је da је cev cetvrtastog poprecnog preseka ispunjena dvofaznom mesavinom, koja sadrzi 50% tecne i 50о/о pame faze. Uzeto је и obzir i delovanje zemljine teze i to u vertikalnom pravcu. Jednacine odГLanja mase, kolicine kretanja i energije su slicne kao i и prethodnom slucaju, jedino se pri modeliranju dodaje strujanje, novi granicni uslov: profil brzina na ulazu i definise se da је zagrejacki zid na bocnoj stranici cevi. q - -- q g - -- q - -- q у u х z Slika 4.9 Primer vertikalnog strujanja dvofazne mecavine u cevi cetvrtastog poprecnog preseka Granicni uslovi su dati na slici 4. l О. 42 г р = О za fluid О~ х ~ О,038т а 1 == 0.5, а 2 == 0.5 z , • у ... -----·-···--·-·-~ а о • х t ulaz z-0 Wi ... ~ •0,5 · sin(0,025·л"y) 4. Razvoj modela krize kljucanja h za zagrejaeki zid О, 038т ~ х ~ 0.04 т ат х=О.04. -=О дх дТ у ... О, у = а, - = О ду дТ zaO z=h -=О . . az Slika 4.1 О Granicni uslovi и slucaju kvadratnog poprecnog preseka 4.3. Strujanje u vertikalnoj cevi okruglog poprecnog preseka 4.3.1. Opis modela Na slici 4.11 је dat prikaz cevnog isecka koji је modeliran. Pretpostavljeno је da se и cevi nalazi dvofazna mesavina, koja sadrZi 50% tecne faze i 50% parne faze. Zagreva se spoljni zid cevi uniformno elektricnim zagrejacem. ZaoЬljeni deo cevnog isecka predstavlja deo na kome se generisu mehurovi. Podeljen је na 4х4 zone i и njima se na slucajan nacin generisu mehшovi. Na ulazu је zadat granicni uslov za brzinu strujanja tecne faze i jednak је sledecem izrazu: v = О, 5 · sin (О, 025 · я ·у) (4.33) 4.3.2. Jednacine oddanja mase, kolicine kretanja i energije za slucaj polarno- cilindricoih koordinata Dinarnika fluida: Jednacina odrlanja mase: 43 4. Razvoj modela krize kljucanja q g - -- q - -- q z I i i I Slika 4. 11 Primer strujanja dvofazne mesavine u cevi ok:ruglog poprecnog preseka Jednacina odrZanja kolicine kretanja u r, rp i z pravcu: и~ др =-а*-+ r дr (4.35) (4.36) 44 4. Razvoj modela krize kljucanja (4.37) Jednacina odrZanja energije: (4.38) gde је k = 1 za tecnost i k = 2 za рати. Energetska jednacina za zid: (4.39) 45 5. Numericko resenje sistema diferencijalnih jednacina 5. Numericko resenje sistema diferencijalnih jednacina Metoda kontrolnih zapremina se koristi za resavanje parcijalnih diferencijalnih jednacina (4.8, 4.9, 4.10) predstavljenih и prethodnom poglavlju [34]. Diskretizacija se odvija integracijom promenljive velicine и kontrolnoj zapremini и ЗD ortogonalnom koordinatnom sistemu. Jednacine odrZanja mase za tecnu fazu i entalpije za tecnu i parnu fazu su diskretizovane и kontrolnoj zapremini, slika 5.1. Jednacina oddanja mase za tecnu fazu је data u sledecem oЬliku: о о) (a1,q - а1 ,q i,j,kЛx;Лyjдzk ~~~~~~~~~~-+ Лt [ (a1.Piut ) ;+112,j,k - (а1,qщ);-1 1 2,j,k ]лу jдzk + [ (a1,qv1)i,j+112,k - (a1.Pivt);,j - 112,k ЈдzkЛх; + [ О if (u1) ;+112,j,k < О лх (5.2) Slika. 5.1 Tipicna kontrolna zapremina koja se koristi za integraciju jednacina odrZanja. Zapreminski udeo, entalpija i pritisak se racunaju za tacku (ij,k) и sredini kontrolne zapre1nine Dimenzije kontrolne zapremine se racunaju ро sledecim izrazima: ~ = xi+l / 2 -Хј-11 2, Луј =У j+l / 2 - У j-l / 2, i Лzk = zk+I / 2 - zk-1/ 2. Uvodenjem ovih izraza ujednacinu (5.1) dоЫја se: 46 5. Numericko resenje sistema diferencijalnih jednacina Щ,j,k (ai)i,j,k = ai+112,j,k(a1)i+1,j,k + ai- 1/2,j,k(a1 )i- 1,j,k + ai,j+1!2,k(a1)i,j+1,k + ai,j-1/ 2,k(a1)i,j- l ,k + (5.3) ai,j ,k+l/2 (а1 )i,j ,k+1 + Щ,ј,k-112 (а1 )i,j,k- 1 + ьi,j,k gde su g koeficijenti koji se racunaju: a;,j,k = max [ (Piut);+11 2,j,kЛyjЛzk, ОЈ+ max[ -(PiЩ);-i / 2,j,k ЛујЛzk, ОЈ + max[ (PiVi);,j+l / 2,kЛzkЛx;,O ]+ max[ -(PiVi);,j-11 2,kЛzkЛx; , O] + max[ (РЈ Wi)i,j,k+112 Лx;ЛYj , О] + max[ -(Piwt);,j ,k- 112 ЛхiЛУ ј•О ]+ (Pi)i,j,k ЛхiЛу jЛzk Лt й;,ј+1!2,k = max [-у max.(x,y) = у, for у>х (5.4а) (5.4Ь) (5.4с) (5.4d) (5.4е) (5.4f) (5.4g) (5.4h) (5.5) 47 5. Numericko resenje sistema diferencijalnih jednacina --1--------11------+-- ј+1/2 г------ ------, Щ,k : i+1,j,k . ;. . ' ( Up~+1/2,j,k i L----- - ______ Ј -+-------+-----~- ј-1/2 YLi-1/2 х i+1/2 i+З/2 Slika. 5.2 Tipicna kontrolna zapremina (predstavljena isprekidanom linijom) koja se koristi za integraciju jednacine odrzanja kolicine kretanja. Brzina faze (up)i+112J,k (tecna faza: p=l , para: р=2) se racuna za tacku (i+ 1/2j,k). Jednacine odтZanja entalpije za tecnu i za parnu fazu (poglavlje 4) se diskretizuju u sledecem oЫiku: о о [(h ) · . -(ho) . . Ј (apPp)i,j,kЛx;ЛyjЛzk + р 1,1,k р 1,1,k Лt { (арррир );+112,j,k [ (hp)i+l/2,j,k - (hp)i,j,k] + (аррри р );-112,j,k [ (hp)i,j,k - (hp)i-112,j,k ]} Лу jЛzk + {(а рРр v р );,j+l/2,k [ (hp )i,j+l f 2,k - (hp );,j,k] + (а рРр v р );,j-112,k [ (hp );,j,k - (hp)i,j-112,k ]} ЛzkЛх; + { (арРр wp)i,j,k+l /2 [ (hp );,j,k+l / 2 - (hp)i,j,k] + (а рРр wp )i,j,k-112 [ (hp );,j,k -(hp)i,j,k-112 ]} Лх;Лу ј = (Se,p)i,j,kЛx;ЛyjЛzk (5.6) gde је р=Ј za tecnu i р=2 za pamu fazu. Izvorni clan S e,p na desnoj strani jednacine (5.6) se definise preko: (5.7а) (5.7Ь) Posle ubacivanja izraza 5.7 i sredivanja, jednacina (5.6) postaje: 48 5. Numericko reSenje sistema diferencijalnih jednaeina ai,j,k(hp)i,J,k = a;+l/2,J,k (hp);+l,J,k +a;-112,J,k(hp ) i-I,J,k + Щ,j+l/ 2,k (hp)i,j+I,k +ai,j-112,k (hp)i,j - l,k + (5.8) Щ,ј,k+11 2 (hp)i,j,k+I + Q;,j,k- 112 (hp ) i,j,k-l + bi,j,k gdesu а,· 1· k = ai+l l 2 1· k + ai- 112 1· k + ai 1·+112 k + ,, '' )) ' ' о а·· 112 k+a · ·k 1 12+ а · · k 1 12+а · ·k 1,1- ' 1,ј, + l ,j, - l ,j, Ь· · k = (S ) · · kЛх·Л11 ·Лzk +а9 · k (h0 ) . · k 1,1, е,р 1,1, 1 :.r 1 1,1, р z,J, (5 .9а) (5.9Ь) (5 .9с) (5.9d) (5 .9е) (5.9±) (5 .9g) (5 .9h) (5 .9i) Entalpija, zapreminski udeo faze i gustina и kontrolnoj zapremini su definisani na sledeci nacin: h . - {(hp)i,j,k• if (up )i+l/2,j,k >О ( p)z+l / 2,j,k - (h ) · . if ( ) о р 1+l,1,k• ир i+l /2,j,k < · (5.lOa) { (ap)i,j,k• if (up)i+l /2,j,k >О (ap)i+l/2,j,k = ( ) if ар i+l,j,k • (up)i+l/2,j,k < О· (5.1 оь) 49 5. Numericko resenje sistema diferenciialnih jednacina { (Pp )i,J,k• if (up)i+l/2,J,k > О (pp )i+l /2,j,k = ( ) if ( ) о. Рр i+l ,J,k • ир i+l / 2,J,k < (5.lOc) Jednacina odrZaпja kolicine kretanja tecne i рате faze su diskretizovane na nacin kako је to predstavlj eno na slici 5.2. Na primer, diskretizovan oЫik jednacine odrZanja kolicine kretanja u х pravcu za kontrolnu zapreminu se defшise: [( ) ( )о Ј (a~p~)i+Il 2,J,k Лxi+1 12ЛyJдzk ир i+l /2,j ,k - ир i+1! 2,j ,k Лt + {(а рРрИ р )i+l,J,k [{ир ) ;+1,J,k - (и p )i+l/ 2,J,k Ј + (а рр ри p)i,J,k [(ир );+l/2,J,k - (up)i,J,k ]} Лу Jдzk + { (арРр vp )i+ll2,j +1/ 2,k [ (up)i+1/ 2,J+1/ 2,k -(up)i+l /2,J,k] + (а рРр vp)i+l /2,j - 112,k [(ир );+112,J,k - (up)i+l/2,J-112,k ] } ЛzkЛxi+l / 2 + { (арРр wp );+l / 2,j ,k+l / 2 [(и p )i+l /2,J,k+l /2 -(up )i+l / 2,J,k ] + (арРр wp)i+l /2,j,k- 112 [ (up )i+l f 2,j ,k -(и p )i+l /2,J,k- 112 ] } Лхi+112ЛУ Ј = -(ар );+\/ 2,J,k (Pi+l ,J,k - Pi,J,k )Лу JЛzk + (Sт,р );+112,J,k Лх;+112ЛУ Jдzk (5.11) gde se dimenzija dui х ose racuna kao Лx;+uz = xi - х;_1 • Izvorni clan za tecnu fazu је: (5 . 12а) а za pamu fazu: (5 . 1 2Ь) Parametri za odgovarajucu tacku kontrolne zaprernine (xi+112"yj, zk ) se racunaju kao srednja vrednost kontrolne zapremine, na primer: ( ) (ap)i+l ,J,k (x;+1- X;+112 ) + (ap )i,J,k(x;+112 - х;) ар i+l /2,J,k =----------------- xi+l -Х; ( ) (Pp )i+l ,j,k(X;+1- Xi+I /2 ) +(p p )i,j,k(x;+112 - х;) Рр i+ll 2,J,k = --''--_...;.;'-'-- --------=--____c:;_--- -- x;+l -Х; (5. 1 3а) (5. 1 3Ь) Komponente brzine na granici kontrolne zaprernine, za koje ne vafe izrazi za unutracnju tacku kontrolne zapremine, na primer u-brzina za tacku (х;+1, уј, z,J, se racuna kao srednja vrednost komponente susedne brzine, tj.: 50 5. Numericko resenje sistema diferencijalnih jednacina (5.14) Maseni fluks na granici kontrolne zapremine koja prolazi kroz reprezentativnu unutra5nju tacku kontrolne zapremine se racuna ро sledecoj serni: (5.15) (5.lба) (5 . lбЬ) Posle sredivanjajednacina (5.11) se pretvara u sledeci oЫik: а· 112 ·k(u ) · 1/ 2·k-a·1 ·k(u )· 1 ·k+ a· ·k(u )· ·k+ z+ ,ј, р L+ ,ј , - l+ ,ј , р L+ ,], L,j, р l ,j, ai+l / 2,J+l / 2,k (ир )i+l / 2,J+l,k + ai+I / 2,j- 112,k (u р )i+l/2,J-1,k + ~+112,J,k+l/2 (up)i+I / 2,J,k+I + ~+112,J,k-l / 2 (up)i+l/2,J,k- 1 + bi+l / 2,j,k gdesu ai+l/2,J,k = ai+l,J,k + ai,J,k + ai+l /2,}+112,k + ai+l / 2,J-112,k + о ai+l 12,J,k+l / 2 + ai+I I 2,J,k-112 + ai+112,J,k (5.17) (5. 1 8а) (5. 1 8Ь) (5.18с) (5.18d) (5.18е) (5.18f) 51 5. Numericko resenje sistema diferencijalnih jednacina (5.18g) (5.18h) b;+l / 2,J,k = - (ap)i+l / 2,J,k(P;+i,J,k - Pi,J,k)ЛyJЛzk + (Sт,р );+112,j,k Лх;+112Лу )Лzk + aP+l f 2,J,k (u~)i+l / 2,J,k (5.18i) Parcijalna diferencijalna jednacina za kolicinu kretanja sadrZi izvomi clan gradijent pritiska, koji је potrebno odrediti. То se postize algoritmom za resavanje strujnog polja, poznatim pod nazivom SIМPLE (Semi Implicit Method for Pressure Linked Equations), (Patankar, 1980, [34]), uzimajuCi и obzir prisustvo dve faze-tecnu i parnu. Pritisak se racuna na osnovu jednacine za korekciju pritiska. Resenje ove jednacine se dоЫја iterativnim postupkom. Komponenta brzine za n-tu iteraciju se racuna: /n) ::::: v(n- 1) + v' w(n) = w<п-1) + w' (5.19) а pritisak: (5.20) Komponente za korekciju brzinetions и: v: w' i korekciju pritiska su povezane preko smanjene forme и diferencijalnoj jednacini odr:Zanja kolicine kretanja, na primer jednaeina (5.17) se pise u sledecem obliku: 1 1 ' ai+112,j,k(up)i+l / 2,j,k ::::: - (ap)i+ll2,J,k(P;+1,J,k - Pi,J,k)ЛyJЛzk (5.21) ili (5.22) (5.23) Komponenta brzine (up)~=~/2,J,k za п-tи itemciju se racuna na osnovu vrednosti za prethodnu (п-1) -и iteraciju: ( )(п) _ (n-1) • • ир i+l / 2,J,k - (up)i+l / 2,J,k - (dp)i+l / 2,j,k (P;+1,j,k - Pi,j,k) (5.24) Ovo pokazuje da se brzina za (n-1)-u iteraciju koriguje u odnosu na korekciju pritiska da bi se doЬila komponenta brzine za п-tи iteraciju. 52 5. Numericko resenje sistema diferencijalnih jednacina I.zrazi za korekciju komponenti brzina u ostalim pravcima mogu da budu napisani u sledecem oЫiku: (5.25) (5.26) Zamenom izraza za brzinu u jednacini (5.20) u sumu jednacina odrLanja mase za tecnu i pamu fazu Gednacina 5.1 i odgovarajuca jednacina za parnu fazu - indeks р=2), jednacina za korekciju pritiska se izvodi na sledeci nacin: 1 ' 1 ' а· · kP· · k =а· 1 · kP· 1 · k +а· 1 · kP· 1 · k +а · · 1 kP· · 1 k + а · · 1 kP · · 1 k t,J, 1,Ј , 1+ ,Ј, 1+ ,Ј, 1- ,1, 1- ,1, 1,1+ , 1,1+, 1,1-, 1,1-, • +а· ·k 1Р: ·k 1+а· ·k 1Р: ·k l+b· ·k 1,1, + 1,1, + 1,1, - z,J, - 1,1, (5.27) gde Ь· . k = МFWЛх·Лу .Лzk i,1, z Ј Лt [ ( {п-1) (п-1)) ( (п-1) (п-1)) Ј а1Р!Иi +а2Р2и2 i+l/2,j,k - а1Р!И1 +а2Р2И2 i-112,j,k ЛујЛzk - [ ( (п-1) (п-1)) ( (п-1) (п-1)) Ј Лz Лх a1P!V1 +a2P2V2 i,j+l/2,k - a1P!V1 +a2P2V2 i,j-112,k k ; - [ ( (п-1) (п-1) ) ( (п-1) (п-1)) Ј a1P!W1 +a2P2W2 i,j,k+l/2 - a1P!W1 +a2P2W2 i,j,k- 112 Лх;ЛУј (5.28) (5.29а) (5.29Ь) (5.29с) (5.29d) (5.29е) (5.29f) а · ·k = а· 1·k+a·1 ·k+ a· · tk+a· · tk+a· ·k 1 + а· ·k 1 1,1, 1+ ,1, 1- ,1, 1,1+' 1,1-' 1,1, + 1,1, - (5.29g) Resenje sistema algebarskih jednacina se doЬija metodom iteracUe. Unutrasnje iteracije se koriste za resavanje sistema algebarskih jednacina koje su doЬijene 53 5. Numericko resenje sistema diferencijaLnih jednacina prirnenom jednaCina (5.1), (5.6), (5.:.1) i~ (5.~7) na osnovu .kojih se .dobija polj~ zapreminskog udela tecne faze, enta1p1Je tecne i рате faze, brzшe tecne i pame faze 1 korekcija pritiska, respektivno. Konvergencija ovog procesa se doЬija na sledeci nacin: 1. Izracunava se zapreminski udeo tecne faze u kontro1noj zapremini pomocu jednacine (5.1). 2. Izracunava se polje enta1pije za tecnu i parnu fazu и kontro1noj zapremini na osnovujednacine (5.6). з. Izracunava se polje brzina tecne i рате faze и kontro1noj zapremini na osnovu (5. 11). 4. Izracunava se polje korekcije pritiska u kontro1noj zapremini na osnovu (5.27). 5. Izracunava se polje brzina na osnovu (5.24 - 5.26). 6. Izracunava se polje pritiska (5.20). 7. Ponavljaju se posupci iz tacaka 1-6 dok se maseni udeo odredenjednacinom (5.28) ne dostigne odgovarajucu tacnost u okviru kontrolne zapremine (za vreme konvergencije iterativnog postupka desna strana jednacine (5.28) dostize vrednost nula). 8. Каdа se uspostavi odrZanje mase, poveca se vreme, zameni se novo izracunata vrednost kao pocetna i ponovo izracunaju termofizicke osobine na osnovu novih vrednosti zavisno promenljive. 9. Proracun se nastavlja do kraja prelaznog procesa. Resenje sistema algebarskih jednacina se vr8i pomocu metode linija ро linija (line-by- line) i tri-dijagonalnog matricnog algoritma Tree-Diagonal-Matrix-Algorithm (ТDМА), koji su predstavljeni и Patankar (1980) i Tannehill ( 1997). 5.1. Polarno cilindricoi koordinatni sistem Za polarno cilindricni koordinatni sistem usvajaju se sledece koordinate: r, <р i z. Susedne celije se oznaeavaju sa: W (West)-E (East), S (South)-N (North), L (Low)-H (Њgh). 5. 1. 1. Koeficijenti u diskretizovanim jednacinama Diskretizovanajednacina dobijena iz generalizovane diferencijalne jednacine је data и sledecoj formi: gde su: аЕ = DeA( IPeel )+max. [-Р: ,О] aw = DwA( IPewl )+ max. [Fw 'о] aN =DпА( IPenl )+max. [- Fn , О] as = DsA( IPesl ) + max. [;: 'о] (5.30) (5 .31) (5.32) (5.33) (5.34) 54 5. Numericko resenje sistema diferencijalnih jednacina ат =D1A( IPe1I )+max [- F, ,О] ао= DьА( IPeьl )+ max [Fь 'о] 0 р~ЛV а = Р Лt Ь = ScЛV + а~Ф~ ~ =~+~+~+~+~+~+~-~лv Slika 5.3 Polarno-cilindricne koordinate Konvektivni i difuzni clanovi se sada racunaju: F, = (ри )1 А,, D = F,A, , (Diff), F. ( ) А D = F;,Аь ь = ри ь ь• ь (Diff)ь Ае = rеЛ&Лz, Ап = ЛrЛz, А1 = 0,5 ·(t;, + rw ) Л&Лr (DifJ)e = (бr)е, (DifJ)n =rп (б&)п , (DifJ)1 =(бу)1 (5.35) (5.36) (5.37) (5.38) (5.39) (5.40) (5.41) (5.42) (5.43) (5.44) (5.45) (5.46) (5.47) 55 5. Numericko resenje sistema diferencijalnih jednacina 5. 2. Primenjene numericke mreze Sistem Ьilansnih jednacina se resava koriScenjem metode kontrolnih zapremina. sIМPLE (Semi Implicit Method for Pressure-Linked Equations) algoritam se koristi za resavanje jednacine za korekciju pritiska iz jednacina odr"Lanja mase i kolicine kretanja. Tridimenzionalno strujno polje је disk:retizovano u Kartezijanskim koordinatama. Na slici 5.7 је data trodimenzionalna uniforrnna mreZa, kojaje koriscena pri simulaciji k:rize razmene toplote и bazenskom kljucanju. Numericka mreia se sastoji od 40х80х40 kontrolnih zapremina. Donji zid koji se zagreva је podeljen u 4 zone, i svaka zona ima razlicita mesta na kojima se generisu mehurovi, definisani koriscenjem metode slucajnih brojeva. " " ,, ' о. 07 " ·' о. 06 " ,, 1 " 05 ' '· " . 1 1 04 ; 1 1 р 03 о. ocetna pozicija nivoa tecnosti m :: . 1 о. 1 Zo 01 о. ne u kojima se generise parna faza ' ' Zid о 0.01 0.02 0.03 0.04 (40х40х40) Sirina(m) z Slika 5.4 Numericka mre:Za za slucaj bazenskog kljucanja Za vertikalno strujanje na gore и cevi eetvrtastog poprecnog preseka tridimenzionalno strujno polje је diskretizovano и Cartesianskim koordinatarna. Na slici 5.8 је data trodimenzionalna uniformna mreZa, koja је koriscena и ovom modelu. Numericka mre:Za se sastoji od 40х60х40 kontrolnih zapremina. Bocni zid se zagreva i podeljen је и 4 zone, i и svakoj zoni se generisu mehurovi na mestima koji su definisani koriseenjem metode slucajnih brojeva. Zagrejacki zid је predstavljen pomo6u sloja koji se sastoji od l о celija. 56 5. Numericko resenje sistema diferencijalnih jednacina х ~· у Slika 5.5 Numericka mreza za slucaj kljucanja u cetvrtastoj cevi vertikalno na gore Slika 5.6 Numericka mreza za slucaj kljucanja u cevi okruglog poprecnog preseka vertikalno na gore Za slucaj verikalne cevi okruglog poprecnog preseka strujno ро\је је diskretizovano u polarno-cilindricnim koordinatama. Na slici 5.6 је data trodimenzionalna mrel.a, kojaje koriscena и ovom modelu. Simulacija је radena za cevni isecak. Numericka mreZз se sastoji od 40х40х40 kontrolnih zapremina razlicitih dimenzija. Zid koji se zagreva је podeljen u 4 zone, i svaka zona ima razliCita mesta na kojiшa se generisu mehurovi, а koja su definisana koriscenjern metode slucajnih brojeva. 57 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote 6. Prikaz i analiza dobljenih rezultata simulacije krize razmene toplote Jedan od glavnih proЫema и konstrukciji i analizi sigumosti razlicitih vrsta generatora pare i razmenjivaca toplote u energetskim, termotehnickim ili procesnim postrojenjima је kriza razrnene toplote pri kojoj moze doci do ostecenja - pregorevanja zagrejackog zida optere6enog visokirn toplotnim fluksom usled naglog pada koeficijenta prelaza toplote pri mehurastom kljueanju. Кlasican CFD pristup ne moze и potpunosti da objas ni kompleksan proces kljucanja u uslovima velikih vrednosti toplotnog fluksa. Ovaj fenomen zahteva matematicko modeliranje i numericku simulaciju dvofazne mesavine na makro nivou i takode modeliranje i simulaciju na mikro nivou rasta mehura i napustanja zagrejacke povrsine. U ovom radu је razvijen trodimenzijski numericki model za simtllaciju procesa kljucanja u horizontalnoj posudi kvadratnog poprecnog preseka (bazensko kljucanje), u vertikalnoj cevi kvadratnog poprecnog preseka kroz koju struji dvofazna mesavina vertikalno navise i и cilindricnoj vertikalnoj cevi, kroz koju struji dvofazna mesavina vertikalno navise. Dvofazno strujanje је modelirano pomocu modela dva fluida, tj. model se sastoji od jednacina od.Гanja mase, kolicine kretanja i energije za svaku od faza pojedinacno, kao i dodatnih konstitutivnih korelacija kojima se opisuje medufazno delovanje. Takode, proracunava se i nestacionarno temperatursko polje u zagrejackom zidu sa uniforrnnim unutra5nj im izvorom toplote. Da bi se pojava kriticnog toplotnog fluksa mogla predvideti na odgovarajuci nacin, bilo је neophodno uzeti u obzir i rnikro nivo, tj . predvidanje mesta klijalista mehurova, rast mehura i odvajanje od zagrejacke povrsine. Korisceni su odgovarajuCi parametri pomocu kojih su opisani kljucni parametri za generaciju mehurova na zagrejackoj povrsini, kao sto su gustina centara nukleacije mehurova, vreme rasta mehurova na zagrejackoj povrsini i odredeni nivo slucajnosti polozaja mehura и okviru odgovarajuce zone. Toplotni fluks sa zagrejanog zida па dvofaznu mesaviпu se rasporeduje neuniformno па zagrejackoj povrsini sa pikovima na mestima gde se geпerise mehur. Polozaj mehura и odgovarajucoj zoni se defшise pomocu funkcije slucajnih brojeva, dok broj zona и kojima se mehur stvara zavisi od karakteristike materijala i hrapavosti zagrejacke povrsine. Pomocu ovog modela је moguce modelirati i glatke i hrapave povrsine. Nwnericka simulacija је sprovedena za razlicite vrednosti toplotnog fluksa, ali dovoUno velike i Ыiske pojavi zasusenja i nastanku kriticпog toplotnog fluksa tako da је Ьilo moguce dovoljno tacno odrediti njegove vrednosti. Rezultati numericke simulacije koji se odnose na bazensko kljucanje su uporedeni sa odgovarajucim eksperimentalnim rezultatima, koji su dostupni u literaturi i postignuto је zadovoljavajuce slaganje izmedu eksperimentalnih rezultata i rezultata modeliranja. 58 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote Mehurasto kljucanje i nastanak krize razmene toplote su numericki simulirani u dvofaznoj mesavini vode i vodene pare, а zagrejacki zid је od bakra. Polinomima odgovarajuceg stepena su aproksimirani specificni toplotni kapacitet i toplotna provodnost bakra u zavisnosti od temperature zida. Pretpostavljeno је da se dvofazna mesavina vode i vodene pare nalazi na liniji zasicenja. 6.1. Bazensko kljucanje Najpre је dat primer bazenskog kljuCa.nja и kome је zagrejacki zid predstavljen samo jednim slojem celija [35]. То znaci da nema raspodele temperatura ро zidu, vec se samo detektuje temperaturski skok i na taj nacin se potvrduje pojava kriticnog toplotnog fluksa. Simulacija је radena za veliku vrednost toplotnog fluksa od 1 ОО W/cm2. Model predstavljen и poglavlju 4 је resavan numericki. Predvideno је ponasanje zagrejacke povrsine za vreme procesa krize kljucanja i pokazan је nastanak krize kljucanja. Na slici 6.1 је prikazana raspodela zapreminskog udela pame faze na pocetku procesa kljucanja, dok је na slici 6.2 prikazana raspodela zapreminskog udela parne faze za period razvijenog kljucanja. Poredenj em slika 6.1 i 6.2 se moze primetiti povecanje nivoa dvofazne mesavine i povecanje zapreminskog udela рате faze. Takode se na slici 6.2 vidi formiranje vecih parnih mehurova, kao posledica intenzivnog isparavanja za visoke vrednosti toplotnog fluksa od 1 ОО W /cm2, а takode i zbog spajanja nekoliko manjih mehurova u jedan veci. Veci delovi sa parnim mehurovima i visi nivo dvofazne mesavine se formira u slueaju viseg toplotnog fluksa. U slucaju nizih vrednosti toplotnog fluksa nivo dvofazne mesavine је skoro ravan. 0.07 0.06 0.05 -Е -~ 0.04 FI 0.94 .(ј) 0.77 5 0.6 О.ОЗ 0.43 0.26 0.09 0.02 0.01 0.02 0.04 0.06 0.08 Sirina (m) Slika 6. l Zapreminski udeo рате faze na pocetku procesa kljucanja ( top lotni fluks 1 ОО W/cm2, gustina nukleacije 1 cm-2, vreme rasta mehura 5·10- s) 59 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote u slucaju nizih vrednosti toplotnog fluksa manja povrsina zagrejackog zida је prekrivena parom, dok se sa povecanjem fluksa ta povrsina uvecava. U slucaju vecih vrednosti toplotnog fluksa se formiraju vece zaprernine pare и dvofaznoj mesavini zbog intenzivne generacije parne faze i zbog spajanja manjih mehurova и ve6e. Formiranje velikih mehurova sprecava obnavUanje zagrejacke povrsine sa tecnos6u sto dovodi do us]ova zasusenja. Е' ........ ~0.04 :;· Q) 0.07 0.06 0.05 ~ 0.04 'iЂ > О.ОЗ 0.02 0.01 0.02 у z~ х FI 0.79504 0.59628 0.39752 0.19876 FI 0.864741 0.807091 0.749442 0.691792 0.634143 0.576494 0.518844 0.461195 0.403546 0.345896 0.288247 0.230597 0.172948 0.115299 0.0576494 0.04 0.06 Sirina (m) 0.08 Slika 62 Zapreminski udeo рате faze za razvijeno k1jucanje za topJotni fluks od 1 ОО W/cm2 gustina nukleacye 1 cm·2, vreme rasta mehura 5·10·3 s, 0,15 s od pocetka grejanja) 60 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmeпe toplote Na slici 6.3 је prikazano temperatursko polje и toku procesa kljucanja. Za vreme procesa isparavanja zapreminski udeo parne faze na zagrejackom zidu se poveeava, sto dovodi do pojave zasиSeпja i роvеСапја temperature zida па delovima zagrejacke povcline. Prikazana је raspodela temperatura za tri razlicita vremenska treпutka za vreme procesa kljueanja. Slika sa najmanjim delovima sa visokom temperaturom odgovara poeetku procesa, dok veci udeo delova sa visokom temperaturom odgovara kasnijem periodu krize kljucanja. Temperatura dvofazne me5avine је konstantna i jednaka temperaturi saturacije, kako је i prikazano na slici 6.4. Temperaturske promene na jednoj lokaciji na zagrejackom zidu gde se javlja kriza kljucanja su prikazane na slici 6.4. Na pocetku procesa kljucanja temperatura zida је priЫiZno jednaka temperaturi zasicene tecnosti. Kada mehurovi pocinju da rastu i kada se pojavi zasuseпj e zida zbog formiranja velikih mehurova, temperatura zida se poveca za oko 50 К. Kako se proces kljucanja nastavi, delovi na kojima se javlja zasusenje, tj. delovi па kojima se javlja temperaturski skok, menjaju svoja mesta na zagrejackom zidu i povrsine koje oni zauzimaju se povecavaju. Toplotni fluks koji dovodi do zasusenja zida i пaglog temperaturskog skoka је kritican toplotni fluks (СНF). 0.07 О.Об 005 :[ :! 0.04 "iii 5 0.03 0 .02 0.01 о с ст ~·0.02 3 ....... Sirina (m) Т1 i 4 17.911 406.681 395.451 384.221 0.06 Т1 • 417.897 406.652 395.407 384.163 0.08 о с ст о с: О' ~-0.021-F-~-fjf-,,;--:tl ]: Sirina {m) ~- О.02~51'Е--"'.5~~~ ~ Sirina (m) z х_Ј Т1 417,899 406.656 395.414 384.171 z x_d Т1 • 417.91 1 406.681 395 .451 384.221 Slika 6.3 Temperatursko polje pri pojavi zasuseпja u mehurastom kljucanju za toplotni fluks od 1 ОО W /cm2 (О, 14 s-gore levo, О, 12 s-gore desno, О, 14 s -dole levo i posle О, 15 s-dole desno) 61 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote 430 +-* 420 41 0 -~ -1- 400 390 380 * 370 о.оо 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 t(s) Slika 6.4 Temperatura zagrejackog zida ujednoj karakteristicnoj celiji и vremenskom periodu Slede rezultati simulacija dobijeni sa konaenom deЫjinom zagrejackog zida od 2 mm [36]. DеЫј ini zida је pridruZena numericka mreZз sa 1 О slojeva ро visini. Simulacija Је Vl'Sena za eetiri razlicite vrednosti toplotnog fluksa, za 10 W/cm2, 40 W/cm2, 70 W/cm i 85 W/cm2• U zavisnosti od vrednosti toplotnog fluksa kome је izlofena zagrejaCka rmina koristi se i razlicito vreme odvЩanja mehura od zagrejacke poVI'Sine i to 2,s ·1 о· s za 10 W/cm2, 10·2s za 40 W/cm2, 7,0·10·3s za 70 W/cm2 i 5,0·10·3s za 85 W/cm2 (slika 4.2). l u ovom slucaju је predvideno ponaSaпje zagrejacke povrsine sa pojavom krize kJjucanja. Rezultati sumulacije koji se odnose na dvofaznu mesavinu su isti kao i и prethodnom slucaju. Na slikama 6.5, 6.6, 6.7, 6.8 i 6.9 је prikazan zapreminski udeo рате, odnosno tecne faze za razvijeno kljucanje za razlicite vrednosti toplotnog fluksa u razliCitim vremenskirn trenucima. Primeceno је povecanje nivoa dvofazne mesavine i povecanje zapreminskog udela рате faze u dvofaznoj mesavini sa povecanjem vrednosti toplotnog tluksa. Vidi se formiranje vecih раmЉ mehurova na zagrejackom zidu zbog intenzivnog procesa isparavanja, zbog visoke vrednosti toplotnog fluksa, kao i zbog spajanja nekoliko parnih mehurova u jedan veci. Velicina mehurova zavisi i od vrednosti toplotnog fluksa, ~· vece su zapremine sa pamom fazom za veee vrednosti toplotnog t1uksa. Za vreme procesa isparavanja, zapreminski udeo parne faze na zagrejackom zidu se poveeava sto dovodi do pojave zasиSenja zagrejaekog zida i poveeanja ternperature zida. Temperatura dvofazne me5avine је skoro konstantna tokom procesa isparavanja i jednaka је temperaturi zasieenja. 62 6. Prikaz i analiza doЬijenih rezultata simulaci je krize razmene toplote у 0.09 z~ 0.08 х 0.о7 FI 0.06 ! 0.8 е 0.6 ~ О.05 0.4 "' с0.2 ~ 0.04 о.оз 0.02 0.01 о 0.025 0.05 0.075 0.1 Sirina (m) у ':Ф 0.09 х 0.06 FI g: O.Q7 0823438 08818'/S FI 0800313 0.06 о 7387$ ~ 0.8 0077188 0.6 :[ 0.05 001~ 0.4 0!5640е3 :! 0 •92S 0.2 ~ 004 о •ЭОS138 0388375 о 307813 0.03 о 2•625 0 18'1688 о 12312$ 0.02 00015025 0.01 о 0.05 Slrlne(m) Slika 6.5 Zapreminski udeo parne faze za razvijeno kljueanje za toflotni fluk.s 10 W/cm2 za razlicite vremenske periode 0,1 s i 2 s (gustina nukleacije 1 cm·, vreme rasta mehura 2,s·10·2 s) Na slici 6.6 su dati temperaturski profili na poVI'Sini zagrejacke ploce na strani koja је и dodiru sa dvofaznom meSa.vinom, kao i profil temperatura и poprecnom preseku ploce i dvofazne me5avine za razlicite vremenske intervale za toploni fluks od 1 О W/cm2 • 63 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote Sirina (m) Sirina (m) z х _Ј Т1 380 'Л7.6fЈ1 '315.333 'ЛЗ Т1 390 384.333 378.667 373 Slika 6.6 Temperatursko polje za toplotni fluks 10 W/cm2 za razlicite vremenske periode 0,1 s i 2 s (gustina nukleacije 1 cm-2, vreme rasta mehura 2,s·10-2 s) Na slici 6.7 је prikazan zapreminski udeo pare и razlicitim presecima dиZ posude sa kljucalom dvofaznom mesavinom u razlicitim vremenskim trenucima. Slike u levoj koloni prikazuju zapreminski udeo pare u vertikalnom preseku na polovini rastojanja izmedu prednjeg i zadnjeg zida posude, а slike u desnoj koloni isecke u trodimenzionalnoj zapremini dvofazne mesavine. Rezultati pokazuju da и pocetnom periodu zagrevanja (О, 1 s nakon poeetka zagrevanja) mehurovi nastaju na zagrejanoj povrsini, ali njihova zapremina је manja i povclina dvofazne mesavine је jos uvek ravna. Sa daljim tokom procesa, и 0,25 s, zapreminski udeo pare је poveCa.n, а povclina dvofazne meSa.vine postaje ta.lasasta usled parnih mehurova koji isplivavaju iz zapremine dvofazne mesavine и parni prostor. U 0,8 s proces kljueanja је veoma intenzivan, povrsina dvofazne mesavine је razЬijena, а delovi tecnosti sa parnom fazom su noseni skoro do vrha posude. Istovremeno znacajan deo zagrejacke povrsine је prekriven parom. Rezultati numerickih simulacija procesa kljucanja za jos veci toplotni fluks od 70 W /cm2 su prikazani na slici 6.8 i 6.9. Rezultati zapreminskog udela pare na slici 6.8 ukazuju na izrazito dinamican proces kljueanja sa mlazevima pare koja istice iz posude. Na slici 6.9 је prikazano temperatursko polje na povтSini zagrejacke povrsine u kontaktu sa dvofaznom meSa.vinom. Veci deo povтSine grejaea је za skoro 1 ОО К pregrejan u odnosu na temperaturu saturacije od 373 К. 64 0 09 0 08 0.07 0 06 g 0.05 .., с ~ 0.04 о.оз 0.02 0.01 о 0.09 0 .08 0 .07 0 .06 §: 0 .05 <11 с ~ 004 0 03 0 .02 O.Q1 о 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote 0.025 0 .025 F1 0.836563 0.874125 0.811688 о 7.&925 0.886813 0.824375 ом1т 0.4995 0.437063 0.374625 0.312188 024975 0.187313 О.12•875 0.0624375 0.05 0.о75 Sirina (m) FI 0.778125 0.72625 0674375 0.8225 0.570825 0.51875 О."66875 0 .415 о 363125 о 31125 0259375 02075 0.155625 0.10375 0.051875 0 .05 0 .075 Sirina (m) FI 0.841875 0.78575 0.729625 0.6735 0.617375 0.56125 0.505125 0.449 0392875 0.33675 0280825 0.2245 1 0.168375 0.11225 0.056125 0 .05 Sirina (m) 0.075 у .08 z~ .06 х g .04~ FI ёЂ ! о.е 5 0.6 04 0.2 0 .1 .08 у .ов z~ х g .04:!! 'ii А 5 0.7992 05994 0.3996 .02 0.1996 0 .1 у .08 z~ х 0 .1 Slika 6.7 Zapreminski udeo рате faze za razvijeno kljucanje za toplotni fluks 40 W/cm2 za razlicite vremenske periode 0,1 s, 0,25 s, i 0,8 s (gustina nukleacije 1 cm·2, vreme rasta mehura 10·2 s) 65 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote 0 .09 у 0.08 z~ 0.07 FI 0-375 0-75 х О.Об 0.749125 08915 §: е о.еззе75 ~О.05 057625 FI "' о 518е25 ! 0.8 с 0~1 0.6 ~ 0.04 о 403375 0.4 о 34575 0.2 о 288125 0.03 02305 0.172875 0.02 о 11525 о 057625 0.01 о 0.025 0.05 0.075 0.1 Sirina (m) у .08 z~ 0.09 0.08 FI х 0.07 0816438 0783875 о 701i313 §: О.Об 0.65475 FI 0600188 е 0.545825 ! 0.8 ~О.05 0.491083 0.6 "' 0 • 385 с 0381838 0.4 ~ 0.04 0.2 о 327375 0272813 О.ОЗ 0.21825 0.163$88 0.109125 0.02 0.0545825 0.01 о 0.05 0.075 0.1 Sirina (m) Slika 6.8 Zapreminski udeo рате faze za razvUeno kljucanje za toplotni fluks 70 W /cm2 za razlicite vremenske periode 0,1 s i 0,7 s (gustina nukleacije 1 cm-2, vreme rasta mehura 7' 10-3 s) 66 6. Prikaz i analiza doЬijenih rezultata sirnulacij e krize razmene toplote о с: СЈ ~· 0.02 3 - Sirina (m) z x__d Т1 450 424.333 398.667 373 Slika 6.9 Temperatursko polje za toplotni fluks 70 W/cm2 и 0,7 s (gustina nukleacije 1 ст -2, vreme rasta mehura 7' 10·3 s) z z x_d x__d 0.04 0.04 о с ll о ; 0.02 с ст 3 Т1 ~·0.02 v 378 ]: о 376.333 374.667 0.04 0 .02 о 373 Slrina(m) Sirina (m) t=O, ls t=0,7s 500 -tr zona 1 475 -+- zona 2 - -в- zona З ~ 450 "'*"" zona 4 n:s ..... * zona 5 ::::i ~ 425 -+-zona 6 ф а. -+- zona 7 ~ 400 - zona 8 1- 375 350 ---1 о 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Vreme (s) Slika 6.10 Promena temperature zida tokom vremena za osam karakteristicnih celija na gomjoj povr~ini ploce za п= 1 cm-2 i q=85 W /cm2(gore t=O, 1 s i t=O, 7s) 67 6. Prikaz i analiza dobljenih rezultata simulacije krize razmene toplote Pri toplotnom fluku od 85 W/cm2 dolazi do naglog temperaturskog skoka и zagrejackom zidu, predstavljeno slici 6.1 О, sto odgovara pojavi krize kljucanja. U taЬeli 6.1 su uporedno prikazani rezultata numericke simulacije za razlicite vrednosti toplotnog fluksa. Data је zavisnost srednje temperatura zida i srednjeg koeficijenta prelaza toplote u funkciji od vrednosti toplotnog fluksa. Кriza razmene toplote se javlja za vrednost toplotnog fluksa od priblimo 85 W/cm2• Rezultati iz tabele 6.1 su prikazani na slikama od 6.11 do 6.14. Dati su dijagrami zavisnosti srednje temperature zagrejacke ploce i srednjeg koe:ficijenta prelaza toplota u zavisnosti od toplotnog fluksa, kao i zavisnosti srednje temperature zagrejacke ploce i srednjeg koeficijenta prelaza toplota u zavisnosti od razlike srednje ternperature zida i ternperature saturacije za konstantnu vrednost gustine nukleacije od n=l cm·2. Slika 6.15 predstavlja krivu kljucanja. Као sto se moze primetiti sa povecanjem vrednosti toplotnog fluksa se povecava i srednja temperatura zagrejacke ploce i srednji koeficijent prelaza toplote. Tabela 6.1 Srednja temperatura zida i srednji koeficijent r,relza toplote za razlicite vrednosti toplotnog fluksa (za gustinu nukleacije od n=lcm· ) toplotni fluks srednja temperatura srednji koeficijent prelaza q (W/cm2) zida Twsr (К) toplote hsr (W/m2K) 10 391,1 5524,86 40 405,2 12422,36 70 410 18918,92 85 411,6 22020,73 1-- в l 420 410 390 380 +----.~--,-~.....--...,.-~~-т----.г----т-~.....--..., о 20 40 60 80 100 Slika 6.11 Zavisnost srednje temperature zagrejacke ploce od toplotnog fluksa 68 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote 24000 20000 16000 5Z "' .е з: 12000 .._; .t::. 8000 4000-t--~..------..~-.-~-.-~..---.-~-.-~.....---.~-. о 20 40 60 80 100 q(W/cm2) Slika 6.12 Zavisnost srednjeg koeficijenta prelaza toplote na zagrejackoj ploci od toplotnog fluksa l -11-дl 10 20 зо 40 50 Slika 6.13 Zavisnost toplotnog fluksa od д Т sat 69 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote Q" N .€ 3:: -;; .s= 1---cl 24000 20000 16000 12000 8000 4000-<--~~~~~~~~~~~~~~~~ 10 20 30 дТ нt(К) 40 50 Slika 6.14 Zavisnost srednjeg koeficijenta prelaza toplote hsr od ЛТsа1 Do sada su prikazani rezultati za gustinu nukleacije od п=1 cm-2• U daljem radu се Ьiti analiziran uticaj hrapavosti povrsine na pojavu kriticпog toplotnog fluksa, kao i па promenu srednje vrednosti temperature zida zagrejaca i srednje vrednosti koeficijenta prelaza toplote. Na slici 6.15 је data promena temperature zida tokom vremeпa za osam karakteristicnih celija па gornjoj povrsini ploce za drugu hrapavost zagrejacke povrsine, za п=39 cm-2 i q=l 10W/cm2• Primecuje se veca temperatura zagrejecke povrsine, kao i vece povrsine zida sa poviSenom temperaturom za vecu hrapavost zagrejacke povrsine. Kada dode do naglog povecanja vrednosti temperature zida javlja se kriza razmene toplote. Toplotni fluks pri kome se ovo desava је kriticni toplotni fluks. Poredeцje sa odgovarajucim eksperimentalnim rezultatima је dato na slici 6.16 i zadovoljavajuce slaganje rezultata ovde predlozenog modela i eksperimentalnih rezultata је postigrшto. 70 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote о с: cr ~·0 .02 Sirina (m) z х _Ј ~ I~ 391 382 373 t=0,05s 500 т- 1 -h- zona 1 475 -- zona 2 1 -e- zona З ~...... 450 i ""*-- zona 4 --zona 5 о с: 0.04 ~ 0.02 Q) о 0 .04 0.02 Sirina (m) t=0,22s ~ 425 ј --zona 6 Ф -+-zona 7 ! 400 1 - zona8 --~-:;;Ј..--;Б;-....1~------ 375 350 -- 0 0.05 0.1 0.15 Vreme (s) 0.2 z x_d ~ ~ ~ 373 о 0.25 Slika 6.15 Promena temperature zida tokom vremena za osam karakteristicnih celija na gomjoj povrsini ploce za п=39 cm-2 iq=l10W/cm2(gore t=0,05s i t=0,22s) 71 е а С1 i- 11 :х: Е .§. !С с: ·џ; "!> 70 60 50 40 зо 20 10 о о 6. Prikaz i analiza doЫjeniЬ rezultata simulacije krize razmene toplote у 09 ь__х 0.8 0.7 0.6 ~ (Ј) 0.5 ::;· (\) ОА ~ о.з FI § nв 0.2 nб n• 0. 1 n2 о.о Х = 35 ПllП Sirina (m) \ 70 г ~" ео ,.. / / ~ ' ~ • 50 е "° §. " с ~ , :iii эо > 3 2.0 ~ 1111; ~ 10 --- о 0.2 0.4 0.6 0.8 1.0 о 02 О.• 0.6 0.8 zapreшinski 11deo ра111е faze zapreminski udeo parne faze Slika 6.16 Poredenje sa eksperimentalnim rezultatima Ovde predstavljani rezultati su pokazali da је moguee izvrsiti numericku simulaciju krize kljucanja i da је moguce predvideti temperaturski skok и zagrejackom zidu и slucaju da se kritiean toplotni tluks pojavi [ 18). Struktura povr8ine i karakteristike materijala od koga је izraden zid imaju znacajan uticaj na koeficijent prenosa toplote i vrednosti kriticnog toplotnog fl.uksa. Dobijeni rezultati се obezbediti korisne informacije koje su od znacaja za projektovanje i proizvodnju efikasne zagrejacke povrsine sa stanovista prenosa toplote, koje mogu dovesti do vecih vrednosti kriticnog toplotnog fl.uksa, odnosno obezbedenja bladenja zagrejacke povrsine bez pojave zasusenja na njoj. 72 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote 6.1.1. Uticaj obradenosti zagrejacke povrsine na pojavu kriticnog topotnog fluksa lspitivan је uticaj hrapavosti z.agrejacke povгSine na kritican toplotni fluks. Pokazano је daje obradenost povrsine и direktnoj vezi sa brojem centara nukleacije: 1 n=-b2 (6.1) Pri cemu је n-broj centara nukleacije, а Ь2 povrsina proracunske zone na zagrejackom zidu и kojoj se javlja mehur. Sprovedena је analiza uticaja obradenosti zagrejacke povrsine, па pregrevanje zida i prisustvo pare na povrsini zida. Numericke simulacije su sprovedene za razlicite dimenzije zona odgovarajuce gustine nukleacije, kao sto је prikazano и tabeli 6.2. ТаЬеlа 6.2 Dimenzije zona i odgovarajuce gustine nukleacije primenjene u analizi uticaja hrapavosti povrsine па vrednosti pregrejanja zagrejackog zida i dostitanje kriticnog toplotnog fluksa zona dimenzija Ь gustina nukleacije п broj zona primenjen u (cm) (cm-2) simulaciji I 0,3 11 4х4 п 0,2 25 4х4 IП 0.16 39 4х4 IV 0,14 51 4х4 DoЬijeni rezultati su prikazani и taЬelama 6.3 do 6.6. Tabela 6.3 Srednja temperatura zida zagrejaca i srednja vrednost koeficijenta prelaza toplote za razlicite vrednosti toplotnog fluksa za gustinu nukleacije n= l lcm·2 toplotni fluks (W/cm") Tsr ·~о с: а. ro ·~ > о ro ~ N U ..... ro "U·~ ro Ф N ._ О> Ф ro Е N ф ..... > 2.SOE-02 1\ 2 .ООЕ-02 ' ' ' 1 .50Е-02 1 .ООЕ-02 ~ 5.ООЕ-03 - О.ООЕ+ОО 200 ' . 1 -kontaktni ugao = 40 stepeni ~ - - kontaktni ugao = 90 stepeni -- ' r- , ...... 1-----;.....--- :--__ 1 --- --- - - 400 600 800 1000 1200 1400 1600 Toplotni fluks, kW/m2 Slika 6.21 Uticaj kontaktnog ugla Rezultati na slici 6.21 pokazuju da vreme Ьoravka mehura na mestu klijanja ne zavisi od ugla kva5enja pri vrednostima toplotnog fluksa veeim od priЫiZпo 100 W/cm2, tako da је numericka analiza uticaja ovog vremena sprovedena za toplotni fluks cije su vrednosti 10 W/cm2, 40 W/cm2 i 70 W/cm2 na pregrejanje zagrejackog zida (Тz-Тsш) i odgovarajuCi koeficijent prelaza toplote (6.3) 77 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote Dobijeni rezultati su prikazani u taЬeli 6.8 i na slikama 6.22 i 6.23. ТаЬеlа 6.7 Vrednosti vremena Ьoravka mehura na mestu njegovog rasta na zagrejackom zidu su sra.Cunata prema jednacini (4.7) za razlicite vrednosti toplotnog fluksa i ugla kva$enja q (W/cm2) t ь (s) kontaktni ugao 40° kontaktni ugao 65° kontaktni ugao 90° 10 1 75·10-2 ' 2 5·10-2 ' 3·10-2 40 8 3·10-3 ' 10-2 1,25·10-2 70 6 т10-з ' т10-з 1 5·10·3 ' 100 5·10·3 5·10·3 5·10·3 Nakon izvrsenih numerickih eksperimenata dob~eni su sledeci rezultati: ТаЬеlа 6.8 Srednja vrednost temperature zida zagrejaca i srednja vrednost koeficij enta prelaza toplote za razlicite vrednosti kontaktnog ugla q(W/cm2) 10 40 70 800 700 N- 600 Е ~ 500 .::.:. ....... U) -3 400 <;:: :s зоо о i5.. о 1- 200 100 о ·'-- -<--- 10 kontaktni ugao 40° kontaktni ugao 65° kontaktni ugao 90° Twsr(К) hsr(W/m2K) Twsr (K) hsr (W/m2K) Twsr (К) hsr (W/m2K) 391,1 5524,9 391 ,1 5524,9 391 ,9 5291,0 405,2 12422,4 405,9 12158,1 406,2 12048,2 410 18918,9 410,3 18766,8 410,6 18617,0 - - kontaktni ugao 40 stepeni -+- kontaktni ugao 65 stepeni -е- kontaktni ugao 90 stepeni д ~ , / л ~ -- 50 100 Temperatura zida - temperatura saturacije (К) Slika 6.22 Zavisnost srednje vrednosti toplotnog fluksa od razlike temperatura zida i temperatura saturacije za razlicite kontaktne uglove 78 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote ..-.. ::.::: "' ~ 17000 +---- ~ .--~~~~---'~~~, Ф ~ kontaktni ugao 40 stepeni1 о а. -- kontaktni ugao 65 stepeni · ~ 13000 N -&- kontaktni ugao 90 stepeni ro 1 ф .... а. ё .~ 9000 т-~~~~~~-г~~~/---1Г--~~r-~--r-~--r-~r--t--t-----i u ..,, ф о ::.::: 5000 +-~~~~~----"'"-'-~~~~~~~~~_._~_,_~~~~---1 10 50 Temperatura zida - temperatura saturacije (К) 100 Slika 6.23 Zavisnost srednjeg koeficijanta prelaza toplote od razlike temperatura zida i temperatura saturacije za razlicite kontaktne uglove Na osnovu rezultata na slikama 6.22 i 6.23 moze se zakljuCiti da kontaktni ugao ne utice znaeajno na temperaturu zida zagrejaca i koeficijent prelaza toplote. 6.1.3. Racunanje koeficijenta prelaza toplote kori.Scenjem razlicitih empirijskib izraza i poredenje sa sopstvenim modelom Kolicina toplote ро jedinici vremena koja se preda fluidu је jednaka: (6.4) Odavde se dobija izraz za koeficijent prelaza toplote: (6.5) Na slici 6.24 su dati rezultati numericke simulacije temperature gomje povr5ine zagrejackog zida na kojoj kljuca voda na pritisku od 1 bar i pri vrednosti toplotnog fluksa od 10 W/cm2. Za ovde prikazani primer osrednjavanjem temperature na povr5ini zagrejackog zida na slici 6.24 se doЬija premajednacini 6.5 srednja vrednost koeficijenta prelaza toplote koja iznosi h = 5524, 9 W/m2K. 79 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote о с ст ~· 0.02 -3 - Sirina (m) z х_Ј Т1 390 384.333 378.667 373 Slika 6.24 Temperatura zida koji se zagreva (gornja povrsina) U literaturi postoje razliciti empirijski izrazi koji se koriste za odredivanje koeficijenta prelaza toplote. Ti izrazi su dati za specificne uslove и kojima se odvija proces, kao i za razlicite vrste fluida i vrste materijala i stepena obradenosti zagrejacke р1осе. Ovde su navedeni samo neki od njih i moze se primetiti da se doЬijaju znacajno raz)icite vrednosti koeficijenta prleza toplote и zavisnosti od toga koji se izraz koristi. То је i oeekivano, jer uslovi и numerickom eksperimentu nisu identicni uslovima и eksperimentima pomocu kojih su doЬijeni ovi empirijski izrazi. Mostinski (17] i Starczewski [17] su predlozili sledecu korelaciju za racunanje koeficijenta prelaza toplote. Ovi izrazi ne uzimaju и obzir uticaj zagrejacke povrsine. !( ) ] 8 0.17 4 1.2 10 10 PR = . PR + PR + PR gde је Рс (bar) kritican pritisak, а PR = 1!_. Рс (6.6) Srednja vrednost koeficijenta prelaza toplote pri kljucanju vod.e za p=l bar doЬijena primenom ovog izrazaje h=5296,183 W/m2K. Foster-Zuber-ova korelacija : (6.7) za koeficijent pralaza toplote za dvofaznu mesavinu voda-vodena para pri pritisku p=lbar. Na osnovu nje se doЬija h=8999,572 W/m2K. Srednje vrednosti koeficijenta prelaza toplote za bazensko kljucanja se mogu odrediti primenom razlicitih empirijskih izraza. Jedan od najjednostavnijih, ali i najpoznatijih i najcesce koriscenihje i korelacija koju su razvili Shoш·ki i Judd [17]: 80 6. Prikaz i analiza doЬijenih rezu\tata s imulacije krize razmene toplote (6.8) Primenom ovog izraza za male vrednosti toplotnog fluksa od 1 О W/crn2 se doЬija da је h= 9486,833 W/m2K. Rohsenow [17] је razvio izraz za koeficijent prelaza toplote za bazensko kljucanje, gde Је: ih=...!L ЛУ'ь (6.9) gde је C.rJF0,015 [18] konstanta, koja zavisi od vrste zagrevanja i vrste zagrejacke povrsine i fluida, а n= l ,7 [18] је koeficijent koji se doЬija eksperimentalno. Ovaj izraz se koristi za kljueanje vode. Na ovaj nacin se doЬija h=6627,477 W/m2K. Sledeca korelacija za racunanje koeficijenta prelaza toplote koja је predloiena је Pioro- va korelacija [17]. 2 { } ~ - 3 т hl - q т Nu-CifK Pr , - -Csf 025 Pr k hљP~s [ o-g(p-pg )] . (6. 1 О) gde је C,r0,015 konstanta [18], koja zavisi od vrste zagrevanja i vrste zagrejacke povrsine (Ьakama ploca) i fluida (dvofazna mesavina voda-vodena para), а т је koeficijent koji se doЬija eksperimentalno [18]. Primenom ovog izraza se dobUa h= 3174,974 W/m2K. Veoma su eesto koriscene Kutateladze-ove korelacije [17]. ''Nova" Kutateladze-ova korelacija za kljueanje vode i pritisak od р= 1 bar: o-g м4 p-pg . =(:.)' 1 N -3 37 10-9 к-2м-4 d · N - hьl и - . х , g е Је и - , k DoЬija se h= 8213" 104 W/m2K. "Stara" Kutateladze-ova korelacija za kljucanje vode i pritisak od p=l bar: Nu = 0.44 Ко.1 Proзs ili h l = 0.44 1 х10-4 qp р Proзs ( ЈОЈ k gh1gpgµ р-Pg h= 5206,991 W/m2K. (6.11) (6.1 2) 81 6. Prikaz i analiza doЬijenih rezultata simulaci je krize razmene toplote Labuntsov korelacija [ 17] za bazensko kljueanje vode i pritisak od р= 1 bar: [ ( ) 0.67 ]( 2 ЈО.33 h = 0.075 1+10 Pg k q o.61 p-pg va (~ + 273.15) (6.1 3) h= 1712,954 W/m2K. Кruzilin korelacija [17] za dvofaznu me5avinu voda-vodena рага: (6. 14) h=8983,821 W/m2K. U taЬeli 6.9 је dat uporedni prikaz poredenja rezultata modela sa empirijskirn izrazima za razlicite vrednosti toplotnog fluksa. Rezultati modela su dati i za razliCite hrapavosti zagrejacke povгSine. Napominjem da u empirijskim izrazima ne postoje clanovi koji se odnose na hrapavost povr5ine, neki uzimaju u obzir karakteristike zagrejacke povrsine i dvofazne mesavine, а neke sarno vrednosti toplotnog fluksa, jednacina 6.8. 10000 9000 ....... ~ 8000 N Е ~ 7000 8 Mostinski 1 Starcweski ~ 6000 - о Foster-ZuЬer а. з О Shourki i Judd <О 5000 N 8 Rohsenow <О 11 а. 4000 O Pioro i: 111 :ff ..: 3000 8 Nova Kutateladze 8 О Starз Kutateladze ~ 2000 8 LaЬunstsov 1000 D Кruzilin о izraz Slika 6.25 Poredenje eksperimentalnih rezultata sa razlicitim empirijskim izrazima za odredivanje koeficijenta prelaza toplote za toplotni fluks 10 W/cm2 i n=l cm·2 82 h(WlmzК) Model n=l n=ll n=25 n=39 n=50 Mostinski 6. Prikaz i analiza dobijenih rezultata simulacije k:rize razmene toplote Tabela 6.9 Poredenje rnodela i empirijskih izraza za koeficijent prelaza toplote za kljucanje vode, zagrejacki zid od bakra i pri pritisku р= 1 bar IOW/cmz 40W/cm" 70W/cm2 85W/cm2 80 90 100 110 W/cm2 W/cm2 W/cm2 W/cm2 5524,9 12422,4 18918,9 22020,7 7407,4 12945,0 21621,0 23746,7 8620,7 13513,5 20000,0 24064,2 31746,0 27343,7 33742,3 67307,7 45045,0 i 5296,2 13976,7 20679,1 23689,5 22705,2 24656,6 26543,8 28375,1 Starczewski Foster-Zuber 8999,5 10333,8 10684,3 10793,4 12268,4 12757,8 13398,6 12430,5 Shourki i 9486,8 25035,9 37041,6 42434,0 40670,9 44166,З 47546,8 50827;2 Jttdd Rohseno\v 6627,5 16777,5 24409,8 27800,9 26694,3 28886;2 30999,0 33043,1 Pioro 3174,9 8000,47 11618,3 13223,8 23002,1 50269,7 233333 65406,5 "Nova" 8213,1 20676,6 30015,3 34158,6 32806,9 35484,0 38063,4 40557,8 K utateladze "Stara" 5206,9 13741,3 20330,8 23290,6 22322,9 24241 ,3 26096,8 27897,3 Kutateladze Labuntsov Кruzilin 1712,9 4336,4 6309,0 7185,5 6899,5 7466,О 8012,1 8540,4 8983,8 23709,4 35077,6 40184,1 38514,5 41824,5 45025,8 48132,2 6.2. Strujanje u vertikalnoj cevi kvadratnog poprecnog preseka U ovom delu је data numericka simulacija za slucaj vertikalnog stU]anJa u cevi kvadratnog poprecnog preseka. Zagrejacki zid је predstavJjen samo jednim slojem сеЩа. Ostali zidovi su adijabatski izolovani i za njih nije sprovedena simulacija, niti је odredivano temperatursko polje. Na ulazu u cev se zadaje strujanje preko profila brzine strujanja tecne faze. Pretpostavljeno је da se pre pocetka procesa zagrevanja и cevi kvadratnog poprecnog preseka nalazi 50% tecne faze i 50% рате faze. Zagrejacki zid је podeljen na cetiri zone, svaka zonaje povгSine 2 crn2, sto odgovara gustini nukleacije od п=О,5 crn-2• Moze se primetiti talasasti profil koji razdvaja tecnu i parnu fazu. Javlja se kao posJedica intenzivnog procesa isparavanja, delovanja sile zemJjine gravitacije i kao posledica uvodenja tecne faze na dnu isparivacke cevi. Temperatursko polje u zagrejackom zidu је slicno kao i za slueaj bazenskog kljucanj a, javljaju se ternperaturski skokovi na rnestirna gde se generisu mehurovi i gde zapreminski udeo рате faze prelazi 0,9. Ternperaturski skok је i u ovorn slucaju oko 50 К, slika 6.28. Ternperatura dvofazne mesavine ostaje nepromenjena i pribliZno jednaka temperaturi saturacije. 83 130 W/cm2 43918,9 31895,1 13805,6 57132,3 36956,5 176935 45330,8 31357,9 9551,9 54103,0 6. Prikaz i analiza dobijenih rezultata simulacije krize razmene toplote 0.08 0.06 s; (/) :;· ~0.04 ~ 0.02 0.04 ~ 0.02 ftl~ ·~(;:' q х ~у FI 0.776848 0.582636 0.388424 о 194212 Slika 6.26 Zapreminski udeo рате faze pri vertikalnom stтujanju na gore za toplotni tluks od 1 ОО W/cm2, gustinu nukleacij e n=0,5cm·2, и trenutku t=O,Ol s 0.07 0.06 FI 0.815241 0.760891 0.05 0.706542 е 0.652192 0.597643 ~.04 0.543494 ·с: 0.469144 СЂ 0.434795 о.оз 0.360446 0.326096 0.271747 0.217397 0.02 0.163048 о. 106699 0.0543494 0.01 0.02 0.04 0.06 0.08 Visina (m) Slika 6.27 Zapreminski udeo parne faze pri vertikalnom strujanju na gore za toplotni fluks od 100 W/cm2, gustinu nukleacUe n=0,5cm·2, u trenutku t=O,Ols, 2D prikaz 84 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote -3 ......... 001 Dv.1... . 0.02 V/f). о.оз q (!ђ) 0.04 0.03 0.02 1«') 0.01 • э. ~,. S'f..'f' Нi====*===~ о.ов 0.02 о Sirina (m) х ~z у Т1 417.91 1 406.68 395.45 384.219 Т1 417.911 406.68 395.45 384.219 Slika 6.28 Temperaturska raspodela unutar zagrejackog zida pri vertikalnom strujanju па gore za toplotni fluks od 1 ОО W/cm2, gustinu nukleacije n=0,5cm·2, и trenutku t=O,O 1 s Zatim је sprovedeno strujanje u cevi kvadartnog poprecnog preseka, gde se jedna strana cevi zagreva, а zid koji se zagreva је predstavljen pomocu sloja od 10 celija. Rezultati simulacije su dati па slikama 6.29 i 6.30. 85 6. Prikaz i anaJiza dobijenih rezultata simulacije krize razmene toplote < ~: ::::1 0.08 0.06 Q) 0.04 3 - 0.02 о о.оз 0.02 0.01 DuЫna(m) 0.02 0.04~ .§ х ~у FI • 0.78972 0.59229 0.39486 0.19743 О.О~~ ·~ о С<) 0.04 Sirina (m) FI 0.925312 0.863625 0.801937 0.74025 0.678562 0.616875 0.555187 0.4935 0.431812 0.370125 0.308437 0.24675 0.185062 0.123375 0.0616875 0.06 0.08 Slika 6.29 Zapreminski udeo tecne faze za toplotni fluks od 1 ОО W/cm2, gustinu nukleacije n=0,5cm·2, и trenutku t=0,2 s PretpostavJjeni su isti uslovi kao i и prethodnom slueaju. Rezultati simulacije su isti kao i и slucaju vertikalnog strujanja па gore kada se zagrejacki zid predstavlja jednim slojem celija. Jedina је razlika и temperaturskom polju и zidu. Temperaturski skok је ist i kao i и prethodnim slucajevima, oko 40 К, samo se primecuje preciznije temperatursko polje unutar zagrejackog zida. Temperatura dvofazne mesavine је nepromenjeпa i priblifuo jednaka temperaturi saturacije. Na samom ulazu и kanal se ne javlja povecanje temperature, jer па ulazu postoji dotok tecne faze, Cija је temperatura jednaka temperaturi saturacije i па taj пacin se vrsi hladenje zagrejackog zida па ulazu. 86 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote < ёi)" ~- 0.04 i.&- ---'-------" 1 3 0.02 о Sirina (m) х z_d Т1 417.911 406.68 395.449 384.218 Slika 6.30 Temperatursko polje za toplotni fluks od 100 W/cm2, gustinu nukJeacije п=О,5 cm·2, u trenutku t=0,2 s 6.3. Strujanje u cilindricnoj vertikalnoj cevi Za ovaj slucaj simulacije Ьilo је neophodno uvesti polarno-cilindricni koordinatni sistem. Matematicki model i algoritam za resavanje jednacina odrZa.nja u dati и poglavlju 4 i poglavlju 5. Simulacija је radena za cevni isecak, pri cemu se ugao menja od 0° do 45°. Pretpostavljeno је da se zid zagreva unifomшo i zid је podeljen na 4х4 zone od kojih se u svakoj zoni generise jedan mehur. Povrsina zone је 3,5325 cm2, sto odgovara gustini nukleacije 0,28 cm·2. PoloZa.j mehura unutar zone se odreduje stohasticki, а broj mehurova zavisi od karakteristika materijala i od hrapavosti zagrejacke povrsine. U ovom slucaju је jedan mehur u jednoj zoni na zagrejackom zidu. z ~ 0.11 01 х 0.1 009 F1 0.06 0.560212 о S2285S 0485517 ЕОО7 044817 0•10822 ';' 0 .06 0.373475 FI с 0.336127 i '"'"' ~ 0.05 0.29878 0.261•32 0.597906 0.04 0.224065 0.398604 ~ 0.186737 L- 0.14939 0.199302 0.03 0.1 12042 ~ О.07<4696 0.02 00373475 0.01 оо 0.05 0.1 Sirina(m) Slika 6.31 Zapreminski udeo ~arne faze za strujanje u cilindricnoj vertikalnoj cevi za toplotni fluks od 100 W/cm , gustinu nukleacije п=О,28 cm·2, и trenutku t=0,01s 87 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote Prvo је uradena simulacija samo za slucaj dvofazne mesavine, slika 6.31 , а и nastavku је и model uvrsten i zagrejacki zid i doЬijeni su rezultati simulacije i u tom slueaju. Na slikama 6.32-6.35 su dati rezultati simulacije za slucaj strujanja dvofazne mesavine u cilindricnoj cevi kada је zagrejacki zid kruZпi luk modeliran jednim slojem 6elija. Primecuje se pove6anje temperature dvofazne mesavine tokom vremena. Zapreminski udeo рате faze se ne menja znacajno, sto је posledica stalnog ubacivanja tecne faze па ulazu u isparivacku cev. 0. 1 ~ 0.075 5· Q) ~ 0 05 0.11 0.1 0 .09 0.08 I O.D7 (\! 0.06 с ·с;; 5 0.05 0.04 о.оз 0.02 0.01 0.05 А 0.612441 0.5716 11 0.530782 0.489953 0.449123 0.408294 0.367464 0.326635 0.285806 0.244976 0.204147 0.163318 0.122488 0.0816588 0.0408294 Precnik (m) z ~ FI 0.799776 0.599832 0.399888 0.199944 0 .1 х Slika 6.32 Zapreminski udeo pame faze za strujanje u cilindricnoj vertikalnoj cevi za toplotni fluks od 100 W/cm2, gustinu nukleacije п=О,28 cm-2, и trenutku t=O,Ols Moze se zakljuciti da se model koji је predlozen и ovom radu moze koristiti i za simulaciju strujanja dvofazne mesavine u cilindricnoj vertikalnoj cevi и koju se sa donje strane ubacuje dvofazna mesavina uz male modifikacije modela, poglavlja 4 i 5. 88 6. Prikaz i analiza dobiienih rezultata s imulacije krize razmene toplote z 0.1 k_ х < 11?: 0.075 - Т1 ::i !!) i 373 ......... 3 0.05 372.9 ......... 372.8 372.7 0.025 Slika 6.33 Raspodela temperature za strujanje u cilindricn~~ vertikalnoj cevi za toplotni fluks od 100 W/cm2, gustinu nukleacije п=О,28 cm-, и trenutku t=O,Ol s z k_ 0.11 х 01 F1 0.58556'2 009 0.5'652S 0.507487 0.08 0-S 0.429412 ~0.о7 03903'/S s. Q.351337 FI • 0.06 0.3123 i '"~ ~ 0.05 0.273282 0.234225 0.59922 ...._ 0.195187 0.39948 r- 0.151515 0.04 1 о 117112 0.19974 0.078075 о.оз 0.0390375 0.02 0.01 0.1 $Ма(m) Slika 6.34 Zapreminski udeo рате faze za strujanje u cilindricnoj vertikalnoj cevi za toplotni fluks od 100 W/cm2, gustinu nukleacije п=О,28 cm-2, и trenutku t=O, l s 89 6. Prikaz i analiza doЬijenih rezultata simulacije krize razmene toplote < !!!: :::i Q) 0.1 30.05 ......... о О. Precnik (m) z к_ Т1 ii 450 424.333 398.667 373 х Slika 6.35 Raspodela temperature za strujanje и cilindricnoj vertikalnoj cevi za toplotni fluks od 100 W/cm2, gustinu nukleacije п=О,28 cm·2, и trenutku t=O,l s 90 7. Zakljucak 7. Zakljucak U okviru doktorske disertacije је modelskim pristupom i numerickim simulacijama istraZivan mehanizam pojave krize kljucanja. Razvijen је trodimenzionalni model dvofaznog strujanja tecne i parne faze sa granicnim uslovima na zagrejackoj povrsini kojima se odreduju uslovi generacije parnih mehurova koji dovode do krize kljucanja. Dvofazna mesavina је predstavljena modelom dva tluida, koji se sastoji od jednacina odrZa.nja mase, kolicine kretanja i energije za svaku od faza. Medufazno dejstvo је uzeto и obzir primenom odgovarajucih konstitutivnih korelacija. Dinamika generacije mehura na zagrejackom zidu је modelirana pomocu gustine ceлtara nukleacije i vremena rasta mehura na zidu. Zagrejacka povrsina је podeljena na zone, ciji је broj ро jedinici povrsine jednak gustini centara nukleacije, а lokacija centra nukleacije u zoni је odredena funkcijom slucajnih brojeva. Matematicki model је numericki resen metodom konacnih zapremina i koriscenjem SIМPLE (Semi-Implicit Method for Pressure-Linked Equations) algoritrna koji је modifikovan za potrebe spregnutog resavanja Ьilansnih jednacina dvofaznog strujanja. Sprovedena је trodimenzionalna numericka simulacija bazenskog kljucanja za zasicenu tecnost pri atmosferskim uslovima sa ciljem predvidanja toplotnog fluksa. Rezultati pokazuju vecu okvasenu povrsinu na zagrejackom zidu pri manjim vrednostima toplotnog fluksa i povrsinsko zasusenje za kriticne vrednosti toplotnog fluksa. Pokazano је da smanjenje gustine centara nukleacije vodi ka smanjenju vrednosti kriticnog toplotnog fluksa. DoЬijeni rezultati kriticnog toplotnog fluksa se dobro sla.Zu sa raspolozivim eksperimentalnim rezultatima. Ovaj pristup је originalan kako zbog primene dvofluidnog dvofaznog modela za predvidanje krize kljueanja, tako i zbog definisanja novih granicnih uslova na zagrejackoj povrsini. Pokazano је da hrapavost povrsine znacajano utice na vrednosti kriticnog toplotnog fluksa. DoЬijeni numericki rezultati pokazuju da se sa povecanjem hrapavosti zagrevane povrsine povecava kritiCni toplotni fluks, jer se sa povecanjem hrapavosti povecava gustina klijalista mehurova. Parametarski је analiziran uticaj hrapavosti zida zagrejacke cevi, ugla kva5enja tecnosti na zagrejackoj povrsini i vrednosti toplotnog fluksa na pojavu krize kljucanja. Na osnovu ove analize је moguce odredivanje kriticnih vrednosti toplotnog fluksa za razlicite uslove kljucanja, kao i predvidanje pojave mehanizma pregrevanja isparivackih cevi. Na osnovu predlozenog modela i izvrsene numericke simuJacije Ьilo је moguce izracunati koeficijent prelaza toplote sa zagrejackog zida na dvofaznu mesavinu. Rezultati ovog proracuna su uporedeni sa empirijskim izrazima iz literature koji su korisceni za izracunavanje koeficijenta prelaza toplote za razlicite fluide i razliCite materijale i karakteristike zagrejacke povrsine. Rezultati numerickih simulacija su uporedeni sa raspolozivim eksperimentalnim rezultatima. Dati su dijagrami na kojima su uporedeni rezultati doЬijeni primenom modela i eksperimenta na kojima је predstavljena zavisnost broja klijaliSta mehurova (zavisi od hrapavosti povcline) i vrednosti kriticnog toplotnog fluksa, kriva kljucanja (zavisnost toplotnog fluksa i razlika temperature zida i temperature saturacije) i zavisnost koeficijenta prelaza toplote i razlike temperature zida i temperature saturacije. Eksperimenti su radeni pod razlicitim uslovima za razlicite hrapavosti zagrejackih povrsiлa, ali se moze zakljuCiti da su rezuJtati poredenja zadovoljavajuci. 91 7. Zakljucak Na osnovu rezultata ovog istra.Zivanja moze se zakljuciti da se ovde predstavljeni numericki postupak moze koristiti za predvidanje procesa koji se odvijaju u isparivackim cevima za razlicite parametre procesa. Takode se moze upotrebiti za doЬijanje optimalnih vrednosti geometrijskih karakteristika isparivackih cevi, razlicite materijale i kvalitete obrade zagrejacke povrsine od kojih se izraduju cevi, kao i strujno-termicke procese koji su promenljivi za razlicite uslove procesa kljucanja. Ispitivanje uslova pri kojima dolazi do pregrevanja isparivacke cevi, do njenog oStecenja, topljenja ili pucanja је od izuzetnog znacaja za uslove sigurnosti i pouzdanosti svakog generatora pare u procesu proizvodnje. Sa stanovista inzenjerske prakse razvijeni numericld postupak omogucava pouzdano predvidanje kriticnog toplotnog f1uksa, sto је veoma znacajno za analize sigumosti razlicitih vrsta generatora pare, kao i za razvoj novih generatora pare u cilju postizanja vece vrednosti kriticnog toplotnog fluksa, odnosno vece termicke pouzdanosti. Dalja istrazivanja koriscenjem numerickog postupka razvijenog и ovoj disertaciji mogu Ьiti usmerena ka predvidanju uslova nastanka krize kljucanja i vrednosti kriticnog toplotnog fluksa u strujnim kanalima razlicitih oblika i dimenzija poprecnog preseka i pri strujanjirna razlicitih fluida, odnosno za uslove koji su od znacaja za pouzdan rad i termicku sigumost generatora pare. 92 Literatura Literatura 1. S. Nukiyama, Maximurn and Мinimum Values of Heat Transmitted from Metal to Boiling Water Under Atmospheric Pressure, Journal о/ the Japanese Society о/ Mechanical Engineers, Vol. 37, (1934), рр. 367-374, (in Japanese), lnternational Journal o/HeatandMass Transjer, Vol. 9, (1966), рр. 1419-1433. (in English) 2. К. Ма, С. Рал, Тhе effect of heated wall thickness and materials on nucleate boiling at high heat tlux, lnternational Communications in Heat and Mass Transfer, Vol. 26, Issue 8, (1999), рр. 1103-1114 3. G. Р. Celata, М. Cumo, У. Katto, А. Mariani, Prediction of the critical heat flux in water subcooled flow boiling using а new mechanistic approach , lntemational Joumal ој Heat and Mass Transjer, Vol. 42, Issue 8, (1999), рр. 1457-1466 4. G. Р. Celata, К. Mishima, G. Zummo, Critical heat flux prediction for saturated tlow boiling of water in vertical tubes, lnternational Joumal ој Heat and Mass Transfer, Vol. 44, (2001), рр. 4323-433 1 5. А. Inoue, S. Lee, Influence of two-phase flow characteristics on critical heat flux in low pressure, Experimental Тhermal and Fluid Science, Vol. 19, (1999), рр. 172-181 6. Ј . С. Sturgis, 1. Mudawar, Critical heat tlux in а long, rectangular channel subjected to one-side heating-П. Analisis of critical heat flux data, lnternational Journal ој HeatandMass Transjer, Vol. 42, (1999), рр. 1849-186 7. 1. L. Pioro, W. Rohsenow, S. S. Doerffer, Nucleate pool-boiJing heat transfer. 1: review of pararnetric effects of boiling surface, Intemational Journal о/ Heat and Mass Transfer, Vol. 47, (2004), рр. 5033-5044 8. I. L. Pioro, W. Rohsenow, S. S. Doerffer, Nucleate pool-boiling heat transfer. П: assessment of prediction methods, lntemational Joumal of Heat and Mass Transfer, Vol. 47, (2004), рр. 5045-5057 9. I. L. Pioro, Experimental evaluation of constants for the Rohsenow pool boiling correlation, Jntemational Joumal ој Heat and Mass Transjer, Vol. 42, (1999), рр. 2003-20013 10. У. Не, М. Shoji, S. Maruyama, Numerica1 study of high heat flux pool boiling heat transfer, lnternational Joumal ој Heat and Mass Transfer, Vol. 4, (2001), рр. 2357- 2373 11. М. Shoji, Studies ofboiling chaos: а review, lnternational Journal ој Heat and Mass Transjer, Vol. 47, (2004), рр. 1105-1128 12. Т. G. Theofanous, Ј. Р. Ти, А. Т. Dinh, Т. N. Dinh, Тhе boiling crisis phenomenon Part I : nucleation and nucleate boiling transfer, Experimental Тhermal and Fluid Science, Vol. 26, (2002), рр. 775-792 13. Т. G. Theofanous, Т. N. Dinh, Ј. Р. Tu, А. T.Dinh, The boiling crisis phenomenon Part ll : dryout dynamics and burnout, Experimental Тhermal and Fluid Science, Vol. 26, (2002), рр. 793-810 14. I. 1. Gogonin, 1.А. Kutateladze, Critical heat flux as а function of heater size for а liquid boiling in а large enclosure, Ј Eng. Phys. , Vol. 33 (5), (1977), рр. 1286-1289 15. 1. Vaquila, L. I. Vergara, М. С. G. Passeggi, R. А Vidal, Ј. Ferron, Chernical reactions at swfaces: titanium oxidation, Surf Coat. Technol, Vol. 122, (1), (1999), рр. 67-71 16. М. Kureta, Т. Нibik:i, К. Мishima, Н. Akimoto, Study on point of net vapor generation Ьу neutroo radiography in subcooled boiling flow along oarrow rectangular channels with short heated length, Intemational Journal ој Heat and Mass Transfer, Vol. 46, (2003), рр. 1171-1181 93 Literatura 17. М. Kureta, Н. Akimoto, Critical heat flux correlation for subcooled boiling flow in nazтow channels, International Joumal ofHeat and Mass Transfer , Vol. 45, (2002), рр.4107-4115 18. Z. Stosic, V. Stevanovic, Тhree--dimensional numerical simulation of burnout on horizontal surface in pool boiling, Proceedings of FEDSM'OЗ, 4 th АSМЕ!ЈSМЕ Joint Fluids Engineering Conftrence (Advances in CFD), Honolulu, Hawaii, USA, July 6- 11, 2003 19. Т. Hibiki, M. lshii, Active nucleation site density in boiling systems, International Journal of Heat and Mass Transfer Vol. 46, (2003), рр. 2587- 2601 20. L.A. Payan-Rodriguez, А. Gallegos-Muiioz , G.L. Porras-Loaiza , М. Picon- Nuiiez, Critical heat flux prediction for water boiling in vertical tubes of а steam generator, International Journal of Тhermal Sciences , Vol. 44, (2005), рр. 179- 188 21. J.Мitrovic, How to create ап efficient surface for nucleate boiling, Intemational Journal ofТhermal Sciences, Vol. 45, (2006), рр. 1-15 22. B.R. Vijayarangan, S. Jayanti, A.R. Balakrishnan, Studies on critical heat flux in flow boiling at near critical pressures, International Journal of fleat and Mass Transfer, Vol. 49, Issues 1-2, (2006), рр. 259-268 23. Ј. Le Corre, S. Уао, С. Н. Amon, А mechanistic model of critical heat flux under subcooled flow boiling conditions for application to опе- and three-dimensional computer code, Nuclear Engineering and Design, 2009. 24. N. Ivanov Kolev, Herzogenaurach, Germany, Letter to the Editor, Nuclear Engineering and Design, Vol. 239, (2009), рр. 187-192 25. N. Afgan, Prenos toplote kljucanjem i toplotni jluks pregrevanja smesa etilalkohol-benzol,_ Doktorska disertacija, Institut za nuklearne nauke "Boris Кidric"~ decembar 1965, Univerzitet u Beogradu, Beograd 26. S. Lazic, Pregrevanje tecnih kapi, Magistarski rad, Institut za nuklearne nauke "Boris Кidric", maj 1966, Univerzitet и Beogradu, Beograd 27. D. Spasojevic, Minimalno toplotno opterecenje isparivackog kanala па granici stabllnosti parametara, Magistarski rad, Masinski fakultet, Univerzitet и Beogradu, avgust 1970, Beograd 28. V. Valent, Dinamika rasta mehura i analiza prenosa toplote i mase p1·i kljucanju smesa etilalkohol-voda, Magistarski rad, Tehnolosko-metaluski fakultet, septembar 1971, Univerzitet и Beogradu, Beograd 29. М. Stefanovic, Analizajluktuacije temperature и dvofaznom toku pri pothlatlenom kljucanju, Doktorska disertacija, MзSinski fakultet, Unjverzitet u Beogradu, 1975, Beograd 30. V. Jovic, Prilog proucavanju nestacionarnog strujanja dvofazne smese voda- vazduh и paralelnim kanalima, Dok:torska disertacija, MзSinski fakultet, Univerzitet и Beogradu, 1992, Beograd 31. Р.В. Whalley, Two-Phase Flow and Heat Transfer, Oxford University Press, 1996. 32. М. Pezo, V. D. Stevanovic, Z. Stevanovic, А two-dimensional model of the kettle reboiler sbell side theпnal-hydraulics, Intemational Journal of Heat and Mass Transfer , 49, (2006), рр.1214-1224 33. S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, 1980. 34. V. D. Stevanovic, Thermal-hydraulics of steam generators - mode/ling and numerical simulation, University of Belgrade, Faculty ofMechanical Engineering, Belgrade, 2006 94 Literatura 35. М. Pezo, V. Stevanovic, Mathematical Modelling and Numerical Simulation of Burnout in Pool Boiling, Dubrovnik 4th Conference оп SustainaЫe Development ој Energy, Water and Environment Systems, Dubrovnik, Croatia, June 2007 36. М. Pezo, V. Stevanovic, Numerical prediction of critical heat flux in pool boiling with the two-fluid model, International Journal ој Heat and Mass Transfer, Vol. 54, (2011), рр. 3296-3303 37. W. Н. McAdams, Heat Transmission, 2d ed., McGraw Company, Inc., New York, 1942 95 Biografija: lme i prezime: Daturn rodenja: Mesto rodenja: Porodicno stanje: Skolovanje: 1980.-1988. 1988.-1992. 1992.-1997. 5.12.1997. 1999.-2004. 8.6.2004. Кretanje u poslu: 1998.-1998. 1998.-1999. 1999.-2004. 2004.-2011 . Мilada Pezo 24.10.1973. Sarajevo, Jugoslavija jedno dete Osnovna skola "UZicka RepuЫika" u Beogradu IX Beogradska gimnazija "Мihailo Petrovic-Alas" u Beogradu Stuctije na Masinskom fakultetu u Beogradu, odsek Automatsko upravljanje Odbranjen diplomski na Masinskom fakultetu u Beogтadu PosledipJomske studije na Masinskom fakultetu u Beogradu, odsek Terrnoenergetika Odbranjen magistarski rad na Masinskom fakultetu u Beogradu sa temom "Numericka simulacija i analiza dvodimenzionog dvofaznog strujanja и horizontalnom isparivacu" Fabrika za proizvodnju metalnih pтoizvoda Intersilver, Beograd Stipendista Мinistarstva za nauku i tehnoloski razvoj u Laboratoriji za termotehniku i energetiku lnstituta za nuklearne naukeVINCA Istra.Zivac pripravnik u Laboratoriji za termotehniku i energetiku Instituta za nuklearne nauke VINCA Istra.Zivac saradnik u Laboratoriji za termotehniku i energetiku Instituta za nuklearne nauke VINCA ~3)". \ . ~ • 96 lnl,m.lr:tonal j()llJ'Tl.11 of Hear .md Mass Tr.uisfer S4 {2011) 3296-3303 Contents lists available at ScienceDirect International journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ ijhmt Numerical prediction of critical heat flux in pool boiling with the two-fluid model Mjlada Pezoa. Vladimir Stevanovicb.• •Laboratory ofThennal Engineering and Enetgy. "VINCA'" lnstirure of Nudeor Sdences. P.O. Box 522. I 1001 Belgrade, Serbia •Faculty of Mtchanlcal Erlgineerlng, University of Belgrade, KroUice Marije 16, 11120 Belgrade 35, Serbia ART ICLE I N FO ABSTRACT Artlde hisrory: Received 14 March 2011 Accepted 23 March 2011 Available online 23 April 2011 Keywords: Pool boiljng Critical heal Oux Modeling Three-dimensional numerical simulations of the atmospheric saturated pool boiling are performed with the aim of predicting the critical heat flux. The two-phase mixture in pool boiling is described with the transient two-fluid model. The transient heat conduction in the horizontal heated wall is also solved. Dynamics of vapor generation on the heated wall is modeled through the density or nucleation sites and the bubble residence time on the wall. The hearer's surface is divided into zones, which number per unit area equals the density of nucleation sites, while the location of nucleation sire within each zone is determined by a random function. The results show a replenishment of the heater's surface with water and surface wetting for lower heat fluxes, while heater's surface dry-out is predicted at critical heat flux values. Also. it is shown that the decrease or nucleation site density leads co the reduction or critical hear flux values. Obtained results of critical heat flux are in good agreement with available measured data. The presented approach is original regarding both the application of the two-Ou id two-phase model for the prediction or boiling crisis in pool boiling and the defined boundary conditions at the heated wall surface. 1. Introduction The critical heat flux (CHF) occurs in boiling when generated va- por covers a cerrajn area of the heated surface. which leads to the abrupt rise of the heated sutface temperature. This phenomenon is called the boiling crisis since the vapor layer prevents the contact of liquid with the heated surface and the formation of bubbles. Also. the temperature rise can lead to a heated wal l destruction. which is characterized as the burnout. The phenomenon is also known under the names departure from nucleate boiling (DNB) and dry-out The boiling crisis occurs both in a pool boiling and a heated channel with a coolant flow and convective boiling. This pa- per is related to the CHF in pool boiling. i.e. to the boiling of a stag- nant buJk of liquid that is in contact with a heated surface at a temperature higher than the saturation temperature of the liquid. Prediction of the CHF is very important for safery and efficiency of thermal equipment. such as steam generators or heat exchangers due to a potential damage of a heated sutface or a severe deterio- ration of the heat transfer at heat exchanger's surfaces. Due to its complexity and the importance for thermal safety, the boiling crisis has been attracting researchers for a long time. The boiling crisis was experimentally obseived by Nukiyama [1 ). A classical model of boiling crisis in pool boiling is based on the works of Kutateladze (2,3) and Zuber [4 [. The model assumes that • CorttSpondlngauthor. Tel.: +381 113370561 : fax: +38111 3370 364. E·rnai/ addresses: m.itadaOvinca.rs (M. Pezo). vstevmo111c-.nas.bg.ac.rs (V. St:evanovic). 0017-9310/S - see front m..u:rer o 2011 Elsevier Ltd. All rlghrs reserved. dol: l 0. 1016/j.ljht'~rmaSStTansfer .2011.03.057 © 2011 Elsevier Ltd. AU rights reseived. the CHF is determined by the Taylor instability of the vapor layer fonned above the heared surface and the Kelvin- Helmholtz insta- bility of the vapor stems through which the generated vapor es- capes from the heated sutface into the bulk of two-phase mixture [SI. Based on this hydrodynamic model the correlation for the critical hear flux was developed: (1) where a value of the constant C was determined as 0.131 by Kuta- teladze (21 for che boiling on horizontal cylinders. and as rt/ 24 -"I 0.1308 by Zuber [4 [ and 0.149 by Lienhard 151 forthe flat plate geometry of che heated surface. Theofanous and co-workers (6,7J revised chls model. Their experimental investigation showed that the CHf value is influenced not only by the hydrodynamic condi- tions in the vicinity of the heated surface. but also by the heater's surface micro-conditions, which they characterized as fresh. med- ium aged and heavily aged surfaces. For each age type of heater's surface they determined the nucleation site density prior to CHF. Depending on the heater's sutface micro conditions. measured CHF values were in the range of±40% of the values calculated by Eq. (1 ). The other conclusion in (6,7) is that there is no structured pattern of the boiling two-phase mixture at the occurrence of boil- ing crisis, hence, there is no validation for the existence of the vapor seems above the vapor blanket at a surface exposed co high heat flux. Bang et al. (SI also reported visual observations of the rwo- phase mixture pacrem at the instance of the boiling crisis. Again, M. Pezo, V. Srevanovic/ Jnremational journal of Hear and Mass Transfer 54 (2011 ) 3296-3303 3297 Nomenclature a thermal diffusivity (m2/s) b width of the nucleation zone (m) Cp specific heat Ufkg K) Co interfacial drag coefficient D diameter (m) Db bubble departure diameter (m) F force per unit volume (N/m3) g gravitational acceleration (m/s2) h enthalpy Ufkg) h12 latent heat of evaporation U/kg) k thermal conductivity (W/mK) Le water capillary length (m) n density of nucleation sites ( m- 2 ) p pressure (Pa) q heat flux (W/m2 ) qh volumetric heat rate (W/m3 ) qb volumetric heat source for bubble generation on the heater's surface (W/m3) t time (s) T temperature (K) u velocity (m/s) x coordinate ( m) these results did not show the existence of the vapor stems and the Kelvin-Helmholtz instability as the reason for the CHF. Instead, fast camera photographs showed that CHF occurs when a liquid film on a heated surface is evaporated beneath a vapor blanket. where the vapor blanket is formed by the coalescence of bubbles in intensive boiling. The existence of vapor stems was also not observed in the investigation presented by Chung and No [9]. Jnstead, the existence of a local nucleate boiling under the large vapor film and at the edge of it was observed. Consequently, in [9] the CHF is determined by the boiling activity and dry spots originating from nucleate. boiling phenomena. A number of models have been developed for the CH.F prediction. and a few are mentioned here. Dhir and Liaw (10] pre- dicted CH.F under the assumption of the existence of the liquid lay- ers surrounding stationary vapor stems, Kolev [ 11) a.nd Zhao et al. [12] presented models of nucleate boiling and CHF based on the individual bubble behavior, while He et al. [13) developed a model based on a two-phase pattern formed by the liquid layer on the heated surface and the vapor stems. while CHF occurs due to the depletion of the liquid layer. Hence, presented models of CHF are based on assumed local patterns of two-phase mixture on the heated surface. But all these assumptions are subject to discussion since the.re is no firrn agreement that structured pattern exists at a heated surface at the instance of CHF [6,7). The approach for the boiling crisis prediction was developed in (14) without the pre- assumptions about the two-phase mixture structure on the heater's surface. The two-fluid model of two-phase mixture in pool boiling was applied, but without the calculation of the heated wall tran- sient conduction and temperature increase at the CHF occurrence. Instead, the CHF was indicated indirectly by the steam void close to 1 in the two-phase mixture adjacent to the heater's surface. for a time period that equals or exceeds experimentally observed per- iod between the incipient of temperature rise and the instant of burnout. In the present work the boiling crisis in pool boiling is investi- gated by the application of the two- fluid model for the prediction of water-steam two-phase mixture conditions on the heated hor- izontal wall. Hence, no assumptions are introduced about the two- phase mixture pattern on the heater's surface in pool boiling and under CHF conditions. Instead, the generated steam upward flow from the heated wall surface and the water circulation in the pool Greek ex void fraction r phase transition rate (kg/m3s) 8 wetting contact angle (degree) p density (kg/m3) (J surface tension (N/m) ' phase change relaxation time (s) Indexes c condensation e evaporation k phase index (k = 1, 2) p particle sat saturation w wall parameter 1 water 2 ·steam 21 interfacial saturated liquid II saturated vapor are predicted by solving the mass, momentum and energy balance equations for each (vapor and liquid) phase. The interface transfer processes between the liquid and vapor phase, i.e. the interface friction and the evaporation rate are calculated by appropriate clo- sure laws. The two-fluid model is coupled with the model of tran- sient three-dimensional heat conduction in the heated wall. The boundary conditions between rhe heated wall and the boiling two-phase mixture above it are modeled with the bubble nucle- ation site density, the bubble residence time on the heated wall and with a certain level of randomness to the location of bubble. nucleation sites. Numerical simulations have been performed with the developed model for high heat flux conditions, which lead to the boiling crisis. The boiling crisis and corresponding CHF are de- tected by a rapid increase of the calculated wall surface tempera- ture. This condition occurs when the heater's surface is covered with a vapor blanket that hinders the replenishment of the surface by the pool water. The calculated CHF values are compared with available experimental data and good agreement is obtained for a range of nucleation site density values. The applied numerical and modeling method has shown robustness by providing stable calculations for wide ranges of applied modeling parameters of pool boiling, such as the density of nucleation sites and the bubble residence time. 2. Problem statement Pool boiling is simulated in a square vessel, initially filled with saturated stagnant water up to 0.02 m, Fig. 1. The vessel is open to the atmosphere. The water level corresponds to the collapsed two- phase mixture level of the physical experiments presented in (6,7). The bottom copper wall of 2 mm thickness is heated by a uniform volumetric heating source. At the beginning of the simulation, the temperature of the wall is equal to the saturated temperature of water. At the initial state, the. coUapse.d and swell levels coincide, while later on, during the vapor generation, the swell level dynamically moves, and its position is predicted with the two-fluid model pre- sented in Section 4. It was found that the inclusion of the swell le- vel, as the upper boundary for the liquid flow domain, is necessary 3298 M. Pezo, V. Stevanavic / International joumal of Heat and Mass Transfer 54 (2011) 3296-3303 generation ¥aU 0.25a 0.5a 0.15a a Wid!h(m) 2 3 4 5 6 7 8 Zones 9 10 II 12 l3 14 15 16 fig. 1. Numerical grid used in the simulation of pool boiling. in order to reliably predict the water recirculation in the pool boil- ing. The volume above the swell level, up to the vessel exit is filled with vapor. There is no water feeding during the boiling, which leads to the constant water depletion due to the evaporation. For the evaporation of the total water mass within the vessel under very high heat flux of 1500 kW/m2• a period of approximately 27 s is needed. Since numerically predicted pool boiling steady- state conditions or the CHF are reached only a few seconds after the initial state. the depletion of water mass has no influence on here presented computational results. Dynamics of pool boiling is modeled according to the following algorithm. The surface of the bottom wall is divided into 16 equal zones as shown in Fig. l (4 >< 4 square zones are used). The width of zone is varied according to the density of nucleation sites. Each zone represents one nucleation site, and each zone comprises ten or more control volumes. The grid refinement tests have shown that with ten or more control volumes within a nucleation zone there is no change in the calculated critical heat flux or two-phase flow structure in the pool boiling. A location of bubble generation is randomly chosen among the control volumes comprised within the zone. where only one bubble is generated within each zone. The density of nucleation sites and corresponding width of the zone for random bubble generation are input parameters for here presented simulations. The relation between the density of nucle- ation sites n and the zone width b can be derived from the simple geometric condition that one square meter is covered with n nucle- ation sites, i.e. (2) The period of vapor generation at a randomly chosen location (nucleation site) is equal to the bubble residence time on the hea- ter's surface. It is also assumed that the heat transfer from the hea- ter's surface towards the fluid takes place at the nucleation sites and it is totally consumed for the vapor generation. The heat transfer from the heater's surface towards the liquid or vapor phase in con- trol volumes without vapor generation is neglected, i.e. the sensible heat transferred to the liquid or vapor phase is neglected. The vessel volume and the heater's wall are discretized with the three-dimensional mesh in Cartesian co-ordinates. Fig. 1. Dimen- sions of the control volume are chosen in order to catch the dynamics of the vapor generation at the heater's surface (mode.led by the nucleation zones and the random choice of the nucleation site within the zone). As presented, the micro conditions at the heater's surface, which determine the vapor generation, are pre- scribed with the width of the zone, the randomly chosen nucle- ation site within the zone and the bubble residence time on the heater's surface. The velocity boundary conditions at the pool lateral walls are prescribed as the so-called numeric adiabatic conditions for the velocity components (i.e. there is no change of the velocity compo- nents in the direction perpendicular to the wall), while the velocity is zero on tb.e heater's surface. Also, the adiabatic conditions are prescribed for the enthalpy and temperature on the boundaries of the calculation domain. 3. Mode.ling of vapor generation dynamics The density of nucleation sites has significant influence on the boiling dynamics and conditions under which the boiling crisis is reached. The nucleation site density is determined with the hea- ter's surface roughness. the wetting contact angle. thermo-physical characteristics of the boiling fluid and the wall heat flux. All these boiling conditions can vary the density of nucleation sites in a wide range. Some empirical correlations for the prediction of nucleation sites density are addressed below and assessed with the experi- mental results from [6,7]. If the relation between the heat flux and the wall superheating in pool boiling is known, the nucleation site density can be esti- mated by the following correlation (15]: n = ( q /3)8/3, BHllr;., (3) where BH is the parameter based on the fluid thermo-physical properties: BH = 0.5k1 (-5!......) t/ 4 Prj1/ 12 ( fJ.1 Cp1) t/3 µ1a1 \jj2h12 (4) The wall superheating in Eq. (3), ~Tsat. can be calculated from a nucleate boiling correlation. He.re, the Thom et al. correlation [16) is applied due to its reliabiJity for the water boiling conditions and its simplicity. The following relation is expressed from the Thom correlation: M - ( q )~ sar - 1971.2exp(2p/ 8687 · 103) ' (5) According to Eqs. (3)-(5) and for the heat flux value of 1500 kW/m2 and atmospheric conditions, the nucleation site density is 1.5 x 106 sites/m2 . The other applied approach for the estimation of the nucleation density is based on a theoretical model of bubble growth on the heater's surface [18). The following correlation for the bubble diameter at the moment of separation from the wall was derived in [18): Db = 0.02080L,, (6) where the contact angle 8 is measured in degrees and the capillarity length is given by L, = -.,.---.,.. J (j g(p1 -P2f (7) Eq. (6) was assessed in (19] and its good prediction ability is con- firmed by comparing with experimental results. The bubble depar- ture diameter dependence on the heat flux is not shown by Eq. (6). The weak relation between the heat flux and the diameter of bubble · at detachment is also confirmed in (15). For conditions close to the departure from nucleate boiling, it can be assumed that the whole M hz.o. V S1tvanovlc/lntmiacional]ouma1 of Hear and Moss Transftr 54 (2011) 3296-3303 3299 heater's surfan- is covered with bubbles. which is desrnl>ed with simple geometric relation that one square meter is covered with n bubbles of diameter D,,. as follows: n~ = 1, (8) where n is che density of nudeac1on sites. For a heavily aged surface che comact angle of s• could be assumed I 17). and for atmospheric conditions the bubble departure diameter of 2.6 x 10-4 m is calcu- lated wirh Eq. (6). Then. the density of nucleation sites of 1.5 x 107 m 2 is calculated with Eq. (8). For a fresh heater. a contact angle of 400, ( 171 could be assumed, and the predictions ofEqs. (6)- (8) are Db - 2.08 x I 0~1 m and n • 2.3 x I os sites/m2• Hence, the va- lue of n m 1.5 x 1 a6 sites/m2 predicted wi th Eq. (3) is within here predicted range with Eqs. (6)-(8) for fresh and aged heaters. The experimental Investigation by the high speed infrared imaging of nano-film heaters from behind in the pool boiling in [6,7] showed that the nucleation site densities for aged heaters are not higher than 6 x 105 sites/m2• For fresh heaters. up to 105 sites/m2 of the density of nucleation si tes are measured. These values are lower than values calculated with Eqs. (3)-(5) and Eqs. (6)-(8). Also. experimental observations showed the increase of the density of nucleation sites with the heat flux increase. It should be also noted that greater differences of density of nucleation sites are indicated under nearly the same values of wetting contact an- gles f 6. 7). This contrasts with the previous belief that contact angle is the key, or even sole parameter that determines nucleation site density [11 .21 .221. ln this paper rhe density of nudeat:ion sites is the input param- eter that depends on the state of the heater surface roughness and the condition of its aging. Since the nucleation site density is mod- eled with che zones shown m Fig. 1. whose width is calculated with Eq. (2~ and it is assumed that the complete beater surface is cov- ered with bubbles at the instant of CHF (represenred with Eq. (8)), the width of the zone is equal to the bubble departure diam- eter. i.e. b • D,,. The second parameter that determines the pool boiling dynam- ics is the bubble residence time on the beater's surface. This is the time of bubble growth up to the bubble departure diameter. It can be predicted from the relation between bubble diameter and bub- ble growth time (20): Db = 2 ( yja + y2ja2 + 2/Jja ./(ii, (9) where Ja is the Jacob number: (10) In (20) it is reported that for contact angles between 40" and 900 the parameter y varies from 0.1 up to 0.49. respectively. The empirical parameter P is equal co 6. The correlation (9) is applicable to the general case when the heal 1s transferred LO the growing bubble boch from the heater's surface into the root of a bubble and from the superheared liquid layer around lhe bubble.. By substiruting Eq. (6) into Eq. (9~ the bubble residence time on the heater's surface is derived: (0 02080)2L: ""--~ 4a1 (}!fa+ ../rJi' + 2p)a) 2 (11) as the function ofche wetting contact angle 8 and the Jacob number. The wall superheating in Eq. (10) for the Jacob number is predicted with Eq. (5) derived from the Thom correlation. In this way the bub- ble residence time in Eq. ( 11) is directly dependent on the heat flux, through the Jacob number. and the contact angle. The relation be- tween the bubble residence Lime on the heater's surface and the -com.d atlgle "40 degrees ·, - - coni.:s en 119 s 90 dognles ' ' ' ' ~ r-, r-:.::: r----- ~---- --- O.OE.OO 200 llOO llOO 1000 1200 1400 1500 HNI b. kW""a Fig. 2.. Relauon ~rwttn the bubble RStdence lime on the he;iter's surfact .md the heat Oux. wall heat flux is shown in Fig. 2 for two wetting contact angles. A weak dependence of the bubble residence time on the wetting con- ract angle is shown for higher heat fluxes. In thi.s paper the bubble residence time is determined according to Eq. ( 11) (the values pre- sented in Fig. 2). To summarize. in the numerical investigation presented in this paper, the boundary condilions for the bubble generation are determined with the density of nucleation sites, as a spatial param- eter. which is introduced through the width of lhe zones presented in Fig. 1, and with lhe bubble residence time on the heater's wall, as a temporal parameter. 4. Governing equ~tioos Three-dimensional liquid and vapor two-phase flow in pool boiling is modeled by the cwo fluid model (231. Mass. momentum and energy conservation equations are wnnen for each phase. while interface transfer processes are modeled by uclosure laws". This approach implies non-equilibrium thermal and flow conditions. Conservation equations take the following form in the indicial notation: Mass balance: (12) Momentum balance: a(a.,p.uu) a(lXtPtUt1Ut1) ap ar + ~ - ci:taicj +r1.•p•g1 .. c- 11•cr. - r c)u.tJ+(-1)•+1 F21.1- (13) Energy balance: Parameters u. p. hand Tare nme averaged instantaneous veloc- ity. pressure. enthalpy and temperarure. respectively. The index k is 1 for water and 2 for steam. The source terms for mass. momen- tum and thennal energy conservation are written on the r.h.s. of Eqs. {12}-(14). The intensity of phase transition. which is the mass of evaporation or condensation per unit volume and time, are de- noted with re and re. respectively. The force of vapor and liquid interfadal drag per unit volume In I Cartesian direction is denoted with F21.1· The term qb represents a volumetric heat rate from the 3300 M. Pno. V. Stevonovlc/ lnrt>matronal Journal of Heat and Moss Tran$ftr 54 (201 I ) 3296- 3303 heater's surface co the liquid per unit volume m the randomly cho- sen control volume within each zone or vapor generation if the va- por void fraction is lower than 1. The volume fraction balance is added to the above system of equations as follows: l'.1 + !%2= 1 , (15) Energy equation for lhe heated wall is written as: arw a (k arw) . . pc- =- w - + qh - qb. at ax, ax, (16) where Qh is the volumetric heat source in the wall (for instance elec- tric heating), while Qb is the volumetric heat sink in the control vol- umes on the wall surface due co the bubble rise (ac the same time this value represents the heac source in the fluid control volumes on the wall when rhe bubble rise occurs). Calculation of the interfacial drag force (F21) is a crucial step for the proper prediction of the relative velocities between the vapor and liquid phase, and consequently the void fraction. The interfa- dal drag force per unit volume of computational cell is calculated as 123): (17) where Co is che interfacial drag coefficient. and Dp is the diameter of the dispersed particle. The correlation for the interfacial drag coef- ficient Co is proposed in the following form: (gt:.p) 1/2 Co - 1.487Dp T (1 - oc2)3(1 - 0.7Soc2)2, (18) where the dependence on the mixture void fraction !%2 has rhe same function form as the CAlliARE code correlation (24( for the inter- face friction in the transitional rwo-phase flow patterns. Other interfacial forces, such as the lift force and the virtual mass force are neglected. The intensity of evaporation rate in rwo-phase mixture is calcu- lated with the empirical model that cakes into account the phase change relaxation time t. The intensity of evaporation rate 1s great- er than zero when the liquid enthalpy is greater than the liquid sat- uration enthalpy, i.e. h 1 > h'. and it is predicted as: r. =11.,p1 h, - h' T, h12 ' (19) s. Numerical method The set of balance Eqs. ( 1 2H 14~ (16) is solved by using the control volume based finite difference method. The semi implicit method for pressure-linked equations (SIMPLE) numerical mechod (25) is used for solving the pressure-correction equation derived from the momentum and mass balance equations for the rwo- phase flow conditions ( 26). Three-dimensional flow field is discret- ized in Cartesian co-ordinates. Numerical grid is made from 40 x 50 x 40 control volumes. Numerical grid consists of two pans: the heated wall (40 x 10 x 40 control volumes) and the two phase mixcure (40 x 40 x 40 control volumes). Heated bottom wall is divided into four by four zones. where the locations of vapor generation in the conrrol volumes ac the heater's surface are deter- mined with different values of the random function for each zone. A discretization of partial differential equations is carried out by their integration over control volumes of basic and staggered grids. Fully implicit time integration and the power law numerical scheme are applied. The set of algebraic equations is solved iteratively with the line-by-line three-diagonal-matrix algorithm (TOMA) 1251. The calculation error for every balance equation and every control volume is kept within limits of 10 5 by iterative solution of sets of linear algebraic equations. 6. Results of three-dimensional numerical simulations and discussion Fig. 3 shows three-dimensional view at the void fraction distri- bution in pool boiling for two different values of heat fluxes under atmospheric conditions. The larger void lumps and the higher swell level position are formed in case of the higher heat flux. In case of the lower heat flux the swell level is almost plane. while in case of the higher heat flux a typical bursting of large vapor lumps is observed at the swell level. Also. the formation ofvapor layers on certain parts of the hearer's surface is observed for 1000kW/m2• Fig. 4 shows vapor blankets formation on the heater's surface with the increase of the heat flux. In case of lower heat flux the va- por spots are formed at the nudeation sites, which locations are randomly chosen within hundred ( 10 x 10} cells at the heater sur- face, while wider vapor blankets are formed with the heat flux in- crease. These vapor blankets at the higher heat flux are formed due ro the intensive vapor generation at the randomly chosen nude- y ~z .... 0 ll .§. 08 >-0.04 07 0.6 0.5 0.02 0.4 0.3 0.2 0 I 00 '~"~ .. .'.fJt~4 . I I • I I I t • • y ~z ....... 0.9 .§. 08 >-0.04 0.7 0.6 0.5 0.02 0.4 0.3 02 01 0.0 f1s. l. Cues through the botl1ng pool (bc~l nux 200 ltW/rrr - top, I OOO ltW/rrr2 - bcmom. number ol nudr.itlon site density IS 10' s1tes/m1• bubble residen« time is 13 ms In c.ue of 200 kW/m1 ~nd 5 ms in ase of 1000 kW/m1}. M. Pezo. V. Slevonovic/lnremmionaJ}GUmal of Hear and MW Trorufrr 54 (2011) 3296-3303 3301 0 0 0 Fig. 4 . Void fraction ar the heaters surface (cop) and 0.01 m abovt' the surface (boa:om). under 200 kW/m1 - left and 1000 kW/m?- right (number of nucleation sites 10 both cases is 10" sites/m1, bubble residence time i5 13 ms in use of200 kW/m1 and 5 ms in G1se of 1000 kW/rn1). ation sites, as well as due to the merging of several vapor Jumps into a larger one. The formatted vapor blanket hinders the replen- ishment of the heater's surface with the liquid, which leads to the burnout condition. Experimentally observed [7] and calculated two-phase mixture structures and void fractions in pool boiling prior to burnout are shown in Fig. 5. The pictures at the cop of figure show the forma- tion of vapor blanket at the heater's surface, wit.h vapor plumes ris- ing from the bottom. Both experimental image and numerical prediction show highly distorted swell level. The diagrams at the bottom of Fig. 5 show that the averaged void value at the heater's surface is approximately 0.9., than it drops to the value of nearly 0.4 within the pool boiling mixture and rises above the height of 0.03 m. The results in Fig. 5 demonstrate that the applied numeri- cal methodology is able to predict the two-phase mixture structure in pool boiling. Fig. 6 shows the measured and calculated CHFs versus the nudeation site density. In the experimental investigation pre- sented In 16.7) the burnout is identified through the heater over- heating with its subsequent mechanical failure. In the calculation the CHF is identified by the abrupt wall surface temperature rise. It is shown that here applied numerical methodology is able to pre- dict the experimental observation char the increase of nucleation site density leads to the increased resistance to burnout The bub- ble residence times used in the calculations are chosen according to Fig. 2. Calculated CHF values are within the range of measured ones, showing good agreement for both fresh and aged heaters. The performed calculations are based on the experimental infor- mation about the actual nucleation sire density. This is the neces- sary input parameter that must be known in order to define adequate boundary conditions for tbe bubble generation. The re- sults in Fig. 6 show that CHF increases with the nucleation site den- sity increase. This finding is in contrast with several previous models of boiling crisis I 11 .21 ). which predict the increase of CHF with the decrease of nucleation site density. Obviously the CHF value strongly depends on the surface roughness and aging, i.e. on the nucleation site density, the parameter that is taken into account by the present modeling approach. The Kutateladze- Zuber correlation (Eq. (1 )) gives the CHF value of 1200 kW/mi, the value that is wi thin the range presented in Fig. 6 and corresponds to the heater with the bubble nucleation site density of approximately 45 si tes/cm2• The heat flux dependence on the difference between tbe hea· ter's surface mean temperature and saturation temperature is shown in Fig. 7 for different nucleation site densities up to the point of CHF occurrence. The points on the dashed lines present s teady-state conditions of pool boiling. The relation between the CHF values and the temperature difference is presented with the full line. As presented with a dashed line for a certain adopted con- stant nucleation site density, the temperature dJfference increases M. Puo. v. SlMJ110\llC/ln1tmarional}oumal of Heal and Mass Trans/~ 54 (2011) 3296-3303 E E 0 !" II ::i: 70 60 50 I 40 ...... .... ::x: 30 ~ w ::x: 20 10 0 X • 35 mm 0 0.2 0.4 0.6 0.8 VOID FRACTION 0.9 0.8 0.7 0.6 o.s 0.4 0.3 0.2 0.1 0.0 1.0 'E ~0.04 I- :c C> ii:i 0 .03 :c 0 0 .07 0 .06 0 .05 e ;: 0 .04 :c C> w 0 .03 :c 0 .02 0 .01 0 0 0.9 0.8 0.7 0.8 0.5 0.4 0.3 0.2 0.1 o.o 0 .01 0.02 0.03 0.04 WIDTH (m) ,,,... u / .JI" \ 7 ) ""- ~ 0 .2 0 .4 0 .6 0 .8 VOID FRACTION 1 ~ s. R.ld1oanpl11c ·~ of void fraction 171 (top left) .ind c.ilcu~!rd ~ues ( top right) and corrrsponclms •~• i1wnged void fr.lctJons u a funellon of lletgh1 allow ~ beat~{-•su~ bottom left c.tlcula!rd - bottom nghl). Calculanon is petfonned with 20" 10' nudeahon sates/m'. S x 10-• s ofbubblt' ~ldt'ntt nme and bear Hux of 1000kW/m2• 1eDO· ? 1:100· ~ l; 1100. Ii: ~ .. 900· :r ~ c 700 0 ~ • model 0 0 O experlmenl·fresh heaters 0 experiment-heavily aged heaters 0 • Ii 0 0 0 500+-~~-,....-~~T,~~---......-~~-r-~~~,~~--t 0 w ~ ~ ~ ~ 60 Nucleation Site Oensrty (s1teslcm1) Fig. 6. Comparl~on of measured (71 and calculated CllF values for different nucleation site densities. with the critical heat flux increase. The point on the full line pre- sents the boiling crisis. i.e. the transient condition when the hea- ter's surface temperature abruprly rises. 7. Conclusion The three-- \ 1 --- ~ - I -<>-25 I -=-~~~ 200 ~ $-----;. - --391 ---s1 0 10 20 30 Temperature Difference. TW·Tsat (K) Fig. 7. Critical hear Rux versus excess wall temperature. ~ 50 sure is presented. Two-phase mixture now in pool boiling is simu- lated with the developed two-fluid model and it is coupled with the solving of the transient beat conduction in the heated wall. Ob- tained numerical results are compared with the available experi- mental observations. lr is shown that the vapor blanket forms on the horizontal heater's top surface in cases of high heat nux values. which leads ro the heated wall surface remperarure rise and boiling crisis conditions. The comparison of numerical and experimental results of two-phase mixture void fraction tn pool boiling and oit- M. P-. V Srrwlnovlc/ lntmulliooolJoumal of Htat 1111d Mass Transfer S4 (2011) l296--l3W 3303 1cal heat nwc values shows that the developed modeling approach and numerical procedure are able to simulate complex two-phase Bow conditions of pool boiling and to predict the condition of boil- ing crisis. Dynamics of the bubble generation is mode led through the den- sity of nucleation sites and the bubble residence time on the hea- ter's surface, while the nudeation sites are randomly distributed on the heater's surface. The numerical results show that the in- crease of density of nudeatioo sites provides higher critical heat flwc values. Regarding previous models of boiling crisis published in the literature, which were based on the assumption about the two-phase m1Xture panem at the heater surface in pool boiling at CHF conditions. the approach presented in this paper caJculares the two-phase mixture structure based on the two-fluid model. i.e. based on the solving of the fiow equations. without a need for any assumptions about the two-phase mixture pattern. Acknowledgement This work was supported by the Ministry of Science and Tech- nological Development of the Republic of Serbia (Grant No. 174014). References ( 1 I S. Nulaymu. Mu.mum and minimum v~u~ of be.at transm1n~ from ~ral to bcnhns water under .atmospheric prcssurr. j.Jpn. Soc. Mech. Eng.37 (1934) 367-374. in JaJWlese: S. Nuluy.ama. Maximum and minimum values of he.at 11.ansmined from metal to boiling water under aanospheric pressurr. lnL J. Heat Mass Transfer 9 ( 19GG) 1419- 1433. in English. (21 S.S. Kumcladze. On the transition to film boiling under natural convection, Kotlotufbomenie 3 ( 1948) 10. (3( S.S. K11t.1reladzc. Hydrodynamic model of heat lraMfcr crisis In free-convection bolling.j. Tech. Phys. 20 (11)(!950) 1389- 1392, 141 N. Zuber. On the sabilicy of boiling heal transfer, ASMEJ. llcar Transfer 80 (2) ( 1958)711- 720. ISI J.H. Ucnh.ird. A HCilt TransferTexlbook. Pblop.uon Prtss. Clmbridgc. 2002 PP- 466--474. (6( T.C. TheOl'•nous. J.P. Tu. A.T. Dinb. T.N. Dinh.~ bo1hn1 crisis phenomenon Put I: nucleation and nodeare boiling transfer. Exp. Therm. Fluid Sci. 26 (2002) 775- 792. (71 T.C. lbcofanous. T.N. Dinb. J.P. Tu. A.T. Dinh. The bo1lin1 cns1s phflM>menon Part 11; dryout dyrumics and burnout. Exp. Therm. Auld Sa. 26 (2002) 793- 810. fBl 1.C. 8.l11J, S.H. Chang. W.P. Baek. Visu.aliuuon o( .I pnno ple mtth.anism of a uic.ll be.ir nux in pool boiling. lnL J. Heal Mass Transfer 48 (2005) 5371- 5385. . . (91 HJ. Chung. ll.C. No. A nucleate boilmg 1im11.itlon model for the pred.1coon of pool boiling CHF. lnLJ. Heat Mass Transfer SO (2007 ) 2944- 2951. . . 11 OI V.K. Ohlr, S.P. Liaw, Framework for a unified model for nucleate and trans111on boiling. ASMEJ. Heat Transfer 111 (1989) 739- 746. ( 11 ( N.I, Kolev. How accurately can we predict nucleate boll Ing?, Exp Therm. Fluid Sd. I 0 ( 1995) 370- 378. (121 Y.H. Zhao, T. Masuoka. T. Tsurut:i, Unified thcore1lcal prediction of ful~y developed nucleate boiling and critical he.it flux based on a dynanuc microl.iyer model. lnt. J. Heat Mass Tr.insfer 4 5 (2002) 3189- 3197. ( 131 Y. He. M. Shoji. s. M.iruymu. Numerical S1udy of he.it nux pool boiling but transfer. lnL J. Heat Mass Transfer 44 (2001 ) 2157-2373. f14( z. StOSK. V. Stevanovic. Three-