Modelling of core melt scenarios in nuclear reactors

Modelling of core melt scenarios in nuclear reactors

6

Arun K. Nayak*,†, Parimal P. Kulkarni‡, Pradeep Pandey*, Sumit V. Prasad* *Reactor Engineering Division, Bhabha Atomic Research Centre, Mumbai, India, † Homi Bhabha National Institute, Mumbai, India, ‡Bhabha Atomic Research Centre, Mumbai, India

Nomenclature d g h h hbed hfg hfs hfus hin hlv k l p q00 q000 r S t tcr u v w z A Ach Cp Dn E F G

particle diameter (m) gravitational acceleration (m/s2) heat transfer coefficient (W/m2 K) enthalpy (kJ/kg) bed height (m) latent heat of vaporization (J/kg) liquid saturation enthalpy (J/kg) latent heat of fusion (J/kg) inlet enthalpy (J/kg) latent heat of vaporization (J/kg) thermal conductivity (W/m K) crack length (m) pressure (Pa) heat flux (W/m2) volumetric heat generation rate (W/m3) radial direction (m) source term (W/m3) time (s) crust thickness (m) displacement in radial direction (m) velocity (m/s) displacement in axial direction (m) axial direction (m) cross section area (m2) area of openings (m2) specific heat capacity at constant pressure (J/kg K) diameter of nozzle (m) Young’s modulus of elasticity (N/m2) interfacial drag force per unit volume (N/m3) mass flux (kg/m2 s)

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment https://doi.org/10.1016/B978-0-08-102337-2.00006-7 © 2019 Elsevier Ltd. All rights reserved.

488

Gr J K Lb M Nn Nu Pr Q T V Vp

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Grashof number (dimensionless) superficial velocity (m/s) permeability (m2) length of the opening (m) mass (kg) number of inlet nozzles Nusselt number (dimensionless) Prandtl number (dimensionless) volumetric heat generation (W/m3) temperature (K) volume (m3) volume in porous region (m3)

Greek letters α αl β ε εs γ γs η κ μ ν ρ σ σb σ max σ mc σ th σ ST σ1 ξ ψ τ Γ

void fraction linear thermal expansion coefficient (K1) liquid fraction of corium bed porosity emissivity density of openings per unit area (m3) surface energy (N/m) bed passability (m2) bed permeability (m2) dynamic viscosity (Pa s) Poisson’s ratio (dimensionless) density (kg/m3) normal stress (N/m2) bending stress (N/m2) fracture stress(N/m2) surface tension between melt and steam (N/m) thermal stress (N/m2) surface tension (N/m) maximum stress (N/m2) critical Taylor wavelength (m) number of eruption channels shear stress (N/m2) vapour generation rate (kg/s)

Subscripts c c d f fb g

crust critical dryout fluid film boiling gas

Modelling of core melt scenarios in nuclear reactors

i in l m nb p r rad rel s sat sat st sub sup v w z 0

6.1

489

interfacial inlet liquid melt nucleate boiling particle radial radiation relative solid saturation saturated steam subcooled superheated vapor water axial inlet

Introduction to severe accident phenomena and progression in water-cooled reactors

The water-cooled reactors have engineered safety systems to restore the reactor to safe states in the event of any accidents. Most of the safety systems have diverse and redundant features, for example, the shut down system (SDS 1 and 2), emergency core cooling system (ECCS), double containments, etc. The purpose is to enhance the availability of these safety systems so that failure of one of them does not lead to progression of accidents. Besides, the reactors have several multiple barriers to arrest the release of activity if a hypothetical accident progresses with the failure of the safety systems. In spite of all these, an accident of further low probability can be cast beyond the acceptable design basis envelope, which can progress to severe accident. Severe accidents are beyond design basis accidents, which arise due to multiple failures of safety systems and lead to significant core degradation [1]. Whatever could be the initiating events, a severe accident occurs only when the ratio of heat generation rate to heat removal rate exceeds a certain threshold. The fuel is unable to get cooled and gets damaged leading to its melt down. The severity of such an accident depends on the nature and extent of core damage during the accident progression. Of course, these accidents are very low probable ones, but certainly they have high consequences in terms of risk to the plant personnel and public. In view of this, most of the newgeneration reactors are designed with built-in severe accident management strategies. The nuclear utilities also plan to extend or build additional safety features (SAMG) in existing plants in order to enhance the defence in depth of these plants so as to cope with such severe accidents.

490

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Main initiating events, which can progress into severe accident, can be postulated as (i) an unmitigated prolonged station black out (SBO), (ii) LOCA with loss of emergency core cooling (ECC), and (iii) reactivity-induced accidents. LWRs have two independent and fast-acting shutdown systems; hence, reactivity-induced transients without SCRAM are not envisaged. In case of an unmitigated SBO or LOCA with loss of ECCS, the reactor SCRAM occurs. However, the decay heat is a concern in absence of safety systems. If the safety systems fail, water in the core is boiled off by the decay heat and the reactor water level gradually decreases. The steam generated in the core is relieved via the safety relief valves (SRV) into the containment; if an operator recognizes a drop in the reactor water level, depressurization of RPV and low-pressure injection of water from alternative water injection system are expected to keep reactor water level high enough to cool the core. However, if these fail, reactor water level falls below the top of active fuel and fuel is exposed above the water surface causing temperature to rise in fuel pins. If fuel temperature rises, the channel box and control rod temperatures rise due to the radiation from the fuel. At high temperatures, metal-water reaction occurs and hydrogen is generated by the steam oxidation reaction of the Zr in the cladding and used in other components. At around 1500 K, the heat release rate due to oxidation of metallic zirconium into ZrO2 is comparable to decay heat, and this rate increases with rise in temperature. In fact, at above 1800 K, this rate is almost 10 times higher [1]. This causes rapid rise in fuel and core structures to temperature of melting point of Zircalloy (2100 K) and even above 3000 K to cause melting of fuel. The oxidation rate is limited by availability of steam in the core; diffusion of steam through the hydrogen blanketed cladding; diffusion rate of oxygen through the already-oxidized zirconium on the surface of clad material. The total amount of hydrogen released is mostly dependent on the oxidation of zircalloy, although oxidation of core structures such as stainless steel and B4C can contribute to hydrogen release. Oxidation of B4C may contribute to release of other gases such as CO, CO2, aerosols in addition to hydrogen. Deflagration and detonation of hydrogen may occur if this gas mixes with air above certain concentrations. The loss of core geometry can occur over a period of minutes to hours as experienced in TMI-2 and Fukushima LWRs, starting from temperature above 1500 K to 3000 K [1]. Chemical reactions play a key role in formation of low melting point alloys, responsible for core degradation. For low-pressure accident conditions, Zircalloy clad starts to balloon and rupture at temperatures between 1000 K and 1200 K depending on the internal pressure of fuel rods due to fission gases and mechanical behaviour of clad material. For high-pressure accident conditions, the failure of Zircalloy clad is delayed and occurs at temperature above 1500 K due to clad collapse on fuel. Chemical interactions between Zircalloy clad with other fuel materials can cause clad failure locally due to formation of low melting temperature alloys. At temperatures between 1500 K and 1700 K, depending on the intensity of transient, chemical interactions occur among fuel clad, structural materials, and control rod materials, i.e. Fe-Zr, B4C-Zr, Fe-B4C, etc. This causes early liquefaction and relocation of core structures. Re-criticality becomes a major concern when control rod materials are segregated from fuel with re-flooding. The whole relocation process may take only a few minutes to several hours depending on the accident progression. Above

Modelling of core melt scenarios in nuclear reactors

491

2000 K, melting of Zircalloy clad can happen and may fall into lower core regions. For fast transients with heating rates in excess of 0.3–0.5 K/s, the zircalloy clad can melt and drain into lower core regions with very short time delay preventing dissolution of UO2 with molten Zircalloy [1]. The high temperature of the fuel pellet causes fission gas release from the pellet to the pellet-cladding gap. Due to eutectic reactions, the control rod (B4C) and fuel pellet (UO2) melts and moves to the lower regions of the core due to gravity. Mixture of molten fuel and core structure materials is called corium, and solidified corium is called crust. With time, the core support plate fails due to core melt that moves to the lower regions along with core melt during relocation. If the core support plate is damaged, corium influxes to the lower plenum as a jet via the damaged opening of the core support plate. In the lower head, water may be present during the relocation of core debris as in TMI-2 accident. In such situation, molten high-temperature fuel coolant interactions occur resulting in a steam explosion, or the hot melt may be solidified as chunks with broken debris and cause violent boiling. The occurrence of high yield steam explosion is very dangerous, which has the potential for mechanical damage of reactor structures or/and even failure of lower head due to addition of large amount of energy release in very small time period. This energy is maximized if lower head does not fail and may induce high loads on the bolts that hold the upper head causing their rupture. The upper head is accelerated upward to the containment roof like a missile and may cause failure of containment wall. This mode of failure of containment is known as ‘alpha mode failure’. At the same time, hydrogen is produced due to the steam oxidation reaction at the surface of the corium with water. Thus, the corium accumulated in the lower plenum forms the layer of particulate corium (particulate debris bed), metal layer with the low density, and corium, which is surrounded by the solidified corium (crust). In this condition, heat transfer between the RPV lower head and the adjacent crust is impeded by the gap between them. The gap is expanded by creep deformation of the RPV lower head by over-heating from the crust. As a result, corium cooling is maintained due to the surrounding water and temperature rise in the RPV lower head is suppressed by the flow of cooling water into the gap formed between the deformed RPV lower head and the crust. This creates a continuous heat removal path from the molten corium and the debris located in the lower head. With time, the water in the lower plenum is depleted. Radiation heat transfer becomes the dominant mode of heat removal to cool the corium. Since this mode of heat transfer is poor, the temperature of the corium, the RPV lower head, instrumentation tubes and housing, etc. continuously increases. In TMI-2 accidents, it was observed that in the lower head, dense layers of oxidic materials involving UO2, ZrO2 and metallic species of U, O, Zr and Fe mixtures such as steel were present [1]. Because of poor heat removal, the structures in the lower plenum may be damaged causing failure of the RPV lower head resulting into discharge of corium through the failed opening. The failed opening is gradually expanded by the high-temperature outflow of the corium (ablation effect), and corium flow rate may increase. In case of dry lower head, which is unlikely, molten corium at high temperatures (more than 2500°C) strikes the lower head having a melting temperature of 1500°C

492

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

(since made of steel). There are possibilities that the hot corium may melt through the lower head or the lower head may fail by creep. However, many tests suggest that since steel is a good conductor of heat, the molten corium soon form the crust, which behaves as insulating material due to its poor thermal conductivity. Thus, the hot corium does not directly interact with the steel of lower head. However, direct impact of high superheat molten corium may ablate the lower head and melt through the lower head. Several tests conducted demonstrate that it is possible to cool the lower head and retain the melt by cooling the lower head externally as that adopted for VVER-440, AP600, AP1000, APR-1400, etc. With failure of lower head, the molten corium moves into the reactor cavity. If RPV failure occurs at high-pressure condition (if depressurization fails), high-pressure melt ejection (HPME) could occur. The melt will thus be forcibly ejected out of the vessel into the reactor cavity, and the corium jet is torn by high-velocity gas flow and fragmented corium droplets are generated. Because the corium droplets have a large surface area per volume, instantaneous heat transfer from the corium to the gas increases. Steam oxidation reaction at the surface of the corium droplet may occur. Combustion of hydrogen produced during melt dispersion into the containment may also contribute to containment pressurization. As a result, containment may be pressurized and early failure may occur due to Direct Containment Heating (DCH) resulting in release of large amount of fission products and risk to public is very high. If water is present in reactor cavity (collected from several sources during accident progression or some LWRs keep large amount of water in lower head), steam explosion is a concern due to interaction of molten corium with water. The in-vessel and ex-vessel steam explosion phenomena are the same, except the yield could be different depending on the thermal hydraulic conditions such as water depth, melt velocity, subcooling of water, etc. and the melt composition prevailing in in-vessel and ex-vessel scenarios (ex-vessel could have more steel and other structural materials in corium). ‘Stratified steam explosion’ may occur if melt is flooded with water. These explosions are less energetic than those with melt jet falling in water, because of less amount of melt-water interaction and large steam space for expansion. The fuel coolant interaction in the reactor cavity forms a particulate debris bed, which is still generating decay heat; the size and porosity of debris depends on various factors such as melt composition, height of water pool, melt velocity, occurrence of steam explosion, etc. The coolability of such debris bed is a major issue; otherwise, it results into re-melting and further progression of the accident. If containment is dry or when the debris is dried-up with evaporation and boiling of water present in the containment cavity, molten corium concrete interaction (MCCI) occurs. The MCCI ablates the concrete and produces different gases such as CO, CO2, water vapour, which pressurize the containment. The amount of gases released depends on type of concrete such as siliceous concrete or lime stone concrete. With heat loss from top of melt pool due to radiation heat transfer, a crust is formed at top, which behaves like an insulation for the containment and preventing rapid containment heating. The heat generating melt pool, however, has the possibility of seeping through the ground after ablation of basemat depending on the thickness of the basemat.

Modelling of core melt scenarios in nuclear reactors

6.2

493

Current status of application of CFD codes for severe accidents

6.2.1 Applicability of CFD to simulate severe accident phenomena Since the occurrence of severe accident in TMI in 1979, the emphasis on numerical simulation of severe accident phenomena has been limited to development of integral system codes (or software systems) to simulate the entire core melt accident, from the initiating event to radioactive release outside the containment with simplistic approach and with several assumptions in modelling. Simulation of entire phenomena using CFD is still not a reality, because severe accidents phenomena involve multiphase, multicomponent heat, mass, and momentum transfer with phase change with a wide variation of properties during the progression of accidents, which is beyond the capability of existing CFD codes. Only recently, with more and more information available from experiments and development of computational tools, use of CFD for specific phenomena during a given accident phase has taken place. CFD codes are generally used to simulate part of an accident scenario such as hydrogen distribution in containment, natural convection of core melt in calandria vessel or lower head of PWR vessel, cooling of debris bed formed in water-filled containment cavity, etc. Their main goal has been to provide a more detailed understanding of physical phenomena in severe accidents, which is also partly used to validate the predictions of system codes. Table 6.1 gives the scope of CFD codes to simulate different phenomena of severe accident. While the approach for simulation of core damage scenario is mostly by system codes, very little efforts have been done to simulate the fuel melting by CFD. Similar is the case for fission product transport and aerosol behaviour in containment. The major focus is on lumped parameter approach. Although Martı´n-Fuertes et al. [2] did simulation of aerosol behaviour with MELCOR-CFX combination and presented detailed calculations concerning aerosol deposition in the steam generator tube. The simulations were obtained by means of an in-house code application, named COCOA, as well as with CFX computational fluid dynamics code, in which several models for aerosol deposition were implemented and tested. For in-vessel retention of corium in lower head of the PWRs, a lot of CFD simulations have been carried out in the past. When the melt falls in the lower head, due to the high-Rayleigh-number, turbulent natural convection occurs in the volumetrically heated corium melt pools. It is difficult to capture the heat transfer characteristics in detail. Very recently, Zhang et al. [3] did Large Eddy Simulation (LES) to analyse the complicated turbulent and multi-phase flows. They employed Wall-Modelled LES (WMLES) method for the simulation of three typical corium pool heat transfer experiments: BALI, LIVE, and COPRA experiments. The numerical simulation for the transient formation of two-layer corium pool was performed with prototypic corium to obtain the heat transfer data in reactor situation. There was a distinctive gap between the temperatures in two layers. The heat flux showed a decrease at the two-layer interface and a huge jump at the metal layer side. The thin metal layer

494

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Table 6.1

Scope of CFD modelling in severe accident phenomena

Sr No

Severe accident phenomena

1 2

Core damage Fission product release and transport In-vessel retention

3

Specific phenomena involved Fuel heat-up, candling, and melting Solid-gas interaction: Aerosol behaviour Heat transfer from heat generating melt to water l

l

l

4 5

Vessel failure Steam explosion

Thermo-mechanical analysis of vessel Corium fragmentation in water and heat transfer l

l

l

l

6

Debris bed coolability

7

Molten corium concrete interaction (If no core catcher is present) Corium coolability in ex-vessel core catcher

8

9 10

Iodine chemistry in containment Hydrogen combustion and containment loading

Turbulent natural convection of melt pool with heat generation and phenomenon of solidification Boiling on external surface Bubble crowding and DNB

Melt jet breakup Fragmentation Heat transfer Vapour explosion

Natural convection in heat-generating debris bed with single phase and boiling heat transfer, pressure drop, counter current flooding limit, local dryout Interaction of molten core with concrete basemat, concrete ablation, gases generation, changes in corium properties, phase inversion Molten core sacrificial material ablation, hightemperature thermochemical reactions, phase inversion, heat transfer during top flooding, bottom flooding and indirect cooling, crust formation and rupture, water ingression, steam eruption, corium fragmentation and coolability Iodine transport, chemical reaction, and scrubbing Hydrogen distribution, hydrogen combustion, containment loading

on the top at the early stage of the two-layer formation is a concern and may threaten the vessel integrity. In this regard, Horvat and Mavko [4] performed the transient calculations with the CFX 5.7 fluid dynamic software based on the SST turbulence model, and the results showed good comparison with the experimental results from Asfia and Dhir [5]. Buck et al. [6] used the DIVA module of the ASTEC V1.3 code and the melt pool model of

Modelling of core melt scenarios in nuclear reactors

495

ATHLET-CD late-phase modules to simulate the LIVE experiment. Tran and Dinh [7] developed a phase-change effective convectivity model (PECM) to simulate the dynamics of the melt pool heat transfer in the lower head. The PECM method bypassed the Naiver-Stokes equations, and only to solve energy conservation equation using the Fluent solver. The convective terms were included in the energy equation to describe the turbulent natural convection heat transfer (Tran and Dinh, [7]). The effectiveness of in-vessel retention (IVR) by external reactor vessel cooling (ERVC) strongly depends on the critical heat flux (CHF). As long as the local CHF does not exceed the local heat flux, the lower head of the pressure vessel can be cooled sufficiently to prevent from failure. For determination of CHF, the approach so far has been by empirical correlations based on experiments only. Until recently, Zhang et al. [8] carried out CFD simulation to investigate the CHF of ERVC. This simulation was performed by a CFD code Fluent coupled with a boiling model by UDF (user-defined function). The experimental CHF of ERVC obtained by State Nuclear Power Technology Research and Development Centre (SNPTRD) was used to validate this CFD simulation. Buck et al. [6] designed LIVE experimental facility for studying severe accident in PWRs. This was aimed at investigating the prototypic core melt behaviour in the lower plenum of RPV and influence of water cooling on vessel outer surface. Experiments were conducted under varying conditions such as initial melt temperature (maximum 350°C), initial external air/water cooling, central melt pouring, noneutectic and eutectic melts, etc. for in-vessel retention. Experiments were performed in a 1:5 scaled semi-spherical lower head for PWR. To simulate the corium melt, a noneutectic binary mixture of NaNO3 and KNO3 was used. The information obtained from the LIVE experiments includes the melt temperature evolution during different stages of the test, the heat flux distribution along the reactor pressure vessel wall in transient and steady-state conditions, the crust growth velocity, and the influence of the crust formation on the heat flux distribution along the vessel wall. Experimental results of LIVE test were used for validation of calculations and applying models with stand-alone DIVA module of the ASTEC V1.3 code. Tests were conducted in LHF [9] to examine the lower vessel deformation and failure behaviour due to corium relocation and pool formation. Wideranging experimental results have been obtained regarding heat flux distribution, crust thickness, and influence of initial melt pool temperature on the coolability. Caroli et al. [10] investigated the thermal and mechanical behaviour of a PWR vessel in consequence of a severe accident with core melt and flooding of the reactor pit by refined finite element analyses. Willsch€utz et al. [11] estimated the failure time of the RPV of an LWR by numerical analysis during severe accident scenarios and validated the predictions with FOREVER results. In FOREVER [12], tests were carried out on scaled vessels lower head to understand the molten pool convection, the timing and mode of the vessel failure. The vessel instrumentation included thermocouples measuring the temperatures at inner wall, outer wall and inside melt pool, a pressure transducer, and a set of linear position transducers. During the progression of the accident, if external reactor vessel cooling fails, the corium falls into a pool of water and may cause fuel coolant interaction. A lot of experimental and numerical studies were carried out to study the FCI phenomena. The first

496

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

step of FCI is jet breakup. The mechanism for jet breakup during FCI was studied by Thakre et al. [13] using FLUENT 14.0. The jet breakup mechanism was found to be dependent on the injection velocity, which thus affects the size distribution of the droplet. Also, the jet break up length was found to be a function of coolant density ratio. Analysis of identifying potential factors, which may affect the jet instability, determining the scaling laws, and predicting the jet behaviour for severe accident conditions was carried out by Dinh [14,15] using MELT 3D. The major factors affecting jet dynamics were found out to be relative (jet-pool) velocity and the transient flow pattern, the density of melt, and the effective density of coolant. Metal jet breakup, cooling and solidification in water was studied by Zhou et al. [16]. Phase change during jet penetration and breakup was successfully modelled by integrating solid liquid phase change in MH code. Potential effects of prototypic FCI conditions on film boiling heat transfer and flow structure were studied by Dinh [14,15] using Local Homogeneous slip model in NARAL code and CFX-F3D code. It highlighted the significant uncertainties in predicting film boiling heat transfer associated with physical properties of the emitting surface and vapour at the prototypic FCI conditions. 3-D simulation of fuel coolant interaction while plunging jet penetration into a denser liquid pool was carried out by Park [17] using Rigid body Dynamics coupled MPS (RD-MPS). A new model comprising rigid body dynamics and moving particle semi-implicit method was incorporated and the model could very well predict the jet penetration depth of plunging water jet in denser fluid pools. Estimation of pressure loadings on the typical PWR cavity structures during an ex-vessel steam explosion was studied by Cizelj et al. [18] using CFX-5.7.1. The CFD simulations showed that the pressure loads on the cavity walls were the highest during the pre-mixture high-pressure relief, when also pressure shocks propagate to the walls. This process is completed within the first 3–5 ms. Molten droplet preconditioning (deformation and pre-fragmentation) during FCI was studied by Thakre et al. [19] using ANSYS-FLUENT. It inferred that VOF method could capture Rayleigh Taylor instability of two-phase flow. The deformation rate of the droplet increases with increasing Weber number. Deformation of droplet of different sizes was also studied under pressure pulse triggered steam explosion. Post TMI accident, coolability of heat-generating debris bed has become an important issue to be addressed. While heat transfer in porous medium is widely studied in open literature, in severe accident context, the main importance is prediction of dryout heat flux and two-phase pressure drop inside volumetrically heated porous medium that has not been studied in deep. To achieve long-term coolability of the configuration, all evaporated water has to be replaced by water inflow due to natural forces. At the same time, the produced steam must escape the porous structure driven by buoyancy forces. If the heat generation is too high and the coolant is unable to take out the heat, the bed reaches dryout condition where the temperature of the bed increases sharply. The heat flux corresponding to this is called as dryout heat flux (DHF). If there is an excessive pressure drop across the bed, water may not reach the entire length of the bed with the available gravity head and dryout may occur. Hence, in order to study the coolability of the debris bed, it is very important to know the DHF and pressure drop characteristics of the bed.

Modelling of core melt scenarios in nuclear reactors

497

The hydrodynamics and heat transfer behaviour in debris bed is complex one influenced by many factors like l

l

l

l

l

l

the mean size of the particles, porosity, which is a function of the size, and the shape of the particles, which constitute the debris bed, operating conditions such as water entry from the top or the bottom of the debris bed, water temperature, magnitude of noncondensable gases generated during MCCI, and spatial distribution of the bed porosity.

These dependencies make modelling a complex issue. While there have been many experiments followed by empirical correlations based on Ergun equation, solution of detailed mass, momentum and energy equations for prediction of dryout heat flux has not been reported at large. One such attempt was made in MC3D-REPO code [20]. MC3D-REPO is a special application for multiphasic thermal-hydraulic calculations in porous media of the more general software MC3D (for multicomponent-3D) [21]. As it uses the usual assumption of thermal equilibrium in the Representative Elementary Volume, it can be used for nuclear debris bed cooling for the prediction of dry out. MC3D uses mixture momentum equation when some components are assumed to have the same velocity; and/or for energy equation, if some components are assumed to have the same temperature; and only one basic pressure field for all the components. Another example is the WABE-2D model [22] aimed at solving the problem of coolability of degraded core material during a severe accident in a light water reactor (LWR) to study the transient boil-off and quenching behaviour of debris beds. It has been developed in the frame of the KESS code system, which is considered to describe the processes of core heat up, melting, degradation, and relocation to the lower plenum as well as the subsequent behaviour. WABE-2D was later upgraded to MEWA-3D to consider the melt water interaction also in complex geometries. The detailed modelling has been presented in Section 6.5. The COUPLE model [23] of SCDAP/RELAP5 is a simplistic model, which calculates the heat-up of the debris and surrounding structures in the reactor lower head vessel after slumping of the core material to the lower head occurs. This model takes into account the decay heat and internal energy of the newly fallen or formed debris and then calculates the transport equation by conduction of this heat in the radial and axial directions to the wall structures and the water surrounding the debris. The most important use of this model is to calculate the heat-up of the vessel wall so that the time at which the vessel may rupture can be determined. In case of dry cavity, the molten corium falls on the concrete basemat. Similarly, if core catcher is present at the bottom, the molten corium falls on what is termed as a sacrificial material. Concrete or sacrificial material gets ablated due to intense heat content of the corium as well as decay heat. This sacrificial material performs several functions like: l

l

reduction of the heat flux by a considerable increase in the molten corium volume and heat transfer surface, oxidation of free zirconium without hydrogen release resource,

498 l

l

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

elimination of the heat flux focusing effect by lowering density of the post-interaction oxidic melt below the density of molten steel, which guarantees the inversion of metallic and oxidic melts, the inversion of metallic and oxidic melts excludes the risk of steam explosion, during water supply on the molten metal surface and restricts the hydrogen generation rate during steammetal reactions.

The major challenges in simulation of these phenomena are l

l

l

l

heat and mass transfer in the stratified molten pool with the diffusion of major components of the mixture, dependence of density on composition of the ablated material to simulate the inversion of layers, chemical reactions between corium and sacrificial material/concrete, in case the molten core-concrete interaction, additional complexity of generation of noncondensable gases due to concrete decomposition.

This makes the CFD simulation of molten core-concrete interaction a challenging task. Bolshov and Chudanov [24] developed the model for the SOCRAT severe accident system of codes for simulation of the heat and mass transfer problem in the pools containing multicomponent mixtures. Their mathematical model was based on the Navier–Stokes equations with variable properties to account for a density change after concrete melting together with the heat transfer equation, which describe the diffusion of components between phases. They validated their code with MASCA-RCW test. Chai et al. [25] developed a new computational code using multi-physics models to simulate MCCI phenomena based on the moving particle semi-implicit (MPS) method. Momentum and energy equations were used to solve the velocity and temperature fields, and multi-physics models were developed on the basis of the basic MPS method. For validation, the CCI-2 experiment was simulated by applying the developed code. As per their simulation, with respect to sidewall ablation, good agreement was observed between the simulation and experimental results. However, axial ablation was found to be slower in the simulation, which may be due to the underestimation of the enhancement effect of heat transfer provided by the moving bubbles at the bottom. The simulation results confirm the rapid erosion phenomenon observed in the experiment, which, in the numerical simulation, is explained by solutal convection provided by the liquid concrete at the corium/concrete interface. Similar model has been used by Li et al. [26]. They carried out 3-D simulation of 2-D MCCI experiment VULCANO VB-U7 with improved MPS method. Heat conduction, phase change, and corium viscosity models have been developed and incorporated into MPS code developed by them. The influence of thermally stable silica aggregates was investigated by setting up different simulation cases for analysis. The simulation results suggested reasonable models and assumptions to be considered in order to achieve best estimation of MCCI with prototypic oxidic corium and siliceous concrete. Once the desired ablation has taken place in the core catcher, water is flooded into the core catcher to remove the decay heat according to accident management strategy of some LWRs. Top flooding is the simplest strategy to tackle the core melt situation

Modelling of core melt scenarios in nuclear reactors

499

where ample amount of melt is flooded on the top with water. The coolability of molten pool under top flooding strongly depends on the water ingression inside the crust. In literature, few experiments have been reported on melt coolability under top flooding condition. From these experiments, it was observed that in some cases water ingression occurred, in some cases it did not occur, and in some cases it stopped midway. However, the reason behind this was unclear. Besides, it is observed that there are very few efforts on modelling of water ingression phenomena in melt pool when flooded with water. Originally, the motivation behind modelling of water ingression phenomenon did not aim to study the quenching of molten corium pool during a postulated severe accident condition in an LWR. In fact, models were developed for simulation of cracking behaviour of hot rocks in geological reservoirs. In this context, Lister [27] did pioneering work in modelling the penetration of water into hot rocks by considering the simplest possible one-dimensional model based on the concept of crack front propagation. Epstein [28] used Lister’s models of bulk permeability of cracked rock and developed a model for water penetration into initially molten, heat generating rock like material at low pressure, which resembles the water ingression phenomena into molten corium pool. However, the applicability of the model on melt coolability in reactor condition has not been established previously. A detailed modelling of melt coolability under top flooding condition has been presented in Section 6.5. A few experiments on top flooding confirm the fact that top flooding may not be a reliable strategy to achieve complete coolability in case of core melt accidents. Hence, efforts were made to study the melt coolability behaviour under an alternative strategy, i.e. bottom flooding. While a few small-scale experiments showed evidence of better coolability under bottom flooding condition, there is lack of understanding of the phenomenology of melt coolability under bottom flooding. There are no mechanistic models in literature for predicting coolability under bottom flooding. To some extent, Paladino et al. [29] and Widmann et al. [30] attempted modelling of melt coolability under bottom flooding, but their focus was to predict the porosity of the melts formed during bottom injection and its effect on coolability. No attempt is made yet to model entire phenomenology of bottom flooding, which is very essential in order to study coolability of molten core in ex-vessel core catchers of nuclear reactors. Similar to top flooding, phenomenon of melt coolability under bottom flooding was postulated and a model has been developed considering the heat transfer from water at bottom to the melt, steam formation and pressurization at the bottom, steam eruption through melt, porous zone formation, and subsequent coolability of the melt. The model is presented in Section 6.5

6.3

Simulation of in-vessel retention of corium in PHWRs using CFD

The pressurized heavy-water reactors (PHWRs) have several inherent and engineered safety systems to restore the reactor to safe states in the event of an accident. Most of the safety systems have diverse and redundant features. There are two independent

500

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

fast-acting shutdown systems of diverse nature. The purpose is to enhance the availability of these safety systems so that failure of one of them does not lead to progression of accidents. In spite of all these, as said before, a low probability accident can be cast beyond the acceptable design basis envelope, which can progress to severe accident by failure of multiple safety systems. In PHWRs, severe accident is termed as ‘Severe core damage accident (SCDA)’, which generally starts with limited core damage accident (LCDA) with failure of some channels where cooling is insufficient [31]. The residual coolant inside the channel boils off and fuel temperature starts rising and this may result in fuel failure. These failed fuel channels discharge large amount of enthalpy in the form of steam and hot fuel nodules into the moderator and pressure starts building inside calandria. The initial pressure spike can damage the neighbouring channels by collapsing the calandria tube on to the pressure tube. In addition to this, pressure tube may balloon at high pressure and temperature conditions, and may touch the calandria tube. Once the annular gap between calandria tube and pressure tube becomes negligible, thermal resistance to the heat transfer from hot channel to relatively cold moderator decreases and rate of deposition of decay heat in the moderator significantly increases. Large amount of heavy-water moderator present in calandria vessel provides a huge heat sink. This can significantly delay accident progression and can even terminate further escalation if moderator cooling can be restored. However, if moderator cooling and make up is not available, temperature of moderator starts rising and eventually boiling takes place (Fig. 6.1A). This will pressurize the calandria vessel leading to opening of the Over Pressure Rupture Device (OPRD). This causes moderator to gradually boil off and channels start uncovering. These uncovered channels start heating up leading to rapid temperature excursion and fuel melting. Since the fuel channels are located horizontally in the calandria, their failure mode will be different from those for the vertical fuel assemblies in PWRs. The long fuel channels will sag when they are heated and dislocate to the lower rows (Fig. 6.1B). Gradually more and more channels get uncovered due to boiling of moderator. Unable to bear the load of disassembled channels, covered channels also collapse. This collapsed core will form a debris bed at the bottom of the calandria vessel, which is generating decay heat and is being cooled by moderator (Fig. 6.1C). Boiling of moderator continues and after certain period, entire moderator is lost from calandria vessel. The hot core debris consisting of broken pressure tubes and calandria tubes, fragmented fuel, control rods and structural material starts heating up in the absence of moderator cooling and forms a mushy zone. This semi-solid corium mass gets heated up by decay heat and is cooled by relatively cooled calandria vessel wall through indirect cooling by calandria vault water. This phenomenon of simultaneous heating and cooling results in formation of molten pool surrounded by a thick crust of the corium inside the calandria vessel (Fig. 6.2) Thus, it is interesting to evaluate whether calandria vessel (CV) of PHWR can act as inherent ‘core catcher’ having large amount of coolant available in the vault around CV. While we are evaluating the capability of calandria to act as inherent core catcher, it would be essential to know the modes, which can lead to calandria vessel failure.

Modelling of core melt scenarios in nuclear reactors

501

Fig. 6.1 Progression of severe accident in PHWRs. (A) Fuel heatup. (B) Channel uncovery and (C) disassembly.

6.3.1 Necessity of in-calandria retention in PHWR In PHWR type reactors, the in-calandria retention of corium is probably the best option for managing a core melt accident. Calculations have shown that if corium breaches the calandria vessel and enters the calandria vault, due to MCCI, large

502

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.2 Terminal debris configuration.

amount of hydrogen (>2000 kg) and other fission gas will generate [32]. Due to this, there is a possibility of containment failure. In addition, managing such large amount of hydrogen by passive autocatalytic recombiners (PARs) is a concern. Hence, it is essential to contain the corium inside the calandria vessel and cool it from outside by calandria vault water.

6.3.2 Past studies in the field of severe accident The design parameters of PHWR, PWR, and BWR plants are quite different. Some striking dissimilarities between design features of PHWR plant from other light-water reactors are: the plant adopts a dual primary heat transport system having heavy water in horizontal pressure tubes under forced circulation in a figure-of-eight loop and has an additional amount of cooling water in the calandria vessel and calandria vault. Another feature is that the Calandria vessel is always submerged in water during normal operation. The design pressure of containment is lower than that of the typical PWR or BWR. In PHWRs, severe accidents leading to core damage may occur due to failure of multiple safety systems, for example a prolonged ‘station blackout’ (SBO) scenario with multiple failures of cooling systems [33]. Numerical analysis was performed to evaluate the thermal loads and stress in calandria vessel in severe accident scenario

Modelling of core melt scenarios in nuclear reactors

503

by using a series of codes [33]. Thermal loads were evaluated from the temperature distribution in the vessel wall. The Chaboche model was used for calculation of plastic strain. A maximum strain of 9.5% was found which is less considering the critical strain value at failure. Hence, it was concluded that the calandria does not undergo creep failure because of temperature and structural loading as long as the calandria vessel wall is cooled by reactor vault water. Detailed analysis for CANDU severe accident was recently published by Dupleac [34]. He studied molten pool and crust formation during a CANDU late-phase severe core damage accident employing the ANSYS-FLUENT code. The analysis used the high-resolution computational fluid dynamics–large-eddy simulation (CFD–LES) methods for CANDU configuration and materials aiming to determine the potential local effects on flow and heat transfer mechanisms during terminal debris melting phenomena. LES was used for modelling of the basic turbulent flow features inside the melt pool with the same resolution level as in Reynolds-averaged Navier–Stokes (RANS) modelling. The distinguishing two scales in flow fluctuations; i.e. (1) large-scale fluctuations that are calculated directly and (2) small, grid-scale fluctuations, which are modelled indirectly. A slice from the 3D calandria geometry was used, in order not to exceed the computing resources. The thickness and position on calculation domain of the crust were calculated in the twophase (solid and liquid) simulations. As per the simulation, due to low conductivity of the crust (porous solid material), the temperature dropped from the molten pool temperature to below 600 K on the inner surface of the calandria wall. Also, one important observation was that the variation of crust thickness in time has different evolution for the upper part than that for the lower one, which means that more energy is transferred through the free upper surface by radiation and convection to the steam. In another study conducted by Prasad and Nayak [35], experiments were carried out in a scaled test facility. The main objective of this exercise was to determine the extent of heat transfer from the melt pool to the outside vault water through calandria vessel. In the very first experiment of this kind, a scaled calandria vessel was immersed in a scaled vault containing water at ambient temperature. Molten glass at 1200°C was used as a simulating material for corium, which was poured in the calandria vessel submerged in the calandria vault water. The temperatures profiles were recorded until ambient temperatures were attained in the system. High-temperature gradient was observed within molten pool along the radial direction. Within about 100 s from the start of the experiment, crust of about 20-mm thickness was formed and it took almost 2.5 h to completely solidify the melt. Maximum average inner surface temperature of vessel was found to be 249°C and maximum average outer surface temperature of vessel was 94°C. The steam generation in the vault water was negligible. When the vessel was opened, no gap between the crust and vessel was observed. A thick crust had formed on the inner side of the vessel during the initial phase of experiment, which considerably insulated the vessel and water in the calandria vault did not boil. The heat flux at the outer wall of the vessel was found to be less than critical heat flux. This exercise was repeated with in situ heating of the melt by electrical heaters to quantify the effects of decay heat on the heat transfer behaviour and the thermal load on the vessel wall and calandria vault. Decay heat of 0.7 MW/m3 in the melt pool was

504

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

simulated using high-watt cartridge-type heaters. The temperature distributions inside the molten pool, across the vessel wall thickness, and vault water were measured. After 3 h, melt temperature attained steady state at 800°C; however after 4 h, the temperature of melt started rising. With heat generation in melt, it was observed that initially crust thickness grew up slowly to 40mm and then it stabilized at that thickness. Maximum average inner surface temperature of vessel was found to be 265°C and maximum average outer surface temperature of vessel was 98°C. Maximum heat flux was still less than CHF and the steam generation in the calandria vault was negligible. When the vessel was opened, no gap between the crust and vessel was observed. Thermal stratification in water tank was witnessed and water below 200 mm of tank level almost remained at ambient temperature and did not take part in heat removal process. No steam generation was measured in the steam flow metre. This finding was also confirmed from the level of the water tank that remained practically constant. Experimental results obtained from this study were compared with previous results (without decay heat).

6.3.3 CFD simulation of in-vessel retention CFD simulation of severe accident in calandria vessel in PHWRs is important to have better understanding of the progression of the accident, detailed thermal-hydraulic mapping of core, thermal loads on calandria vessel, and development and improvement of accident management strategies. To achieve the goal, an approach to effective use of CFD has been made to guide and support the development of models suitable for accident analysis. The following section of the numerical investigation involves simulation of tests carried out for melt coolability in simulated PHWR calandria vessel by Prasad and Nayak [35].

6.3.3.1 Geometric modelling and meshing The cylindrical vessel simulating calandria vessel in experiment has 300-mm diameter and 460-mm length, which is made of SS304L. The wall thickness of cylindrical test section is 25 mm. Two-dimensional geometry of same dimension has been used for numerical simulation. The calandria vault water has not been considered as a part of the computational domain; rather its effect is taken by applying appropriate boundary condition on calandria vessel outside surface. Scaled test facility and corresponding model of computational domain is shown in Figs 6.3 and 6.4. The grid resolution required for independent results for molten pools is a function of pool Rayleigh number. In the case of PHWR, the Rayleigh criterion (109–1012) is smaller than for LWR systems (due to lower volumetric heat generation rate and dispersion of the debris in the large calandria vessel); thus, it is possible to have larger grid cells. In this analysis, unstructured tri-grid method of meshing is adopted. Each grid point (or node) is uniquely identified by indices i, j, and the corresponding Cartesian coordinates as xi, j, yi, j respectively. The numerical solution provided by the model tends to be a unique solution as we increase the mesh density. But, a fine mesh also leads to an increase in computational time and resources required. The mesh is said to be

Modelling of core melt scenarios in nuclear reactors

505

Inlet for melt Water tank Water Cylindrical vessel

(A)

500.0 200.0

350.0

900.0 650.0

Molten pool

(B)

Fig. 6.3 (A) Actual test facility; (B) Schematic of the facility.

Fig. 6.4 (A) Schematic of computational domain for case without decay heat; (B) Tri Mesh for case without decay heat case; (C) Schematic of computational domain for decay heat case; (D) Tri mesh for case with decay heat case.

506

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

refined when there is no significant change in the values of solution results with further refinements of mesh. Hence, a grid independence test has been performed on the present 2-D model to determine the accuracy of the results and to minimize errors due to meshing. Four different grid systems having element counts of 2946, 3937, 5802, and 22099 were used for grid test. The mesh having 5802 elements is selected for the current numerical investigation, as the relative deviation of temperature is almost constant beyond this cell count.

6.3.3.2 Physical Model The computational domain contains three regions; vessel wall, molten glass, and air. The study with decay heat employs heat source inside the molten pool. Pressure-based transient solver is used for computation. Pressure-based solvers are generally used for low-speed incompressible flows. In this method, flow field is obtained by solving x, y, z momentum equations and pressure field is obtained by solving a pressure correction equation. Time-varying temperature field is obtained by solving transient energy equation. For transient simulations, the governing equations are discretized in both space and time. The spatial discretization for the time-dependent equations is identical to the steady-state case. Temporal discretization involves the integration of every term in the differential equations over a time step. Gravitation effect with corresponding gravitational acceleration of 9.81 m/s2 is considered. Effect of solidification and melting, energy, viscous (laminar) and radiation is included for modelling. The radiation model is used to account for the radiation exchange in an enclosure of gray-diffuse surfaces. The energy exchange between two surfaces depends on their size, separation distance, and orientation. The present study is based on the following assumptions/simplifications: a) the volume change upon phase change (solidification/melting) is ignored, b) molten glass in the liquid phase is considered to be an incompressible and Newtonian fluid, c) both the solid as well as the liquid phase is homogeneous and isotropic, and the solidification process is symmetric within a segment, d) the solid is homogeneously distributed in the mushy region (homogeneous mixture), e) the complex stratified temperature fields of the pools are not used in the present analysis, f ) radiation heat transfer between inner surface of the calandria vessel and the top surface of the molten pool is considered, and g) conduction through liquid glass is considered as the only mode of heat removal from the debris to the vessel.

6.3.3.3 Material Properties The calandria vessel wall is made of SS 304L. Sodium borosilicate glass forms the molten pool and remaining volume of vessel is filled by air. All the thermo-physical properties of participating media are considered to be temperature dependent. This temperature dependency has been provided in the form of polynomial function or table. Brief description of properties is given here.

Modelling of core melt scenarios in nuclear reactors

507

1. Calandria vessel; SS304L (at steam saturation temperature) Density: 7800 kg/m3 Specific Heat: 500 J/kg K Thermal conductivity: 16 W/m K 2. Molten Sodium Borosilicate glass (at 1200°C) Density: 2200 kg/m3 Viscosity: 0.0025 kg/m sec Thermal Conductivity: 1 W/m K Specific Heat: 800 J/kg K Solidus Temperature: 600°C Liquidus Temperature: 850°C Latent heat of fusion: 100 kJ/kg l

l

l

l

l

l

l

l

l

l

6.3.3.4 Boundary and initial conditions The boundary conditions used for analysis are shown in Fig. 6.5. The inner and outer surfaces of cylindrical calandria vessel are stationary wall with no slip condition and coupled thermal condition. Molten pool and air interface is taken as wall with no slip condition. The upper surface of the molten pool/debris bed exchange heat with the vessel inner wall by thermal radiation and with the air above the surface by natural convection and conduction as shown in Fig. 6.5. The emissivity of the upper surface of the debris bed is taken as 0.7 (for T > 1073 K) and the emissivity of the calandria wall is 0.3 [36]. Heat transfer from the molten pool to the water in calandria vault takes place by conduction through calandria vessel wall. For the case with decay heat, conditions prevailing in the experiment are replicated. Heated surface temperature is maintained at 900°C corresponding to volumetric heat generation of 0.7 MW/m3. This is equivalent to deposition of 9.2 kW power in the melt pool. The calandria vessel is submerged in water. As mentioned earlier, the vault water is not part of computational domain and hence its effect in the heat transfer has been accounted by applying appropriate convective heat transfer coefficient. At the outer surface of calandria vessel, heat transfer is due to natural convection and by partial nucleate boiling. Heat transfer coefficient on the outer surface of calandria vessel wall is supplied as a function of vault water temperature as an input to the solver. This is implemented to provide the local heat transfer values with change in temperature with time. The Nusselt criterion for natural convection at outer surface of calandria is calculated according to Churchill and Chu correlation [37], i.e. 8 > > > > > <

91=6 > > > > > =

Gr Pr 1=2 NuAvg ¼ 0:6 + 0:387 "
 #16=9 > > > > 0:559 9=16 > > > > > > ; : 1+ Pr

(6.1)

For fully developed nucleate boiling, the heat transfer coefficient can be reproduced well by the simple empirical equation according to Fritz [38]. The analysis starts after

508

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.5 (A) Pictorial representation; (B) Boundary condition without decay heat; (C) Boundary condition with decay heat.

Modelling of core melt scenarios in nuclear reactors

509

the simulant melt has been poured into the scaled calandria vessel. At t ¼ 0 s, the temperature field (without decay heat case) for all domains is known: corium simulant/ debris temperature is 1100°C, and calandria wall temperature is 32°C. The calandria vault water temperature at the beginning is 32°C. For simulating 30-s duration, it takes around 2 min clock time on Intel i5 3.6 GHz computer with 64-GB RAM. Same initial conditions are used for analysis with decay heat except that the initial heater surface temperature is maintained at 800°C.

6.3.3.5 Modelling and solution algorithm The CFD calculations are carried out on the computational domain by solving the conservation equations of mass, momentum, and energy for transient conduction and natural convection with appropriate initial and boundary conditions using finite volume method. For the pressure-velocity coupling, the SIMPLE scheme has been used. Second-order upwind scheme is used for the discretization of convection and diffusion terms in the conservation equations. Second-order implicit formulation is used for the transient calculations. Convergence is assumed when the residual of u, v velocity is less than 103, and for energy, it is less than 106. This numerical simulation also involves solidification and/or melting taking place over a range of temperatures (e.g. in binary alloys). This has been modelled with enthalpy-porosity model for tracking of liquid-solid front. The liquid-solid mushy zone is treated as a porous zone with porosity equal to the liquid fraction, and appropriate momentum sink terms are added to the momentum equations to account for the pressure drop caused by the presence of solid material. For solidification/melting problems [39], the energy equation is written as ∂ ðρH Þ + ΔðρvH Þ ¼ ΔðkΔT Þ + S ∂t

(6.2)

where, H ¼ enthalpy, ρ ¼ density, ν ¼ fluid velocity, and S ¼ source term.

The enthalpy of the material is computed as the sum of the stored heat, h, and the latent heat, ΔH: H ¼ h + ΔH

(6.3)

where, ðT h ¼ href +

CpdT Tref

(6.4)

510

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

where, href ¼ reference enthalpy, Tref ¼ reference temperature, and Cp ¼ specific heat at constant pressure.

The liquid fraction, β, can be defined as β¼0

if T < solidus temperature

β ¼1

if T > liquidus temperature

(6.5)

The latent heat content can now be written in terms of the latent heat of the material, L: ΔH ¼ βL

(6.6)

The latent heat content can vary between zero (for a solid) and L (for a liquid). The momentum sink due to the reduced porosity in the mushy zone can be written as a function in the following form: s ¼ f ðβ, v, mushy zone constantÞ

(6.7)

where β is the liquid volume fraction, ν is the solidified material pull velocity out of the domain (here v ¼ 0).

6.3.3.5.1 Analysis without decay heat The simulation starts with melt pool being at 1100°C and vessel being at 32°C. Heat is convected away from the outer surface of the vessel by vault water at 32°C by natural convection and nucleate boiling. The temperature is monitored at several locations throughout the computational domain, which includes inner and outer surfaces of the vessel at different orientations and molten glass pool. These points of observation are shown in Fig. 6.6. Fig. 6.7 shows the scheme of nomenclature for different angles and heights in the vessel. Fig. 6.8A and B shows the predicted variation of temperature on the inner and outer surface of the simulated calandria vessel at different angles as shown in Fig. 6.7. It is clear from the figures that the maximum temperature on the vessel inner wall as well as outer wall occurs at the bottom most (at 0° location). The maximum temperature on inner surface of vessel instantly rises to nearly 180°C and then gradually falls to about 100°C and then stabilizes at that temperature. The curve for temperature at bottom most point (0° location) and at 45° location on inner surface is almost coinciding. Similar trend is obtained for vessel outer temperature at bottom and at 45° location. From the outer wall of the calandria vessel, heat is transferred mainly by natural convection. At bottom most point, velocity being lowest, the heat transfer is on the lower side. On the other hand, at 90° location, the buoyancy is high which enhances the heat transfer. Hence, temperature at bottom most point (0° location) is higher than

Modelling of core melt scenarios in nuclear reactors

511

Fig. 6.6 Radial location on test section.

120 mm 100 mm 80 mm 60 mm 40 mm Inner surface Outer surface

0° Ref. line 180°

Fig. 6.7 Circumferential location on test section.

90°

45°

0°

that of 90° location. The temperature of the topmost point (180°) is lowest on both inner and outer wall as the melt pool does not touch the vessel at that point. The rise in temperature is due to conduction and radiation heat transfer. Fig. 6.9 shows comparison of simulation results with experimental findings of scaled facility of Prasad and Nayak [35]. It is evident that the physical modelling

512

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.8 Temperature distribution at various locations: (A) vessel inner surface; (B) vessel outer surface.

Fig. 6.9 (A) Comparison of temperature variation on the inner and (B) outer surface at 0° location.

and simulation technique used for the solution captures the actual physics reasonably well. The minor differences between the two findings are mainly due to the uncertainties in the experimental procedure and approximation in imposition of boundary condition at the vessel outer surface as calandria vault water is not taken in to computational domain. Simulation results predict the vessel inner surface quite well but over-predict the vessel outer surface by around 12%. However, both predict the temperature to stabilize near 100°C on the inner surface of vessel and close to 80°C on the outer surface. Fig. 6.10 shows variation of temperature at different locations in the molten pool along the vertical direction from the vessel inner surface. The temperature shows a gradual increase from near wall points, which clearly indicate the efficient cooling on the outer surface. However, molten glass pool closer to interface record higher temperature as heat transfer by radiation and natural convection by air is not sufficient to remove the heat effectively from free surface of the pool. The crust formed at the interface further hampers the rate of heat transfer out of the molten pool. The formation of crust on the inner surface of the vessel and marching of solid-liquid interface is shown in Fig. 6.10B. It is clearly visible from the figure that the crust starts forming immediately after coming into contact with relatively cooler vessel surface.

Modelling of core melt scenarios in nuclear reactors

513

Fig. 6.10 Predicted and measured. (A) Temperatures in the melt pool at location 0° at different radial locations; (B) Crust thickness variation with time.

6.3.3.5.2 Numerical analysis of scaled experimental facility with decay heat As mentioned earlier, the foregoing experiment [35] relating to coolability of scaled calandria vessel under severe accident was repeated to quantify the effects of decay heat. The decay heat was simulated by cartridge-type high watt heaters whose surface was maintained at 900°C. The decay heat was equivalent to 0.7 MW/m3. The same was simulated in numerical calculation and constant surface boundary condition was applied to the heater surfaces. Fig. 6.11 shows the numerical results of temperature variation of inner and outer vessel surface at various circumferential locations. At inner and outer surface, maximum temperature of vessel surface was observed at 0°. It can be seen that, at initial time step, temperature of vessel inner and outer surface rises suddenly to a peak temperature and then decreases. This was due to the formation of crust inside melt, which offers resistance to heat transfer. As visible from these results, the temperature of the calandria inner wall surface does not saturate, rather behaves like quasi steady state after 50 min. However, the temperature does not reach alarming values. Similar trend is obtained for the outer vessel surface, where the temperature at 45° position reaches the boiling point. This is due to the effect of decay heat generation inside the melt.

514

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.11 (A) Temperature distribution on inner and (B) outer surface of calandria vessel at various locations under decay heat.

Fig. 6.12A shows predicted and experimental temperature distribution in molten pool as function of time at different radial locations away from CV inner wall. The corium simulant (glass) is poured into scaled calandria vessel at nearly 1100°C. All the thermocouples immediately read similar temperatures. The rate of temperature

Modelling of core melt scenarios in nuclear reactors

515

Fig. 6.12 Predicted and measured (A) temperature distribution in the molten pool inside test section at 0° location; (b) Crust thickness.

drop is maximum near the CV wall and it decreases towards the free surface of the simulant pool. Temperature near the wall; at 20 mm and 40 mm sharply drops to 500°C and 600°C and then settles at around 500°C and 400°C respectively. This steep temperature gradient at near wall locations, just after pouring indicates rapid crust formation even in the presence of decay heat. The crust thus formed creates significant resistance to further heat transfer. This leads to huge temperature gradient in radial direction near the wall. Temperature of melt near centre and free surface experience small temperature drop and stabilizes at slightly elevated temperatures. Fig. 6.12A also validates the numerical results against experimental findings. This gives enough confidence over numerical models to simulate prototype conditions. Fig. 6.12B shows the crust formation with time from the bottom of the cylindrical vessel. Initially the crust thickness grew and after that it almost remained constant. The thickness of crust does not increase linearly as in the case without decay heat. Here the heat generated by heaters immersed in the pool prevents formation of crust away from the wall and solidification front does not march monotonically. Moreover, the crust near the 45° location does not grow more than 20 mm due to presence of heaters in the vicinity, which maintains a constant temperature of 900°C unlike that in experiment.

516

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

6.3.3.5.3 Comparison of results of analysis with and without decay heat Fig. 6.13A shows the comparison of the test section inner surface temperature variation with and without decay heat. It is observed that in case of decay heat, higher temperature of the inner surface of the test section occurs initially. This was due to presence of both stored and decay heat inside the molten pool. After 110 min, the effect of decay heat has been observed in the temperature distribution profile. In case of no decay heat, temperature reduce gradually while in presence of decay heat, it does not saturate and increases gradually after 110 min. It is evident from Fig. 6.14 that in case no decay heat, higher temperature of the outer test section with slight nucleate boiling occurs initially. The interface between crust and melt (mushy region) offers resistance to heat transfer in case of decay heat. In this figure, both the inner and outer Calandria vessel surface temperatures are found to reach a quasi-steady-state temperature after 25 min. Fig. 6.15 shows the chronological contours of temperature in the pool from the beginning of simulation for scaled model. These pictures provide an overall progression of events taking place during cooling of molten glass with and without presence of decay heat. The isotherms inside the lower pool are studied at different times for understanding the way by which the high-temperature fluid is directed for passive heat removal. In Fig. 6.15, the molten pool temperature distributions are presented at 2-min, 30-min, 60-min, 120-min, and 180-min moments of time for both; no decay and with decay consideration. Temperature distribution in the molten pool and the crust thickness variation with time are output data of the models. Fig. 6.15A–F shows the formation of crust near the inner wall surface of test section. It can be observed that crust is formed within a short time. Temperature of glass monotonically decreases when decay heat is absent (Fig. 6.15A–E) which causes solidification front to march towards centre. The crust thickness goes on increasing in this case as seen in Fig. 6.15B–E. However when decay heat is supplied through simulated immersed heaters, the heat loss from glass is compensated by heaters. This does not allow

Fig. 6.13 Comparison of temperature distribution on vessel inner surface at location 0°.

Modelling of core melt scenarios in nuclear reactors

517

Fig. 6.14 Crust thickness (A) comparison between CFD simulation of two cases (with and without decay heat); (B) comparison of crust thickness with decay heat between simulation and experimental finding.

the solidification front to move beyond approximately 40 mm. The contour plots after 180 minutes clearly show that without decay heat, temperature near melt centre gradually decreases to nearly 700 K while with decay heat present, temperature remains above 900 K.

6.3.4 CFD simulation of melt pool coolability in calandria vessel in prototype plant In previous sections, we have discussed results obtained from CFD simulation of scaled experimental facility for severe accident with and without decay heat. These results were compared with corresponding experimental results. The analysis of

518

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Temperature 1372.95 1296.87 1220.78 1144.70 1068.62 992.54 916.46 840.37 764.29 688.21 612.13 536.04 459.96 383.88 307.80 [K] At 2 min

Temperature 1372.27 1296.17 1220.07 1143.97 1067.88 991.78 915.68 839.58 763.49 687.39 611.29 535.20 459.10 383.00 306.90 [K] At 2 min

(A)

(F) Temperature 1372.27 1296.17 1220.07 1143.97 1067.88 991.78 915.68 839.58 763.49 687.39 611.29 535.20 459.10 383.00 306.90

Temperature 1372.95 1296.87 1220.78 1144.70 1068.62 992.54 916.46 840.37 764.29 688.21 612.13 536.04 459.96 383.88 307.80

[K]

[K]

At 30 min

(B)

(G)

Temperature 1372.27 1296.17 1220.07 1143.97 1067.88 991.78 915.68 839.58 763.49 687.39 611.29 535.20 459.10 383.00 306.90 [K]

Temperature 1372.95 1296.87 1220.78 1144.70 1068.62 992.54 916.46 840.37 764.29 688.21 612.13 536.04 459.96 383.88 307.80 [K]

At 60 min

(C)

At 60 min

(H) Temperature 1372.95

Temperature 1372.27

1296.87

1296.17

1220.78

1220.07

1144.70

1143.97

1068.62

1067.88

992.54

991.78

916.46

915.68

840.37

839.58

764.29

763.49

688.21

687.39

612.13

611.29

536.04

535.20

459.96

459.10

383.88

383.00 [K]

At 30 min

[K]

306.90

At 120 min

(D)

307.80

At 120 min

(I)

Temperature 1372.27

Temperature 1372.95

1296.17

1296.87

1220.07

1220.78

1143.97

1144.70

1067.88

1068.62

991.78

992.54

915.68

916.46

839.58

840.37

763.49

764.29

687.39

688.21

611.29

612.13

535.20

536.04

459.10

459.96 383.88

383.00

307.80

306.90 [K]

(E)

[K]

At 180 min

At 180 min

(J)

Fig. 6.15 Temperature distribution contours for no decay heat (left) (A)–(E) and with decay heat (right) (F)–(J) at t ¼ 2 min, 30 min, 60 min, 120 min, and 180 min.

Modelling of core melt scenarios in nuclear reactors

519

results obtained from simulation of scaled experimental facility clearly indicates that the formulation used for the CFD solution closely captures the physical phenomena and is able to simulate the experimental behaviour reasonably well. The following section deals with application of CFD models to prototypic condition. For prototypic condition, actual calandria vessel was modelled in 2-D. The Indian PHWR 700 MWe calandria vessel has diameter close to 8 meters and length is 6 meters. It is made of SS 304L having 32 mm wall thickness. The simulation starts when the moderator has already boiled off and the debris bed is above solidus temperature ( 2950 K) and being indirectly cooled by calandria vault water, which is at 60°C. The phenomena considered are natural convection in melt pool, solidification and melting in molten debris, conduction in vessel, convection on outer side of vessel. The thermo-physical properties of all participating media are considered to be temperature dependent. This temperature dependence can be provided either in the form of polynomial function or table. Brief description of material properties is given as follows.

6.3.4.1 Material Properties Calandria vessel; SS304L (at steam saturation temperature) l

l

l

Density: 7800 kg/m3 Specific heat: 500 J/kg K Thermal conductivity: 16 W/m K

Corium (properties at 2950 K) l

l

l

l

l

l

Density: 9800 kg/m3 Thermal conductivity: 3 W/m K Specific heat: 630 J/kg K Solidus temperature: 2900 K Liquidus temperature: 3100 K Latent heat of fusion: 260  103 J/kg

6.3.4.2 Boundary and initial condition Boundary and initial conditions for simulation of the prototype are supplied in similar fashion as that of scaled facility. The inner and outer surfaces of cylindrical calandria vessel are stationary wall with no slip condition and conducting vessel wall. Molten pool at air/steam interface is taken as wall with no-slip condition. The upper surface of the molten pool/debris bed exchanges heat with air/steam above the free surface by natural convection and conduction. It exchanges heat with vessel inner wall by thermal radiation. The emissivity of the upper surface of the debris bed εdb is taken as 0.67 and the emissivity of the calandria wall, εcs is 0.3. Heat transfers from the molten pool to the water in calandria vault take place by conduction through calandria vessel wall. The calandria vault water is not part of computational domain and therefore appropriate convective heat transfer coefficient is supplied by a user-defined function based on

520

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.16 Nodalization of the prototypic geometry.

temperature of calandria vault water. The heat is thus conducted from corium crust to the calandria vessel wall and then gets convected by vault water. For decay heat analysis, volumetric heat generation rate of 1.0 MW/m3 is supplied to the melt pool. Fig. 6.16 shows tetra mesh with minimum quality 0.2 used for simulation.

6.3.4.3 Corium coolability behaviour inside the calandria As mentioned earlier, simulation starts when the debris are at liquidous temperature ( 2950 K). Fig. 6.17A–D shows the temperature contours inside corium chronologically from beginning to the end of 3 h. The debris is initially at 2950 K inside the calandria, which is being cooled by calandria vault water. As simulation proceeds, the temperature of debris increases due to decay heat and melting starts forming a molten pool. It is observed that convection is taking place inside the melt. However, debris near the calandria wall does not melt and remains at relatively lower temperature. Thus, the molten corium is completely surrounded by the solid corium (crust), which is at low temperature. This crust effectively insulates the hot corium to directly interact with vessel wall. The thermal conductivity of the corium is very low and this causes the overall heat transfer coefficient to decrease significantly. Therefore, although the corium pool temperature near the centre is very high (close to 3500 K) after roughly 2 hours, still bulk of the water in the calandria vault does not start boiling. Fig. 6.18A shows the vessel inner and outer temperature variation with time. Partial nucleate boiling on the outer surface of calandria starts around 1.5 hours. This is marked by fluctuation of calandria vessel outer surface temperature as shown in Fig. 6.18A. Nearly after 2.5 h, water inside the calandria vault attains the saturation

3.507e+003 3.280e+003 3.053e+003 2.827e+003 2.600e+003 2.373e+003 2.146e+003 1.920e+003 1.693e+003 1.466e+003 1.240e+003 1.013e+003 7.861e+002 5.594e+002 3.327e+002

(A)

1800 s

(B)

Temperature vessel wall

Temperature prototype 3.507e+003 3.280e+003 3.053e+003 2.827e+003 2.600e+003 2.373e+003 2.146e+003 1.920e+003 1.693e+003 1.466e+003 1.240e+003 1.013e+003 7.861e+002 5.594e+002 3.327e+002

3.338e+003 3.137e+003 2.935e+003 2.734e+003 2.532e+003 2.331e+003 2.129e+003 1.928e+003 1.727e+003 1.525e+003 1.324e+003 1.122e+003 9.208e+002 7.194e+002 5.179e+002 3.165e+002 [K]

(C)

3600 s

Temperature prototype 3.507e+003 3.280e+003 3.053e+003 2.827e+003 2.600e+003 2.373e+003 2.146e+003 1.920e+003 1.693e+003 1.466e+003 1.240e+003 1.013e+003 7.861e+002 5.594e+002 3.327e+002

7200 s

Temperature Corium 3.338e+003 3.137e+003 2.935e+003 2.734e+003 2.532e+003 2.331e+003 2.129e+003 1.928e+003 1.727e+003 1.525e+003 1.324e+003 1.122e+003 9.208e+002 7.194e+002 5.179e+002 3.165e+002

[K]

10800 s

(F)

7200 s (enlarged view)

Fig. 6.17 Temperature contours in the debris bed (A) 1800 s, (B) 3600 s; (C) 5400 s; (D) 7200 s; (E) 10800 s, and (F) 7200 enlarged view.

521

(E)

(D)

5400 s

Modelling of core melt scenarios in nuclear reactors

Temperature prototype 3.507e+003 3.280e+003 3.053e+003 2.827e+003 2.600e+003 2.373e+003 2.146e+003 1.920e+003 1.693e+003 1.466e+003 1.240e+003 1.013e+003 7.861e+002 5.594e+002 3.327e+002

Temperature prototype

522

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

1100

CV inner CV outer

1000

Temperature (K)

900 800 700 600 500 400 300

2000

4000

(A)

6000

8000 Time (s)

10,000

12,000

14,000

3500

3000

Temperature (K)

2500

2000

1500

1000

500

0 −3800

(B)

−3600

−3400

−3200

−3000

Y [mm]

Fig. 6.18 (A) Temperature profile on the calandria vessel inner and outer wall at the bottom; (B) Temperature of corium in vertical direction from bottom of calandria wall.

Modelling of core melt scenarios in nuclear reactors

523

300,000

Heat flux (W/m2)

250,000

200,000

150,000

100,000

50,000

2000

4000

6000

8000

10,000

12,000

14,000

Time (s) Fig. 6.19 Heat flux on the calandria vessel outer surface.

temperature and causes the temperature on outer side of CV to be steady. The temperature on the inner side of the CV drops sharply at the beginning of simulation owing to formation of crust. This temperature profile becomes nearly constant after 2.5 hrs. It is seen that the vessel outer temperature remains close to 400 K and inner temperature remains around 620 K. Fig. 6.18B shows temperature of calandria vessel and molten pool after 3 hrs. The profile has been plotted along the line starting from the outer surface of the calandria vessel and ending at the solidified interface between air/steam and pool at zero degree location as shown in Fig. 6.18B. It can be observed from the profile that the CV wall temperature rises from nearly 400 K on outer surface to 620 K on inner wall. The temperature profile suggests that the crust thickness at this time is roughly 50 mm. The maximum temperature of the pool is close to 3500 K. Fig. 6.19 shows the outer wall surface heat flux variation with time. It is observed that the peak wall heat flux is around 190 kW/m2. The literature values show that the CHF occurs above 280 kW/m2 [40] for cylindrical geometries. This shows that there is sufficient margin against CHF.

6.3.5 Summary CFD model was benchmarked with various experimental data, which have given lot of confidence and insights regarding in-vessel retention in PHWR. The model was extended to prototypic condition for predicting the vessel behaviour during severe accident. CFD model has brought out that, vessel temperature remains low

524

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

during severe accident condition. Formation of crust on the inner surface of calandria vessel plays an important role in retaining the corium safely inside CV. The thick crust has very poor thermal conductivity and thus effectively insulates the calandria vessel from getting exposed to very high temperature. This crust acts as protecting layer for CV wall and also decreases rate of heat transfer to the calandria vault water. This prolongs the presence of water in the vault and thus arrests the further progression of accident. The heat flux on the outer surface of CV remains well below the CHF observed in experiments. Thus, calandria vault in PHWR has the potential for cooling the high-temperature debris inside the calandria vessel during severe accident.

6.4

CFD simulation of debris bed coolability

One of the major issues in coolability of the particulate debris bed is whether the decay heat can be completely removed from the bed by the coolant supplied from the top and/or the bottom of the bed. Coolability of particulate debris is to be considered under conditions of quenching of dry hot debris as well as boil-off of a water-filled bed in presence of decay heat. The latter case is mainly related to long-term cooling and former is during the initiation of re-flooding action as an accident management strategy. Coolability of debris beds is characterized by the dryout heat flux (DHF), which is the maximum heat flux that can be removed from the debris bed through the upper surface, without occurrence of dry spots inside the bed [41]. Bed porosity and particle sizes are the key parameters determining the dryout heat flux when operating conditions like system pressure and water supply (from the top or from the bottom) remain constant. In the top flooding configuration, a locally homogeneous inflow of water from above results in the spatially averaged two-phase flow against the upward steam flow and dryout occurs when the escaping steam prevents water inflow into the bed. Investigations by Sch€afer et al. [42] with water fed from the bottom of the bed indicated a significant increase of the power, which can be removed from the bed before reaching dryout. Most previous investigations of debris bed coolability have been carried out under strongly idealized conditions, assuming the debris beds as homogeneous (porosity and particles are uniformly distributed throughout the bed and the bed height is laterally constant), thus providing practically 1D behaviour. But, in reality, debris beds have a heap-like shape and have a nonuniform particle distribution or porosity stratification.

6.4.1 Modelling of debris bed coolability Very few CFD models are reported in literature for simulation of debris bed coolability under severe accident scenarios. The MEWA-2D (Melt and WAter) code (B€ urger et al. [22]) has been developed at IKE, University of Stuttgart, for modelling of core degradation, melt progression, and coolability during the late phase of severe accidents with core melting. In the MEWA code, the debris bed is modelled by applying a quasicontinuum approach and using a two-dimensional description in either

Modelling of core melt scenarios in nuclear reactors

525

cylindrical or cartesian coordinates. Three separate phases, solid particles, liquid coolant (water), and gas (vapour) are considered as described in Rahman et al. [41]. The solid matrix is assumed to be fixed. The mass conservation equations of the fluids are as given in Eqs (6.8), (6.9) !

∂t εαρg + — ρg j g ¼ Γ

(6.8)

!

∂t εð1  αÞρl + — ρl j l ¼ Γ

(6.9)

where ε is the porosity, α is the void fraction, jg and jl are the superficial velocities of gas and liquid, respectively, and Γ is the phase change rate, representing either evaporation or condensation. Temporal and spatial derivatives of the velocities can be neglected, which leads to the following simplified momentum conservation equations ! F pg

!

— pg ¼ ρg g +

εα

(6.10)

εα

! F pl

!

— pl ¼ ρl g +

+

! Fi

εð1  αÞ



! Fi

ε ð1  α Þ

(6.11)

For the particle-fluid drag forces, the classical Ergun [43] law (originally proposed for single-phase flow through porous media) is adapted to two-phase flow by introduction of relative permeabilities and passabilities [44] as follows: ! F pg

μg ! ρg ! ! ¼ εα jg+  j g j g KK rg ηηrg

! F pl

¼ ε ð1  α Þ

!

 μl ! ρ ! ! j l + l  j l j l KK rl ηηrl

(6.12)




(6.13)

where the basic permeabilities (K) and passabilities (η) are derived from Ergun formula: K¼

ε 3 dp 2 ε 3 dp and η ¼ 150ð1  εÞ 1:75ð1  εÞ

(6.14)

where dp is the particle diameter. To consider the drag between liquid and gas, different approaches are used. In some models, it is considered implicitly in the Krg and Krl terms as mentioned in Eqs (6.12), (6.13). Table 6.2 gives the details of the models [44]. Another way to consider interfacial drag is by adding explicit drag term in the momentum equation as given here

526

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Table 6.2 Relative permeabilities and relative passabilities of different models Model

Krel,v

ηrel,v

Krel,l

ηrel,l

Reed Lipinski Hu and theophanous

α3 α3 α3

α5 α3 α6

(1  α)3 (1  α)3 (1  α)3

(1  α)5 (1  α)3 (1  α)6

— pg ¼

μ ρ   Fi jg + jg jg + ρg g + α KK rel, g ηηrel,g

(6.15)

— pl ¼

μ ρ Fi jl + jjl jjl + ρl g  ð1  αÞ KK rel, l ηηrel,l

(6.16)

Schulenberg and M€ uller developed the correlation for drag force term as 2 
j ρl K  jl g  Fi ¼ 350ð1  αÞ α ρ  ρg g α ð1  αÞ ησ l 7

(6.17)

The expressions for permeability and passability are given as  3

5

Krel, v ¼ α , Krel,l ¼ ð1  αÞ ,ηrel, l ¼ ð1  αÞ ηrel, v ¼ 3

0:1α4 if α  0:3 α6 if α > 0:3

(6.18)

Tung and Dhir model also used a different approach to model the interfacial drag force based on the flow pattern in the porous medium. Based on visual observations in airwater flow experiment, they defined flow patterns ranges for bubbly, slug and annular flow. The expressions for relative permeability and passability are defined as


1ε Krel, v ¼ 1  εα


1ε Krel, v ¼ 1  εα

4=3




1ε α , ηrel, v ¼ 1  εα

4=3

4

4=3




1ε α , ηrel, v ¼ 1  εα 3

α4 for 0  α  0:6

(6.19)

α3 for 0:6  α  1

(6.20)

4=3

and Krel, l ¼ ð1  αÞ4 , ηrel, l ¼ ð1  αÞ4

(6.21)

Modelling of core melt scenarios in nuclear reactors

527

If the particles are having smaller diameter, capillary forces also play important role. The pressure in liquid phase and gas phase is different and is given as [41] rffiffiffi ε pc ¼ pg  pl ¼ cos θ  σ l J ðs l Þ κ

(6.22)

where σl is the surface tension, θ is the contact angle, ε is the porosity, sl is the saturation (liquid fraction in the pore), and κ is the permeability. The Leverett function depends on the liquid saturation. The correlation of Hofman and Barleon gives [41]: J ðsl Þ ¼ 0:38  ðsl + 0:014Þ0:27

(6.23)

In the energy conservation equations of gas and liquid, mechanical work due to friction and pressure forces is generally neglected. Radiation heat transfer is implicitly considered in the solid energy conservation equation through an effective thermal conductivity. For the gas, energy conservation yields:   ! 

∂ εαρg eg + r  ρg j g hg ¼ r  λeff :g rTg + Qsg  Qg, sat + Γhg, sat ∂t  ! 

∂ ðεð1  αÞρl el Þ + r  ρl j l hl ¼ r  λeff :l rTl + Qsl  Ql, sat  Γhl, sat ∂t

(6.24) (6.25)

The solid-state energy equation is given as

∂ ðð1  εÞρs es Þ ¼ r  λeff :s rTs + Qdecay  Qsg  Qsl  Qs, sat ∂t

(6.26)

6.4.1.1 Pressure drop simulation in debris bed A comparative study for prediction of pressure drop in two-phase flow inside heated debris bed with different models as mentioned earlier was carried out. For this, experimental data from Debris facility at IKE were used. Experiments were carried out in DEBRIS facility as shown in Fig. 6.20 [44]. The experimental set-up consists of a pressure vessel designed for pressures up to 40 bar in which the crucible filled with particles is mounted. The pressure vessel is connected to a storage tank filled with demineralized water and a pumping system, which allows performing boiling experiments with feeding water to the crucible at the bottom (bottom-flooding) or at the top (top-flooding). The debris bed is heated via an oil-cooled 2-winding induction coil by an RF generator. The RF generator operates at a frequency of 200 kHz and has a nominal output power of 140 kW. For two-phase pressure drop and dryout experiments, a crucible made of PTFE is used. It has a total height of 870 mm and an inner diameter of 125 mm. It is equipped with 60 thermocouples (1 mm, Type N), of which 51 are located in the debris bed on 25 levels. The thermocouples measure the temperature in the voids between the

528

Pressure vessel

B

Hot water storage tank

B

Water return line TH5

A

A

(3%)

max TL05

MS03 TH4 (intern)

RFgenerator

Electric heater

TH3

max TH2

TH1

min

60TE 7P

TL01

MS07 TL02

Fig. 6.20 Debris test facility [44].

TL03

TL04

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

MS02

Modelling of core melt scenarios in nuclear reactors

P = 3 bar

0

−1000

-dp/dz − rlg (Pa/m)

-dp/dz − rlg (Pa/m)

1000

P = 1 bar

0

−2000 −3000

Reed Lipinski Hu & Theophanus Schulenberg & Mueller Tung & Dhir Experimental

−4000 −5000 −6000 0.0

(A)

529

0.1

0.2

0.3

0.4

0.5

Jg m/s

0.6

0.7

0.8

−1000 −2000 −3000 Reed Lipinski Hu & Theophanus Schulenberg & Mueller Tung & Dhir Experimental

−4000 −5000

0.9

−6000 0.0

(B)

0.1

0.2

0.3 Jg (m/s)

0.4

0.5

Fig. 6.21 Predicted pressure drops under top flooding condition at (A) 1 bar and (B) 3 bars [44].

particles, which are filled by liquid, vapour, or a mixture of both. Pressure drop was measured with DPTs located at 8 locations along the height of the bed. Fig. 6.21 shows the comparison of experimental and predicted results with different models. It was seen that the models without consideration of explicit drag term fail to predict the two-phase-flow pressure drop at lower gas velocities for 1 bar pressure (Fig. 6.21A). Better prediction is given by the Schulenberg and Mueller model, which is able to predict the pressure drop characteristics with typical lying S-shaped curve. For 3-bar pressure (Fig. 6.21B), the prediction by Schulenberg and Mueller is much better. In two-phase pressure drop in a porous medium, there are different types of pressure drops components like liquid particle, vapour particle, and vapour liquid drag. In case of top-flooding and low inflow bottom flooding, the vapor velocities are low; thus, the interfacial drag between vapour and liquid plays an important role in the total pressure drop, which is taken into account in the Schulenberg and Mueller model. Hence, it predicts the pressure drop well. However, as the liquid velocities increase, the liquid particle drag becomes more dominant. The coefficients of Reed model are obtained from air–water flow rate at high inflows. Hence, Reed model predicts the pressure drop well. Although the variations in the predicted dryout heat fluxes among all the models are small, there is a significant difference in pressure drop characteristics under bottom flooding condition (Fig. 6.22A–D). At a low inflow velocity of 0.70 mm/s, no model seems to be able to predict the pressure gradients accurately, except the Schulenberg and Mueller model, which captures the trend and predicts reasonably well in very low vapour velocities range. For the inflow velocity of 2.8 mm/s (Fig. 6.22C), the Reed model gives the best prediction. The models of Hu and Theofanous and Lipinski give reasonably good predictions for the low vapour velocities, for a vapour velocity higher than 0.4 m/s the deviations become bigger. The predictions from the models of Schulenberg and M€uller, and Tung and Dhir are much smaller than the experimental values. When the flow velocity increases to 7.1 mm/s (Fig. 6.22D), the particle liquid drag becomes more dominant, the pressure gradients increase consistently with vapour velocity. All the models predict this tendency correctly; however, the deviations are generally very high except for the model of Tung and Dhir.

530

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment 6000

2000

Reed Jl0 = 0.5 mm/s Lipinski Hu & Theofanus Schulenberg & Mueller Tung & Dhir Experimental

5000 4000 -dp/dz − rlg (Pa/m)

-dp/dz − rlg (Pa/m)

4000

0 −2000

3000

Reed Jl0 = 0.69 mm/s Lipinski Hu & Theofanus Schulenberg & Mueller Tung & Dhir Experimental

2000 1000 0 −1000 −2000

−4000

−3000 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1

(A)

Jg (m/s)

14000 -dp/dz − rlg (Pa/m)

12000 10000

70000 60000

6000 4000

0.6

0.8

1.0

1.2

50000

Reed Jl0 = 7.1 mm/s Lipinski Hu & Theofanus Schulenberg & Mueller Tung & Dhir Experimental

40000 30000 20000 10000

2000 0

(C)

0.4

Jg (m/s)

Jl0 = 2.8 mm/s

8000

0.0

0.2

80000 Reed Lipinski Hu & Theofanus Schulenberg & Mueller Tung & Dhir Experimental

-dp/dz − rlg (Pa/m)

16000

0.0

(B)

0 0.2

0.4

0.6

Jg (m/s)

0.8

1.0

0.0

(D)

0.2

0.4

0.6

0.8

1.0

Jg (m/s)

Fig. 6.22 Predicted pressure drops under bottom flooding condition at 1 bar at (A) jl0 ¼ 0.5 mm/s, (B) jl0 ¼ 0.69 mm/s, (C) jl0 ¼ 2.8 mm/s, and (D) jl0 ¼ 7.1 mm/s [44].

6.4.2 Simulation of dryout: Top flooding conditions The code MEWA was applied on the tests conducted with top flooding condition at pressure of 1 bar in DEBRIS test facility [45]. Experimentally, dryout heat flux of 604 kW/m2 was observed. Fig. 6.23 shows the experimental temperature contours in debris bed at dryout condition. Fig. 6.24 shows the MEWA-2D code predictions for these experimental conditions. The dryout behaviour was predicted quite close to experiments by MEWA-2D code using Modified Tung and Dhir model. Of course, the location of the dryout was not predicted correctly. This can be possible because of nonuniformity in the particle sizes of the bed whereas the code uses uniform size.

6.4.3 Simulation of dryout heat flux: Bottom flooding Experiments were carried out for bottom flooding for dryout heat flux in DEBRIS facility. For studying the dryout heat flux under bottom flooding condition, a 10-mm id downcomer was placed at the centre of a polydispersed bed (Fig. 6.25) composed of spheres of 2-, 3-, and 6-mm diameter [45]. It was observed that because of the additional downflow generated by the downcomer, not only the dryout heat flux increases, but there was also a phenomenon

Modelling of core melt scenarios in nuclear reactors

531

T (°C) 0.6

165.0 158.5 152.0 145.5 139.0 132.5 126.0 119.5 113.0 106.5 100.0

0.4

0.2

−0.625

0.000

0.625

Fig. 6.23 Experimental Dryout heat flux for top flooding (DHF ¼ 604 kW/m2) [45].

0.7

Time t = 2134.6s

0.6

Height [m]

0.5 0.4 0.3

TP [K] 392 390 388 386 384 382 380 378 376

0.2 0.1 0.0 0.0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

Radius [m]

Fig. 6.24 MEWA calculation of top flooding with Modified Tung and Dhir model (DHF ¼ 614 kW/m2) [45].

532

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.25 Test section with installed downcomer in the bed [45].

observed, that of rewetting. It can be seen from Fig. 6.26 that, initially, there was dryout observed at a location, which spread downwards, but it was again rewetted by the water flowing up from the downcomer. The simulation results of MEWA code using Modified Tung and Dhir Model as well as Reed model are given in Figs 6.27 and 6.28. It can be seen that both the models not only are unable to capture the rewetting phenomena but also over-predict the enhancement of the dryout heat flux due to bottom flooding. A summary of the experiments and the simulation is given in Table 6.3.

6.4.3.1 Simulation of dryout heat flux: Effect of capillary forces in a stratified porous debris bed The debris bed can have a nonuniform particle distribution or porosity stratification. Tests conducted in the FARO series [46], as well as in KROTOS [47], with prototypic material, showed that the debris bed formed due to melt fragmentation in water, constitutes particles of various sizes ranging from close to 0 mm to slightly more than 6 mm. In the KROTOS tests, most of the particles were found to remain in the range of 0.045–2.0 mm. In the FARO tests (L-24 and L-31), the particles were found to be in the size range between 0.7 and 6.0 mm. Similarly, in the CCM tests conducted at the Argonne National Laboratory [48], most of the particles were found to lie between 0.25 and 10 mm. It is not only the size of the particles but also the porosity of the debris bed, which influences the coolability during the flooding. Similar to the variation of the particle size, the debris bed formed is found to have a spatially non-uniform porosity distribution. The FARO tests found that there were cake regions formed in the

T (°C)

0.6 111.0 110.0 109.0 108.0 107.0 106.0 105.0 104.0 103.0 102.0 101.0

0.5

0.4

0.3

t=9560 0.6

155.0 149.4 143.8 138.2 132.6 127.0 121.4 115.8 110.2 104.6 99.00

0.4

0.2

0.2

0.1

−0.625

0.000 0.625 Diameter [m]

T (°C)

−0.625

0.000 0.625 Diameter [m]

t=9594 0.6

T (°C) 170.0 165.3

t=9610 0.6

155.8

0.5

146.3

0.5

136.9 127.4

0.4

0.4

117.9 108.5

0.3

99.00

0.3

0.2

0.2

0.1

0.1

−0.625

0.000 0.625 Diameter [m]

T (°C) 175.0 167.6 160.2 152.8 145.4 138.0 130.6 123.2 115.8 108.4 101.0

Modelling of core melt scenarios in nuclear reactors

t=9490

−0.625 0.000 0.625 Diameter [m]

Fig. 6.26 Experimental dryout heat flux ¼ 896 kW [45].

533

534

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

0.7

Time t = 8099.6s JLRC (-)

0.6

0.1

Height (m)

0.5 0.4

TP (K) 388 387 386 385 384 383 382 381 380 379 378 377 376

0.3 0.2 0.1 0.0 0.0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

Radius (m)

Fig. 6.27 MEWA calculation of bottom flooding with Modified Tung and Dhir model [45] (DHF ¼ 1370 kW/m2).

0.7

Time t = 3274.5s JLRC (-) 0.08

0.6

Height (m)

0.5 0.4 0.3

TP (K) 406 402 398 394 390 386 382 378 374

0.2 0.1 0.0 0.0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

Radius (m)

Fig. 6.28 MEWA calculation of bottom flooding with REED model (DHF ¼ 1145 kW/m2) [45].

Modelling of core melt scenarios in nuclear reactors

Table 6.3

535

Summary of the experiments with downcomers Experimental

MEWA

Scenario

DHF kW/m2

MTD

REED

Top flooding Bottom flooding

604 896

614 1370

715 1145

middle of the bed where the particles not completely quenched had coalesced together, resulting in a radially stratified bed with lower porosity at the centre of the bed and higher porosity at the periphery. The stratification can happen not only in radial directions but also in axial directions. In the COTELS tests [49], the debris bed formed with prototypic material was found to be like a sand bed. Most of the tests conducted so far have revealed that the porosity of the debris bed lies between 0.2 and 0.4. Numerous experimental studies were conducted in small-scale facilities in the 1970s and 1980s to estimate coolability of particulate beds, e.g. by Hardee and Nilson [50], Dhir and Catton [51], Catton et al. [52], and Reed et al. [53]. They determined the effects of bed height, pressure, particle size distribution, and porosity on the coolability of debris beds. Cho and Bova [54] found that during the top flooding process, the quench front propagates faster in the middle of the bed than in the outer regions. Schrock et al. [55] studied the flooding of an isothermal particle bed with steam and water. Hu and Theofanous [56] investigated the coolability behaviour of deep volumetrically heated particle beds with non-spherical particles. Tung and Dhir [57] conducted experiments in debris beds with both radial and axial porosity distributions. They used steel particles for the debris bed and heated the bed inductively. They found that the quench front velocity is controlled by the counter current flooding limit during top flooding of an axially stratified bed. During top flooding of a radially stratified porous bed, they observed very little crossflow at the boundary between the layers. The quenching was found to occur first in the high porosity region while the low porosity region was found to get quenched later by the bottom flooding mode due to inflow of liquid from the bottom of the quenched high porosity region. However, the size of particles they used in the tests, were much larger than the particle sizes observed in corium fragmentation tests. Tung et al. [58] conducted coolability experiments with top flooding and gas injection from bottom and investigated the bed quenching behaviour. Sehgal et al. [59] studied the quenching characteristics of homogenous and axially stratified sand particle beds with top and bottom flooding using downcomers of different size. They found that the quenching period is decided by the low porosity layer and its mean particle size. The quenching period was also found to be reduced significantly with bottom flooding caused by the water circulating from the water overlayer through the downcomers. Nayak et al. [60] investigated the influence of radial stratifications of the bed on the quenching behaviour using smaller-size particles with prototypic porosity distributions. For this, they performed experiments in a facility named as POMECO

536

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Steam Downcomers Low porosity (0.26) Water Side pipe injection

Large downcomer Small downcomers

High porosity (0.38)

High porosity (0.38)

48 Thermocouples Sand bed

Air/argon injection

Fig. 6.29 POMECO facility [60].

(POrous MEdia COolability). It basically consisted of two stainless steel tanks with the lower tank being the test section and the upper tank used for the water overlayer providing the top flooding. Both had the same cross-section of 350-mm square. The test section had a height of 500 mm and the height of upper tank was 900 mm. Seven downcomers were placed in the test section in order to bring water from the top tank to the bottom of the bed for bottom flooding purpose (Fig. 6.29). The corium porous debris bed was simulated by using river sand particles with proper size distribution. The size and porosities of the sand particles used in the experiments were close to those of the corium debris bed formed in an accident situation. The radially stratified bed was prepared by putting the high porosity sand with porosity of 0.38 at the periphery and the low porosity sand of porosity 0.26 at the centre of the bed. The sizes of sand particles ranged from 0 to 5 mm. The sand bed was heated volumetrically with closely spaced electrical heaters and was flooded with water from the top, and from both top and bottom when downcomers are employed. The quenching time periods for different conditions were measured and compared with those of homogeneous and axially stratified beds. For the case with top flooding only, the DHF value measured in dryout experiments was 112 kW/m2, which is between the values of 45 kW/m2 and 222 kW/m2 measured in earlier experiments for homogeneous beds with porosity 0.26, particle size 0.8 mm and porosity 0.38, particle size 0.92 mm, respectively. In the experiment dryout first was detected in the region with higher porosity (Fig. 6.30). Fig. 6.31 shows the results of computation where the capillary force has been neglected. The dryout first occurs at the upper left corner in the low porosity layer. A kind of natural circulation loop establishes, with water inflow from top into the layer of higher porosity and supply of water to the lower porosity layer mainly via the side and the bottom. The calculated dryout heat flux is 72 kW/m2, which is too small compared to the experimental value of 112 kW/m2. Further, the calculated dryout location contradicts the experimental findings, where dryout was first detected in the layer with higher porosity.

Modelling of core melt scenarios in nuclear reactors

537

Fig. 6.30 Temperature history in high porosity (left) and low porosity (right region) [41].

Fig. 6.31 Dryout prediction without capillary forces for top flooding conditions [41].

Fig. 6.32 shows the result with capillary forces included. In this calculation, the Leverett function according to Hofman and Barleon (Eq. 6.23) was applied. As the wetting properties of the debris in the experiment was not known, the cosine of the contact angle was set to 1, i.e. the maximum wetting conditions were assumed in the calculations in order to maximize the capillary effects. Here, as in the experiment, a dry zone occurs first in the high porosity zone. The mechanism can be

538

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.32 Dryout prediction with capillary forces for top flooding conditions [41].

understood from Fig. 6.32. Water flowing into the higher porosity layer is directly drawn to the lower porosity layer by capillary forces, due to higher capillary pressure in the lower porosity layer (liquid pressure is lower there) and no loop establishes.

6.4.4 Summary Due to the presence of porous medium, there is a substantial change in the mass, momentum, and energy equations due to additional terms of porosity, permeability, and interfacial drags between solid, liquid, and gas in the conservation equations. However, to a great extent, empirical models have been developed for interfacial drag, which need to be tested with a large variety of experimental data. Nevertheless, these models are found to be capturing the physics. In a laterally stratified bed with strong differences in the solid matrix (porosity, particle diameters), dryout can even occur in the more permeable region (higher porosity region). This phenomenon is explained by the capillary forces moving water from the high porosity to the low porosity region to such an extent that the cooling limit is reached in the high porous layer first. For prediction of pressure drop, the interfacial drag plays an important role. Schulenberg and M€ uller model with explicit consideration of interfacial drag between liquid and vapour phases is more suitable for top-flooding and low inflow rate bottom-

Modelling of core melt scenarios in nuclear reactors

539

flooding, whereas the models without interfacial drag like Reed model can better predict bottom-flooding with relatively high inflow rates.

6.5

Modelling of melt coolability in ex-vessel conditions

During a severe accident, with failure of the reactor vessel, the molten corium can be relocated in the containment cavity forming a melt pool. The melt pool can be flooded with water at the top for quenching it. The effectiveness of heat removal from the melt pool shall depend on the water ingression depth inside the melt pool to cool and quench it. The water ingression into the melt pool is a complex thermal hydraulics problem, which depends on the thermal and mechanical properties of corium such as the viscosity, the melt temperature, and thermal conductivity, which can change drastically due to its interactions with the in-vessel and ex-vessel structures and also with the concrete basemat. Similarly, the modulus of elasticity, tensile strength, and linear thermal expansion coefficient, which are some of the key parameters that determine the crack formation during melt quenching and responsible for crust cracking and water ingression, can change depending on the amount of structural and concrete materials mixing with corium. Originally, the motivation behind modelling of water ingression phenomenon did not aim to study the quenching of molten corium pool during a postulated severe accident condition in an LWR. In fact, models were developed for simulation of cracking behaviour of hot rocks in geological reservoirs. In this context, Lister [27] has done pioneering work in modelling the penetration of water into hot rocks by considering the simplest possible one-dimensional model based on the concept of crack front propagation. Bj€ ornsson et al. [61] found that penetration of water into hot rock is the primary reason for the intense heat release of the subglacial Grimsv€otn geothermal area. Jagla and Rojo [62] presented a model to predict the statistical properties of columnar quasi-hexagonal crack patterns, as observed in the columnar jointing of basaltic lava flows. While Lister’s model considered penetration of water into hot but initially solid rock under high-pressure condition, recently, Epstein [28] used Lister’s models of bulk permeability of cracked rock and developed a model for water penetration into initially molten, heat generating rock like material at low pressure, which resembles the water ingression phenomena into molten corium pool. In addition to difficulty in modelling the water ingression phenomena, there are only a few experiments conducted on quenching behaviour and water ingression in top flooded molten pool. In the MACE experimental programme [63], it was found that a tough crust is formed on the upper surface of the melt pool during top flooding situation, which was found to limit the access of the water overlayer to the melt pool below the crust. Posttest examination of the debris of the MACE M3b-test indicated that the crust thickness was about 10 cm, amounting to 1/2 of the initial mass [64]. After 20 min of melt/water heat transfer in M3b, several independent sources of data indicate that the crust had anchored to the test section side walls and the melt separated from the crust [65]. After separation, melt/water heat transfer rate dropped significantly. Test results from COTELS project [49] indicate that water ingression through cracks/defects in core material interacting with concrete can contribute to debris

540

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

coolability. Water can penetrate into debris at the sidewalls due to erosion at this interface as well as direct penetration of water into channels located in the central regions of the debris. However, the authors did not quantify the relative contributions of water ingression at the sidewall core/concrete interface versus ingression in the central core material region to the overall debris coolability. Melt coolability research was also performed at the Sandia National Laboratory [66] and WETCOR [67] experiments with stainless steel and corium melts, respectively. High-pressure corium quenching tests were performed in FARO facility [68]. A series of low-temperature simulant material experiments were conducted by Theofanous et al. [69] at the University of California, Santa Barbara. Also, the melt coolability research has been pursued in Germany [70] in the COMET project where large-scale (several hundred kilograms) melt pools were cooled and stabilized by injecting water at the bottom of the pool. Farmer et al. [71] have reported the test results of the OECD-sponsored MCCI programme on ex-vessel debris coolability and the two-dimensional molten core–concrete interactions under both wet and dry cavity conditions. They used fully oxidized 400 kg of PWR core melt, initially containing 8 wt.% calcined siliceous concrete in a specially designed two-dimensional siliceous concrete test section. The core–concrete interaction results indicated that radial ablation is a key element of overall cavity erosion process. Late-phase flooding was found to cause a significant increase in upward heat flux from the debris. Crust breach by an insertable lance caused a dramatic increase in debris cooling rate. Recently, Lomperski and Farmer [72] reported the results of tests carried out to assess the significance of water ingression cooling in the quenching of molten corium. They conducted tests for quenching of about 75 kg of melt at 2100°C consisting of UO2, ZrO2 and chemical constituent of concrete in the range 4%–23%. Their results indicated that, for the lower concrete content melts, cooling rates exceeded the conduction-limited rate with the difference being attributed to the water ingression mechanism. In this context, a mathematical model has been developed by Nayak et al. [73] to investigate and clarify the water ingression and melt coolability behaviour when the melt pool is flooded from the top in an ex-vessel situation. The model has been developed with the following postulations:

6.5.1 Postulation of melt coolability phenomena l

l

l

l

As soon as water is introduced in the pool, a vapor film is formed over the melt surface (Fig. 6.33). The heat transfer from melt pool to water overlayer is by convection and radiation. Heat conduction is dominant mode of heat transfer inside the melt pool. As the melt is cooled at the top, a solid crust is formed at the top of molten surface due to cooling. The crust sees low temperature at water side and high temperature at the melt side (Fig. 6.33), resulting in thermal stress development. When the thermal stresses exceed the fracture stress, the crust becomes mechanically unstable, it disintegrates into debris. Water ingresses into the debris by gravity. It takes the heat away by convection. During this heat transfer, there is a counter-current flow of steam and water into the debris bed. At high

Modelling of core melt scenarios in nuclear reactors

Overlying water pool

541

Fig. 6.33 Domain of postulation for melt coolability under top flooding.

Vapor film Debris bed Solid crust Melt pool

l

heat flux, the upward steam flow may limit the water downflow leading to counter current flooding limit (CCFL). This would ultimately cause dryout. Once the dryout condition has occurred, water would not ingress into the debris.

6.5.2 Development of numerical model The model solves the conservation equations of mass, momentum, and energy from the molten pool to the overlaying water (Fig. 6.33) as described here: (i) Governing equations for molten pool

In molten pool, in the containment condition and near solidification temperature, the viscosity of the melt is several orders of magnitude higher than that for normal fluids like water, mostly due to inclusion of concrete and other structural materials in it. Hence, it is assumed that the dominant mode of heat transfer is due to conduction only. With this, the transient 2D heat conduction equation can be written as ρCp


 ∂T 1 ∂ ∂T ∂2 T 0 00 ¼ r k +k 2 +q ∂t r ∂r ∂r ∂z

(6.27)

where q000 is the volumetric heat generation rate due to decay heat. The boundary conditions are shown in Fig. 6.34. (ii) Governing equations in solid crust:

Heat transfer in solid crust is due to conduction only. The energy equation can be written as ρCp


 ∂T 1 ∂ ∂T ∂2 T 0 00 ¼ r k +k 2 +q ∂t r ∂r ∂r ∂z

The boundary conditions for this region are shown in Fig. 6.35.

(6.28)

542

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.34 Boundary conditions for melt pool region.

r T = Tmelt z r

R r

r

Melt pool

0

0 r

z

0 z

0

0

Fig. 6.35 Boundary conditions for crust region. (iii) Governing equations in debris region:

Mass balance equation rν¼0

(6.29)

Momentum balance equation κ v ¼ ðr  p + ρgÞ μ

(6.30)

Energy balance equation: γ1

∂T q00

+ v  rT ¼ αr2 T + ∂t ρCp f

(6.31)

Here q000 is heat generation term in debris region (in debris region water or steam is present between particles so here heat generation can be different from crust) and, γ 1 ¼ ε + ð1  εÞ

ρCp

ρCp



s f

The boundary conditions for this region are shown in Fig. 6.36.

(6.32)

Modelling of core melt scenarios in nuclear reactors

543

Fig. 6.36 Boundary conditions for debris region.

Here ψ is the stream function defined as νr ¼

1 ∂ψ 1 ∂ψ ;νz ¼  r ∂z r ∂r

(6.33)

Calculation of surface heat transfer coefficient: For calculation of heat transfer from top of the debris to the pool of water, appropriate heat transfer correlations are used as given here. If surface heat flux is greater than the debris bed dry out flux, combination of radiation and film boiling heat transfer is used, i.e. for q00surface > q00d ; h ¼ hfb + hrad

(6.34)

The film boiling heat transfer coefficient is calculated using Berenson’s [74] film boiling model as 0 11=4 Bhfg  ρv  g  ðρl  ρv Þ  kv3 C rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C hfb ¼ 0:425B @ A σ ST μv   ΔT g  ð ρ l  ρv Þ The radiation heat transfer coefficient is calculated as

2 ðT + Tsat Þ hrad ¼ σ stefan  εs  T 2 + Tsat

(6.35)

(6.36)

When the surface heat flux falls below the dryout flux, i.e. If q00surface < q00d ; h ¼ hnb

(6.37)

For this condition, nucleate boiling is considered and the heat transfer coefficient is calculated from the nucleate boiling correlation by Rohsenow [75] as given by  3

2 CPl 1:25  105 μl hlv ΔTsup hlv hnb ¼  5:1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi CPl  μl σ ST Kl gðρl  ρv Þ

(6.38)

544

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

The dryout heat flux was calculated using the following relationship based on counter current flooding limitations (CCFL) as given by Sinha et al. [76]. 1=2

qd ¼ 0:539hlv

½ðρl  ρv Þρv gdε3 =ð1  εÞ h i2 1 + ðρv =ρl Þ1=4

(6.39)

Solid crust generation rate: Crust generation rate is calculated from energy balance as given here: ðð

  ! ∂T  ∂T  km  + kc  dA  Δt ∂z pool ∂z crust 00 0

¼ A  Δzcrust  ρs  Cps  ΔT + A  Δzcrust  ρs  hfus + q  A  Δzcrust

(6.40)

There is no heat flux continuity at crust–melt interface. This heat flux discontinuity is responsible for the crust growth. If the left-hand side of the equation is negative, it implies that crust is dissolving in molten pool. This will occur due to high heat generation rate in pool, which leads to increase in pool temperature, whereas top flooding is unable to take the heat away.

6.5.2.1 Stresses in solid crust As cooling initiates, there will be a temperature distribution in the crust. Initially, the crust will be formed only at the melting temperature of the material. At the top of the crust, the crust gets cooled and is at water saturation temperature. At the bottom, the crust is at melting temperature. Hence, thermal stresses are induced inside the crust. The 2-D axi-symmetric equations for stress developed in solid crust are: ∂ ∂σ z ¼0 ðr  τrz Þ + ∂z ∂r

(6.41)

∂ ∂τrz ¼ σθ ðr  σ r Þ + r ∂z ∂r

(6.42)

These equations straightway come from force balance in differential element. The boundary conditions are given as follows (Fig. 6.37): Fig. 6.37 Stresses in the crust and boundary condition.

trz = 0

sz = 0 u=0 ∂w =0 ∂r

u=0 w=0

trz = 0

w

w=0

u

Modelling of core melt scenarios in nuclear reactors

545

Where u is radial and w is the axial displacement. Stress strain relationship is given as 0 10 1 1 1ν ν ν 0 σr εr B ν 1ν ν 0 C B σθ C B εθ C E B C B C¼ C ν 1ν 0 CB @ σ z A ð1 + νÞð1  2νÞ B @ A @ ν 1  2ν A εz 0 0 0 τrz γ rz 2 0

(6.43)

So we can write 0

1 ∂u 0 1 B ∂r  αl ΔT C B u C εr B C B εθ C B r  αl ΔT C B C¼B C @ εz A B ∂w C B  αl ΔT C B C ∂z γ rz @ ∂u ∂w A + ∂r ∂z

(6.44)

From these equations, ultimately two equations in u and w are obtained. At the top of the crust, there will be no shear stress and the axial stress due to self-weight of the debris is considered to be negligible. These equations were discretized in implicit manner and solved with boundary conditions as shown in Fig. 6.37, to obtain stresses distribution in crust.

6.5.2.2 Criteria for fracture of the solid crust The solid crust, being a ceramic mixture of the oxides, behaves like brittle material. Brittle materials have no yield limit and they fracture catastrophically at the peak stress, which is termed as fracture stress. For fracture initiation, Ritchie’s model [77] was used. In this model, fracture initiation occurs only when thermal stress exceeds critical stress in some region. This region is taken generally of the order of eight-grain size. So fracture initiation occurs only when stresses are above critical limits in more than eight-grain size area. Since this requires detailed knowledge of dimensions of critical region (like eight-grain size), it may not be known for some materials. Hence, an alternative model is used to simulate fracture initiation, which is known as Beremin’s [78] model. In this model, the critical stress σc required for the separation of cleavage facets, which can be related to the length ‘l0’ of a micro-crack by the equation as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Eγ s σc ¼ π ð1  ν2 Þl0

(6.45)

where E is Young’s modulus of elasticity, ν is the poison’s ratio, and γ is the surface energy. The length ‘l0’ is maximum size of micro-crack present in some volume. This

546

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

is a probabilistic model and distribution of micro-cracks is required for the longer ones. In each volume, the probability of finding a crack of length between l0 and l0 + dl0 can be written as: Pðl0 Þdl0 ¼

α1 β

l0 1

dl0

(6.46)

where α and β are material constants for a particular value of V0. Hence, in a given volume for a stress level σ, the probability of failure is ∞ ð

Pðσ Þ ¼

Pðl0 Þdl0 ¼ lc0

α1 1

β1  1 lc β1 1

(6.47)

0

So from some stress level σ, the critical length of micro-crack needed is obtained from Eq. (6.45) and the probability of failure of the concerned grid is calculated using Eq. (6.47). The total probability of failure is

 σw PR ¼ 1  exp  σu

(6.48)

where σ u is material constant and is given as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n

X m Vi m σ i1 σu ¼ V0 i¼1

(6.49)

In Eq. (6.49), Vi is the volume of the ith element experiencing the maximum stress σ i. In brittle materials when fracture occurs, it moves catastrophically. In lower points of crust, temperatures are as high as melting point; so crack may be arrested. To model this, the R6 model is used as given by R6 model: ( Kr ¼



 1  0:14L2r 0:3 + 0:7 exp 0:65L6r 0

for Lr > Lmax r

where Lr ¼

Applied load Limit load based on yeild stress

Kr ¼

Stress intensity factor Fracture toughness

forLr  Lmax r

(6.50)

Modelling of core melt scenarios in nuclear reactors

547

1.2 FAL

1

Lr

0.8

B

0.6 0.4 A

0.2 0 O0

0.2

0.4

0.6

0.8 Kr

1

1.2

1.4

1.6

Fig. 6.38 R6 model curve.

The curve is given in Fig. 6.38, if at some time Lr and Kr lie under curve shown (FAL), it indicates that crack is arrested. The governing equations in melt pool and solid crust region are discretized using finite difference method and solved implicitly using Gauss elimination method to obtain temperature distribution in the melt pool. After evaluating the temperatures, the crust growth rate is calculated and subsequently the thickness of the solid crust region is updated. In debris region, pressure term is eliminated in the momentum equation with the help of stream function approach. Then applying Boussinesq approximation, we get momentum equation in ψ and T only. Energy equation is also reduced in terms of ψ and T. Then both the equations are solved implicitly. After the temperature distribution has been obtained, the stress equations are solved implicitly by finite difference technique. With the stresses, the fracture conditions are evaluated. If criteria for fracture are satisfied, the fractured grids are merged into debris region. The calculation procedure for debris surface heat transfer coefficient has already been discussed. When the surface heat flux falls below the dry out flux, water is considered to ingress into debris; otherwise, dryout occurs.

6.5.3 Model validation with experiments For validation of the model, experimental data from COMECO [79] were taken. The COMECO facility consists of a test section (200 mm  200 mm cross-section) with a maximum pool height of 300 mm. Fig. 6.39 shows a schematic of the facility. The test section walls are made of carbon steel plates with thickness of 25 mm. The test section was connected to an upper tank whose height is about 1000 mm and its cross-section same as that of the test section. The upper tank was used for the purpose of water flooding on the top of the melt pool to the desired height. For this, water was preheated to a suitable temperature in the water storage tank before it was delivered to the upper

548

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

Fig. 6.39 Schematic of COMECO test facility. Venturi meter Water tank Water

Melt from furnace

Melt

Heaters

Test section

Thermocouples

tank via the water line. The level of the water in the upper tank was kept at around 700 mm during the tests. The melt pool was heated by heaters, located outside the test section on the four sidewalls. The heaters were made of molybdenum silicate (MoSi2) alloy that can be operated at temperature of about 1700°C. The four heaters could deliver a total power of 16 kW (4 kW/heater) to the melt pool, which corresponds to a heat generation rate of nearly 1.33 MW/m3. This gives an average heat flux of 66.7 kW/m2 in each face of the test section. The power input to the test section was constant throughout the test. Fig. 6.40A shows transient temperature history of the molten pool during quenching, as predicted by our model. Theoretical results also show that temperature of top layer falls suddenly in 500 s, which is closer to measured value. The transient temperature history due to quenching is similar at all the depths as in the experiments. Model predicts water ingression of 0.11 meters, as shown in Fig. 6.40A. It is very close to experimental result, which was around 100 mm (Fig. 6.40B). In Figs 6.41 and 6.42, the total depth of water ingression, crust height, and dry debris have been plotted with respect to time. It can be seen that water ingression rate is faster initially. During this time, the crust height is found to be fluctuating. This is due to the fact that initially whatever crust forms, disintegrates into debris and the debris gets quenched due to water ingression into it. Since brittle fracture is the mode of crust disintegration as assumed in the model, the debris formation occurs in steps, i.e. when crust disintegration stops, depth of crust increases again periodically. There arises a big question ‘why water ingression stops after some depth?’ As depth of water ingression increases, it means depth of debris flooded with water increases. This water takes the heat produced, i.e. volumetric heat generation due to decay heat. However, due to increased debris depth as disintegration of crust proceeds, the path

1200

o

Temperature ( C)

1000 800

Distance from the Top of the melt pool 1 (27 mm) 2 (47 mm) 3 (70 mm) 4 (122 mm) 5 (170 mm) 6 (222 mm) 7 (274 mm)

600 400 200 0 0

250

500

750 1000 1250 1500 1750 2000 Time (s)

(A) 1200 1000 Temperature (°C)

Distance from top of melt pool 800

1 (27 mm) 2 (47 mm) 3 (70 mm) 4 (122 mm) 5 (170 mm) 6 (222 mm) 7 (274 mm) 8 (322 m)

600 400 200 0 400

0

800

(B)

1200

1600

2000

Time (s)

Fig. 6.40 (A) Model predictions; (B) COMECO test results.

0.12 0.10

Height (m)

0.08 0.06

Water ingression Dry debris height Crust height

0.04 0.02 0.00 0

1000

2000 Time (s)

3000

4000

Fig. 6.41 Depth of water in ingression and crust height predicted by model.

550

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

0.30

Water ingression 0.25

Dry debris

Bed height (m)

0.20

Solid crust 0.15

0.10

Melt pool

0.05

0.00 0

1000

2000

3000

4000

Time (s)

Fig. 6.42 Various regions in pool vs time.

resistance increases for water to ingress down to cool the top of the crust, which is at the bottom of debris. So the crust top surface temperature progressively increases as time proceeds. Since bottom part of crust is always at liquidus temperature, the thermal gradient inside the crust becomes lower and lower with time, which reduces the thermal stresses inside the crust. After some time, the stresses become so low to disintegrate the crust further. That is why water ingression stops after some depth. In addition, with reduction in thermal gradient across the crust, the heat absorption rate from molten part also slows down, which reduces the solidification rate. After water ingression stops, crust depth increases slowly. Increasing depth of crust having low conductivity material acts as a blanket for molten material filled below. So heat transfer from molten part decreases. It leads to very low crust formation rate. After some time, when net heat produced in molten pool equals to heat taken away by the overlaying water, crust formation stops. And these lead to an eventual steady state.

Discussions on current model vis-à-vis Lister’s model The model developed by Lister is based on transient creep, which primarily depends upon strain rate and is definitely appropriate for the large scale and long time frame events like geological phenomena wherein the rocks are hot and under pressure due to loading of several tones of cold water lying over them. The hot rocks are considered to crack mainly due to creep and the crack front moves progressively over prolonged

Modelling of core melt scenarios in nuclear reactors

551

period. Whereas in the present system, the crust developed due to quenching of molten material, cracks only when the thermal stresses built-in it exceeds the critical stress for fracture. Lister’s model is based on strain rate. Strain rate is calculated from an initially defined crack propagation velocity, which is continuous function of time. However, from the COMECO experiments and present simulation results, we observed that the crack propagation is discrete and sudden until a depth of 100 mm below which no water ingression is observed due to presence of solid rock-like structure. This is an indication of the material behaving like brittle materials such as glass, which are typical to the crack formation in brittle materials. Further, the COMOCO tests indicated that the crack occurs only a few times during the quenching process and is not continuous over the entire period of quenching. In fact, most of the time, the crack propagation velocity was zero. In view of the different mechanics for crack formation, it can be presumed that Lister’s model is not applicable for water ingression in molten corium pool at ex-vessel situation. In addition, as observed in the present simulations and COMECO tests, during the first 100 s, the temperature drop is relatively higher. Hence, it would be presumed that the stresses will be higher during this period. In the subsequent, 400 s, the temperature drop rate is found to be negligible as compared to initial 100 s. According to Lister’s model, there should be negligible stress in this period as compared to that during the initial period since the stress is calculated from the strain rate, which, in turn, is based on rate of change of temperature. Hence, cracks should have occurred in the crusts during the initial 100 s. However, fracture is found to occur during the later period as indicated by large bursts of steam (Fig. 6.43) when cracks occur in the crusts and water ingresses even up to 1000 s (Fig. 6.40B). So the phenomena contradict the assumptions in the model of Lister.

Fig. 6.43 Measured steam flow rate during quenching of top flooded melt pool in COMECO.

552

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

6.5.4 Numerical simulation of melt coolability under bottom flooding In bottom flooding, the melt is allowed to fall on a porous surface so as to allow the water to be fed from bottom through some nozzles. Mostly, passive means are used to allow water flooding [80] where water is fed from gravity driven water tank. At the melt water interface, a sacrificial material is placed. When the melt falls on the sacrificial surface, it gets ablated and the water is then forced upwards by gravity. In some cases like COMET [79], the melt falls on the porous material itself and then water is fed through the porous material. In the present model, this situation is considered where the melt falls on the porous material and then water is fed from bottom at fixed rate. We postulate that, when melt falls on the porous material, it starts getting cooled and a thin crust starts getting formed. When water supply is continued, water reaches the melt from below and starts further cooling the melt and the crust starts growing. As a result of heat transfer from melt to the water, steam starts getting formed and it expands in volume, which exerts back-pressure to water. As a result, steam pressure builds up in the area below the melt crust. Hence, the crust is subjected to stresses mainly because of the water hydrostatic head and steam pressure on one side, hydrostatic head of melt pool on the other side, and thermal stresses as a result of temperature gradient across the thickness. When the stresses in the crust exceed the fracture stress, the crust breaks resulting in debris formation. The steam passes through the debris continuously pressurizing it, which ultimately causes eruption at selected sites. As a result of eruption, the steam escapes through the melt, cooling it instantly and forming porous passage. The erupted melt forms debris-like structure on the top of the melt. The eruption sites act as passage for water and steam, which brings about a multi-dimensional coolability. Based on this phenomenology, an integrated model has been developed for each phenomenon taking place as postulated here. Fig. 6.44A–D shows the phenomenology of the melt coolability under bottom flooding. Fig. 6.44A shows the scheme for bottom flooding. In Fig. 6.44B, thin crust generation and steam pressure built-up is depicted. When the pressure exceeds the fracture stress, the crust breaks and porous zone is formed and ultimately, melt eruption occurs (Fig. 6.44C). As a result of melt eruption, the coolability is enhanced as shown in Fig. 6.44D. Based on these phenomena, models are developed for each phenomenon.

6.5.4.1 Governing equations (i) Governing equations for molten pool

In molten pool, it is assumed that the dominant mode of heat transfer is due to conduction only. With this, the transient 2-D axi-symmetric heat conduction equation can be written as ρCp


 ∂T 1 ∂ ∂T ∂2 T 0 00 ¼ r k +k 2 +q ∂t r ∂r ∂r ∂z

(6.51)

Modelling of core melt scenarios in nuclear reactors

553

Melt pool Crust formation

Porous material

Steam formation

Water inlet

(A)

(B)

Enhanced coolability

Melt eruption

(C)

(D)

Fig. 6.44 Phenomenology of melt coolability under bottom flooding.

Similarly, in crust layer, energy equation is given as ρCp


 ∂T 1 ∂ ∂T ∂2 T 0 00 ¼ r k +k 2 +q ∂t r ∂r ∂r ∂z

(6.52)

The boundary conditions for these equations are given in Fig. 6.45. At the top, heat loss from radiation is considered. The emissivity of corium is considered to be 0.87 [32]. (ii) Equation in porous zone:  Velocity is predicted by Darcy equation as given below

v¼

dp κ dz μ

(6.53)

Rate of vaporization is given as

V  A  ρ  hin  hfs + h  A  ðTmelt  Tsat Þ Γ¼ hfg

(6.54)

554

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

T

k

T

4

r

T

k

4

Tsat

T

k

0

k

T r

0

r

r

0

k

T

0

r k

T r

h (T Tsat )

Fig. 6.45 Boundary conditions for melt and crust layer.

Steam generated in a given time is Mst ¼ Γ  dt

(6.55)

Specific volume of the steam ϑ is ϑ ¼ Mst =Vp

(6.56)

From equation of state, the pressure developed is calculated as p ¼ f ðϑÞ

(6.57)

(iii) Stresses in the crust:

The crust is subjected to stresses as shown in Fig. 6.46. At the top of the melt, there is hydrostatic head and at the bottom, the pressure is exerted by the steam. In addition to this, the top edge of the crust is at melting temperature, whereas the bottom end is at much lower temperature. This exerts thermal as well as mechanical stresses. The bending stresses on the circular-plate-type crust as a result of clamped edges are given as σ b ¼ 0:75 

ðΔp  r 2 Þ t2cr

(6.58)

Modelling of core melt scenarios in nuclear reactors

s z = − rgH

555

melt

u=0 w=0

u=0 w=0 s z= − p

w=0

Fig. 6.46 Stresses in the crust.

The thermal stresses as a result of temperature gradient are given as σ th ¼

α  ΔT  E 2ð 1  ν Þ

(6.59)

The total strain is additive and the material, being brittle, remains in elastic region. Hence, we can add the individual stresses to obtain the total maximum stress acting on the crust as σ tot ¼ σ th + σ b

(6.60)

The crust will break if the total stress exceeds the strength of the crust, i.e. σ tot > σ max

(6.61)

After the crust is broken, melt eruption takes place at certain locations. The location of melt eruption sites is random in nature and it involves certain empiricism to predict the exact location of eruption sites. However, we can estimate the distribution of the sites in the entire area. Paladino et al. [29] have developed an empirical model for density of the openings and their diameter. The average number of channels is given as Ψ¼

 Γ Cpw  ΔTsub + hfg + Cpst  ΔTsup c hpool  Ach  Nn  ΔTm, w

(6.62)

where ΔTm, w is average temperature drop between melt and water. This parameter is a measure of heat transfer enhancement due to formation of small porous openings, which enhance the area available for heat transfer. The average density of channels per unit area can be given as γ¼

1 ξ2

(6.63)

556

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

where ξ is the spacing between the channels given by Zuber’s [81] modified Critical Taylor Wavelength formula as ξ ¼ 2π

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ mc ðρm  ρst Þg

(6.64)

The diameter of the eruption area is given as d¼2

rffiffiffiffiffi ψ γπ

(6.65)

The resulting porosity is given as: ε¼

Ψ  π  D2n  hpool π  L2b  hpool

(6.66)

Once the crust is broken, the melt layer is considered to be made of two zones, mainly porous zone and nonporous zone. We know the diameter, porosity, and the number of openings. We can predict the coolability by modifying the initial domain and governing equations. The equation for porous medium is modified as
 ∂T ∂T 1 ∂ ∂T ∂2 T 00 0 +v ¼ r keff + keff 2 + q (6.67) ρeff Cp, eff ∂t ∂y r ∂r ∂r ∂z where ν is the velocity of steam flowing through the bed calculated using Eq. (6.53) and the flow rate, which is kept fixed during bottom injection. The effective properties are volume averaged over the void fraction in bed with boundary conditions as keff

∂T ¼ h ð T  T∞ Þ ∂z

(6.68)

where h is evaluated from surface temperature based on heat transfer regime, i.e. film boiling or nucleate boiling. For film boiling, Berenson’s [74] model is used and for nucleate boiling, Rohsenow [75] correlation is used. For the solid zone, the boundary conditions are modified as shown in Fig. 6.47. In addition to top and bottom, now the span of the solid zone has reduced and additional convective boundary condition at one side has been introduced, which makes it coolable from two dimensions.

6.5.4.2 Validation of the model with bottom flooding experiment To validate the model, experimental data from DECOBI tests were used [29]. The details of the tests are: l

l

Simulant material: CaO–B2O3 Melting temperature: 1178 K

Modelling of core melt scenarios in nuclear reactors

k

Porous zone

557

T h T Tsat z

Solid zone

Heat loss T k h T Tsat z k

T h T Tsat z

Fig. 6.47 Modified boundary conditions for the domain.

l

l

l

l

l

l

Volume: 6 L Melt height: 14.2 cm Melt temperature: 1400 K Water injection rate: 0.014 kg/s Inlet temperature: 27°C Decay heat: Not simulated

Fig. 6.48 shows the experimental and predicted temperatures in the central region of the pool and Fig. 6.49 gives the pressure predicted by the model. It can be seen that the pressure and temperature are predicted by the model are in very good agreement with the experimental data. Table 6.4 shows the predicted diameter of the eruption zone and peak pressure. The model was able to predict the diameter of eruption zone with good precision.

Centre region

1200

Temperature (°C)

1000

150 mm predicted 150 mm experimental 100 mm predicted 100 mm experimental

800 600 400 200 0 400

500

600

700 Time (s)

800

Fig. 6.48 Comparison of predicted temperature with experiment.

900

1000

558

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

1.8

1.7

Experimental Predicted

Pressure (bar)

1.6

1.5

1.4

1.3

1.2 500

1000

1500 Time (s)

2000

2500

Fig. 6.49 Predicted pressure history.

Table 6.4

Predicted parameters

Parameter

Experimental

Predicted

Error

Diameter of eruption zone Maximum pressure pulse

130 mm 1.57 bar

126 mm 1.67

3.1% 1.7 %

The progression of coolability under bottom flooding as predicted with this model is shown in Fig. 6.50. It can be observed that, at 6 s after the water is introduced, there is eruption taking place. Once the eruption takes place, the water starts cooling the melt rapidly. Within 2000 s, the melt is quenched which otherwise takes hours to cool if cooled from top alone.

6.6

Closure

Severe accidents in nuclear reactor involves many complex phenomena like core damage, fission product release, fuel meltdown, vessel failure, melt pool formation, molten core corium interaction, steam explosion, hydrogen generation, distribution and deflagration in-vessel and ex-vessel corium coolability and debris bed coolability, etc. It is very difficult to solve the entire phenomenology of severe accidents using CFD, owing to multi-phase and multi-physics nature, and lack of complete understanding of physics behind it. However, individual phenomena have been mechanistically modelled with many simplistic assumptions. In this chapter, simulation of a few of these phenomena has been discussed. CFD simulation of in-vessel retention

Modelling of core melt scenarios in nuclear reactors

t=5 s

t = 100 s

t=6 s

t = 1000 s

559

t = 10 s

t = 2000 s 1473 1373 1273 1173 1073 973 873 773 673 573 473 373

Fig. 6.50 Progression of coolability under bottom flooding.

of corium in PHWRs has been presented. In-vessel retention simulation includes heat transfer by convection, conduction, and radiation associated with solidification of corium and properties changing with temperature. These phenomena were in-general well captured using CFD. The results are benchmarked with experimental data, which validates the models used in the simulation for the specific case. However, in order to generalize the simulation for a wider range of simulation, improvement in few areas is needed. For example, the simulation predominantly uses boiling boundary condition on outer wall of the vessel. The correlations used for obtaining the heat transfer coefficient needs to be validated with experimental data on large vessels. In addition, presently in this chapter, initially solid debris has been considered which

560

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

ultimately melts and forms a convective pool. For initial period, when the viscosity of the melt is very high, laminar model is applicable. However, as time progresses, the melt becomes superheated and causes turbulent natural convection. Simulation of turbulent natural convection in such large geometries is a challenge. Further, the density variation of corium with temperature is very important as the natural convection is involved. The data for change in properties with temperature need to be established properly to give accurate predictions. Coolability of heat generating debris bed is important in context of wet containment when the molten pool fragments and forms a particulate bed. Due to the presence of porous medium, addition of solid particle– vapour, and solid particle–liquid drag, adds to the complexity. However, empirical models for drag forces, validated with experimental data, predict the most important parameters in coolability of debris bed; dryout heat flux and two phase pressure drop. However, there is not a single model which can predict the pressure drop and dryout heat flux over entire range of vapour flow rates. There is a scope for complete 3 phase (debris, water, and steam) CFD simulation of debris bed coolability. Ex-vessel coolability of molten corium is the most complex phenomenon in severe accident. Due to multi-physics (thermal hydraulics and mechanics) and multi-phase (solid, liquid, and vapour) nature of the problem, solving it with CFD makes it extremely difficult. On the other hand, postulation of phenomena with simplistic assumptions, using governing equation wherever applicable and plugging empirical correlations helps it bringing out physics of the phenomena. The above method was applied for modelling of melt coolability in ex-vessel conditions in this chapter as examples. Validated models based on postulated phenomena have been presented for corium coolability under top and bottom flooding. The model for top flooding considers the heat transfer behaviour in axial and radial directions from the molten pool to the overlaying water, crust generation and growth, thermal stresses built-in the crust, disintegration of crust into debris, natural convection heat transfer in debris and water ingression into the debris bed. Assumptions such as when the stresses in the crust exceed the strength, the crust breaks and forms debris, have been considered in the simulation. In reality, the prediction of crack formation due to thermal stresses, and its growth is really challenging. Further, modelling of water ingression into such crack and occurrence of dryout is beyond the capabilities of present CFD models. Similarly, in modelling of melt coolability under bottom flooding, simplistic assumption were considered like when steam eruption occurs, the melt is converted into porous bed with a certain porosity and uniform particle diameter obtained from experiment. In reality, simulation of conversion of melt into solid fragments needs in-depth understanding of physics and rigorous modelling which needs to be done in future. Steam explosion is another such phenomenon which includes multiple stages like melt jet breakup, fragmentation of melt, triggering, expansion, and shockwave propagation. Most of the efforts have been dedicated to simulate the jet breakup and fragmentation and to a good extent it has been simulated by CFD. Very few efforts have been made to simulate triggering, fine fragmentation and subsequent pressurization due to fast release of energy after fine fragmentation. Overall, it can be said that, there is tremendous scope for CFD simulation of severe accident phenomena in future.

Modelling of core melt scenarios in nuclear reactors

561

References [1] B.R. Sehgal, Nuclear Safety in Light Water Reactors: Severe Accident Phenomenology, 1st ed., Academic Press, 2011. [2] F. Martin-Fuertes, R. Barbero, J.M. Martin-Valdepenas, M.A. Jimenez, Analysis of source term aspects in the experiment Phebus FPT1 with the MELCOR and CFX codes, Nucl. Eng. Des. 237 (5) (2007) 509–523. [3] L. Zhang, Y. Zhou, S. Luo, Y. Zhang, G. Su, Z. Ma, L. Pan, Large eddy simulation for the thermal behavior of one-layer and two-layer corium pool configurations in HPR1000 reactor, Appl. Therm. Eng. 145 (2018) 38–47. [4] A. Horvat, B. Mavko, Numerical investigation of natural convection heat transfer in volumetrically heated spherical segments, Proceedings of the International Conference Nuclear Energy for New Europe, vol. 37, (2004) INIS-SI–06-001. [5] F.J. Asfia, V.K. Dhir, An experimental study of natural convection in a volumetrically heated spherical pool bounded on top with a rigid wall, Nucl. Eng. Des. 163 (1996) 333–348. [6] M. Buck, M. Burger, A. Miassoedov, X. Gaus-Liu, A. Palagin, L. Godin-Jacqmin, C. T. Tran, V. Chudanov, The LIVE program – results of test L1 and joint analyses on transient molten pool thermal hydraulics, Prog. Nucl. Energy 52 (2010) 46–60. [7] C.-T. Tran, N. Dinh, Simulation of core melt pool formation in a Reactor pressure vessel lower head using an Effective Convectivity Model, Nucl. Eng. Technol. 41 (2009) 929–944. https://doi.org/10.5516/NET.2009.41.7.929. [8] X. Zhang, et al., CFD simulation on critical heat flux of flow boiling in IVR-ERVC of a nuclear reactor, Nucl. Eng. Des. 304 (2016) 70–79. [9] T.Y. Chu, M.M. Pilch, J.H. Bentz, J.S. Ludwigsen, W.-Y. Lu, L.L. Humphries, Lower head failure experiments and analyses, (1998)NUREG/CR-5582, SAND98-2047. [10] C. Caroli, A. Miliozzi, F. Milillo, Thermo-mechanical behaviour of the reactor pressure vessel and corium molten pool in a severe accident with core melt down, SMiRT 15, Seoul Korea (Paper P02/2), 1999. [11] H.G. Willsch€utz, E. Altstadt, B.R. Sehgal, F.P. Weiss, Coupled thermal structural analysis of LWR vessel creep failure experiments, Nucl. Eng. Des. 208 (2001) 265–282. [12] B.R. Sehgal, Assessment of reactor vessel integrity (ARVI), Nucl. Eng. Des. 235 (2005) 213–232. [13] S. Thakre, L. Manickam, W. Ma, A numerical simulation of jet breakup in melt coolant interactions, Ann. Nucl. Energy 80 (2015) 467–475. [14] T.N. Dinh, Experimental and analytical studies of melt jet-coolant interactions: a synthesis, Nucl. Eng. Des. 189 (1999) 299–327. [15] T.N. Dinh, Investigation of film boiling thermal hydraulics under FCI conditions: results of analyses and a numerical study, Nucl. Eng. Des. 189 (1999) 251–272. [16] Y. Zhou, et al., Numerical simulation of metal jet breakup, cooling and solidification in water, Int. J. Heat Mass Transf. 109 (2017) 1100–1109. [17] S. Park, 3-D simulation of plunging jet penetration into a denser liquid pool by the RD-MPS method, Nucl. Eng. Des. 299 (2016) 154–162. [18] L. Cizelj, B. Koncˇar, M. Leskovar, Vulnerability of a partially flooded PWR reactor cavity to a steam explosion, Nucl. Eng. Des. 236 (2006) 1617–1627. [19] S. Thakre, W. Ma, L. Li, A numerical analysis on hydrodynamic deformation of molten droplets in a water pool, Ann. Nucl. Energy 53 (2013) 228–237.

562

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

[20] G. Berthoud, Models and validation of particulate debris coolability with the code MC3DREPO, Nucl. Eng. Des. 236 (2006) 2135–2143. [21] G. Berthoud, M. Valette, interaction, Development of a multidimensional model for the premixing phase of a fuelcoolant, Nucl. Eng. Des. 149 (1) (1994) 409–418. [22] M. B€urger, M. Buck, W. Schmidt, W. Widmann, Validation and application of the WABE code: investigations of constitutive laws and 2D effects on debris coolability, Nucl. Eng. Des. 236 (2006) 2164–2188. [23] SCDAP/RELAP5-3D Code Manual Volume 2: Modeling of Reactor Core and Vessel Behavior During Severe Accidents, 2003, NEEL/EXT-02-00589. [24] L. Bolshov, V. Chudanov, CFD approach to modeling of Core-Concrete Interaction, SMiRT 19, Toronto, August 2007, Paper # B07/1. [25] P. Chai, M. Kondo, N. Erkan, K. Okamoto, Numerical simulation of 2D ablation profile in CCI-2 experiment by moving particle semi-implicit method, Nucl. Eng. Des. 301 (2016) 15–23 ISSN 0029-5493. [26] X. Li, A. Yamaji, Three-dimensional numerical study on the mechanism of anisotropic MCCI by improved MPS method, Nucl. Eng. Des. 314 (2017) 207–216. [27] C.R.B. Lister, On the penetration of water into hot rock, Geophys. J. R. Astron. Soc. 39 (1974) 465–509. [28] M. Epstein, Dryout heat flux during penetration of water into solidifying rock, J. Heat Transf. 128 (2006) 847–850. [29] D. Paladino, S.A. Theerthan, B.R. Sehgal, DECOBI: investigation of melt coolability with bottom coolant injection, Prog. Nucl. Energy 40 (2002) 161–206. [30] W. Widmann, M. B€urger, G. Lohnert, H. Alsmeyer, W. Tromm, Experimental and theoretical investigations on the COMET concept for ex-vessel core melt retention, Nucl. Eng. Des. 236 (2006) 2304–2327. [31] IAEA Tecdoc 1594, Analysis of Severe Accidents in Pressurized Heavy Water Reactors, International Atomic Energy Agency (IAEA), (2008). [32] IAEA Tecdoc 1727, Benchmarking Severe Accident Computer Codes for Heavy Water Reactor Applications, International Atomic Energy Agency (IAEA), (2013). [33] P.P. Kulkarni, S.V. Prasad, A.K. Nayak, P.K. Vijayan, Thermal and structural analysis of calandria vessel of a PHWR during a severe accident, Nucl. Eng. Technol. 45 (4) (2013) 463–476. [34] D. Dupleac, Analysis of late phase severe accident phenomena in CANDU plant, ASME J. Nuclear Rad. Sci. 3 (2017) 020904. [35] S.V. Prasad, A.K. Nayak, Experimental investigation of heat transfer during severe accident of a pressurized heavy water reactor with simulated decay heat generation in molten pool inside calandria vessel, Nucl. Eng. Des. 303 (2016) 75–87. [36] J.T. Rogers, Thermal and hydraulic behavior of CANDU cores under severe accident conditions, Final report in: Analytical Methods and Results, vol. 1, AECL, Ontario, 1984. [37] J.R. Welty, C.E. Wicks, R.E. Wilson, G.L. Rorrer, Fundamentals of Momentum Heat and Mass Transfer, fifth ed., John Wiley and Sons, Hoboken, NJ, 2008. [38] H.D. Baehr, K. Stephan, Fundamentals of Momentum, Heat and Mass Transfer, second ed., Library of Congress, Stuttgart, 2008. [39] V.R. Voller, Modeling solidification processes, Mathematical Modeling of Metals Processing Operations Conference, American Metallurgical Society, Palm Desert, CA, 1987 Technical report. [40] D. Dupleac, Analysis of late phase severe accident phenomena in CANDU plant, J. Nucl. Eng. Radiat. Sci. 3 (2017) 1 020904.

Modelling of core melt scenarios in nuclear reactors

563

[41] S. Rahman, M. Buerger, M. Buck, G. Pohlner, R. Kulenovic, A.K. Nayak, B.R. Sehgal, Analysis of dryout behaviour in laterally non-homogeneous debris beds using the MEWA-2D code, 13th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-13) Kanazawa City, Ishikawa Prefecture, Japan, September 27– October 2, 2009. [42] P. Sch€afer, M. Groll, R. Kulenovic, Basic investigations on debris cooling, Nucl. Eng. Des. 236 (2006) 2104–2116. [43] S. Ergun, Fluid flow through packed column, Chem. Eng. Prog. 48 (2) (1952) 89–94. [44] P.P. Kulkarni, M. Rashid, R. Kulenovic, A.K. Nayak, Experimental investigation of coolability behaviour of irregularly shaped particulate debris bed, Nucl. Eng. Des. 240 (2010) 3067–3077. [45] A.K. Nayak, P.P. Kulkarni, N. Kumar, R. Kulenovic, P. Schaefer, P. Shroeder, Studies on coolability behaviour of particulate debris beds, (2009)DST-DAAD project report INT/ DAAD/P-158/2007. [46] D. Magallon, I. Huhtiniemi, H. Hohmann, Lessons learnt from FARO-TERMOS corium melt, Nucl. Eng. Des. 189 (1999) 223–238. [47] I. Huhtiniemi, D. Magallon, H. Hohmann, Results of recent KROTOS FCI tests: alumina vs corium, Nucl. Eng. Des. 189 (1999) 379–389. [48] B.W. Spencer, K. Wang, C.A. Blomquist, L.M. McUmber, J.P. Schneider, Fragmentation and quench behaviour of corium melt streams in water, (1994)NUREG/CR-6133, ANL93/32. [49] H. Nagasaka, M. Kato, I. Sakaki, Y. Cherepnin, Y. Vasilyev, A. Kolodeshnikov, V. Zhdanov, V. Zuev, COTELS project (1): Overview of project to study FCI and MCCI during a severe accident, OECD Workshop on Ex-Vessel Debris Coolability, Karlsruhe, Germany, November 15–18, 1999. [50] H. Hardee, R. Nilson, Natural convection in porous media with heat generation, Nucl. Sci. Eng. 63 (1977) 119–132. [51] V. Dhir, I. Catton, Study of dryout heat fluxes in beds of inductively heated particles, (1977)Technical report, NUREG CR-0262. [52] I. Catton, V. Dhir, C. Somerton, An experimental study of debris coolability under pool boiling conditions, (1983)Technical report NP-3094, EPRI. [53] A. Reed, K. Boldt, R. Gorham-Bergeron, R. Lipinski, T. Schmidt, DCC-1 and DCC-2 degraded core coolability analysis, (1985)Technical report CR-439 NUREG SAND851967. [54] D.H. Cho, L. Bova, Formation of dry pockets during water penetration into a hot particle bed, Trans. Am. Nucl. Soc. 41 (1983) 418. [55] V.E. Schrock, C.H. Wang, S. Revankar, L.H. Wei, S.Y. Lee, D. Squarer, Flooding in particulate beds and their role in dryout heat flux, Proceedings of the Sixth Information Exchange Meeting on Debris Coolability, UCLA, USA, 1984. [56] K. Hu, T.G. Theofanous, Scale effects and structure of dryout zone in debris bed coolability, Proceedings of Sixth Information Exchange Meeting on Debris Coolability, UCLA, USA, 1984. [57] V.X. Tung, V.K. Dhir, Quenching of debris beds having variable permeability in the axial and radial directions, Nucl. Eng. Des. 99 (1987) 275–284. [58] V.X. Tung, V.K. Dhir, D. Squarer, Quenching by top flooding of a heat generating particulate bed with gas injection at the bottom, Proceedings of the Sixth Information Exchange Meeting on Debris Coolability, EPRI Report NP-4455, 1986. [59] B.R. Sehgal, M.J. Konovalikhin, Z.L. Yang, I.V. Kazachkov, M. Amjad, G.J. Li, Experimental investigations on porous media coolability, (2001)RIT/NPS research report.

564

Advances of Computational Fluid Dynamics in Nuclear Reactor Design and Safety Assessment

[60] A.K. Nayak, B.R. Sehgal, A.V. Stepanyan, An experimental study on quenching of a radially stratified heated porous bed, Nucl. Eng. Des. 236 (19) (2006) 2189–2198. [61] H. Bj€ornsson, S. Bj€ornsson, Th. Sigurgeirsson, Penetration of water into hot rock boundaries of magma at Grı´msv€otn, Nature 295 (1982) 580–581. https://doi.org/10. 1038/295580a0. [62] E.A. Jagla, A.G. Rojo, Sequential fragmentation: the origin of columnar quasihexagonal patterns, Phys. Rev. E 65 (2002) 026203. [63] B.R. Sehgal, B.W. Spencer, ACE program phase D: melt attack and coolability experiments (MACE) program, in: Second OECD (NEA) CSNI Specialist Meeting on Molten Core Debris-Concrete Interactions, 1992 25/11. [64] M.T. Farmer, B.W. Spencer, D.J. Kilsdonk, R.W. Aeschlimann, Results of MACE corium coolability experiments M0 and M1b, in: Proc. ICONE-8 Conference on Nuclear Engineering, Baltimore, MD, April 2–6, 2000. [65] M.T. Farmer, B.W. Spencer, A review of database pertaining to ex-vessel debris coolability, ACEX-TR-C34, (2001). [66] R.E. Blose, J.E. Gronager, A.J. Suo-Anttila, J.E. Brockmann, SWISS: Sustained heated metallic melt-concrete interactions with overlaying water pools, NUREG/CR-4727, (1987)SAND85-1546 R3, R4, R7. [67] R.E. Blose, D.A. Powers, E.R. Copus, J.E. Brockmann, R.B. Simpson, D.A. Lucero, Coreconcrete water interactions with overlaying water pools – the WETCOR 1 test, NUREG/ CR-5907, (1993). [68] D. Magallon, H. Hohmann, High pressure corium melt quenching tests In FARO, Nucl. Eng. Des. 155 (1995) 253–270. [69] T.G. Theofanous, C. Liu, W.W. Yuen, Coolability and quench of corium-concrete interactions by top flooding, MACE-TR-D15, (1998). [70] H. Alsemeyer, B. Spencer, W. Tromm, The COMET concept for cooling of ex-vessel corium melts, ICONE-6, San Diego, USA, 1998. [71] M.T. Farmer, S. Lomperski, S. Basu, Results of reactor material experiments investigating 2-D core-concrete interaction and debris coolability, in: Proceedings of ICAPP-04, June 13–17, Pittsburg, USA, 2004. [72] S. Lomperski, M.T. Farmer, Experimental evaluation of the water ingression mechanism for corium cooling, Nucl. Eng. Des. 237 (2007) 905–917. [73] A.K. Nayak, R.K. Singh, P.P. Kulkarni, B.R. Sehgal, A numerical and experimental study of water ingression phenomena in melt pool coolability, Nucl. Eng. Des. 239 (2009) 1285–1293. [74] P.J. Berenson, Film boiling heat transfer from a horizontal surface, J. Heat Transf. 83 (1961) 351–357. [75] W.H. Rohsenow, A method of correlating heat transfer data for surface boiling of liquids, Trans. ASME 74 (1952) 969–976. [76] R. Sinha, A.K. Nayak, B.R. Sehgal, Modeling the natural convection heat transfer and dryout heat flux in a porous debris bed, ASME J. Heat Transf. 130 (10) (2008) 104503. [77] R.D. Ritchie, J.F. Knott, J.R. Rice, On the relationship between critical tensile stress and fracture toughness in mild steel, J. Mech. Phys. Solids 21 (1973) 395–410. [78] F.M. Beremin, A local criterion for cleavage fracture of a nuclear pressure vessel steel, Metall. Trans. A. 14A (1983) 2277–2287. [79] A.K. Nayak, B.R. Sehgal, R.S. Rao, Study on water ingression phenomenon in melt coolability, The 11th International Topical Meeting on Nuclear Reactor ThermalHydraulics (NURETH-11) Paper: 449, Avignon, France, 2005.

Modelling of core melt scenarios in nuclear reactors

565

[80] S. Lomperski, M.T. Farmer, Performance testing of engineered corium cooling systems, Nucl. Eng. Des. 243 (2012) 311–320. [81] N. Zuber, On stability of boiling heat transfer, ASME Trans. 80 (1958) 711–720.

Further reading [82] H. Hohmann, D. Magallon, H. Schins, A. Yerkess, FCI experiments in the aluminum oxide/water system, Nucl. Eng. Des. 155 (1995) 391–403. [83] G. Repetto, T. Garcin, M. Rashid, R. Kulenovic, W. Ma, L. Li, Investigation of multidimensional effects during debris cooling, 5th European Review Meeting on Severe Accident Research (ERMSAR-2012) Cologne, Germany, March 21–23, 2012. [84] P.P. Kulkarni, A.K. Nayak, Study on coolability of melt pool with different strategies, Nucl. Eng. Des. 270 (2014) 379–388. [85] N. Singh, P.P. Kulkarni, A.K. Nayak, Experimental investigation on melt coolability under bottom flooding with and without decay heat simulation, Nucl. Eng. Des. 285 (2015) 48–57. [86] S.V. Prasad, A.K. Nayak, Experimental Study on Melt Coolability Capability of Calandria Vault Water During Severe Accident in Indian PHWRs for Prolonged Duration, ASME J of Nuclear Rad Sci 4 (3) (2018) 031014.