UNIVERZITET U BEOGRADU MAŠINSKI FAKULTET Sanja S. Milivojević NUMERIČKA SIMULACIJA PROSTIRANJA TEMPERATURSKIH TALASA PRI STRUJANJU NOSIOCA TOPLOTE U SLOŽENIM CEVNIM MREŽAMA doktorska disertacija Beograd, 2013 UNIVERSITY OF BELGRADE FACULTY OF MECHANICAL ENGINEERING Sanja S. Milivojević NUMERICAL SIMULATION OF HOT WATER TEMPERATURE WAVES PROPAGATION IN COMPLEX NETWORKS Doctoral Dissertation Belgrade, 2013 i Komisija za pregled i odbranu Mentor prof. dr Vladimir Stevanović, redovni profesor Univerzitet u Beogradu, Mašinski fakultet Ĉlanovi komisije prof. dr Aleksandar Gajić, redovni profesor Univerzitet u Beogradu, Mašinski fakultet prof. dr Miroslav Stanojević, redovni profesor Univerzitet u Beogradu, Mašinski fakultet prof. dr Branislav Ţivković, vanredni profesor Univerzitet u Beogradu, Mašinski fakultet dr Milan Rajković, nauĉni savetnik Institut za nuklearne nauke ”Vinĉa” Datum odbrane ii PREDGOVOR Ovaj rad je nastao posle višegodišnjeg istraţivanja na Katedri za termoenergetiku, Mašinskog fakulteta Univerziteta u Beogradu. Tema rada je predviĊanje prostiranja temperaturskih talasa u nosiocu toplote u sloţenim cevnim mreţama. PredviĊanje prostiranja talasa temperature je znaĉajno za rad razliĉitih termoenergetskih i termotehniĉkih postrojenja sa stanovišta sigurnosti, pouzdanosti ili ekonomiĉnosti rada. Motiv za ovaj rad je proistekao iz ţelje da se da doprinos optimizaciji prelaznih reţima rada sistema daljinskog grejanja sa stanovišta povećanja energetske efikasnosti snabdevanja potrošaĉa toplotom. TakoĊe, istraţivanje uslova nastanka i dinamike odvijanja prelaznog procesa pri ispadu parnog bloka termoelektrane je od velikog znaĉaja za termohidrauliĉku sigurnost, pošto je pojava termohidrauliĉkog udara, koji pri tome moţe nastati, opasna po sigurnost postrojenja i bezbednost pogonskog osoblja. Zahvaljujem se mentoru prof. dr Vladimiru Stevanoviću na podršci, strpljenju, korisnim savetima i uputstvima koje mi je pruţio u toku izrade disertacije. TakoĊe se zahvaljujem kolegama Blaţenki Maslovarić i Milanu M. Petroviću na saradnji i doprinosu koji su dali u formiranju konaĉne verzije disertacije. Beograd, oktobra, 2013. Sanja S. Milivojević iii NUMERIČKA SIMULACIJA PROSTIRANJA TEMPERATURSKIH TALASA PRI STRUJANJU NOSIOCA TOPLOTE U SLOŽENIM CEVNIM MREŽAMA APSTRAKT U okviru doktorske disertacije je razvijen numeriĉki postupak za proraĉun prostiranja temperaturskih talasa i razdelne površine izmeĊu stuba teĉnosti i parne faze u cevovodima i sloţenim cevnim mreţama. Postupak je zasnovan na rešavanju energijske jednaĉine jednodimenzionalnog strujanja homogenog fluida. Energijska jednaĉina se rešava numeriĉki duţ karakteristiĉnog pravca odreĊenog kretanjem fluidnog delića u vremensko-prostornom koordinatnom sistemu. Postupak numeriĉkog rešavanja je eksplicitan i vremenski korak integracije je odreĊen minimalnim vremenom potrebnim da fluidni delić preĊe rastojanje od poĉetnog poloţaja, koji se u opštem sluĉaju nalazi izmeĊu ĉvorova mreţe, do najbliţeg susednog numeriĉkog ĉvora u pravcu strujanja, u okviru cele mreţe (Courant-ov kriterijum). Primenom Lagrange-ovog interpolacionog polinoma se odreĊuje vrednost temperature ili entalpije u poĉetnom vremenskom trenutku izmeĊu ĉvorova numeriĉke mreţe, na karakteristici koja predstavlja pravac prostiranja fluidnog delića. Razvijene su bilansne jednaĉine koje omogućuju proraĉun graniĉnih uslova, kao što su spoj tri ili više cevi u ĉvoru, toplotno- razmenjivaĉke podstanice kod potrošaĉa i u izvoru toplote. Razvijena metoda je primenjena na simulaciju prelaznog procesa u okviru sistema daljinskog grejanja toplane Zemun nastalog usled promene snage izvora toplote. Razvijeni numeriĉki postupak je primenjen i za simulaciju i analizu hidrauliĉkog udara izazvanog kondenzacijom pare. Neposredni kontakt pothlaĊene teĉnosti i pare u stanju zasićenja dovodi do intenzivne kondenzacije. Brzina kondenzacije zavisi od specifiĉne razdelne površine teĉne i parne faze i koeficijenta prelaza toplote usled kondenzacije. Razdelna površina teĉnosti i pare ima veoma nepravilan oblik. Za vreme hidrauliĉkog udara izazvanog kondenzacijom pare nestacionarni mlazevi teĉnosti i kapi se odvajaju sa ĉela stuba teĉnosti i ukljuĉuju se u zapreminu pare. Ova pojava ukljuĉivanja znaĉajno povećava razdelnu površinu i brzinu kondenzacije. U cilju predviĊanja brzine kondenzacije u ovakvim sloţenim uslovima, razvijena je korelacija izmeĊu ubrzanja ĉela vodenog stuba i proizvoda specifiĉne razdelne površine teĉne i iv parne faze i koeficijenta prelaza toplote pri kondenzaciji. Dvofazni sistem teĉnosti i pare je opisan zakonima odrţanja mase, koliĉine kretanja i energije, koji formiraju sistem kvazilinearnih parcijalnih diferencijalnih jednaĉina hiperboliĉkog tipa. Dobijeni sistem jednaĉina je rešen numeriĉki primenom metode karakteristika. Korišćena su tri karakteristiĉna pravca, dva odreĊena pravcem prostiranja talasa pritiska i treći odreĊen prostiranjem fluidnog delića. Praćenje fluidnog delića i razdelne površine vode i pare izvršeno je rešavanjem energetske jednaĉine u prostoru sa taĉnošću ĉetvrtog reda, dok je prisustvo vode, dvofazne mešavine ili pare odreĊeno vrednošću termodinamiĉkog stepena suvoće. Validacija razvijenog numeriĉkog postupka je sprovedena poreĊenjem numeriĉkih rezultata sa raspoloţivim eksperimentalnim podacima. Razvijeni numeriĉki postupak omogućava pouzdano predviĊanje uslova pojave hidrauliĉkog udara izazvanog kondenzacijom pare, kao i dinamiĉke promene pritiska izazvane udarom. Ključne reči: strujnotermiĉki prelazni procesi, temperaturski talas, metoda karakteristika, hidrauliĉki udar, kondenzacija Naučna oblast: Mašinstvo Uža naučna oblast: Termoenergetika UDK broj: UDK 532.542:621.17:536.7(043.3) v NUMERICAL SIMULATION OF HOT WATER TEMPERATURE WAVES PROPAGATION IN COMPLEX NETWORKS ABSTRACT A high-order accurate numerical method for prediction of thermal transients and propagation of interfacial area between liquid column and steam in pipelines and complex pipe network is presented. It is based on the numerical solution of the transient energy equation for one-dimensional homogeneous fluid flow. Transient one- dimensional energy balance equation is solved along the characteristic path that is determined by the fluid particle motion in time-space coordinate system. Numerical method of solving is explicit and the time step of the integration is determined by the minimal time needed for fluid particle to cross the distance between the initial point, in general case positioned between two adjacent nodes of numerical mesh, and the closest node in the flow direction, in whole mesh (Courant criterion). The initial value of temperature or enthalpy between nodes of numerical mesh, on characteristic path which represents fluid particle propagation direction, is determined by the application of Lagrange's interpolation polynomial. Balance equations for boundary conditions are developed, e.g. junctions of three or more pipelines, heat exchanger substations at costumers and in heat source. Thermal transients caused by an increase and decrease of the heat power plant load are simulated for real operating conditions of the district heating system Zemun. Developed numerical model is also applied for simulation and analysis of condensation induced water hammer. Direct contact of subcooled liquid and saturated steam leads to intensive condensation. Condensation rate is highly influenced by the interfacial area concentration and condensation heat transfer coefficient. The interfacial area has very irregular shape. During the condensation induced water hammer transient liquid jets are formed and liquid column is dissintegrated so that droplets from the head are entrained in steam volume. This droplet entrainment increases interfacial area and condensation rate. Two-phase system of liquid and steam is described by the mass, momentum and enthalpy conservation equations which form a system of quasi linear partial differential equations of the hyperbolic type. Given system is solved by the vi method of characteristics. The explicit time integration is performed along the three characteristic paths, where two of them are determined with the pressure wave propagation and the third by the fluid particle flow. The tracking of the interface between steam and water column head is achieved by the application of the higher order degree numerical scheme to solving of the energy equation. The thermodynamic quality determines the presence of water, two-phase mixture or steam. Validation of the model is conducted through comparison of results of computer simulations with available experimental measurements in literature. Developed numerical method is capable for reliable predictions of conditions which leads to condensation induced water hammer, as well as dynamic change of pressure caused by the condensation induced water hammer. Key words: thermal hydraulic transient, temperature wave, method of characteristics, water hammer, condensation Scientic discipline: Mechanical engineering Scientic subdiscipline: Thermal power engineering UDC number: UDC 532.542:621.17:536.7(043.3) vii SPISAK OZNAKA UPOTREBLJENIH U RADU Latiniĉna slova a ubrzanje razdelne površine teĉne i parne faze, m/s2 ai specifiĉna razdelna površine, m 2 /m 3 A površina popreĉnog preseka cevi, površina popreĉnog preseka mlaza, m2 C + , C - , C P oznaka karakteristika C toplotni kapacitet, J/kg c lokalna brzina zvuka u fluidu, m/s cp specifiĉni toplotni kapacitet pri konstantnom pritisku, J/(kgK) d preĉnik cevi, m D preĉnik mlaza, m DH hidrauliĉki preĉnik cevi, m e energija fluida, J/kg f koeficijent trenja g ubrzanje Zemljine teţe, m/s2 h specifiĉna entalpija, J/kg koeficijent prelaza toplote, W/(m 2 K) Ja Jakob-ov broj L duţina mlaza, duţina cevi, m m maseni protok, kg/s Nu Nusselt-ov broj p pritisak, Pa Pr Prandtl-ov broj q toplotni fluks, W/m 2 Re Reynolds-ov broj Rg gasna konstanta, J/(kgK) r latentna toplota kondenzacije, J/kg St Stanton-ov broj Su Suratman-ov T temperatura, 0 C t vreme, s Δt vremenski korak integracije, s u brzina strujanja fluida, m/s, ui brzina kretanje razdelne površine teĉne i parne faze, m/s V zapremina, m 3 v specifiĉna zapremina, m3/kg W teţinska funkcija We Weber-ov broj X jednaĉina (7.25) x prostorna koordinata, m xt termodinamiĉki stepen suvoće Δx prostorni korak integracije, m Y jednaĉina (7.26) Z jednaĉina (7.27) viii Grĉka slova  izraz (7.48), zapreminski udeo pare  izraz (7.50) c brzina kondenzacije, kg/(m 3 s)  izraz (7.49)  izraz (7.51)  malo rastojanje, m  ugao nagiba strujnog kanala u odnosu na horizontalu, o  eksponent adijabate  toplotna provodnost, W/(mK) talasna duţina, m   dinamiĉka viskoznost, Pas  kinematska viskoznost, m 2 /s  gustina, kg/m 3   površinski napon, N/m  napon smicanja, N/m 2 * z  sila trenja na zidu po jedinici zapremine strujnog kanala, N/m 3  koeficijenti u jednaĉini (7.28) Indeksi A, B, C, D, L, P, R oznaka pripadnosti odreĊenoj taĉki - Slika 3.1 C cev c, kond kondenzacija DB Dittus - Boelter-ova korelacija i razdelna površina (eng. interface) j, j+1 oznaka cevi u spoju K kap pothlaĊene teĉnosti cr kritiĉna vrednost nt nestacionarno trenje o okolina p para, povrat r razvod P (C P) ĉestica (eng. particle) res stanje u rezervoaru sat saturacija (zasićenje) ST stub (pothlaĊene) teĉnosti t nestacionarnost (eng. transient), teĉnost v voda ’ stanje zasićene vode ’’ stanje zasićene pare 0 poĉetni trenutak 1 teĉna faza 2 parna faza ix SADRŽAJ: 1. UVOD 1 2. DOSADAŠNJA ISTRAŢIVANJA PROSTIRANJA TEMPERATURSKIH TALASA U CEVOVODIMA I CEVNIM MREŢAMA 8 2.1 Prostiranje temperaturskog talasa u sistemu daljinskog grejanja 8 2.2 Hidrauliĉki udar izazvan kondenzacijom pare 19 3. NUMERIĈKE METODE ZA PRAĆENJE PROSTIRANJA FRONTA TALASA 41 3.1 Lax-ova metoda 42 3.2 Lax-Wendroff-ova metoda 43 3.3 MacCormack-ova metoda 45 3.4 Rusanov (Burstein - Mirin) metoda 46 3.5 Warming - Kutler - Lomax metoda 47 3.6 Metode trećeg reda sa podešavanjem parametra  48 3.7 Implicitne metode 49 3.8 Šema Godunov-a 52 3.9 Roe-ova šema 55 4. NUMERIĈKA METODA ZA PRAĆENJE PROSTIRANJA FRONTA SKALARNE PROMENLЈIVE DUŢ KARAKTERISTIĈNOG PRAVCA 60 4.1 Graniĉni uslovi 67 4.2 Hidrauliĉki model 69 5. REZULTATI SNIMANJA PRELAZNOG PROCESA U SISTEMU DALJINSKOG GREJANJA 70 5.1 Prelazni temperaturski reţim u sistemu daljinskog grejanja toplane „Zemun“ 70 6. NUMERIĈKE SIMULACIJE PRELAZNIH REŢIMA RADA SISTEMA DALЈINSKOG GREJANJA 78 7. TERMOHIDRAULIĈKI MODEL ZA PRORAĈUN HIDRAULIĈKOG UDARA IZAZVANOG KONDENZACIJOM PARE 81 7.1 Matematiĉki model jednodimenzionalnog nestacionarnog strujanja homogenog stišljivog fluida 81 x 7.1.1. Formiranje sistema kvazilinearnih parcijalnih diferencijalnih jednaĉina hiperboliĉkog tipa 86 7.2 Numeriĉko rešavanje jednodimenzionalnog nestacionarnog strujanja homogenog stišljivog fluida primenom metode karakteristika 87 7.2.1 Formiranje sistema obiĉnih diferencijalnih jednaĉina sa totalnim diferencijalima – formiranje jednaĉina karakteristika 89 7.2.2 Formiranje sistema diferencnih jednaĉina 91 7.2.3 Rešavanje sistema diferencnih jednaĉina primenom metode karakteristika 91 7.3 Homogeni model dvofaznog strujanja 96 7.4 OdreĊivanje brzine zvuka 97 7.5 Uticaj nestacionarnog trenja 99 7.6 Matematiĉko definisanje graniĉnih uslova 103 7.6.1. Graniĉni uslov na poĉetku cevi 104 7.6.1.1 Rezervoar na poĉetku cevi 104 7.6.2. Graniĉni uslov na kraju cevi 105 7.6.2.1 Zatvoreni kraj cevi na kraju cevi 106 7.6.2.2 Rezervoar na kraju cevi 106 7.6.3. Graniĉni uslov u spoju dve cevi istog popreĉnog preseka 106 8. KOMPJUTERSKI PROGRAM ZA SIMULACIJU HIDRAULIĈKOG UDARA IZAZVANOG KONDENZACIJOM PARE 109 8.1. Opis programa 109 8.2. Modeliranje strujnog prostora 115 9. NUMERIĈKA SIMULACIJA KRETANJA RAZDELNE POVRŠINE I HIDRAULIĈKOG UDARA IZAZVANOG KONDENZACIJOM PARE 117 9.1 Kretanje razdelne površine teĉne i parne faze u oscilujućem manometru 118 9.1.1 Diferencijalna jenaĉina kretanja razdelne površine u oscilujućem manometru 120 9.1.2 Grafiĉki prikaz dobijenih rezultata pri numeriĉkoj simulaciji kretanja razdelne površine u oscilujućem manometru 122 9.2 Hidrauliĉki udar usled kondenzacije pare na stubu pothlaĊene vode 126 xi 9.2.1 Diferencijalna jednaĉina kretanja razdelne površine u stubu teĉnosti promenljive visine u vertikalnoj cevi 127 9.2.2 Grafiĉki prikaz dobijenih rezultata pri numeriĉkoj simulaciji hidrauliĉkog udara usled kondenzacije pare na stubu pothlaĊene vode 130 9.3 Hidrauliĉki udar u vertikalnoj cevi za ispuštanje pare u bazen sa pothlaĊenom vodom 142 9.3.1 Grafiĉki prikaz dobijenih rezultata pri numeriĉkoj simulaciji hidrauliĉkog udara u vretikalnoj cevi za ispuštanje pare u bazen sa pothlaĊenom vodom 143 9.4 Analitiĉko odreĊivanje proizvoda koeficijenta prelaza toplote i razdelne površine pri kondenzaciji pare na pothlaĊenoj vodi pri hidrauliĉkom udaru izazvanom kondenzacijom pare 146 9.4.1 OdreĊivanje koeficijenta prelaza toplote i razdelne površine pri kondenzaciji pare na ĉelu stuba teĉnosti 151 9.4.2 Prelaz toplote pri kondenzaciji pare na kapima pothlaĊene teĉnosti otkinutim sa ĉela stuba teĉnosti i ukljuĉenim u parni prostor 152 9.4.3 Kondenzacija na mlazu pothlaĊene teĉnosti otkinutom sa ĉela stuba teĉnosti 155 10. ZAKLJUĈAK 160 LITERATURA 163 1 1. UVOD Tokom prelaznih procesa u termoenergetskim postrojenjima dolazi do prostiranja talasa temperature u nosiocu toplote ili radnom fluidu. Temperaturski talasi nastaju pri nestacionarnom strujanju fluida i usled promene toplotnog fluksa na površinama za razmenu toplote u izvoru toplote (kao što su vrelovodeni kotlovi, nuklearni reaktori i generatori pare), u ponoru toplote (različite vrste razmenjivača toplote i kondenzatora) ili pri nestacionarnoj razmeni toplote između fluida i okoline kroz zidove cevovoda i izolaciju. Prostiranje talasa temperature je određeno strujanjem fluidnih delića, provođenjem i turbulentnom razmenom toplote u pravcu strujanja i razmenom toplote kroz zidove strujnog kanala. Tačno predviđanje prostiranja talasa temperature je značajno za rad različitih termoenergetskih i termotehničkih postrojenja sa stanovišta sigurnosti, pouzdanosti ili ekonomičnosti rada. Potreba za predviđanjem kretanja temperaturskih talasa u strujnim kanalima postoji, na primer, u složenim cevnim mrežama sistema daljinskog grejanja, u cevovodima kondenzata i napojne vode u termoelektranama, u nuklearnim sistemima za proizvodnju pare itd. Predviđanje prostiranja temperaturskih talasa je značajno za optimizaciju prelaznih režima rada sistema daljinskog grejanja sa stanovišta povećanja energetske efikasnosti snabdevanja potrošača toplotom. Do prelaznih režima u radu toplovoda dolazi usled promene spoljne temperature vazduha, promenljive potrošnje tople vode u sistemima koji obezbeđuju potrošnu toplu vodu i pri zaustavljanju ili pokretanju isporuke toplote. Na efikasnost snabdevanja toplotom u prelaznim režimima utiču vremenska usklađenost promene snage toplotnog izvora sa zahtevima potrošača, potrebni vremenski interval za transport toplote od izvora do pojedinih toplotnih podstanica, akumulacija toplote u okviru cevne mreže i gubici toplote u okolinu. S obzirom na složenost toplovodnih mreža, dinamičke promene potrošnje i karakteristike procesa transporta toplote, predviđanje prelaznih režima je jedino moguće primenom odgovarajućih numeričkih modela i kompjuterskih programa. U složenim cevnim mrežama toplovoda potrebno je predvideti vreme za koje temperaturski talas putuje od izvora toplote do potrošača pri promeni snage izvora ili 2 potrošnje, što značajno utiče na ekonomičnost pogona i kvalitet snabdevanja potrošača toplotom (Gabrielaitiene i dr., 2007; Stevanović i dr., 2007a; Stevanović i dr., 2007b; Stevanović i dr., 2009). Modeli prelaznih procesa u sistemima daljinskog grejanja su zasnovani na stacionarnim hidrauličkim proračunima raspodele protoka i promene pritiska u okviru složenih cevnih mreža i rešavanju jednodimenzionalne nestacionarne energijske jednačine (Kunz i dr., 1991; Gabrielaitiene i dr., 2007). S obzirom na složenost cevnih mreža sistema daljinskog grejanja, razvijeni su i modeli temperaturskih prelaznih procesa kod kojih se veći broj deonica sistema objedinjuje u okviru jedne ekvivalentne deonice (Larsen i dr., 2004). Na taj način se pojednostavljuje model i ubrzava simulacija. Razvijeni modeli su testirani za uslove prostiranja temperaturskih talasa u pojedinim delovima sistema daljinskog grejanja (Kunz i dr., 1991; Gabrielaitiene i dr., 2007; Larsen i dr., 2004). U toku je primena numeričkih modela za optimizaciju prelaznih režima u okviru nekoliko sistema daljinskog grejanja u Danskoj, pri čemu rezultati pokazuju uštede od najmanje 2% u ukupnoj potrošnji energije goriva (7 technologies, 2007). Drugi primer prelaznog procesa u energetskom postrojenju, čija dinamika odvijanja zavisi od prostiranja temperaturskih talasa u fluidu, jeste ispad parnog bloka termoelektrane, pri kome dolazi do složenih termičkih i hidrauličkih pojava u sistemu napojne vode i cevnom sistemu kotla. Kod parnih blokova posebna pažnja se mora posvetiti termohidrauličkoj sigurnosti, pre svega zaštiti od pojave termohidrauličkog udara, pošto je ova pojava posebno opasna po sigurnost postrojenja i bezbednost pogonskog osoblјa (Prica i dr., 2004; Prica i dr., 2009; Prica i dr., 2011; Milivojević i dr., 2013). Do termohidrauličkog udara dolazi pri kontaktu pare i pothlađene vode u cevi ili sudu pod pritiskom. Kontakt pare i pothlađene vode dovodi do intenzivne kondenzacije, praćene padom pritiska u zapremini pare koja se kondenzuje i ubrzanjem okolne mase vode ka zapremini pare, što dovodi do pojave hidrauličkog udara i oštećenja postrojenja. Hidraulički udar izazvan kondenzacijom pare je brza pojava koja može izazvati značajnu materijalnu štetu i ozbiljne povrede pogonskog osoblja u svim sistemima u kojima para i pothlađena tečnost dolaze u neposredan kontakt, kao što su parni blokovi 3 termoelektrana, sistemi daljinskog grejanja, procesna i hemijska industrija, kao i industrija vodosnabdevanja. Hidraulički udar izazvan kondenzacijom pare može nastati i kao posledica lošeg projektovanja, pogrešnog vođenja elektrane pri radu ili pri incidentnim okolnostima. Engleski naziv ove pojave je "condensation induced water hammer" (hidraulički udar izazvan kondenzacijom pare). Sam naziv "water hammer" (vodeni čekić) podrazumeva sve nestacionarne promene uzrokovane promenom pritiska (Fox, 1977). Naziv vodeni čekić potiče od zvuka koji nastaje pri prolasku talasa pritiska kroz cevovod, a koji podseća na zvuk koj nastaje pri udaru čekićem o cevovod. U isparivačkim kanalima generatora pare je potrebno tačno odrediti nestacionarne promene temperature radnog fluida u cilju predviđanja kretanja granice ključanja u prelaznim režimima izazvanim promenom snage ili protoka radnog fluida (Stevanović i Jovanović, 2000). Od položaja granice ključanja zavisi promena koeficijenta razmene toplote duž zagrejačkog i isparivačkog kanala, raspodela toplotnog fluksa, stabilnost strujanja u isparivačkom kanalu, toplotna snaga nuklearnog reaktora s obzirom na povratnu spregu kojom gustina nosioca toplote utiče na reaktivnost (Lahey i Moody, 1979), kao i drugi efekti od značaja za siguran i pouzdan rad generatora pare. Do sada je razvijen veći broj numeričkih postupaka za rešavanje jednačina održanja u fluidnom toku. Jedna od prvih metoda za rešavanje jednačina strujanja je metoda karakteristika. Ona potencijalno daje najtačnije rezultate (Wulf, 1987), ali do sada je uglavnom korišćena za rešavanje izotermskih strujanja koja se opisuju samo jednačinama konzervacije mase i količine kretanja (Streeter i Wylie, 1967). Proširenje primene metode karakteristika i na rešavanje energijske jednačine je urađeno uvođenjem trećeg karakterističnog pravca određenog kretanjem fluidnog delića u prostorno vremenskom koordinatnom sistemu (Stoop i dr., 1985; Stevanović i dr., 1994). Međutim, numeričke šeme prvog reda, kao što je standardna metoda karakteristika, dovode do numeričke difuzije, tako da rešenje talasa temperature prikazuje položaj talasa na znatno većoj dužini od one u realnoj pojavi. U cilju prevazilaženja numeričke difuzije razvijene su numeričke šeme višeg reda tačnosti, kao što su kvadratna uzvodna interpolacija za konvektivnu kinetiku - QUICK šema (Versteeg i Malalasekera, 1995), koja je drugog reda tačnosti. Kao nedostatak 4 numeričkih šema višeg reda tačnosti javljaju se oscilatorne promene u okolini maksimalnih i minimalnih vrednosti sračunate veličine. Iz potrebe za rešavanjem prethodno navedenih problema, u ovoj disertaciji je primenjena nova numerička metoda za simulaciju prostiranja temperaturskih talasa pri strujanju nosioca toplote u složenim cevnim mrežama u jednofaznim i dvofaznim sistemima. Naime, jednačina konzervacije energije se rešava duž karakterističnog pravca u vremensko-prostornom koordinatnom sistemu, određenog prostiranjem fluidnog delića, pri čemu se početna vrednost energetskog parametra (temperature ili entalpije), koja se u opštem slučaju nalazi između fiksnih čvorova numeričke mreže, određuje primenom Lagrange-ovog polinoma trećeg reda. Primenjena metoda je četvrtog reda tačnosti, što omogućava efikasno i tačno rešavanje nestacionarne jednačine konzervacije energije, uz praktično eliminisanje numeričke difuzije. U cilju simulacije strujanja u složenim mrežama razvijeni su odgovarajući modeli graničnih uslova, kao što su spoj dve ili više cevi, isticanje iz cevovoda, spoj sa rezervoarom ili akumulatorom toplote, modeli razmenjivača toplote sa usrednjenim parametrima, pumpi i pumpnih stanica i slično. U cilju rešavanja realnih problema iz inženjerske prakse razvijen je odgovarajući kompjuterski program, sa algoritmom za definisanje konfiguracije mreže i scenarija prelaznih režima. Efikasnost postupka se ogleda u njegovoj stabilnosti pri rešavanju strujanja u složenim cevnim mrežama i zadovoljavajućem vremenu rešavanja koje omogućava sprovođenje simulacija u realnom vremenu ili vremenu kraćem od realnog. Ova disertacija je podeljena u deset poglavlja. U drugom poglavlju su prikazana dosadašnja istraživanja prostiranja temperaturskih talasa kroz složene cevne mreže i istraživanja hidrauličkog udara izazvanog intenzivnom kondenzacijom pare, koji je upravo određen dinamikom kretanja faza različitih temperatura (pothlađene tečnosti i zasićene ili pregrejane pare). U poglavlju tri su opisane numeričke metode višeg reda tačnosti za praćenje kretanja fronta talasa, zajedno sa poređenjem rezultata dobijenih njihovom primenom i analitički dobijenih rezultata. 5 U poglavlju četiri dat je opis numeričke metode višeg reda tačnosti, koja je u ovoj disertaciji primenjena za predviđanje prostiranja talasa temperature u složenim cevnim mrežama. Predstavljeno je numeričko rešavanje energijske jednačine duž karakterističnog pravca definisanog kretanjem fluidnog delića. Postupak numeričkog rešavanja je eksplicitan i vremenski korak integracije je određen minimalnim vremenom potrebnim da fluidni delić pređe rastojanje od početnog položaja između dva susedna numerička čvora do sledećeg susednog čvora unutar cevi u okviru cele mreže (Courant- ov kriterijum). Takođe, prikazana je primena Lagrange-ovog interpolacionog polinoma kojim se određuje vrednost skalarne promenljive veličine u početnom vremenskom trenutku između čvorova numeričke mreže, na karakteristici koja predstavlja pravac prostiranja fluidnog delića. Razvijene su bilansne jednačine koje omogućuju proračun graničnih uslova, kao što je spoj više cevi u čvoru, toplotno-razmenjivačke podstanice kod potrošača i u izvoru toplote. Hidraulički proračun je zasnovan na efikasnoj numeričkoj metodi bilansiranja promene pritiska po prstenovima cevne mreže (Stevanović i dr., 2006b; Stevanović i dr., 2007c). U poglavlju pet su prikazani rezultati snimanja prelaznog procesa u toplotno- razmenjivačkim stanicama kod potrošača u sistemu daljinskog grejanja, koji su korišćeni za validaciju u ovoj tezi razvijenog numeričkog postupka za predviđanje prostiranja talasa temperature u složenim cevnim mrežama. Merenja dinamičkih promena temperatura i protoka su izvršena u okviru sistema daljinskog grejanja toplane ”Zemun“ pri promeni snage izvora toplote (Stevanović i dr., 2006a). Validacija programa za realne uslove rada na osnovu izvršenih merenja je data u poglavlju šest. Dobijeno je dobro slaganje rezultata kompjuterske simulacije prostiranja temperaturskog talasa, usled prelaznog režima rada sa izmerenim vrednostima. U cilju predviđanja prostiranja temperaturskih talasa u dvofaznom toku tečne i parne faze, u poglavlju sedam je opisan matematički model jednodimenzionalnog nestacionarnog strujanja homogenog stišljivog fluida i način numeričkog rešavanja sistema kvazilinearnih parcijalnih diferencijalnih jednačina hiperboličkog tipa primenom metode karakteristika. Zatim je predstavljen homogeni model dvofaznog toka, navedene su korelacije za određivanje lokalne brzine zvuka u fluidu i korelacije 6 kojima se uzima u obzir nestacionarno trenje, određen je prelaz toplote pri kondenzaciji za slučaj termohidrauličkog udara. Takođe, definisani su odgovarajući granični uslovi na krajevima strujnog prostora. Rešavanjem energetske jednačine razvijenim numeričkim postupkom, u kojoj je zavisno promenljiva entalpija, omogućeno je za proračun temperaturskih talasa u dvofaznom toku predviđanje kako temperaturskih talasa, tako i kretanje razdelne površine između faza u dvofaznom toku. Poglavlje osam sadrži opis razvijenog programa sa algoritmima glavnog programa i nekih podprograma. U ovom poglavlju je predstavljen i način modeliranja strujnog prostora i pripreme ulaznih podataka za program. U poglavlju devet su prikazane primene razvijenog programa na simulaciju kretanja razdelne površine tečne i parne faze i simulaciju hidrauličkog udara izazvanog kondenzacijom pare na primeru nekoliko eksperimenata. Dato je poređenje dobijenih rezultata sa raspoloživim izmerenim vrednostima, poređenje dobijenih rezultata sa rezultatima drugih kompjuterskih programa i poređenje sa rezultatima mogućih analitičkih rešenja. Predstavljene su mogućnosti programa i njegova tačnost pri praćenju kretanja razdelne površine tečne i parne faze, određivanju maksimalne vrednosti pritiska i prostiranja talasa pritiska usled hidrauličkog udara u dvofaznom fluidu. Deseto poglavlje je zaključak, u okviru koga je istaknut značaj istraživanja prostiranja temperaturskih talasa u jednofaznom i dvofaznom toku za optimizaciju rada prelaznih režima u sistemu daljinskog grejanja i analize udesnih stanja u termoenergetskim postrojenjima izazvanih hidrauličkim udarom. Navedeni su doprinosi koji su ostvareni izradom ove disertacije. Validacijom razvijenog numeričkog postupka za uslove prostiranja temperaturskih talasa u složenoj cevnoj mreži pokazano je da dobijeni rezultati omogućuju određivanje vremena potrebnog da temperaturski talas pređe rastojanje od izvora toplote do najudalјenijeg potrošača, odnosno vremensko kašnjenje promene snage potrošnje u odnosu na promenu snage izvora toplote. Takođe, rezultati dobijeni primenom razvijenog kompjuterskog programa na simulaciju kretanja razdelne površine i hidrauličkog udara izazvanog kondenzacijom pare u različitim geometrijama strujnih kanala i za različite nivoe pritisaka u sistemu potvrđuju da je 7 sprovedeni proračun stabilan i da razvijeni kompjuterski program omogućava pouzdano predviđanje uslova nastanka hidrauličkog udara izazvanog intenzivnom kondenzacijom pare, kao i dinamičke promene pritiska izazvane udarom. 8 2. DOSADAŠNJA ISTRAŽIVANJA PROSTIRANJA TEMPERATURSKIH TALASA U CEVOVODIMA I CEVNIM MREŽAMA U analizama i pri predviđanju ponašanja termoenergetskih postrojenja u prelaznim pogonskim uslovima javlja se potreba za tačnim proračunom prostiranja talasa temperature. U ovom poglavlju su prikazana dosadašnja istraživanja prelaznih režima rada sistema daljinskog grejanja, koja zavise od prostiranja talasa temperature, kao i značaj optimizacije tih pogonskih uslova u cilju povećanja energetske efikasnosti transporta toplote. Takođe, prikazana su i istraživanja hidrauličkog udara izazvanog kondenzacijom pare u dvofaznom stujanju, koja su značajna za sigurnost termoenergetskih postrojenja, a čija dinamika odvijanja zavisi od prostiranja talasa temperature. 2.1 Prostiranje temperaturskog talasa u sistemu daljinskog grejanja Sistemi daljinskog grejanja predstavljaju jedan od najefikasnijih načina grejanja stambenih i poslovnih objekata i proizvodnje sanitarne tople vode. Njihovom primenom je omogućeno postizanje veće efikasnosti i prema tome nižih troškove grejanja zbog proizvodnje toplote u okviru centralizovanih izvora toplote termoelektrana - toplana i toplana. Efikasan vid snabdevanja toplotom je u kogeneraciji - spregnutoj proizvodnji toplote i električne energije. Sa ekološkog stanovišta kogeneracija ili drugi vid kombinovane proizvodnje toplote omogućuje efikasnije mere prečišćavanja dimnih gasova i efikasniju potrošnju goriva. Distribucijom toplote iz jedinica za centralizovanu proizvodnju toplote smanjuje se zagađenje vazduha, jer se toplota proizvodi u jednom izvoru umesto u više pojedinačnih izvora toplote. Imajući u vidu da je u našoj zemlji ugalj značajan energent, centralizovana proizvodnja toplote poboljšava njegovu racionalnu eksploataciju. Pogonsko osoblje sistema daljinskog grejanja često i po nekoliko puta dnevno mora da donese odluke o načinu vođenja sistema. Ove odluke uključuju određivanje količine toplote koju je potrebno proizvesti, startovanje i zaustavljanje jedinica za proizvodnju toplote, startovanje i zaustavljanje pumpi, vrednost temeperature odlazne 9 tople vode i masenog protoka tople vode. Dok se uticaj donetih odluka na vođenje sistema oseti kod potrošača potrebno je da protekne izvesno vreme, koje može da bude od nekoliko sati pa do više od jednog dana, zavisno od veličine i složenosti konfiguracije mreže sistema daljinskog grejanja. Zbog toga je osoblju vrlo teško da odredi kada i koje radnje treba da primeni. Dakle, osoblju je potrebna informacija o vremenu potrebnom da potrošači osete promenu temperature, protoka, odnosno količine toplote koja im se isporučuje. Zato je pri vođenju procesa od velike pomoći postojanje modela koji bi mogao da simulira posledice promene parametara u toplani. Optimizacija sistema daljinskog grejanja zahteva zadovoljenje potreba potrošača uz minimalno moguće troškove. Ovaj problem je dinamičke prirode. Osnova za optimizaciju vođenja postrojenja je postojanje dobre dokumentacije o procesu, postojanje on-line merenja, pristup neophodnim informacijama i razumevanje različitih podprocesa. Takođe su potrebne metode i strategija za određivanje optimalne radne tačke na osnovu dostupnih podataka i matematičkog opisa procesa. U literaturi postoji nekoliko različitih predloga za smanjenje troškova vođenja sistema daljinskog grejanja, ali se samo u nekoliko referenci mogu naći studije zajedno sa izvršenim testiranjem predloženog modela. Usled velike složenosti problema većina referenci sugeriše metode koje se odnose na posmatranje samo jednog dela celokupnog problema, dok se za ostale delove pretpostavlja da se ponašaju po nekom unapred predviđenom scenariju ili se jednostavno zanemare (Benonysson i dr., 1995). Predloženi modeli, se mogu prema tome podeliti u tri kategorije:  Određivanje optimalnog opterećenja u svim izvorima toplote, toplani u sistemu daljinskog grejanja i toplotno razmenjivačkim stanicama kod potrošača, u nekim slučajevima se razmatra i akumulacija toplote, koristeći konstantnu ili unapred određenu temperaturu razvoda (dinamika mreže grejanja se ne razmatra).  Snižavanje temperature razvoda bez uzimanja u obzir mogućnosti akumulacije toplote u grejnoj mreži.  Potpuna dinamička optimizacija sistema podrazumeva određivanje i optimalne temperature razvoda i optimalnog toplotnog opterećenja istovremeno. Do sada nije pronađeno sveobuhvatno rešenje iako je ovaj problem razmatrao veći broj autora (Benonysson i dr., 1995). 10 Dinamičko ponašanje potrošača i dinamika distribucije u mreži veoma utiču na rad sistema daljinskog grejanja i posledica su velikog vremena odziva mreže, akumulacije toplote u mreži i gubitaka toplote u okolinu. Za planiranje i optimizaciju rada je prema tome neophodno imati odgovarajući model za simulaciju ponašanja potrošača i konfiguraciju mreže sistema grejanja. U (Benonysson i dr., 1995) je predstavljen takozvani metod čvorova. Ovaj metod se može primeniti na simulaciju masenog protoka i temperaturskog talasa u sistemu centralnog grejanja koji se menjaju kao posledica promene potreba potrošača za toplotom i promene temperature odlazne tople vode. Da bi se primetila promena temperature tople vode potrebno je onoliko vremena koliko je potrebno temperaturskom talasu da stigne od toplane do potrošača, a to zavisno od udaljenosti potrošača od toplane može da potraje i nekoliko sati. Vreme prostiranja temperaturskog talasa ne zavisi samo od udaljenosti potrošača već i od toplotnog kapaciteta materijala cevi, najčešće čelika, koji i prima i predaje toplotu vodi u sistemu zavisno od promene njene temperature. Toplotni gubici su direktno proporcionalni razlici temperatura tople vode i temperature tla (imajući u vidu da su cevi sistema daljinskog grejanja uglavnom ukopane u zemlju). Prema tome gubici toplote se povećavaju kako temperatura vode raste ili kako temperatura tla pada. Stvarni toplotni gubici po jedinici dužine zavise i od izolacije cevi i prečnika cevi. Imajući u vidu da se temperatura tla kroz koje ne prolaze cevi menja sporo u toku godine i da prati promenu temperature okolnog vazduha sa izvesnim zakašnjenjem, promene u toplotnim gubicima posmatrane u kratkom vremenskom periodu uglavnom zavise od promena temperature tople vode. Gubici usled trenja u cevima pretvaraju energiju pumpanja u toplotnu energiju težeći da povise temperaturu vode u cevnoj mreži sistema, pojava poznata kao disipacija energije. U većini slučajeva disipacija ima neznatan značaj, ali u slučajevima gde je brzina strujanja vode relativno velika, količina toplote oslobođena na ovaj način može biti iste veličine kao i toplotni gubici sa cevovoda u tlo. Promene pritiska (i masenog protoka) kroz cevnu mrežu se prostiru oko 1000 puta brže nego promene temperature, s obzirom da se talas pritiska prostire brzinom zvuka kroz vodu – oko 1000 m/s, dok se temperaturski talas prostire brzinom bliskoj 11 brzini strujanja tople vode. Ovo navodi na zaključak da je sa stanovišta optimizacije vođenja postrojenja dinamika strujanja u cevnoj mreži od manjeg uticaja u odnosu na dinamiku promene temperatura. U metodi čvorova, mreža sistema daljinskog grejanja je predstavljena određenim brojem čvorova, njihovim vezama i odgovarajućim tehničkim podacima, kao što su toplotni kapaciteti i prečnici cevi. Princip metode čvorova je da prati koliko vremena je masenom deliću potrebno da stigne od jednog do drugog čvora. Na osnovu podataka o vrednosti temperatura u pojedinim čvorovima u toku vremena, izračunata je temperatura vode na ulazu u cev, a na osnovu ovoga i toplotnih gubitaka određena je temperatura vode na izlazu iz cevi. Na ovaj način temperature u svim čvorovima mogu biti određene za svaki vremenski trenutak. Poznavajući vrednosti temperatura u pojedinim čvorovima mreže (uključujući temperaturu odlazne tople vode iz toplane) u bilo kom vremenskom trenutku i maseni protok između čvorova, moguće je odrediti temperaturu u narednom vremenskom trenutku. Ovako je dalje moguće odrediti temperature u svim čvorovima u svim vremenskim trenucima. Sledeći korak je određivanje protoka vode. Pretpostavimo da je potrošač povezan na mrežu daljinskog grejanja preko razmenjivača toplote. Potreba potrošača je zavisna od temperature na sekundarnoj strani, temperature povrata i toplotnog opterećenja. Putem merenja je ustanovljena sledeća veza između odlazne i povratne temperature vode na primarnoj strani i veličina na sekundarnoj strani 1 1 2( ) ,qp s p sTp Tp k Tr Tr k Q     (2.1) gde su pTp i sTp temperature povrata sa primarne i sekundarne strane ( 0 C), pTr i sTr temperature razvoda sa primarne i sekundarne strane ( 0 C), Q toplotno opterećenje (W), a 1, 2,k k q konstante (-). Primenom jednačine (2.1) određuje se temperatura povrata i na osnovu nje maseni protok vode na primarnoj strani. Da bi se primenio metod čvorova korišćena je temperatura razvoda na primarnoj strani za prethodni vremenski trenutak, t- 1 , 1p tTr  , pa je tako maseni protok vode na primarnoj strani u trenutku t određen kao , 1 1 1 2, 1 , , , [(1 ) ] t p t q pv p t s t s t t Q m c k Tr Tp k Tr k Q       (2.2) 12 gde je vp c specifični toplotni kapacitet vode (kJ/(kgK)), a t vreme (s). Koristeći relaciju (2.2) uz poznate dimenzije cevi i vrednosti temperatura na primarnoj i sekundarnoj strani se mogu izračunati maseni protok i brzine u strujnom toku. Model mreže i potrošača može da se predstavi koristeći sledeće pretpostavke: da je poznata geometrija mreže, da je za svakog potrošača u svakom vremenskom trenutku poznato toplotno opterećenje na sekundarnoj strani, temperatura vode koja ide ka potrošaču (odnosno temperaturu razvoda) i temperatura povratne vode, da je u vremenskom trenutku t=1, poznata temperatura i vrednosti temperatura u prethodnim vremenskim trenucima u svim tačkama mreže (izuzev u toplani) kao i vrednosti masenog protoka u svim cevima mreže. Na osnovu ovih pretpostavki se bira temperatura odlazne tople vode toplane tTr , 1...t T . Onda je moguće odrediti temperature i masene protoke u svim tačkama u bilo kom vremenskom koraku koristeći gore navedeni metod i na osnovu izračunatih protoka izračunati pad pritiska. Ovaj model se može primeniti na različite sisteme koji služe centralizovanom snabdevanju toplotom: toplane, postrojenja sa protivpritisnim turbinama ili postrojenja sa turbinama sa oduzimanjem pare. Na Slici 2.1 je šematski predstavljen iterativni postupak određivanja optimalne temperature odlazne tople vode iz toplane. Dinamičke simulacije prostiranja temperaturskog talasa i akumulacije toplote u sistemu daljinskog grejanja omogućuju istraživanja različitih radnih režima i razvoj modela optimizacije, koji će omogućiti efikasnije vođenje postrojenja. Dinamičkom kontrolom temperature u mreži sistema daljinskog grejanja moguće je ostvariti ekonomske uštede. Postoje mnoge razvijene metode za numeričku simulaciju dinamičkog ponašanja složenih cevnih mreža. Problem koji se javlja je validacija primenjenih modela. Validacija modela toplotne dinamike sistema daljinskog grejanja se obično vrši za jednu cev ili za sistem sa ograničenim brojem informacija o dinamičkom ponašanju potrošača. Validacije na celim sistemima su ograničene posebno u slučaju kada se posmatraju stvarni podaci od potrošača. Najčešće su dostupni podaci o vremenskoj promeni toplotnog zahteva potrošača za samo jedan deo potrošača. Za preostale potrošače vremenske promene su veštački formirane ili je primenjena srednja 13 godišnja vrednost potrošnje toplote, pri čemu nedostaje dinamička promena temperature kod potrošača za posmatrani vremenski period. Slika 2.1 Algoritam određivanja optimalne odlazne temperature tople vode iz toplane korišćenjem simulacija prelaznih procesa. Vremenska promena ponašanja potrošača iz grada Naestved u Danskoj je iskorišćena za istraživanje mogućnosti nekih modela da predvide prostiranje talasa temperature kroz čitavu mrežu daljinskog grejanja (Gabrielaitiene i dr., 2007). Sistem daljinskog grejanja u Naestved-u je modeliran primenom dva pristupa, metodom čvora razvijenom na Tehničkom univerzitetu Danske i primenom komercijalnog softvera TERMIS. Rezultati dobijeni primenom ovih modela su poređeni sa izmerenim vrednostima. Podsistem daljinskog grejanja u Naestved-u obuhvata 10%, odnosno 6 MW, ukupnog toplotnog opterećenja celog sistema daljinskog grejanja. Ovaj podsistem snabdeva 41 potrošača, čime su obuhvaćeni stambeni blokovi, administrativni i industrijski objekti. U administrativnim i industrijskim objektima 35-40% njihovog toplotnog zahteva otpada na pripremu sanitarne tople vode, dok je ovaj procenat u 14 stambenim zgradama manji i iznosi 25%. Na Slici 2.2 je prikazana distributivna mreža podsistema daljinskog grejanja u Naestved-u. Temperatura razvoda u podsistemu daljinskog grejanja je merena i korišćena je kao ulazni podatak za model. Tačka merenja ukupne količine toplote predate podsistemu je tretirana u modelu kao izvor toplote (Slika 2.2). Potrošači su na Slici 2.2 predstavljeni malim krugovima. Svaki krug predstavlja jednu podstanicu u okviru stambene zgrade. Slika 2.2 Podsistem daljinskog grejanja u Naestved-u, Danska. (Gabrielaitiene i dr., 2007) Merenja spoljašnje temperature, temperature razvoda i povrata na izlazu iz toplotnog izvora, temperature razvoda i povrata na primarnoj strani podstanica su izvršena u toku dva radna dana 21. i 22.01.2004. godine, za koje se ispostavilo da su bili najhladniji dani te godine sa spoljašnjom temperaturom od -9 0C. Vremenska promena potrošnje toplote je merena kod svih potrošača u okviru ovog podsistema. Podaci korišćeni u radu (Gabrielaitiene i dr., 2007) su mereni u vremenskim intervalima od sat 15 vremena, osim očitavanja na 5 min za šest najvećih potrošača, koji predstavljaju 50% ukupnog toplotnog opterećenja. Izmereni podaci su omogućili validaciju modela, kao i ispitivanje uticaja primene srednjih godišnjih vrednosti potrošnje toplote na rezultate modeliranja. Prilikom modeliranja strujanja u složenoj cevnoj mreži sistema daljinskog grejanja u Naestved-u posmatrano je idealizovano jednodimenzionalno strujanje fluida pri kome su promene temperatura uzrokovane brzinom strujanja i akumulacijom toplote u cevima. Uticaj turbulentnih fluktuacija i pojava sekundarnog strujanja su zanemareni zbog manjeg uticaja na dinamiku procesa. Metoda čvora i program TERMIS se zasnivaju na kvazi-dinamičkom pristupu, gde se temperatura procenjuje dinamički, dok se pritisak i protok računaju na osnovu stacionarnog modela. Jednačina (2.3) predstavlja osnovu procene dinamičke promene temperature   0v vp v o T T c C h T T t x         (2.3) gde je cp specifični toplotni kapacitet, C toplotni kapacitet, h koeficijent prelaza toplote, Tv temperatura vode, a To temperatura okoline. Razlika modela se ogleda u određivanju razmene toplote između fluida i okoline, odnosno u određivanju toplotnog kapaciteta materijala cevi i akumulacije toplote u njima, koja je zanemarena u TERMIS-u. Koeficijent prelaza toplote određen metodom čvora je primenjen i u TERMIS-u. Isti vremenski korak od 5 min, koji se poklapa sa rezolucijom merenja, je primenjen u oba programa kako bi se osigurali slični uslovi modeliranja. Diskretizacija cevovoda je urađena automatski i zavisi od vremenskog koraka. Vrednosti dobijene primenom modela su upoređene sa izmerenim podacima temperature povrata u toplotnom izvoru (Slike 2.3 i 2.4) i temperatura razvoda kod potrošača (Slike 2.5, 2.6 i 2.7). Niža temperatura povrata je poželjna kada postoji hladnjak dimnih gasova da bi se što bolje iskoristila otpadna toplota dimnih gasova. Da bi se zadovoljio toplotni zahtev potrošača moguće je povećati maseni protok nosioca toplote uz manju temperatursku razliku razvoda i povrata ili je moguće smanjiti potreban maseni protok, a time i snagu i troškove pumpanja, a povećati temperatursku razliku, snižavanjem temperature povrata. Sračunate temperature povrata toplotnog 16 izvora su prikazane na Slici 2.3 zajedno sa izmerenim vrednostima. Slaganje rezultata je zadovoljavajuće, ali se primećuje da su sve sračunate vrednosti niže od izmerenih. Slika 2.3 Izmerene i predviđene temperature povrata toplotnog izvora dobijene metodom čvora i primenom komercijalnog programa TERMIS. (Gabrielaitiene i dr., 2007) Slika 2.4 Izmerene i predviđene temperature povrata kada udeo potrošača ima konstantnu osrednjenu vrednost toplotnog opterećenja. (Gabrielaitiene i dr., 2007) 17 Izvršeno je i ispitivanje uticaja primene srednjih godišnjih vrednosti potrošnje toplote na rezultate modeliranja. Podaci o vremenskoj promeni potrošnje toplote su često poznati samo za najveće potrošače, dok su osrednjene vrednosti potrošnje toplote korišćene za preostale potrošače. Da bi se sagledalo kako ovo utiče na predviđanje temperature povrata, merene vrednosti vremenske promene potrošnje toplote korišćene su samo za 14 najvećih potrošača kojima odgovara 70% ukupnog toplotnog opterećenja i u drugom slučaju za samo 6 najvećih potrošača kojima odgovara 50% ukupnog toplotnog opterećenja, a za sve ostale potrošače primenjene su konstantne osrednjene vrednosti potrošnje toplote. Na osnovu podataka o godišnjoj potrošnji toplote je procenjeno da potrošači sa konstantnim toplotnim opterećenjem rade na 80% maksimalnog opterećenja. Takođe je pretpostavljeno da je razlika temperatura razvoda i povrata kod ovih potrošača konstantna i odgovara srednjoj vrednosti u toku godine. Slika 2.4 prikazuje promene izmerenih vrednosti temperatura povrata i temperatura sračunatih modelom sa vremenom u navedenim slučajevima, takođe je predstavljen i referentni slučaj sa osrednjenim vrednostima toplotnog opterećenja i razlike temperatura razvoda i povrata u svim podstanicama. Sračunata temperatura povrata u ovom slučaju odstupa do 5 0C od izmerenih vrednosti. Što je manji udeo potrošača sa vremenski promenljivom potrošnjom razmatran, veće je odstupanje sračunate od izmerene temperature povrata. Amplituda odstupanja izmerenih i predviđenih temperatura nije ravnomerna u toku simuliranog perioda (Slika 2.4). Najveća odstupanja se javljaju ujutro prvog i drugog dana (4 - 9h i 29 - 33h) i naveče drugog dana (39 - 45h), u periodima kada su najveća povećanja toplotnog zahteva. Validacija oba pristupa je izvedena i poređenjem rezultata dobijenih primenom modela sa izmerenim rezultatima temperature razvoda kod potrošača. Promene temperature razvoda i brzine sa vremenom su date na Slikama 2.5 - 2.7 za potrošače koji se nalaze na različitim udaljenostima od toplotnog izvora. Najbliže toplotnom izvoru je potrošač EE041 (Slika 2.5), koji se nalazi na rastojanju od izvora 0,1 km. Podstanica N1263 se nalazi na rastojanju 1,173 km od izvora, duž prilično ravnog dela cevovoda (Slika 2.7). Podstanica EE020 se nalazi na rastojanju 1,126 km u udaljenom delu cevovoda sa mnogo kolena i armature, odnosno gde su veći lokalni otpori. Na 18 Slikama 2.5-2.7 se može videti da je trend predviđenih i izmerenih temperaturskih profila sličan. Slika 2.5 Promena temperature razvoda i brzine sa vremenom za potrosača EE041 koji se nalazi blizu izvora toplote. (Gabrielaitiene i dr., 2007) Slika 2.6 Promena temperature razvoda i brzine sa vremenom za potrosača EE020 koji se nalazi na dalekom kraju cevovoda sa kolenima. (Gabrielaitiene i dr., 2007) 19 Poređenjem je utvrđeno da su lokalna odstupanja između rezultata dobijenih primenom modela i izmerenih temperatura razvoda značajnija za potrošače u udaljenim delovima mreže, a naročito tamo gde ima veći broj kolena i cevne armature, odnosno gde su veći lokalni otpori. Mogući razlog za duže prostiranje talasa je što talas putuje kroz veći broj kolena, lukova i armature u prethodnim cevovodima. Prenošenje toplote u takvim cevovodima je intenzivnije zbog razbijanja profila brzine i pojave sekundarnog strujanja. Brzine termalne difuzije su naglašenije a kao rezultat temperaturski talas se prostire brže. Ovi efekti nisu uzeti u obzir korišćenim modelima, pa je zbog toga odstupanje predviđenih i izmerenih temperatura veće. Slika 2.7 Promena temperature razvoda i brzine sa vremenom za potrosača N1263 koji se nalazi duž prilično ravnog cevovoda. (Gabrielaitiene i dr., 2007) 2.2 Hidraulički udar izazvan kondenzacijom pare Hidraulički udar izazvan kondenzacijom pare je pojava koja se dešava vrlo brzo i ima veliku razornu moć. Mnoga istraživanja se sprovode kako bi se bolje upoznala ova opasna pojava i odredili uslovi koji dovode do njenog nastanka. Preciznijim uvidom u mehanizam nastanka hidrauličkog udara mogu se poboljšati sistemi zaštite postrojenja i tako sprečiti njegove posledice. 20 Hidraulički udar izazvan kondenzacijom pare nastaje neposrednim kontaktom pothlađene tečnosti i pare, koji dovodi do intenzivne kondenzacije. Specifična zapremina pare je znatno veća od specifične zapremine tečnosti koja nastaje kondenzacijom ove pare, usled čega dolazi do pada pritiska u oblasti koju je zauzimala para. Razlika pritisaka u delovima cevovoda ispunjenim tečnošću i parom dovodi do kretanja stuba tečnosti ka parnom prostoru. Krećući se kroz cevovod stub tečnosti se ubrzava sve dok na svom putu ne naiđe na neku prepreku poput ventila, zatvorenog kraja cevi ili drugog stuba tečnosti. U blizini prepreke stub tečnosti usporava zatim udara o prepreku i odbija se od nje. Na mestu udara tečnosti dolazi do porasta pritiska, čiji se talas dalje prostire kroz masu tečnosti. U slučaju intenzivne kondenzacije pare masa tečnosti može dostići znatno ubrzanje pre udara, tako da porast pritiska pri udaru može biti razoran. Dakle, amplitude talasa pritiska indukovanog pojavom kondenzacije pare u direktnom kontaktu sa tečnom fazom mogu biti vrlo velike i prouzrokovati mehanička oštećenja opreme, kao što su oštećenja zida cevi, cevne armature, ovešenja ili sudova pod pritiskom, koja mogu biti takva da ugroze bezbednost zaposlenih u pogonu. U literaturi je navedeno da je do pojave hidrauličkog udara izazvanog intenzivnom kondenzacijom pare dolazilo u različitim termohidrauličkim sistemima: u sistemu napojne vode i generatoru pare nuklearne elektrane (Serkiz, 1983; Beuthe, 1997), u napojnoj pumpi parnog kotla u termoelektrani (de Vries i Simon, 1985), u sistemu daljinskog grejanja (Kirsner, 1999) i u rashladnom sistemu sa amonijakom (Martin, 2009). Do termohidrauličkog udara u sistemu napojne vode u termoelektrani dolazi nakon ispada napojne pumpe ili pri ponovnom uspostavlјanju punog protoka napojne vode nakon prestanka napajanja ili značajnog smanjenja u dotoku napojne vode (Stevanović i dr., 2008; Stevanović i dr., 2011). Naime, tokom prestanka napajanja kotla, smanjuje se protok kroz ekonomajzer i isparivač, pritisak u cevnom sistemu opada, a voda dodatno adijabatski isparava, tako da se znatno povećava udeo pare u cevima. Nakon ovog perioda prekida napajanja, pri ponovnom uspostavlјanju protoka, pothlađena napojna voda iz linije napojne vode dolazi u dodir sa parom što može 21 dovesti do intenzivne kondenzacije pare, naglog ubrzavanja stuba tečnosti i pojave termohidrauličkog udara. Takođe, nakon ispada napojnih pumpi, do mešanja pothlađene vode i dvofazne mešavine može doći i usled gravitacionog dejstva, kao što je, na primer, slučaj u usisnom cevovodu napojne pumpe, kojim se pothlađena voda dovodi do napojne pumpe iz napojnog rezervoara, koji se nalazi na višoj koti u odnosu na pumpu (Slika 2.8). Slika 2.8 Ključanje zagrejane vode u napojnoj pumpi i usisnom cevovodu pumpe nakon ispada bloka. Razlog pojave ovog termohidrauličkog udara je niz prelaznih procesa u turbopostrojenju. Zatvaranje stop ventila ispred turbina visokog i srednjeg pritiska nakon ispada bloka dovodi do zaustavljanja protoka pare kroz stupnjeve turbine i kroz linije oduzimanja kojima se greje linija kondenzata i voda u napojnom rezervoaru. Protok kondenzata iz kondenzatora turbine niskog pritiska ka napojnom rezervoaru se i dalje održava, ali temperatura kondenzata opada, s obzirom da je prestalo njegovo zagrevanje u zagrejačima kondenzata parom iz turbinskih oduzimanja. Uticanje hladnijeg kondenzata dovodi do opadanja pritiska u napojnom rezervoaru. Napojni rezervoar je ispunjen zasićenom vodom i parom, tako da sniženje temperature napojne vode dovodi do kondenzacije pare i opadanja pritiska. Ovo opadanje pritiska dovodi do ključanja zagrejane vode koja ispunjava usisni cevovod napojne pumpe i samu pumpu, Para iz turbine - ZAUSTAVLJENA Kondenzat iz ZVP - ZAUSTAVLJEN Rezervoar napojne vode Glavna napojna pumpa - ISKLJUČENA Buster pumpa - ISKLJUČENA ka kotlu - ZAUSTAVLJENO Kondenzat iz ZNP 22 kao što je prikazano na Slici 2.8. Protok kroz glavnu napojnu pumpu je zaustavljen, a nepovratni ventil na potisu napojne pumpe je zatvoren da ne bi došlo do povratnog strujanja iz cevnog sistema kotla na visokom pritisku u napojni rezervoar na nižem pritisku. Nakon desetak minuta od ispada pumpe, ohlađena voda iz napojnog rezervoara dolazi u dodir sa parom u usisnom cevovodu i napojnoj pumpi, para se brzo kondenzuje, opada pritisak u turbonapojnoj pumpi i usisnom cevovodu pumpe, tako da se vodeni stub iz napojnog rezervoara ubrzava ka pumpi i dolazi do termohidrauličkog udara. Termohidraulički udar u nuklearnim elektranama nastaje zbog pojave parnog prostora u cevima koje su u normalnim uslovima rada ispunjene vodom, zbog kondenzacije pare u cevovodima koji pri normalnim uslovima rada sadrže vodu i vodenu paru ili, kao što je najčešći slučaj, prilikom brzog otvaranja ili zatvaranja ventila. Prema podacima iz literature Liu i dr. (1997) preko 100 hidrauličkih udara izazvanih kondenzacijom pare je zabeleženo u nuklearnim elektranama u SAD u periodu od 1969 – 1981. Stoga su sprovedena eksperimentalna istraživanja hidrauličkog udara čiji je osnovni cilj određivanje graničnih uslova, odnosno praga, kada dolazi do ove pojave, kao i ilustracija uticaja hidrauličkih i termičkih efekata na ovu nestacionarnu pojavu. Osnovni parametar od kog zavisi predviđanje ponašanja sistema pri hidrauličkom udaru je brzina kondenzacije. Na brzinu kondenzacije utiče koeficijent prelaza toplote pri kondenzaciji, čija vrednost prema eksperimentalnim istraživanjima (Liu i dr., 1997) može varirati i do nekoliko redova veličine, npr. od 1 kW/(m2K) u raslojenom toku do oko 2000 kW/(m 2 K) u mehurastom toku. Liu i dr. (1997) su na osnovu svojih istraživanja došli do zaključka da je intenzitet indukovanog porasta pritiska funkcija pritiska vode u rezervoaru eksperimentalne aparature, temperature vode i temperature pare, kao i da postoj prag u uslovima koji dovode do pojave hidrauličkog udara izazvanog kondenzacijom pare. Primetili su takođe, da postoji značajno rasipanje eksperimentalnih rezultata za iste početne uslove eksperimenta. Do termohidrauličkog udara u sistemu daljinskog grejanja došlo je u Fort Wainwright-u, Alaska (Kirsner, 1999) tokom remonta sistema, odnosno uklanjanja azbestne izolacije sa glavnog parovoda. U normalnim uslovima temperatura, u podzemnim tunelima u kojima su parovodi smešteni, je 70 0C, što je previsoka 23 temperatura za boravak radnika. Pošto nije bilo moguće prekinuti strujanje pare kroz parovode dok se ne završi remont odlučeno je da se temperatura pare snizi u toku noći tako da u jutarnjim časovima temperatura u tunelima bude podnošljivih 50 0C. U toku noći se u parovodima nakupljao kondenzat. Uvođenjem pare pod visokim pritiskom u liniju glavnog parovoda, sa koje je uklonjena azbestna izolacija, svi uslovi za hidraulički udar izazvan kondenzacijom pare su ostvareni. Kada je jedan od radnika otvorio ventil, toliko da samo podigne sedište ventila, čuo se prasak i para je eksplodirala. Vrela para i kondenzat su ispunili tunel sa toplotnim talasom temperature 93 0 C, koji je uzrokovao teške povrede radnika i otežavao njihovo kretanje kroz postrojenje u pokušaju spasavanja jer je sve bilo klizavo i vidljivost otežana. Hidraulički udar izazvan kondenzacijom pare takođe može da se javi i u rashladnim sistemima sa amonijakom (Martin, 2009). Ovi sistemi rade na znatno nižim pritiscima nego što su radni pritisci u napojnim linijama u termo ili nuklearnim elektranama pa se ne očekuju visoki pikovi pritiska od 5 MPa i teške povrede radnika. Nasuprot tome, često se hidraulički udar izazvan kondenzacijom pare javlja u povratnim linijama rashladnih sistema, kada se topao gas uvodi preko pothlađene tečnosti kao deo periodičnog ciklusa odmrzavanja, i uzrokuje teška oštećenja imovine i stradanje ljudi iz pogona. Istraživanja hidrauličkog udara izazvanog intenzivnom kondenzacijom pare u suprotnosmernom strujanju pothlađene vode i pare u horizontalnim cevima i cevima sa malim nagibom su zbog svog velikog značaja za sigurnost rada generatora pare i sistema napojne vode u nuklearnim elektranama vršena i ranije (Bjorge i Griffith, 1984; Chun i Yu, 2000). Date su preporuke za prevenciju od hidrauličkog udara izazvanog kondenzacijom pare u takvim (sličnim) uslovima dvofaznih strujanja. Takođe, hidraulički udar izazvan uticanjem hladne vode u horizontalnu cev ispunjenu parom je eksperimentalno istraživan i numerički simuliran (Barna i dr., 2010). Numeričke simulacije su izvršene primenom programa WAHA3, koji se zasniva na jednodimenzionalnom modelu dva fluida za nestacionarno strujanje pare i vode. U ovoj studiji hidraulički udar izazvan intenzivnom kondenzacijom pare je uzrokovan brzom kondenzacijom parnog mehura koji je zarobljen između čepova pothlađene vode, koji se 24 formiraju zbog Kelvin-Helmholtz-eve nestabilnosti sloja vode u raslojenom dvofaznom toku u horizontalnoj cevi. Još jedan mehanizam koji dovodi do hidrauličkog udara izazvanog kondenzacijom pare je pojava poznata pod nazivom vodeni top. Može se javiti prilikom ispuštanja pare u zapreminu hladne vode. Para može ostati zarobljena u cevi kada se ventil za ispuštanje pare zatvori. Kontakt pare i hladne vode uzrokuje brzu kondenzaciju i pothlađena voda biva uvučena u cev. Formirani vodeni čep se zaustavlja na ventilu uz veliki porast pritiska. Simulacija vodenog topa primenom termohidrauličkog programa za predviđanje sigurnosti rada nuklearnog reaktora RELAP5/MOD3 je prikazana u (Yeung i dr., 1993). Hidraulički udar izazvan kondenzacijom pare je pojava koja se javlja kada se u direktnom kontaktu nađu para i pothlađena voda, stoga je za njenu numeričku simulaciju primenom kompjuterskih programa neophodno razviti model direktne kondenzacije. Jedan od primera gde se pothlađena voda i para u stanju zasićenja nađu u kontaktu preko velike površine pri čemu dolazi do trenutne kondenzacije je hlađenje jezgra reaktora u udesnim situacijama, gde se para hladi ubrizgavanjem hladne vode. Do sada ne postoji bezbedan i potvrđen metod koji dovoljno dobro predviđa termohidrauličko opterećenje za vreme hidrauličkog udara izazvanog kondenzacijom pare. Na Tehničkom univerzitetu u Hamburgu je, u okviru istraživačkog saveza CIWA (Condensation Induced Water Hammer), koji je osnovan 2010. godine od strane Nemačkog federalnog ministarstva obrazovanja i istraživanja (BMBF), sa ciljem pronalaska efikasne metode za predviđanje hidrauličkog udara izazvanog kondenzacijom pare, napravljena eksperimentalna aparatura kojom je moguće vršiti istraživanje parametara koji određuju hidraulički udar izazvan kondenzacijom pare. Eksperimentalna aparatura sastoji se od test dela sa 5 merača pritiska koji velikom brzinom snimaju promenu pritiska (Slika 2.9). Generator pare proizvodi paru u stanju zasićenja na pritisku od 0,4 do 1 MPa, koja ulazi u test deo sa leve strane kroz separator 1. Zapreminski protok pare kontroliše se regulacionim ventilom sa gornje desne strane eksperimentalne aparature. Voda ulazi u test deo sa desne strane ili kroz ulaznu mlaznicu ili kroz separator 2, formirajući tako suprotnosmerni režim strujanja. 25 Pothlađena voda struji u zatvorenom krugu kroz razmenjivač toplote i pumpu. Svi podaci se beleže mernim akvizicionim sistemima. Slika 2.9 Eksperimentalna aparatura na Hamburško tehnološkom univerzitetu. (Urban i Schlüter, 2013) Za rad (Urban i Schlüter, 2013) je sprovedeno 36 eksperimenata hidrauličkog udara. U predstavljenom eksperimentu početni pritisak u sistemu je 0,4 MPa, a zapreminski protok pare u stanju saturacije je manji od 0,2 m 3 /h. Eksperiment počinje ubrizgavanjem pothlađene vode zapreminskim protokom 1,4 m3/h i temperaturske razlike 27 K u odnosu na paru u stanju saturacije. 30 s posle početka eksperimenta je detektovan i snimljen hidraulički udar. Na Slici 2.10 je prikazana promena pritiska u 5 merača pritiska u test delu, u trenutku udara. Merač pritiska p1 se nalazi blizu mesta ulaza pare, a p5 je sa desne strane blizu mesta ubrizgavanja vode. Samo su 4 merača pritiska detektovala porast pritiska usled hidrauličkog udara. Razlog tome je raspodela faza unutar test dela. Merač pritiska p1 nije bio u kontaktu sa tečnom fazom u trenutku udara. Na osnovu rezultata predstavljenih na Slici 2.10 može se zaključiti da se kolaps velikog parnog mehura dogodio u delu levo od p2. Promena pritiska tokom kondenzacije velikog parnog mehura se može posmatrati na krivim linijama p2-p5. 26 Na Slici 2.10 svaki merač pritiska je snimio zaravan maksimalnog pritiska od oko 70 bar, trajanja 200μst  . Dužina ubrzanog vodenog čepa može da se odredi na osnovu vremenske razlike između rastućeg i opadajućeg nagiba krive linije. Opadajući nagib je rezultat odbijanja talasa visokog pritiska. Rezultujući talas podpritiska se prostire kroz tečnu fazu istom brzinom kao i talas visokog pritiska, nastao nakon sudara dva vodena čepa, prosečnom brzinom 1429eksa  m/s, što je u okviru teorijskih granica (brzina zvuka u vodi na pritisku od 0,4 MPa i na temperaturi 117 0 C je 1523teoa  m/s), pri čemu kasni za vremensku razliku t . Ako se uzme konstantna brzina od 1429 m/s tada se dužina vodenog čepa može izračunati kao 2 eksa tL   . (2.4) Slika 2.10 Promena pritiska u meračima pritiska od p1-p5 u trenutku nastanka udara. (Urban i Schlüter, 2013) Uzimajući u obzir izračunatu dužinu vodenog čepa, lokacija hidrauličkog udara može biti sužena na prostor dužine 343 mm levo od merača pritiska p2 za posmatrani eksperiment. Eksperimenti pokazuju da se većina hidrauličkih udara javlja između merača pritiska p1 i p2. Nakon lociranja oblasti u kojoj se dešava hidraulički udar, u ovom delu se može ugraditi providna cev tako da je moguće videti stanje strujanja u cevi i sprovođenje optičkih istraživanja. Optičko istraživanje se može sprovesti 27 primenom kamere sa velikom brzinom snimanja do 5000 fps, pozadinskog osvetljenja i laserskog optičkog merenja temperatura i brzinskog polja tečne faze. U (Urban i Schlüter, 2013) su izvršeni testovi kako bi se potvrdila mogućnost ugradnje providne cevi i primene laserski indukovane fluorescencije (LIF) na određivanje raspodele temperatura. Posmatran je rezervoar vode, u kom se razdelna površina zagreva primenom toplog vazduha, čime nastaje slojevita raspodela temperature u rezervoaru. Na levoj strani na Slici 2.11 je prikazana kvantitativna raspodela temperature u rezervoaru vode dobijena primenom LIF, a sa desne strane je prikazana ista raspodela sa temperaturama merenim termoparovima. Temperaturski profil se dobro slaže sa predviđenom raspodelom temperatura. Referentna merenja temperature termoparovima, koja su na rastojanju 13 i 38 mm od granice, dobro se slažu sa temperaturama izmerenim primenom optičke metode na ovim rastojanjima (Urban i Schlüter, 2013). Slika 2.11 Temperaturski profil u zapremini vode u rezervoaru. (Urban i Schlüter, 2013) Kao deo istraživačkog saveza CIWA TÜV NORD SysTec GmbH & Co. KG je razvio program za simulaciju brzih nestacionarnih promena pritiska u jednofaznom i dvofaznom strujanju u cevnim sistemima DYVRO. Blömeling i dr. (2013) su simulirali hidraulički udar izazvan kondenzacijom pare sa 1D programom DYVRO. Sproveli su proračune sa konstantnom vrednošću koeficijenta prelaza toplote i izvršena je verifikacija programa DYVRO poređenjem dobijenih rezultata sa rezultatima proračuna primenom dvofaznog čepastog modela. 28 Posmatrani slučaj, na kome je izvršeno poređenje modela, se sastoji od cevi, postavljene pod uglom  u odnosu na vertikalu, povezane sa rezervoarom u kome vlada konstantan pritisak p0 i temperatura T0 (Slika 2.12). Cev je delimično ispunjena i vodom i parom. Voda koja se nalazi između rezervoara i pare u domenu obeleženim sa slug predstavlja čep. Domen koji sadrži mešavinu vode i pare je obeležen sa .mix Cev može da se puni vodom iz pravca nasuprot rezervoaru. Brzina i temperatura napojne vode obeležene su na Slici 2.12 sa FWu i FWt . Oblik razdelne površine između vode i pare je stalan i jednak dvostrukoj vrednosti poprečnog preseka cevi, osim u slučaju kada nema napojne vode nasuprot čepa. Tada je razdelna površina jednaka poprečnom preseku cevi. Slika 2.12 Konfiguracija čepastog modela. (Blömeling i dr., 2013) DYVRO i čepasti model su primenjeni na isti scenario hidrauličkog udara sa jednakim početnim i graničnim uslovima uzetim iz literature (Blömeling i dr., 2013):  prečnik cevi D = 4,3 cm,  pritisak u rezervoaru p0 = 21,7 bar,  temperatura vode T0 = TFW = 27 0 C,  dužina čepa 4D,  dužina zapremine sa zasićenom parom 2D,  brzina napojne vode uFW = 0,9125 m/s  pritisak pare pm = p0 i  koeficijent prelaza toplote  2200 kW/ m K  . 29 Odlično slaganje oba modela (Slika 2.13) je na početku simulacije, što znači da su oba modela izračunala istu vrednost za brzinu kondenzacije. Kasnije dolazi do odstupanja rezultata. Ovo odstupanje proističe iz činjenice da DYVRO mora da odredi razdelnu površinu između faza. Za razliku od čepastog modela, DYVRO razvuče razdelnu površinu preko nekoliko kontrolnih zapremina zbog numeričke difuzije, tako da se razdelna površina s vremenom pomalo povećava, što dovodi do većih brzina kondenzacije i odstupanja na Slici 2.13. Slika 2.13 Promena pritiska pare. (Blömeling i dr., 2013) Blömeling i dr. (2013) su takođe izvršli i validaciju DYVRO-ovog modela i čepastog modela sa eksperimentalnim podacima za vodeni top. Eksperimentalna aparatura vodenog topa prikazana na Slici 2.14 se sastoji od vertikalne cevi dužine 0,7 m i prečnika 0,038 m, koja je u početnom trenutku ispunjena vodom i čiji je otvorеni kraj uronjen u veliki rezervoar sa pothlađenom vodom na pritisku 1 bar i temperaturi 18,3 0C. Vodena para u stanju zasićenja (1 bar, 99,6 0C) se uvodi u cev sa gornje strane vertikalne cevi čime se zagreva voda u cevi i zidovi cevi i potiskuje voda naniže ka rezervoaru. Kada temperatura vode na razdelnoj površini pare i vode unutar cevi dostigne vrednost temperature zasićenja, para počinje da zauzima mesto vode u cevi. Kada para dođe do donjeg kraja cevi, ona dolazi u kontakt sa hladnijom vodom i vrlo brzo se kondenzuje. Zbog razlike pritisaka voda se uvlači u cev i udara o njen gornji zatvoreni kraj, pri čemu dolazi do kondenzacije pare i nastanka talasa pritiska. 30 Slika 2.14 Eksperimentalna aparatura za vodeni top. (Blömeling i dr., 2013) Eksperimentalni i numerički rezultati, su dati na Slici 2.15. Ceo proces se može podeliti na tri faze:  kondenzacija pare na vodenoj površini. Pritisak u zapremini pare opadne 0,3 bara i uzrokuje ubrzanje vodenog stuba naviše u cev,  u drugoj fazi pritisak osciluje oko vrednosti od 0,7 bara što znači da kretanje vodenog čepa u većoj ili manjoj meri kompenzuje kondenzaciju pare i  udar vodenog stuba o gornji kraj cevi sa maksimalnom vrednošću pika pritiska pri udaru od 76 bara za 135 ms. Ovo je stohastički eksperiment, u 80 ponavljanja eksperimenta sa istim početnim uslovima dolazi do pikova pritiska između 34 - 90 bara. Kada se uzme u obzir stohastička priroda eksperimenta slaganje rezultata numeričkih simulacija i eksperimenta je dobro. 31 Slika 2.15 Pritisak pare - eksperimentalne vrednosti i rezultati simulacija. (Blömeling i dr., 2013) Određivanje koeficijenta prelaza toplote je veoma složen zadatak jer on zavisi od velikog broja faktora pa mora da bude određen za svaki oblik strujanja i za svaku geometriju posebno. Svi prethodno sprovedeni proračuni su imali konstantnu vrednost koeficijenta prelaza toplote. Blömeling i dr. (2013) su razvili novi pristup za određivanje koeficijenta prelaza toplote za određena strujanja i geometrije, koji se zasniva na teoriji obnavljanja površine, pri čemu nisu potrebne intervencije od strane korisnika. Teorija obnavljanja površine zasniva se na konceptu turbulentnih vrtloga koji mogu narušiti laminarni podsloj tečnosti na razdelnoj površini u kontaktu sa parom i zameniti tečnost u stanju saturacije, koja je u kontaktu sa parom, sa svežom 32 pothlađenom tečnošću iz dubljih zona. Na ovaj način kondenzacija je kontrolisana vremenom kontakta na razdelnoj površini, koje se još naziva period obnavljanja površine i predstavlja kontakt vrtloga sa tečne strane i parne faze u stanju saturacije, a računa se kao odnos karakteristične dužine Lt i brzine vrtloga Vt. Ako su karakteristične vrednosti poznate koeficijent prelaza toplote se može izračunati kao 0.5 1Pr n n t t t t t V L            (2.5) gde su  i n konstante, t , t i Prt su toplotna provodljivost, kinematska viskoznost i Prandtl-ov broj tečne faze, sledstveno. Izbor konstanti i karakterističnih turbulentnih veličina zavisi od odabranog modela turbulencije. Pouzdan način predviđanja direktne kondenzacije je neophodan zato što ova pojava pod određenim okolonostima može dovesti do hidrauličkog udara izazvanog kondenzacijom pare i time dovesti do oštećenja komponenata u cevnom sistemu. Direktna kondenzacija je opasnija pri nižim pritiscima u sistemu zbog većih razlika u gustinama između dve faze. Ceuca i dr. (2013) su poredili različite pristupe za izračunavanje koeficijenta prelaza toplote zasnovane na teoriji obnavljanja površine. Razvijen je hibridni model za koeficijent prelaza toplote na osnovu dva modela teorije obnavljanja površine. Osnovna razlika između dva modela za određivanje koeficijenta prelaza toplote, zasnovana na teoriji obnavljanja površine, je u načinu na koji određuju kako je karakteristična dužina vrtloga odgovorna za odvođenje toplote sa razdelne površine na zapreminu tečnosti. Model Hughes i Duffey (Ceuca i dr., 2013) razmatra najmanje veličine u turbulentnom strujanju, Kolmogorov-e mikro veličine, i smatra da su one odgovorne za odvođenje toplote što rezultira kraćim periodom obnavljanja i dozvoljavanjem da više toplote bude preneseno. Shen-ov model (Ceuca i dr., 2013) smatra da veći vrtlozi transportuju više energije sa razdelne površine i tako produžuju period obnavljanja površine. U (Ceuca i dr., 2013) je dat prikaz hibridnog modela za proračun koeficijenta prelaza toplote zasnovanog na teoriji obnavljanja površine. Hibridni model obuhvata i velike i male vrtloge koji odvode toplotu sa razdelne površine i predaju je zapremini tečnosti (Slika 2.16). Model se zasniva na pretpostavci da se energija prenosi pretežno 33 malim vrtlozima kada turbulentni Reynolds-ov broj prelazi prag od 500, dok se ispod ove vrednosti za izračunavanje koeficijenta prelaza toplote koriste veliki vrtlozi. Slika 2.16 Shematski prikaz čepastog toka sa rasporedom vrtloga odgovornih za prenos toplote sa razdelne površine na zapreminu tečnosti. (Ceuca i dr., 2013) Eksperimentalna merenja su izvršena na postrojenju PMK 2 na Istraživačkom institutu za atomsku energiju u Budimpešti, Mađarska. Ovo eksperimentalno postrojenje ima sličnu konstrukciju kao glavna grana hladioca jezgra reaktora tipa VVER. Dve od tri numeričke simulacije sprovedene primenom CFD-a se dobro slažu sa merenjima. Na postrojenju je izvršeno 35 simulacija hidrauličkog udara. Shematski prikaz postrojenja je dat na Slici 2.17. Voda temperature 25 0 C, masenog protoka 1,01 kg/s dovodi se kroz cev prečnika 73 mm u test deo koji sadrži paru u stanju saturacije na pritisku od 1,45 MPa. Slika 2.17 Shematski prikaz eksperimentalne aparature PMK 2. Lokacija sonde koja sadrži termoparove i senzore za zapreminski udeo pare je označena sa 8, merači pritiska sa 10, dok je senzor od žičane mreže broj 9. (Ceuca i dr., 2013) 34 Dobijeni oblici strujanja pri CFD simulaciji primenom hibridnog modela su predstavljeni na Slici 2.18. Voda ulazi u horizontalni test deo i struji sa leve na desnu stranu. Kada voda dođe do desnog kraja, formira se talas male amplitude. Talas se ubrzava pri svom kretanju nalevo zbog opadanja pritiska u levom delu cevi što je rezultat direktne kondenzacije pare i na taj način se razvija čep. Kada vodeni talas dostigne visinu cevi, voda zarobljava parni mehur, koji se kondenzuje dozvoljavajući vodi da zauzme njegovo mesto i uzrokuje ubrzanje vodenog talasa ka levom kolenu, istovremeno menjajući oblik strujanja od horizontalnog raslojenog toka u čepasti tok. Strujanje se nastavlja sa leve strane ka desnoj, sve dok voda konačno ne ispuni test deo i izađe iz proračunskog domena. Ubrzanje vode zbog kolapsa mehura je pokretač hidrauličkog udara izazvanog kondenzacijom pare. Slika 2.18 Promena zapreminskog udela parne faze i oblici strujanja dobijeni prilikom simulacije PMK 2 eksperimenta primenom CFD-a i hibridnog modela u karakterističnim vremenskim trenucima 6,8 s, 7 s, 7,2 s, 7,25 s, 7,3 s, 7,6 s, 8,4 s i 8,5 s. (Ceuca i dr., 2013) Validacija rezultata je urađena poređenjem rezultata CFD-a i izmerenih rezultata na eksperimentalnom postrojenju. Dostupni rezultati obuhvataju 4 lokalno merene temperature i 4 lokalno merena zapreminska udela pare, obeležena brojem 8 na Slici 2.17, nazvana VT1 do VT4 u smeru strujanja. Posmatranjem rezultata CFD simulacije na Slici 2.19, primenom samo modela koji uzima u obzir uticaj velikih vrtloga (Shen-ov model), može se zaključiti da se 35 temperature razlikuju od izmerenih vrednosti, iako je dobijeni oblik strujanja prilično realan. Razvijeni hibridni model postiže dobro slaganje sa eksperimentalnim merenjima. Promene temperatura imaju sličan karakter kao promene merenih rezultata. Niža temperatura vode na kraju CFD simulacije je uzrokovana graničnim uslovom - adijabatski izolovanim zidom i činjenicom da je čitava cev na početku simulacije bila ispunjena parom u stanju saturacije što nije slučaj u eksperimentu. Slika 2.19 Poređenje promena temperature u toku vremena na mestu postavljanja merne sonde primenom CFD modela sa Shen-ovim i hibridnim modelom i sa izmerenim rezultatima. (Ceuca i dr., 2013) Kod kotlova, koji rade sa visokim radnim parametrima, posebna pažnja se mora posvetiti termohidrauličkoj sigurnosti, pre svega zaštiti od pojave termohidrauličkog udara u sistemu napojne vode, ekonomajzeru, spusnim cevima od ekonomajzera i isparivačkim cevima kotla. Protočni kotlovi su podložni brzim promenama termohidrauličkih parametara tokom prelaznih procesa, kako zbog visokih toplotnih opterećenja, tako i zbog relativno manje količine vode u cevnom sistemu. Zbog manje količine vode u sistemu, u slučaju poremećaja u protoku napojne vode kod protočnih 36 kotlova, u kraćem periodu može doći do poremećaja u hlađenju toplotno visokoopterećenih isparivačkih ili pregejačkih površina. Prestanak efikasnog hlađenja ovih površina dovodi do naglog porasta temperature metala zida cevi, što uzrokuje promene u strukturi metala, dodatna termomehanička naprezanja, a može dovesti i do pucanja cevi. Ekstremni slučajevi termičkih naprezanja nastaju kada u cevnom sistemu isparivača nema protoka, a u ložištu je prisutna vatra (usled prisustva zaostalog uglјa koji sagoreva ili rada mazutnog gorionika). Korišćenjem kompjuterskog programa za simulacije i analize prelaznih termohidrauličkih procesa u cevovodima i sudovima pod pritiskom, razvijenog na Mašinskom fakultetu u Beogradu, izvršena je simulacija pojave termohidrauličkog udara u sistemu napojne pumpe Termoelektrane "Nikola Tesla B" (TENT B) (Stevanović i dr., 2008, Stevanović i dr., 2011). Slika 2.20 Promena pritiska u napojnom rezervoaru pri ispadu bloka TENT B dana 07.11.2006. god. u u 16:15 časova (sa delovanjem sigurnosno-zaštitnog sistema). Para koja nastaje u usisnom cevovodu i napojnoj pumpi pri ispadu bloka i opadanju pritiska u napojnom rezervoaru (Slika 2.20), struji naviše ka napojnom rezervoaru usled dejstva sile potiska u dvofaznoj mešavini vode i pare (Slika 2.8). Para pri suprotnosmernom strujanju ka napojnom rezervoaru delimično sa sobom povlači i vodu dejstvom međufaznog trenja na razdelnim površinama pare i vode, tako da dolazi do intenzivnog izbacivanja dvofazne mešavine iz usisnog cevovoda napojne pumpe u 37 napojni rezervoar. Ovo izbacivanje dvofazne mešavine iz usisnog cevovoda napojne pumpe uzrokuje usisavanje hladnije vode iz napojnog rezervoara u usisni cevovod. Za to vreme horizontalna linija usisnog cevovoda i napojna pumpa su većim delom ispunjeni parom. Pri ovim uslovima ohlađena napojna voda koja utiče iz napojnog rezervoara u vertikalni usisni cevovod turbonapojne pumpe dolazi u dodir sa parom i nastaju uslovi za intenzivnu kondenzaciju pare. Para se kondenzuje na pothlađenoj vodi, prostor koji je zauzimala para preuzima voda, tako da se stub vode ubrzava ka napojnoj pumpi i zatvorenom nepovratnom ventilu na potisu pumpe, udara o kolo pumpe i zatvoreni kraj cevovoda i dolazi do intenzivnog hidrauličkog udara. Ovoj pojavi pogoduje i velika visinska razlika između napojnog rezervoara i napojne pumpe od 50 m. Na Slici 2.21 trenutak udara je prikazan impulsnim porastom pritiska preko 200 bara na mestu napojne pumpe u 550 s. Ova promena pritiska razarajuće deluje na strukturu cevovoda i napojne pumpe. Slika 2.21 Promena pritiska u napojnoj pumpi nakon ispada sa pojavom termohidrauličkog udara. Metode za zaštitu od termohidrauličkog udara u sistemu napojne pumpe su prikazane na Slikama 2.22 - 2.24. Metoda predložena za sprečavanje termohidrauličkog udara, prikazana na Slici 2.22, sastoji se od odzračne linije čiji je zadatak da paru koja bi nastala nakon ispada pumpe sprovede u napojni rezervoar. Ova metoda nije efikasna jer se ne formiraju dovoljne uzgonske sile za uspostavljanje strujanja pare sa potisa napojne pumpe ka napojnom rezervoaru. Drugo rešenje, prikazano na Slici 2.23, podrazumeva ugradnju cirkulacione pumpe koja bi se uključila nakon ispada napojne 38 pumpe, a njenim radom bi se obezbedila cirkulacija napojne vode iz napojnog rezervoara i hlađenje pumpe. Ovo rešenje nije praktično i iziskuje znatne troškove za ugradnju cirkulacione pumpe i njen pogon. Treće rešenje, prikazano na Slici 2.24, predloženo je i uspešno primenjeno početkom devedesetih godina na TENT B (Jovanović i Jocić, 1996). Primenom ovog rešenja obezbeđuje se stalan protok kroz turbonapojnu pumpu putem otvaranja rasteretne linije ka atmosferskom ekspanderu, što dovodi do hlađenja pumpe i sprečava ključanje vode u sistemu napojne pumpe. Slika 2.22 Sprečavanje termohidrauličkog udara u sistemu napojne pumpe pomoću rasteretne linije za isticanje pare sa potisa pumpe ka napojnom rezervoaru. Efekti mera zaštite od termohidrauličkog udara nakon ispada turbonapojne pumpe, koji se postižu hlađenjem turbonapojne pumpe putem ispuštanja vode is sistema turbonapojne pumpe ka atmosferskom ekspanderu, kroz odzračne linije na TENT-u B, prikazani su na Slici 2.25. Na Slici 2.25 je prikazana promena stepena suvoće u turbonapojnoj pumpi u slučaju bez delovanja sistema zaštite, kada dolazi do termohidrauličkog udara, kao i u slučaju kada deluje sistem zaštite. U slučaju bez delovanja zaštite primećuje se porast stepena suvoće u turbonapojnoj pumpi iznad nule, do vrednosti kada nastaje hidraulički udar (porast pritiska na Slici 2.21). Vidi se da u slučaju delovanja sistema zaštite ne dolazi do ključanja vode u turbonapojnoj pumpi, pa samim tim ni do stvaranja uslova za pojavu termohidrauličkog udara. Para iz turbine Kondenzat iz ZNP Kondenzat iz ZVP Rezervoar napojne vode Motorni ventil Glavna napojna pumpa Buster pumpa ka kotlu 39 Slika 2.23 Različite metode za sprečavanje termohidrauličkog udara u sistemu napojne pumpe: pomoćna cirkulaciona pumpa omogućava hlađenje usisnog cevovoda i napojne pumpe. Slika 2.24 Različite metode za sprečavanje termohidrauličkog udara u sistemu napojne pumpe: princip rešenja primenjenog na TENT B sa isticanjem vode ka atmosferskom ekspanderu. Pomoćna cirkulaciona pumpa Para iz turbine Kondenzat iz ZNP Kondenzat iz ZVP Rezervoar napojne vode Glavna napojna pumpa Buster pumpa ka kotlu Atmosferski ekspander Para iz turbine Kondenzat iz ZNP Kondenzat iz ZVP Rezervoar napojne vode Motorni ventil Glavna napojna pumpa Buster pumpa ka kotlu 40 Slika 2.25 Promena stepena suvoće u turbonapojnoj pumpi sa i bez delovanja zaštite od termohidrauličkog udara. 41 NUMERIČKE METODE ZA PRAĆENJE PROSTIRANJA FRONTA 3. TALASA U ovom poglavlju su predstavljene osnovne numeričke metode za rešavanje transportnih jednačina kojima se opisuje kretanje diskontinuiteta talasa temperature. Više različitih numeričkih metoda je moguće primeniti za rešavanje istog problema. Razlike u rešenjima dobijenim primenom različitih modela su često veoma male, čime se otežava postupak izbora optimalne tehnike rešavanja. U odabiru metode pomaže stečeno iskustvo u programiranju primenom različitih metoda. Uobičajeni problemi koji se javljaju u mehanici fluida su veoma nelinearni. Sistem nelinearnih parcijalnih diferencijalnih jednačina koji opisuje neki problem mora da bude rešen za nepoznate pritiske, gustine, temperature i brzine. Rešavanje ovakvog sistema je komplikovano. Burgers je 1948. godine (Tannehill i dr., 1997) uveo jednostavnu nelinearnu jednačinu koja ima sve članove koji blisko preslikavaju fizičke osobine jednačina fluida, odnosno ova jednačina ima nestacionarni, konvektivni i viskozni član. Burgers-ova jednačina ima sledeći oblik 2 2 nestacionarni član konvektivni član viskozni član u u u u t x x          . (3.1) Ovo je jednačina paraboličkog tipa kada se uključi viskozni član, a ako se on zanemari dobija se jednačina hiperboličkog tipa 0 u u u t x       . (3.2) Zbog jednostavnosti posmatramo skalarnu jednačinu gde su u opštem slučaju obe nepoznate u i F(u) vektorske veličine 0. u F t x       (3.3) Ako se uvede Jacobian-ova matrica A=A(u), koja je jednaka /i jF u  u opštem slučaju i dF/du za pojednostavljenu jednačinu, jednačinu (3.3) možemo napisati i kao 0 u u A t x       . (3.4) 42 Ova jednačina je hiperboličkog tipa, što znači da su sve sopstvene vrednosti matrice A realni brojevi. Slede neke najčešće korišćene metode za rešavanje Burgers-ove jednačine sa tipičnim rezultatima i diskusijom uticaja nelinearnih članova (Tannehill i dr., 1997). 3.1 Lax-ova metoda Lax-ova metoda (1954) je metoda prvog reda za rešavanje parcijalnih diferencijalnih jednačina hiperboličkog tipa. Metode prvog reda tačnosti se retko koriste. Ova metoda je tipičan primer za pokazivanje primene na nelinearne jednačine i disipativnog karaktera rezultata. Osnovna parcijalna diferencijalna jednačina u formi zakona održanja (3.3) se za Lax-ovu metodu razvija u Taylor-ov red oko tačke  ,x t pri čemu se zadržavaju samo prva dva člana     , , , ... x t u u x t t u x t t t            (3.5) zamenom izvoda po vremenu dobija se     , , , ... x t F u x t t u x t t x            . (3.6) Primenom centralnih razlika i osrednjavanjem prvog člana sledi 1 1 1 11 2 2 n n n n j j j jn j u u F Ft u x          . (3.7) Primenom Lax-ove metode na pomeranje diskontinuiteta udesno, sa u=1 pre diskontinuiteta i u=0 posle diskontinuiteta, dobija se rešenje prikazano na Slici 3.1. Lokacija diskontinuiteta je tačno predviđena, ali se vidi disipativna priroda ovog metoda u razmazivanju diskontinuiteta preko nekoliko rastojanja između čvorova mreže. Ovo razmazivanje fronta postaje sve izraženije kako opada Courant-ov broj. 43 Slika 3.1 Numeričko rešenje Burgers-ove jednačine primenom Lax-ovog metoda. (Tannehill i dr., 1997) 3.2 Lax-Wendroff-ova metoda Lax-Wendroff-ova metoda (1960) je jedna od prvih metoda konačnih razlika drugog reda za rešavanje hiperboličkih parcijalnih diferencijalnih jednačina. Jednačina (3.3) se razvija u Taylor-ov red oko tačke  ,x t pri čemu se zadržavaju tri člana       2 2 2 , , , , ... 2x t x t tu u u x t t u x t t t t                    (3.8) Izvod prvog reda se može direktno zameniti iz diferencijalne jednačine. Sledi razvoj smene za drugi izvod. Ako se u jednačini (3.3) prebaci izvod po x sa desne strane znaka jednakosti i uradi izvod po vremenu jednačine, pri čemu je redosled diferenciranja F zamenjen sledi 2 2 2 2 u F F F t t x x t x t                       . (3.9) Pošto je F=F(u) sledi u F F u u A F Ft x u x x A F F u u t x A t u t t                                  (3.10) tako će drugi izvod u jednačini (3.8) posle smene (3.10) glasiti 2 2 . u F A t x x            (3.11) u x 44 Jacobian A sadrži samo jedan element za Burgers-ovu jednačinu, a A je matrica kada su u i F vektori. Zamenom prvog i drugog izvoda po u u jednačini Taylor-ovog reda (3.8) dobija se       2 , , ... 2 tF F u x t t u x t t A x x x                 (3.12) Primenom centralnih razlika sledi     2 1 11 1/2 1 1/2 1 1 2 2 n n j jn n n n n n n n j j j j j j j j F Ft t u u A F F A F F x x                       .(3.13) Jacobian-ova matrica je određena za polovinu intervala 1 1 1/2 1/2 . 2 2 j j j j j j u u u u A A i A A                   (3.14) Slika 3.2 Numeričko rešenje Burgers-ove jednačine primenom Lax-Wendroff-ovog metoda. (Tannehill i dr., 1997) Primenom Lax-Wendroff-ove metode je tačno i precizno pozicioniran položaj diskontinuiteta (Slika 3.2). Disperzna priroda ove metode je potvrđena prisutstvom oscilacija u blizini diskontinuiteta. Iako metoda koristi centralne razlike asimetrija se javlja zbog kretanja talasa. Veće oscilacije nastaju kada je Courant-ov broj jednak 0,6, nego za vrednost 1, odakle sledi da će kvalitet rešenja biti sve manji sa opadanjem vrednosti Courant-ovog broja. u x 45 3.3 MacCormack-ova metoda Metoda MacCormack-a (1969) predstavlja izmenjenu metodu Lax-Wendroff-a u dva koraka. Ovu metodu je mnogo lakše primeniti u odnosu na Lax-Wendroff-ovu metodu jer se ne javlja Jacobian. Primenom na neviskoznu Burgers-ovu jednačinu MacCormack-ova metoda je definisana sa sledeća dva koraka:     1 1 1 1 1 1 1 korak predviđanja, 1 korak popravke. 2 n n n n j j j j n n n n n j j j j j t u u F F x t u u u F F x                     (3.15) Član 1n ju  je privremeno pretpostavljena vrednost zavisno promenljive u u vremenskom trenutku n+1. Jednačina popravke daje konačnu vrednost zavisno promenljive u u trenutku n+1. U prvoj jednačini (3.15) uzeta je naredna razlika (u tačkama j+1 i j), a u drugoj jednačini (3.15) prethodna razlika (j i j-1). Ove razlike mogu da zamene mesta, što je kod nekih problema, posebno kod problema sa prostiranjem diskontinuiteta, povoljno uraditi. Slika 3.3 Numeričko rešenje Burgers-ove jednačine primenom MacCormack-ove metode. (Tannehill i dr., 1997) Rezultat primene ove metode je prikazan na Slici 3.3. Položaj fronta je dobro definisan. Za isti Courant-ov broj se dobijaju različita rešenja u odnosu na Lax- Wendroff-ovu metodu. Ovo je posledica zamene mesta razlika u koracima predviđanja i popravke i nelinearne prirode parcijalnih diferencijalnih jednačina održanja. u x 46 3.4 Rusanov (Burstein - Mirin) metoda Rusanov (1970) i Burstein i Mirin (1970) su istovremeno razvili eksplicitnu metodu u tri koraka.                       1 1/2 1 1 2 1 1 1/2 1/2 1 2 1 1 2 2 2 1 1 2 1 1 2 1 1 Prvi korak , 2 3 2 Drugi korak , 3 1 Treći korak 2 7 7 2 24 3 4 6 4 . 8 24 n n n n j j j j j n j j j j n n n n n n j j j j j j n n n n n j j j j j j j t u u u F F x t u u F F x t u u F F F F x t F F u u u u u x                                              (3.16) Poslednji član u trećem koraku, koji predstavlja četvrti izvod   4 4 4 u x x    , je dodat kako bi šema postala stabilna. U ovom članu se javlja slobodni parametar . Tačnost trećeg reda metode nije narušena pošto je greška zaokruživanja   4 O x    . Stabilnost Rusanov-e metode je osigurana ako je 1  , gde je t c x     Courant-ov broj, i 2 44 3     . Ako je 0  , odnosno ako ne postoji poslednji član u trećem koraku, ne može da bude zadovoljen uslov 0 1  . Modifikovana jednačina ove metode je       3 3 4 2 4 4 24 5 4 15 4 ... 120 t x xxxx xxxxx c x u cu u c x u                           (3.17) Umanjenje disipacije ove metode može se ostvariti izjednačavanjem člana koji množi četvrti izvod sa nulom, odnosno 2 44    , a umanjenje disperzione greške izjednačavanjem koeficijenta koji množi član petog izvoda sa nulom   2 24 1 4 / 5     . Primena Rusanov-e metode na Burgers-ovu jednačinu kretanja diskontinuiteta udesno daje rezultate prikazane na Slici 3.4. Veličina i položaj diskontinuiteta su tačno predviđeni, ali rezultati pokazuju veće vrednosti (prebacivanje) sa obe strane udarnog fronta. 47 Slika 3.4 Numeričko rešenje Burgers-ove jednačine primenom Rusanov-e metode. (Tannehill i dr., 1997) 3.5 Warming - Kutler - Lomax metoda Warming, Kutler i Lomax (1973) su razvili metodu trećeg reda koristeći razlike koje nisu centralne. Ova metoda koristi MacCormack-ovu metodu za prva dva koraka sračunata u 2 3 t , a treći korak je isti kao kod metode Rusanov-a. Prednost u odnosu na Rusanov-u metodu je što su samo vrednosti u integralnim čvorovima mreže potrebne za račun. Metoda Warming - Kutler - Lomax (WML) primenjena na jednačinu (3.3) postaje                       1 1 2 1 1 1 1 1 2 1 1 2 2 2 1 1 2 1 1 2 2 Prvi korak , 3 1 2 Drugi korak , 2 3 1 Treći korak 2 7 7 2 24 3 4 6 4 . 8 24 n n n j j j j n j j j j j n n n n n n j j j j j j n n n n n j j j j j j j t u u F F x t u u u F F x t u u F F F F x t F F u u u u u x                                             (3.18) Rezultati numeričkog rešavanja Burgers-ove jednačine dobijeni primenom WKL metode su prikazani na Slici 3.5. Rešenje je skoro isto kao kod Rusanov-a. Na osnovu dobijenih rezultata može se zaključiti da bilo koja metoda trećeg reda može da se koristi sa približno istom tačnošću. u x 48 Slika 3.5 Numeričko rešenje Burgers-ove jednačine primenom WKL metode. (Tannehill i dr., 1997) 3.6 Metoda trećeg reda sa podešavanjem parametra  Parametar , koji se javlja u trećem koraku metoda Rusanov i WKL, se može proizvoljno izabrati dok god se ne narušava stabilnost rešenja. Odabrani parametar  na početku proračuna, zadržava istu vrednost kroz celu mrežu. Ako se član numeričkog prigušivanja napiše u formi zakona održanja za treći korak 3 3 u x x          tada se  može menjati od tačke do tačke pri proračunu i održanje fluksa u mreži je osigurano. Primenom ovog pristupa poslednji član u metodama Rusanov ili WKL se može napisati kao    1/2 1/22 1 1 1 1 23 3 3 3 24 24 n n j jn n n n n n n n j j j j j j j ju u u u u u u u                 .(3.19) Vrednosti 1/2 n j  variraju na osnovu efektivnog Courant-ovog broja mreže. Warming i dr. (1973) su predložili da se ovi parametri računaju u svakoj tački mreže kako bi se smanjila disperzivna ili disipativna greška. Rezultati dobijeni primenom metode podešavanja dati su na Slici (3.6). Obe metode trećeg reda obezbeđuju zadovoljavajuće rešenje za slučaj minimalne disperzije. Nešto veće prebacivanje vrednosti nastaje sa leve strane diskontinuiteta, a skoro tačno rešenje se dobija na desnoj strani. u x 49 Slika 3.6 Numeričko rešenje Burgers-ove jednačine primenom metode podešavanja parametra . (Tannehill i dr., 1997) 3.7 Implicitne metode Vremenski centrirana implicitna metoda (trapezna metoda) se zasniva na jednačini       1 31 . 2 n nn n j j t t t u u u u O t               (3.20) Ako u jednačini (3.20) zamenimo izvode po vremenu primenom model jednačine, dobićemo 1 1 . 2 n n n n j j t F F u u x x                        (3.21) Ovo je nelinearni problem i potrebno je primeniti neku vrstu linearizacije ili neku iteracionu tehniku. Beam i Warming (1976) su predložili    1 1 1 n n n n n n n n nFF F u u F A u u u             (3.22) pa je  1 12 . 2 n n n n n j j j j t F u u A u u x x                      (3.23) Ako se izvodi po x zamene centralnim razlikama drugog reda sledi u x 50 1 11 1 1 1 1 1 1 1 1 1 1 4 4 . 2 4 4 n n j jn n n j j j n n n n n j j j j jn n j j tA tA u u u x x F F tA u tAt u u x x x                               (3.24) Jacobian A jednak je u za Burgers-ovu jednačinu, pa je moguće dalje pojednostavljenje desne strane jednačine. Primenjena linearizacija Beam i Warming-a vodi do linearnog sistema algebarskih jednačina u narednom vremenskom trenutku. Ovo je tridijagonalni sistem za čije je rešenje potrebno primeniti Thomas-ov algoritam. Ova metoda je stabilna za bilo koji vremenski korak. Modifikovana jednačina sadrži parne izvode članova. Veštačko ispravljanje je dodato šemi. Uobičajeno se razlika četvrtog reda  2 1 1 24 6 4 8 n n n n n j j j j ju u u u u          (3.25) može dodati jednačini (3.24) pri čemu se formalna tačnost metode neće promeniti. Prema Beam-u i Warming-u implicitna formula (3.24) sa dodatim eksplicitnim prigušivanjem je stabilna ako važi 0 1  . Na Slici (3.7) su prikazani rezultati primene vremenski centrirane implicitne formule na kretanje diskontinuiteta udesno. Rešenje bez prigušivanja je neprihvatljivo. Kada se doda eksplicitno prigušivanje dato jednačinom (3.25) dobijaju se bolji rezultati. Slika 3.7 Numeričko rešenje Burgers-ove jednačine primenom Beam Warming-ove (trapezoidalne) metode. (Tannehill i dr., 1997) u x 51 Kao dodatak trapezoidne formule Beam i Warming (1976) su razvili tehnike tri tačke unazad i Euler-ovu implicitnu metodu kao deo familije tehnika. Beam i Warming verzija Euler-ove implicitne šeme sledi iz Euler-ove formule unazad 1 1 n n n uu u t t           (3.26) koja je za nelinearnu jednačinu (3.3) 1 1 n n n Fu u t x           . (3.27) Ako se primeni isti postupak linearizacije dobija se 1 11 1 1 1 1 1 1 1 1 1 1 2 2 . 2 2 2 n n j jn n n j j j n n n n n j j j j jn n j j tA tA u u u x x F F tA tA ut u u x x x                               (3.28) To je tridijagonalni sistem koji se lako rešava. Ova šema je bezuslovno stabilna, ali se mora dodati prigušenje kao u jednačini (3.25) kako bi se osigurao upotrebljiv numerički rezultat. Jednostavniji oblik implicitnog algoritma se može dobiti ako se koristi zapis u delta obliku. Ovaj oblik ima prednost kod višedimenzionalnih problema jer obezbeđuje stacionarno rešenje koje je nezavisno od vremenskog koraka u problemima koji imaju stacionarno rešenje. Primenom delta oblika trapezoidalna formula jednačine (3.21) se može napisati u obliku 1 1 . 2 n n n n j j j t F F u u u x x                           (3.29) Lokalna linearizacija je primenjena za dobijanje 1n n n j j j jF F A u     . Konačni oblik diferencne jednačine postaje  1 11 1 1 1 . 4 4 2 n n j j n n j j j j j tA tA t u u u F F x x x                      (3.30) Ovaj oblik je mnogo jednostavniji od jednačine (3.24). Tridijagonalna forma je i dalje zadržana, ali desna strana ne zahteva množenje originalnog algoritma. To može biti od važnosti za sisteme jednačina kod kojih je veliki broj operacija. Rezultati 52 dobijeni primenom delta oblika na udarni talas koji se kreće udesno su dati na Slici 3.8. Rešenja sa i bez prigušivanja su u suštini identična onima dobijenim primenom razvijenog oblika. Rešenja neviskozne Burgers-ove jednačine primenom implicitne šeme su uglavnom lošija u odnosu na ona dobijena eksplicitnim tehnikama i zahtevaju više kompjuterskog vremena po koraku integracije. Kada su prisutni diskontinuiteti, rezultati dobijeni eksplicitnim metodama su bolji od rezultata implicitnih metoda koje koriste centralne razlike, stoga se eksplicitne metode preporučuju za rešavanje nestacionarne Burgers-ove jednačine. Slika 3.8 Numeričko rešenje dobijeno primenom vremenski centralne implicitne formule na kretanje diskontinuiteta udesno, delta oblik. (Tannehill i dr., 1997) 3.8 Šema Godunov-a Numeričke metode primenjene na rešavanje Burgers-ove jednačine, koje su prethodno izložene, koriste Taylor-ov red za dobijanje izraza za određivanje vrednosti zavisno promenljive u narednom vremenskom trenutku. Razlike u prostornom pravcu se takođe određuju primenom aproksimacija redovima uz zadovoljenje odgovarajuće tačnosti. Razvoj u Taylor-ov red daje dobre rezultate kada su zadovoljeni uslovi konvergencije redova. U slučaju metoda konačnih razlika smatramo da je razvoj reda prikladno sredstvo za dobijanje aproksimacije razlike i da je funkcija neprekidna, kao i da ima neprekidne izvode najmanje do reda aproksimirane razlike. To sigurno nije zadovoljeno u slučaju udarnog talasa ili nekog drugog diskontinuiteta. Godunov (1959) je prepoznao ovaj osnovni problem i predložio da se izbegne uslov diferencijabilnosti primenom aproksimacije konačnim zapreminama pri rešavanju jednačina održanja i u x 53 određivanja članova fluksa na granicama kontrolne zapremine rešenjem Riemann-ovog problema. Slika 3.9 Izgled kontrolne zapremine u x, t koordinatnom sistemu. (Tannehill i dr., 1997) Za eksplicitnu metodu kontrolna zapremina se prostire po osi t, od t do t +t, a po osi x, od x-x/2 do x+x/2. Ako posmatramo kontrolnu zapreminu sa centrom u(j,n+1/2) (Slika 3.9), rezultujuća numerička aproksimacija zavisno promenljive može biti napisana kao 1 1 1 2 2 n n j j j j t u u f u f u x                      . (3.31) U ovoj jednačini je vrednost promenljive u osrednjena po elementu zapremine tj.   2 2 1 , , x x j x x u u x t dx x        (3.32) a član fluksa je dobijen osrednjavanjem fluksa na granici kontrolne zapremine po vremenu 1 . t t t f fdt t     (3.33) Godunova metoda rešava lokalni Riemann-ov problem na svakoj granici kontrolne zapremine kako bi se dobila vrednost fluksa neophodna za obezbeđenje rešenja. Riemann-ov problem za Burgers-ovu jednačinu je   2 1 0 0, sa početnim uslovima ,0 . 02 j j u xu u u x u xt x                (3.34) 54 Geometrija problema je predstavljena na Slici 3.10. Osrednjene vrednosti zavisno promenljive imaju oblik stepenaste promene raspodele, koja vodi do pojave diskontinuiteta na svakoj granici kontrolne zapremine. Na svakoj granici kontrolne zapremine nastaje udar ili ekspanzija, koji propagiraju sa vremenom. Slika 3.10 Talasni dijagram za Godunov metod. (Tannehill i dr., 1997) Primenom zapisa 1 1/2 1/2 2 j j j j u udx c dt           (3.35) rešenje Riemann-ovog problema za ovu jednačinu može biti napisano za slučaj udarnog i ekspanzionog talasa u tabeli 3.1. Tabela 3.1 Slučaj 1: Udarni talasi 1 1/2 1 1/2 2 1/2 1/2 2 1 1/2 / / 1 / 2 0 1 / 2 0 j j j j j j j j j j j u u u x t c u u x t c u c f u c                   Slučaj 2: Ekspanzioni talasi 1 1 1 1 1 2 1/2 1/2 1 2 1 1/2 1 / / / / 0 0 1 / 2 0 0 1 / 2 0 0 j j j j j j j j j j j j j j j j j j j u u u x t u u x t u x t u u x t u u u f u c u u u c u u                                 55 Usvojena je pretpostavka da nema uzajamnog dejstva talasa iz dve susedne kontrolne zapremine. Ova pretpostavka je potrebna zbog pisanja jednostavnog rešenja za stanje promenljivih na granicama kontrolne zapremine i osigurana je samo onda kada talas može da se kreće najviše za rastojanje polovine kontrolne zapremine. Kao posledica ovoga uslov stabilnosti Godunov-e šeme je max 1 . 2 t u x    (3.36) Primenom jednačine (3.31) možemo dobiti rešenje Burgers-ove jednačine gde je fluks određen rešavanjem Riemann-ovog problema. Tipičan rezultat je predstavljen na Slici 3.11. Rešenje Burgers-ove jednačine primenom Godunov-e metode se lako postiže i dobijaju se odlični rezultati. Kada se ova metoda primeni na rešavanje jednačina održanja, koje opisuju strujanje fluida, neophodno je upotrebiti iterativni postupak rešavanja koji može biti vremenski zahtevan. Slika 3.11 Primena Godunov-e metode na problem udarnog talasa. (Tannehill i dr., 1997) 3.9 Roe-ova šema U rešavanju nelinearnih sistema, umesto da se primeni tačna nelinearna iterativna šema, uspešno je primenjeno rešavanje približnog Riemann-ovog problema. Najpopularniju aproksimaciju Riemann-ovog rešavanja je predložio Roe (1980,1981). On je predložio rešavanje linearnog problema 0, u u A t x       (3.37) u x 56 gde je A matrica koeficijenata koja zavisi od lokalnih uslova. U slučaju Burgers-ove jednačine matrica A je jedan skalarni element. Primenom Roe-ove šeme na Burgers- ovu jednačinu, u je konstanta i označava osrednjenu vrednost od A . Tada se posmatra problem 0, u u u t x       (3.38) gde je u brzina određena za kontrolne zapremine j i j+1 iz prvog uslova. Ovo daje 1 1/2 1 , j j j j j F F u u u u        (3.39) što se redukuje za Burgers-ovu jednačinu u obliku  1 1 1/2 1 / 2 . j j j j j j j j u u u u u u u u           (3.40) Primenom (3.40) može se razviti numerički fluks na osnovu Rankine - Hugoniot-ove relacije koja daje direktno vezu između skoka fluksa i skoka zavisno promenljive u preko talasa  1/21 1 .jj j j jF F u u u    (3.41) Približno Riemann-ovo rešavanje prepoznaje samo diskontinuitete, ali ne može da napravi razliku između ekspanzije i udarnog talasa. Potrebno je uvesti popravku. Slika 3.12 Talasni dijagram za Roe-ovu metodu primenjenu na Burgers-ovu jednačinu. (Tannehill i dr., 1997) 57 Približno Riemann-ovo rešavanje Burgers-ove jednačine predstavljeno je na Slici 3.12, gde vidimo da samo jedan talas proizilazi iz granice kontrolne zapremine. Ovaj talas se prostire bilo u pozitivnom ili u negativnom smeru u zavisnosti od 1/2 /ju dx dt  . Primenom definicije skoka preko talasa sledi  1/21/2 1jj j j jf F u u u      (3.42) ili  1/21 1/2 1 .jj j j jF f u u u       (3.43) Numerički fluks se može napisati u simetričnom obliku   1 1/2 1/21/2 11 . 2 2 j j j jj j j F F f u u u u           (3.44) Potrebno je razmotriti pojedinačne doprinose koji se javljaju kada se talas prostire ili sa desne ili sa leve strane. Ako se svaki slučaj proceni pojedinačno, numerički fluks se može predstaviti jednom jednačinom u obliku  1 1/21/2 1 1 . 2 2 j j jj j j F F f u u u        (3.45) Slika 3.13 predstavlja proračun propagacije diskontinuiteta primenom Roe-ove šeme i rezultati pokazuju odlično slaganje sa analitičkim rešenjem. Slika 3.13 Propagacija diskontinuiteta primenom Roe-ove šeme. (Tannehill i dr., 1997) Roe-ova šema ima uslov stabilnosti, kao i sve eksplicitne metode, da Courant-ov broj mora da bude manji od 1. Roe-ova formulacija propagira razliku zavisno promenljivih između dve tačke kao da diskontinuitet postoji između te dve tačke. Kao posledica ovoga mogu se javiti ekspanzioni talasi koji nisu posledica fizičke pojave. U u x 58 Burgers-ovoj jednačini ovo je problem kada 1/2ju  nestaje. To se naziva zvučni prelaz. Klasični primer je rešenje dobijeno primenom Roe-ove šeme sa ulaznim podacima 0 0 1 0 . 1 x x u x x L         (3.46) Analitičko rešenje je centralna ekspanzija oko x=x0. Rešenje primenom Roe-ove šeme dato je na Slici 3.14. Verno reprodukovani početni uslovi pokazuju netačan stacionarni ekspanzioni udar. Slika 3.14 Ekspanzioni talas primenom Roe-ove šeme. (Tannehill i dr., 1997) Ovo nefizičko ponašanje nastaje zbog toga što šema ne može da napravi razliku između ekspanzionog i kompresionog udara. Svaki od njih je važeće rešenje. Navodi se da postojanje ekspanzionog udara narušava stanje entropije dozvoljavajući netačno fizičko ponašanje. Oleinik (1957) i kasnije Lax (1973) su razvili uslove koji moraju da budu zadovoljeni od strane diskontinualnih rešenja hiperboličkih jednačina. Najjednostavnija formulacija ovih uslova primenjena na Burgers-ovu jednačinu se može napisati kao 0 ,R L x dx u u dt        (3.47) gde su uR i uL vrednosti zavisno promenljivih levo i desno od diskontinuiteta. Harten i Hyman (1983) su predložili sledeću modifikaciju numerike (modifikaciju 1/2ju  ) kako bi se dobila prihvatljiva rešenja koja obuhvataju odgovarajuće fizičko ponašanje. Neka je 1 max 0, 2 j ju u         (3.48) u x 59 tada je 1/2 1/2 1/2 1/2 . j j j j u u u u             (3.49) U slučaju kompresije ( 0  ) koristi se nepromenjena šema za propagaciju diskontinuiteta, dok se u slučaju ekspanzije  1 / 2j ju u   zahteva modifikacija. Slika 3.15 predstavlja rezultate ekspanzije (-1,1) primenom Roe-ove metode sa popravkom entropije. Slaganje numeričkih i analitičkih rešenja je dobro i približno je onome koje se očekuje od metode prvog reda. Slika 3.15 Ekspanzija primenom Roe-ove šeme sa korekcijom entropije. (Tannehill i dr., 1997) u x 60 NUMERIČKA METODA ZA PRAĆENJE PROSTIRANJA FRONTA 4. SKALARNE PROMENLЈIVE DUŽ KARAKTERISTIČNOG PRAVCA U ovom radu je primenjena nova numerička metoda za rešavanje transportne jednačine skalarne strujne promenljive veličine (entalpije, temperature, zapreminskog udela parne faze, koncentracije) u jednodimenzionalnom prostoru. Metoda se zasniva na transformaciji transportne jednačine skalarne promenljive u oblik sa materijalnim izvodom i aproksimacije materijalnog izvoda sa konačnim razlikama duž karakterističnog pravca određenog kretanjem fluidnog delića. Početne vrednosti skalarnog parametra u tačkama iz kojih polaze karakteristike su određene primenom Lagrange-ovog interpolacionog polinoma trećeg reda tačnosti. Ova metoda efikasno smanjuje numeričku difuziju i omogućava praćenje kretanja razdelne površine u dvofaznom strujanju primenom standardne metode. Tačnost metode je demonstrirana na testovima predviđanja nestacionarne granice ključanja i problema prostiranja diskontinuiteta (Stevanović i Jovanović, 2000). Analizirana je zavisnost hibridne metode od prostornog koraka integracije i stepena Lagrange-ovog interpolacionog polinoma. Lagrange-ov interpolacioni polinom 3. stepena je preporučen jer se dobija praktično tačno rešenje sa prihvatljivim brojem čvorova. Jednodimenzionalno jednofazno ili homogeno višefazno strujanje je opisano zakonima održanja mase i energije i zakonom promene količine kretanja u sledećem obliku:  zakon održanja mase   0, u t x       (4.1)  zakon promene količine kretanja  2 *( ) sin ,w uu p g t x x               (4.2)  zakon održanja energije   ( ) sin , e ue p gu q t x t               (4.3) gde je ukupna energija fluidne struje određena sa 61 2 . 2 u e h  (4.4) Konvektivni članovi u jednačinama zakona održanja mase i energije (drugi članovi sa leve strane jednačina (4.1) i (4.3)) su napisani u tzv. konzervativnom ili divergentnom obliku, kao skalarni proizvod gradijenta jediničnog vektora i vektora brzine. Za nestišljivo strujanje transformacijom ovih članova u nekonzervativni oblik i smenom (4.1), (4.2) i (4.4) u zakon održanja energije (4.3), uz zanemarivanje uticaja gravitacije, sledi * 1 ,z V h h p p u u u q t x t x                     (4.5) uvođenjem materijalnog izvoda , D u Dt t x       (4.6) dobija se *1 .z V Dh Dp u q Dt Dt           (4.7) Jednačina (4.7) se primenjuje na strujanja tople vode u sistemima daljinskog grejanja, kao što su startovanje postrojenja i kretanje temperaturskih talasa ka potrošačima, promena snage toplotnog izvora ili potrošnje ili zaustavljanje rada izvora toplote. Predviđanje propagacije fronta talasa entalpije u složenim cevnim mrežama sistema daljinskog grejanja je neophodno za postavku operacionih procedura, koje će omogućiti vremenski tačno i optimalno, prema toplotnom zahtevu, snabdevanje potrošača toplotom (Stevanović i dr., 2006a). Pri prinudnom strujanju fluida prenošenje toplote mehanizmom prelaženja preovladava nad provođenjem toplote kroz fluid, kao i nad mehanizmom turbulentne difuzije. Za takve strujne uslove uticaj promene pritiska na entalpiju i disipacija strujne energije usled trenja na zidovima cevovoda (prvi i drugi član na desnoj strani jednačine (4.7) sledstveno) se mogu zanemariti u odnosu na promenu entalpije usled hlađenja ili zagrevanja fluida kroz zidove cevovoda (treći član na desnoj strani jednačine 4.7). 62 Promena gustine može da nastane kao posledica isparavanja ili kondenzacije u struji homogenog fluida. Tako energijska jednačina (4.7) postaje talasna jednačina prvog reda ,V Dh q Dt   (4.8) koja opisuje transport tople vode entalpije h duž strujnog kanala, odnosno cevovoda, brzinom prostiranja fluidnih delića u. U opštem obliku možemo napisati transportnu jednačinu skalarne promenljive , D S Dt    (4.9) gde je  neka skalarna promenljiva (entalpija, temperatura, zapreminski udeo parne faze, koncentracija), a S izvorni član. Jednačina (4.9) je parcijalna diferencijalna jednačina hiperboličkog tipa. Uzimajući u obzir fizičko značenje jednačine (4.9), direktan postupak njenog rešavanja je integracija duž karakterističnog pravca određenog kretanjem fluidnog delića, koji je u x-t koordinatnom sistemu određen sa 1 . dt dx u  (4.10) Ovakav postupak rešavanja je svojstven metodi karakteristika, pri čemu se on uglavnom primenjuje pri rešavanju bilansa količine kretanja i mase duž karakterističnih pravaca određenih prostiranjem talasa pritiska. Materijalni izvod D/Dt u jednačini (4.9) duž karakterističnog pravca definisanog jednačinom (4.10) se transformiše u izvod po vremenu pa jednačina (4.9) dobija oblik . d S dt    (4.11) Izvod po vremenu u jednačini (4.11) se aproksimira konačnim razlikama, tako da se dobija algebarska jednakost  ( ) ( ) ,i P Pt t t S t     (4.12) gde indeks i označava i-ti čvor na Slici 4.1, dok je P tačka preseka karakterističnog pravca i x koordinate u vremenskom trenutku t. Koordinata xP tačke P je određena na osnovu nagiba karakterističnog pravca i linearne interpolacije brzine između čvorova i i i-1 za strujanje u pozitivnom smeru uP0, odnosno između čvorova i i i+1 za strujanje u 63 negativnom smeru uP0. Sledeće relacije se mogu napisati ,i P P x x u t    (4.13) 1 1 1 0,P i P i i i i i i u u x x dt za u u x x dx u         (4.14) 1 1 1 0.P i P i i i i i i u u x x dt za u u x x dx u         (4.15) Slika 4.1 Numerička mreža sa obeleženim čvorovima (i-2), (i-1), i, (i+1) koji se koriste za formiranje Lagrange-ovog interpolacionog polinoma 3. stepena u slučaju strujanja u smeru x-ose Iz jednačina (4.13), (4.14) i (4.15) dobija se izraz za određivanje koordinate tačke P , 1 i P i u t x x a     (4.16) gde je 1 1 1 za 0,i i i i i u u dt a t x x dx u         (4.17) 1 1 1 za 0.i i i i i u u dt a t x x dx u         (4.18) 64 Za sračunavanje skalarne promenljive u čvoru i u vremenskom trenutku t+t, to jest i(t+t) pomoću jednačine (4.12), potrebno je poznavati i njenu vrednost u tački P, P(t). Za određivanje ove vrednosti razvijen je poseban postupak, koji obezbeđuje numeričku šemu višeg reda tačnosti i u velikoj meri eliminiše numeričku difuziju (Stevanović i dr., 2007c). Navedeni postupak omogućuje pouzdano i efikasno predviđanje prostiranja talasa entalpije u okviru cevne mreže sistema daljinskog grejanja. Skalarna promenljiva P(t) se određuje primenom Lagrange-ovog interpolacionog polinoma 3. stepena. U (Stevanović i Jovanović, 2000) su primenjene i linearna interpolacija i Lagrange-ov polinom 2. stepena kako bi se utvrdila zavisnost tačnosti modela i stepena polinoma, i pokazana je prikladnost primene polinoma 3. stepena. Za formiranje Lagrange-ovog interpolacionog polinoma 3. stepena potrebno je poznavati vrednosti skalarne promenljive veličine u 4 čvora numeričke mreže. Izbor ovih čvorova zavisi od smera strujanja. Koristi se jedan čvor nizvodno od posmatranog čvora, dva čvora uzvodno i posmatrani čvor. U slučaju strujanja u smeru x-ose koriste se vrednosti promenljive u čvorovima numeričke mreže (i-2), (i-1), i, (i+1), a u slučaju strujanja u negativnom smeru (i-1), i, (i+1), (i+2). Kada se proračun radi na granicama sistema, za odgovarajući granični uslov, uzimaju se ona četiri čvora koja u tom slučaju postoje u fizičkom domenu cevi. Npr. ako je i2 i strujanje u pozitivnom smeru, i-2 ne postoji u cevovodu pošto se čvorovi numerišu rastućim nizom prirodnih brojeva. U tom slučaju se uzimaju čvorovi 1, 2, 3 i 4. Slično se dešava i na kraju cevi kada je i  n[j] i strujanje u negativnom smeru, pri čemu se uzimaju čvorovi n[j]-3, n[j]-2, n[j]-1 i n[j]. Primenjen je sledeći algoritam proračuna pogodan za kompjutersko programiranje Lagrange-ovog interpolacionog polinoma. Predstavljena je interpolacija između čvorova (i-2), (i-1), i, (i+1). 1 0 ( ) ( ) m j m P m P j j x x D        , (4.19) gde je parametar m=3 stepen Lagrange-ovog polinoma, a j=0, 1, … , m broj čvorova između kojih se vrši interpolacija. Proizvod razlika x koordinata tačke P i m+1 čvora između kojih se vrši Lagrange-ova interpolacija je 65 1 0 1( ) ( )( ) ... ( )m P P P P mx x x x x x x       , (4.20) 0 1 1 1( )( ) ... ( )( )( ) ... ( ), j 0,1,...,m.j j j j j P j j j j mD x x x x x x x x x x x x             (4.21) Greška interpolacije Lagrange-ovog interpolacionog polinoma je data u (Demidovitch i Maron, 1987) 1 1( ) ( ) ( ) ( 1)! m m m M x x x m       , (4.22) gde je    11 max , m m a x b M x       (4.23) a [a,b] je interval interpolacije po x-osi. Prema jednačinama (4.22), (4.20) i (4.23) greška zaokruživanja za numeričku diskretizaciju skalarne promenljive  duž ose x primenom Lagrange-ovog interpolacionog polinoma na mreži predstavljenoj na Slici 4.1 je O[(x)m+1] (u slučaju uniformne mreže kada je prostorni korak konstantan). Dakle, za polinom 3. stepena greška je četvrtog reda, O[(x)4]. Integracija jednačine (4.11) po vremenu je sprovedena duž karakterističnog pravca primenom Euler-ove eksplicitne metode (4.12). Dakle, integracija po vremenu ima grešku prvog reda tačnosti, O(t). Analize stabilnosti rešavanja pokazuju da metoda karakteristika zahteva vremenski korak integracije koji zadovoljava Courant-ov kriterijum (Tannehill i dr., 1997) min , 2,3,..., 1. i x t i n u          (4.24) Ako je kretanje fluidnog delića u toku jednog vremenskog koraka integracije jednako prostornom koraku 1i i ix x x    , tada karakteristika povezuje tačke (i-1,t), i (i,t+t) u numeričkoj mreži, odnosno tačka P se poklapa sa tačkom (i-1,t), Slika 4.1. Zbog toga nema potrebe da se interpolacijom traži vrednost promenljive u tački P pošto je P(t)= i-1(t). U ovom slučaju metoda karakteristika predviđa vrednost i(t+t) sa greškom samo zbog vremenske integracije jednačine (4.11). 66 Ako je entalpija fluida skalarna promenljiva  , onda zapreminski toplotni fluks Vq , koji je deo izvornog člana, može da ima zadatu vrednosti ili pošto je vrednost entalpije u tački P poznata može da se odredi primenom zakona prostiranja toplote kroz zidove strujnog kanala. Gustina fluida u tački P se može odrediti iz jednačine stanja u obliku  ,h p  . Iako su jednačine zakona održanja izvedene za nestišljivo strujanje uzimanjem u obzir zavisnosti gustine od pritiska u jednačini stanja omogućeno je predviđanje efekata stišljivog strujanja za strujanja sa malim Mach-ovim brojem (Patankar, 1980). Slika 4.2 Predviđanje propagacije diskontinuiteta (Stevanović i Jovanović, 2000) Sledi primer rešavanja problema prostiranja diskontinuiteta. Strujni parametar h ima vrednost h=0 u početnom trenutku, a za 0t  h=1 na ulazu u strujni kanal. Ovaj problem je rešen primenom metode SIMPLE (Patankar, 1980), i hibridne metode (Stevanović i Jovanović, 2000). Hibridna metoda predstavlja kombinaciju metode SIMPLE, primenjene na rešavanje jednačine zakona održanja mase i zakona promene količine kretanja, i metode karakteristika primenjene na rešavanje zakona održanja energije. Rezultati prikazani na Slici 4.2 su dobijeni sa prostornim korakom 1 mx  i konstantnom brzinom 1 m/su  , a desna strana jednačine (4.5) je jednaka nuli. Na Slici 4.2 se vidi da numerička šema metode SIMPLE nije adekvatna za predviđanje prostiranja diskontinuiteta. To je numerička šema prvog reda koja dovodi do izravnavanja fronta diskontinuiteta zbog grešaka numeričke disipacije (Tannehill i dr., 67 1997). Hibridna metoda, koja jednačinu (4.5) rešava primenom metode karakteristika ima mogućnost predviđanja propagacije diskontinuiteta. Tačno rešenje je dobijeno primenom vremenskog koraka /t x u   . Ovo je moguće zato što je prostorni korak x konstantan i rešenje jednačine (4.5) za sve čvorove je duž karakterističnog pravca koji počinje u čvoru (i-1,t), Slika 4.1, i završava u čvoru (i,t+t) , tj. tačka P se poklapa sa čvorom (i-1,t). U ovom slučaju nema potrebe za interpolacijom parametra h u trenutku t. U inženjerskim proračunima ovakav uslov ne može da bude zadovoljen zbog obe promenljive i prostornog koraka x i brzine u. U ovakvim situacijama tačka P na Slici 4.1 se nalazi između čvorova (i-1,t) i (i,t) za ui0, ili između (i,t) i (i+1,t) za ui<0. U cilju zadovoljenja ovih opštih uslova propagacija diskontinuiteta je računata i sa vremenskim koracima  0,75 /t x u   i  0,5 /t x u   . U ovom slučaju diskontinuitet je izobličen zbog disperzione greške, kao što je i očekivano za numeričke šeme drugog reda (Tannehill i dr., 1997). Oblik diskontinuiteta je sačuvan i njegova lokacija je adekvatno predviđena. 4.1 Granični uslovi S obzirom da je svrha dinamičkog modela simulacija i analiza prelaznih termičkih procesa u složenim cevnim mrežama sistema daljinskog grejanja, potrebno je definisati i granične uslove u spoju u čvoru cevne mreže. Veze potrošača sa razvodnim i povratnim cevovodima unutar mreže sistema daljinskog grejanja su ostvarene spajanjem tri ili više cevovoda, kao što je prikazano na Slici 4.3. Cevi sa ili bez pumpi i toplotno-razmenjivačkih stanica u kojima fluid teče ka čvoru su obeležene sa Ji (i = 1,2,......,n), dok su sa Ij (j = 1,2,......,m) obeležene cevi kojima fluid teče od čvora. Sa D je obeležena tačka spoja, a sa D' tačka ispred ili iza toplotno-razmenjivačke stanice, odnosno pumpe. Promena entalpije pri prolasku kroz toplotno-razmenjivačku stanicu ili pumpu je određena sa hJi, odnosno hIi, gde indeksi označavaju uticanje ili isticanje iz spoja cevi. Promena entalpije duž karakterističnog pravca koji stiže u tačku D' u trenutku t+t u cevi Ji je određena sa 68 Slika 4.3 Spoj više cevi u čvoru ', , , ( ) ( ) 1,2,......, . i i i D J P J P J q h t t h t t i n             (4.25) Entalpija fluida na mestu spoja, D, određena je mešanjem entalpija fluida koji utiču u čvor cevima Ji (i = 1,2,......,n), dakle stacionarni energetski bilans strujanja fluida od tačaka D' u svim cevovodima do tačke mešanja D, uzimajući u obzir promene entalpije u pumpama i razmenjivačima toplote daje jednačinu  , , , 1 , , 1 , i i i i i i i i n D J D J D J J J i D n D J D J J i u h h A h u A               (4.26) gde imenilac predstavlja ukupan maseni protok fluida ka spoju D. Površina protočnog preseka cevi Ji je označena sa iJ A , a gustina fluida u tački D' u cevi Ji sa , iD J  . Energetski bilans za svaku toplotno-razmenjivačku stanicu ili pumpu u cevi Ij, daje , 1,2,..., .i iD J D Jh h h j m    (4.27) Pomoću jednačine (4.25) određuje se entalpija ispred toplotno-razmenjivačke stanice ili pumpe u cevima koje vode fluid ka čvoru, a sračunavanjem (4.26) i zatim (4.27) određuje se entalpija neposredno iza razmenjivača toplote ili pumpe u cevima koje odvode fluid. J1 Ji Jn I1 Ij Im (D’,J1) (D’,Ij) D 69 4.2 Hidraulički model Nestacionarno prostiranje temperaturskih talasa u sistemu daljinskog grejanja je spor proces, koji traje od nekoliko desetina minuta do nekoliko sati. Za vreme ovog perioda voda u cevima se može smatrati nestišljivom u odnosu na promene pritiska zbog brzine zvuka kojom se prostire talas pritiska koja je reda veličine 1000 m/s (promene pritiska u sistemu će ovom brzinom biti transportovane i prigušene unutar mreže u mnogo kraćim vremenskim intervalima nego što je period nestacionarnog prostiranja temperaturskih talasa). Ova pretpostavka omogućava sprovođenje hidrauličkog proračuna kao niza stacionarnih stanja primenom stacionarnog hidrauličkog modela strujanja u mreži. Stacionarni hidraulički model primenjen za određivanje pritisaka i masenih protoka unutar cevne mreže se zasniva na prstenastom modelu mreže koji je dat u (Stevanović i dr., 2006b; Stevanović i dr., 2007c). 70 5. REZULTATI SNIMANJA PRELAZNOG PROCESA U SISTEMU DALJINSKOG GREJANJA Validacija programa za realne uslove rada sistema daljinskog grejanja u prelaznim režimima, izvedena je poređenjem dobijenih numeričkih rezultata sa rezultatima merenja vremenskih promena temperature, pritiska i protoka tokom karakterističnih prelaznih režima u odabranom sistemu daljinskog grejanja, kao što je uspostavljanje cirkulacije u sistemu i prostiranje temperaturskog talasa iz izvora toplote ka potrošačima. Merenja dinamičkih promena temperatura i protoka obavljena su u okviru sistema daljinskog grejanja toplane “Zemun“ pri promeni snage izvora toplote (Stevanović i dr., 2006a). 5.1 Prelazni temperaturski režim u sistemu daljinskog grejanja toplane “Zemun“ Snimanje prelaznog procesa je sprovedeno u okviru magistrale 3 sistema daljinskog grejanja toplane “Zemun“ (Slika 5.1). Sistem daljinskog grejanja toplane “Zemun“ se sastoji od tri magistrale. Magistrala 1 (Slika 5.2) i magistrala 2 (Slika 5.3) snabdevaju potrošače toplotom za grejanje, dok magistrala 3 obezbeđuje i toplotu za pripremu sanitarne tople vode. Ukupna snaga toplane je 65 MW, pri čemu se oko 2/3 toplote predaje potrošačima u okviru magistrale 3. Magistrala 3 obuhvata preko pedeset podstanica, većim delom u stambenim zgradama. Podstanice imaju regulacione ventile pomoću kojih se reguliše protok shodno potrebama potrošača. Merene su temperature u razvodnim i povratnim granama u toplani i u tri podstanice PS A (Banijska 30), PS B (Banijska 1) i PS C (Mostarska 12, SC “Partizan“) (označene brojevima 117, 108 i 206 na Slici 5.1). Podstanice PS A i PS C su u većim stambenim zgradama, dok je PS B u okviru sportskog centra. Dužina toplovoda od toplane do PS A je 729 m, do PS B 1615 m, a do PS C je 1718 m. Takođe, izmereni su i protoci kroz podstanice i ukupni protok iz toplane ka magistrali 3. Podstanica PS A je potrošač toplote samo za grejanje, dok PS B i PS C troše toplotu za grejanje i pripremu potrošne tople vode. 71 Slika 5.1 Magistrala 3 sistema daljinskog grejanja toplane “Zemun“. 72 Slika 5.2 Magistrala 1 sistema daljinskog grejanja toplane “Zemun“. Slika 5.3 Magistrala 2 sistema daljinskog grejanja toplane “Zemun“. Merenje promena temperatura i protoka sa porastom snage izvora toplote u podstanicama i u toplani je obavljeno dana 02.03.2007. u vremenu od 930-1200 časova. U navedenom vremenskom intervalu praćene su i promene vrednosti pritisaka, temperatura i protoka u toplani, dok je merenje prelaznog režima sa smanjenjem snage izvora toplote obavljeno dana 22.03.2007. u vremenskom periodu od 1100-1400 časova u istim podstanicama kao prilikom prethodnog merenja (Banijska 1, Banijska 30 i Mostarska 12) i u toplani. Merenja su sprovedena korišćenjem postojeće instrumentacije industrijske tačnosti (Stevanović i dr., 2006a). Protok je meren ultrazvučnim meračima čija je relativna greška merenja ±2%, dok je greška merenja temperature manja od 1 0 C. Merenjima su dobijene dve grupe rezultata. I.25 I.19 I.18 I.17 I.16 I.15 I.13 I.12I.11 I.10 I.9 I.8 I.7 I.6 I.5 I.4 I.3 I.2 I.1 I.24 I.23 I.22 I.21 I.20 TO II.8 II.7 II.6 II.5 II.17 II.16 II.15 II.14 II.13 II.12 II.11 II.10 II.9 II.4 II.3 II.2 II.1 TO 73 Prva grupa rezultata pokazuje promene temperature i protoka u podstanicama sa porastom snage izvora toplote. Merenja su započeta prepodne u 9:30 časova pri visokoj spoljnoj temperaturi od +15 0 C. Rezultati merenja temperatura u razvodnim i povratnim vodovima su prikazani na Slici 5.4. Slika 5.4 Izmerene temperature u razvodnom i povratnom vodu u toplani i tri podstanice PS A, PS B i PS C za vreme i nakon porasta snage izvora toplote. U periodu od 14. do 32. minuta raste snaga toplotnog izvora u toplani (Slika 5.4) i shodno tome iz toplane ka potrošačima kreće temperaturski talas sa maksimalnom temperaturom vode od 70 0 C. Nakon vremenskih intervala koji su uslovljeni brzinom strujanja vode, talas stiže najpre u najbližu podstanicu PS A, a zatim u PS B i najdalju PS C. Maksimalna vrednost talasa temperature koji stiže do podstanica je niža za oko 2 0C. Vrednosti temperatura u povratnim granama praktično ne pokazuju uticaj povećanja snage toplotnog izvora, što ukazuje na veliku akumulacionu sposobnost objekata koji se greju. Jedino u PS A dolazi do porasta povratne temperature u 100-tom minutu, ali to je posledica potpunog zaustavljanja protoka kroz podstanicu, s obzirom da pri visokoj spoljnjoj temperaturi nije bilo više potrebe za grejanjem. Dinamičke promene temperature na izlazu iz PS C su posledica potrošnje tople vode. 74 Protok kroz podstanice za vreme prelaznog temperaturskog režima nastalog usled porasta snage izvora toplote je prikazan na Slici 5.5. Ukupan protok iz toplane ka pojedinim magistralama je prikazan na Slici 5.6. U 90-tom minutu je obustavljeno grejanje i protok kroz magistrale 1 i 2. U ovim magistralama se ne priprema sanitarna topla voda. Slika 5.5 Protoci kroz podstanice za vreme i nakon porasta snage izvora toplote. Slika 5.6 Protok kroz magistrale 1, 2 i 3 meren na izlazu iz toplane za vreme i nakon porasta snage izvora toplote. 75 Druga grupa rezultata pokazuje promene temperature i protoka u podstanicama sa smanjenjem snage izvora toplote. Merenja su započeta prepodne u 11:00 časova pri spoljnoj temperaturi od +12 0 C. Spoljna temperatura je u toku merenja rasla, tako da je do 14:00 časova kada su merenja završena iznosila 14 0C. Rezultati merenja temperatura u razvodnim i povratnim vodovima su prikazani na Slici 5.7. Slika 5.7 Izmerene temperature u razvodnom i povratnom vodu u toplani i tri podstanice PS A, PS B i PS C za vreme i nakon smanjenja snage izvora toplote. U periodu od 58. do 100. minuta pada snaga toplotnog izvora u toplani (Slika 5.7). Pad temperature tople vode u posmatranim podstanicama javlja se prvo u PS A i to 24 min po snižavanju temeprature vode u toplani. U podstanicama PS C i PS B se pad temperature primećuje za 48 min, odnosno 52 min. Temperatura vode u toplani je u 86 min 42 0 C. Nakon vremenskih intervala koji su uslovljeni brzinom strujanja vode, talas stiže najpre u najbližu podstanicu PS A, a zatim u PS B i PS C. Vrednosti temperatura u povratnim granama iz podstanica sa potrošnjom sanitarne tople vode se snižavaju u roku 2-3 min od sniženja temperature u dolaznoj grani, dok se u podstanici iz koje se distribuira samo topla voda za grejanje (PS A) pad povratne temperature primećuje posle 26 min. Do porasta odlazne temperature dolazi u 100-tom minutu. Iste vrednosti vremenskih intervala su potrebne da bi se u navedenim podstanicama primetio porast temperature. 76 Protok kroz podstanice za vreme prelaznog temperaturskog režima nastalog usled smanjenja snage izvora toplote je prikazan na Slici 5.8. Do naglog porasta protoka u podstanici PS A, koja snabdeva potrošače samo toplotom za grejanje, dolazi u 92 min (10 min pošto je temperatura dolazne tople vode počela da pada). U 110 min je dostignuta maksimalna vrednost protoka kroz ovu podstanicu. Temperatura dolazne tople vode pada do 116 min, ali dalje povećanje protoka nije moguće jer se ventil koji služi za regulaciju protoka već potpuno otvorio. U 148 min je temperatura dolazne tople vode 58 0C (vrednost koja se zadržava u narednih 20-tak min), a u 158 min počinje pritvaranje ventila. Značajnije promene protoka u preostalim dvema podstanicama se ne primećuju, što se može objasniti i smanjenom potrošnjom sanitarne tople vode zbog njene snižene temperature. Slika 5.8 Protoci kroz podstanice za vreme i nakon smanjenja snage izvora toplote. Ukupan protok iz toplane ka pojedinim magistralama je prikazan na Slici 5.9. U 106-tom minutu se primećuje porast protoka u magistrali 3. U magistralama 1 i 2 nema značajnijeg porasta protoka. Ove magistrale snabdevaju znatno manji broj potrošača od magistrale 3 i na njima se ne nalaze značajniji i veći potrošači kao na magistrali 3 (Sportski centar, Škola za decu oštećenog vida, Škola za decu oštećenog sluha, Dom zdravlja), pa je njihov toplotni zahtev manji. Zbog ovoga je i manji pad temperature 77 tople vode u povratnom vodu u ovim dvema magistralama, pa je i manja promena protoka tople vode potrebana da bi se zadovoljile potrebe potrošača. Slika 5.9 Protok kroz magistrale 1, 2 i 3 meren na izlazu iz toplane za vreme i nakon smanjenja snage izvora toplote. 78 6. NUMERIČKE SIMULACIJE PRELAZNIH REŽIMA RADA SISTEMA DALЈINSKOG GREJANJA U ovom poglavlju su prikazani rezultati numeričke simulacije prostiranja temperaturskog talasa pri prelaznim režimima rada sistema daljinskog grejanja, koji su prikazani u poglavlju 5. Izmerene vrednosti promene pritiska, temperature, protoka i snage u toku prelaznih režima u tri toplotno razmenjivačke podstanice sistema daljinskog grejanja "Zemun", i promene pritiska, temperature i protoka u toplani su date u (Stevanović i dr., 2006a). Poređenje rezultata numeričke simulacije sa izmerenim vrednostima pokazuje pouzdanost razvijenog modela. Izmereni pogonski uslovi su simulirani pomoću sopstvenog razvijenog numeričkog modela i kompjuterskog programa. Model je zasnovan na rešavanju nestacionarnog bilansa energije za jednodimenzionalno strujanje duž karakterističnog pravca, koji je određen kretanjem fluidnog delića u prostorno-vremenskom koordinatnom sistemu, što je opisano u poglavlju 4. ove disertacije. Postupak numeričkog rešavanja je eksplicitan i vremenski korak integracije je određen minimalnim vremenom potrebnim da fluidni delić pređe rastojanje od početnog položaja između dva susedna numerička čvora do sledećeg susednog čvora u pravcu strujanja unutar cevi u okviru cele mreže (Courant-ov kriterijum). Razvijene su bilansne jednačine koje omogućuju proračun graničnih uslova, kao što su spoj tri ili više cevi u čvoru, toplotno-razmenjivačke podstanice kod potrošača i u izvoru toplote. Hidraulički proračun je zasnovan na efikasnoj numeričkoj metodi bilansiranja promene pritiska po prstenovima cevne mreže (Stevanović i dr., 2006b; Stevanović i dr., 2007c). Poređenje sračunatih i izmerenih vrednosti prostiranja temperaturskih talasa od toplane do tri podstanice u različitim delovima mreže toplovoda i na različitim rastojanjima od izvora toplote pokazuje zadovoljavajuće slaganje. Za potpuno sagledavanje dinamičkog ponašanja sistema toplovoda u prelaznim režimima rada pomoću sprovođenja kompjuterskih simulacija i analiza, definisani su scenariji koji određuju prelazne procese u realnom odabranom sistemu, kao što su redosled i vremena uključivanja pumpi i pumpnih stanica, vremenska promena snage 79 toplotnog izvora i konzuma, dejstvo upravljačkih sistema i slično. Na osnovu ovako određenih scenarija definisani su početni i granični uslovi za proračune prelaznih režima rada. Prikazani su rezultati merenja i kompjuterske simulacije prelaznih temperaturskih procesa u sistemu daljinskog grejanja toplane ”Zemun“. Rezultati kompjuterske simulacije prostiranja talasa temperature su upoređeni sa izmerenim vrednostima i dobijeno je zadovoljavajuće slaganje. Slika 6.1 Poređenje rezultata kompjuterske simulacije prostiranja talasa temperature u sistemu daljinskog grejanja sa izmerenim vrednostima. Prelazni režim sa porastom snage toplane. Rezultati kompjuterske simulacije prostiranja temperaturskog talasa, usled prelaznog režima rada sa porastom snage toplane, su prikazani na Slici 6.1. Dobijene vrednosti pokazuju dobro slaganje sa izmerenim vrednostima. Vremenski intervali između položaja talasa na izlazu iz toplane i na ulazu u podstanice su određeni brzinom strujanja vode u razvodnim cevovodima. U toku prelaznog režima su radile tri paralelno povezane pumpe, čije su karakteristike date u (Stevanović i dr., 2006a). 80 Prostiranje temperaturskog talasa od toplotnog izvora do toplane je praćeno i simulirano, usled prelaznog režima rada sa smanjenjem snage izvora toplote, za spoljašnju temperaturu od 120C. Poređenje vrednosti temperature tople vode dobijene primenom razvijenog softvera i rezultata merenja u pojedinim postanicama prikazano je na Slici 6.2. Slika 6.2 Prostiranje temepraturskog talasa od toplane do potrošača. Poređenje numeričkih i izmerenih vrednosti. Prelazni režim sa smanjenjem snage toplane. 81 7. TERMOHIDRAULIČKI MODEL ZA PRORAČUN HIDRAULIČKOG UDARA IZAZVANOG KONDENZACIJOM PARE Razvijeni numerički postupak za predviđanje prostiranja talasa temperature je primenjen za simulaciju hidrauličkog udara izazvanog intenzivnom kondenzacijom pare. Ova pojava je određena dinamikom prostiranja razdelne površine između pare i pothlađene tečnosti na nižoj temperaturi i intenzivnom kondenzacijom pare na razdelnoj površini. Stoga je neophodno pouzdano predviđanje temperaturskih talasa na razdelnoj površini da bi se dobili pouzdani rezultati. 7.1 Matematički model jednodimenzionalnog nestacionarnog strujanja homogenog stišljivog fluida Za opisivanje hidrauličkog udara izazvanog kondenzacijom pare u cevovodu korišćen je homogeni model jednofaznog i dvofaznog stišljivog strujanja u kome su strujni i termodinamički parametri dobijeni usrednjavanjem parametara tečne i parne faze po jedinici mase. Ovaj model se može primeniti kako na jednofazno stujanje tečnosti ili pare u cevovodu, tako i na strujanje dvofazne mešavine. U slučaju strujanja dvofazne mešavine pretpostavljene su termička i strujna ravnoteža. S obzirom da su pri strujanju u cevovodu dve dimenzije strujnog prostora zanemarljive u odnosu na treću, strujanje je posmatrano kao jednodimenzionalno. Jednačine koje opisuju jednodimenzionalno nestacionarno strujanje homogenog stišljivog fluida u strujnom kanalu konstantnog poprečnog preseka su:  zakon održanja mase D 0, D u t x       (7.1)  zakon promene količine kretanja D 1 sign( ) sin 0, D 2 nt H fu uu p u f u g t x D t            (7.2)  zakon održanja energije 82   i i 2 D 1 D " sign( ) ' 0. D D 2 x c nt H x fu uh p u h f u u sign h h dx t t D t                 (7.3) U sistemu jednačina (7.1) do (7.3) zavisno promenljive veličine su brzina u, pritisk p i specifična entalpija h, a nezavisno promenljive veličine su prostorna koordinata x i vreme t. Termodinamičko stanje fluida i odgovarajuće gustine su određene na osnovu termodinamičkog stepena suvoće dvofazne mešavine ' . " ' t h h x h h    (7.4) Za 1tx  imamo stanje pare, za 0tx  stanje tečnosti, a za 0 1tx  stanje dvofazne mešavine. Gustina fluida se izračunava kao recipročna vrednost specifične zapremine dobijene primenom jednačine stanja za pothlađenu tečnost  ,tv p h , za pregrejanu paru  ,pv p h , a za tečnost i paru u stanju zasićenja, sledstveno  'v p i  ''v p , u obliku polinoma izvedenih na osnovu tablica za tečnost i paru (Wagner i Kretzschmar, 2007)   1/ ( , ), ako je 0 1/ '( ) ( ''( ) '( )) , ako je 0 1. 1 / ( , ), ako je 1 t t t t p t v p h x v p x v p v p x v p h x            (7.5) Treći članovi sa leve strane jednačina zakona promene količine kretanja i zakona održanja energije, jednačine (7.2) i (7.3), opisuju pad pritiska pri stacionarnom strujanju fluida kroz stujni kanal, gde su korišćene veličine DH hidraulički prečnik strujnog kanala i f Darcy-jev koeficijent trenja, a četvrti članovi sa leve strane jednačina (7.2) i (7.3) uzimaju u obzir uticaj nestacionarnog trenja (Prica, 2006) prema (Vardy i Brown, 2003) kao  10nt 1.1844 0.0567log Re 12.86 2 . Re f   (7.6) Poslednji član na levoj strani jednačine (7.3) predstavlja toplotnu snagu po jedinici mase fluida na razdelnoj površini između faza usled faznog prelaza, u ovom 83 slučaju kondenzacije. Toplota se prenosi sa pare ili dvofazne mešavine tečnosti i pare na stub tečnosti. Položaj razdelne površine u strujnom kanalu je označen sa ix , a sa  je obeleženo malo rastojanje od razdelne površine ka pari ili dvofaznoj mešavini i ka stubu tečnosti. U numeričkoj simulaciji  je jednako rastojanju između dva susedna čvora numeričke mreže. Brzina kondenzacije (kg/(m3s)) je određena kao ,c ic q a r   (7.7) pri čemu je cq toplotni fluks usled kondenzacije po jedinici razdelne površine (W/m2), ia specifična razdelna površina tečne i parne faze (m 2 /m 3 ), r latentna toplota kondenzacije (kJ/kg). Toplotni fluks usled kondenzacije po jedinici razdelne površine se računa kao  ,c c sat tq h T T  (7.8) gde je ch koeficijent prelaza toplote usled kondenzacije (W/(m 2 K)),  sat tT T razlika temperature pare u stanju zasićenja i temperature pothlađene tečnosti (K). Slika 7.1 Razdelna površina tečne i parne faze i prenos toplote usled kondenzacije. Funkcija  'sign h h određuje smer prostiranja toplote pri kondenzaciji, odnosno prenos toplote se odvija sa pare ili dvofazne mešavine kroz razdelnu površinu na stub tečnosti. 84 Glavni zadatak kod određivanja brzine kondenzacije pri hidrauličkom udaru izazvanom kondenzacijom pare je predviđanje specifične razdelne površine i koeficijenta prelaza toplote usled kondenzacije. Razdelna površina tečnosti i pare ima veoma nepravilan oblik (Slika 7.1). Za vreme hidrauličkog udara izazvanog kondenzacijom pare nestacionarni mlazevi tečnosti i kapi se odvajaju sa čela stuba tečnosti i uključuju se u zapreminu pare. Ova pojava uključivanja značajno povećava razdelnu površinu i brzinu kondenzacije. Uvrštavanjem izraza (7.8) u jednačinu (7.7) dobija se proizvod c ih a čija vrednost se računa kao zbir proizvoda koeficijenta prelaza toplote i specifične razdelne površine za kondenzaciju na čelu stuba tečnosti (ST) i za kondenzaciju uključenih kapi (K)     .c i c i c iST Kh a h a h a  (7.9) Koeficijent prelaza toplote pri kondenzaciji na čelu stuba tečnosti je aproksimiran prelazom toplote pri turbulentnom opstrujavanju cevi primenom poznate Dittus-Boelter (1930) korelacije   0.8 0.4, / 0.023Re Pr ,c ST t H t th D (7.10) pri čemu indeks t označava tečnu fazu, Re je Reynolds-ov broj, Pr je Prandtl-ov broj, a  toplotna provodnost. Specifična razdelna površina za čelo stuba tečnosti je određena na osnovu idealizovane pretpostavke da se razdelna površina poklapa sa površinom poprečnog preseka cevi Ai , 1 ,i ii ST i i A A a V A x x      (7.11) čime razdelna površina biva jednaka recipročnoj vrednosti prostornog koraka integracije, odnosno rastojanju između dva susedna čvora x . Koeficijent prelaza toplote za direktnu kondenzaciju na pothlađenim kapima se računa na osnovu teorijski predviđene i eksperimentalno potvrđene konstantne vrednosti Nusselt-ovog broja (Kronig i Brink, 1950; Kuznecov, 1989; Ghiaasiaan, 2008). , .c K K t h D Nu C    (7.12) 85 Vrednost konstante C=17,9 je usvojena prema (Kronig i Brink, 1950). Prečnik kapi se računa kao 4 4 4 4 4 4 10 , 10 , 10 5 10 , 5 10 , 5 10 K m Y D Y Y Y                  (7.13) gde je   2 ,kr p p t We Y u u     (7.14) a vrednost kritičnog Weber-ovog broja je 0,799crWe  (Saito i dr., 1978).  je površinski napon, a indeks p označava parametre parne faze. Veća neodređenost postoji kod predviđanja razdelne površine između parne faze i kapi uključenih sa čela stuba tečnosti. Pretpostavljeno je da uključivanje kapi i odgovarajuća specifična razdelna površina zavise od ubrzanja čela stuba tečnosti   33 , ,0 4 10 ,i K K ia a Du Dt    (7.15) gde je 2 3 ,0 0 40 m /mKa   . Parametri na desnoj strani jednačine (7.15) su određeni poređenjem rezultata sprovedenih numeričkih simulacija hidrauličkog udara sa eksperimentalnim merenjima dostupnim u literaturi (Zaltsgendler i dr., 1996; Liu i dr., 1996; Yeung i dr. 1993). Prvi član sa leve strane jednačine (7.15) uzima u obzir uticaj hidrodinamičkih stanja prostiranja stuba tečnosti ka zapremini pare, koje je uglavnom određeno početnim stanjem nestacionarnog procesa odmah posle otvaranja ventila ili uklanjanja neke druge pregrade između zapremina tečnosti i pare. Otvaranje ventila koji u početnom trenutku razdvaja tečnost i paru uzrokuje stohastičke vrednosti oblika razdelne površine, a na koje utiče formiranje mlazeva tečnosti i uključivanja kapi u zapreminu pare. Zbog ove stohastičke prirode formiranja razdelne površine, vrednost parametra ,0Ka varira unutar propisanog opsega od nule do maksimalne vrednosti 2 340 m /m . Vrednost nula znači da nema uključenih kapi koje doprinose povećanju razdelne površine za direktnu kondenzaciju. Drugi član u jednačini (7.15) uzima u obzir uticaj ubrzanja stuba tečnosti na razbijanje čela stuba tečnosti i uključivanje kapi sa čela stuba u zapreminu pare. 86 7.1.1 Formiranje sistema kvazilinearnih parcijalnih diferencijalnih jednačina hiperboličkog tipa Jednačina stanja, napisana u obliku  ,p h  , je diferencijabilna po vremenu i po prostornoj koordinati. , ph p h t p t h t                      (7.16) . ph p h x p x h x                      (7.17) Transformacijom jednačine zakona održanja mase (7.1), tako što se materijalni izvod gustine razvije prema definiciji (7.18) D , D u t t x         (7.18) i u dobijenu jednačinu uvrste jednačine (7.16) i (7.17), a zatim ponovo primeni definicija materijalnog izvoda na pritisak i entalpiju dobija se D D 0. D Dph p h u p t h t x                   (7.19) Smenom jednačine (7.19) u jednačinu zakona održanja energije (7.3), pri čemu se brzina zvuka u jednofaznom fluidu, a time i u homogenom modelu dva fluida, prema definiciji   1 2 s c p    u funkciji pritiska, gustine i entalpije određuje kao 1 , 1 ph c p h                (7.20) dobija se   i i 2 2 2D "sign( ) ' 0. D 2 x c nt p H x fu up u u h c c f u u sign h h dx t x h D t                              (7.21) U jednačinama (7.21), (7.2) i (7.3) razvijanjem materijalnog izvoda i prebacivanjem na levu stranu jednačine svih onih veličina koje se diferenciraju po 87 vremenu i koordinati dobijaju se tri jednačine koje formiraju sistem kvazilinearnih parcijalnih diferencijalnih jednačina hiperboličkog tipa. 2 , p p u u c X t x x           (7.22) 1 , u u p u Y t x x          (7.23) 1 , h h p p u u Z t x t x                (7.24) gde je   i i 2 2 "sign( ) sin ' , 2 x c nt p H x fu u u h X c f u u gu sign h h dx h D t                            (7.25)  sign sin , 2 nt H fu u u Y f u g D t        (7.26)   i i 2 " sign( ) ' . 2 x c nt H x fu u u h Z f u u sign h h dx D t              (7.27) Da bi se ovaj sistem jednačina mogao rešiti potrebno je zadati odgovarajuće početne i granične uslove. Početni uslovi su definisani strujnotermičkim parametrima fluida u početnom vremenskom trenutku, pre dejstva poremećaja. Granični uslovi su definisani na osnovu stanja fluida na početku i na kraju deonica, kao i na granicama posmatranog sistema. Analitičko rešenje ovog sistema nije moguće dobiti pa se nekom numeričkom metodom određuje partikularni integral. 7.2 Numeričko rešavanje jednodimenzionalnog nestacionarnog strujanja homogenog stišljivog fluida primenom metode karakteristika Metoda karakteristika se koristi za numeričko rešavanje sistema parcijalnih diferencijalnih jednačina hiperboličkog tipa. Primena ove metode na rešavanje višefaznih strujanja je ograničena na jednodimenzionalna strujanja. Metoda karakteristika daje najtačnije rezultate i zbog toga je veoma pogodna za referentne 88 proračune i benchmark testove (testove na osnovu kojih se određuju mogućnosti kompjuterskih programa da simuliraju neku pojavu) (Wulff, 1987). Visoka tačnost metode karakteristika potiče iz činjenice da ona prevodi parcijalne diferencijalne jednačine u obične diferencijalne jednačine ili u integralne jednačine bez uticaja numeričke difuzije. Metoda karakteristika je jedina metoda koja tačno prati prostiranje diskontinuiteta primenom izvoda prvog reda (Wulff, 1987). Za takve diskontinuitete kao karakteristične koordinate koriste se Lagrange-ove koordinate. Primenom metode karakteristika se rešava sistem od tri kvazilinearne parcijalne diferencijalne jednačine (7.22) do (7.24), sa tri zavisno promenljive veličine (pritisak, brzina i entalpija) i dvema nezavisno promenljivima (vreme i prostorna koordinata). Postupak se sastoji u sledećem; traže se linearne kombinacije ove tri parcijalne diferencijalne jednačine, takve da sadrže izvode nepoznatih funkcija u samo jednom pravcu (ovde u pravcu prostorne koordinate - x). Za tri parcijalne diferencijalne jednačine postoje tri karakteristična pravca duž kojih važe izvedene jednačine. Na ovaj način se metodom karakteristika sistem kvazilinearnih parcijalnih diferencijalnih jednačina (7.22) do (7.27) prevodi u sistem diferencijalnih jednačina sa totalnim diferencijalima, pri čemu se određuju familije krivih u prostorno-vremenskoj ravni duž kojih važe izvedene transformacije. Totalni diferencijali se zatim zamenjuju konačnim razlikama čime se dobijaju tri diferencne jednačine. Rešavanjem ovih jednačina po zavisno promenljivim veličinama dobijaju se vrednosti strujno-termičkih parametara fluida duž strujnog prostora tokom vremena. Familije krivih u prostorno-vremenskoj ravni fizički predstavljaju prostiranje talasa pritiska i entalpije fluidnih delića u svakoj tački strujnog polja. Korak integracije je određen pomoću Courant-ovog kriterijuma. Courant-ov kriterijum povezuje prostorni i vremenski korak integracije. Na osnovu njega se određuje diferencna mreža u prostorno-vremenskom koordinatnom sistemu i za svaki čvor se rešavaju po tri diferencne jednačine. 89 7.2.1 Formiranje sistema običnih diferencijalnih jednačina sa totalnim diferencijalima - formiranje jednačina karakteristika Množenjem jednačine (7.22) koeficijentom 1, jednačine (7.23) koeficijentom 2, jednačine (7.24) koeficijentom 3 i sabiranjem nastalih jednačina uz odgovarajuće grupisanje članova dobija se   3 2 3 1 2 3 1 2 1 2 3 1 2 3 . p u h u p u t t t x u h c u u X Y Z x x                                                    (7.28) Zavisno promenljive veličine, uopšteno označene nekom funkcijom  , ,f p u h su diferencijabilne po vremenu i po prostornoj koordinati . x t f f df dt dx t x                (7.29) Smenom jednačine (7.29) u (7.28) dobija se     3 2 3 2 3 1 1 1 2 2 2 1 2 1 2 3 3 3 1 2 3 . u dt p u dp u u dx t dx dt u du c u c u dx t dx dt h dh u u dx t dx X Y Z                                                                             (7.30) Koeficijenti 1, 2 i 3 se određuju iz uslova da su izrazi u jednačini (7.30) uz parcijalni izvod zavisno promenljivih veličina, p, u ili h, po vremenu jednaki nuli. Grupisanjem po koeficijentima 1, 2 i 3 dobija se sistem od tri linearne homogene jednačine 1 2 3 2 1 2 3 1 1 1 0, 1 0, 1 0. dt dt u dt u dx dx dx dt dt c u dx dx dt u dx                                          (7.31) 90 Rešenja ovog sistema biće netrivijalno ako i samo ako je determinanta sistema jednaka nuli 2 1 1 1 1 0 0. 0 0 1 dt dt u dt u dx dx dx dt dt c u dx dx dt u dx                  (7.32) Rešenja determinante po dt dx predstavljaju karakteristične pravce, tj. jednačine karakteristika 1 1 1 , , . dt dx u u c u c         (7.33) Zamenom 1dt dx u  u jednačine (7.31) dobija se 1 2 3 3 0, 0,       (7.34) odnosno za 1dt dx u c   u jednačine (7.31) 1 3 2 2 1 , 0, c        (7.35) i konačno za 1dt dx u c   u jednačine (7.31) dobija se 1 3 2 2 1 , 0. c         (7.36) Zamenjujući jednačine (7.34), (7.35) i (7.36) u jednačinu (7.30) gube se parcijalni diferencijali zavisno promenljivih i dobijaju se obične diferencijalne jednačine sa njihovim totalnim diferencijalima. Ove jednačine se nazivaju jednačinama karakteristika. Za 1dt dx u  , CP karakteristika . u dp dh u Z dx dx    (7.37) Za 1dt dx u c   , C + karakteristika   1 . u dp du X c u Y c dx dx c             (7.38) Za 1dt dx u c   , C - karakteristika   1 . u dp du X c u Y c dx dx c                (7.39) 91 7.2.2 Formiranje sistema diferencnih jednačina Smenom jednačina (7.33) u jednačine (7.37), (7.38) i (7.39) dobijaju se za 1dt dx u  , CP karakteristika 1 ,dh dp Zdt    (7.40) za 1dt dx u c   , C + karakteristika   ,dp cdu X cY dt    (7.41) za 1dt dx u c   , C - karakteristika   .dp cdu X cY dt    (7.42) Jednačina (7.40) opisuje strujanje fluidnog delića, koje je povezano sa prostiranjem fronta entalpije, a jednačine (7.41) i (7.42) prostiranje (propagaciju) talasa pritiska. Zakoni održanja mase i promene količine kretanja rešavaju se duž C+ i C- karakteristika, a zakon održanja energije se rešava duž CP karakteristike. U sistemu jednačina (7.40) do (7.42) se totalni diferencijali aproksimiraju konačnim razlikama. Konačne razlike se uzimaju duž karakterističnih pravih linija. Na taj način se dobija sistem diferencnih jednačina. Koeficijenti u jednačinama (7.40) do (7.42) se smatraju konstantnim u toku vremenskog koraka integracije, a njihove vrednosti se dobijaju linearnom interpolacijom rezultata prethodnog računskog koraka. 7.2.3 Rešavanje sistema diferencnih jednačina primenom metode karakteristika Na Slici 7.2 je predstavljen vremenski trenutak t, bilo da je to početni ili prethodni vremenski trenutak, i trenutak t+∆t kao naredni vremenski trenutak. Tačke A, B i C predstavljaju proizvoljno izabrana tri uzastopna čvora, unutar posmatranog strujnog prostora, u kojima se računaju vrednosti zavisno promenljivih veličina u trenutku t. Tačka D je mesto do kog stigne poremećaj za vreme ∆t iz tačaka R i L. U tački D se seku sve tri karakteristične linije: C+ koja prolazi kroz tačku R, C- koja prolazi kroz tačku L i CP koja prolazi kroz tačku P, koje polaze iz različitih tačaka (negde između čvorova A, B i C) u prethodnom vremenskom trenutku t. Posledično, tačka D predstavlja stanje u strujnom prostoru koje nastaje u narednom vremenskom trenutku t+∆t kao posledica prostiranja poremećaja koji je nastao u vremenskom 92 trenutku t. Vremenski korak se određuje iz Courant-ovog kriterijuma stabilnosti prema kome poremećaj koji krene iz tačke R, krećući se brzinom u+c, i poremećaj koji krene iz tačke L, krećući se brzinom u-c, ne smeju da pređu tačku D jer bi to izazvalo nestabilnost rešenja. Rastojanja duž x-ose između čvorova A i B, i B i C su međusobno jednaka i konstantna tokom vremena, tj. AB BC . Zavisno promenljive veličine su poznate u svim čvorovima u trenutku t i potrebno je izračunati njihove vrednosti u trenutku t+∆t. Slika 7.2 Prostorno-vremenska ravan. Aproksimirajući totalne diferencijale iz jednačina (7.40), (7.41) i (7.42) konačnim razlikama duž karakterističnih pravaca dobija se sistem za 1 R R t x u c     , C + karakteristika ( ) ( ) ( ) ,R RD R D R R R R R c c p p u u X Y t v v       (7.43) za 1 L L t x u c     , C - karakteristika ( ) ( ) ( ) ,L LD L D L L L L L c c p p u u X Y t v v       (7.44) za 1 P t x u    , C P karakteristika ( ) ( ) .D P P D P Ph h v p p Z t     (7.45) Rešavanjem algebarskih jednačina (7.43) i (7.44) po pD i uD dobija se ,Dp        (7.46) ,Du        (7.47)  t x x A R B D L CP C P C - C + t+ t t t 93 gde su ,R R c v   (7.48) ( ) ,R R R Rp u X Y t       (7.49) ,L L c v   (7.50) ( ) .L L L Lp u X Y t       (7.51) Uvrštavanjem pD iz jednačine (7.46) u jednačinu (7.45) dobija se hD ( ) .D P P D P Ph h v p p Z t     (7.52) Jednačinama (7.46), (7.47) i (7.52) su određene vrednosti zavisno promenljivih veličina u trenutku t+∆t u čvoru D samo preko vrednosti zavisno promenljivih veličina u trenutku t u tački P. Da bi se odredile ove veličine potrebno je prethodno odrediti vrednosti zavisno promenljivih u tačkama R i L linearnom interpolacijom vrednosti promenljivih u čvorovima A i B, i B i C. Tačka R ,B R B R B A B A x x u u x x u u      (7.53) ,B R B R B A B A x x c c x x c c      (7.54) ,B R B R B A B A x x p p x x p p      (7.55) 1 . B R R R t x x u c     (7.56) Rešavajući jednačine (7.53) do (7.56) po uR, cR i pR dobija se  1 , 1 B B R b u ac u a b      (7.57)  1 , 1 B B R a c bu c a b      (7.58)   , 1 B B R B B A u c t p p p p a b x         (7.59) kao i 94   , 1 B B R B B A u c t v v v v a b x         (7.60) gde je    , .B A B A t t a u u b c c x x         (7.61) Tačka L ,B L B L B C B C x x u u x x u u      (7.62) ,B L B L B C B C x x c c x x c c      (7.63) ,B L B L B C B C x x p p x x p p      (7.64) 1 . B L L L t x x u c     (7.65) Rešavajući jednačine (7.62) do (7.65) po uL, cL i pL dobija se  1 , 1 B B L d u ec u e d      (7.66)  1 , 1 B B L e c du c e d      (7.67)   , 1 B B L B B C u c t p p p p e d x         (7.68) kao i   , 1 B B L B B C u c t v v v v e d x         (7.69) gde je    , .B C B C t t e u u d c c x x         (7.70) Tačka P Za 1 0 ; B dt dx u   (7.71) 95 ,B P B P B A B A x x u u x x u u      (7.72) ,B P B P B A B A x x h h x x h h      (7.73) ,B P B P B A B A x x p p x x p p      (7.74) 1 . B P P t x x u    (7.75) Rešavajući jednačine (7.71) do (7.75) dobija se , 1 B P u u a   (7.76) ( ) , 1 B P B B A u t h h h h a x       (7.77)   , 1 B P B B A u t p p p p a x       (7.78) kao i   . 1 B P B B A u t v v v v a x       (7.79) Za 1 0; B dt dx u   (7.80) ,B P B P B C B C x x u u x x u u      (7.81) ,B P B P B C B C x x h h x x h h      (7.82) ,B P B P B C B C x x p p x x p p      (7.83) 1 . B P P t x x u    (7.84) Rešavajući jednačine (7.71) do (7.75) dobija se , 1 B P u u c   (7.85) ( ) , 1 B P B B C u t h h h h c x       (7.86) 96   , 1 B P B B C u t p p p p c x       (7.87) kao i   . 1 B P B B A u t v v v v c x       (7.88) Kada se odrede sve zavisno promenljive veličine (p, u, h), kao i specifične zapremine, u čvorovima R, L i P pristupa se određivanju pritiska, brzine i entalpije u čvoru D prema jednačinama (7.46), (7.47) i (7.52). Specifična zapremina fluida u čvoru D određuje se iz jednačine stanja, a lokalna brzina zvuka u fluidu iz korelacije za brzinu zvuka. U ovoj metodi je prostorni korak, tj. rastojanje između čvorova, konstantan. Vremenski korak integracije se određuje Courant-ovim kriterijumom stabilnosti rešenja min , 1,2, ... , , J J x t J n c u         (7.89) pri čemu je minimalni vremenski korak, za određenu vrednost prostornog koraka, definisan maksimalnom vrednošću zbira brzine zvuka i apsolutne vrednosti brzine strujanja fluida. Treba napomenuti da se traži maksimalna vrednost zbira brzine zvuka i apsolutne vrednosti brzine od svih n čvorova u posmatranom sistemu. Kako bi se umanjila numerička difuzija pri prostiranju talasa entalpije, a koja predstavlja oštar (jasan) diskontinuitet između entalpija hladne tečnosti i vrele pare, početna vrednost entalpije u čvoru P je određena primenom Lagrange-ovog interpolacionog polinoma trećeg reda. 7.3 Homogeni model dvofaznog strujanja Homogeni model posmatra strujanje dva fluida kao strujanje jednog fluida čiji su strujni i termodinamički parametri dobijeni usrednjavanjem parametara tečne i parne faze po jedinici mase dvofazne mešavine. 97 Homogeni model omogućava jednostavnije posmatranje strujanja dva fluida, vode i vodene pare, mešavine vode i pare ili bilo koje tečne i gasne faze, u odnosu na model dva fluida, koji posmatra svaku fazu odvojeno. Iako je primena homogenog modela ograničena zbog pretpostavke koju u sebi sadrži, da su brzine dveju faza jednake, on naročito olakšava proračun strujanja dve faze pri kojem dolazi do promene faza. Specifične zapremine pothlađene vode, zasićene vode, zasićene pare i pregrejane pare se određuju polinomima oblika: , 0 0 , K N n k n k k n v C p h    (7.90) a specifične entalpije zasićene vode i pare polinomima oblika: 0 . N n n n h C p   (7.91) Jednačine stanja su aproksimirane polinomima koji su preuzeti iz programa PIPES (Choi, 1983), a određeni su metodom najmanjih kvadrata vrednosti iz parnih tablica standarda ASME. Navedene korelacije su date u tekstu odgovarajućih podprograma. 7.4 Određivanje brzine zvuka Brzina zvuka u jednofaznom fluidu se prema definiciji određuje kao 2 . s p c         (7.92) Korišćenjem veze između gustine i specifične zapremine, ( 1 v   ) brzina zvuka se može izraziti u funkciji pritiska i specifične zapremine kao 2 2 . s p c v v         (7.93) 98 Za potrebe razvijenog programa brzinu zvuka jednofaznog fluida je potrebno izraziti u zavisnosti od pritiska i entalpije. Diferenciranjem jednačine stanja ( vv(p,h) ) dobija se . ph v v dv dp dh p h              (7.94) Pošto definicija brzine zvuka glasi da je kvadrat brzine zvuka jednak parcijalnom izvodu pritiska po gustini pri konstantnoj entropiji, prema drugom zakonu termodinamike ( q ds T   ) sledi q  0. Daljim uvrštavanjem u jednačinu koja opisuje prvi zakon termodinamike ( q du pdv   ), uz korišćenje definicije entalpije ( ( )dh du d pv  ) sledi: .dh vdp (7.95) Zamenom jednačine (7.95) u jednačinu (7.94) dobija se , ps h v v v v p p h                      (7.96) čijim se uvrštavanjem u jednačinu (7.93) konačno dobija jednačina za određivanje brzine zvuka u jednofaznom fluidu u funkciji pritiska i entalpije i specifične zapremine ili gustine 2 1 1 ili , 1 1 1 p ph h c c v v v p v h p h                                  (7.97) Da bi se uštedelo vreme pri računanju moguće je koristiti relacije za brzinu zvuka u pothlađenoj i zasićenoj vodi   ,t tc c h (7.98) a u pregrejanoj i zasićenoj pari  , .p pc c p h (7.99) U program, kojim se računa brzina prostiranja talasa pritiska u dvofaznoj mešavini vode i vodene pare, ugrađena je i korelacija dobijena teorijskim modelom 99 kojim se pretpostavlja da ne dolazi do faznih prelaza, (Choi, 1983, Tong i Weisman, 1979)   2 2 1 , 1 '' 'p t c c c        (7.100) gde je  zapreminski udeo parne faze u dvofaznoj mešavini i određuje se kao: '' , V V   a  gustina dvofazne mešavine i određuje se kao:  '' 1 ' .      Brzina zvuka u dvofaznoj mešavini se u funkciji entalpije i pritiska određuje prema korelaciji za čepasti tok kao   , 1 p t t p c c c c c     (7.101) gde se zapreminski udeo parne faze u dvofaznoj mešavini računa prema 1 , 1 1 t t t p x v x v     (7.102) pri čemu je tx termodinamčki stepen suvoće. 7.5 Uticaj nestacionarnog trenja U opštem slučaju uticaj trenja možemo posmatrati kao uticaj njegovih komponenata, stacionarnog i nestacionarnog trenja. Ove dve komponente trenja imaju različit uticaj na strujanje stoga imaju i različite matematičke formulacije pa ih je potrebno odvojeno računati. U praktičnoj primeni u inženjerskoj praksi potrebno je uglavnom dati ukupan uticaj trenja na neko strujanje, dok je analiza njegovoh pojedinih komponenata manje značajna. Prilikom numeričke simulacije hidrauličkog udara primećeno je da uzimanje u obzir samo pojave stacionarnog trenja nije dovoljno za postizanje odgovarajuće tačnosti pri modeliranju hidrauličkog udara. Pri hidrauličkom udaru se javljaju veliki skokovi 100 (padovi) pritisaka koji rezultuju skokovitim promenama brzine, pa se u matematički model koji opisuje ovakva strujanja mora uvesti nestacionarno trenje. Tačnost simulacije prostiranja fronta talasa pritiska zavisi od mogućnosti primenjenog modela da modelira pojave kao što su kavitacija, međudelovanje fluida i strukture i trenje o zidove cevi, koje utiče na razvijanje oblika talasnog fronta, na njegovu amplitudu i brzinu prostiranja (Vardy i Brown, 2003). Trenje fluida o zidove cevi se najčešće u nestacionarnom strujanju modelira primenom kvazi stacionarnog modela, u kome se pretpostavlja da je otpor nestacionarnom strujanju jednak otporu koji bi postojao u stacionarnom strujanju, pri čemu to strujanje ima istu srednju vrednost brzine kao i posmatrano nestacionarno strujanje. Ova pretpostavka daje numeričke rezultate koji se dobro slažu sa eksperimentalnim vrednostima u slučajevima kada su ubrzanja mala i kada se posmatra strujanje daleko od talasnog fronta. U blizini fronta uticaj nestacionarnog trenja na trenje o zidove cevi je značajan, zato što od njega zavisi brzina prostiranja talasnog fronta i razvoj njegovog oblika. Modeli za određivanje nestacionarnog trenja se mogu podeliti u tri grupe. Prva grupa se zasniva na trenutnom ubrzanju, druga koristi brzine/ubrzanja iz prethodnog vremenskog perioda, a treća se zasniva na nepovratnoj termodinamici (Vardy i Brown, 2003). Modeli zasnovani na trenutnom ubrzanju su najpopularniji za upotrebu u numeričkoj simulaciji. Najčešće korišćen model je Brunone-ov model. Nestacionarno trenje je proporcionalno koeficijentu čija se vrednost određuje empirijski. Druga grupa modela se zasniva na Zielke-ovom pristupu za laminarna strujanja. Vardy i Brown su primenili Zielke-ov model na idealizovani model turbulentnog strujanja u kome je pretpostavljeno da se turbulentna viskoznost linearno menja u blizini zida cevi. Ovo je omogućilo da se utvrdi analitička veza napona smicanja na zidu koji prati iznenadnu promenu brzine. Dvodimenzionalni modeli za određivanje nestacionarnog trenja su u mogućnosti da daju više detalja i bolje predviđanje od jednodimenzionalnih modela ali nisu pogodni za uopštene analize. Vardy i Brown (2003) su razvili model nestacionarnog turbulentnog trenja u jednodimenzionalnim glatkim cevima. Primena modela je ograničena Reynolds-ovim 101 brojem 10 8 , ali se granica primene može proširiti po potrebi. Oni su primenili idealizovanu raspodelu viskoznosti u poprečnom preseku. Poprečni presek je podeljen na dva regiona: jezgro i anularni prostor. U anularnom prostoru je pretpostavljeno da se viskoznost menja linearno od laminarne vrednosti na zidu do maksimalne viskoznosti jezgra na granici ovih regiona, koja je dalje uniformna u celom regionu jezgra. Numerički rezultati su upoređeni sa eksperimentalnim merenjima za stacionarno strujanje. Pretpostavljeno je da je viskoznost u oba regiona konstantna u toku kratkog perioda nestacionarnog strujanja, što je potvrđeno eksperimentalnim i teorijskim studijama. Ova pretpostavka se ne može primeniti za duže periode. Primenom ove metode određena je težinska funkcija turbulentnog strujanja koja zavisi od Reynolds- ovog broja. Dobijeni rezultati su analitički definisani i važe za ceo opseg vremena. Metoda je primenjena za određivanje koeficijenta koji opisuje nestacionarno trenje u modelu u kom je pretpostavljeno da je napon smicanja na zidu proporcionalan trenutnoj vrednosti ubrzanja. Ovaj pristup je primenjen za izračunavanje nestacionarnog trenja u ovoj disertaciji, i pokazao se kao dovoljno dobar za prikazivanje uticaja ove pojave na promene nastale u strujanju usled hidrauličkog udara. Prema Literaturi (Vardy i Brown, 2003) komponenta nestacionarnog trenja se određuje na osnovu izraza * 0 2 , T lam wu U W dt R t      (7.103) gde je wu nestacionarno trenje (N/m 2 ), lam kinematska viskoznost fluida u laminarnom toku (m 2 /s), R poluprečnik cevi (m), T vreme računato od početka nestacionarne pojave (s), t * vreme mereno unazad od trenutka T, ( *t T t  ) (s), /U t  lokalno ubrzanje fluida u trenutku t (m/s 2 ), W težinska funkcija koja zavisi od t* (-). Množilac 2 sa desne strane jednačine (7.103) je uveden zbog kompatibilnosti sa težinskom funkcijom koju je razvio Zielke za nestacionarno laminarno strujanje. Laplace-ovom transformacijom jednačine (7.103) dobija se 102 ' 2 ' ' ,lamwu U W R t          (7.104) gde je 'W funkcija od odnosa /w lam  . Vrednost ove funkcije zavisi od vremena nastanka poremećaja u strujnom toku. Za izrazito male vrednosti vremena ova funkcija parametarski zavisi od Reynolds-ovog broja dok se za veće vrednosti ovaj uticaj Reynolds-ovog broja ne oseća. Radi lakše upotrebe pri numeričkim proračunima Zielke je dao sledeću aproksimaciju težinske funkcije 'W ' .a A W s B   (7.105) Parametri A i B su određeni tako da obezbede tačnu aproksimaciju asimptota funkcije pri velikim i malim vrednostima koordinate s (Laplace-ova transformacija vremena), i računaju se kao 2 w lam A R    i 2 ' 0 , A B W        (7.106) gde je 0 'W asimptotska vrednost 'W pri malim vrednostima s, a zavisi od Rejnolds- ovog broja. Uvrštavanjem izraza za A i B u jednačinu (7.105) i primenom inverzne Laplace- ove transformacije dobija se * , 2 w lamC aW e            (7.107) pri čemu se uzima da je w lam   =1, i gde je * 2 lamt R    bezdimenziono vreme, a * 2 lamC BR   koeficijent kojim se uzima u obzir opadanje intenziteta trenja. Radi lakše praktične primene izraz za 0 'W (odnosno B) nije izveden, već je sprovođenjem velikog broja proračuna za veliki opseg Reynolds-ovih brojeva ustanovljena relacija koja povezuje C * i Re i važi za 2000 < Re < 108 , i ima oblik 103 * 12,86 , Re C   gde je 10 0,0567 15,29 log . Re         (7.108) Pri malim vrednostima Reynolds-ovih brojeva C * se jako menja dok pri velikim vrednostima Re koeficijent C * manje menja svoju vrednost. Pad pritiska duž strujnog kanala, pri proračunu hidrauličkog udara izazvanog kondenzacijom pare, određuje se prema   ,nt t p u f sign u x t          (7.109) gde je *2 ,ntf C odnosno  10nt 1.1844 0.0567log Re 12.86 2 , Re f   (7.110) a u brzina kretanja razdelne površine tečne i parne faze. 7.6 Matematičko definisanje graničnih uslova Pri izračunavanju parametara fluida na početku ili na kraju cevi neke od karakteristika potrebnih za izračunavanje nije moguće nacrtati jer one izlaze iz posmatranog strujnog prostora. Zbog toga je na tim mestima potrebno napisati dopunske hidrauličke uslove koji ujedno predstavljaju granične uslove, koji su pored početnih uslova potrebni za rešavanje sistema kvazilinearnih parcijalnih diferencijalnih jednačina hiperboličkog tipa, jednačine (7.25) do (7.27). U programu se pravolinijske deonice, odnosno cevi, obeležavaju opadajućim nizom prirodnih brojeva u pretpostavljenom smeru proticanja fluida, a čvorovi unutar cevi rastućim nizom. Dakle, granične uslove je potrebno definisati za prvi i poslednji čvor deonice (cevi). Granične uslove, prema mestu na kome se nalaze, delimo na granične uslove na početku, granične uslove na kraju cevi i granični uslov u spoju dve cevi. U posmatranim slučajevima na početku i na kraju cevi se javljaju: zatvoreni kraj cevi ili rezervoar. 104 7.6.1 Granični uslov na početku cevi Na Slici 7.3 je prikazan početak cevi u prostorno-vremenskoj ravni sa nacrtanom (levom) C - karakteristikom. Ovo je jedina karakteristika koja se može nacrtati u ovom slučaju. Jednačina koja važi duž C- karakteristike je ( ) ( ) ( ) ,L LD L D L L L L L c c p p u u X Y t v v       za nagib 1 . L L t x u c     (7.111) Veličine uL, cL, pL,vL, XL i YL se određuju na već opisani način redom jednačinama: (7.66), (7.67), (7.68), (7.69), (7.25) i (7.26). Slika 7.3 Početak cevi u prostorno-vremenskoj ravni sa nacrtanom C- karakterisitkom. 7.6.1.1 Rezervoar na početku cevi Kada je cev svojim prvim čvorom spojena sa rezervoarom i ako su zaustavni pritisak i entalpija fluida u rezervoaru zadate funkcije vremena; pritisak i entalpija fluida u tački D (Slika 7.3) određuju se jednačinama   2 0, , 2 B D rez B u p p t v   (7.112)   2 0, , 2 B D rez u h h t  (7.113) x  t B D L C - t t+ t t 1 2 3čvor 105 pri čemu su ulazni gubici zanemareni. Kako bi se smanjio broj iteracija, a time skratilo vreme računanja, dinamički pritisak se umesto preko brzine u trenutku t+∆t (brzine u tački D), računa koristeći poznatu brzinu fluida u čvoru B koja je rezultat prethodnog računskog koraka, brzina u trenutku t. Brzina fluida u tački D određuje se iz jednačine (7.111) kao   ,L LD L D L L L L L v v u u p p Y X t c c            (7.114) a specifična zapremina iz jednačine stanja  ,D D Dv v p h . 7.6.2 Granični uslov na kraju cevi Na Slici 7.4 je prikazan kraj cevi u prostorno-vremenskoj ravni sa nacrtanim C + i C P karakterisitkama. U ovom slučaju nije moguće nacrtati samo C- karakteristiku. Jednačina koja važi duž C+ karakteristike je ( ) ( ) ( ) .R RD R D R R R R R c c p p u u X Y t v v       za nagib 1 . R R t x u c     (7.115) Veličine uR, cR, pR,vR, XR i YR se određuju na već opisani način redom jednačinama: (7.57), (7.58), (7.59), (7.60), (7.25) i (7.26). Slika 7.4 Kraj cevi u prostorno-vremenskoj ravni sa nacrtanim C + i C P karakterisitkama.  t x A R B D P C P C + t+ t t t N-2 N-1 Nčvor 106 7.6.2.1 Zatvoreni kraj cevi na kraju cevi U neposrednoj blizini zatvorenog kraja cevi fluid miruje, njegova brzina jednaka je nuli. Uslov koji mora biti zadovoljen je 0.B Du u  (7.116) Uvrštavanjem uslova (7.116) u jednačinu (7.115) dobija se ( ) .R RD R R R R R R c c p p u X Y t v v      (7.117) Entalpija fluida u tački D određuje se iz jednačine (7.52) pri čemu se uzima da se tačka P poklapa sa čvorom B čime je: 0 , , ,P P B P B P Bu p p h h v v    i sledi ( ) .D B B D B Bh h v p p Z t     (7.118) Specifična zapremina fluida u tački D određuje se iz jednačine stanja  ,D D Dv v p h . 7.6.2.2 Rezervoar na kraju cevi Pritisak i entalpija fluida u tački D, kada je cev svojim krajem spojena sa rezervoarom, se određuju kao i u slučaju kada je rezervoar spojen sa početkom cevi jednačinama (7.112) i (7.113). Dinamički pritisak se takođe, umesto preko brzine u tački D, računa koristeći poznatu brzinu fluida u čvoru B koja je rezultat prethodnog računskog koraka. Brzina fluida u tački D određuje se iz jednačine (7.115)   ,R RD R D R R R R R v v u u p p Y X t c c            (7.119) a specifična zapremina iz jednačine stanja  ,D D Dv v p h . 7.6.3 Granični uslov u spoju dve cevi istog poprečnog preseka Pošto su u svim eksperimentima, koji su numerički simulirani, cevi u spoju imale međusobno jednake prečnike navedene relacije važe samo za taj slučaj. Kada su 107 poprečni preseci cevi koje se nadovezuju jedna na drugu različiti, relacije koje definišu strujanje su znatno komplikovanije zbog pojave lokalnih gubitaka energije strujanja fluida na mestu spoja, pa se pritisci i brzine strujanja na mestu prelaza iz jedne cevi u drugu međusobno razlikuju. Slika 7.5 Spoj dve cevi istog poprečnog preseka u prostorno-vremenskoj ravni sa nacrtanim C + , C P i C - karakterisitkama. U spoju dve cevi, za jednu cev se mogu nacrtati samo neke karakteristike, pri čemu one karakteristike koje se ne mogu nacrtati za tu cev, mogu se nacrtati za cev sa kojom je ona u spoju. Na Slici 7.5 nalazi se spoj dve cevi prikazan u prostorno-vremenskoj ravni. Obeležavanje cevi je takvo, da broj cevi opada u smeru pretpostavljenog strujanja fluida tako da je J+1. cev svojim N-tim čvorom spojena sa prvim čvorom J-te cevi. Jednačina C+ karakteristike napisana za J+1. cev (7.120) i jednačina C- karekteristike napisane za J-tu cev (7.121) , 1 , 1( ) ( ) ( ) , R R D J R D J R R R R R c c p p u u X Y t v v        (7.120) , ,( ) ( ) ( ) , L L D J L D J L L L L L c c p p u u X Y t v v       (7.121)  t x x A R BJ+1 DJ+1 L CP C P C - C + t+ t t t N-1 N 2čvor BJ DJ cev (J+1) cev (J)1 u 108 uz korišćenje izraza (7.48) do (7.51) dobijaju oblik , 1 1 , 1 1,D J J D J Jp u      (7.122) , , .D J J D J Jp u   (7.123) Pošto u slučaju spoja dve cevi istog poprečnog preseka nema gubitaka strujne energije fluida u spoju, sledi da je , 1 ,D J D Jp p  i , 1 , .D J D Ju u  Iz jednačina (7.122) i (7.123) dobija se 1 1, , 1 1 ,J J J JD J D J J J p p              (7.124) 1, , 1 1 .J JD J D J J J u u            (7.125) Entalpija fluida u tački D se određuje iz jednačine (7.52), tj. jednačine CP karakteristike, uz pretpostavku da je entalpija fluida jednaka u tački D kao delu J+1. cevi i tački D kao delu j-te cevi     , 1 , 1 , 1 , 1 , 1 , 1 , , 1 , , , , , , ; za 0 , ; za 0 P J P J D J P J P J D J D J D J P J P J D J P J P J D J h v p p Z t u h h h v p p Z t u                          (7.126) pri čemu se isključuje slučaj , , 10 0 .D J D Ju i u   Specifična zapremina u tački D se računa primenom jednačine stanja  ,D D Dv v p h . 109 8. KOMPJUTERSKI PROGRAM ZA SIMULACIJU HIDRAULIČKOG UDARA IZAZVANOG KONDENZACIJOM PARE Kompjuterski program nestacionarnog kompresibilnog strujanja, razvijen u cilju ispitivanja prostiranja temperaturskih talasa u fluidu, je napisan u kompjuterskom jeziku Turbo Pascal. Sastoji se od glavnog programa i nekoliko podprograma koji izvršavaju odreĎene operacije redom kojim je definisano u glavnom programu. 8.1 Opis programa Glavni program započinje učitavanjem ulaznih podataka. Iz ulazne datoteke se učitavaju: vremenski kraj računanja, ukupan broj cevi, dužina cevi i broj čvorova unutar svake cevi, zatim unutrašnji prečnik, koeficijent trenja, nagib cevi u odnosu na horizontalu, granični uslovi na početku i na kraju cevi za svaku cev; početne vrednosti pritiska, entalpije i brzine u svim čvorovima duž cevi; entalpija i pritisak fluida u rezervoaru, minimalne dozvoljene vrednosti entalpije i pritiska fluida, odnosno granice u kojima važe primenjene jednačine stanja. Ovakvom koncepcijom ulazne datoteke omogućeno je zadavanje željene konfiguracije cevne mreže i termohidrauličkih uslova. Nakon učitavanja ulaznih podataka se na početku svake iteracije, odnosno u svakom novom vremenskom trenutku, odreĎuje vremenski korak integracije, dt, prema Courant-ovom kriterijumu stabilnosti rešenja, jednačina (7.89), pri čemu se vremenski korak računa za svaku cev u posmatranom sistemu i onda se traži njegova najmanja vrednost u skupu svih dobijenih rezultata. Zatim je potrebno proveriti da li su parametri na početku cevi izračunati. Ovako se skraćuje vreme računanja jer kod spoja dve cevi, ako su izračunati parametri na kraju J+1 cevi, onda su izračunati i parametri na početku J-te cevi. Ako su na početku cevi parametri već izračunati, proverava se da li su izračunati parametri na kraju cevi. Ukoliko nisu izračunati parametri u oba slučaja, pristupa se njihovom izračunavanju iz odgovarajućih podprograma. Programom su obuhvaćeni granični uslovi: početak cevi spojen sa rezervoarom, kraj cevi spojen sa rezervoarom, zatvoreni kraj cevi na kraju 110 cevi i spoj dve cevi istog poprečnog preseka. U podprogramima „Rezervoar“ i „Rezervoar_na_kraju“ se računaju pritisak i entalpija onom čvoru kojim je cev spojena sa rezervoarom (u prvom, odnosno u poslednjem čvoru cevi), korišćenjem jednačina: (7.113) za pD i (7.114) za hD, dok se brzine fluida u prvom, odnosno u poslednjem čvoru cevi, izračunavaju primenom (7.115) i (7.120) za uD, sledstveno. Pritisak i entalpija fluida u rezervoaru su poznate veličine koje se učitavaju iz ulazne datoteke. Podprogram „KrajCevi_nizvodno“ se koristi za izračunavanje pritiska, entalpije i brzine fluida u poslednjem čvoru cevi sa zatvorenim krajem koristeći jednačine: (7.117) za uD , (7.118) za pD i (7.119) za hD. U podprogramu „Spoj2Cevi“ se izračunavaju pritisak, entalpija i brzina fluida na mestu gde se spajaju dve cevi istog poprečnog preseka bez gubitka strujne energije. Pri tome se koriste jednačine: (7.125) za pD, (7.126) za uD i (7.127) za hD. Algoritam izvršavanja podprograma „Spoj2Cevi“ prikazan je Slikom 8.1. Sada, kada su poznate sve potrebne veličine na granicama cevi, izračunavaju se parametri fluida: pritisak, entalpija i brzina, u svim čvorovima unutar svih cevi, ne uzimajući u obzir granične slučajeve, tj. prvi i poslednji čvor. U podprogramu „Cev“ se rešavaju jednačine (7.46) do (7.52) i pozivaju se podprogrami za izračunavanje strujnih i termodinamičkih parametara fluida u tačkama R, L i P: „RR“, „LL“ i „PP“, zatim podprogrami za izračunavanje desnih strana kvazilinearnih parcijalnih diferencijalnih jednačina hiperboličkog tipa: „XXX“, „YYY“, „ZZZ“. U podprogramu „RR“ se izračunavaju: brzina fluida u tački R, uR iz jednačine (7.57), brzina zvuka u tački R, cR iz (7.58), pritisak u tački R, pR iz (7.59) i specifična zapremina fluida u tački R, vR iz (7.60). U podprogramu „LL“ se izračunavaju: brzina fluida, brzina zvuka, pritisak i specifična zapremina fluida u tački L redom jednačinama (7.66) do (7.69). Podprogram „PP“ se koristi za odreĎivanje pritiska i entalpije fluida u tački P pomoću Lagrange- ovog interpolacionog polinoma 3. stepena. Prvo se odreĎuje brzina u tački P (Slika 7.2), jednačinama (7.76) ili (7.85), zavisno od toga da li je strujanje u čvoru B pozitivno ili negativno (smer strujanja jeste ili nije u pravcu ose, koja je orjentisana pri zadavanju ulaznih veličina). OdreĎuje se i x-koordinata tačke P. Zatim se prema znaku brzine uP odreĎuju koordinate četiri tačke (čvora) - jedan čvor nizvodno i dva uzvodno od posmatranog čvora uključujući i posmatrani čvor. Svakom odabranom čvoru se dodeljuju odgovarajuće vrednosti entalpija i pritisaka (veličine čiju je vrednost potrebno 111 POČETAK Da li se ra ? čunaju parametri na kraju cevi J=J+1 I=N J=J I=1 NE DA NE DA u <0B hD,J+1 za J=J+1,I=N se poziva Podprogram “PP” KRAJ za J=J+1,I=N se poziva Podprogram “RR” i računaju i iz jednačina (3.37) i (3.38)   za J=J,I=1 se poziva Podprogram “LL” i računaju i iz jednačina (3.39) i (3.40)   u =uD,J D,J+1 p =pD,J D,J+1 za J=J,I=1 se poziva Podprogram “PP” hD,J u ,uD,J D,J+1 D,J D,J+1 D,J D,J+1 p ,p h ,h Slika 8.1 Algoritam izvršavanja podprograma „Spoj2Cevi“. 112 odrediti u tački P). Tako definisane te četiri tačke se koriste pri formiranju Lagrange- ovog interpolacionog polinoma 3. stepena. Na kraju se poziva podprogram „Flag“ u kome se za odreĎenu koordinatu xP i definisane četiri tačke sa svojom apscisama i ordinatama računa suma (4.19). Lagrange-ovim interpolacionim polinomom 3. stepena se izmeĎu ove četiri tačke odreĎuje željeni parametar u tački xP, a to je u jednom slučaju pritisak, a u drugom slučaju entalpija u tački P. Ako se tačka P nalazi u blizini nekog od četiri izabrana čvora, ne primenjuje se Lagrange-ova interpolacija već entalpija i pritisak u tački P postaju jednaki vrednostima entalpije i pritiska u čvoru u čijoj se blizini tačka P nalazi. Na kraju se odreĎuje i specifična zapremina u tački P korišćenjem jednačine stanja. Podprogrami „XXX“, „YYY“ i „ZZZ“ služe za izračunavanje jednačina (7.25), (7.26) i (7.27). Iz podprograma „XXX“ i „ZZZ“ se poziva podprogram „KTR“ koji odreĎuje vrednost razlomka "ch   , iz jednačine zakona održanja energije (7.3), u funkciji ubrzanja razdelne površine tečne i parne faze. U programu je proces kondenzacije tako definisan da se odvija do trenutka kada se u neposrednoj blizini mesta o koje udara talas pritiska naĎe samo tečna faza, a to je ustvari znak da se sva para (ili njen najveći deo) kondenzovala i da je razdelna površina udarila u kraj cevi. Nestacionarni proces kondenzacije se ovde modelira na taj način što se u podprogramima „XXX“ i „ZZZ“ odvodi konstantna vrednost količine toplote od parne faze ili dvofazne mešavine, kojom je predstavljen intenzitet kondenzacije, a definisana je u podprogramu „KTR“. Ova količina toplote se odvodi od oblasti pare ili dvofazne mešavine koja se nalazi u kontaktu sa pothlaĎenom vodom i koja tokom vremena menja svoj položaj, a zatim se ista ta količina toplote dovodi uzvodno tečnoj fazi čime se postiže očuvanje zakona o održanju energije. Na Slici 7.1 je šematski prikaz modeliranja kondenzacije. Podprogrami za izračunavanje veličina stanja i termodinamičkih parametara obuhvataju podprograme za odreĎivanje temperature zasićenja, entalpija i specifičnih zapremina zasićene vode i zasićene pare, specifične zapremine vode, pare i mešavine 113 vode i pare, parcijalnog izvoda specifične zapremine po entalpiji pri konstantnom pritisku i brzina zvuka u zasićenoj vodi, zasićenoj pari, vodi, pari i mešavini vode i pare. Na kraju svake linije programa u kojoj se izračunavaju vrednosti zavisno promenjlivih veličina (p, h, u), pozivaju se podprogrami koji služe za proveru da li se dobijena vrednost nalazi unutar opsega u kom važe jednačine stanja koje se koriste u programu. Podprogram „ProvU“ se koristi za proveravanje izračunate vrednosti brzine fluida. Izračunata vrednost brzine fluida se poredi sa brzinom zvuka u toj istoj tački. Ako je izračunata vrednost brzine fluida veća od brzine zvuka (za vrednosti parametara u toj tački), brzina fluida u toj tački uzima vrednost brzine zvuka pri čemu zadržava svoj znak. Ovo se može dogoditi ako se simulira strujanje sa naglim promenama, pa vrednosti parametara u posmatranoj tački dobiju vrednosti na granici realnih vrednosti. Ovaj podprogram omogućava da se premosti moguća nastala numerička nestabilnost i održi fizikalnost rešenja, s obzirom da namena programa nije da proračunava nadzvučno strujanje. Podprogram „ProPri“ vrši poreĎenje izračunatih vrednosti pritiska na mestu odakle se ovaj podprogram poziva sa minimalnom dozvoljenom vrednosti pritiska koja se zadaje u ulaznoj datoteci. Ako je izračunata vrednost manja od zadate pritisak fluida na tom mestu postaje jednak zadatoj minimalnoj dozvoljenoj vrednosti. Minimalna vrednost pritiska se odreĎuje iskustveno, s obzirom na strujanje koje se proračunava, ali se prvenstveno mora imati u vidu opseg pritisaka za koji su napisani podprogrami koji na osnovu jednačina stanja izračunavaju ostale potrebne veličine (specifičnu zapreminu, brzinu zvuka...). Funkcija podprograma „ProEnt“ je ista kao kod podprograma „ProPri“ samo što se u njemu proverava izračunata vrednost entalpije. Na kraju svake iteracije, pošto se izračunaju vrednosti svih veličina u posmatranom vremenskom trenutku, vrši se ispisivanje dobijenih rezultata u izlaznu datoteku. Redosled ispisivanja i format u kome će se prikazati rezultati odreĎuje se u podprogramu „Izlaz“. 114 t=0 t=t+ t Podprogram “VREMENSKI KORAK” Podprogram “ ”CEV Podprogram “ ”IZLAZ u = u p = p h = h t t+ t t t+ t t t+ t    J=1,R0 J=1,R0 Podprogram “ULAZ” KRAJ POČETAK G[1,J]=3 G[1,J]=9 G[1,J]=11 G[1,J]=12 Podprogram “Spoj 2 cevi” Podprogram “Rezervoar” Podprogram “Rezervoar na kraju” Podprogram “Kraj cevi nizvodno” DA DA DA DA NE NE NE NENE Da li su i ? zračunati parametri na kraju cevi G[2,J]=3 G[2,J]=9 G[2,J]=11 G[2,J]=12 Podprogram “Spoj 2 cevi” Podprogram “Rezervoar” Podprogram “Rezervoar na kraju” Podprogram “Kraj cevi nizvodno” DA DA DA DA NE NE NE NENE DA DA Da li su i ? zračunati parametri na početku cevi NE DA t t< max Slika 8.2 Algoritam izvršavanja glavnog programa. 115 Dobijene vrednosti u ovom koraku postaju početne vrednosti za sledeći korak. Na kraju se vrši provera da li je računsko vreme dostiglo zadati vremenski kraj računanja i ako nije postupak se ponavlja za novi vremenski trenutak t+dt. Algoritam po kome se izvršava glavni program dat je na Slici 8.2. Glavni program sadrži i deo u kome se odreĎuje položaj razdelne površine (mesto u cevi gde se nalaze u dodiru tečna i parna faza) i definišu čvorovi koji odreĎuju oblasti iz kojih se odvodi i kojima se dovodi toplota. Položaj razdelne površine se odreĎuje na osnovu termodinamičkog stepena suvoće. Nije moguće precizno locirati položaj razdelne površine, nego samo odrediti dva čvora izmeĎu kojih se ona nalazi. Ovaj deo glavnog programa nije naznačen na algoritmu glavnog programa zbog jednostavnijeg snalaženja u njemu. 8.2 Modeliranje strujnog prostora Za potrebe programa kojim se simulira prostiranje temperaturskih talasa u fluidu potrebno je strujni prostor podeliti na konačan broj deonica (u programu J1,…,R0), koje su sastavljene od pravih cevi konstantnog poprečnog preseka. Da bi se uštedeo memorijski prostor, poželjno je da deonice budu približno jednake dužine jer se nizovi, koji se koriste u programu npr. p, u i h, dimenzionišu prema maksimalnom broju čvorova, odnosno prema cevi najveće dužine. Rastojanje izmeĎu čvorova je jednako za sve deonice, da bi se povećala tačnost pri izračunavanju, a odreĎuje se kao     0 , 1,..., 1 L J x J R n J     gde je [ ]L J dužina J- te cevi (m), a [ ]n J broj čvorova unutar J-te cevi. Ako prostorni korak integracije ima konstantnu vrednost (xconst. za svaku deonicu), onda je poželjno odrediti ga za najkraću cev, a zatim prilagoditi dužine ostalih deonica tom odreĎenom prostornom koraku, u onolikoj meri koliko to realni uslovi strujanja dozvoljavaju. Sledeći uslovi moraju da budu ispunjeni pri modeliranju nekog strujnog prostora da bi proračun mogao da se sprovede prema izloženom algoritmu 116 1. Deonice se obeležavaju opadajućim nizom prirodnih brojeva u pretpostavljenom smeru strujanja. Brojevi deonica idu od 1 do n, pri čemu se nijedan broj iz niza ne sme izostaviti. Jedna deonica ne može da ima dve različite oznake, niti dve različite deonice mogu imati istu oznaku. 2. Čvorovi se obeležavaju rastućim nizom prirodnih brojeva u pretpostavljenom smeru strujanja, pri čemu se nijedan broj u nizu ne sme preskočiti. Mesto prvog čvora u deonici se naziva početkom deonice, a mesto čvora označenog najvećim brojem krajem deonice. 3. Granični uslovi na početku i na kraju deonice se zadaju nizovima     01, , 2, , 1,..., ,G J G J J R pri čemu elementi prvog niza  1,G J označavaju vrstu graničnog uslova na početku cevi, a drugog niza  2,G J na kraju cevi. Vrsta graničnog uslova se odreĎuje tako što se elementima niza dodeljuju prirodni brojevi iz Tabele 8.1. Tabela 8.1 Vrste graničnih uslova Vrsta graničnog uslova Broj koji se dodeljuje elementima niza   0, , 1,2; 1,..., .G I J I J R  Spoj dve cevi 3 Rezervoar 9 Rezervoar na kraju cevi 11 Kraj cevi nizvodno 12 117 9. NUMERIČKA SIMULACIJA KRETANJA RAZDELNE POVRŠINE I HIDRAULIČKOG UDARA IZAZVANOG KONDENZACIJOM PARE Razvijeni modeli i kompjuterski program nestacionarnog kompresibilnog strujanja primenjeni su za proračune strujanja kada se u neposrednom kontaktu nalaze pothlađena voda i vodena para i pri kojima dolazi do hidrauličkog udara usled intenzivne kondenzacije pare. Validacija razvijenog programa izvršena je poređenjem dobijenih numeričkih rezultata sa dostupnim rezultatima eksperimentalnih merenja iz literature (Yeung i dr., 1993) i poređenjem sa numerički dobijenim rezultatima primenom komercijalnih programa TUF (Liu i dr., 1997) i RELAP5/MOD3 (Yeung i dr., 1993). U cilju provere mogućnosti modela da predvidi kretanje razdelne površine vodenog stuba i parnog prostora u cevovodu, simulirani su idealizovani uslovi strujanja za koje je moguće naći analitička rešenja (Jun, 1987). U cilju predviđanja navedenih termičkih i hidrauličkih uslova, razvijeni numerički postupak obuhvata i modele kondenzacije i nestacionarnog trenja, kao i model za praćenje kretanja razdelne površine tečne i parne faze. Izvršene su tri grupe proračuna:  numerička simulacija kretanja razdelne površine u oscilujućem manometru koji se sastoji od U-cevi. Masa vode u manometru se izvodi iz ravnoteže zadavanjem početne brzine (Jun, 1987),  numerička simulacija kretanja razdelne površine i hidrauličkog udara izazvanog kondenzacijom pare u jednostavnoj eksperimentalnoj aparaturi koja se sastoji od rezervoara, horizontalne i vertikalne cevi. Horizontalna cev je jednim svojim krajem spojena sa rezervoarom. Rezervoar i horizontalna cev su ispunjeni pothlađenom vodom. U vertikalnoj cevi, koja je spojena brzo delujućim ventilom sa horizontalnom cevi, a na svom drugom kraju je zatvorena, nalazi se para u stanju zasićenja. Otvaranjem ventila dolaze u kontakt pothlađena voda i para i nastaje hidraulički udar. Proračuni su sprovedeni za dva različita nivoa pritisaka u sistemu (Liu i dr., 1997) i 118  numerička simulacija hidrauličkog udara izazvanog kondenzacijom pare u eksperimentalnoj aparaturi sačinjenoj od vertikalne cevi, ispunjene parom u stanju zasićenja i zatvorene na jednom kraju, a drugim krajem uronjene u bazen ispunjen pothlađenom vodom. Direktnim kontaktom pothlađene vode i pare nastaje hidraulički udar. 9.1 Kretanje razdelne površine tečne i parne faze u oscilujućem manometru Aparatura, prikazana na Slici 9.1, se sastoji od U-cevi spojene na vrhu tako da je formiran zatvoren sistem (Jun, 1987). Ukupna dužina manometra je 20 m, a prečnik U- cevi je 1 m. U početnom trenutku sistem je ispunjen vodom i vazduhom, pri čemu voda ispunjava i levi i desni krak manometra do visine 5 m od dna, brzina fluida u manometru je 2,1 m/s, a ubrzanje je jednako nuli. Pomeranje razdelne površine x se meri od početnog položaja kada su nivoi vode u oba kraka jednaki. Cilj proračuna je testiranje numerički dobijenih rezultata očuvanja mase u sistemu, koja je konstantna, kao i praćenje kretanja razdelne površine vode i vazduha i modeliranje periode oscilovanja, čije je analitičko rešenje poznato. U početnom vremenskom trenutku temperatura je u svakoj tački sistema jednaka, odnosno voda i vazduh se unutar U-cevi nalaze na istoj temperaturi, koja iznosi 50 0C. Pritisak na razdelnoj površini vode i vazduha, kao i u prostoru iznad nje, je jednak atmosferskom pritisku. U delu U-cevi ispunjenom tečnošću pritisak ima različite vrednosti u skladu sa odgovarajućim vrednostima hidrostatičkog pritiska, dok je pritisak u delu koji je ispunjen vazduhom konstantan. Pošto razvijeni program nema u sebi ugrađene odgovarajuće jednačine za određivanje termodinamičkih veličina stanja vazduha već samo relacije za određivanje stanja vode i vodene pare, potrebno je izvršiti određeno prilagođavanje parametara ulaznih veličina (Prica i dr., 2008). U programu je izvršen proračun sa vodenom parom, umesto sa vazduhom, na pritisku na kom para ima istu gustinu kao vazduh u posmatranom eksperimentu. Uzeta je nešto veća vrednost entalpije pare od entalpije pare u stanju zasićenja, kako bismo bili sigurni da vodena para neće preći u oblast 119 vlažne pare prilikom računa. Zbog gore navedenih razloga usvojeno je da je vrednost pritiska na razdelnoj površini 0,1936 MPa. Pošto u ovom eksperimentu stišljivost vode ne dolazi do izražaja, ova razlika pritisaka između pritiska pare u programu i pritiska vazduha u eksperimentu, koja u konkretnom slučaju iznosi  0,1 MPa (0,0936 MPa), praktično ne izaziva promenu gustine vode. 1 0 m x 0 5 m vazduh voda 1 m Slika 9.1 Shematski prikaz aparuture oscilujućeg manometra. Tabela 9.1 Početni uslovi kretanja razdelne površine u oscilujućem manometru: stanje vode p1  0,1936 MPa h1  209,3 kJ/kg stanje pare p2  0,1936 MPa h2 > h’’(p2)  2800 kJ/kg 120 9.1.1 Diferencijalna jednačina kretanja razdelne površine u oscilujućem manometru U oscilujućem manometru na kretanje razdelne površine utiče samo Zemljina gravitacija (posmatra se slučaj bez trenja), a izazvano je zadavanjem kretanja vode i vazduha konstantnom i međusobno jednakom brzinom u početnom trenutku. Slika 9.2 predstavlja shematski prikaz stanja u U-cevi u nekom vremenskom trenutku, pošto je započelo kretanje vode. Diferencijalna jednačina kretanja razdelne površine pod dejstvom gravitacione sile je ,um x mg  (9.1) gde je mu ukupna masa vode u sistemu 2 u L m A (kg), m masa vode koja predstavlja razliku masa vode u jednom i u drugom kraku manometra 2m A xg (kg), A površina poprečnog preseka U-cevi 2 4 d A   (m2), x položaj razdelne površine u posmatranom trenutku (m). U konačnom obliku diferencijalna jednačina kretanja razdelne površine ima oblik 4 0, g x x L   (9.2) a početni uslovi za njeno rešavanje su: x|t0 0, x |t0  dx dt  u0 2,1 m/s. Jednačina (9.2) je homogena linearna diferencijalna jednačina 2. reda sa konstantnim koeficijentima čije opšte rešenje (Mamuzić, 1981) glasi      1 1 2 2x t C x t C x t  . (9.3) Pošto su koeficijenti u jednačini (9.2) realne konstante njena rešenja se traže u obliku   , 1,2itix t e i   , (9.4) 121 pa opšte rešenje (9.3) dobija oblik   1 21 2 t t x t C e C e    . Uvrštavanjem   tx t e u (9.2) se dobija kvadratna jednačina koja se zove karakteristična jednačina diferencijalne jednačine (9.2) 2 4g L    . (9.5) Diskriminanta karakteristične jednačine (9.5) je manja od nule pa su njena rešenja kompleksni brojevi 1,2 ,a ib   gde je a realni deo kompleksnog broja, a b imaginarni deo kompleksnog broja. U ovom slučaju je realni deo kompleksnog broja jednak nuli, a imaginarni 4g b L  , pa se dobija opšte rešenje jednačine (9.2) 4 4 1 2 . g g it it L Lx C e C e    (9.6) L L /2 vazduh voda d t=0, x0 t, x x 2 x Slika 9.2 Shematski prikaz stanja u U - cevi u posmatranom vremenskom trenutku. 122 Koristeći Euler-ovu formulu cos sin ,ie i    (9.7) i određujući konstante iz početnih uslova dobija se analitičko rešenje kretanja razdelne površine u konačnom obliku 0 4 sin( ), 4 L g x u t g L  (9.8) a brzina kretanja razdelne površine u toku vremena je jednaka 0 4 cos( ). dx g u u t dt L   (9.9) 9.1.2 Grafički prikaz dobijenih rezultata pri numeričkoj simulaciji kretanja razdelne površine u oscilujućem manometru Zbog jednostavnosti diferencijalne jednačine pomeranja razdelne površine u slučaju oscilujućeg manometra, rešenja dobijena primenom metode karakteristika su poređena sa analitičkim rešenjem ove diferencijalne jednačine. Izvršeno je i poređenje rezultata dobijenih primenom metode Runge-Kutta i analitičkog rešenja kako bi se pokazala izuzetna tačnost ove metode i opravdala njena upotreba, gde nije moguće odrediti rešenje analitičkim putem, kao metode koja služi za proveru tačnosti numeričkog rešenja diferencijalne jednačine dobijenog drugim metodama. Krećući se brzinom od 2,1 m/s u smeru x koordinate (Slika 9.1) vodeni stub gura vazduh u levom kraku manometra naviše, zatim se pod dejstvom Zemljine gravitacije kreće vertikalno naniže i gura vazduh u desnom kraku prema vrhu manometra. Postupak se neprestalno ponavlja, pošto se posmatra slučaj kada nema trenja o zidove manometra. Rešavanjem diferencijalne jednačine kretanja razdelne površine (9.2) primenom metode karakteristika i analitičkim putem dobija se promena položaja razdelne površine vode i vazduha i masenog protoka na razdelnoj površini u toku vremena, što je prikazano na Slikama 9.3 i 9.5. Vrednosti dobijene na ova dva načina su međusobno upoređene i dobijeno je zadovoljavajuće slaganje. Dobijeni period oscilovanja je isti primenom obe metode pri promeni položaja razdelne površine, dok se u slučaju 123 promene masenog protoka analitičkim rešavanjem dobija za 4% veći period oscilovanja. Zbog toga se u toku vremena povećava vremenski razmak između odgovarajućih amplituda dobijenih primenom različitih metoda, Slika 9.5. Amplitude oscilovanja, dobijene primenom metode karakteristika i rešavanjem diferencijalne jednačine kretanja analitičkim putem, su jednake pri promeni masenog protoka na razdelnoj površini vode i vazduha. Pomeranje razdelne površine u toku vremena dobijeno primenom metode karakteristika ima veće amplitude od analitičkog rešenja pomeranja (Slika 9.3). Dijagrami na Slikama 9.4 i 9.6 prikazuju promene istih veličina kao na dijagamima na Slikama 9.3 i 9.5. Ovde je izvršeno poređenje vrednosti dobijenih rešavanjem diferencijalne jednačine (9.2) numeričkim putem primenom metode Runge- Kutta i analitičkim putem. Pošto je poklapanje rezultata idealno, svaka vrednost dobijena numeričkim putem metodom Runge-Kutta nalazi se na liniji koja predstavlja analitičko rešenje. Ovim je pokazana mogućnost primene metode Runge-Kutta kao metode kojom će se vršiti poređenje dobijenih rezultata primenom neke druge numeričke metode kada je teže, ili čak nemoguće dobiti analitičko rešenje posmatrane, tj. rešavane, diferencijalne jednačine. Da bi se ilustrovao značaj primene interpolacionog polinoma višeg reda, pri određivanju početnih vrednosti entalpije i pritiska između čvorova, sproveden je i proračun sa linearnom interpolacijom. Na Slici 9.7 se jasno vidi da primena linearnog interpolacionog polinoma daje znatno veću grešku pri određivanju položaja razdelne površine. 124 Slika 9.3 Promena položaja razdelne površine u toku vremena u U – cevi manometra. U problemu nije uzeto u obzir trenje o zidove cevi. Poređenje analitičkog rešenja i rezultata dobijenih numeričkim rešavanjem primenom metode karakteristika. Slika 9.4 Poređenje analitičkog rešenja promene položaja razdelne površine u toku vremena u U - cevi manometra i rezultata dobijenih numeričkim rešavanjem metodom Runge-Kutta. U problemu nije uzeto u obzir trenje o zidove manometra. 125 Slika 9.5 Promena masenog protoka na razdelnoj površini dveju faza u toku vremena u U – cevi manometra. U problemu nije uzeto u obzir trenje o zidove cevi. Poređenje analitičkog rešenja i rezultata dobijenih numeričkim rešavanjem primenom metode karakteristika. Slika 9.6 Poređenje analitičkog rešenja promene masenog protoka na razdelnoj površini dveju faza u toku vremena u U - cevi manometra i rezultata dobijenih numeričkim rešavanjem metodom Runge-Kutta. U problemu nije uzeto u obzir trenje o zidove manometra. 126 Slika 9.7 Poređenje rešenja promene položaja razdelne površine u toku vremena u U – cevi manometra dobijenih numeričkim rešavanjem primenom metode karakteristika uz primenu Lagrange-ovog polinoma 3. stepena i linearne interpolacije na određivanje početnih vrednosti entalpije i pritiska između čvorova numeričke mreže. Trenje o zidove manometra nije uzeto u obzir. 9.2 Hidraulički udar usled kondenzacije pare na stubu pothlaĎene vode Eksperimentalno istraživanje hidrauličkog udara izazvanog intenzivnom kondenzacijom pare je izvršeno na aparaturi prikazanoj na Slici 9.8 (Liu i dr., 1997). Sistem se sastoji od rezervoara, horizontalne cevi dužine 5,5 m i vertikalne cevi dužine 5 m. Obe cevi imaju jednake prečnike od 0,0921 m. Rezervoar i horizontalna cev, koja je povezana sa dnom rezervoara, u početnom vremenskom trenutku su ispunjeni pothlađenom vodom, a u vertikalnoj cevi se nalazi para u stanju zasićenja. Vertikalna cev je svojim početkom povezana sa horizontalnom cevi preko brzo delujućeg loptastog ventila, a na svom kraju je zatvorena. Merači pritiska su postavljeni neposredno ispred zatvorenog kraja vertikalne cevi. Usled brzog otvaranja loptastog ventila pothlađena voda dolazi u direktan kontakt sa parom pri čemu dolazi do njene intenzivne kondenzacije na čelu vodenog stuba što uzrokuje nastanak hidrauličkog udara. Izvršene su dve grupe proračuna za različite vrednosti pritisaka u sistemu. Početni uslovi 127 proračuna su dati u Tabeli 9.2. U oba slučaja sistem miruje u početnom trenutku, nema kretanja ni vodenog stuba ni pare, u  0. p 5.5m 5 m Slika 9.8 Shematski prikaz eksperimentalne aparature za ispitivanje hidrauličkog udara. Tabela 9.2 Početni uslovi hidrauličkog udara izazvanog intenzivnom kondenzacijom pare: 1. slučaj stanje vode: p1  5,51 bar t1  22 0 C h1  92,802 kJ/kg stanje pare: p2  3,82 bar t2  tsat(p2) h2  h”(p2) 2. slučaj stanje vode: p1  6,54 bar t1  23 0 C h1  97,0 kJ/kg stanje pare: p2  4,76 bar t2  tsat(p2) h2  h”(p2) 9.2.1 Diferencijalna jednačina kretanja razdelne površine u stubu tečnosti promenljive visine u vertikalnoj cevi Na Slici 9.9 je predstavljen shematski prikaz stanja u vertikalnoj cevi u nekom vremenskom trenutku t, kada je vertikalna cev ispunjena do visine x pothlađenom vodom, a iznad vode se nalazi para u stanju zasićenja. U ovom slučaju kretanje razdelne 128 površine je izazvano razlikom pritisaka u pari i u vodi, odnosno razlikom pritisaka ispred i iza razdelne površine posmatrano u smeru strujanja. L -x1 x L1 p(x) pres m=m(x) m Ar d 1 Slika 9.9 Shematski prikaz stanja u vertikalnoj cevi u posmatranom vremenskom trenutku. Voda u horizontalnoj cevi se nalazi na pritisku jednakom pritisku u rezervoaru (Slika 9.8), a para u vertikalnoj cevi je na nižem pritisku. Usled ovoga dolazi do uticanja vode iz horizontalne u vertikalnu cev i potiskivanja pare, odnosno pomeranja razdelne površine. Pritisak u horizontalnoj cevi je sve vreme konstantan, dok se pritisak u vertikalnoj cevi menja. Da bi se odredila ova promena pritiska, uzeta je pretpostavka da se vodena para u posmatranim uslovima (nema kondenzacije) ponaša kao idealan gas pa važi jednačina stanja idealnog gasa ,gpV mR T (9.9) gde je Rg gasna konstanta vodene pare (J/(kgK)), m masa vodene pare (kg), V zapremina vodene pare (m 3 ). Pri adijabatskom sabijanju gasa važi jednačina adijabate .pv const  (9.10) 129 Kao i pritisak, i zapremina vodene pare je promenljiva veličina i zavisi od nivoa vode u cevi x, odnosno položaja razdelne površine 1( ),rV A L x  (9.11) gde je L1 dužina vertikalne cevi (m), Ar razdelna površina, 2 1 4 r d A   (m2), d1 prečnik vertikalne cevi (m). Diferenciranjem jednačine (9.9) dobija se .gpdV Vdp mR dT  (9.12) Pošto se rezultat traži u obliku pp(x), a i VV(x), potrebno je iz jednačine (9.12) eliminisati T. Diferenciranjem jednačine (9.10) i uvrštavanjem specifične zapremine v=V/m, iz jednačine (9.9) dobija se 1 , T dT dp p      (9.13) odnosno, jednačina (9.12) , uz smenu (9.13) i / ,gV mR T p postaje . dV dp V p   (9.14) Da bi se dobio rezultat u obliku pp(x), u (9.14) se uvodi smena dobijena diferenciranjem jednačine (9.11) ,rdV A dx  (9.15) i 1( )rV A L x  iz jednačine (9.11), tako da jednačina (9.14) postaje 1 . dp dx p L x   (9.16) Integracijom jednačine (9.16) od početnog do posmatranog stanja sledi promena pritiska pare u funkciji pomeranja fronta vodenog talasa   10 1 , L x p x p L         (9.17) gde je p0 – pritisak pare u početnom trenutku (Pa). 130 Sada se može napisti diferencijalna jednačina kretanja razdelne površine pod dejstvom promene pritiska i gravitacione sile (posmatra se slučaj bez trenja)      ,res rMx p p x A m x g   (9.18) gde je M ukupna masa vode (kg), koja je jednaka zbiru mase vode koja ulazi u vertikalnu cev m(x)Arx (kg) i mase vode u horizontalnoj cevi m2A2L2 (kg), pri čemu je A2Ar pošto su obe cevi istog prečnika. Sređivanjem jednačine (9.18) dobija se jednačina pomeranja razdelne površine u obliku pogodnom za rešavanje primenom metode Runge-Kutta.     2 2 . resp p x xg x x L x L      (9.19) Početni uslovi su sledeći:   00 0 0 m 0, 0 , 3,82 bar st t t dx x x p x p dt        u prvom slučaju, a   00 4,76 bartp x p   u drugom posmatranom slučaju. 9.2.2 Grafički prikaz dobijenih rezultata pri numeričkoj simulaciji hidrauličkog udara usled kondenzacije pare na stubu pothlaĎene vode Izvršeni su proračuni za dve grupe početnih vrednosti pritisaka u sistemu. Prva grupa dijagrama predstavlja rezultate kada je pritisak u rezervoaru jednak 5,51 bar, a u vertikalnoj cevi 3,82 bar, dok drugu grupu čine rezultati dobijeni u slučaju kada je pritisak u rezervoaru jednak 6,54 bar, a u vertikalnoj cevi 4,76 bar. Na Slikama 9.10 i 9.11 su prikazane promene položaja razdelne površine tečne i parne faze i promene pritiska u aparaturi u toku vremena, pri čemu su zanemareni uticaji kondenzacije i trenja. Ovi rezultati su dobijeni rešavanjem diferencijalne jednačine kretanja razdelne površine primenom dveju numeričkih metoda, metode karakteristika i metode Runge–Kutta, za obe navedene grupe početnih uslova. Usled razlike pritisaka u rezervoaru i u vertikalnoj cevi, pri čemu je pritisak u rezervoaru viši od pritiska u vertikalnoj cevi, razdelna površina tečne i parne faze započinje kretanje ka zatvorenom kraju vertikalne cevi. Krećući se ka zatvorenom kraju stub tečnosti sabija 131 paru i pritisak pare raste. Nastali talas pritiska prostire se do zatvorenog kraja cevi, odbija se od njega i kreće ka rezervoaru. Amplitude oscilovanja sve vreme ostaju konstantne jer je uticaj trenja o zidove cevi zanemaren. U rezultatima dobijenim primenom metode karakteristika, na rešavanje promene položaja razdelne površine, primećuje se mali pad vrednosti amplitude što je posledica numerike. Postignuto je zadovoljavajuće slaganje rezultata dobijenih metodom Runge-Kutta i metodom karakteristika. Na dijagramima je vidljivo odstupanje amplituda i perioda oscilovanja rešenja dobijenog metodom karakteristika od rešenja dobijenog metodom Runge-Kutta. Mogući razlog ovog odstupanja je u tome što se prilikom izvođenja diferencijalne jednačine kretanja razdelne površine, (9.19), uvodi gruba aproksimacija da se vodena para u vertikalnoj cevi ponaša kao idealan gas. Otvaranjem brzo delujućeg ventila (Slika 9.8) u trenutku t=0 u neposredni kontakt su dovedene pothlađena voda i para u stanju zasićenja usled čega dolazi do intenzivne kondenzacije i vodena masa počinje da se kreće. Usled kondenzacije pritisak u parnom prostoru naglo opada i stub tečnosti se ubrzava ka zatvorenom kraju vertikalne cevi (Slika 9.8). Kada se sva para iz vertikalne cevi kondenzuje, čelo vodenog stuba udara u zatvoreni kraj cevi u 1,34 s i izaziva porast pritiska od 86 bar (Slika 9.12a). Posle udara o zatvoreni kraj cevi, talas pritiska se prostire ka rezervoaru, gde se odbija od vodene mase usled čega mu slabi intenzitet. Ovakav usporeni talas se ponovo kreće ka zatvorenom kraju vertikalne cevi i udara o njega, ali ovaj put oslabljenog intenziteta. Ovaj se proces periodično ponavlja. Amplitude talasa pritiska s vremenom slabe zbog trenja o zidove cevi i elastičnih deformacija cevi. Numerički dobijeni rezultati se dobro slažu sa izmerenim vrednostima (Slika 9.12b). 132 a) b) Slika 9.10 Pomeranje položaja razdelne površine (gore) i promena pritiska (dole) u toku vremena u vertikalnoj cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8). Slučaj kada je pritisak u rezervoaru 5,51 bar, a u vertikalnoj cevi 3,82 bar, kada su zanemareni uticaji kondenzacije i trenja. Poređenje numeričkih rešenja dobijenih metodom karakteristika i metodom Runge-Kutta. 133 a) b) Slika 9.11 Pomeranje položaja razdelne površine (gore) i promena pritiska (dole) u toku vremena u vertikalnoj cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8). Slučaj kada je pritisak u rezervoaru 6,54 bar, a u vertikalnoj cevi 4,76 bar, kada su zanemareni uticaji kondenzacije i trenja. Poređenje numeričkih rešenja dobijenih metodom karakteristika i metodom Runge-Kutta. 134 a) b) Slika 9.12 Promena pritiska u toku vremena u blizini zatvorenog kraja vertikalne cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8) dobijena primenom razvijenog programa (gore) i rezultati eksperimentalnih merenja (Liu i dr., 1997). Slučaj kada je pritisak u rezervoaru 5,51 bar, a u vertikalnoj cevi 3,82 bar, kada su uzeti u obzir uticaji kondenzacije i trenja. 135 Na Slici 9.13 je prikazano prostiranje čela vodenog stuba za vreme hidrauličkog udara. Na početku eksperimenta vertikalna cev je bila ispunjena parom u stanju zasićenja, a čelo vodenog stuba se nalazilo na ulazu u vertikalnu cev. U trenutku potpune kondezacije pare u 1,34 s, vodeni stub dolazi do zatvorenog kraja vertikalne cevi i udara o njega. Intenzivno prostiranje talasa pritiska kroz eksperimentalnu aparaturu uzrokuje lokalni pad pritiska ispod vrednosti pritiska u stanju zasićenja i posledično adijabatsko isparavanje vode usled čega nastaju male zapremine pare na vrhu vertikalne cevi. Prema sprovedenim proračunima posle 2,7 s (Slika 9.13) nema više adijabatskog isparavanja. Slika 9.13 Pomeranje položaja razdelne površine u toku vremena u vertikalnoj cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8) dobijeno primenom razvijenog programa. Slučaj kada je pritisak u rezervoaru 5,51 bar, a u vertikalnoj cevi 3,82 bar, kada su uzeti u obzir uticaji kondenzacije i trenja. Na Slici 9.14 je prikazan nagli porast brzine čela vodenog stuba na početku nestacionarne promene zbog intenzivne kondenzacije pare na pothlađenoj vodi posle otvaranja ventila. Brzina čela vodenog stuba dostiže vrednost od 5 m/s i osciluje oko ove vrednosti. Na dijagramu su jasno uočljiva dva perioda neznatnog opadanja brzine. U trenutku udara čela vodenog stuba o zatvoreni kraj cevi se dešava iznenadni pad brzine sa promenom smera brzine. Do udara čela vodenog stuba o zatvoreni kraj cevi se 136 vodeni stub kretao vertikalno naviše, u trenutku udara ima brzinu jednaku nuli, a posle udara se kreće u suprotnom smeru, vertikalno naniže. Slika 9.14 Promena brzine čela vodenog stuba u toku vremena dobijena primenom razvijenog programa. Slučaj kada je pritisak u rezervoaru 5,51 bar, a u vertikalnoj cevi 3,82 bar, kada su uzeti u obzir uticaji kondenzacije i trenja. Pad brzine čela vodenog stuba se može iskoristiti za procenu pika pritiska prema korelaciji Joukowsky (1904) p c u    . (9.20) Za gustinu vode od 1000 kg/m3, brzinu zvuka 1500 m/s, i promenu brzine -5 m/s, dobija se pik pritiska 75 bar što je blisko predviđanju numeričkog modela, 84 bar (Slika 9.12a). Na Slici 9.15 je predstavljena promena pritiska u toku vremena, u blizini zatvorenog kraja vertikalne cevi, u aparaturi za simulaciju hidrauličkog udara (Slika 9.8), bez uzimanja u obzir uticaja nestacionarnog trenja u proračunskom modelu. Nedostatak nestacionarnog trenja dovodi do smanjenog prigušivanja oscilacija talasa pritiska. Izmerene vrednosti pritiska (Slika 9.12b) pokazuju prigušivanje talasa pritiska posle 4 s. Izračunate vrednosti primenom metode karakteristika uz uzimanje u obzir 137 nestacionarnog trenja (Slika 9.12b), pokazuju prigušivanje posle približno 3,7 s, dok proračuni bez uključivanja nestacionarnog trenja ne predviđaju prigušivanje oscilacija talasa pritiska čak ni posle 5 s. Slika 9.15 Promena pritiska u toku vremena u blizini zatvorenog kraja vertikalne cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8). Pritisak u rezervoaru je 5,51 bar, a u vertikalnoj cevi 3,82 bar. Rezultati dobijeni primenom razvijenog programa pri čemu je zanemaren uticaj nestacionarnog trenja. Uticaj kondenzacije je uzet u obzir. Na Slikama 9.16 i 9.17 je ilustrovan značaj primene numeričke šeme višeg reda tačnosti (Lagrange-ov interpolacioni polinom trećeg reda) za određivanje početnih vrednosti entalpije fluidnog delića, na karakteristici koja predstavlja prostiranje fluidnog delića. Primena linearne interpolacije rezultuje dužim vremenskim periodima između pojave pikova pritiska, nižom tačnošću predviđenih amplituda pikova i nedostatkom prigušivanja talasa pritiska (Slika 9.16). Posledica primene linearne interpolacije za određivanje entalpije na karakteristici fluidnog delića je i značajno pomeranje čela vodenog stuba posle nastanka prvog pika pritiska, dok primena Lagrange-ove interpolacije trećeg reda ne predviđa postojanje pare u vertikalnoj cevi posle približno 2.7 s (Slika 9.17). 138 Slika 9.16 Promena pritiska u toku vremena u blizini zatvorenog kraja vertikalne cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8). Pritisak u rezervoaru je 5,51 bar, a u vertikalnoj cevi 3,82 bar. Rezultati dobijeni primenom razvijenog programa. Poređenje rezultata dobijenih primenom Lagrange-ovog interpolacionog polinoma trećeg reda i linearne interpolacije na određivanje početne vrednosti entalpije na karakteristici fluidnog delića. Uticaj kondenzacije i nestacionarnog trenja je uzet u obzir. Test kojim je proveren uticaj promene numeričke mreže na rezultate proračuna pokazuje da su slični rezultati dobijeni primenom koraka 0,1 m i 0,05 m između dva susedna čvora proračunske mreže (Slika 9.18). Broj čvorova u vertikalnoj cevi eksperimentalne aparature, prikazane na Slici 9.2, je sledstveno 101 i 51. Dalje smanjenje rastojanja između čvorova ne utiče na rezultate, ali utiče na vremensko trajanje proračuna s obzirom da je metoda karakteristika eksplicitna numerička metoda u kojoj je odnos vremenskog i prostornog koraka integracije definisan Courant-ovim kriterijumom (7.78). Takođe je izvršen proračun za eksperimentalne uslove sa nešto višim početnim pritiskom vode od 6,54 bar. Izmerene vrednosti, predstavljene na Slici 9.19b, pokazuju 139 da je pik pritiska 162 bar. Rezultati dobijeni primenom razvijenog programa (Slika 9.19a), pokazuju dobro slaganje sa izmerenim vrednostima. Ako se uporede oba slučaja, može se zaključiti da porast pritiska od 1,03 bar (sa 5,51 bar na 6,54 bar, pri čemu su ostali početni uslovi približno jednaki), udvostručuje porast pritiska pri udaru vodenog stuba (Slike 9.12 i 9.19). Slika 9.17 Pomeranje položaja razdelne površine u toku vremena u vertikalnoj cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8). Pritisak u rezervoaru je 5,51 bar, a u vertikalnoj cevi 3,82 bar. Poređenje rezultata dobijenih primenom razvijenog programa i Lagrange-ovog interpolacionog polinoma trećeg reda i linearne interpolacije na određivanje početne vrednosti entalpije na karakteristici fluidnog delića. Uticaj kondenzacije i nestacionarnog trenja je uzet u obzir. Autori eksperimenata hidrauličkog udara izazvanog intenzivnom kondenzacijom pare (Zaltsgendler i dr., 1996) su objavili da je pri istim početnim uslovima eksperimenta dobijeno značajno rasipanje izmerenih rezultata. Ovo rasipanje je posledica stohastičke prirode formiranja razdelne površine tečne i parne faze, koje je uzeto u obzir preko parametra aK,0 u jednačini (7.15) poglavlje 7. Oba eksperimenta koja su numerički simulirana i čiji su rezultati dati u ovom poglavlju su simulirani sa vrednošću ,0 40Ka  m 2 /m 3 . 140 a) b) Slika 9.18 Provera uticaja promene prostornog koraka numeričke mreže na rezultate proračuna dobijene primenom razvijenog programa. Promena pritiska u toku vremena u blizini zatvorenog kraja vertikalne cevi (gore) i pomeranje položaja razdelne površine u toku vremena (dole), u vertikalnoj cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8). Pritisak u rezervoaru je 5,51 bar, a u vertikalnoj cevi 3,82 bar. Uticaj kondenzacije i nestacionarnog trenja je uzet u obzir. 141 a) b) Slika 9.19 Promena pritiska u toku vremena u blizini zatvorenog kraja vertikalne cevi u aparaturi za simulaciju hidrauličkog udara (Slika 9.8) dobijena primenom razvijenog programa (gore) i rezultati eksperimentalnih merenja (Liu i dr., 1997) Slučaj kada je pritisak u rezervoaru 6,54 bar, a u vertikalnoj cevi 4,76 bar, kada su uzeti u obzir uticaji kondenzacije i trenja. 142 9.3 Hidraulički udar u vertikalnoj cevi za ispuštanje pare u bazen sa pothlaĎenom vodom Hidraulički udar u vertikalnoj cevi za ispuštanje pare u bazen sa pothlađenom vodom se naziva još vodenim topom, s obzirom da voda u zatvoreni kraj udara velikom brzinom, što dovodi do naglog porasta pritiska. Ovu pojavu su eksperimentalno istraživali Yeung i dr. (1993). Eksperimentalna aparatura, prikazana na Slici 9.20, se sastoji od vertikalne cevi dužine 0,7112 m i unutrašnjeg prečnika 0,0381 m, izrađene od metala i velikog rezervoara sa vodom u koji je cev uronjena svojim dnom do dubine od nekoliko centimetara. Merači pritiska postavljeni su na vrhu cevi. 0 .7 1 1 2 m 0.0381 m p para Slika 9.20 Shematski prikaz eksperimentalne aparuture za simulaciju vodenog topa. Vodena para u stanju zasićenja na pritisku 1,023 bar se uvodi u cev konstantnom brzinom kroz cevčicu, koja se nalazi na vrhu cevi, a donji kraj cevi je uronjen u veliki rezervoar vode temperature 49 0 C, na istom pritisku 1,023 bar. Posmatran je problem od trenutka kada se u vertikalnoj cevi nalazi samo para, a ventil na vrhu cevi je zatvoren. Početna brzina u svim delovima sistema je jednaka 0 m/s. Ovakva situacija se može javiti na izlazu iz parne turbine, prilikom ispuštanja pare u veću zapreminu hladne vode. Direktan kontakt pare i hladne vode vodi do intenzivne kondenzacije pare i nastanka 143 hidrauličkog udara. Pothlađena voda utiče u cev u obliku čepa, koji udara u zatvoreni ventil za ispuštanje pare, pri čemu nastaje veliki porast pritiska na ventilu. Ovakav hidraulički udar može da dovede do oštećenja opreme i obustave rada elektrane. Tabela 9.3 Početni uslovi za simulaciju vodenog topa: stanje vode: p1  1,023 bar t1  49 0 C h1  209,35 kJ/kg stanje pare: p2  1,023 bar t2  tsat(p2) h2  h’’(p2) 9.3.1 Grafički prikaz dobijenih rezultata pri numeričkoj simulaciji hidrauličkog udara u vertikalnoj cevi za ispuštanje pare u bazen sa pothlaĎenom vodom Pri kontaktu pare i pothlađene vode dolazi do intenzivne kondenzacije pare. Pritisak u pari opada i stub tečnosti počinje da ulazi u cev. Porastom nivoa vode u cevi, usled razlike pritisaka, proces kondenzacije se prenosi dalje ka njenom vrhu gde na kraju vodeni stub i udara. Maksimalni intenzitet pritiska nastaje pri prvom udaru. Na Slici 9.21 su prikazani rezultati promene pritiska u toku vremena u blizini zatvorenog kraja vertikalne cevi izazvane hidrauličkim udarom. Ovi pikovi pritiska su izazvani udarima vodenog stuba o zatvoreni kraj cevi. Rezultati dobijeni primenom razvijenog programa (Slika 9.21a) su poređeni sa rezultatima programa RELAP (Yeung i dr., 1993) (Slika 9.21b), koji je razvijen za sigurnosne termohidrauličke analize procesa u nuklearnim elektranama. Intenzitet i vremenski trenutak nastanka prvog pika pritiska, koji predstavlja opasnost od izazivanja mehaničkih oštećenja cevovoda, poklapaju sa rezultatima iz literature sa kojima su poređeni. Ostali pikovi predviđeni primenom programa RELAP (Slika 9.21b) su izazvani akustičkim prostiranjem talasa pritiska u vertikalnoj cevi između bazena ispunjenog vodom i zatvorenog kraja cevi, dok razvijeni program (Slika 9.21a) predviđa samo četiri udara izazvana kretanjem vodenog stuba u cevi. 144 a) b) Slika 9.21 Promena pritiska u toku vremena u blizini zatvorenog kraja vertikalne cevi aparature za simulaciju vodenog topa. Rezultati dobijeni primenom razvijenog programa (gore), rezultati dobijeni primenom programa RELAP5/MOD3 (Yeung i dr., 1993) (dole). Kondenzacija pare na pothlađenoj vodi počinje u 0 s. 145 Na Slici 9.22 prikazano je pomeranje čela vodenog stuba unutar vertikalne cevi od trenutka kad počne kondenzacija do trenutka potpunog kondenzovanja pare unutar cevi. Posle prvog udara vodenog stuba u 0,11 s talas pritiska se odbija od kraja cevi i obara pritisak, krećući se naniže ka bazenu. Obaranjem pritiska ispod vrednosti pritiska zasićenja, određenog temperaturom vode, voda isparava i nastaje dvofazna mešavina. Zapremina pare približne dužine od 0,15 m se formira od vrha cevi. Vodeni stub ponovo kreće naviše i opet udara o zatvoreni kraj (oko 0,3 s), ovaj put slabijeg intenziteta. Vreme između nastanka dva pika pritiska predstavlja vreme između dva udara vodenog stuba o zatvoreni kraj. Drugi i treći udar vodenog stuba o zatvoreni kraj cevi i odbijanje talasa pritiska izazivaju manje isparavanje vode i mali četvrti pik pritiska posle kog sledi potpuna kondenzacija pare, nema više pare u vertikalnoj cevi. Slika 9.22 Pomeranje položaja razdelne površine u toku vremena u vertikalnoj cevi u aparaturi za simulaciju vodenog topa (Slika 9.20). Rezultati dobijeni primenom razvijenog programa. Slika 9.23 predstavlja promenu brzine u toku vremena čela vodenog stuba. Brzina raste do 8,4 m/s, zatim polako opada i u trenutku udara vodenog stuba o zatvoreni kraj vertikalne cevi u 0,11 s menja svoj znak. Ovaj proces se dešava za veoma kratko vreme, pa otuda i naziv vodeni top. Prvi pik pritiska ima maksimalnu vrednost i 146 najopasniji je jer može prouzrokovati teška oštećenja opreme. Predstavljeni rezultati su dobijeni sa sa vrednošću ,0 2Ka  m 2 /m 3 (jednačina 7.15), što odgovara malom broju uključenih kapi sa čela vodenog stuba i maloj specifičnoj razdelnoj površini između vodene pare i kapi, mnogo manjoj nego što je kod eksperimenata hidrauličkog udara usled kondenzacije pare na stubu pothlađene vode, poglavlje 9.2. Slika 9.23 Promena brzine čela vodenog stuba u toku vremena u vertikalnoj cevi u aparaturi za simulaciju vodenog topa (Slika 9.20) dobijena primenom razvijenog programa. 9.4 Analitičko odreĎivanje proizvoda koeficijenta prelaza toplote i razdelne površine pri kondenzaciji pare na pothlaĎenoj vodi pri hidrauličkom udaru izazvanom kondenzacijom pare Hidraulički udar izazvan intenzivnom kondenzacijom pare nastaje usled neposrednog kontakta pothlađene tečnosti, koja je na nižoj temperaturi, i pare na višoj temperaturi. Formiranje razdelne površine, odnosno kontaktne površine pare i tečnosti, je stohastičke prirode. Usled intenzivne kondenzacije pare razdelna površina se naglo povećava. Što je veća razdelna površina, biće veća brzina kondenzacije. Određivanje razdelne površine je veoma složeno, jer se samo na početku procesa može pretpostaviti 147 da je razdelna površina jednaka površini poprečnog preseka strujnog kanala u kome se dešava hidraulički udar, a nakon toga dolazi do otkidanja nestacionarnih mlazeva tečnosti i kapi sa čela stuba tečnosti u parni prostor. Geometrija mlazeva tečnosti i otkinutih kapi, kao i broj kapi, se menjaju za vreme hidrauličkog udara. Ove podatke je analitičkim putem nemoguće odrediti. U Tabeli 9.4 je dat pregled početnih vrednosti i nekih ulaznih podataka za proračun kondenzacije na stubu tečnosti za tri slučaja, na koje je primenjen razvijeni program za proračun kondenzacije pare na stubu pothlađene tečnosti. Simulacije su sprovedene u dva geometrijski različita strujna prostora:  u eksperimentalnoj aparaturi za ispitivanje hidrauličkog udara (Slika 9.8), za dva nivoa pritiska u sistemu, 5,51 bar (slučaj 1 - Tabela 9.4) i 6,54 bar (slučaj 2 - Tabela 9.4) i  u eksperimentalnoj aparaturi za simulaciju vodenog topa (Slika 9.20) (slučaj 3 - Tabela 9.4). AC(m 2) V2,0 (m 3) L ( m ) p1,0 (bar) T1,0 ( 0C) p2,0 (bar) T2,0 ( 0C) p2,0 ≈ p1,0 T2,0 = Tsat(p2,0) Slika 9.24 Stanje u vertikalnoj cevi, ispunjenoj parom, na početku eksperimenta u trenutku kada pothlađena voda počinje da utiče u cev. 148 Tabela 9.4 Početne vrednosti za proračun kondenzacije na stubu tečnosti i određivanje koeficijenta prelaza toplote i razdelne površine Slučaj 1 - 5,51 bar Slučaj 2 - 6,54 bar Slučaj 3 - vodeni top Geometrija (C - cev) L (m) 5,000 5,000 0,711 dC (m) 0,092 0,092 0,038 x (m) 0,100 0,100 0,024 AC (m 2 ) 6,662 310 6,662 310 1,140 310 Fizičke osobine vode i pare (0 – početni vremenski trenutak 0 s, 1 – voda, 2 - para) p1 (bar) 5,510 6,540 1,023 t1 ( 0 C) 22,0 23,0 49,0 1 ( 3kg/m ) 997,977 997,791 988,497 1 (  W/ mK ) 0,602 0,604 0,642 1 (  kg/ ms ) 9,543 410 9,320 410 5,562 410 cp1 (  J/ kgK ) 4182,0 4181,0 4179,0  ( N/m ) 4,761 210 4,614 210 5,892 210 Pr1 (-) 7,013 6,717 3,658 u1 ( m/s ) 4,659 5,322 3,895 h1 ( J/kg ) 656179,0 97079,0 205234,0 p2 (bar) 5,510 6,540 1,023 t2 ( 0 C) 155,000 162,000 100,000 2 ( 3kg/m ) 2,924 3,438 0,603 2 (  W/ mK ) 0,032 0,033 0,025 2 (  kg/ ms ) 1,418 510 1,441 510 1,228 510 cp2 (  J/ kgK ) 2447,0 2515,0 2078,0 Pr2 (-) 1,023 1,030 0,999 u2 ( m/s ) 4,659 5,322 3,895 h2 ( J/kg ) 2752411,0 2759858,0 2675955,0 r ( J/kg ) 2096232,0 2662779,0 2470721,0 p2,0 (bar) 3,824 4,761 1,023 t2,0 ( 0 C) 142,0 150,0 100,0 2,0 ( 3kg/m ) 2,073 2,548 0,603 Ukupna zapremina (9.21) i ukupna masa pare (9.22) u cevi na početku eksperimenta V2,0 ( 3m ) 3,331 210 3,331 210 8,108 410 m2,0 ( kg ) 6,905 210 8,487 210 4,889 410 U Tabeli 9.4 je nulom obeležen početni vremenski trenutak, brojem 1 je označeno stanje tečne faze, a brojem 2 stanje parne faze. 149 Za prva dva slučaja, kada do hidrauličkog udara dolazi pri direktnom kontaktu parne i tečne faze, ostvarenom otvaranjem brzodelujućeg ventila, neposredno posle otvaranja ventila, (Slika 9.24) smatramo da se pritisak u parnom prostoru veoma brzo izjednači sa pritiskom tečnosti u rezervoaru i u horizonatalnoj cevi (p2,0  p1,0 = 5,51 bar za prvi slučaj odnosno 6,54 bar za drugi slučaj). Na tom pritisku dolazi do kondenzacije pare. U trećem slučaju je pritisak pare i tečnosti jednak p2,0 = p1,0 = 1,023 bar. Ukupna zapremina pare u vertikalnoj cevi pre otvaranja ventila, odnosno početka eksperimenta je 2 2,0 4 C C d V A L L    ( 3m ), (9.21) gde je dC prečnik vertikalne cevi (m), L visina vertikalne cevi (m), a ukupna masa pare u cevi je 2,0 2,0 2,0m V  ( kg ). (9.22) Za brzinu tečne faze u1 i parne faze u2 je usvojena srednja vrednost brzine čela vodenog stuba, dobijena kao rezultat numeričkog proračuna primenom razvijenog kompjuterskog programa, od početka kondenzacije u 0 s do trenutka kada se sva para koja se nalazila u vertikalnoj cevi na početku eksperimenta iskondenzovala  (s). Ovaj trenutak se poklapa sa trenutkom nastanka prvog pika pritiska, odnosno prvim udarom stuba tečnosti o zatvoreni kraj cevi. Neposredno ispred čela vodenog stuba, posmatrano u smeru kretanja, imamo parnu fazu, a iza čela vodenog stuba imamo tečnu fazu. Vremenski trenutak  je određen na osnovu dobijenih numeričkih rezultata (Dijagrami 9.12, 9.19 i 9.21). Analitičkim putem se proizvod koeficijenta prelaza toplote i specifične razdelne površine pri kondenzaciji pare na stubu pothlađene tečnosti može odrediti polazeći od jednakosti 2,0 c m r r xA     ( W ), (9.23) 150 gde je količina toplote koja se oslobodi pri kondenzaciji ukupne mase pare u vertikalnoj cevi za vreme , jednaka količini toplote koja se preda pri kondenzaciji, u jedinici vremena, zapremini jednog proračunskog elementa, čiji je prostorni korak x . Uvrštavanjem izraza za brzinu kondenzacije c ic q a r   ( 3 kg m s ), (9.24) gde je toplotni fluks koji se pri kondenzaciji na čelu stuba tečnosti sa pare preda tečnoj fazi  2,0 1,0c cq h T T  ( 2 W m ), (9.25) iz jednačine (9.23) može da se izrazi   2,0 2 2,0 1,0 4 c i C m r h a x d T T         ( 3 W m K ). (9.26) Numerički dobijene vrednosti proizvoda koeficijenta prelaza toplote pri kondenzaciji pare na stubu pothlađene tečnosti slede iz zavisnosti primenjene u razvijenom kompjuterskom programu 3 i c i num a Du h a b x Dt    ( 3 W m K ). (9.27) Empirijski koeficijenti a i b u jednačini (9.27) su određeni poređenjem rezultata dobijenih numeričkim proračunom sa eksperimentalnim rezultatima. Za vrednost ubrzanja razdelne površine ( i Du Dt ) uzeta je apsolutna vrednost srednje vrednosti ubrzanja čela vodenog stuba, dobijena kao rezultat numeričkog proračuna primenom razvijenog kompjuterskog programa, od početka eksperimenta do udara stuba tečnosti o zatvoreni kraj cevi (od 0 s - ). Poređenje vrednosti proizvoda c ih a sračunatih analitičkim i numeričkim putem je dato u Tabeli 9.5. U Tabeli 9.5 vrednost c i numh a predstavlja srednju vrednost numerički dobijenih rezultata od početka eksperimenta do udara stuba tečnosti o zatvoreni kraj cevi (od 0 s - ). 151 Tabela 9.5 Proizvod koeficijenta prelaza toplote i specifične razdelne površine pri kondenzaciji pare na pothlađenoj vodi pri hidrauličkom udaru izazvanom kondenzacijom pare. Poređenje analitički i numerički dobijenih vrednosti. Slučaj 1 - 5,51 bar Slučaj 2 - 6,54 bar Slučaj 3 - vodeni top  (s) 1,284 1,220 0,103 hcai (W/(m 3 K)) 1,272 610 2,001 610 8,543 610 a 4849940 5238061 457555 b 4250 2982 6518 D D iu t (m/s 2 ) 10,609 12,462 115,825 ich a num (W/(m 3 K)) 5,498 710 6,978 710 7,464 1010 Stohastička priroda formiranja razdelne površine u posmatranim eksperimentima je uzrok zbog koga se javlja razlika u redu veličine (Tabela 9.5) proizvoda koeficijenta prelaza toplote i specifične razdelne površine pri kondenzaciji pare na stubu pothlađene tečnosti dobijenih analitičkim i numeričkim putem. 9.4.1 OdreĎivanje koeficijenta prelaza toplote i razdelne površine pri kondenzaciji pare na čelu stuba tečnosti Za određivanje koeficijenta prelaza toplote pri kondenzaciji pare na čelu stuba tečnosti koristi se Dittus-Boelter-ova korelacija za određivanje Nusselt-ovog broja pri prinudnoj konvekciji, pri turbulentnom strujanju u cevi, koja glasi 0.8 0.4 1Nu 0,023Re PrDB C , (9.28) gde je Reynolds-ov broj za strujanje tečnosti kroz cev prečnika dC, brzinom u1 jednak 1 1 1 Re CC u d    (-). (9.29) Koeficijent prelaza toplote pri kondenzaciji na čelu stuba tečnosti jednak je 1NuDB cST C h d   ( 2 W m K ), (9.30) gde je 1 toplotna provodljivost tečne faze (W/mK). Toplotni fluks koji se pri kondenzaciji na čelu stuba tečnosti sa pare preda tečnoj fazi kroz razdelnu površinu jednak je 152  2 1c ST cSTq h T T  ( 2 W m ), (9.31) gde su T1 i T2 temperature tečne i parne faze sledstveno (K). Ukupna površina za kondenzaciju pare na čelu stuba tečnosti je jednaka u idealnom slučaju površini poprečnog preseka cevi 2 4 C i ST C d A A    ( 2m ). (9.32) U Tabeli 9.6 su navedene sračunate vrednosti koeficijenta prelaza toplote i razdelne površine pri kondenzaciji na čelu stuba tečnosti. Tabela 9.6 Prelaz toplote pri kondenzaciji pare na čelu stuba tečnosti (DB - Dittus- Boelter-ova korelacija, c – kondenzacija, ST – čelo stuba tečnosti) Slučaj 1 - 5,51 bar Slučaj 2 - 6,54 bar Slučaj 3 - vodeni top ReC (-) 448751,748 524718,214 279291,124 NuDB (-) 1666,074 1855,833 878,781 hc,ST (W/(m 2 K)) 10893,809 12170,698 14817,258 qc,ST (W/m 2 ) 1448876,549 1691726,984 755680,150 Ai ST (m 2 ) 6,662 310 6,662 310 1,140 310 9.4.2 Prelaz toplote pri kondenzaciji pare na kapima pothlaĎene tečnosti otkinutim sa čela stuba tečnosti i uključenim u parni prostor Naglo ubrzanje stuba tečnosti, uzrokovano intenzivnom kondenzacijom pare na pothlađenoj tečnosti, uzrokuje otkidanje kapi sa čela stuba tečnosti. Posledica većeg ubrzanja stuba tečnosti je veći broj otkinutih kapi. Na taj način se povećava razdelna površina, na kojoj se dešava fazni prelaz, pa je veća i brzina kondenzacije. Posmatramo sfernu kap pothlađene tečnosti unutrašnjeg prečnika d0 i uniformne temperature T0, koja dolazi u direktan kontakt sa parom u stanju zasićenja uniformnog pritiska p, i temperature T, pri čemu da bi došlo do kondenzacije mora da bude   0satT T p T  . Pretpostavlja se da je kap sferno simetrična i da nema unutrašnjeg 153 strujanja u kapi. Usled kondenzacije kap raste. Za uobičajene fluide promena prečnika kapi, koja nastaje zbog kondenzacije, je veoma mala pa se usvaja pretpostavka da je prečnik kapi konstantan tokom kondenzacije 0 2 const.d d r   Na taj način se može dobiti analitičko rešenje, koje se svodi na provođenje toplote u krutoj sferi 21 1 1 2 1T T r t r r r             . (9.33) Početni i granični uslovi za ovakav problem su 0 1 1 0 za 0 i 0 za =0 za / 2sat T t T T r T r d r       . (9.34) Kronig i Brink (1950) (Ghiaasiaan, 2008) su rešili problem nestacionarne difuzije unutar sfernog fluidnog delića sa unutrašnjom cirkulacijom. Prema njihovom rešenju kada t  , što predstavlja rešenje blisko ravnotežnom rešenju, sledi Nusselt- ov broj za kondenzaciju na pothlađenim tečnim kapima Nu 17,9.K  (9.35) Koeficijent prelaza toplote pri kondenzaciji na pothlađenim tečnim kapima otkinutim sa čela stuba tečnosti i uključenim u parni prostor je 1 NuK cK K h d   ( 2 W m K ), (9.36) pri čemu je prečnik tečnih kapi jednak 2 2 1 Wecr Kd u    (m), (9.37) gde je  površinski napon (N/m), We 0.799cr  kritični Weber-ov broj prema (Saito i dr. 1978) (-), a 2 gustina pare (kg/m 3 ). Toplotni fluks koji se pri kondenzaciji na kapi pothlađene tečnosti sa pare preda tečnoj fazi preko razdelne površine jednak je  2 1c K cKq h T T  ( 2W/m ), (9.38) gde su T1 i T2 temperature tečne i parne faze sledstveno. 154 Ukupna površina za kondenzaciju pare na kapi tečnosti je jednaka proizvodu ukupnog broja kapi i površine jedne sferne kapi i K K KA N A ( 2m ), (9.39) gde je NK ukupan broj kapi (-), a 2 4 2 K K d A         površina jedne sferne kapi (m2) . Ukupan broj kapi se određuje na osnovu ukupne razdelne površine pri kondenzaciji na kapima, dobijenoj iz energetskog bilansa, koji glasi: masa pare koja se kondenzuje u jedinici vremena jednaka je količniku ukupne količine toplote koja se odvede od pare preko stuba tečnosti i svih kapi u jedinici vremena i latentne toplote  2,0 1 cST i ST c K i K m q A q A r   ( kg s ), (9.40) odavde sledi da je ukupna površina za kondenzaciju pare na kapi tečnosti jednaka 2,0 c ST i ST i K c K m r q A A q    ( 2m ). (9.41) Poznavanjem ukupne razdelne površine može se konačno odrediti broj uključenih kapi NK. U Tabeli 9.7 su navedeni rezultati proračuna određivanja koeficijenta prelaza toplote i razdelne površine pri kondenzaciji na pothlađenim kapima otkinutim sa čela stuba tečnosti i uključenim u parni prostor. Tabela 9.7 Prelaz toplote pri kondenzaciji na pothlađenim kapima vode otkinutim sa čela stuba tečnosti i uključenim u parni prostor. (K – kap, Wecr (Saito 1978, Sami 1988), hc, K (Ghiaasiaan (Kronig i Brink 1950)) Slučaj 1 - 5,51 bar Slučaj 2 - 6,54 bar Slučaj 3 - vodeni top Wecr (-) 0,799 0,799 0,799 dK (m) 5,994 410 3,787 410 4,590 310 hc,K (W/(m 2 K)) 17985,228 28551,764 2505,363 qc,K (W/m 2 ) 2392035,268 3968695,218 127773,497 Ai1 K (m 2 ) 1,129 610 4,505 710 6,618 510 Ai K (m 2 ) 4,310 210 4,386 210 8,540 210 N K (-) 38187 97356 1290 Odnos površina svih kapi i čela vodenog stuba Ai K/Ai ST (m 2 ) 6,469 6,583 74,903 155 9.4.3 Kondenzacija na pothlaĎenom mlazu tečnosti otkinutom sa čela stuba tečnosti Za analitičko definisanje kondenzacije na pothlađenom mlazu tečnosti potrebno je poznavati geometriju formiranih mlazeva. Formiranje mlazeva zavisi od oblika razdelne površine tečne i parne faze i intenziteta kondenzacije pare, čijim se povećanjem razdelna površina naglo povećava. Geometrija mlazeva tečnosti i njihov broj se menjaju za vreme hidrauličkog udara, pa je stvarne dimenzije mlaza veoma teško odrediti. Kako bi se sproveo proračun kondenzacije pare u stanju zasićenja na pothlađenom mlazu tečnosti u ovom poglavlju sprovedene su dve grupe proračuna, koje se razlikuju prema usvojenoj pretpostavci o početnom prečniku mlaza (D0). U prvoj grupi je pretpostavljeno da je početni prečnik mlaza jednak polovini talasne dužine najnestabilnijeg talasa, koji nastane pri formiranju mlaza, a u drugoj grupi je početni prečnik mlaza jednak prečniku cevi u kojoj dolazi do kondenzacije pare na pothlađenoj tečnosti. Takođe je usvojeno da je dužina mlaza jednaka početnom prečniku mlaza (L/D0=1), kao i da su površina poprečnog preseka mlaza na ulazu i ukupna površina mlaza međusobno jednake (A0/AS=1). Uzeto je da je odnos faktora trenja hrapave i glatke mlaznice jednak jedinici (F=1). U Tabeli 9.8 su navedene neke korelacije za kondenzaciju pare u stanju zasićenja na mlazu pothlađene tečnosti (Ghiaasiaan, 2008). U ovim korelacijama se javlja Stanton-ov broj koji je jednak 1 1 1 10 St ,i p h c u  (9.42) gde je 1ih koeficijent prelaza toplote na razdelnoj površini sa strane tečne faze (W/(m 2 K)), 1pc specifični toplotni kapacitet tečne faze pri konstantnom pritisku (J/(kgK)), a 10u brzina mlaza na ulazu (m/s). Weber-ov broj za paru u stanju zasićenja je 2 2 2 0 2We , u D   (9.43) 156 a za tečni mlaz 2 1 10 0 10We . u D   (9.44) Jakob-ov broj za tečnu fazu se određuje primenom izraza  1 2 10 1Ja , pc T T r   (9.45) gde je 2T temperatura pare u stanju zasićenja (K), 10T temperatura mlaza tečnosti na ulazu (K), a r latentna toplota kondenzacije (J/kg). Suratman-ov broj jednak je 1 0 2 1 Su . D    (9.46) Tabela 9.8 Empirijske korelacije za kondenzaciju pare u stanju zasićenja na mlazu pothlađene tečnosti (Ghiaasiaan, 2008). Isackenko i dr.(1971) 0,42 0,17 0,09 0,13 0,35 10 1 1 2 0 St 0,0335 Re Pr Ja We L D          gde je L dužina mlaza  m , 0D početni prečnik mlaza  m . Benedek (1976) 0,06 0,0840 1St 0,00286 Ja S A A       gde je 0A površina poprečnog preseka mlaza na ulazu  2m , S A ukupna površina mlaza  2m . Sklover i Rodivilin (1975) 0,6 0,4 0,55 0,11 0,4 10 1 1 2 0 St 2,7 Re Pr Ja We L D          De Salve i dr.(1986) 0,52 0,38 0,52 0,19 10 1 1 0 St 3,25 Re Pr Ja L D         Kim i Mills (1989) 0,57 0,2 0,7 0,19 0,18 10 1 0 St 3,2 Re Pr Su F L D           gde je F odnos faktora trenja hrapave i glatke mlaznice Sam i Patel (1984) 0,35 10 2 0,35 0,3 10 10 2 St 0,075Re za 1< We 13 St 0,075Re We za We 1       157 Ako se usvoji da je početni prečnik mlaza jednak maksimalnom prečniku mlaza, koji je jednak polovini talasne dužine najnestabilnijeg talasa 0 max / 2TD d   (Whalley, 1996), gde je talasna dužina najnestabilnijeg talasa jednaka   1/2 1 2 T C g            (m), (9.47) pri čemu je 2 / 6C  za dvodimenzionalne talase, dobija se vrednost 0 max 0,017 m.D d  Primenom različitih empirijskih korelacija uz pretpostavke 0 1 L D  , 0 1 S A A  i F=1, dobijaju se analitičke vrednosti za Stanton-ov broj i koeficijent prelaza toplote pri kondenzaciji pare, na pothlađenom mlazu tečnosti, na razdelnoj površini sa strane tečne faze, predstavljene u Tabeli 9.9. Sračunate vrednosti bezdimenzionih brojeva su: We2 = 22,757, We10 = 7767,170, Ja1 = 0,265, Su = 888568,745, Re10 = 83076,256. Tabela 9.9 Vrednosti Stanton-ovog broja i koeficijenta prelaza toplote pri kondenzaciji pare na pothlađenom mlazu tečnosti, na razdelnoj površini sa strane tečne faze, kada se usvoji da je maksimalni prečnik mlaza jednak polovini talasne dužine najnestabilnijeg talasa 0 max / 2.TD d   Korelacije St 1ih (W/(m 2 K)) Isackenko i dr. (1971) 0,42 0,17 0,09 0,13 0,35 10 1 1 2 0 St 0,0335 Re Pr Ja We L D          0,014538 283410,00 Benedek (1976) 0,06 0,0840 1St 0,00286 Ja S A A       0,003197 62327,77 Sklover i Rodivilin (1975) 0,6 0,4 0,55 0,11 0,4 10 1 1 2 0 St 2,7 Re Pr Ja We L D          0,040231 784278,98 De Salve i dr. (1986) 0,52 0,38 0,52 0,19 10 1 1 0 St 3,25 Re Pr Ja L D         0,012392 241572,71 Kim i Mills (1989) 0,57 0,2 0,7 0,19 0,18 10 1 0 St 3,2 Re Pr Su F L D           0,006293 122683,18 Sam i Patel (1984) 0,35 10 2 0,35 0,3 10 0 2 St 0,075Re za 1< We 13 St 0,075Re We za We 1L       0,001423 27743,34 0,020909 407603,22 158 Kada se usvoji da je maksimalni prečnik mlaza jednak prečniku cevi 0 max 0.0921 mD d  , za sračunate vrednosti bezdimenzionih brojeva: We2 = 123,244, We10 = 42063,867, Ja1 = 0,265, Su = 4812130,814, Re10 = 449907,580, dobijaju se analitičke vrednosti Stanton-ovog broja i koeficijenta prelaza toplote pri kondenzaciji pare, na pothlađenom mlazu tečnosti, na razdelnoj površini sa strane tečne faze date u Tabeli 9.10. Takođe su usvojene pretpostavke 0 1 L D  , 0 1 S A A  i F=1. Tabela 9.10 Vrednosti Stanton-ovog broja i koeficijenta prelaza toplote pri kondenzaciji na pothlađenom mlazu tečnosti, na razdelnoj površini sa strane tečne faze, kada se usvoji da je maksimalni prečnik mlaza jednak prečniku cevi 0 max 0.0921 m.D d  Korelacije St 1ih (W/(m 2 K)) Isackenko i dr.(1971) 0,42 0,17 0,09 0,13 0,35 10 1 1 2 0 St 0,0335 Re Pr Ja We L D          0,019704 384124,06 Benedek (1976) 0,06 0,0840 1St 0,00286 Ja S A A       0,003197 62327,77 Sklover i Rodivilin (1975) 0,6 0,4 0,55 0,11 0,4 10 1 1 2 0 St 2,7 Re Pr Ja We L D          0,040231 784278,98 De Salve i dr.(1986) 0,52 0,38 0,52 0,19 10 1 1 0 St 3,25 Re Pr Ja L D         0,006522 127134,15 Kim i Mills (1989) 0,57 0,2 0,7 0,19 0,18 10 1 0 St 3,2 Re Pr Su F L D           0,003256 63483,80 Sam i Patel (1984) 0,35 10 2 0,35 0,3 10 10 2 St 0,075Re za 1< We 13 St 0,075Re We za We 1       0,000788 15359,69 0,019215 374589,23 Primenom različitih korelacija za određivanje Stanton-ovog broja, a time i koeficijenta prelaza toplote toplote pri kondenzaciji na pothlađenom mlazu tečnosti, na razdelnoj površini sa strane tečne faze, dobijaju se vrednosti koje su u rasponu od 15.359,69 - 784.278,98 W/(m 2 K). Zbog ovih velikih odstupanja je umesto analitičke zavisnosti u razvijenom programu primenjena korelacija za određivanje proizvoda koeficijenta prelaza toplote pri kondenzaciji i specifične razdelne površine u funkciji ubrzanja razdelne površine, pri čemu su empirijski koeficijenti određeni poređenjem 159 rezultata dobijenih numeričkim proračunom sa eksperimentalnim rezultatima. Ovom korelacijom je obuhvaćena i kondenzacija na mlazevima tečnosti i na kapima otkinutim sa čela stuba tečnosti. 160 10. ZAKLJUČAK U termohidrauličkim analizama prelaznih režima rada ili u okviru sigurnosnih analiza termoenergetskih postrojenja potrebno je predvideti prostiranje temperaturskih talasa u jednofaznom ili dvofaznom toku. U okviru ove disertacije je razvijen numerički postupak za proračun prostiranja temperaturskih talasa i razdelne površine između stuba tečnosti i parne faze u cevovodima i složenim cevnim mrežama. Predviđanje prostiranja temperaturskih talasa je značajno za optimizaciju rada prelaznih režima u sistemu daljinskog grejanja i analize udesnih stanja u termoenergetskim postrojenjima poput intenzivne kondenzacije sa hidrauličkim udarom. Razvijeni postupak je primenjen za termohidrauličku analizu prelaznih režima rada u sistemu daljinskog grejanja nastalih usled promene snage toplotnog izvora i predviđanje uslova nastanka hidrauličkog udara u dvofaznom toku izazvanog intenzivnom kondenzacijom pare. Energijska jednačina jednodimenzionalnog strujanja homogenog fluida se rešava numerički duž karakterističnog pravca određenog kretanjem fluidnog delića u vremensko-prostornom koordinatnom sistemu. Na taj način se parcijalna diferencijalna jednačina, koja predstavlja bilans energije, prevodi u običnu diferencijalnu jednačinu, u kojoj se zatim totalni diferencijali aproksimiraju konačnim razlikama prvog reda i dobija se diferencna jednačina čijim se rešavanjem po zavisno promenljivoj, uz odgovarajuće početne i granične uslove, dobija njena vrednost duž strujnog prostora tokom vremena. Rešenje se nalazi u okviru numeričke mreže sa konstantnim korakom između čvorova tako da se početak karakterističnog pravca u opštem slučaju nalazi između čvorova mreže, a završetak je u najbližem susednom čvoru u pravcu strujanja. Zavisno promenljiva u energijskoj jednačini je temperatura ili entalpija. Rešavanjem energijske jednačine sa entalpijom kao zavisno promenljivom ostvaruje se praćenje razdelne površine između tečne i parne faze, pri čemu se prisutstvo faze određuje na osnovu vrednosti termodinamičkog stepena suvoće. U cilju eliminisanja numeričke difuzije početna vrednost temperature ili entalpije se određuje primenom Lagrange- ovog interpolacionog polinoma trećeg stepena na osnovu poznatih vrednosti, koje su 161 zadate kao početni uslovi ili su rezultat prethodnog koraka integracije. Na taj način je dobijen numerički postupak četvrtog reda tačnosti u prostoru. Validacija razvijenog numeričkog postupka za uslove prostiranja temperaturskih talasa u složenoj cevnoj mreži izvršena je poređenjem rezultata numeričke simulacije prelaznih procesa u odabranom sistemu daljinskog grejanja sa izmerenim vrednostima. Poređenje sračunatih i izmerenih vrednosti prostiranja temperaturskih talasa od toplane do tri podstanice u različitim delovima mreže topovoda i na različitim rastojanjima od izvora toplote pokazuje zadovoljavajuće slaganje. Oblik i amplitude sračunatih talasa su očuvani i ne dolazi do deformacije talasa usled numeričke difuzije. Dobijeni rezultati omogućuju određivanje vremena potrebnog da temperaturski talas pređe rastojanje od izvora toplote, tj. toplane do najudalјenijeg potrošača, odnosno vremensko kašnjenje promene snage potrošnje u odnosu na promenu snage izvora toplote. Na osnovu ove vremenske konstante sagledava se vreme kada treba menjati snagu izvora da bi se pravovremeno zadovolјile potrebe potrošača. Obavljena istraživanja su deo verifikacije metodologije koja su podloga za optimizaciju režima rada sistema daljinskog grejanja u cilju povećanja energetske efikasnosti snabdevanja potrošača toplotom. Za potrebe predviđanja hidrauličkog udara izazvanog kondenzacijom pare bilo je potrebno termohidraulički model za rešavanje jednodimenzionalnog strujanja homogenog fluida proširiti sa modelom za praćenje kretanja razdelne površine, modelom kondenzacije i modelom nestacionarnog trenja. Za predviđanje položaja razdelne površine primenjen je Lagrange-ov interpolacioni polinom 3. stepena umesto linearne interpolacije koja se standardno koristi. Na taj način se postiže tačnije predviđanje kretanja razdelne površine stuba tečnosti i pare, koje je određeno duž karakterističnih pravaca kretanja fluidnih delića u prostorno-vremenskom koordinatnom sistemu, pri čemu je prisustvo tečne ili parne faze određeno na osnovu vrednosti entalpija faza. Značajan uticaj na prostiranje talasa pritiska, koji je nastao usled pojave hidrauličkog udara izazvanog kondenzacijom pare, ima pojava nestacionarnog trenja i vrednost brzine, odnosno intenziteta kondenzacije. Pri hidrauličkom udaru se javljaju skokovite promene pritiska velikog intenziteta koje dovode do naglih promena brzine. Uticaj nestacionarnog trenja se oseća usled promene brzine, pa je za postizanje 162 odgovarajuće tačnosti pri modeliranju hidrauličkog udara uzimanje u obzir samo pojave stacionarnog trenja nedovoljno. Intenzitet kondenzacije je određen odgovarajućim neravnotežnim modelom na osnovu razlika temperatura tečne i parne faze u okolini razdelne površine. Simulacija procesa kondenzacije je složena zbog nemogućnosti da se tačno definiše oblik razdelne površine, koji je nepravilan i nestalan usled odvajanja nestacionarnih mlazeva tečnosti i kapi sa čela stuba tečnosti. Pored toga, potrebno je odrediti i koeficijent prelaza toplote usled kondenzacije. Stoga je uticaj razdelne površine i koeficijenta prelaza toplote pri kondenzaciji uzet u obzir preko nove korelacije kojom je definisan proizvod ove dve veličine kao funkcija ubrzanja razdelne površine. Empirijski koeficijenti u ovoj korelaciji su određeni poređenjem rezultata dobijenih numeričkim proračunom sa eksperimentalnim rezultatima. Validacija rezultata dobijenih primenom razvijenog kompjuterskog programa na simulaciju kretanja razdelne površine i hidrauličkog udara izazvanog kondenzacijom pare u različitim geometrijama strujnih kanala i za različite nivoe pritisaka u sistemu izvršena je poređenjem sa analitičkim rešenjem kretanja razdelne površine, raspoloživim izmerenim vrednostima i rezultatima drugih kompjuterskih programa. Postignuto je zadovoljavajuće slaganje rezultata. Dobijeni rezultati pokazuju da je sprovedeni proračun stabilan i da razvijeni kompjuterski program omogućava pouzdano predviđanje uslova nastanka hidrauličkog udara izazvanog intenzivnom kondenzacijom pare, kao i dinamičke promene pritiska izazvane udarom. 163 LITERATURA 1. Barna, I.F., Inre, A.R., Baranyai, G., Ezsol, G., Experimental and Theoretical Study of Steam Condensation Induced Water Hammer Phenomena, Nuclear Engineering and Design, 240, 2010, str. 146-150. 2. Benonysson, A., Bohm, B., Ravn, H.F., Operational Optimization in a District Heating System, Energy Conversion and Management, Vol. 36, No. 5, 1995, str. 297-314. 3. Beuthe, T.G., Review of Two-Phase Water Hammer, Proceedings of the 18th Annual Canadian Nuclear Society Conference, Toronto, Canada, 1997. 4. Bjorge, R.W., Griffith, P., Initiation of Water Hammer in Horizontal and Nearly Horizontal Pipes Containing Steam and Subcooled Water, ASME J. Heat Trans., 106, 1984, str. 835-840. 5. Blömeling, F., Neuhaus, T., Schaffrath, A., Development and Validation of DYVRO for the Simulation of Condensation Induced Water Hammer: Comparison with a Two-Phase Slug Model and with Experimental Data of the Water Cannon Experiment, The 15 th International Topical Meeting on Nuclear Reactor Thermal - Hydraulics (NURETH-15), Pisa, Italy, May 12-17, 2013, NURETH15-307 6. Ceuca, C.S., Macian-Juan, R., Benchmark of Surface Renewal Theory Based Heat Transfer Coefficients for the Simulations of Direct Contact Condensation In Pipes Using The 1D and 3D Approach, The 15 th International Topical Meeting on Nuclear Reactor Thermal - Hydraulics (NURETH-15), Pisa, Italy, May 12-17, 2013, NURETH15-441 7. Choi, D.K., PIPES – A Computer Code for Analysis of Dynamic Hydraulic Response in Plant Piping Systems, Windsor, 1983. 8. Chun, M.H., Yu, S.O., A Parametric Study and a Guide Chart to Avoid Condensation/Induced Water Hammer in a Horizontal Pipe, Nuclear Engineering and Design, 201, 2000, str. 239-257. 9. Demidovitch, B., Maron, I., Elements de calcul numerique, Editions Mir, Moscou, 1987. 10. De Vries, M., Simon, A., Suction Effects on Feedpump Performance: A Literature Survey, Sulzer Brothers Limited, EPRI Report CS-4204, Palo Alto, USA, 1985. 11. Dittus, F.W., Boelter, L.M.K., Heat Transfer in Automobile Radiators of the Turbulent Type. University of California Publications, Vol. 2, 1930, str. 443-461. 164 12. Fox, J.A., Hydraulic Analysis of Unsteady Flow in Pipe Networks, Macmillan Press Ltd., First Ed., 1977. 13. Gabrielaitiene, I., Bohm, B., Sunden, B., Modelling temperature dynamics of a district heating system in Naestved, Denmark – A case study, Energy Conversion & Management, Vol. 48, 2007, str. 78-86. 14. Ghiaasiaan, S.M., Two-phase flow, boiling, and condensation. Cambridge University Press, 2008, str. 482. 15. Joukowsky, N., Water Hammer. Transleted by O. Simon, Proc. American Water Works Association, 1904, str. 24. 16. Jovanović, M., Jocić, Lj., Instalacija i postupak za sprečavanje hidrauličkog udara u postrojenju pumpe za napajanje vodom generatora pare, Patent br. 47840, Zavod za intelektualnu svojinu, Beograd, 1996. 17. Jun, L., Oscillating Manometer, Numerical Benchmark Test No. 2.2, DOE/EPRI Workshop on Two-Phase flow Fundamentals, Rensselear Polytechnic Institute, Troy, NY USA, 1987, str. 52-57. 18. Kirsner, W., Waterhammer, HPAC Heating/Piping/Air/Conditioning, January 1999, str. 113-122. 19. Kronig, R., Brink, J.C., On the theory of extraction from falling droplets. Applied Science Research, A2, 1950, str. 142-154. 20. Kunz, V.J., Haldi, P.A., Sarlos, G., Dynamic behaviour of district heating networks, Fernwaerme International, Vol. 20, 1991, str. 104-119. 21. Kuznecov, Y.N., Heat transfer in safety problems of nuclear reactors, Energoatomizdat, Moskva, 1989, str. 67. 22. Lahey, R.T., Moody, F.J., The thermal hydraulics of a boiling water nuclear reactor, American Nuclear Society, LaGrange Park, 1979. 23. Larsen, H.V., Bohm, B., Wigbels, M., A comparison of aggregated models for simulation and operational optimization of district heating networks, Energy Conversion & Management, Vol. 45, 2004, str. 1119-1139. 24. Liu, W.S., Tahir, A., Zaltsgendler, E., Kelly, W., Leung, R.K., Development Status of TUF Code, Proceedings of the 17 th Annual CNS Conference, Canadian Nuclear Society, Fredericton, Canada, 1996. 25. Liu, W.S., Hanna, B., Zaltsgendler, E., Advances in Modelling of Condensation Phenomena, Proceedings of the Annual Meeting of the Nuclear Society, Toronto, Canada, 1997. 165 26. Mamuzić, Z. P., Izabrana poglavlja iz oblasti običnih i diferencijalnih jednačina, Mašinski fakultet, Beograd, 1981. 27. Martin, C.S., Condensation – Induced Water Hammer in a Horizontal Pipe, Proceedings of the 3 rd IAHR International Meeting of the Workgroup on Cavitation and Dynamic Problems in Hydraulic Machinery and Systems, Brno, Czech Republic, 2009, str. 397-407. 28. Milivojevic, S., Stevanovic, V., Maslovaric, B., Numerical Simulation of Condensation Induced Water Hammer, The 15 th International Topical Meeting on Nuclear Reactor Thermal - Hydraulics (NURETH-15), Pisa, Italy, May 12-17, 2013, NURETH15-171. 29. Patankar, S.V., Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980. 30. Prica, S., Stevanović, V., Maslovarić, B., Numerical simulation of condensation induced waterhammer, Proceedings of the 12 th International Conference on Nuclear Engineering, April 25-29, 2004, Arlington, Virginia USA, ICONE12-49404. 31. Prica, S., Numerical Simulation of Condensation Induced Water Hammer, MSc. Thesis, Faculty of Mechanical Engineering, Belgrade, 2006. 32. Prica, S., Stevanović, V., Maslovarić, B., Numerical Simulation of Condensation Induced Water Hammer, FME TRANSACTIONS, New Series, Vol. 36, No. 1, 2008, str. 21-26. 33. Prica, S., Stevanović, V., Maslovarić, B., Vapour-Liquid Interface Tracking And Condensation Induced Water Hammer Predictions, 2nd International Congress of Serbian Society of Mechanics, Palić, Srbija, 2009. 34. Prica, S., Maslovarić, B., Stevanović, V., Numerical Prediction of Temperature Waves in Complex Pipeline Networks, III International Symposium Contemporary Problems of Fluid Mechanics, Beograd, Srbija, 2011. 35. Rapid Return on Investment on Temperature Optimization, 7 Technologies, Izveštaj dostupan na www.7t.dk, 2007. 36. Saito, T., Hughes, D. D., Carbon, M. W., Multi-fluid modeling of annular two-phase flow, Nuclear Engineering and Design, 50, 1978, str. 225-271. 37. Serkiz, A.W., An Evaluation of Water Hammer in Nuclear Power Plants, Proceedings of the 2 nd International Topical Meeting on Nuclear Reactor Thermalhydraulics, Santa Barbara, California, USA, 1983. 166 38. Stevanovic, V., Studovic, M., Bratic, A., Simulation and Analysis of a Main Steam Line Transient with Isolation Valves Closure and subsequent Pipe Break, International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 4, No. 5, 1994, str. 387-398. 39. Stevanovic, V., Jovanovic, Z., A hybrid method for the numerical prediction of enthalpy transport in fluid flow, International Communications in Heat and Mass Transfer, Vol. 27, No. 1, 2000, str. 23-34. 40. Stevanovic, V., Prica, S., Maslovaric, B., Zivkovic, B., Todorović, M., Galić, R., Optimizacija rada sistema daljinskog grejanja primenom numeričkih modela za simulaciju transporta toplote u složenim toplovodnim mrežama u stacionarnim i prelaznim režimima, Elaborat 1. godine rada na projektu Ministrastva nauke i zaštite životne sredine republike Srbije, broj 242008, Mašinski fakultet, Beograd, 2006. 41. Stevanović, V., Živković, B., Nikodijević, S., Maslovarić, B., Prica, S., Todorović, M., Galić, R., Hidraulički proračun složenih cevnih mreža sistema daljinskog grejanja, KGH – klimatizacija, grejanje, hlađenje, Godina 35, Broj 2, 2006, str. 27- 32. 42. Stevanović, V., Živković, B., Maslovarić, B., Prica, S., Todorović, M., Galić, R., Mandić, D., Dragojević, D., Nikodijević, S., Trkulјa, V., Merenje i simulacija prelaznih temperaturskih procesa u sistemu dalјinskog grejanja, KGH – klimatizacija, grejanje, hlađenje, Godina 36, Broj 2, 2007, str. 21-25. 43. Stevanović, V., Živković, B., Maslovarić, B., Prica, S., Todorović, M., Galić, R., Mandić, D., Dragojević, D., Nikodijević, S., Trkulјa, V., Termohidraulički proračuni sistema dalјinskog grejanja u cilјu povećanja energetske efikasnosti transporta toplote, Zbornik radova 38. međunarodni kongres o grejanju, hlađenju i klimatizaciji (KGH), Beograd, Srbija, 2007. 44. Stevanovic, V., Prica, S., Maslovaric, B., Zivkovic, B., Nikodijevic, S., Efficient numerical method for district heating system hydraulics, Energy Conversion and Management, Vol. 48, No. 5, 2007, str. 1536-1543. 45. Stevanovic, V., Jovanovic, M., Prica, S., Maslovaric, B., Condensation induced water hammer in thermal plants, Proceedings of the 11 th International Conference on Multiphase Flow in Industrial Plants, Palermo, Italy, September 7-10, 2008, str. 783-790. 46. Stevanovic, V., Živković, B., Prica, S., Maslovaric, Karamarkovic, V., Trkulja, V., Prediction of thermal transients in district heating systems, Energy Conversion and Management, Vol. 50, No. 9, 2009, str. 2167-2173. 167 47. Stevanovic, V., Prica, S., Maslovaric, Waterhammer in Pipelines of Steam Boilers, Proceedings of the 4 th IAHR International Meeting on Cavitation and Dynamic Problems in Hydraulic Machinery and Systems, Belgrade, Serbia, October 26-28, 2011, str. 57-65. 48. Stoop, P.M., Van der Bogaard, J.P.A., Koning, H., CHARME-01, А Thermo- Hydraulic Code for The Calculation of Fast Transients Inside Piping System, 8 th Int. Conf. on Structural Mechanics in Reactor Technology, Brussels, 1985. 49. Streeter, V.L., Wylie, E.B., Hydraulic Transients, McGraw Hill, New York, 1967. 50. Tannehill, J.C., Anderson, D.A., Pletcher, R.H., Computational Fluid Mechanics and Heat Transfer, Second Edition, Taylor & Francis, 1997. 51. Tong, L.S., Weisman, J., Thermal Analysis of Pressurized Water Reactors, ANS monograph, 1979. 52. Urban, C., Schlüter, M., Optical Investigations on Key Phenomena of Condensation Induced Water Hammers, The 15 th International Topical Meeting on Nuclear Reactor Thermal - Hydraulics (NURETH-15), Pisa, Italy, May 12-17, 2013, NURETH15-351 53. Vardy, A.E., Brown, J.M.B., Transient Turbulent Friction In Smooth Pipe Flows, Journal Of Sound And Vibration, 259(5), 2003, str. 1011–1036. 54. Versteeg, H.K., Malalasekera, W., An introduction to Computational Fluid Dynamics, Longman Group Ltd., Harlow, 1995. 55. Wagner, W., Kretzschmar, H.J., International Steam Tables, Springer-Verlag, Berlin, 2007. 56. Whalley, P.B., Two-Phase Flow and Heat Transfer, Oxford University Press, 1996. 57. Wulff, W., Computational Methods for multiphase flow, A Review prepared for the Second International Workshop on Two-Phase Flow Fundamentals, Rensselaer Polytechnic Institute Troy, New York, USA, 1987, str. 38-42. 58. Yeung, W.S., Wu, J., Fernandez, R.T., Sundaram, R.K., RELAP5/MOD3 Simulation of the Water Cannon Phenomenon, Nuclear Technology, Vol. 101, Feb. 1993, str. 244-251. 59. Zaltsgendler, E., Tahir, A., Leung, R.K., Condensation-Induced Waterhammer in a Vertical Upfill Pipe, Transactions of the American Nuclear Society, 74, 1996, str. 346. BIOGRAFIJA Sanja Milivojević, devojačko Prica, rođena je 10.05.1979. godine u Rijeci, Republika Hrvatska. U Rijeci je završila osnovnu školu i krenula u gimnaziju (prirodno- matematički smer). Srednješkolsko obrazovanje završila je 1997. godine u Zemunu, u “Zemunskoj gimnaziji”. Iste godine upisala je Mašinski fakultet Univerziteta u Beogradu. Diplomirala je 26.07.2002. godine na Katedri za termoenergetiku, odbranivši diplomski rad na temu “Numerička simulacija mehurastog toka u pravougaonom kanalu sa preprekom”, pod mentorstvom prof. dr Vladimira Stevanovića. Poslediplomske (magistarske) studije je upisala 2002. godine na Mašinskom fakultetu Univerziteta u Beogradu, na Katedri za Termoenergetiku, gde je i magistrirala 02.06.2006. godine odbranivši tezu pod naslovom "Numerička simulacija hidrauličkog udara izazvanog kondenzacijom pare". Mentor magistarskog rada bio je prof. dr Vladimir Stevanović. U periodu od 2003. do 2006. godine je bila stipendista-istraživač Ministarstva nauke i zaštite životne sredine Republike Srbije. Od aprila 2006. godine radi kao istraživač-saradnik na Katedri za termoenergetiku Mašinskog fakulteta u Beogradu, a od juna 2006. godine radi u Inovacionom centru Mašinskog fakulteta kao saradnik. Od 2007. godine je bila angažovana u nastavi na Mašinskom fakultetu u Beogradu iz predmeta Termodinamika B, Kompjuterske simulacije strujnotermičkih procesa i CFD, Generatori pare, Nuklearni reaktori i Dvofazna strujanja sa faznim prelazom. Učestvovala je u izradi 6 projekata finansiranih od strane Ministarstva prosvete, nauke i tehnološkog razvoja Republike Srbije i 11 projekata saradnje sa privredom. Autor je i koautor 34 naučno-stručna rada, od čega je 5 radova objavljeno u vodećim međunarodnim časopisima sa SCI liste, 5 radova je objavljeno u vodećim nacionalnim naučno-stručnim časopisima, 9 u zbornicima međunarodnih naučno-stručnih skupova, a 15 u zbornicima domaćih naučno-stručnih skupova. Koautor je jednog malog patenta zaštićenog kod Zavoda za intelektualnu svojinu Republike Srbije. Recenzent je studentskih naučnih radova na međunarodnoj konferenciji ICONE.