© Strojni{ki vestnik 46(2000)9,607-621 © Journal of Mechanical Engineering 46(2000)9,607-621 ISSN 0039-2480 ISSN 0039-2480 UDK 532.5:519.61/.64:621.039 UDC 532.5:519.61/.64:621.039 Izvirni znanstveni ~lanek (1.01) Original scientific paper (1.01) Vpliv natan~nosti numeri~nih metod na rezultate simuliranj me{alne faze parne eksplozije The Influence of the Accuracy of Numerical Methods on Steam-Explosion Premixing-Phase Simulation Results Matja` Leskovar - Jure Marn - Borut Mavko Parna eksplozija je pojav, do katerega lahko pride, kadar se zmešata dve kapljevini, od katerih je temperatura prve kapljevine višja od temperature vrelišča druge. Z raziskavami parnih eksplozij se resno ukvarjajo tudi v jedrski tehniki, ker so pri nekaterih scenarijih težkih nezgod v jedrskih elektrarnah, ko pride staljena sredica v stik s hladilno vodo, izpolnjeni pogoji za parno eksplozijo. Največjo pozornost posvečajo mešalni fazi parne eksplozije, saj je moč parne eksplozije odvisna predvsem od premesanosti faz. Da bi ugotovili, kako vpliva natančnost numeričnega reševanja večfaznih enačb na rezultate simuliranj mešalne faze parne eksplozije, smo razvili lasten simulirni program ESE (Evaluation of Steam Explosions). S programom ESE smo opravili številna primerjalna simuliranja različnih preskusov, pri katerih so hladne ali vroče kroglice različnih premerov in materialov spuščali v posodo, napolnjeno z vodo. Vsak preskus smo simulirali z metodo “upwind” prvega reda natančnosti in z metodo velike ločljivosti, ki je drugega reda natančnosti. Rezultate simuliranj smo primerjali med seboj in z razpoložljivimi preskusnimi podatki. © 2000 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: eksplozija pare, tok večfazni, metode numerične, simuliranje) A steam explosion is a physical event which can occur when two fluids are mixing and the temperature of one fluid is higher than the boiling point of the other. Steam explosions are an important area of study in nuclear engineering because the conditions for a steam explosion are fulfilled during some scenaria of severe nuclear reactor accidents, when the molten core comes into contact with the coolant water. Research is mainly focused on the steam-explosion premixing-phase since it determines the extent of the steam explosion. In order to discover how the steam-explosion premixing-phase simulation results are influenced by the accuracy of the numerical methods used for solving the multiphase-flow equations a simulation code named ESE (Evaluation of Steam Explosions) has been developed. With ESE a number of simulations of premixing-phase experiments, where different jets of cold or hot spheres are injected into a water pool, have been performed. Each experiment has been simulated with a first order accurate upwind method and the high-resolution method, which is second-order accurate. The simulation results and available experimental data have been compared. © 2000 Journal of Mechanical Engineering. All rights reserved. (Keywords: steam explosion, multiphase flow, numerical methods, simulation) 0 UVOD Pri mešanju dveh kapljevin, od katerih je temperatura prve kapljevine višja od temperature vrelišča druge, lahko pride do parne eksplozije (PE). Pogoji za PE so tako izpolnjeni pri nekaterih scenarijih težkih nezgod v jedrskih elektrarnah, ko pride staljena sredica po različnih poteh v stik s hladilno vodo. Potek PE lahko glede na dogajanja razdelimo v štiri faze, to je v fazo mirnega mešanja taline in hladiva, fazo sprožitve eksplozije, fazo stopnjevanja in širjenja 0 INTRODUCTION A steam explosion (SE) can occur when two fluids are mixing and the temperature of one fluid is higher than the boiling point of the other. The conditions for a SE are fulfilled during some scenaria of severe nuclear reactor accidents, when the molten core comes into contact with the coolant water. A SE can be divided into four stages: the formation of a premixed region, the occurrence of a trigger, the propagation of the trigger and the expansion. To find gfin^OtJJlMlSCSD 00-9 stran 607 |^BSSITIMIGC M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy eksplozije ter fazo raztezanja in opravljanja dela. Da bi ugotovili, ali in kako je varnost jedrskih elektrarn ogrožena zaradi PE, izvajajo po svetu številne preskuse, povezane s pojavom PE, in razvijajo različne računalniške programe za simuliranje PE [1]. Upoštevajoč aktualnost in pomembnost pojava PE na potek nekaterih scenarijev težkih nezgod v jedrskih elektrarnah smo se lotili razvoja lastnega simulirnega programa ESE (Evaluation of Steam Explosions) [2]. Zaradi zapletenosti pojava PE smo se usmerili le na modeliranje prve faze PE, to je na modeliranje mirnega mešanja taline in hladiva pred izbruhom prave eksplozije. Mešalni fazi PE posvečajo največjo pozornost tudi drugi raziskovalci, saj je moč PE odvisna predvsem od premešanosti faz. Glavni namen naših raziskav je bil spoznati vpliv natančnosti uporabljenih numeričnih metod na rezultate simuliranj in tako ugotoviti, ali bi bilo sedanje simulirne programe mešalne faze PE, ki vsi temeljijo na numeričnih metodah prvega reda natančnosti, primerno nadgraditi z natančnejšimi numeričnimi metodami drugega reda natančnosti [3]. 1 MATEMATIČNI MODEL Vsako od faz p v večfaznem toku taline f, vode w in pare v smo opisali z verjetnostno kontinuitetno: out if and how the safety of nuclear power plants is threatened by SEs, a number of experiments related to SEs are carried out in the world, and different SE-simulation computer codes are being developed [1]. After considering the real posibility and affects of a SE on some severe accident scenaria in nuclear power plants we started to developed our own computer code called ESE (Evaluation of Steam Explosions) [2]. Since a SE is a very complex process we focused only on modeling the premixing stage of the SE, which covers mixing of the molten core debris with water prior to any SE. The SE premixing stage is also the main interest of other researchers because it determines the extent of the steam explosion. The primary purpose of our research was to analyze how the accuracy of the numerical methods used influences the simulation results, and to find out in this way if it would be reasonable to upgrade existing SE premixing-phase simulation codes (which are all based on first-order accurate numerical methods) with second-order accurate methods [3]. 1 MATHEMATICAL MODEL Each phase p in the multiphase flow of melted fuel f, water w and vapor v was described using the probabilistic continuity: dt (apPp) + V-(apPpvp) gibalno: momentum: 3 r dt ( apPpvp ) +V-( apppvp ®vp ) = -apVp+apppg r+apMp + vipG in energijsko enačbo: and energy equation: d dt(apPphp)^(apPp ki smo jih izpeljali s statističnim povprečenjem enofaznih enačb ([4] in [5]). V verjetnostnih večfaznih enačbah (1) do (3) se pojavljajo poleg vseh osnovnih členov, ki jih brez fazne verjetnosti a srečamo tudi v nepovprečenih enofaznih enačbah, še številni dodatni členi, ki jih pridelamo med povprečenjem in jih moramo ustrezno modelirati. Zaradi preglednosti so vsi ti členi združeni v spremenljivke G, M in E. Pri večfaznem toku imajo najpomembnejšo vlogo tisti členi, ki opisujejo medfazni prenos snovi, gibalne količine in energije. Modelirali smo jih v odvisnosti od strukture toka (preglednica 1), ki smo jo določevali s kriterijem fazne verjetnosti taline af in razmerja fazne verjetnosti pare 9v = av/(av +aw) kakor je opisano v [6]. Fazna verjetnost a (r,t) določa verjetnost, da bomo včasu t v točki r naleteli na fazo p. Model večfaznega toka smo prevzeli iz [6], kjer so tudi podane vse konstitutivne enačbe za izračun medfaznega prenosa za vsako izmed hpvp apEp+ G ip p (1), (2) (3), obtained by ensemble averaging the single-phase equations ([4] and [5]). In these probabilistic multiphase equations (1) to (3) different types of terms appear. Basic terms without the phase probability a , which also appear in anaveraged single-phase equations, and a number of additional terms produced during the averaging process which have to be adequately modeled. For the sake of clarity these terms were lumped together in variables G, M and E. In multiphase flow the terms describing the interfacial exchange of mass, momentum and energy are the most important. They were modeled according to the flow regime (table 1) determined by the fuel phase presence probability af and the gas void fraction 5 =a l(a +a ) criteria as described in [6]. The phase presence probability a (r,t) determines the probability that phase p will occur in place r at time t. The multiphase-flow model has been taken from [6], where the interfacial exchange laws for each VH^tTPsDDIK stran 608 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy Preglednica 1. Kartogram tokovnih struktur, ki jih upoštevamo pri medfaznem prenosu (prevzeto po [6]) Table 1. Diagram of flow regimes considered in characterizing interface transfers (taken from [6]) tok okoli razpršenih delcev taline flow over immersed fuel particles af<0,3 mehurčkasti tok bubbly flow Jv <0,3 čepasto-turb. tok churn-turb. flow 0,30,3 talina / fuel para / steam I— upoštevanih struktur toka. Zaradi preobsežnosti konstitutivnih enačb ne bomo navajali, ampak se bomo omejili le na opis glavnih značilnosti večfaznega modela. Za af<0,3 smo obravnavali večfazni tok taline, vode in pare kot gibanje delcev taline v dvofaznem toku vode in pare, za katerega smo določili tokovne strukture na podlagi razmerja fazne verjetnosti pare: Jv <0,3 - mehurčkasti tok, 0,3 0,7 - kapljičasti tok. Za af > 0,3, ko so delci taline že gosto posejani, smo obravnavali večfazni tok kot tok vode in pare skozi porozno plast delcev taline. Medfazni prenos smo računali s pomočjo obrazcev, veljavnih za dvofazni tok, navzočnost tretje faze pa smo upoštevali z zmanjšano koncentracijo površine med preostalima dvema fazama za faktor fij = aj /( a + ak), ki upošteva vpliv dodatne faze k na koncentracijo površine medfazne ploskve faze i za interakcijo s fazo j. 2 NUMERIČNA OBRAVNAVA Diskretizirane verjetnostne večfazne enačbe smo reševali v dvodimenzionalnem valjnem koordinatnem sistemu na premaknjeni mreži z metodo končnih razlik. Da bi spoznali vpliv natančnosti numeričnih metod na rezultate simuliranj, smo parcialne diferencialne enačbe reševali primerjalno z metodo “upwind” (UW) prvega reda natančnosti in z metodo velike ločjivosti (VL - HR), ki je drugega reda natančnosti [7]. 2.1 Metoda velike ločljivosti Dobra lastnost metod, ki temeljijo na diskretizaciji prvega reda natančnosti, je njihova robustnost, njihova glavna pomanjkljivost pa je numerična difuzija. Ta difuzija je posledica napak pri diskretizaciji odvodov in se kaže v glajenju strmih gradientov koncentracije, hitrosti in temperature. Numerični difuziji se izognemo z natančnejšimi diskretizacijami drugega reda natančnosti, pri katerih pa se kot posledica numerične disperzije zaradi napak - voda / water of the considered flow regimes are also described. Since the exchange laws are too extensive they will not be presented here and only the basic features of the model will be described. For af <0.3 we considered the fuel particles immersed in a two-phase gas-liquid flow, whose own flow regimes are defined by the value of the void fraction: Jv <0.3 - bubbly flow, 0.3 < Jv <0.7 -churn-turbulent flow, and Jv > 0.7 - droplet flow. For af > 0.3, as the fuel particles are densely packed, we considered a flow of gas and liquid through a porous bed of fuel particles. The interfacial exchange was calculated using modified two-phase flow correlations. The influence of the third phase was considered with a reduced area concentration between the other two phases by a factor of fij = aj/(aj + ak), which represents the effect of the additional phase k on the area concentration of phase i for its interaction with phase j. 2 NUMERICAL TREATMENT The discretised multiphase equations were solved in the 2D cylindrical coordinate system on a staggered grid with the finite-differences method. To find out how the numerical method’s accuracy influences the simulation results, the partial differential equations were solved, for comparison, using the first-order accurate upwind method (UW) and the high-resolution method (HR), which is second-order accurate [7]. 2.1 High-resolution method A good property of methods based on first order accurate discretisations is their robustness, whereas their main drawback is the numerical diffusion. This diffusion is a consequence of the discretisation error, and reflects in the smoothing of steep gradients of concentration, velocity and temperature. Numerical diffusion can be avoided using second-order accurate discretisations, but as a consequence of numerical dispersion, arising from sec- stran 609 SUMEČ M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy drugega reda natančnosti v rešitvah pojavijo oscilacije v okolici strmih gradientov. Te oscilacije, ki so značilne za vse diskretizacije višjih redov natančnosti, so v splošnem bolj neprijetne od numerične difuzije, ker dajejo nefizikalne rešitve, reševanje enačb pa lahko postane nestabilno. Pomanjkljivostim osnovnih metod prvega in drugega reda natančnosti se izognemo z metodami velike ločljivosti (VL). Metode VL s primernim kombiniranjem diskretizacij prvega in drugega reda natančnosti minimizirajo numerično difuzijo diskretizacij prvega reda natančnosti in eliminirajo oscilacije diskretizacij drugega reda natančnosti, pri tem pa ostanejo natančnosti drugega reda [8]. Glavne značilnosti metode VL si oglejmo na preprostem primeru 1D kontinuitetne enačbe brez izvirnih členov za nestisljivo tekočino, ki se giblje z nespremenljivo pozitivno hitrostjo v > 0 : ond-order errors, oscillations occur in the vicinity of steep gradients. These oscillations, which are characteristic of all higher order accurate discretisations, tend to be more problematic than numerical diffusion since they produce unphysical results and the equation solving procedure may become unstable. High-resolution (HR) methods combine first- and second-order accurate numerical methods in such a way that the weaknesses of both methods - the large amount of numerical diffusion of first-order accurate methods and the oscillations of second-order accurate methods - are suppressed. Despite this HR methods remain second-order accurate [8]. The main characteristics of the HR method will be presented using the simple example of a 1D continuity equation without source terms for an incompressible fluid moving with constant positive velocity v>0: da d(av) ----+— = 0 dt dx (4). Konvekcijski člen zapišemo v diskretizirani obliki s fluksi velike ločljivosti kot: The convection term may be written in a discretised form with HR fluxes as: d[avj dx F i+1/2 F i-1/2 Dx Fhr /2 =F +1/2 F2.red + (1-F +1/2 )F1. 1.red i+1/2 i+1/2 (5), kjer omejitveni faktor F določa delež prispevkov where the limiter F defines the contribution of first- fluksov prvega in drugega reda natančnosti k fluksu and second-order accurate fluxes to the HR flux. The velike ločljivosti. Omejitveni faktor, ki je funkcija limiter, which is a function of the smoothness of the gladkosti rešitve J, definirane kot razmerje sosednjih solution J, defined as the ratio of adjacent solution gradientov rešitve: gradients: mora izpolnjevati določene pogoje [8], če želimo dobiti metodo VL. Pri vseh izračunih, ki smo jih opravili z metodo VL, smo uporabili omejitveni faktor “van Leer” [8]: F (J) flukse prvega in drugega reda natančnosti pa smo določevali z “upwind” F1 red = Fuw oz. metodo Lax-Wendroff (LAX) F2red=Flax [8]. Omenjene razlike med posameznimi metodami si lahko nazorno ogledamo na primeru izračuna gibanja skoka fazne verjetnosti, z a = 1 za xe[0,1;0,3] in a = 0 zunaj tega območja, s hitrostjo v = 1. Na sliki 1 je prikazana izračunana porazdelitev fazne verjetnosti po času 0,2, 0,4 in 0,6 skupaj z analitično rešitvijo (anal). Vidimo, da je numerična difuzija najbolj izrazita pri metodi UW, da LAX metoda povzroča oscilacije, in da metoda VL združuje le dobre lastnosti obeh metod. Ker metoda LAX ustvarja fizikalno nesprejemljive rezultate, ni primerna za simuliranje mešalne faze PE. (6), has to fulfill special conditions [8] if we want to ob-tain a HR method. In all calculations which were per-formed with the HR method, the van Leer flux-limiter function has been used [8]: J\ + J 1 + JI (7), and the first- and second-order accurate fluxes have been calculated with the upwind F1.red = Fuw or Lax-Wendroff (LAX) method F2.red = Flax [8]. The differences between the methods are presented for the example of phase presence probability shock propagation with a = 1 for x e [0.1,0.3] and a = 0 elsewhere, moving with velocity v = 1. In Figure 1 the calculated phase-presence probability distribution after steps 0.2, 0.4 and 0.6 is presented together with the analytical solution (anal). It can be seen that the numerical diffusion is best expressed with the UW method, that the LAX method produces oscillations and that the HR method combines only the good properties of both methods. Since the LAX method produces physically unacceptable results, it is not suitable for SE premixing-phase simulations. VBgfFMK stran 610 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy 1.2 0.8 0.6 0.4 0.2 0.2 0.4 0.6 razdalja / distance 1.2 1 0.8 0.6 0.4 0.2 0 0 0.8 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 razdalja / distance 0.2 0.4 0.6 0.8 razdalja / distance Sl. 1. Gibanje skoka fazne verjetnosti na mreži velikosti 160 točk, računano z metodami UW (zgoraj, levo), LAX (zgoraj, desno) in VL (spodaj) Fig. 1. Phase-presence probability shock convection on a mesh of 160 points calculated with the UW (top, left), LAX (top, right) and HR (bottom) methods 2.2 Tlačna enačba S časovno odvisnimi verjetnostnimi večfaznimi enačbami ((1) do (3)) lahko za vsako fazo posebej izračunamo fazno verjetnost a, hitrost v in entalpijo h. Za tlak p nimamo ustrezne eksplicitne enačbe, saj se pojavlja tlak le implicitno v vektorski gibalni enačbi (2) [9]. Tlak ima pri večfaznem toku pomembno vlogo, saj mora oblikovati hitrostno polje posameznih faz tako, da je ves čas izpolnjena tudi trivialna enačba o vsoti faznih verjetnosti: 2.2 Pressure equation With the time-dependent probabilistic multiphase equations ((1) to (3)) for each phase the phase-presence probability a , velocity vr and en-thalpy h can be calculated. For pressure p there is no explicit equation since it occurs only implicitly in the vector-momentum equation (2) [9]. Pressure has an important role in multiphase flow since it has to influ-ence the velocity field in such a way that the trivial phase-presence probability sum equation is fulfilled all the time: p (8). Tlačno enačbo za večfazni tok smo izpeljali iz vsote kontinuitetne enačbe (1) po vseh fazah ob upoštevanju pogoja (8): The multiphase-flow pressure equation has been derived from the sum of the continuity equation (1) over all phases considering condition (8): an+ 1=Z r rn + Gn+1Dt - V • ( an+1r n+1/2 vn+p1 ) Dt + + ( Dt2 )V- \an+1rn+1/2 1Wpn r 1 (pn+1-pn)drn n+ n+1/2 n+1/2 p r r dpn rn + Gn+1 Dt - V • ( an+1rn+1/ 2 v^+p1 ) Dt + + ( Dt2 )V- \an+1rn+1/2 1Wpn r (9), r Dt an+1dp[v?+p1- rVpn\(pn+1-pn)\Dt gfin^OtJJlMlSCSD 00-9 stran 611 |^BSSITIMIGC 0 1 1 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy kjer je in where vn+1 = vn+1/2 + Dvp = vn+1/2 - Dt (W+1 -Vpn) = v:+p1 - DtVpn+ rn+1=rn+1/2 + Dr rn+1/2 + dr dpn Zaradi preglednosti so indeksi faz izpuščeni, indeks p pa sedaj označuje vpliv tlaka. Tlačna enačba (9), ki upošteva tudi temperaturno raztezanje faz, se inherentno prilagodi naravi toka in je zato primerna za obravnavo vseh vrst večfaznih tokov, tako stisljivih kakor tudi nestisljivih tekočin [10]. Tlačno enačbo smo reševali s kvadratno konjugirano gradientno metodo CGSTAB [10], ki smo jo izbrali zato, ker konvergira hitro tudi za velika razmerja gostot posameznih faz. Hitrost konvergence metode CGSTAB smo testirali za različna razmerja gostot faz preprostega dvofaznega toka (dviganje velikega mehurja) in jo primerjali s hitrostjo konvergence dveh drugih metod, ki se pogosto uporabljata za reševanje tlačne enačbe, to je metode ADI in metode SIP [10]. S slike 2 je razvidno, da najhitreje konvergira uporabljena metoda CGSTAB in da je hitrost njene konvergence skoraj neodvisna od razmerja gostot posameznih faz, medtem ko postaneta metodi ADI in SIP neuporabni pri razmerju gostot, večjem od približno 1000. 100000 10000 1000 100" 10 1------- 1 and (pn+1 -pn ) , rn+1/2 =r(hn+1,pn ) (10) (11). For the sake of clarity the phase indexes are omitted. Index p now denotes the pressure influence. The pressure equation (9), which also considers the phase’s temperature expansion, automatically adjusts to the local nature of the flow and therefore can be applied to all types of compressible or incompress-ible multiphase flow [10]. The pressure equation has been solved us-ing the squared-conjugate-gradient method (CGSTAB) [10], which was chosen since it also con-verges quickly for high phase-density ratia. The con-vergence rate of the CGSTAB method was tested for different ratia of a two-phase flow (rise of a big bubble) and was compared with other proposed pressure-equation solvers: the Alternate Direction Implicit solver (ADI) and the Stones Strongly Implicit Procedure (SIP) [10]. From Figure 2, we observe that the CGSTAB method has the fastest convergence rate, which is nearly independent of the density ratio. In contrast, the ADI and SIP methods become useless if the density ratio is greater than about 1000. 10000 1000 100 10 0.1 10 100 1000 razmerje gostot / density ratio 10000 10 100 1000 razmerje gostot / density ratio 10000 Sl. 2. Hitrost konvergence metod ADI, SIP in CGSTAB za reševanje tlačne enačbe v odvisnosti od razmerja gostot faz Fig. 2. Convergence rate of ADI, SIP and CGSTAB pressure equation solvers for different phase-density ratia 3 SMULIRANJA PRESKUSOV S programom ESE, ki smo ga poprej natančno testirali [3], smo simulirali številne preskuse mešalne faze PE (pregl. 2), ki so jih izvedli na napravi QUEOS [11]. Pri teh preskusih so spuščali valjast curek hladnih ali vročih kroglic različnih premerov in materialov v posodo z vodo, ki je imela temperaturo vrelišča. 3 EXPERIMENT SIMULATIONS Using ESE, which has been extensively tested previously [3], a number of premixing experiments (table 2) performed at the QUEOS facility [11] have been simulated. In these experiments a cylindri-cal jet of cold or hot spheres of variing diameter and material was injected into a pool filled with water at boiling point. VBgfFMK stran 612 1 1 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy Preglednica 2. Vrednosti osnovnih parametrov simuliranih preskusov QUEOS Table 2. Conditions of simulated QUEOS premixing experiments Preskus Experiment material material gostota density kg/m3 premer diameter mm masa mass kg prostornina volume cm3 temperatura temperature K Q01 jeklo / steel 7800 4,7 20 4440 300 Q02 jeklo / steel 4,7 10 2220 300 Q05 ZrO2 5750 5,0 7,0 1830 300 Q09 ZrO2 5,0 7,0 1880 1300 Q06 ZrO2 10,0 7,0 1900 300 Q07 ZrO2 10,0 7,0 1900 1300 Q08 Mo 10200 4,2 10 1800 300 Q04 Mo 4,2 10 1800 1300 Q11 Mo 4,2 5,7 1025 1800 Q12 Mo 4,2 6,9 1183 2300 3.1 Modeliranje preskusov Simuliranja smo opravili v dvodimenzionalnem valjnem koordinatnem sistemu ob predpostavki osnosimetričnosti preskusov. Ta predpostavka je le delno upravičljiva. Začetna porazdelitev kroglic namreč predvsem na čelu curka zaradi sprostilnega mehanizma kroglic, pri katerem se simetrično odpreta dve ravni premični plošči, ki zadržujeta kroglice, ni osno simetrična. To se dobro vidi na posnetkih preskusa Q08, kjer je bil omenjeni pojav najbolj izrazit, iz dveh med seboj pravokotnih smeri (sl. 3). V modelu smo predpostavili, da so kroglice v curku na začetku simuliranja enakomerno porazdeljene. To je prav tako zelo idealizirano dejansko stanje, saj je curek kroglic zaradi načina odpiranja zaklopke, pri “vročih” preskusih pa še dodatno zaradi potiska segretega plina, nekoliko razpršen že pred vstopom v vodo (sl. 4). Ker poudarek naših raziskav ni bil na modeliranju vseh posebnosti posameznih preskusov, ampak predvsem na modeliranju najpomembnejših pojavov med mešalno fazo PE, in to predvsem z namenom spoznati vpliv natančnosti numeričnih metod na rezultate simuliranj, teh in tudi nekaterih drugih posebnih pomanjkljivosti našega modela nismo poskušali odpraviti. Pri simuliranjih smo obravnavali le z vodo napolnjeni spodnji del posode in plast pare nad vodno gladino, vse druge značilnosti preskusov, vključno s padajočim curkom kroglic, pa smo v model vnesli z ustreznimi robnimi pogoji (R.P) (slika 5). Pri “hladnih” preskusih je zgornji rob simulirnega območja segal 10 cm nad začetni nivo vodne gladine, pri “vročih” preskusih pa še 10 cm višje, saj bi sicer pri “vročih” preskusih, pri katerih je zaradi uparjanja dvig gladine vode večji, voda iztekala iz obravnavanega območja. Robne vrednosti hitrostnega polja kroglic in polja fazne verjetnosti kroglic na zgornjem robu simulirnega območja smo 3.1 Experiment modeling The simulations were performed in 2D in the cylindrical coordinate system assuming axial symmetry of the experiments. This presumption is only partially valid; because of the spheres’ release mechanism with two sliding doors holding the spheres opening symmetrically, the initial sphere distribution is not axially symmetric, especially at the jet’s head. This can best be seen in experiment Q08, where the described phenomenon was the strongest, with images taken from two perpendicu-lar sides (Figure 3). In the model a uniform initial sphere distribution in the jet was assumed. This is also a very idealized condition since the jet was, because of the spheres’ release mechanism and due to gas burst durig “hot” experiments, somewhat dis-persed before the water entry (Figure 4). Since the emphasis of our research was not in modeling all the specialties of these experiments, but mainly to model the most important phenomena appearing during the SE-premixing phase with the intention of determining the influence of the numerical meth-ods’ accuracy on the simulation results, no attempt has been made to remove these and other specific drawbacks of our model. In the simulations only the lower part of the vessel filled with water and a vapor layer over the water level have been treated. All other experimental features including the falling-spheres jet have been incorporated into the model with corresponding boundary conditions (B.C.) (Figure 5). In “cold” experiments the simulation region upper boundary ex-tended 10 cm over the initial water level and in “hot” experiments 10 cm higher since otherwise in “hot” experiments, where the water level swell is higher, water would flow out of the simulation region. The sphere velocity and phase-presence probability boundary conditions on the upper edge of the simula- griTM5)tJJ[M]S[]SD 00-9 stran 613 |^BSSITIMIGC M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy L^dU Sl. 3. Posnetka preskusa Q08 s severne strani (osvetlitev od zadaj) in z zahodne strani (črno ozadje) (posnetki vzeti iz [11]) Fig. 3. Images of experiment Q08 from the north side (back-lighting) and from the west side (black background) (images taken from [11 ]) curek vročih kroglic / jet of hot spheres curek hladnih kroglic / jet of cold spheres para / vapor 10 cmT simetrijska os symmetry axis voda water Sl. 4. Posnetki curkov kroglic ob stiku z vodo za “hladen” preskus Q08 ter za “vroča” preskusa Q04, pri katerih je razpršitev kroglic najbolj izrazita, in Q11 (od leve proti desni) (posnetki vzeti iz [11]) Fig. 4. Images of sphere jets at water contact for the “cold” experiment Q08 and for the “hot” experiment Q04, where the sphere dispersion is most expressive, and Q11 (from left to right) (images taken from [11]) R.P./B.C: curek / jet: vz rob(t), arob(t) tlak / pressure: prob(t) 20 cm 100 cm 82 cm Sl. 5. Modela “hladnih” (leva polovica) in “vročih” (desna polovica) QUEOS preskusov Fig. 5. Models of “cold” (left half and “hot” (right half QUEOS experiments določili posredno iz znanih povprečnih vrednosti na vodni gladini (vzi = 5,12 m/s, ai = 0,17). 3.2 Rezultati simuliranj 3.2.1 “Hladni” QUEOS preskusi Na sliki 6 so prikazani rezultati simuliranj značilnega “hladnega” preskusa (Q08), opravljenih z tion region have been determined indirectly from known values at the water level ( vzi = 5.12 m/s , ai = 0.17 ). 3.2 Simulation results 3.2.1 “Cold” QUEOS experiments In Figure 6 the simulation results of the typical “cold” experiment (Q08) performed with the first-order VH^tTPsDDIK stran 614 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy metodo UW prvega reda natančnosti in z metodo HR. Intenzivnost fazne verjetnosti kroglic, ki je definirana kot integral fazne verjetnosti kroglic v vodoravni smeri, je prikazana v kartezičnem koordinatnem sistemu z obrisi vrednosti 0,002; 0,005; 0,01; 0,02; 0,03 itn. Fazna verjetnost vode je prikazana v valjnem koordinatnem sistemu z ekvidistančnimi obrisi vrednosti 0,1; 0,2; 0,3 itn. Številke na robu slik fazne verjetnosti vode pomenijo indekse mrežnih točk (vključno z robnimi točkami) v obeh koordinatnih smereh, številke na robu slik intenzivnosti fazne verjetnosti kroglic pa označujejo oddaljenost v cm od levega oz. spodnjega roba obravnavanega območja. accurate UW method and the HR method are presented. The spheres’ phase-presence probability intensity, defined as the integral of the phase presence probability in the horizontal direction, is presented in the cartesian coordinate system with contours showing 0.002, 0.005, 0.01, 0.02, 0.03, etc. values. The water phase-presence probability is presented in the cylindrical coordinate system with equidistant contours showing 0.1, 0.2, 0.3, etc. values. The numbers on the sides of the water phase-presence probability pictures represent mesh indexes (including boundary points) and the numbers on the sides of the spheres’ phase-presence probability intensity pictures represent the distance in cm from the left and lower boundary of the simulation region. 20 40 60 Sl. 6. Od leve proti desni si sledijo posnetek preskusa Q08 po 0,22 s [11] ter izračuni intenzivnosti fazne verjetnosti kroglic in fazne verjetnosti vode, opravljeni z metodo UW (leva stran) in metodo VL (desna stran) Fig. 6. From left to right are the image of experiment Q08 after 0.22 s [11] and the calculations of the spheres’ phase-presence probability intensity and the water phase-presence probability performed with the UW (left side) and the HR methods (right side) Vidimo, da se rezultati simuliranj kakovostno precej dobro ujemajo s fotografijo preskusa. Če primerjamo izračune, opravljene z metodama HR in UW, pridemo do naslednjih splošnih ugotovitev: - Gradienti fazne verjetnosti so pri izračunih z metodo VL v splošnem večji kakor pri izračunih z metodo UW. Tako je prehod med vodo in paro pri metodi UW zelo razmazan tudi na območjih, kjer bi morali imeti ostro vodno gladino (npr. dvig nivoja vodne gladine, ujetje parnega mehurja). - Pri izračunih z metodo VL se kroglice razpršijo v obliko narobe obrnjene gobe s širšim in bolj topim klobukom. - Pri metodi VL se oblak razpršenih kroglic proti koncu simuliranja popolnoma loči od parnega stebra, kakor je razvidno tudi iz posnetkov preskusov, medtem ko ostane pri metodi UW prehod zabrisan. Na sliki 7 so primerjalno prikazani rezultati izračunov lege čela curka in polovične širine klobuka curka z obema numeričnima metodama skupaj s preskusnimi meritvami. Vidimo, da se pri metodi UW zaradi izrazitejše numerične difuzije in posledične večje The simulation results agree qualitatively quite well with the image of the experiment. The comparison of the calculations performed with the HR and UW methods reveals the following general findings: - Phase-presence probability gradients tend to be higher for calculations performed with the HR method than with the UW method. So for the UW method the transition between water and vapor is also extensively spread in regions where a sharp water-vapor interface should be (water level swell, entrapped vapor bubble). - For HR-method calculations the spheres form an inverted mushroom structure with a wider and more blunt hat. - For the HR method the sphere cloud separates completely from the vapor chimney as seen in the experiments, whereas for the UW method the transition is blurred. In Figure 7 the jet penetration depth and the plume half-width calculated with both numerical methods are presented for comparison with the expe-rimental measurements. The UW method predicts a gfin^OtJJIMISCSD 00-9 stran 615 |^BSSITIMIGC M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy razpršitve čela curka v smeri gibanja napove hitrejše dejansko prodiranje čela curka kakor pri metodi VL. Pričakovali bi, da se bo iz istega razloga povečala tudi širina curka, vendar se izkaže, da delovanje numerične difuzije spremeni tudi geometrijsko obliko glave curka kroglic (sl. 6) in zato se vrtinci v vodi, ki obdaja curek, tvorijo drugače kakor bi se sicer. Glava curka zaradi tega ne zavzame tako izrazite oblike narobe obrnjene gobe in zato je širina curka manjša. Kakor je razvidno s slike 7 ležijo meritve širine curka pri večini preskusov v območju, ki ga omejujeta napovedi obeh numeričnih metod. Razlog, zakaj metoda VL napove nekoliko prehitro širitev curka, tiči verjetno v dejstvu, da curka nismo modelirali dovolj stvarno. faster effective jet-front propagation than the HR method because the numerical diffusion is more expressive and consequently the jet’s head dispersion is greater. One would expect that for the same reason the jet’s plume width would also increase, but it turned out that the numerical diffusion changes the jet’s shape (Figure 6) in such a way that the formation of the water vortices around the jet is different. As a result, the jet’s head does not form such an expressive “inverted mushroom” structure and the jet’s plume width is smaller. From Figure 7, the measured plume width lies mostly between both the calculated curves. The reason why the HR method slightly overpredicts the jet’s spreading could be the unrealistic jet modeling. 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 0.05 --------hr • Q01 0.1 0.15 0.2 čas / time (s) 0.3 0.25 0.2 0.15 0.1 0.05 0 0.25 0.3 0.3 0.25 0.2 0.15 0.1 0.05 0 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 ^^ ^-"^*—" —- ^ ^^*^5s^ ^^ ^:fc^ uw ^*6^ hr ^^>* ¦ Q02 ^S; 0.3 0.25 0.2 0.15 0.1 0.05 0 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 0.3 0.25 0.2 0.15 0.1 0.05 0 0.1 0.2 čas / time (s) 0.3 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 X^ *^* ^-^^ uw hr ¦ Q08 0.3 0.25 0.2 0.15 0.1 0.05 0 0 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 Sl. 7. Primerjava lege čela curka (padajoča krivulja) in polovične širine razpršenega curka kroglic (rastoča krivulja) pri izračunih z metodama UW in VL s podatki preskusov Q01, Q02, Q05, Q06 in Q08 Fig. 7. Comparison of spheres’ jet penetration depth (falling curve) and plume half width (rising curve) between UW- and HR-method calculations and experiments Q01, Q02, Q05, Q06 and Q08 3.2.2 “Vroči” preskusi QUEOS 3.2.2 “Hot” QUEOS experiments Na sliki 8 so prikazani rezultati simuliranj In Figure 8 the simulation results of the typical značilnega “vročega” preskusa (Q11), opravljenih z “hot” experiment (Q11) performed with both numerical VBgfFMK stran 616 0 0 0 0 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy obema numeričnima metodama. Oznake na slikah in vrednosti obrisov so enake kakor pri “hladnem” preskusu na sliki 6. Vidimo, da je parni steber zaradi uparjanja bolj izrazit in ostane dalj časa odprt kakor pri “hladnih” preskusih in se ob navpični simetrijski osi razvije izboklina kot posledica Rayleigh-Taylorjeve (RT) nestabilnosti. Podobne RT nestabilnosti se pojavljajo tudi pri preskusih, vendar zaradi 3D geometrije običajno ne ob osi. Kažejo se kot razpad curka v več manjših curkov (Q09, Q11) ali v oblikovanju ostre konice na čelu curka (Q12). S programom ESE razpada curka v več manjših curkov ne moremo modelirati, ker smo zaradi 2D obravnave v valjnem koordinatnem sistemu omejeni le na osno simetrične pojave. methods are presented. The figure denotations and the contour values are the same as in the “cold” experiment in Figure 6. Because of vapor generation the vapor chimney is more pronounced and remains open for a longer timer than in the “cold” experiments, and a small embossment forms at the vertical symmetry axis as a consequence of the Rayleigh-Taylor (RT) instability. Similar RT instabilities also appear in experiments, but because of the 3D geometry they are usually not near the axis. They are reflected in the jet breakup into a greater number of smaller jets (Q09, Q11) or in the formation of a sharp leading nib (Q12). With ESE the jet breakup into more jets of smaller size cannot be modeled because in the 2D treatment with a cylindrical coordinate system its application is limited only to axially symmetric phenomena. Sl. 8. Od leve proti desni si sledijo posnetek preskusa Q11 po 0,26 s [11] ter izračuni intenzivnosti fazne verjetnosti kroglic in fazne verjetnosti vode, opravljeni z metodo UW (leva polovica) in metodo VL (desna polovica) Fig. 8. From left to right are the image of experiment Q11 after 0.26 s [11] and the calculations of the sphere phase-presence probability intensity and the water phase-presence probability performed with the UW (left side) and HR methods (right side) Na sliki 9 so prikazani izračuni lege čela curka in širine curka skupaj s preskusnimi podatki. Vidimo, da se izračuni tudi pri simuliranju “vročih” preskusov relativno dobro ujemajo z meritvami. Razlike med izračuni z obema numeričnima metodama so podobne kakor pri “hladnih” preskusih, le hitrost prodiranja čela curka je pri metodi VL v drugi polovici simuliranja večja kakor pri metodi UW, kar je drugače, kakor je bilo pri izračunih “hladnih” preskusov. Do tega pride zaradi izbokline (sl. 8), ki nastane ob simetrijski osi zaradi RT nestabilnosti in je bolj izrazita pri simuliranjih z metodo VL. Zanimivo je, da se pri preskusu Q04, pri katerem je opisani pojav razpršitve kroglic med odpiranjem zaklopke največji (sl. 4), rezultati simuliranja z metodo UW bolje ujemajo z eksperimentalnimi podatki kakor rezultati, dobljeni z metodo VL. Poglavitni razlog za tako dobro ujemanje izračunov, opravljenih z metodo UW, je numerična difuzija metode, ki curek kroglic razprši v smeri gibanja, podobno kakor potisk plina po odpiranju zaklopke, ki pa ga nismo modelirali. In Figure 9 the results of the penetration-depth and plume-width calculation are presented together with the experimental measurements. Also for the “hot” experiments the results of the simulations are in relatively good agreement with the experimental data. The differences be-tween the calculations performed with both numerical meth-ods are similar to those of the “cold” experiments. The exception being the speed of the spheres’ front in the second part of the simulation which is higher with the HR method than the UW method which is in opposite to the calculations involving the “cold” experiment. The reason is the formation of the embossment (Fig. 8) at the symme-try axis due to the RT instability, which is more pronounced in the simulations using the HR method. It is interesting that for the Q04 experiment, where the described spheres-spreading effect caused by the gas burst is the strongest (Fig. 4), the UW method gives better agreement with the experimental data than the HR method. The reason may be that the numerical diffusion of the UW method produces a similar spreading of the spheres to the gas burst, which was not modeled. gfin^OtJJlMlSCSD 00-9 stran 617 \^BSSITIMIGC M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 v ^*5^ J? >L uw hr Q09 ^ """ 0.25 0.2 - 0.15 0.1 ¦ o 0 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 ^^— -—"^ ¦^^ uw hr *- ¦ Q04 0.3 0.25 0.2 0.15 0.1 0.05 0 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 v \. uw hr Q07 0.05 0.1 0.15 0.2 čas / time (s) 0.25 ^* ------uw -------hr ¦ Q11 -----r 0.35 0.3 0.25 0.2 0.15 0.1 •»•55,; 0.05 -----0 0.3 0.3 0.25 0.2 0.15 0.1 : 0.05 0 0 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.4 j 0.2 . N 0 -0.2 -0.4 - -0.6 - -0.8 - -1 0 0.3 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 V v uw hr Q12 0.1 ^ 0.25 0.2 0.15 0.1 0.05 0 0.2 0.3 čas / time (s) Sl. 9. Primerjava lege čela curka (padajoča krivulja) in polovične širine razpršenega curka kroglic (rastoča krivulja) pri izračunih z metodama UW in VL s podatki preskusov Q09, Q07, Q04, Q11 in Q12 Fig. 9. Comparison of spheres jet penetration depth (falling curve) and plume half-width (rising curve) between UW- and HR-method calculations and experiments Q09, Q07, Q04, Q11 and Q12 Na sliki 10 so primerjalno prikazani rezultati izračunov pretoka pare in skupne prostornine pare za vse simulirane “vroče” preskuse. Vidimo, da se rezultati kakovostno ujemajo z dinamiko uparjanja med preskusi, kolikostno pa prihaja do odstopanj pri napovedi pretoka pare in pri napovedi skupne prostornine pare. Pri vseh simuliranjih je pretok pare nič, dokler prve kroglice ne pridejo v stik z vodo v posodi, medtem ko je pri preskusih pretok pare znaten že prej. Razlog za to je razširjanje pare, ki je v cevi in delu posode nad vodno gladino in se segreva med prostim padom vročih kroglic, česar nismo modelirali, v kombinaciji z uparjanjem na vodni gladini zaradi toplotnega sevanja. Precej k pretoku pare prispevajo tudi kroglice, ki priletijo v vodo pred glavnim čelom curka. Sevanje v programu ESE upoštevamo le s preprostim znotraj-celičnim modelom, v katerem predpostavimo, da se vsa izsevana energija absorbira že znotraj ene mrežne celice. Ker model ne upošteva medceličnega vpliva sevanja, z njim ni mogoče opisati uparjanja na vodni gladini zaradi sevanja na daljavo. In Figure 10 the comparison of the steam flow rate and total steam volume calculation for all simulated “hot” experiments is presented. The re-sults agree qualitatively with the dynamics of the steam production. However, there are some devia-tions in the absolute value of the steam flow rate. In all simulations the steam flow rate is zero until the first hot spheres enter the water, whereas in the experiments there exists significant steam flow also before. The reason is heating of the atmosphere above the water during the hot spheres’ free fall-dashes which was not modeled-combined with the water surface boiling due to radiation. In addition, the spheres which enter the water before the bulk cause significant early steam flow. Radiation was taken into account in ESE only with a simple intracell model presuming that the entire radiated energy is absorbed within one mesh cell. Since the model does not consider the intercell radiation influence it can-not describe the water surface boiling due to radiation over a distance. At the end of the simulations VBgfFMK stran 618 0 0 0 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy Proti koncu simuliranj pa izračuni odstopajo od meritev prav v nasprotno smer, kakor so na začetku, saj je pri nekaterih preskusih pretok pare že nič (Q07, Q04), medtem ko simuliranja še vedno napovedo nezanemarljivo uparjanje. Razlog, da pri preskusih Q07 in Q04 pretok pare tako skokovito pade na nič, je povezan z dejstvom, da jim pri preskusih ni uspelo vzpostaviti popolnoma enovite začetne porazdelitve temperature vode [11], ampak se je na dnu posode pojavil nekaj decimetrov debel pas z nekaj stopinj hladnejšo vodo, ki v trenutku zaduši uparjanje. 500 450 400 350 300 250 200 150 100 50 0 400 350 300 250 200 150 100 50 0 the calculations deviate from the measurements in the opposite direction to the beginning, since in some experiments the steam flow is zero (Q07, Q04) whereas the simulations still predict significant steam produc-tion. The reason for the abrupt steam flow termina-tion in experiments Q07 and Q04 is related to the fact that in the experiments they could not achieve a com-pletely uniform initial water temperature distribution [11]. So at the bottom of the vessel a layer of a few decimeters with water a few degrees colder was formed, which instantly suppresses boiling. 250 200 150 100 50 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 90 80 70 60 50 40 30 20 10 0 600 500 400 300 200 100 0 120 100 80 60 40 20 0 0 0.05 0.1 0.15 0.2 čas / time (s) 0.25 0.3 900 800 700 600 500 400 300 200 100 0 0.05 0.1 0.15 0.2 0.25 0.3 čas / time (s) 180 160 140 120 100 80 60 40 20 0 0 0.05 0.1 0.15 0.2 0.25 čas / time (s) Sl. 10. Primerjava pretoka pare (nazobčana krivulja) in skupne prostornine pare (gladka krivulja) pri izračunih z metodama UW in VL s podatki preskusov Q09, Q07, Q04, Q11 in Q12 Fig. 10. Comparison of steam flow rate (toothed curve) and total steam volume (smooth curve) between UW- and HR-method calculations and experiments Q09, Q07, Q04, Q11 and Q12 V prvem delu simuliranj se izračuni nastajanja pare z obema metodama skoraj ne razlikujejo, pozneje pa je nastajanje pare pri metodi UW bistveno večje. Razloga za večje nastajanje pare pri metodi UW v drugem delu simuliranja sta vsaj dva, obema pa botruje velika numerična difuzija metode: - Zaradi velike numerične difuzije metode UW je razpršitev kroglic v navpični smeri večja kot pri In the first part of the simulations both meth-ods give similar results for the steam generation, whereas in the second part the steam generation is much higher for the UW method. There are at least two reasons for this, both a consequence of the high numerical diffusion of the UW method: - Because of the high numerical diffusion of the UW method the spheres’ vertical spreading is higher than for the HR method. Since the water’s stran 619 glTMDDC 0 0 0 0 M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy metodi VL. Ker je temperatura vrelišča vode odvisna od tlaka, ta pa se zvišuje z globino vode, je pri večji navpični razpršitvi kroglic večja površina kroglic, ki je v stiku z manj podhlajeno vodo, zato je uparjanje večje. - Ker metoda UW gladi strme gradiente fazne verjetnosti, voda hitreje prodira v parni mehur, ki obdaja vroče kroglice, in tako se uparjanje poveča zaradi večje koncentracije vode v območju vročih kroglic. 4 SKLEP Rezultati simuliranj, opravljenih z razvitim programom ESE, so pokazali, da so izračuni, opravljeni z metodo VL, ki je drugega reda natančnosti, v splošnem boljši od izračunov, opravljenih z metodo UW prvega reda natančnosti. Glede na to, da so odstopanja rezultatov simuliranj od preskusnih meritev zaradi pomanjkljivosti uporabljenih fizikalnih in matematičnih modelov istega velikostnega reda, kakor so razlike med izračuni, opravljeni z metodama prvega in drugega reda natančnosti, je s tega vidika uporaba natančnejših metod drugega reda natančnosti za sedanje modele mešalne faze PE vprašljiva. Če pa upoštevamo dejstvo, da so sedanji programi za simuliranje mešalne faze PE v razvojni fazi, saj jih, da bi dobili boljše ujemanje izračunov z eksperimentalnimi podatki, še vedno dopolnjujejo z različnimi fizikalnimi in matematičnimi modeli, je uporaba natančnejših numeričnih metod priporočljiva. Z uporabo natančnejših numeričnih metod se namreč zmanjša nezanesljivost simuliranja zaradi numeričnih napak in je tako iz rezultatov simuliranj lažje razbrati vpliv in ustreznost fizikalnih in matematičnih modelov. Dejstvo, da je nekatere opazovane fizikalne količine pri nekaterih preskusih bolje napovedala manj natančna metoda UW, nas ne sme zavesti. Dober simulirni program namreč ne sme temeljiti na naključjih, ko v nekaterih okoliščinah nenatančnost numerične metode odtehta pomanjkljivosti fizikalnega in matematičnega modela. Numerične difuzije metode UW, ki včasih sicer ugodno vpliva na rezultate, ne moremo nadzirati in je odvisna tako od časovnega koraka kakor od mrežne razdalje. Če se že izkaže, da je model zaradi neupoštevanja nekaterih fizikalnih pojavov (odboj kroglic, turbulentno mešanje itn.) premalo difuziven in omenjenih pojavov iz različnih razlogov ne moremo ali nočemo modelirati, je veliko bolje, da difuzijo nadzorovano vložimo v matematični model z dodatnimi difuzijskimi členi, kakor pa da upamo, da bo za to v pravšnji meri poskrbela numerična difuzija manj natančne metode prvega reda natančnosti. ^BSfiTTMlliC | stran 620 boiling point temperature depends on pressure, which increases with water depth, at higher vertical spreading the spheres’ surface area, which is in contact with less subcooled water, is greater, as is the steam generation. - Since the UW method smoothes steep phase-presence probability gradients, water penetrates faster in the depleted steam region around the hot spheres, and so the steam generation increases because of the higher water concentration in the hot-spheres region. 4 CONCLUSION The simulation results performed with the developed code called ESE showed that the calcula-tions performed with the HR method, which is sec-ond-order accurate, are in general better than the cal-culations performed with the first-order accurate UW method. The deviations between the simulation re-sults and the experimental measurements are, because of deficiencies in the physical and mathematical models, of the same order as the differences between the calculations performed with the first- and second-order accurate numerical methods. From that point of view the use of second-order accurate numerical methods for existing SE premixing-phase models can-not be justified. But if we take into account that exist-ing SE premixing-phase simulation codes are in the development stage, since they are still supplemented with different physical and mathematical models to get a better agreement of calculation results with ex-perimental data, the use of more accurate numerical methods is recommended. With the use of more ac-curate numerical methods the simulation uncertainty due to numerical errors is reduced and so the influ-ence and adequacy of physical and mathematical models can be more easily determined from the simulation results. The fact that some physical quantities for some experiments were better predicted by the less accurate UW method should not mislead us. A good simulation code must not be based on chances, where in some circumstances the numerical method’s inaccuracy compensates for the models deficiencies. The numerical diffusion of the UW method, which sometimes positively influences the results, cannot be controlled and depends on both the time step and the mesh size. If it turns out that the model is not sufficently diffusive because some physical phenomena (spheres repulsion, turbulent mixing, etc.) are not considered, and we are not able or do not want to model these phenomena for different reasons, it is better to intro-duce diffusion in the mathematical model controlled with additional diffusive terms, than to hope that the numerical diffusion of the first-order accurate method will take care of it. M. Leskovar - J.Marn - B. Mavko: Vpliv natan~nosti - The Influence of the Accuracy 5 LITERATURA 5 REFERENCES [I] Turland, B.D., G.P. Dobson, G.P. (1996) Molten fuel coolant interactions: A state of the art report. European Commission. Nuclear Science and Technology, Luxembourg, ISSN 1018-5593. [2] Leskovar, M., B. Mavko (1998) ESE a 2D compressible multiphase flow code developed for MFCI analysis - code validation. Nuclear Energy in Central Europe ‘98, Čatež, Slovenija, Proceedings, 319-326. [3] Leskovar, M. (1999) Model drugega reda natančnosti za opis večfaznega mešanja tekočin. Doktorska disertacija, Fakulteta za matematiko in fiziko, Univerza v Ljubljani. [4] Leskovar, M., J. Marn (1996) Evaluation of steam explosions: Probabilistic approach. Zeitschrift für Angewandte Mathematik und Mechanik, Vol. 76, 291-292. [5] Marn, J., M. Leskovar (1996) Simulation of steam explosion premixing phase using probabilistic multiphase flow equations. Fluid Mechanics Research, Vol. 22 (1), 44-55. [6] Amarasooriya, W.H., TG. Theofanous (1991) Premixing of steam explosions: A three-fluid model. Nucl Eng& Design, Vol. 126, 23-39. [7] Leskovar, M., B. Mavko, J. Marn (1997) High resolution method applied to premixing phase of steam explosion. Transactions of ANS, Albuquerque, New Mexico, USA, Vol. 77, 433-434. [8] Leveque, R.J. (1992) Numerical methods for conservation laws. Lectures in Mathematics, ETH Zürich, Birkhauser Verlag, Basel. [9] Leskovar, M. (1998) ESE a 2D compressible multiphase flow code developed for MFCI analysis - code description. Nuclear Energy in Central Europe ‘98, Čatež, Slovenija, Proceedings, 311-319. [10] Ferziger, J.H., M. Peric (1996) Computational methods for fluid dynamics. Springer-Verlag, Berlin, Heidelberg. [II] Meyer, L., G. Schumacher (1996) QUEOS, a Simulation-experiment of the premixing phase of a steam explosion with hot spheres in water, base case experiments. FZKA Report 5612, Forschungszentrum Karlsruhe. Naslov avtorjev: dr. Matjaž Leskovar profdr. Borut Mavko Institut “Jožef Stefan” Odsek za reaktorsko tehniko Jamova 39 1000 Ljubljana prof.dr. Jure Marn Fakulteta za strojništvo Univerze v Mariboru Smetanova 17 2000 Maribor Authors’ Addresses: Dr. Matjaž Leskovar ProfDr. Borut Mavko “Jožef Stefan” Institute Reactor Engineering Division Jamova 39 1000 Ljubljana, Slovenia ProfDr. Jure Marn Faculty of Mechanical Eng. University of Maribor Smetanova 17 2000 Maribor, Slovenia Prejeto: Received: 12.6.2000 Sprejeto: Accepted: 20.12.2000