Modeling and verification of three-dimensional simulation for BWR in-vessel core degradation

Annals of Nuclear Energy 101 (2017) 182–195

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Modeling and verification of three-dimensional simulation for BWR in-vessel core degradation Tsuyoshi Okawa ⇑, Tetsuo Nakajima Regulatory Standard and Research Department, Secretariat of Nuclear Regulation Authority (S/NRA/R), 1-9-9 Roppongi-First Bldg., Roppongi, Minato-ku, Tokyo 106-8450, Japan

a r t i c l e

i n f o

Article history: Received 7 June 2016 Received in revised form 29 August 2016 Accepted 25 September 2016

Keywords: BWR core degradation Severe accident Three-dimensional simulation Tight coupling calculation

a b s t r a c t A Multifunction Model has been developed for describing BWR in-vessel core degradation from initial cladding temperature rise to eventual molten materials accumulation in the lower part of reactor vessel as a fundamental safety research by Regulatory Standard and Research Department, Secretariat of Nuclear Regulation Authority (S/NRA/R). The Multifunction Model aims to reduce the uncertainties of severe accident assessment related to physical models on core damage progression, melt formation and material relocation processes in the reactor vessel. The Multifunction Model consists of a thermalhydraulic module, a fuel pin behavior module and a neutronic module with three-dimensional geometry. For reducing the uncertainties, the thermal-hydraulic module incorporated multi-phase, multi-component and multi-velocity field with multiple chemical and eutectic reactions. Four fuel failure models were incorporated into the fuel pin behavior module in order to simulate the failures of cladding rupture, brittle fracture, cladding melting and pellet dissolution. In order to decrease the calculation time, the data table of macroscopic cross-section was prepared for neutronic diffusion calculation of three fuel conditions in advance. The verification results revealed that the Multifunction Model demonstrated a reasonable capability to simulate the BWR in-vessel core degradations. The validation will be carried out for existing severe accident experiments. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction A large amount of fuel debris would be generated during coremelt/degradations in the severe accidents at Fukushima Daiichi Units 1–3. The fuel debris location and volume are still unclear due to high radiation condition in these reactors. It is difficult to assess these degraded in-vessel conditions precisely, since current severe accident codes have a certain amount of uncertainty to estimate core degradation processes. The degraded core processes are a key factor in the progression of severe accidents, since the processes provide the initial conditions for in- and ex-vessel phenomena. Regulatory Standard and Research Department, Secretariat of Nuclear Regulation Authority (S/NRA/R) has investigated the Fukushima accident progression. The S/NRA/R started development of a time-dependent three-dimensional detailed model, ‘‘Multifunction Model” to reduce the uncertainties of severe accident assessment related to physical models of in-vessel core degradation of Boiling Water Reactor (BWR) as one of the fundamental ⇑ Corresponding author. E-mail addresses: [email protected] (T. Okawa), tetsuo_nakajima01@ nsr.go.jp (T. Nakajima). http://dx.doi.org/10.1016/j.anucene.2016.09.042 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.

safety research activities (Okawa and Nakajima, 2016). The Multifunction Model simulates behaviors of fuel and structural materials under severe accident conditions as well as normal operational and transient conditions. The current severe accident codes (Sehgal, 2012) can be divided into two groups: (i) integral code group: MAAP (FAUSKE, 2016) and MELCOR (U.S. Nuclear Regulatory Commission, 2015) developed in the United States, and ASTEC (Chatelard et al., 2014) developed in France and Germany, and (ii) detailed code group: SCDAP/RELAP5 (U.S. Nuclear Regulatory Commission, 2001) developed in the United States, ICARE/CATHARE (Seiler et al., 2008) developed in France and IMPACT/SAMPSON (Satoh et al., 2000) developed in Japan. The Multifunction Model is classified as the detailed code. The Multifunction Model consists of a thermal-hydraulic module, fuel pin behavior module and neutronic module with a control junction. Under an agreement made with Japan Atomic Energy Agency (JAEA), the neutronic and the fuel pin behavior modules applied the existing codes: SKETCH (JAEA, 2013) and FEMAXI-6 (Suzuki and Saito, 2006) respectively which are well verified and validated by the JAEA. However the thermal-hydraulic module was newly developed, and therefore, the thermal-hydraulic module has been verified carefully over two years by the S/NRA/R. For reducing the uncertainties, the thermal-hydraulic module incorporated

183

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

multiple melting points, multiple eutectic reaction growth rates, an improved candling model and a complicated structural model. The fuel pin behavior module added four failure models of cladding rupture, brittle fracture, cladding melting and pellet dissolution. In order to decrease the calculation time of the neutronic module, the data table of macroscopic cross-section was considered for three fuel conditions in advance. The verification results revealed that the Multifunction Model demonstrated a reasonable capability to simulate the BWR in-vessel core degradations. The validation will be carried out for using existing severe accident experiments. The present paper describes development of the prototype of the Multifunction Model, the findings through the verifications for each calculation module and tight coupling calculation results with all of the modules. 2. Uncertainty in severe accident simulation for BWR in-vessel core degradation There are a multitude of uncertainties in simulating of BWR invessel core degradations with current severe accident codes. The tiny differences (or uncertainties) between calculations and facts in the initial condition lead to certain changes in the course of subsequent severe accident progressions. Hence the S/NRA/R identified that the key issues of the current severe accident codes for in-vessel core degradation assessment are uncertainties for physical models: (i) multi-eutectic reaction, (ii) multi-phase change and (iii) mass and energy transfer in the vessel and (iv) fine structural geometry, particularly complicated BWR lower structures (tie plate, fuel support piece, core support plate and control rod guide tube) under the core zone. The concerns should be resolved for enhancing simulation of severe accident. 3. Engineering challenge for severe accident model development The purpose of the Multifunction Model is to reduce and decrease the foregoing uncertainties of the physical models described in Section 2, as an engineering challenge. For reducing uncertainties of the physical models, it is required for the thermal-hydraulic module to incorporate multi-phase, multicomponent and multi-velocity field with multiple chemical and eutectic reactions. Furthermore, in order to simulate BWR invessel core degradations in the normal operational conditions (setting initial conditions) and severe accident conditions, the Multifunction Model should include neutronic diffusion and fuel pin behavior analyses. Additionally the fine three-dimensional geometry is required for the complicated BWR lower structures. The Multifunction Model should simulate the following core degradation progress: (i) cladding temperature rise, oxidation and embrittlement due to lower water level, (ii) collapse of control rods, (iii) melting of fuel pins and assemblies, (iv) molten materials flow to the lower part of fuel assemblies and formation of debris beds, (v) molten materials flow to the core support plate and rupture of the plate and (vi) molten materials flow along the lower part of reactor vessel. 4. Modeling and verification of thermal-hydraulic, fuel pin behavior and neutronic modules The present work is an attempt to model and verify the thermal-hydraulic, fuel pin behavior and neutronic modules. All of the modules must reduce principally the uncertainties mentioned in Section 2. Each of the modules will be presented in Sections 4.1,4.2 and 4.3. The Multifunction Model programming languages are FORTRAN. The model size is about 200,000 lines, dis-

tributed in 1700 subroutines. The Multifunction Model target is grid computers with Linux operating systems. 4.1. Thermal-hydraulic module The fundamental idea of the thermal-hydraulics module (3-dimensional, 3-velocity field, 3-phase and multi-component) developed by the S/NRA/R is based on the models of the SIMMER code (Bohl and Luck, 1990) and other multiphase ideas (Ishii and Hibiki, 2011). The convective components are liquid (UO2, Zr, B4C, Fe, eutectic, coolant, etc.), solid particle/chunk (UO2, Zr, B4C, Fe, eutectic, etc.) and gas /vapor (UO2, SUS, steam, fission product, etc.). The total number of the liquid and gas components is twenty and ten, respectively (see Table 1). The ten eutectic components are introduced to the module to evaluate more realistic core degradation. The six structural components are assumed to evaluate the complicated core structural effect for simulating the molten materials relocation as shown in Table 2. The velocity fields are currently three: q1:Coolant, q2:Corium, q3: Gas/steam mixed gas to reduce the calculation time. Almost all of the core components are already prepared in the data base of the thermal-hydraulic module. The several phase change models are assumed for simulating material relocation, practically achievable. After fuel pin failures, the data of liquid/solid fuel pellet and liquid cladding are transferred from the fuel pin behavior module to the thermalhydraulic module. The liquid/solid fuel and fuel particle on the fuel cladding outer surface from the fuel pin behavior module are conserved in the thermal-hydraulic module. The data of heat flux into coolant are transferred from the thermal-hydraulic module to the fuel pin behavior module. The data of cladding surface tempera-

Table 1 Liquid (with Eutectic) and Gas Components. Number

Components

Materials

Velocity field

Le1 Le2 Le3 Le4 Le5 Le6 Le7 Le8 Le9

Liquid fuel Particle fuel Chunk fuel Liquid Zr Particle Zr Liquid steel Particle steel Liquid control rod Particle control rod

1) UO2 2) U3O8

q2: Corium velocity field

Le10

Coolant (Water)

H2O

q1: Coolant velocity field

Le11 Le12 Le13 Le14

Liquid B4C-Fe eutectic Particle B4C-Fe eutectic Liquid UO2-Zr eutectic Particle UO2-Zr eutectic

q2: Corium velocity field

Le15 Le16

Liquid B4C-Fe-Zr eutectic Particle B4C-Fe-Zr eutectic

Le17

Le19 Le20

Liquid B4C-Fe-Zr-UO2 eutectic Particle B4C-Fe-Zr-UO2 eutectic Liquid Zr-Inconel eutectic Particle Zr-Inconel eutectic

1) 2) 1) 2) 3) 1) 2) 3) 1) 2) 3) 4) 1) 2)

Ge1 Ge2 Ge3 Ge4 Ge5 Ge7 Ge8 Ge9 Ge10

Fuel vapor Steel vapor Coolant vapor Fission gas Hydrogen gas Nitrogen gas Oxygen gas Argon gas Helium gas

UO2 SUS H2O Xe H2 N2 O2 Ar He

Le18

1) Zr 2) ZrO2 1) SUS 2) Fe2O3 1) Ag-InCd 2) B4C

Fe B4C U Zr O Zr Fe B4C Zr Fe B4C UO2 Zr Inconel

q3: Gas velocity field

184

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

Table 2 Structural components. Number

Name of components

Application structures

Materials

S1(1) S1(2) S2

Cylindrical structure surface Cylindrical structure surface Steel wall structure

– – 1) SUS 2) Fe2O3

S3

Zr wall structure

Fuel cladding surface Control rod surface Control rod blade Core support plate Penetration parts of reactor pressure vessel (RPV) bottom Upper part of fuel assembly Channel box

S4(1) S4(2) S4(3)

Crust wall structure Cylindrical crust structure Cylindrical crust structure

Crust on S2 and S3 structure Crust on S1(1) structure Crust on S1(2) structure Control rod model Penetrations of reactor pressure vessel (RPV) bottom

S5

Water rod structure

Water rod

S6

Penetrations of reactor pressure vessel (RPV) bottom-detail model

Penetrations of reactor pressure vessel (RPV) bottom

tures are transferred from the fuel pin behavior module to the thermal-hydraulic module. A control rod is modeled by annulus heat interactive structure in the thermal-hydraulic module. The control rod changes to liquid/solid particle when its temperature reaches the eutectic temperature, 1523 K (which is much lower than the melting point of B4C, 2723 K). The thermal-hydraulic module calculates both the inside and outside of the control rod behavior (mainly eutectic reaction). The conservation equations of mass, momentum and energy for liquid phase can be written as the following equations, respectively: Mass conservation:


m @q
m v q Þ ¼ Cm þ r  ðq @t

m ¼ 1; . . . ; NLd ;

ð1Þ

where N Ld ¼ 20 is the number of liquid-field density components subscribed bym,
m is the macroscopic density of liquid-field density component q m, v q is the liquid-field velocity vector, Cm is the total mass transfer rate per unit volume from liquidfield density component m. The first term of the left-hand side is time-dependent term, and the second one represents the convection term. The right-hand side is the source term of liquid phase. Momentum conservation:

X
qvq @q
q v q v q Þ þ aq rp  q
q g þ K qS v q  K qq0 ðv q0  v q Þ  VM q þ r  ðq @t q0 X ¼  ðCqq0 v q  Cq0 q v q0 Þ  CqS v q  CSq v q ; ð2Þ q0

where
q is the summation of liquid-field densities which have the q liquid-field velocity vq, aq is the summation of liquid-field volume fractions which have the liquid-field velocity vq, p is the local pressure, g is the acceleration of gravity, K qq0 is the momentum exchange function between the velocity fields vq and vq0 , K qS is the momentum exchange function between the velocity fields vq and the structure field,

1) Zr 2) ZrO2 1) UO2 2) U3O8 3) Zr 4) ZrO2 5) SUS 6) Fe2O3 7) Ag-In-Cd 8) B4C 1) Zr 2) ZrO2 Inner SUS liner Welding materials

Cqq0 is the mass transfer rate per unit volume from the velocity field vq to the velocity field vq0 , CqS is the mass transfer rate per unit volume from the velocity field vq to the structure field, VM q is the virtual mass force of the velocity field vq. The first and the second terms of the left-hand side are timedependent and convection terms, respectively. The third term of the left-hand side represents the pressure gradient term. The right-hand side is the source term of liquid phase momentum. Energy conservation:

 
M eM @q @ aM
M eM v q Þ þ p þ r  ðq þ r  ðaM v q Þ @t @t ¼ Q IH þ Q N þ Q MF þ Q VC þ Q HT

M ¼ 1; . . . ; NLe ;

ð3Þ

where N Le ¼ 20 is the number of liquid-field energy components subscribed by M, eM is the specific energy of liquid-field energy component M, aM is the volume fraction of liquid-field energy component M, Q IH is the energy dissipation rate per unit volume due to friction between velocity fields, Q N is the energy generation rate per unit volume due to decay heat, Q MF is the energy transfer rate per unit volume due to Melting/ Freezing, Q VC is the energy transfer rate per unit volume due to Vaporization/Condensation, Q HT is the energy transfer rate per unit volume due to heat transfer from fuel pins, structure fields, etc. The first term of the left-hand side is the time-dependent term, and the second one represents the convection term of liquid phase energy. The third term of the left-hand side is pressure dependent. The right-hand side represents the source term of liquid phase energy. In the module, the number of liquid-field energy components NLe is twenty as shown in Table 1. The first velocity field q1 corresponds to coolant, and all of the components other than coolant belong to the second velocity field q2. Similarly, the conservation equations of mass, momentum and energy for gas phase can be

185

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

written. The number of gas-field energy components is ten. All of the gas components have the third velocity field q3. 4.1.1. Important physical models for BWR In-Vessel core degradation After stability limits of materials in BWR core structures are reached, melting and liquefaction occurs and the core degradation proceeds through core melt and relocation. To enhance the calculation model, the following three effects: a) melting point and eutectic reaction zone growth rate, b) candling effect and c) complicated structures under core zone must be incorporated into the Multifunction Model from the knowledge of CORA experiment and inherent characteristics of BWR in-vessel structures. 4.1.1.1. Melting point and eutectic reaction zone growth rate. During BWR in-vessel core degradations, it is highly significant to consider various melting points and eutectic reactions for several materials of the BWR in-vessel components based on the CORA experimental knowledge (Hering and Hofmann, 2007). Hence, a MATPRO material database (A Library of Materials Properties for Light-WaterReactor Accident Analysis) (INEEL, 2003) was introduced for key single structural core materials (not compound) to simulate the detail solidus/liquidus temperatures, for instance, UO2, ZrO2, B4C, Zr and SUS. Compound materials reactions are divided into two groups: for intact core geometries and for un-intact core geometries. One of the key phenomena in the BWR severe accident is eutectic reaction in the intact core geometry before molten materials relocation, for example: UO2 fuel-Zr cladding interaction, B4C neutron absorberSUS cladding interaction and Zr-Inconel grid spacer spring interaction. Hofmann, FZK (Hofmann, 1999) gave the physic-chemical material behaviors in the BWR core with increasing temperature and proposed the reaction zone growth rates in Arrhenius diagram against reciprocal temperature for SUS-Zr, UO2-Zr, B4C-Zr, B4C-SUS, Zr-Inconel interactions etc. With intact core geometries, the Hofmann correlations therefore were introduced to the module. After eutectic reactions in the intact core geometries, chemical reactions start locally between liquid molten eutectic materials and intact structures. With these complicated reactions, the Hofmann correlations are applied to three occasions as follows:  Molten B4C-stainless steel eutectic attack to the Zr channel box,  Molten B4C-stainless steel-Zr eutectic attack to the Zr fuel pin,  Molten B4C-stainless steel-Zr-UO2 eutectic attack to the core structures. The concept of eutectic reaction zone growth rate is applied to the foregoing un-intact structures as shown in Fig. 1. The liquid molten material starts to travel from above on the wall structure. The boundary between the liquid molten material and wall structure is heated up by the decay heat or the chemical reaction of the liquid molten materials. The time-depend eutectic reaction zone parabolic growth rate is defined as,

 1=2 nðt þ DtÞ ¼ nðtÞ2 þ K  Dt ;   Q ; K ¼ K 0 exp  RT

ð4Þ ð5Þ

where n(t) is the thickness of contact surface between the liquid molten material and the wall structure, Dt is the time step used for the calculation, K is the reaction zone growth rate in m2/s, K0 is the constant for the reaction zone growth rate in m2/s, Q is the activation energy in J/mol, R is the gas constant which is 8.314 in J/mol/K, T is the temperature in K. A Binary Contact Area (BCA) is analyzed by the flow regime subroutine. Introducing the Dn into

the following equation, the volume and thickness of eutectic reaction zone in the wall structure is obtained:

Dv ¼ BCA  Dn;

ð6Þ

where Dv is the volume of eutectic reaction zone in m , BCA is the binary contact area in m2, S is the surface area of the wall structure in the calculation mesh in m2, Dx is the eutectic width in m, 3

Dx ¼

BCA  Dn: S

ð7Þ

When the eutectic width expands to a thin wall layer, the equivalent eutectic layer becomes an eutectic particle and the particle falls down to the next calculation mesh. Quantitative analysis of the eutectic reactions however cannot be estimated exactly in the foregoing zones, since there is no database for these complicated molten material combinations. Hence in calculation, heaviest materials will represent these complicated molten materials in a calculation cell. 4.1.1.2. Candling effect. Based on the CORA experimental experience, Hofmann and Hagen et al. (Hagen et al., 2015) suggested that during initial material relocation, it is essential to simulate relocation of molten materials by droplet/rivulets, what is known as ‘‘Candling” as shown in Fig.2(a). The previous Multifunction Model (see Fig. 2(b)) was capable of simulating the dispersed flow (metal/ water particle). This is however different description from the actual phenomenon of Fig. 2(a). A bubbly flow for core degradations is applied to the film flow effect of ‘‘Candling,” as shown in Fig.2(c). By applying the bubbly flow regime (void fraction: ag 6 0.3) to molten material instead of the dispersed flow regime (void fraction: ag P 0.7), ‘‘Candling” phenomenon including heat transfer and momentum exchange can be simulated adequately because molten material interacts directly with structure components of fuels and control rods, etc. From the CORA experimental results, it is very important to understand the local molten materials pool/blockage effect (Veshchunov and Palagin, 1998; Veshchunov and Shestak, 2008; Veshchunov et al., 2013) at the grid spacers in relocation of initial eutectic/melting materials during ‘‘Candling”. Hence very thin three grid spacers are arranged in an axial direction as a major impediment to the downward movement of melts. Additionally, Zr-Inconel grid spacer spring interaction must be also considered since Inconel produces a large amount of eutectic material with Zr. 4.1.1.3. Complicated structures under core zone. There would be huge uncertainties for molten materials relocation in the BWR lower head, since the BWR structures (tie plate, fuel support piece, core support plate, control rod guide tube) under the core zone are quite complicated geometry. It is desirable to model the complicated structures three dimensionally. The fine three-dimensional geometry was made for the complicated zone and also the molten materials flow in the following four paths as shown in Fig. 3 (Gauntt and Humphries, 1997) to simulate realistically the material traveling: (i) main flow path of the molten control rods (which flow out first), (ii) main flow path of the molten fuel, (iii) sub-flow path of the molten fuel to the control rod guide tube via the blade channel, (iv) sub-flow path of the molten control rods to the fuel bundle. 4.1.2. Verification calculation for high temperature molten and particle fuel into saturated Water The purpose of preliminary verification was to observe phase changes of the components during cooling phase in a typical pressure vessel. Fig. 4 shows the initial calculation condition with fine meshes (top) and the calculation results (bottom) for verification calculation. The top figure shows that the highly-heated liquid and particle fuels are dropped into the saturated water of 559 K.

186

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

Fig. 1. Traveling liquid molten materials downward along with eutectic reaction.

The water was filled up to the core support plate. The temperatures of liquid and particle fuels were 3114 K (just above the melting point of UO2) and 3000 K, respectively. The volume of the fuels was assumed to be 35% equal to the volume fraction of fuel pellets in the single fuel assembly. The pressure of the calculation zone was 7 MPa (assuming normal operational condition). The three dimensional calculation model of single BWR fuel assembly equivalent model was from top to bottom of a typical BWR reactor pressure vessel. The axial and radial calculation dimensions were same as a Japanese single 9x9 BWR fuel assembly. The bottom figures revealed the time dependence of the transient volume fraction in each axial position in the reactor pressure vessel.

Immediately after the fuel reached to the surface of saturated water (at 1 s), the coolant boiled and the volume fraction of coolant decreased suddenly near the core support plate. Simultaneously, the liquid fuel (red) became to the fuel particle (orange) due to freezing by contact with the saturated water. From two to three seconds after the fuel drop, the particle fuel started to accumulate on the surface of bottom-end of reactor pressure vessel and the coolant started to boil at the same time under the core zone. Four seconds after the fuel drop, most water was disappeared by boiling due to drop of the hot fuel particle. At nine seconds, most of the coolant area was occupied by the fuel particle. The foregoing results were confirmed

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

Fuel rod, Control rod, Nuclear instrumentaon etc.

Water droplet Gas phase

Molten material (Candling)

(a) Actual “Candling” phenomenon Dispersed flow

Solid metal parcle Water droplet Gas phase

Fuel rod, Control rod, Nuclear instrumentaon etc.

Liquid metal droplet

(b) Dispersed flow model - Previous model Bubbly Dispersed flow flow

Water droplet Gas phase Molten material Fuel rod, Control rod, Nuclear instrumentaon etc.

187

Atomic Energy Agency (JAEA) is selected for the fuel pin behavior module. The FEMAXI-6 code is for analysis of thermal and mechanical behaviors by using Finite Element Method (FEM) for Light Water Reactor (LWR) fuel rods during steady-state and transient conditions. The coupled analysis makes it possible to evaluate thermodynamic, mechanical and material transfer behaviors at severe accidents like the Fukushima accident. The thermodynamic analysis subroutine calculates the temperature and gas pressure distributions in a one-dimensional axisymmetric problem in the radial direction. The fuel rod is divided into several segments in the axial direction considering the axial power distribution. The mechanical analysis subroutine calculates the interactions of pellets and claddings using the two-dimensional axisymmetric finite element method. The fuel rod region is divided into several segments in axial direction, same as the above thermodynamic analysis. The each segment is divided into several ring elements in the radial direction. Newton-Raphson method is applied to evaluate the stress-strain relationship of fuel and cladding mechanical behaviors. It is required to improve the FEMAXI-6 in severe accident conditions for example: (i) material properties at very high temperature, (ii) Zr oxidation model at high temperature, (iii) hydrogen production model, (iv) molten fuel relocation model, (v) Zr-UO2 eutectic reaction and melting model. Because of adding the foregoing points, the module is capable of describing the oxidation thickness of cladding outer surface (original Zr to ZrO2) and the weight of hydrogen generated by high temperature steam. Furthermore, the module simulates that the inner UO2 reacts with Zr and forms the UO2-Zr eutectic which expands to the outer ZrO2 cladding. The code evaluates the weight of molten eutectic and residual UO2. 4.2.1. Failure criteria of fuel pin The original FEMAXI-6 analysis however is capable of simulating up to transient conditions; hence the code must be improved for fuel failures. The models were incorporated into the FEMAXI6 in order to simulate the following cladding failures: cladding rupture, brittle fracture, cladding melting and pellet dissolution. 4.2.1.1. Cladding rupture. Cladding rupture occurs below 1470 K which is one of the design basis events due to over plastic deformation of the Zr cladding. The criterion of plastic deformation is performed using the three equations for different temperature conditions from the FRAPTRAN criteria (Geelhood et al., 2011), as follows: 940 < T < 1200 K

efail ¼ 1:587979  109 T 4  6:692798  106 T 3 þ 1:053049  102 T 2  7:331051T þ 1906:22;

ð8Þ

1200 < T < 1700 K

efail ¼ 1:67939  108 T 3 þ 6:23050  105 T 2  7:360497

(c) Candling model - New model Fig. 2. Candling phenomenon and calculation model.

 102 T þ 28:1199; 1700 K < T

efail ¼ 0:544589; not quantitatively but qualitatively by considering proven eventadvancements.

4.2. Fuel pin behavior module It is desirable to evaluate cladding and fuel deformation up to normal operations, since the fuel condition affects to severe accident progress. The proven FEMAXI-6 code developed by Japan

ð9Þ

ð10Þ

where efail is the plastic deformation at the rupture (m/m), T is the cladding temperature (K). 4.2.1.2. Brittle fracture. Cladding brittle failure occurs at quench phases after a high temperature condition. The failure is governed by the amount of oxygen concentration in b-Zr phase. After disappearance of the b-Zr phase, the oxidized cladding tube is no longer mechanically stable. The cladding is damaged in the two cases, (i) the cladding temperature exceeds 1973 K and then is quenched

188

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

1. Flow path of molten control rods (flow out first) 2. Flow path of molten fuel 3. Flow path of molten fuel to control rod guide tube channel 4. Flow path of molten control rods to fuel bundle

Fig. 3. Complicated structure under core zone and four paths of molten materials relocation.

and (ii) the cladding temperature is below 1973 K and the oxygen concentration is high in the b-Zr phase. These criteria were applied to brittle failure. The following criteria in the b-Zr phase are proposed by Pawel of Oak Ridge National Laboratory (ORNL) (U.S. Nuclear Regulatory Commission, 2001) based on PBF in-pile reactor experiment of EG&G Idaho, Inc:  High speed cooling: >100 K/s – Solid solubility of oxygen concentration at starting of quench >90% – Average oxygen concentration >0.65 wt% – Cladding temperature >1973 K  Low speed cooling: – Oxygen concentration >0.3 mm.

>1.0 wt%

with

b-Zr

thickness

point: 2963 K and prevents the relocation of molten Zr cladding at high temperature conditions. Hence the ZrO2 thickness was selected for the cladding melting parameter as an engineering viewpoint. From the irradiation experience of BWR high burnup 9x9 fuel (Hirano et al., 2005), the average oxidized ZrO2 thickness is about 10 lm. At 2963 K, the ZrO2 thickness was assumed then 10 lm from the irradiation result. Also no oxide layer was assumed for the low temperature criterion. The inner side of the ZrO2 layer is a-Zr(O) layer which has the low melting temperature, 2248 K. The temperature is the starting point of liquid phase of Zr-ZrO2 phase diagram (INEEL, 2003). The criterion of cladding melting failure therefore was interpolated between 2248 K and 2963 K:

f ox ¼ 1:02  103 ðT K  2248Þ1:4 ; f melt ¼

4.2.1.3. Cladding melting. After a normal operation, it is generally acknowledged that the Zr cladding is cover with an oxidized Zr, ZrO2 uniformly. The oxide layer on the surface has a high melting

f ox ; Thox

ð11Þ ð12Þ

where Thox is the evaluated outer oxidized thickness (micrometer), Tk is the cladding temperature (K), fox is the failure thickness

189

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

1 cell

Top View

1 cell

Side View

20 cells

Liquid fuel Fuel particle 16 cells

Core support plate Control rod guide tube

5 cells

Coolant (Water)

9 cells

Control rod drive housing

6 cells

Bottom-end of reactor pressure vessel

56 cells in the axial direction

(a) Calculation geometry

Initial condition Fuel particle : 30% vol. Fuel liquid : 5% vol.

Freezing of fuel liquid to fuel particle

Relocation of fuel particle downward due to gravity

Vaporization of water due to the interaction with hot fuel particle

Progress of fuel particle deposition on the surface of bottom-end of reactor pressure vessel

Relocation of fuel particle downward

0 sec

1 sec

Deposition of fuel particle on the surface of bottom-end of reactor pressure vessel

2 sec

3 sec

4 sec

9 sec

(b) Verification results Fig. 4. Verification geometry and results for high temperature molten and particle fuel into water.

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

120

2,500

2.0 Oxygen concentration (%)

Failure

Cladding plastic strain (%)

Limit of brittle fracture: 0.65 (wt%)

100

Cladding temperature (K)

Oxygen concentration (wt%)

Cladding strain (%)

Limit of cladding burst strain (%) 80

60

40

1.5

2,000

Failure 1,500

1.0 1,000

0.5

Cladding temperature (K)

190

500

20

0

0.0 0

1000

2000

3000

4000

0

1000

2000

0 4000

3000

Time(s)

Time (s)

(a) Cladding rupture failure

(b) Cladding brittle fracture 3000

3,500

Failure

Failure

3,000

2500

Temperature (K)

Temperature (K)

2,500

2,000

1,500

2000

1500 Pellet average temperature (K) 1000

Limit of pellet dissolution failure: 2500(K)

1,000

500

500

Cladding surface temperature (K) Limit of cladding melting failure (K)

0

0

0

100

200 300 Time(s)

400

500

0

1000

2000

3000

4000

5000

6000

7000

Time (s)

(c) Cladding melting failure

(d) Pellet dissolution failure Fig. 5. Verification results of fuel pin failures.

(micrometer), if fmelt P 1, the cladding melting failure will occur. The fitting coefficient between 2248 K and 2963 K was assumed as the liquidus of the diagram. For this reason, if the oxidized layer is very thin, the brittle oxide layer will collapse/fail due to the internal pressure by the inner molten Zr. 4.2.1.4. Pellet dissolution. In the case of low-level heating rate, the temperature distribution in a pellet is almost flat; UO2 and Zr form the liquid eutectic inside the cladding. From the PHEBUS knowledge (Ikeda et al., 2003), the cladding failure occurred around 2500 K. Hence when an averaged volume fuel temperature reaches 2500 K, all the cladding and fuel relocates to a coolant flow channel in the calculation mesh with low-level heating, like decay heat. Additionally in the case of high-level heating rate, the foregoing UO2-Zr eutectic would not mature and hence another criterion is introduced the melting point (3123 K) of UO2 pellet.

4.2.2. Verification calculation Fig. 5 shows the results of verification calculation about four failure cases. For the all calculations, the fuel pin was assumed to be the 9x9 BWR fresh fuel pin. The initial pressures of coolant and fuel plenum were 0.1 MPa and 3 MPa respectively. From the following verification, it was verified that the introduced four failure criterion works properly for fuel pin failure conditions. 4.2.2.1. Cladding rupture failure (see Fig. 5(a)). The plenum pressure increased up to 9.5 MPa with cladding temperature rise. The figure shows the plastic strain and burst strain profiles. At the cladding rupture failure (3940 s), the cladding plastic strain reached the criterion. The cladding plastic strain was then found from calculation to be 53.8% which was high enough for coolant flow blockage. 4.2.2.2. Cladding brittle fracture (see Fig. 5(b)). The figure shows the cladding brittle failure after quench in the case of high speed

191

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

cooling. After heating the cladding to 1600 K, the average oxygen concentration in the b-phase was over 0.65wt% which corresponded to the criterion of cladding brittle failure. Hence the cladding brittle failure occurred immediately after quench (3500 s). 4.2.2.3. Cladding melting failure (see Fig. 5(c)). The cladding surface temperature profile is shown in figure. The cladding melting failure occurred at 415 s, since the cladding temperature reached 2963 K (limit of cladding melting failure). 4.2.2.4. Pellet dissolution failure (see Fig. 5(d)). The eutectic failure between UO2 and Zr occurred at 6567 s during cladding temperature rise in the figure, since the average fuel temperature attained 2500 K. 4.3. Neutronic module The input/output system was programed for a nodal neutron diffusion code, SKETCH developed by Japan Atomic Energy Agency (JAEA). In order to reduce the Multifunction Model’s calculation time, the data table of macroscopic cross-section was prepared

for the SKETCH code for several fuel conditions in advance. The SRAC code (Okumura et al., 2007) with JENDL-4 nuclear library database (Nuclear Data Center, 2014) generates the macroscopic cross-section table for intact cores, degraded cores and debris beds. The following three conditions are considered as SRAC calculations: (i) intact fuel condition in a core zone, (ii) un-intact fuel condition in a core zone, (iii) fuel debris condition under a core zone (at a bottom of reactor pressure vessel). 4.3.1. Calculation method The macroscopic cross-section data cover all the fuel/debris compositions for a typical BWR. The arbitrary fuel and debris cross sections for the SKETCH code are evaluated by interpolating and extrapolating the foregoing macroscopic cross section data table. The cross section table is 2-energy-group: thermal group (105– 1.8554 eV) and fast group (1.8554–107 eV). The maximum simulation geometry is a 1/4 core geometry with a fuel assembly-wise. An infinite medium multiplication factor, kinf of volume-averaged in a calculation cell is used as the eigen value indication, instead of the effective multiplication factor to reduce the neutronic calculation time. If infinite medium multiplication factors kinf of volumeaveraged are below 0.90 in every cell of the vessel, the neutronic

100 (%)

10.0

8.0

6.0

4.0

2.0

Core zone (Radial direction)

Ratio of the calculation result of SKETCH to MVP (%)

0.0 0

20

40

60

80

100

Coolant void ratio (%)

Ratio of the calculation result of SKETCH to MVP (%)

10.0

100 (%)

Core zone (Axial direction)

8.0

6.0

4.0

In-Vessel (RPV)

Near core support plate

2.0

Ratio of the calculation result of SKETCH to MVP (%)

(Radial direction) 0.0 0

1000

2000

3000

Fuel temperature (K) Ratio of the calculation results of SKETCH to MVP (%)

(a) Three dimensional calculation geometry of Monte Carlo Code, MVP.

(b) Comparing SKETCH calculation results with that of MVP for core degradation conditions.

Fig. 6. Verification results of neutronic calculation between SKETCH and MVP codes.

192

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

module is not active and indicates the value of infinite medium multiplication factor kinf alone as an engineering judgment. Moreover, a first-order perturbation method is applied to reactor power analyses. The infinite medium multiplication factor kinf timedepend is defined as,

kinf ¼

ðmRf 1 /1 þ mRf 2 /2 Þ ; ðRa1 /1 þ Ra2 /2 Þ

ð13Þ

where m is the average number of fission neutrons released as the result of fission, Rfi is the macroscopic fission cross-section in m2,

Upper head

k = 36 43 (7 cells)

Dryer

k = 34 36 (3 cells)

Separator

k = 29 33 (5 cells) k = 26 28 (3 cells)

Shroud head

k = 22 25 (4 cells)

Top of active fuel k = 14 21 (8 cells) Bottom of active fuel Control rod guide tube

k = 8 13 (6 cells) k = 5 7 (3 cells)

Control rod drive housing

k = 1 4 (4 cells)

Lower head

(a) Axial calculation mesh in a reactor pressure vessel (RPV).

(b) Radial calculation mesh in a reactor pressure vessel (RPV). Fig. 7. In-vessel calculation geometry of the Multifunction Model.

193

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

Rai is the macroscopic absorption cross-section in m2, ui is the neutron flux in neutrons/m2s, i is the energy group. 4.3.2. Verification calculation For verification of the SRAC and SKETCH codes of the module, the calculated effective multiplication factor keff (coolant void reactivity and Doppler effects) for the 1/4 core geometry was verified by MVP (Nagaya et al., 2005): General Purpose Monte Carlo Codes for Neutron and Photon Transport Calculations developed at Japan Atomic Energy Research Institute (JAERI, at the present JAEA) with a full core geometry. A very fine MVP calculation model inside a BWR reactor pressure vessel has been made as shown in Fig. 6 (a). The model includes the core zone, the core support plate zone and also the upper and lower head of reactor pressure vessel. For verification, the annular core geometries (the degradation ratio of 70% for the inner fuel assemblies and control rods) were selected for simulating a core degradation phase; the central part of core (fuel assemblies and control rods) was missed as a hypothetical assumption. The results of the coolant void and Doppler reactivities are shown in Fig. 6(b). The coolant void ratio and fuel temperature range were from 0 to 100% and from 293 to 3273 K, respectively. The ratio of the coolant void reactivity of SKETCH to MVP was from 5% to 8% overall range of the coolant void ratio. The ratio of the Doppler reactivity of SKETCH to MVP was constantly about 6% throughout a range of the fuel temperatures. The diffusion calculation method of the SKETCH code was slightly difficult to solve the unusual annular core geometry such as severe accident conditions. No significant differences between the SKETCH and MVP however were observed in the coolant void and Doppler reactivities for the degraded geometry. 5. Tight coupling method and geometry, and verification 5.1. Tight coupling method and geometry The Multifunction Model consists of the neutronic diffusion calculation module (based on SKETCH), the thermal-hydraulic calculation module and the fuel pin behavior calculation module Initial condition

Coolant boiling

Fuel melt

(based on FEMAXI-6) in order to evaluate steady-states, initial fuel deformation phases and fuel melting phases. The three modules communicate with each other through a ‘‘control module” and the data are exchanged among the three modules in optimized time steps at the state of the in-vessel. The main transferred data are spatial distributions of materials and temperatures, power density and so on. The Multifunction Model covers a BWR core zone and complicated structures under the core zone (the core support plate, the control guide tubes and the bottom of reactor pressure vessel). The calculation meshes (U.S.NRC, 2016) are same in all the modules as shown in Fig. 7 for the maximum simulation geometry: 1/4 core dimensions of a typical BWR reactor pressure vessel. The axial direction divided into 43 meshes and the relatively-fine mesh is applied to the core and the lower head of reactor pressure vessel, since these zones are essential for evaluating core degradation. For the radial direction, the five to nine fuel assemblies are lumped together as a unit fuel assembly cell as shown by the red lines. 5.2. Verification of tight coupling simulation for single BWR fuel assembly equivalent model using three modules All of the modules were connected by the control module. The preliminary calculation was executed using a BWR single fuel assembly equivalent model to verify tight coupling calculation in a hypothetical condition. The simulation covered wide ranges from a water flooding condition to a relocation condition of molten materials. The calculation zone of the single BWR fuel assembly equivalent model was from the bottom to the top of a typical BWR reactor pressure vessel. The number of calculation cells was 56 in the axial direction same as Fig. 4(a) calculation geometry. In the three-dimensional calculation geometry, the single fuel assembly alone based on the 9x9 BWR fuel assembly was arranged without control rods. The seventy-two fuel pins were lumped together as a single fuel pin. The heat power of the fuel rod was assumed to be a decay heat level calculated by the widely known Way-Wigner equation (Way and Wigner, 1948). The assumed decay heat was relatively higher than that of a normal severe Rupture of fuel rod

Melt relocation downward

Melt relocation to lower zone

Separator

Coolant boiling

Coolant

Fuel rod Core support plate

Melt relocation downward

Rupture of fuel rod

Fuel melt

Center of the active core

Control rod guide tube

Molten materials on the surface of bottom-end of reactor pressure vessel

On the core support plate

Bottom-end of reactor pressure vessel

The bottom of RPV

0.1 sec

301 sec

1552 sec

1553 sec

1605 sec

1607 sec

Volume fractions of components at this position are shown below.

Fig. 8. Core degradation progress: volume fraction of materials at the important calculation steps.

194

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

Melt components generated 1-cell above the core center at 1552 sec are passing.

Melt components generated at the core center at 1605 sec are passing. Melt components generated 2-cell above the core center at 1598 sec are passing.

Volume fraction of the materials at the central part of fuel

Melt components generated 1-cell above the core center at 1552 sec are passing.

Melt components generated at the core center at 1605 sec are passing.

Melt components generated 2-cell above the core center at 1598 sec are passing.

Volume fraction of the materials on the core support plate

Volume fraction of the materials on the surface of bottom-end of reactor pressure vessel Fig. 9. Volume fraction of materials at the important zones.

T. Okawa, T. Nakajima / Annals of Nuclear Energy 101 (2017) 182–195

accident condition, since the reactor was shut down at the starting time of this calculation. The axial power distribution was assumed ‘‘as chopped cosine.” The initial pressure and temperatures of gas, structure and coolant were 7 MPa and 559 K, respectively. The initial water level was at the Top of Active Fuel (TAF). Fig. 8 illustrates, the result of volume fraction of each component, showing the representative degradation progress of the single fuel assembly equivalent model. After fuel rupture at 1552 s due to meet the criterion of pellet dissolution failure in the fuel pin behavior module (2500 K), the molten fuel materials flowed downward owing to gravity after 1553 s. The figure indicated that the fuel elements were properly transferred from the fuel pin behavior module to the thermal-hydraulic module. The calculation result suggested that the calculation progress was rational and qualitative as coupling of thermal-hydraulic and fuel pin behavior modules. Fig. 9 shows the volume fraction profile at the central part of the fuel, on the core support plate and on the surface of bottom-end of reactor pressure vessel, respectively. In these figures, the volume of structural materials was not shown for clarity to reveal core material relocation. In the central part of fuel, at 1553 s after inputting decay heat, a small amount of fuel chunk (orange), Zr particle (dark green) and eutectic liquid (light blue) appeared. The relocating fuel elements from the central part of fuel induced by decay heat appear at 1598 and 1605 s intermittently. The volume fraction of the fuel rod (black) suddenly decreased due to relocation of all of the fuel elements in the cell and the volume of fuel rod was substituted by that of gas (gray) at 1605 s. On the core support plate, relocation of the very hot fuel elements induced rapid boiling from 1550 to 1560 s and the volume of coolant (blue) decreased. Finally, the traveling fuel elements reached on the surface of bottom-end of reactor pressure vessel. The fuel chunk and Zr particle gradually accumulated on the surface of bottom-end of reactor pressure vessel. The results are not quantitative but qualitative by proven event-advancements. In the near future, quantitative evaluations for larger or more complicated geometries: (i) four fuel assemblies with unit fuel support piece and (ii) 1/4 core geometries will be carried out. 6. Conclusions For the purpose of reducing uncertainties in severe accident simulation, the Multifunction Model has been developed for simulating BWR in-vessel core degradations since 2012. The Multifunction Model permits to decrease the uncertainties of assessment that are key issues of current severe accident codes for the in-vessel core degradation, since the advanced physical models were incorporated. The thermal-hydraulic module was evaluated for high temperature molten and particle fuel into water. The fuel pin behavior module was evaluated in four fuel pin failure modes for the Japanese 9x9 BWR fuel pin. The neutronic module was evaluated in the hypothetical core geometry, the annular core geometry. The tight coupling analysis of the neutronic, thermal-hydraulic and fuel pin behavior modules was assessed in the three-dimensional, single fuel assembly equivalent geometry. The verification results revealed that the Multifunction Model demonstrated a reasonable capability to simulate the BWR in-vessel core degradations. As an engineering challenge, the validation of Multifunction Model will be undertaken for quantitative simulations by means of the existing severe accident tests. The papers related to foregoing validation experiments will be published in the near future. Acknowledgements The Multifunction Model analysis has been performed with O. Kitamura of Mizuho Information & Research Institute (MHIR,

195

Japan). The researchers for the collaboration work are S. Hagen, P. Hofmann, W. Hering, J. Stuckert and A. Miassoedov of Karlsruhe Institute of Technology (KIT). The authors appreciate their important suggestions for the material relocation during initial core degradation. References Bohl, W.R., Luck, L.B., 1990. SIMMER-II: A Computer Program for LMFBR Disrupted Core Analysis. LA-11415-MS. Los Alamos National Laboratory. Chatelard, P. et al., 2014. ASTEC V2 Severe Accident Integral Code Main Features, Current V2.0 Modelling Status, Perspective. Nucl. Eng. Des. 272, 119–135. Fauske & Associates, LLC, 2016. MAAP-Modular Accident Analysis Program. Available from: . Gauntt, R.O., Humphries, L.L., 1997. Final Results of the XR2-1 BWR Metallic Melt Relocation Experiment. SAND-97-1039. Sandia National Laboratories (SNL). Geelhood, K.J. et al., 2011. FRAPTRAN 1.4: A Computer Code for the Transient Analysis of Oxide Fuel Rods. NUREG/CR-7023, Vol. 1. United States Nuclear Regulatory Commission (U.S.NRC). Hagen, S. et al., 2015. Private Communication. Karlsruhe Institute of Technology (KIT). Hering, W., Hofmann, C., 2007. Degraded core reflood: present understanding and impact on LWRs. Nucl. Eng. Des. 237, 2315–2321. Hirano, Y., et al., 2005. Irradiation Characteristics of BWR High Burnup 9x9 Lead Use Assemblies. In: Proc. 2005 Water Reactor Fuel Performance Mtg. Kyoto, Japan. Hofmann, P., 1999. Current knowledge on core degradation phenomena, a review. J. Nucl. Mater. 270, 194–211. Ikeda, T., Terada, M., Karasawa, H., 2003. Analysis of core degradation and fission products release in phebus FPT1 test at IRSN by detailed severe accidents analysis code, IMPACT/SAMPSON. J. Nucl. Sci. Technol. 40, 591–603. Idaho National Engineering and Environmental Laboratory (INEEL), 2003. SCDAP/ RELAP5-3D Code Manual Volume 4: MATPRO – A Library of Materials Properties for Light-Water-Reactor Accident Analysis. INEEL/EXT-02-00589. Ishii, M., Hibiki, T., 2011. Thermo-Fluid Dynamics of Two-Phase Flow. Springer. Japan Atomic Energy Agency (JAEA), 2013. SKETCH-N 1.0: Solve Neutron Diffusion Equations of Steady-State and Kinetics Problems. Available from: . Nagaya, Y. et al., 2005. MVP/GMVP II: General Purpose Monte Carlo Codes for Neutron and Photon Transport Calculations based on Continuous Energy and Multigroup Methods. JAERI 1348. Japan Atomic Energy Research Institute (JAERI). Nuclear Data Center, 2014. JENDL-4.0. Available from: . Okawa, T., Nakajima, T., 2016. Multifunction Model Features and Current Status for BWR Core Degradation. In: 2016 International Congress on Advances in Nuclear Power Plants (ICAPP 2016). San Francisco, CA, USA. Okumura, K. et al., 2007. SRAC2006: A Comprehensive Neutronics Calculation Code System. JAEA-Data/Code 2007-004. Japan Atomic Energy Agency (JAEA). Satoh, N. et al., 2000. Development of Molten Core Relocation Analysis Module MCRA in the Severe Accident Analysis Code SAMPSON. J. Nucl. Sci. Technol. 37, 225–236. Sehgal, B.R., 2012. Nuclear Safety in Light Water Reactors – Severe Accident Phenomenology. Academic Press. Seiler, N. et al., 2008. Investigations on Boron Carbide Oxidation for Nuclear Reactors Safety-General Modelling for ICARE/CATHARE Code Applications. Nucl. Eng. Des. 238, 820–836. Suzuki, M., Saito, H., 2006. Light Water Reactor Fuel Analysis Code FEMAXI-6 (Ver.1); Detailed Structure and User’s Manual. JAEA-Data/Code 2005-003. Japan Atomic Energy Agency (JAEA). U.S. Nuclear Regulatory Commission, 2001. SCDAP/RELAP5/MOD 3.3 Code Manual, Code Architecture and Interface of Thermal Hydraulic and Core Behavior Models. NUREG/CR-6150. U.S. Nuclear Regulatory Commission, 2015. MELCOR Computer Code Manuals. SAND2015-6692R, Vol. 2, Ver. 2.1.6840 2015. U.S. Nuclear Regulatory Commission, 2016. Part II Introduction to Reactor Technology-BWR. Available from:. Veshchunov, M.S., Palagin, A.V., 1998. Modeling of Chemical Interactions of Fuel Rod Materials at High Temperatures II. Investigation of Downward Relocation of Molten Materials. J. Nucl. Mater. 252, 110–120. Veshchunov, M.S., Shestak, V.E., 2008. Model for Melt Blockage (Slug) Relocation and Physico-Chemical Interactions during Core Degradation under Severe Accident Conditions. Nucl. Eng. Des. 238, 3500–3507. Veshchunov, M.S. et al., 2013. Modelling the formation and oxidation of molten pools. Ann. Nucl. Energy 61, 54–62. Way, K., Wigner, E.P., 1948. The rate of decay of fission products. Phys. Rev. 73, 1318–1330.