Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 UDK - UDC 536.423.1:621.039.534 Kratki znanstveni prispevek - Short scientific paper (1.03) Model stenskega uparjanja za popis podhlajenega vrenja toka pri nizkih tlakih Wall-Evaporation Model for a Description of Sub-Cooled Flow Boiling under Low- Pressure Conditions Boštjan Končar - Borut Mavko - Ivo Kljenak V prispevku je predstavljen model stenskega uparjanja za popis podhlajenega vrenja pri nizkih tlakih. Model temelji na lokalnih mehanizmih nukleacije mehurčkov in obsega vpliv debelitve toplotne mejne plasti ob greti steni. Model smo vključili v enorazsezni termo-hidravlični program RELAP5. Spremenjeni program smo testirali na številnih nizkotlačnih preizkusih podhlajenega vrenja iz literature ([1] do [4]). V nasprotju s sedanjim programom smo dosegli zelo dobro ujemanje izračunov z meritvami. © 2005 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: uparjanje stensko, vrenje toka podhlajeno, programi termohidravlični, RELAP5) In this paper we propose a wall-evaporation model for sub-cooled flow boiling under low pressures. The model is based on local mechanisms of bubble nucleation and captures the effect of the developing thermal boundary layer near the heated wall. The model was incorporated into the one-dimensional thermal-hydraulic code RELAP5. The modified code was validated against a number of published low-pressure sub-cooled boiling experiments ([1] to [4]), and in contrast to the existing code it shows good agreement with the experimental data. © 2005 Journal of Mechanical Engineering. All rights reserved. (Keywords: hot wall evaporation, sub-cooled flow boiling, thermal-hydraulic code, RELAP5) 0 UVOD Podhlajeno vrenje toka se pojavi, ko ob greti steni kanala nastaja para, kljub temu da je povprečna temperatura toka kapljevine po prerezu toka nižja od temperature vrelišča. Pojav je posebej pomemben v jedrskih reaktorjih hlajenih z vodo, pri katerih parni mehurčki v sredici vplivajo na obnašanje reaktorskega sistema, tako med normalnim obratovanjem kot tudi v primeru nezgodnih scenarijev. Za celovite analize prehodnih pojavov v hladilnih sistemih jedrskih reaktorjev se večinoma uporabljajo enorazsezni termohidravlični programi (RELAP5, TRAC, CATHARE, ATHLET). V zadnjih letih se je zaradi potreb po varnostnih analizah raziskovalnih jedrskih reaktorjev, ki obratujejo pri nizkih tlakih, in zaradi raziskav koncepta hlajenja zbiralnika zadrževalnega hrama prihodnjih 0 INTRODUCTION Sub-cooled flow boiling occurs when vapour is produced as a cold liquid flows along a heated channel, even though the average liquid temperature over the channel cross-section is lower than the liquid saturation temperature. The phenomenon is especially important in water-cooled nuclear power reactors, where the presence of vapour bubbles in the core influences the behaviour of the reactor system, either during normal operation or in the case of accident scenarios. To perform integral simulations of the transient phenomena in nuclear-reactor coolant systems one-dimensional thermal-hydraulic codes (such as RELAP5, TRAC, CATHARE and ATHLET) are widely used. In recent years, the interest in numerically simulating sub-cooled flow boiling at low pressures (1 to 3 bar) has increased, driven by the need to perform 646 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 lahkovodnih reaktorjev (PLVR - ALWR) povečalo zanimanje za numerične simulacije podhlajenega vrenja toka pri nizkih tlakih (1 do 3 bar). Vendar je večina konstitucijskih modelov v današnjih termohidravličnih programih razvita in preverjena pri visokih tlakih (zaradi uporabnosti v jedrskih elektrarnah), zato jih ne moremo brez sprememb uporabiti za preračune nizkotlačnih sistemov. Eden od glavnih ciljev razvijalcev programov je razširitev uporabnosti celovitih programov tudi v območje nizkih tlakov. V zadnjem času je bilo z različnimi termohidravličnimi programi izvedenih več simulacij podhlajenega vrenja pri nizkih tlakih. Rezultati analiz modelov podhlajenega vrenja ([5] do [8]), ki se uporabljajo v sedanjih termohidravličnih programih, so pokazali, da je izračunani delež parne faze pri nizkih tlakih mnogo nižji kakor tisti v preizkusih. Pri nizkih tlakih je razlika med gostotama kapljevine in pare bistveno večja kakor pri visokih tlakih (1:1590 pri 1 bar, 1:6.2 pri 150 bar), zato enaka količina nastale pare zavzame bistveno večjo prostornino. Pri visokih tlakih so parni mehurčki razmeroma majhni, medtem ko so pri nizkih tlakih le-ti bistveno večji. Večanje deleža parne faze vzdolž toka je tako močno odvisno od spremembe velikosti mehurčkov. V splošnem na količino proizvedene pare med pojavom podhlajenega vrenja vplivajo različni fizikalni mehanizmi, kakor so stensko uparjanje, kondenzacija parnih mehurčkov v podhlajeni kapljevini, koncentracija medfazne površine (določena z velikostjo mehurčkov) in medfazno trenje. V tem prispevku smo posebej obravnavali modeliranje stenskega uparjanja. Predlagali smo novi posplošeni model stenskega uparjanja za uporabo v termohidravličnih programih in ga vključili v program RELAP5. 1 MODEL STENSKEGA UPARJANJA Model stenskega uparjanja popisuje hitrost generacije parnih mehurčkov na greti površini stene. Predstavljen je enorazsežni model stenskega uparjanja za uporabo v termohidravličnih programih. Pri razvoju modela smo upoštevali naslednje predpostavke: 1. V razmerah nasičenja, v katerih sta kapljevina in para v toplotnem ravnovesju (Tl=T =T ), se celotni dovedeni toplotni tok porabi samo za generacijo pare. 2. Model uparjanja mora temeljiti na osnovnih safety analyses of research nuclear reactors operating at low pressures and to investigate the sump-cooling concept for advanced light-water reactors (ALWR). However, as most of the constitutive models in present-day thermal-hydraulic codes have been developed and validated for high-pressure conditions (due to the relevance for nuclear power plants), they cannot be straightforwardly applied for calculations involving low-pressure systems. One of the main goals of code developers is to extend the applicability of current integral codes to the low-pressure region. Recently, several simulations of sub-cooled boiling under low-pressure conditions using different thermal-hydraulic codes have been performed. The results of analyses of the sub-cooled boiling models ([5] to [8]) used in the current thermal-hydraulic codes have shown that, at low pressures, the calculated void fraction is significantly lower than in experiments. At low pressures, the difference between the liquid and vapour densities is much higher than at high pressures (1:1590 at 1 bar, 1:6.2 at 150 bar); therefore, the same amount of generated vapour occupies a significantly larger volume. The vapour bubbles are relatively small at high pressures, whereas they are much larger at low pressures. The increasing of the void fraction along the flow therefore strongly depends on the change of bubble size. In general, the amount of vapour produced during the sub-cooled flow boiling process is governed by different physical mechanisms, such as wall evaporation, the condensation of vapour bubbles in a sub-cooled liquid, the interfacial area concentration (determined by the bubble size) and the interfacial drag between the phases. In this paper we focus on the modelling of wall evaporation. A new generic wall-evaporation model for application in thermal-hydraulic codes is proposed and incorporated into the RELAP5 code. 1 THE WALL-EVAPORATION MODEL The wall-evaporation model describes the generation rate of vapour bubbles on a heated wall surface. A one-dimensional wall evaporation model for implementation in thermal-hydraulic codes is presented. When developing the model, the following assumptions were taken into account: 1. Under saturation conditions, where the liquid and vapour are in thermal equilibrium (Tl =Tg=Tsat ), the entire wall heat flux has to be used for vapour generation. 2. The model of wall evaporation has to be based on the basic parameters that describe bubble nucleation: Model stenskega uparjanja - Wall-Evaporation Model 647 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 parametrih, ki opisujejo mehurčkasto vrenje: velikost mehurčka, gostota nukleacijskih jeder N in frekvenca nukleacije mehurčkov f. 3. Pri nizkih tlakih mora biti delež toplotnega toka za uparjanje bistveno večji kakor v modelih stenskega uparjanja, ki se uporabljajo v sedanjih termohidravličnih programih. Kakor je navedeno v literaturi, termohidravlična programa RELAP5 ([7] in [8]) in ATHLET [5] napovesta prenizko hitrost stenskega uparjanja. Hitrost stenskega uparjanja lahko določimo z uporabo t.i. modela razdelitve toplotnega toka, ki predpostavlja, da se celotni dovedeni toplotni tok porabi delno za uparjanje pregrete kapljevine v mejni plasti tik ob greti površini in delno za segrevanje podhlajene kapljevine v jedru toka. Za osnovo pri razvoju modela stenskega uparjanja smo uporabili model razdelitve dovedenega toplotnega toka Kurula in Podowskega [9]. Po njunem modelu je vsaka enota grete površine razdeljena na dva dela: prvi del obsega vplivno območje nukleacije mehurčkov A b, medtem ko na preostalem delu grete površine A1, poteka konvektivni prenos toplote s stene na kapljevino. Abb in A1, v brezrazsežni obliki pomenita deleže celotne grete površine: bubble size, the number of active nucleation sites per unit surface and the nucleation frequency of the bubbles. 3. Under low-pressure conditions, the fraction of wall heat flux used for wall evaporation must be significantly higher than in the wall-evaporation models used in the existing thermal-hydraulic codes. According to the literature, the thermal-hydraulic codes RELAP5 ([7] and [8]) and ATHLET [5] predict a wall-evaporation rate that is too low. The wall-evaporation rate can be determined from the so-called heat-flux partitioning model, which assumes that the total applied heat flux is consumed partially for the evaporation of the superheated liquid in the boundary layer adjacent to the heated surface and partially for the heating of the sub-cooled liquid in the core flow. As a basis for the development of the wall-evaporation model, the heat-flux partitioning model of Kurul and Podowski [9] was used. According to their model, each unit of the heated surface consists of two parts: the first part of the heated area occupies the influence area of bubble nucleation, Abub, whereas single-phase convection from the wall to the liquid takes place on the remaining part of the heated area, A1f. In non-dimensional form, Abub and A1f represent the fractions of the total heated area: Abub + A1f =1 Celotni dovedeni toplotni tok iz stene na dvofazni vrelni tok sestoji iz treh različnih komponent (sl. 1): (1). The total heat flux transferred from the wall to the two-phase boiling flow consists of three different components (Figure 1): qw = q1f +qQ +qe (2), kjer: . q,, označuje toplotni tok zaradi enofazne konvekcije, ki se vrši zunaj vplivnega območja nukleacije mehurčkov; . qQ označuje toplotni tok zaradi površinskega hlajenja znotraj vplivnega območja mehurčkov Abub (to je prenos toplote s stene na podhlajeno kapljevino, ki periodično zapolnjuje izpraznjen prostor mehurčkov, ki se trgajo s stene med nukleacijskim krogom); . q označuje uparjalni toplotni tok, ki se porablja za neposredno generacijo parnih mehurčkov. Vplivna površina mehurčkov na enoto površine stene Ab b je definirana kot: where: . qu denotes the single-phase convection heat flux that takes place outside the influence area of the bubble nucleation, . qn denotes the surface-quenching heat flux within the bubble influence area, Abb (i.e., the heat transferred from the wall to the sub-cooled bulk liquid that periodically fills the volume vacated by the departing bubbles during the bubble nucleation cycle), . q denotes the wall-evaporation heat flux that is used for the direct generation of vapour bubbles. The bubble influence area per unit wall area, Abb, is determined as: 1,NK l^L (3), 648 Končar B. - Mavko B. - Kljenak I. Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 tok podhlajene kapljevine subcooled liquid flow Abub Sl. 1. Porazdelitev toplotnega toka na steni Fig. 1. Heat-flux partitioning at the wall kjer sta Na gostota aktivnih nukleacijskih jeder in db premer mehurčkov, nastalih na greti steni. Premer d bw je izračunan z uporabo Unalove korelacije [10], ki upošteva vpliv tlaka, podhladitve kapljevine, dovedenega toplotnega toka in hitrosti toka na velikost mehurčka. Parameter K določa velikost vplivnega območja mehurčka v okolici nukleacijskega jedra. Navadno se vzame nespremenljiva vrednost parametra K = 4 [11]. Pri največji gostoti aktivnih nukleacijskih jeder N vplivno območje mehurčkov pokriva celotno greto površino (Ab = 1). Pri vrednosti K= 4 vplivna površina mehurčkov postane enaka celotni greti površini, kadar povprečni razmik med dvema sosednjima nukleacijskima jedroma doseže velikost dveh premerov mehurčka (2 db). Osnovni model Kurula in Podowskega [9] smo prilagodili za uporabo v enorazsežnih termohidravličnih programih, ki uporabljajo spremenljivke povprečene po prečnem prerezu kanala. Uparjalni toplotni tok je sorazmeren gostoti aktivnih nukleacijskih jeder Na, frekvenci nukleacije mehurčkov f, masi mehurčka in uparjalni toploti hlg: where Na is the density of active nucleation sites and dbw is the diameter of the bubbles generated on the heated wall. The diameter dbw is calculated from the Unal correlation [10], which takes into account the effects of pressure, liquid sub-cooling, heat flux and liquid flow velocity on the bubble size. The parameter K determines the size of the bubble influence area around the nucleation site. A constant value of K = 4 is assumed [11]. At the maximum density of nucleation sites, Na, the bubble influence area covers the entire heating surface (Abub= 1). At the value of K= 4, the bubble influence area becomes equal to the total heating surface if the average spacing between two neighbouring nucleation sites reaches the size of two bubble diameters (2 dbw). The basic model of Kurul and Podowski [9] has been adapted for application in one-dimensional thermal-hydraulic codes, which use variables averaged over the channel cross-section. The evaporation heat flux is proportional to the density of nucleation sites, Na, the mass of a single bubble, the bubble nucleation frequency, f , and the latent heat, hlg: qe=Na f[^dbw \pghlg (4). Frekvenca nukleacije mehurčkov je definirana po Coleu (1960, citirano v Ivey [12]): The bubble nucleation frequency is defined by Cole (1960, cited by Ivey [12]): f = J- 4g(pl-pg) kjer je f odvisen samo od premera mehurčka dbw in od gostot faz rl in rg . Toplotni tok zaradi površinskega hlajenja izračunamo kot: 3 dbwPl (5), where f is affected only by the bubble diameter, dbw, and by the phase densities rl and rg. The quenching heat flux is calculated as: qQ=hQAbub{Tw-Tl\Š (6), Model stenskega uparjanja - Wall-Evaporation Model 649 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 kjer so hQ toplotna prestopnost, Tw temperatura stene in T povprečna temperatura kapljevine. Koeficient hQ izračunamo tako, da upoštevamo časovno odvisen prevod toplote s stene na polneskončno kapljevino [11]: where hQ is the heat-transfer coefficient, Tw is the wall temperature, and Tl is the average liquid temperature. The coefficient hQ is calculated by considering transient conduction from the wall to a semi-infinite liquid [11]: 1,6f Vp fkl lPc pl (7). Komponenta toplotnega toka enofazne konvekcije zunaj vplivnega območja mehurčkov je izračunana takole: The single-phase convection heat-flux component outside the bubble influence area is calculated in the following way: qu=hu-(1-Abub)-(Tw-Tl)-š (8), kjer je toplotna prestopnost h1d za enofazni turbulentni konvektivni tok, izračunana po Kurulu in Podowskem [9]: where the heat-transfer coefficient h1f for the singlephase turbulent convective flow is calculated according to Kurul and Podowski [9]: h1=St-Prcp (9), kjer je St Stantonovo število. Faktor ^smo vključili v en. (6) in (8) za segrevanje kapljevine, da bi zadostili pogoju, da se pri nasičenju celotni dovedeni toplotni tok porabi za generacijo pare: where St is the Stanton number. The inhibiting factor x was introduced in Eqs. (6) and (8) for the liquid-heating components to satisfy the condition that at saturation, the entire wall heat flux is consumed for vapour generation: x ( hlsat l ) j ( hlsat hlcr ) 1 + (hat-hl)l(hat-hlcr) (10). Faktor L je padajoča funkcija podhladitve kapljevine in popisuje dvorazsežno naravo razvijajoče se toplotne mejne plasti ob steni, kjer so temperaturni gradienti največji. Z debelitvijo plasti vzdolž kanala se vedno več kapljevine uparja, medtem ko je vedno manj toplote na voljo za segrevanje kapljevine. Kritična entalpija hlc , ki je izračunana po korelaciji Saha-Zuber [13], označuje mesto, kjer postane vpliv toplotne mejne plasti pomemben. Temperaturo stene T in temperaturo kapljevine Tl v en. (6) in (8) izračunamo z implicitno sklopitvijo energijske enačbe in enačbe za prevod toplote v steni [14], tako da je gostota nukleacijskih jeder N edina preostala neznana veličina. Iz vsote komponent dovedenega toplotnega toka lahko izračunamo neznano vrednost N: The factor x is a decreasing function of the liquid sub-cooling and captures the two-dimensional nature of the developing thermal boundary layer near the heated wall, where the liquid temperature gradient is significant. As the layer gets thicker along the channel, more vapour is being generated, thus less heat is transferred to the liquid. The critical enthalpy, hlcr, calculated from the Saha-Zuber correlation [13], defines the location where the effect of the thermal boundary layer becomes significant. The wall temperature, Tw, and average liquid temperature, Tl, in Eqs. (6) and (8) are calculated by the implicit coupling of the energy equation and the equation for the heat-transfer conduction through the wall [14], so that Na remains the only unknown variable. The unknown value of Na can be calculated from the sum of the heat-flux components: N h1 (Tw-TlK n-d f 6 dbw-hfg-Pg(1 + L) K h1 (Tw-Tl)š (11), kjer je e faktor površinskega hlajenja, ki je definiran kot razmerje med toplotnima tokovoma qQ in q: where e is the quenching factor, defined as the ratio between heat fluxes qQ and qe: qQ 3 \hQ-T(T-Tl) 2 J dbw-p^hlg (12). 650 Končar B. - Mavko B. - Kljenak I. Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 Ko poznamo vrednosti komponent dovedenega toplotnega toka, lahko iz uparjalnega toplotnega toka q izračunamo hitrost stenskega uparjanja Hitrost stenskega uparjanja Gw je definirana kot količina pare na enoto časa in prostornine, ki nastaja v pregreti toplotni mejni plasti kapljevine ob steni: kjer sta A in V greta površina in prostornina računske celice. 2 VKLJUČITEV MODELA V PROGRAM RELAP5 Termohidravlični program RELAP5 se zelo pogosto uporablja za izvajanje varnostnih analiz v jedrskih reaktorjih. Sedanji program RELAP5 uporablja Laheyev model [15] za izračun stenskega uparjanja. Naše predhodne analize [7] so pokazale, da Laheyev model izračuna prenizke hitrosti uparjanja pri preizkusnih podhlajenega vrenja, izvedenih pri tlakih pod 3 bar. Model ne upošteva eksplicitno premera mehurčka, ki pri nizkih tlakih bistveno vpliva na hitrost stenskega uparjanja in posledično na količino pare v vrelnem kanalu. Zato smo v zadnjo razpoložljivo verzijo programa RELAP5/MOD3.2.2 Gamma vključili novi model stenskega uparjanja. Program temelji na dvofluidnem sistemu neravnotežnih enorazsežnih transportnih enačb [14] za prenos: mase: kjer se spodnji indeks k nanasa bodisi na kapljevito fazo (k=l) bodisi na plinasto fazo (k=g). Členi na levi strani enačb (14) do (16) opisujejo spremembo fizikalnih veličin po času in kraju. Člen izmenjave mase na desni strani en. (14) lahko razdelimo na prenos mase na medfaznem stiku pare in kapljevine Gik (kondenzacija) in na prenos mase zaradi Knowing the values of the heat-flux components, the wall-evaporation rate can be calculated from the evaporation heat flux, qe. The wall-evaporation rate, Gw, is defined as the amount of steam per unit time and unit volume generated in the superheated liquid thermal boundary layer near the wall: where Aw and V are the heated area and the volume of the computational cell. 2 INCLUSION OF THE MODEL IN THE RELAP5 CODE The thermal-hydraulic code RELAP5 is frequently used to perform safety analyses of nuclear reactors. The existing code RELAP5 uses Lahey’s model [15] to calculate wall evaporation. Our previous analyses [7] have shown that Lahey’s model calculates evaporation rates that are too low, for sub-cooled boiling experiments performed at pressures below 3 bar. The model does not explicitly take into account the bubble size, which significantly affects the wall-evaporation rate at low pressures, and consequently the amount of vapour in the boiling channel. Thus, the new wall-evaporation model was included in the latest available code version RELAP5/ MOD3.2.2. Gamma. The code is based on the two-fluid system of non-equilibrium one-dimensional transport equations [14] for mass where the subscript k refers either to the liquid phase (k=l) or to the gas phase (k=g). The terms on the left-hand side of Eqs. (14–16) describe the variation of physical variables in time and space. The massexchange term on the right-hand side of Eq. (14) can be divided into the mass-transfer term at the vapour-liquid interface, Gik (e.g., condensation) and the mass- (13), Jt (akrk)+~Adx (akrkvkA)=Gk (14), gibalne količine: momentum ak rk—vk+-akrk—v 2 =-ak—p + akrkg +WFRk+Gk(vik-vk) -Ci |vr|vr -CVMaVM (15), in notranje energije: and internal energy jt (akrk ek) + -—(akrkekvkA) = -p — ak - A—(akvkA) + Qwk+Qik +Gkhk +DISSk (16), Model stenskega uparjanja - Wall-Evaporation Model 651 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 uparjanja v toplotni mejni plasti ob steni Gw (stensko uparjanje): Členi na desni strani gibalne en. (15) po vrsti pomenijo: silo tlaka, silo teže, trenje ob steni, prenos gibalne količine zaradi izmenjave mase med fazama, medfazno trenje in silo navidezne mase. Prva dva člena na desni strani energijske en. (16) pomenita delo tlaka, tretji člen toplotni tok, ki prehaja s stene na fazo k, četrti člen medfazni prenos toplote, peti člen prenos notranje energije zaradi izmenjave mase med fazama in zadnji člen disipacijo notranje energije zaradi stenskega trenja. Prenos mase, energije in gibalne količine med fazama in med dvofaznim tokom in steno je popisan z dodatnimi zapiralnimi enačbami [14]. V primeru podhlajenega vrenja v toku, se medfazna kondenzacija (Gik < 0) in stensko uparjanje (Gw > 0) pojavljata sočasno. Hitrost stenskega uparjanja Gw je definirana z en. (13). 3 REZULTATI IN RAZPRAVA Pri podhlajenem vrenju toka je neto količina pare v kanalu odvisna od interakcije med stenskim uparjanjem, kondenzacijo in relativno hitrostjo mehurčkov glede na kapljevino. Ker je glavni namen prispevka modeliranje stenskega uparjanja, podrobnosti modeliranja preostalih mehanizmov niso obravnavane (najdemo jih lahko v [16]). Nabor zapiralnih enačb za stensko uparjanje, kondenzacijo in relativno hitrost med fazama (model gonilnega toka) smo z uporabo lastnih FORTRAN-skih podprogramov [17] vključili v sedanji program RELAP5/MOD3.2.2 Gamma. Izračune s spremenjenim programom (označenim z RELAP5_new) smo primerjali z izračuni s sedanjim programom (označenim z RELAP5_old) in z nizkotlačnimi preizkusi podhlajenega vrenja iz literature ([1] do [4]). Obratovalne razmere obravnavanih preizkusov so podane v preglednici 1. Preizkusi Zeitouna in Shoukrija ([1] in [2]) ter Donevskega in Shoukrija [3] so bili izvedeni v 306 mm dolgi greti cevi kolobarjastega prečnega prereza z notranjim premerom 13 mm in zunanjim premerom 25 mm. Testna sekcija Dimmicka in Selanderja [4] je bila 308 mm dolga cev z greto steno notranjega premera 12,3 mm. Za vsak preizkus smo razvili vhodni transfer due to evaporation in the wall’s thermal boundary layer, Gw (wall evaporation): ^ (17). The terms on the right-hand side of the momentum Eq. (15) represent, respectively, the pressure gradient, the body force, the wall friction, the momentum transfer due to mass exchange, the interfacial drag force and the virtual mass force. The first two terms on the right-hand side of Eq. (16) represent the pressure work, the third term represents the heat flux from the wall to phase k, the fourth term describes the interfacial heat transfer, the fifth term describes the internal energy transport due to mass exchange and the last term represents the energy dissipation due to wall friction. The transfer of mass, energy and momentum between the phases and between the two-phase flow and the wall is described by additional closure equations [14]. In the case of sub-cooled boiling flow, interfacial condensation (Gik < 0) and wall evaporation (Gw > 0) occur simultaneously. The wall-evaporation rate, Gw, is defined by Eq (13). 3 RESULTS AND DISCUSSION The net amount of vapour in the channel during sub-cooled flow boiling depends on the interaction between the wall evaporation, the condensation and the relative velocity of the bubbles with regard to the liquid. As this paper focuses on wall-evaporation modelling, the details about the modelling of other mechanisms are not provided; however, they can be found in [16]. The set of closure relations for the wall evaporation, the condensation and the relative velocity between phases (drift-flux model) was incorporated into the current code RELAP5/MOD3.2.2 Gamma via our own FORTRAN subroutines [17]. The calculations with the modified code (denoted as RELAP5_new) were compared with the existing code calculations (denoted as RELAP5_old) and against low-pressure sub-cooled boiling experiments from the literature ([1] to [4]). The operating conditions of the experiments are given in Table 1. The experiments of Zeitoun and Shoukri ([1] and [2]) and Donevski and Shoukri [3] were performed in a 306-mm-long heated annulus with inner and outer diameters of 13 mm and 25 mm, respectively. The test section of Dimmick and Selander [4] was a 308-mm-long cylindrical tube with a heated wall, having an inner diameter of 12.3 mm. An input model was developed for each 652 Končar B. - Mavko B. - Kljenak I. Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 Preglednica 1. Obratovalne razmere pri preizkusih Table 1. Operating conditions of the experiments Vir Reference [1] 1 213,6 161,2 1,14 Št. preiz. Exp. No. 2 480 [2] 3 508 [4] 7 805 qw kW/m2 4 596 [3] 5 586,4 315,1 2,11 8 472 G kg/m2s 208,1 264,1 263,8 6 481,4 392,1 1,54 620,2 630,7 pin bar 1,14 1,5 1,2 1,65 1,65 DTsub,in oC 13,1 20 16,8 20,1 23,7 18,5 44,3 27,5 model, kjer smo definirali geometrijsko obliko testne sekcije, robne pogoje in diskretizacijo računskega območja [17]. Na slikah 2 in 3 so predstavljeni značilni parametri novega modela stenskega uparjanja. S spremenjenim programom smo simulirali preizkus št.1 (preglednica 1). Na sliki 2 je predstavljena razdelitev dovedenega toplotnega toka vzdolž kanala. Vidimo lahko, da se enofazni konvekcijski toplotni tok qM zmanjšuje, predvsem na račun večanja uparjalnega toplotnega toka q . Največji delež dovedenega toplotnega toka se prenaša na podhlajeno kapljevino zaradi površinskega hlajenja q Območje »polno razvitega vrenja«, kjer je q = 0 Q. dosežemo približno na polovici gretega dela kanala. V tem območju vplivno območje mehurčkov pokriva celotno greto površino (Ab b =1). Kakor je prikazano na sliki 3, se premer mehurčka dbw zvečuje vzdolž kanala, kar je v 1 experiment, where the test-section geometry, the boundary conditions and the discretization of the computational domain were defined [17]. Figs. 2 and 3 present the characteristic parameters of the new wall-evaporation model. Experiment No.1 (Table 1) was simulated with the modified code. Fig. 2 presents the partitioning of the wall heat flux along the heated channel. As can be seen, the single-phase convection heat flux, q1f, decreases, mostly on account of the increasing of the evaporation heat-flux component, qe. Most of the wall heat flux is transferred to the sub-cooled liquid as a result of quenching, qQ. A “fully developed boiling” region, with q1f = 0, is reached approximately halfway along the heated channel. In this region, the bubble influence area covers the entire heated surface (Abub =1). As shown in Fig. 3, the bubble diameter, dbw, increases along the channel, which is in accordance 0.8 0.6 0.4 0.2 0 0 H- ¦A-q1f/qw ^¦qq/qw ^qe/qw -H«, -H*. "H polno razvito vrenje fully developed boiling 0.1 z (m) 0.2 0.3 Sl. 2. Izračunana porazdelitev dovedenega toplotnega toka vzdolž kanala (primer 1, preglednica 1) Fig. 2. Calculated partitioning of the wall heat flux along the channel (Case 1, Table 1) Model stenskega uparjanja - Wall-Evaporation Model 653 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 0.003 0.0025 0.002 0.0015 0.001 0.0005 polno razvito vrenje fully developed boiling 120000 100000 80000 60000 40000 20000 0.1 z (m) 0.2 0.3 Sl. 3. Izračunani razvoj premera mehurčka dbw in gostote aktivnih nukleacijskih jeder N vzdolž kanala (primer 1, preglednica 1) Fig. 3. Calculated evolution of the bubble diameter, dbw, and the density of the active nucleation sites, Na, along the channel (Case 1, Table 1) skladu s preizkusnimi opažanji [2]. Gostota nukleacijskih jeder N se počasi zvečuje v območju »delno razvitega vrenja « (Ab <1), medtem ko se začne zmanjševati v območju polno razvitega vrenja zaradi nadaljnjega večanja premera mehurčka dbw . Na sliki 4 je za preizkusni primer 1 prikazana primerjava hitrosti stenskega uparjanja, izračunanih s spremenjenim (RELAP5_new) in sedanjim programom (RELAP5_old). Kakor vidimo, so hitrosti stenskega uparjanja bistveno višje v primeru spremenjenega programa. 4 3 -----RELAP5_new - - -RELAP5_old with experimental observations [2]. The nucleation site density, Na, slowly increases in the “partially developed boiling” region (Abub <1), whereas it starts decreasing in the fully developed boiling region due to a further increase in the bubble diameter, dbw. In Fig. 4, the wall-evaporation rates calculated by the modified (RELAP5_new) and the existing codes (RELAP5_old) are compared for the same experimental case 1. As can be seen, the wall-evaporation rates are significantly higher in the case of the modified code calculation. 2 0 0.1 z (m) 0.2 0.3 Sl. 4. Izračunane hitrosti stenskega uparjanja vzdolž kanala (primer 1, preglednica 1) Fig. 4. Calculated wall-evaporation rates along the channel (case 1, Table 1) 0 0 0 1 0 654 Končar B. - Mavko B. - Kljenak I. Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 0.0035 0.003 0.0025 0.002 0.0015 0.001 0.0005 — RELAP5_new i Exp. 3 0 0 0 0.3 0.0035 0.003 0.0025 0.002 0.0015 0.001 0.0005 0 RELAP5_new Exp. 4 0.1 0.2 z (m) 0.3 0.1 0.2 z (m) Sl. 5. Izračunane in izmerjene [2] vrednosti premera mehurčka vzdolž kanala (primera 3 in 4, preglednica 1) Fig. 5. Calculated and measured [2] values of the bubble diameter along the channel (cases 3 and 4, Table 1) Preglednica 2. Izmerjene ([1] in [2]) in izračunane (RELAP5_new) vrednosti povprečne temperature stene Table 2. Measured ([1] and [2]) and calculated (RELAP5_new) values of the average wall temperature Št. preiz. Exp. No. Tw,exp K Tw,calc (RELAP5_new) K (Twcalc - Twexp)x100/Twexp % 1 393,0 391,8 -0,3 2 400,0 400,1 0,025 3 410,5 410,0 -0,12 4 406,3 407,0 0,17 Ker preizkusnih podatkov o gostoti aktivnih nukleacijskih jeder ni na voljo, smo izračun N preverili posredno, s primerjavo izračunanih in izmerjenih vrednosti premera mehurčka in temperature stene. Po en. (11) je gostota aktivnih nukleacijskih jeder odvisna od premera mehurčka in temperature stene. Kakor je prikazano na sliki 5, izračun s spremenjenim programom RELAP5 razmeroma dobro napove izmerjene [2] vrednosti premera mehurčka (povprečene po prerezu kanala). V preglednici 2 so za 4 preizkusne primere (preglednica 1) podane izmerjene in izračunane vrednosti temperature grete stene. Za vsak preizkusni primer je podana ena vrednost, ki pomeni temperaturo povprečeno po celotni notranji površini grete stene. Povprečenje je smiselno, saj so tako izmerjene kakor izračunane temperaturne porazdelitve vzdolž kanala skoraj nespremenjene. Kakor je prikazano v preglednici 2, je tudi ujemanje izmerjenih in izračunanih temperatur sten zelo dobro. Celotni model podhlajenega vrenja sestoji iz posameznih podmodelov, ki vplivajo na izračun deleža pare v gretem kanalu. Poleg modela stenskega uparjanja sta najpomembnejša model kondenzacije in model gonilnega toka. Da bi zaupali izračunu hitrosti stenskega uparjanja, smo preverili As the experimental data on the nucleation site density, Na, are not available, the calculation of Na was validated indirectly, by comparing the calculated bubble diameter and the wall temperature with the experimental values. According to Eq. (11), the active nucleation site density depends on the bubble diameter size and on the wall temperature. As shown in Fig. 5, the measured [2] bubble diameter values (averaged over the channel cross-section) are reasonably well predicted by the modified RELAP5 calculation. Table 2 provides measured and calculated values of the heated wall temperature for four experimental cases (Table 1). For each experimental case, one value is provided, representing the temperature averaged over the entire inner surface of the heated wall. The averaging is reasonable since both the measured and calculated temperature distributions along the channel are almost constant. As presented in Table 2, very good agreement between the measured and calculated wall temperatures is also obtained. The overall sub-cooled boiling model consists of several sub-models that influence the calculation of the vapour content in the heated channel, [16]. Besides the wall-evaporation model, the condensation model and the drift-flux model are the most important. To have confidence in the wall-evaporation rate calculation, Model stenskega uparjanja - Wall-Evaporation Model 655 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 tudi ta dva modela. Na sliki 6 je prikazana primerjava izračunanih in izmerjenih [1] porazdelitev deleža parne faze v gretem in v neogrevanem delu navpičnega kanala. V gretem delu kanala na delež parne faze vplivata stensko uparjanje in kondenzacija, medtem ko v neogrevanem delu kanala poteka le kondenzacija parnih mehurčkov, zato se delež parne faze hitro zmanjšuje vzdolž kanala. Izračuni s spremenjenim programom (RELAP5_new) se dobro ujemajo z izmerjenim manjšanjem deleža parne faze in potrjujejo, da je bil uporabljen ustrezen model kondenzacije. V primeru izračunov s sedanjim programom (RELAP5_old) so vrednosti deleža parne faze bistveno nižje vzdolž celotne dolžine kanala, predvsem zaradi premajhne generacije pare ob steni. V neogrevanem delu kanala je nagib manjšanja deleža parne faze manj strm kakor v preizkusu, kar kaže na prenizko hitrost kondenzacije. Spremenjeni program RELAP5 (RELAP5_new) vključuje tudi novi model gonilnega toka [16], ki določa relativno hitrost med fazama. Primerjava izračunanih relativnih hitrosti z izmerjenimi vrednostmi je pokazala dobro ujemanje v primeru izračunov s spremenjenim programom, medtem ko je sedanji program (RELAP5_old) napovedal prevelike relativne hitrosti [16]. Na sliki 7 je prikazana primerjava izračunanih (RELAP5_old in RELAP5_new) in izmerjenih deležev parne faze vzdolž gretega kanala. Večja hitrost stenskega uparjanja v primeru izračuna RELAP5_new vpliva na večji delež pare v kanalu. Izračuni s sedanjim programom napovejo precej manjši delež parne faze od izmerjenega, medtem ko se izračuni s spremenjenim programom dobro ujemajo z izmerjenimi vrednostmi. 0.5 0.4 0.3 0.2 0.1 0 0 • Exp. 1 [1] ¦ ¦ RELAP5_old — RELAP5_new neogrevani del unheated part these two models were also validated. In Fig. 6. a comparison of the calculated and measured [1] void-fraction distributions in the heated and unheated parts of the vertical channel is presented. In the heated part of the channel both the wall evaporation and the condensation have an influence on the void fraction, whereas in the unheated part of the channel only the condensation of vapour bubbles is present, causing the void fraction to decrease rapidly along the channel. The modified code calculations (RELAP5_new) agree well with the measured void-fraction decrease and confirm that an appropriate condensation model was used. The void-fraction values along the entire length of the channel are significantly lower in the case of the existing code calculations (RELAP5_old), mostly due to the vapour generation near the wall being too low. The slope of the void-fraction decrease in the unheated part of the channel is less steep, as in the experiment, which indicates a condensation rate that is too low. The modified RELAP5 code (RELAP5_new) also includes a new drift-flux model [16], which determines the relative velocity between the phases. A comparison of the calculated and measured values shows good agreement in the case of modified code calculations, whereas the existing code (RELAP5_old) predicted relative velocities that were too high [16]. A comparison of the calculated (RELAP5_old and RELAP5_new) and the measured void fractions in the heated channels is presented in Fig. 7. A higher wall-evaporation rate in the case of the RELAP5_new calculations leads to higher vapour contents in the channel. The existing code calculations significantly under-predict the measured void fraction, whereas the calculations with the modified code show good agreement with the measured values. 0.1 0.2 0.3 0.4 0.5 0.4 0.3 0.2 0.1 0 • Exp. 2 [1] ¦ ¦ RELAP5_old RELAP5_new neogrevani del unheated part 0.5 0 0.1 0.2 0.3 0.4 0.5 z (m) z (m) Sl. 6. Izračunane in izmerjene [1] porazdelitve deleža parne faze v gretem in neogrevanem delu kanala (primera 1 in 2, preglednica 1) Fig. 6. Calculated and measured [1] void-fraction distributions in the heated and unheated parts of the channel (Cases 1 and 2, Table 1) 656 Končar B. - Mavko B. - Kljenak I. Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 Exp. 1 [1] RELAP5_old RELAP5_new 0 0.1 0.2 z (m) 0.5 0.4 0.3 0.2 0.1 0 0 0.05 0.1 0.15 0.2 0.25 z (m) 0.3 0.2 0.1 • Exp. 5 [3] ¦ - RELAP5_old -----RELAP5_new 0.1 0.2 z (m) 0.3 0.2 0.1 Exp. 7 [4] RELAP5_old RELAP5_new 0.3 0.3 0.5 0.4 0.3 0.2 0.1 0 0.5 0.4 0.3 0.2 0.1 0 0 0.3 • Exp. 2 [1] RELAP5_old -----RELAP5_new 0.05 0.1 0.15 0.2 0.25 0.3 z (m) Exp. 4 [2] RELAP5_old RELAP5_new 0.05 0.1 0.15 0.2 0.25 0.3 z (m) 0.2 0.1 • Exp. 6 [3] ¦ - RELAP5_old -----RELAP5_new 0.3 0.1 0.2 z (m) 0.3 ^^* • • 0.2 0.1 Exp. 8 [4] RELAP5_old RELAP5_new -0.2 -0.15 -0.05 -0.2 -0.15 -0.1 z (m) Sl. 7. Primerjava vzdolžnih porazdelitev deležev parne faze Fig. 7. Comparison of the axial void-fraction distributions -0.1 z (m) -0.05 4 SKLEPI 4 CONCLUSIONS 0.3 Z namenom, da bi razširili uporabnost To extend the range of applicability of the termohidravličnih programov v območje nizko- thermal-hydraulic codes to low-pressure sub-cooled tlačnega podhlajenega vrenja, smo razvili model boiling, a wall-evaporation model consistent with 0 0 0 0 0 0 0 Model stenskega uparjanja - Wall-Evaporation Model 657 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 stenskega uparjanja, ki je skladen z osnovnimi mehanizmi nukleacije mehurčkov na greti steni. Model upošteva delitev dovedenega toplotnega toka na delež, ki se porabi za generacijo pare in delež, ki se porabi za segrevanje kapljevine. Model je tudi zmožen napovedati prehod iz delno razvitega v polno razvito območje podhlajenega vrenja. Novi model stenskega uparjanja smo vključili v program RELAP5. Spremenjeni program RELAP5 poleg novega modela stenskega uparjanja vključuje tudi nova modela kondenzacije in gonilnega toka, ki vsi skupaj vplivajo na delež parne faze pri podhlajenem vrenju toka. Spremenjeni program smo preverili na vrsti nizkotlačnih preizkusov podhlajenega vrenja. Dosegli smo zelo dobro ujemanje izračunanih in izmerjenih deležev parne faze vzdolž kanala. Nasprotno je sedanji program napovedal bistveno premajhen delež parne faze pri vseh obravnavanih preizkusih. the basic physical mechanisms of bubble nucleation on a heated wall has been developed. The model considers the decomposition of the wall heat flux into vapour-generation and liquid-heating parts, and is able to predict the transition from the partially developed to the fully developed sub-cooled boiling region. The new wall-evaporation model was incorporated into the RELAP5 code. Besides the new wall-evaporation model, the modified RELAP5 code also includes new models for condensation and drift-flux, which altogether influence on the void fraction during sub-cooled flow boiling. The modified code was verified against a series of low-pressure sub-cooled boiling experiments. Very good agreement between the calculated and measured void fractions along the channel was achieved. In contrast, the existing code significantly under-predicts the void-fraction data in all the experiments. 5 SPREMENLJIVKE 5 NOMENCLATURE pospešek navidezne mase prerez greta površina specifična toplota kapljevine koeficient medfaznega trenja koeficient navidezne mase premer mehurčka, ki nastane na steni hidravlični premer kanala specifična notranja energija težnostni pospešek masni pretok specifična entalpija uparjalna toplota toplotna prevodnost kapljevine Nusseltovo število = q D/kl tlak w e Pecletovo število = GD c /kl Prandtlovo število kapljevine površinski toplotni tok vir toplote na enoto prostornine Reynoldsovo število toka = GD /m Stantonovo število = Nu/RePrl e temperatura hitrost nadzorna prostornina relativna hitrost med fazama Grške črke prostorninski delež faze k VM a A Aw c pl Ci C VM d bw De e g G h h lg kl Nu p Pe Prl q Q Re St T v V vr ak m/s2 m2 m2 J/kgK m m J/kg m/s2 acceleration due to virtual mass cross-section heated area liquid specific heat interfacial drag coefficient virtual mass coefficient bubble diameter generated on the wall channel equivalent diameter specific internal energy acceleration due to gravity kg/m2s mass flux J/kg specific enthalpy J/kg latent heat W/mK liquid thermal conductivity Nusselt number = qw De/kl Pa, bar pressure Peclet number = GDecp,l/kl liquid Prandtl number surface heat flux heat source per unit volume Reynolds flow number = GDe/m Stanton number = Nu/RePrl temperature velocity control volume relative velocity between phases Greek symbols volume fraction of phase k W/m2 W/m3 K m/s m3 m/s 658 Končar B. - Mavko B. - Kljenak I. Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 hitrost prenosa mase G kg/sm3 mass transfer rate dinamična viskoznost kapljevine m Pa s liquid viscosity gostota r kg/m3 density Spodnji indeksi Subscripts od meje faz proti fazi k ik from interface to phase k vstopni pogoji in inlet conditions kapljevina l liquid para g vapour faza k k phase k nasičenje sat saturation podhlajeno sub sub-cooling stena w wall od stene proti fazi k wk from wall to phase k 6 LITERATURA 6 REFERENCES [I] Zeitoun, O., Shoukri, M. (1997) Axial void fraction profile in low pressure subcooled flow boiling. Int. J. Heat Mass Transfer 40 (1997), pp. 869-879. [2] Zeitoun, O., Shoukri, M. (1996) Bubble behavior and mean diameter in subcooled flow boiling. J. Heat Transfer-Trans. ASME 118 (1996), pp. 110-116. [3] Donevski, B., Shoukri, M. (1989) Experimental study of subcooled flow boiling and condensation in an annular channel. Thermofluids report No. ME/89/TF/R, Department of Mechanical Engineering, McMaster University, Hamilton, ON. [4] Dimmick, G.R., Selander, W.N. (1990) A dynamic model for predicting subcooled void: experimental results and model development. EUROTHERM seminar, Pisa, Italy, 1990. [5] Hainoun, A., Hicken, E., Wolters, J. (1996). Modelling of void formation in the subcooled boiling regime in the ATHLET code to simulate flow instability for research reactors. Nucl. Eng. Des. 167 (1996), pp. 175- 191. [6] Devkin, A.S., Posedonov, A.S. (1998) RELAP5/MOD3 subcooled boiling model assessment, NUREG/IA- 0025. U.S. NRC, Washington, DC. [7] Končar, B., Mavko, B., (2000) RELAP5 modeling of subcooled boiling in vertical flow. Proceedings of the EigthInternational Conference ICONE-8, Baltimore, MD, USA, 2000. [8] Hari, S., Hassan, Y.A., (2002) Improvement of the subcooled boiling model for low-pressure conditions in thermal-hydraulic codes. Nucl. Eng. Des. 216 (2002), pp. 139-152. [9] Kurul, N, Podowski, M.Z. (1990) Multidimensional effects in forced convection subcooled boiling. Proceedings of the Ninth International Heat Transfer Conference, Vol. 2, Jerusalem, Israel, 1990. [10] Unal, H.C. (1976) Maximum bubble diameter, maximum bubble-growth time and bubble-growth rate. Int. J. Heat Mass Transfer 19 (1976), pp. 643-649. [II] Victor, H, Del Valle, M., Kenning, D.B.R. (1985) Subcooled flow boiling at high heat flux. Int. J. Heat Mass Transfer 28 (1985), pp. 1907-1920. [12] Ivey, H.J. (1967) Relationships between bubble frequency, departure diameter and rise velocity in nucleate boiling. Int. J. Heat Mass Transfer 10 (1967), pp. 1023-1040. [13] Saha ,P, Zuber, N. (1974). Point of net vapor generation and void fraction in subcooled boiling. Proceedings of the Fifth International Heat Transfer Conference, Tokyo, Japan, 1974. [14] RELAP5 Code Development Team (1999) RELAP5/MOD3 Code Manual, Volume IV: Models and correlations, NUREG/CR-5535. Scientech, Inc., Idaho Falls, ID. [15] Lahey, R.T., (1978) A mechanistic subcooled boiling model. Proceedings of the Sixth International Heat Transfer Conference, Vol. I, Toronto, Canada, 1978. Model stenskega uparjanja - Wall-Evaporation Model 659 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)10, 646-660 [16] Končar, B., Mavko, B., (2003) Modelling of low-pressure subcooled flow boiling using the RELAP5 code. Nucl. Eng. Des. 220 (2003), pp. 255-273. [17] Končar, B. (2003) Model of forced convective subcooled boiling at low pressure conditions, PhD Thesis, Faculty of Mathematics and Physics, University of Ljubljana, Slovenia. Naslov avtorjev: dr. Boštjan Končar Authors’ Address: Dr. Boštjan Končar prof.dr. Borut Mavko Prof.Dr. Borut Mavko dr. Ivo Kljenak Dr. Ivo Kljenak Institut Jožef Stefan Jožef Stefan Institute Odsek za reaktorsko tehniko Reactor Engineering Division Jamova 39 Jamova 39 1000 Ljubljana 1000 Ljubljana, Slovenia bostjan.koncar@ijs.si bostjan.koncar@ijs.si borut.mavko@ijs.si borut.mavko@ijs.si ik@ijs.si ik@ijs.si Prejeto: Sprejeto: Odprto za diskusijo: 1 leto 15.7.2003 25.5.2005 Received: Accepted: Open for discussion: 1 year 660 Končar B. - Mavko B. - Kljenak I.