# Numerical Simulation and Validation for Early Core Degradation Phase under Severe Accidents

Authors: Dekui Zhan; Xinhai Zhao; Shaoxiong Xia; Peng Chen; Huandong Chen
Journal: Science and Technology of Nuclear Installations
 (2020)

---

## Abstract

Early core degradation determines the amount of hydrogen generated by cladding oxidation as well as the temperature, the mass, and the composition of corium that further relocates into the lower head of reactor pressure vessel (RPV), which is essential for the effectiveness analysis of in-vessel retention (IVR) and hydrogen recombiners. In this paper, the mechanisms of controlling phenomena in the early phase of core degradation are analysed at first. Then, numerical models adopted to calculate (1) core heating up, (2) cladding oxidation, (3) dissolution between molten zirconium and fuel pellets, and (4) formation of a molten pool in the core active section are presented. Compared with integral codes for severe accident analysis (such as MAAP and MELCOR), the models in this paper are established at the fuel pin level and the calculation is performed in 3D, which can capture the detail local phenomena during the core degradation and eliminate the average effect due to equivalent rings used in integral codes. In addition, most of the control equations in this paper are calculated by implicit schemes, which can improve the accuracy and stability of the calculation. In the simulation, the calculation oxidation is calculated by using the oxygen diffusion model, while the dissolution is calculated with Kim, Hayward, Hofmann, and IBRAE models to perform uncertainty analysis. For the validation, the cladding oxidation model is verified by Olander theoretical cases in the conditions of both steam-rich and steam-starved. The dissolution models are validated by the RIAR experiment. The code is overall verified by Phebus FPT0 on the integral phase of core early degradation. According to the simulation results, it can be inferred that the dissolution reaction between the molten zirconium and fuel pellets is the main reason for the melting of UO2 at low temperature. In the case of starved steam, part of the fuel pellets can melt down even at 2248 K and relocate to the bottom of the core, which is much lower than the melting point of UO2 (3113 K).


---

## Body

## 1. Introduction

In the Chernobyl and Fukushima accidents, the melted reactor cores caused a large amount of radioactive material to be released into the environment. It is essential to study the severe accident phenomena for the design of dedicated mitigation measures to cease the accident process and reduce its consequences.Accidents that can cause the core to melt down mainly include (1) the station blackout accident (SBO), (2) the large break loss of coolant accident (LBLOCA) superimposed with the failure of safety injection pumps, and (3) the anticipated transient without scram (ATWS) [1]. Although the initial events are different, the phenomena in the core early degradation are approximately identical. The controlling phenomena include core uncovery, core heating up, failure of control rods, dissolution of UO2 and ZrO2 by molten zirconium, failure of the cladding oxide layer, corium relocation and solidification, and formation of a molten pool in core active section.For the cladding oxidation, under the condition that the solid zirconium alloy is thick enough and the steam is sufficient, the oxidation of zirconium alloy is dominated by the oxygen diffusion, in which the chemical reaction rate follows the parabolic law [2]. Correlations of Urbanic and Heidrick [3], Prater and Courtright [4], and Leistikow-Schanz [5] are widely used in the calculation of cladding oxidation.De Luze et al. [6] summarized the failure temperatures of Ag-In-Cd control rods in different experiments, including CORA-9, QUENCH-13, and Phebus FPT1. The failure temperatures of control rods in most experiments are between 1400 and 1700 K. It can be concluded that the eutectic reaction between the solid stainless steel and the zirconium alloy is the reason for the partial melting of the cladding and guide tube at about 1400 K, while the dissolution of zirconium by the molten stainless steel and Ag-In-Cd absorber causes the guide tube to be broken in a wide range at around 1700 K.Crucible experiments found that the dissolution process of UO2 can be divided into two phases, including the saturation phase and the precipitation phase. Hofmann [7], Hayward and George [8], and Kim and Olander [9] established separately dissolution models based on their experiments, respectively. However, significantly different results were discovered during the precipitation period. The dissolution of UO2 reaches saturation in Kim and Hayward experiments while Hofmann found the dissolved UO2 mass continues to increase and follows the parabolic law. The ZrO2 dissolution process is also separated into two stages, including the erosion stage and the corrosion stage [10], in which the correlation of Hofmann is widely used.Based on the experimental data and mechanism research studies, many integral codes for severe accident analysis were developed, including MELCOR, SCDAP/RELAP5, and ASTEC. The temperature criteria and the oxide thickness criteria are widely used in most codes to predict the failure moment of the cladding oxide layer. For the relocation process of melt after the failure of the cladding oxide layer, candling model is established in MELCOR by analogy with the flow of molten droplets of the candle along the wall [11]. LIQSOL model is developed in SCDAP/RELAP5 based on the experimental results of CORA, PBF-SFD, and FLHT [12]. Decanting model and porous media model are adopted in ASTEC, which can simulate the radial and axial relocation of corium simultaneously [13].Due to the consideration of computing efficiency, equivalent rings are widely used in integral codes. However, radial power factors of assemblies in the same ring are different, especially between assemblies with black rods and without black rods. When equivalent rings are adopted in the modelling, the power factors of high-power assemblies are actually averaged, causing the delay of the calculated degradation process. It means that integral codes may underestimate the time of some important events in the reactor, such as the start time of cladding oxidation, the failure of the zirconium cladding (namely, the fission product begins to be released in large quantities), and the downward relocation of the molten materials, which are important for accident analysis. In this paper, the models are established at the fuel pin level and the calculation is performed in 3D which can capture more detailed local phenomena during the core degradation phase and eliminate the average effect due to equivalent rings used in integral codes to achieve more accurate simulation for the early phase of core degradation.


## 2. Mechanism of Early Core Degradation

In this section, the mechanisms of controlling phenomena in the early core degradation phase are presented based on the literature review.


### 2.1. Oxidation of Zirconium Cladding

The reaction of zirconium cladding with steam is the most important phenomenon in the core degradation, which generates not only a large amount of heat but also a great quantity of hydrogen. With neglecting considering the hydrogen absorption in the cladding, the zirconium oxidation can be divided into the following stages [7]: (1) diffusion of steam to the outer surface of cladding; (2) chemical absorption of steam; (3) dissociation of steam molecules into oxygen and hydrogen; (4) diffusion of oxygen from the outer surface of cladding to the inside, forming ZrO2, α-Zr(O), and β-Zr three layers in sequence; (5) recombination of hydrogen; and (6) diffusion of hydrogen molecules into the bulk of the flow channel, causing the hydrogen content in the flow channel to increase.


### 2.2. Dissolution of UO2 and ZrO2

Based on the analysis of the experimental results of Hayward, Kim, and Hofmann, Olander [14] proposed that during the saturation period, the UO2 dissolution is dominated by natural convection mass transfer. Liquid U with a density of 1.7 × 104 kg/m3 is produced by a chemical reaction at the interface between the molten zirconium alloy and the solid UO2. However, the density of the liquid Zr is only 6 × 103 kg/m3. The difference in density between liquids U and Zr causes significant convection in the melt. During this period, the U content in the melt increases continuously and gradually reaches the saturation. Dissolved U and O are quickly transferred to the central region of the melt and mixed well with the molten Zr. Therefore, it is found in the experiments that the ratios of U, Zr, and O are uniformly distributed in the molten zone. In the precipitation stage, Olander does not suppose that UO2 will continue to be dissolved. At the end of the saturation period, the dissolution reaction of UO2 is ceased, and the mass ratio of U and Zr in the melting zone is kept unchanged. During this period, O in the UO2 crucible continues to diffuse into the molten zone, causing the previously saturated (U, Zr, O) melt to be supersaturated to precipitate the (U, Zr)O2-x ceramic phase.The dissolution of ZrO2 is similar to that of UO2. Based on the ZrO2 crucible experiments, Hayward and George [10] speculated that in the erosion period, the dissolution of ZrO2 is caused by a significant convective effect. The newly dissolved ZrO2 is quickly brought into the melt central zone. Therefore, no transition layer was created in the experiment in this period. Thereafter, the oxygen concentration gradient at the interface gradually decreases causing the convective effect to be weakened. Once the mass flux of oxygen entering the melt by diffusion is greater than that by dissolution, the erosion period is ceased. In the following period, oxygen in the ZrO2 crucible continues to diffuse into the melt driven by the difference of oxygen potential. Near the interface region, ZrO2 in the crucible is reduced to ZrO2-x and the metallic zirconium in the melt is oxidized to ZrO2−x. Therefore, it was shown in the experiment that the thickness of the ZrO2−x transition layer increased not only towards the crucible region but also towards the melt central region.


### 2.3. Cladding Oxide Layer Failure

The cladding of fuel rods is oxidized by steam at high temperature. Then, a ZrO2 oxide layer is formed on the outside of the cladding. The oxide layer with a high melting point (2963 K) will hold up the molten (U, Zr, O) in the cladding, preventing the downward relocation of the melt. However, the ZrO2 oxide layer is dissolved by the metallic Zr in the molten (U, Zr, O), causing the oxide layer to become thinner. In addition, the creep effect at high temperatures reduces the ability of the oxide layer to hold up the melt and its yield limit decreases with increasing temperatures. Therefore, when the thickness of the oxide layer is reduced to a certain threshold, the oxide layer can no longer hold up the (U, Zr, O) melt. Then, the (U, Zr, O) melt is released to relocate downwards.


### 2.4. Corium Relocation

According to CORA W1 and W2 experimental results, Veshchunov and Palagin [15] infer that in the first period, the molten cladding materials flow down rapidly in the form of droplets or rivulets at about 0.5 m/s. Owing to a large temperature gradient below the front, the molten materials refreeze in the interrod space forming local blockages (debris) and accumulating the melt in the form of molten pool. The melt mass in the molten pool gradually increases. Once the solidified crust at the bottom cannot support the weight of the melt, the molten pool will relocate downwards. During the relocation process, the molten zirconium alloy in the molten pool will dissolve adjacent fuel pellets and be oxidized by steam. The dissolution reaction and the oxidation reaction cause U in the melt to reach a supersaturated state to produce the (U, Zr)O2−x ceramic phase precipitate. The presence of the (U, Zr)O2−x ceramic phase precipitate increases the viscosity of the melt, thereby preventing the molten pool from rapidly relocating towards the bottom of the bundle. At the same time, the oxidation reaction generates a large amount of heat raising the temperature of the melt. Therefore, the solubility of U in the melt is increased, resulting in the reduction of the (U, Zr)O2−x ceramic phase and the viscosity of the melt. Consequently, the molten pool moves downwards with a slow characteristic velocity 1∼2 × 10−3 m/s. As the molten pool relocates downwards, fresh melt continues to accumulate in the molten pool, causing the molten pool to expand in the radial and axial directions to form a massive slug flow.


## 2.1. Oxidation of Zirconium Cladding

The reaction of zirconium cladding with steam is the most important phenomenon in the core degradation, which generates not only a large amount of heat but also a great quantity of hydrogen. With neglecting considering the hydrogen absorption in the cladding, the zirconium oxidation can be divided into the following stages [7]: (1) diffusion of steam to the outer surface of cladding; (2) chemical absorption of steam; (3) dissociation of steam molecules into oxygen and hydrogen; (4) diffusion of oxygen from the outer surface of cladding to the inside, forming ZrO2, α-Zr(O), and β-Zr three layers in sequence; (5) recombination of hydrogen; and (6) diffusion of hydrogen molecules into the bulk of the flow channel, causing the hydrogen content in the flow channel to increase.


## 2.2. Dissolution of UO2 and ZrO2

Based on the analysis of the experimental results of Hayward, Kim, and Hofmann, Olander [14] proposed that during the saturation period, the UO2 dissolution is dominated by natural convection mass transfer. Liquid U with a density of 1.7 × 104 kg/m3 is produced by a chemical reaction at the interface between the molten zirconium alloy and the solid UO2. However, the density of the liquid Zr is only 6 × 103 kg/m3. The difference in density between liquids U and Zr causes significant convection in the melt. During this period, the U content in the melt increases continuously and gradually reaches the saturation. Dissolved U and O are quickly transferred to the central region of the melt and mixed well with the molten Zr. Therefore, it is found in the experiments that the ratios of U, Zr, and O are uniformly distributed in the molten zone. In the precipitation stage, Olander does not suppose that UO2 will continue to be dissolved. At the end of the saturation period, the dissolution reaction of UO2 is ceased, and the mass ratio of U and Zr in the melting zone is kept unchanged. During this period, O in the UO2 crucible continues to diffuse into the molten zone, causing the previously saturated (U, Zr, O) melt to be supersaturated to precipitate the (U, Zr)O2-x ceramic phase.The dissolution of ZrO2 is similar to that of UO2. Based on the ZrO2 crucible experiments, Hayward and George [10] speculated that in the erosion period, the dissolution of ZrO2 is caused by a significant convective effect. The newly dissolved ZrO2 is quickly brought into the melt central zone. Therefore, no transition layer was created in the experiment in this period. Thereafter, the oxygen concentration gradient at the interface gradually decreases causing the convective effect to be weakened. Once the mass flux of oxygen entering the melt by diffusion is greater than that by dissolution, the erosion period is ceased. In the following period, oxygen in the ZrO2 crucible continues to diffuse into the melt driven by the difference of oxygen potential. Near the interface region, ZrO2 in the crucible is reduced to ZrO2-x and the metallic zirconium in the melt is oxidized to ZrO2−x. Therefore, it was shown in the experiment that the thickness of the ZrO2−x transition layer increased not only towards the crucible region but also towards the melt central region.


## 2.3. Cladding Oxide Layer Failure

The cladding of fuel rods is oxidized by steam at high temperature. Then, a ZrO2 oxide layer is formed on the outside of the cladding. The oxide layer with a high melting point (2963 K) will hold up the molten (U, Zr, O) in the cladding, preventing the downward relocation of the melt. However, the ZrO2 oxide layer is dissolved by the metallic Zr in the molten (U, Zr, O), causing the oxide layer to become thinner. In addition, the creep effect at high temperatures reduces the ability of the oxide layer to hold up the melt and its yield limit decreases with increasing temperatures. Therefore, when the thickness of the oxide layer is reduced to a certain threshold, the oxide layer can no longer hold up the (U, Zr, O) melt. Then, the (U, Zr, O) melt is released to relocate downwards.


## 2.4. Corium Relocation

According to CORA W1 and W2 experimental results, Veshchunov and Palagin [15] infer that in the first period, the molten cladding materials flow down rapidly in the form of droplets or rivulets at about 0.5 m/s. Owing to a large temperature gradient below the front, the molten materials refreeze in the interrod space forming local blockages (debris) and accumulating the melt in the form of molten pool. The melt mass in the molten pool gradually increases. Once the solidified crust at the bottom cannot support the weight of the melt, the molten pool will relocate downwards. During the relocation process, the molten zirconium alloy in the molten pool will dissolve adjacent fuel pellets and be oxidized by steam. The dissolution reaction and the oxidation reaction cause U in the melt to reach a supersaturated state to produce the (U, Zr)O2−x ceramic phase precipitate. The presence of the (U, Zr)O2−x ceramic phase precipitate increases the viscosity of the melt, thereby preventing the molten pool from rapidly relocating towards the bottom of the bundle. At the same time, the oxidation reaction generates a large amount of heat raising the temperature of the melt. Therefore, the solubility of U in the melt is increased, resulting in the reduction of the (U, Zr)O2−x ceramic phase and the viscosity of the melt. Consequently, the molten pool moves downwards with a slow characteristic velocity 1∼2 × 10−3 m/s. As the molten pool relocates downwards, fresh melt continues to accumulate in the molten pool, causing the molten pool to expand in the radial and axial directions to form a massive slug flow.


## 3. Numerical Simulation Method

In this section, the physical models adopted in the simulation are presented for the phenomena mentioned above.


### 3.1. Bundle Heat Transfer

The core heating up is accompanied by heat conduction, convection heat transfer, and radiation heat transfer. For the convective heat transfer, the Dittus–Boelter forced convection correlation is adopted when there are drive pumps and the coolant flow rate is large, while the natural convection correlations are used in the case that the flow of coolant is driven by the buoyancy. The models introduced in the following sections are used to calculate the transient heat conduction in the fuel rods and the radiation heat transfer between rods.


### 3.2. Transient Heat Conduction

The two-dimensional axisymmetric unsteady heat conduction equation is used to solve the transient temperature distribution in the fuel rods and control rods. The governing equation is given by(1)ρcP∂T∂t=1r∂∂rλr∂T∂r+∂∂zλ∂T∂z+S˙,in which S˙ is the volumetric heat generation rate. For the absorbers in control rods, S˙=0, while for the fuel pellets,(2)S˙=HtotNfVfPzPr,in which Htot is the total power of core decay heat. Pz is the axial power factor. Pr is the radial power factor of fuel assemblies. Nf is the number of fuel rods in an assembly. Vf is the volume of a fuel rod in one axial node. The following equation is obtained by integrating equation (1):(3)aPTPn+1=aETEn+1+aWTWn+1+aNTNn+1+aSTSn+1+b.The coefficients in equation (3) are shown in (4)∼(6), where the geometric parameters are defined in Figure 1.(4)aE=rPΔrδzE/λE,aW=rPΔrδzW/λW,aN=rNΔzδrN/λN,aS=rSΔzδrS/λS,(5)aP=aE+aW+aN+aS+aP0,(6)aP0=ρcPΔVΔt,b=S˙ΔV+aP0TPn,ΔV=0.5rN+rSΔrΔz.Figure 1
Schematic of mesh generation for fuel pellets or absorbers.


### 3.3. Radiation Heat Transfer

The radiative heat transfer is calculated at the fuel pin level. Due to the shielding effect, central rodP only radiates heat directly to the surrounding 16 rods as shown in Figure 2. According to the geometric position, the surrounding rods can be divided into three types: 1, 2, and 3. The radiative view factors are X1, X2, and X3, respectively. Then, the effective radiation for central rod P is calculated as follows:(7)JPn+1=εσTP4+1−ε∑k=1,4J1,kn+1X1+1−ε∑k=1,4J2,kn+1X2+1−ε∑k=1,8J3,kn+1X3,where J1, J2, and J3 are the effective radiation at the outer surface of rods 1, 2, and 3, respectively. The radiative heat flux at the outer surface of rod P is given by(8)qP=ε1−εσTP4−JPn+1.Figure 2
Schematic of the radiation model between rods.


### 3.4. Cladding Oxidation

The oxygen diffusion model proposed by Olander [16] is adopted to calculate the oxidation of cladding, in which the cladding is divided into the oxide layer (ZrO2) and the metal layer (including a-Zr(O) and β-Zr). The oxidation rate depends on the amount of steam available, the thickness of the metal layer as well as the diffusion rate of oxygen in the cladding (namely temperature). According to the amount of available steam, the steam state in the flow channel can be divided into two cases: steam-rich and steam-starved. The amount of steam available is calculated by the molar steam balance:(9)VchMgdyidt=Fgyi−1−yi−Aclw˙,in which Mg=ptot/RT is the molar density of gas. yi is the steam mole fraction in the bulk gas in node i. The total gas molar flow rate is denoted by Fg and Acl is the external surface area of the cladding. The steam absorption rate is given by(10)w˙a=kgMgyi−ycl,where ycl is the steam molar fraction at the cladding outer surface. The oxygen diffusion rate in the cladding is given by(11)w˙d=ρZrDoCs−Caδ,in which the molar density of Zr in ZrO2 is denoted as ρZr and Do is the diffusion coefficient of oxygen in the oxide layer. Cs, Ca, and Cb are the O/Zr ratio at the outer surface of the oxide layer, at the interface in the oxide layer side, and at the interface in the metal layer side, respectively, which are determined by Zr-O phase diagram. The oxidation rate w˙ is calculated as the minimum of w˙a and w˙d. In the case of oxidation when both oxide and metal phases are present, the thickness of the oxide layer δ is calculated as follows:(12)δn+1=δn+ΔtGCa−Cb−DMCb4E2ξnc¯nCb+GDoδnCs−Ca,where DM is the diffusion coefficient of oxygen in α-Zr(O). The average oxygen content in the metal layer c¯ is determined by the following equation:(13)c¯n+1=c¯n+CbΔtDM2Eξn2c¯Cb−Gξn1−c¯nCbdδdtn+1.For the case of oxidation in all-metal cladding or in all-oxide cladding, other controlling equations are used to calculate the oxide thickness and the average oxygen content in the metal layer.


### 3.5. Dissolution of UO2 and ZrO2

Hofmann [7], Hayward and George [8], and Kim and Olander [9] established dissolution models based on their experiments, respectively. However, the results of each experiment in the precipitate period are significantly different. In the calculation, all the three models as well as the IBRAE model are used to simulate the dissolution of UO2 and ZrO2 and will be verified by the RIAR test in the next section.


#### 3.5.1. Kim UO2 Dissolution Model

When the Ra number in the (U, Zr, O) melt is greater than 1000, the natural convection model is adopted [9]:(14)CUn+1=CUS−CUS−CUn−1/4+0.54S4VgDU3νd1/4Δt−4.Otherwise, the diffusion model is used [9]:(15)CUn+1=CUn2+4DUπe2CUS2Δt0.5,in which CU is the mass fraction of U in the bulk melt. S, V, and ν are the UO2, (U, Zr, O) contact surface, the volume, and viscosity of (U, Zr, O), respectively. d is the diameter corresponding to S and e is the (U, Zr, O) mixture thickness.


#### 3.5.2. Hayward UO2 Dissolution Model

Based on crucible experiments, Hayward found that the mass fraction of U in the melt can be calculated as follows [8, 17]:(16)CUn+1=CUS−CUS−CUn·exp−kUΔt,in which kU is a constant related to the contact area and the initial oxygen content.


#### 3.5.3. Hofmann UO2 and ZrO2 Dissolution Model

The crucible experiments of Hofmann show that UO2 is rapidly dissolved in the initial period, and the mass fraction of UO2 in the melt rapidly increases to 35.8%. Subsequently, UO2 is continuously dissolved, and the UO2 mass fraction in the melt follows the parabolic law [7]:(17)CUO2n+1=0.358+CUO2n−0.3582+Kp1Δt0.5,(18)Kp1=1.02×1011exp−6.77×105RT.The ZrO2 dissolution model is calculated in a similar form of equation (7):(19)CZrO2n+1=0.233+CZrO2n−0.2332+Kp2Δt0.5,(20)Kp2=1.47×1010exp−6.66×105RT.


#### 3.5.4. Simultaneous Dissolution of UO2 and ZrO2

For the calculation of the simultaneous dissolution of UO2 and ZrO2, the first method is to assume it as two separate processes. Within one time step, the dissolution of UO2 is calculated firstly according to the models mentioned above, and then the dissolution of ZrO2 with Hofmann correlation is calculated until the mixture composition reaches the liquidus in the (U, Zr, O) ternary phase diagram. The second method is to adopt the IBRAE model proposed by Veshchunov et al. [18, 19], which can be used to calculate the simultaneous dissolution of UO2 and ZrO2 as well as the oxidation of (U, Zr, O) mixture. In the IBRAE model, both intact fuel rods and the molten pool can be considered as UO2, (U, Zr, O) mixture, and ZrO2 three-layer structure. The interface movement rate can be derived from the flux match. For the interface of UO2 and (U, Zr, O),(21)dr1dt=−DOUO2∂ρO/∂rI1−kOI1cOmax−cOBρO−g1+g2cMρU/cM.For the interface of (U, Zr, O) and ZrO2,(22)r3dr3dt−r2dr2dt=r2DOZrO2∂φO/∂rI2−kOI2cOmax−cOBφO−g1ρU/cM,where r1, r2, and r3 are positions of boundaries of the solid UO2 and ZrO2 layers. ρi, ci, and φi are the molar densities of different components (i=O, U, and Zr) in UO2, (U, Zr, O) mixture, and ZrO2, respectively. The molar density of the melt is denoted by cM. g1 and g2 are temperature-dependent parameters which determine the position of liquidus line in the ternary phase diagram.The oxygen concentration in the melt in (21) and (22) can be calculated by the mass conservation equations, while the oxygen concentration gradient at interface I1 is calculated by the diffusion equation of oxygen in solid UO2 with moving boundary conditions. For the oxygen concentration gradient at the interface I2,(23)if steam−rich,DOZrO2∂φ0∂rI2≈DOZrO2φ0B−φ0I1r3−r2,if steam−starved,DOZrO2∂φ0∂rI2≈kgMgyi.


### 3.6. Criteria of Cladding Oxidation Layer Failure

In the steam-rich case, the oxidation reaction in the cladding dominates. Oxidation of cladding continues until it is completely oxidized. In the calculation, it is assumed that the fuel pellets and cladding will melt together once the temperature reaches the eutectic point (2673 K).In the steam-starved condition, the dissolution reaction in the cladding dominates. In the calculation, it is assumed that the cladding oxide layer will fail once the (U, Zr, O) mixture temperature exceeds theα-Zr(O) melting point (2248 K) and the oxide layer thickness is less than 300 μm.


### 3.7. Corium Relocation

Based on the explanation of Veshchunov and Palagin on CORA W1 and W2 experimental results [15], the corium relocation downward in the bundle is considered as the slug flow with fixed velocity 1 × 10−3m/s. During the relocation, convective heat exchange with the steam, radiant heat exchange with the surrounding nonadjacent rods, and the heat generated by corium oxidation are considered. The energy equation during relocation is given by(24)ρcpdTmpVmpdt=Qox−AmphTmp−Tg−Qrad.


## 3.1. Bundle Heat Transfer

The core heating up is accompanied by heat conduction, convection heat transfer, and radiation heat transfer. For the convective heat transfer, the Dittus–Boelter forced convection correlation is adopted when there are drive pumps and the coolant flow rate is large, while the natural convection correlations are used in the case that the flow of coolant is driven by the buoyancy. The models introduced in the following sections are used to calculate the transient heat conduction in the fuel rods and the radiation heat transfer between rods.


## 3.2. Transient Heat Conduction

The two-dimensional axisymmetric unsteady heat conduction equation is used to solve the transient temperature distribution in the fuel rods and control rods. The governing equation is given by(1)ρcP∂T∂t=1r∂∂rλr∂T∂r+∂∂zλ∂T∂z+S˙,in which S˙ is the volumetric heat generation rate. For the absorbers in control rods, S˙=0, while for the fuel pellets,(2)S˙=HtotNfVfPzPr,in which Htot is the total power of core decay heat. Pz is the axial power factor. Pr is the radial power factor of fuel assemblies. Nf is the number of fuel rods in an assembly. Vf is the volume of a fuel rod in one axial node. The following equation is obtained by integrating equation (1):(3)aPTPn+1=aETEn+1+aWTWn+1+aNTNn+1+aSTSn+1+b.The coefficients in equation (3) are shown in (4)∼(6), where the geometric parameters are defined in Figure 1.(4)aE=rPΔrδzE/λE,aW=rPΔrδzW/λW,aN=rNΔzδrN/λN,aS=rSΔzδrS/λS,(5)aP=aE+aW+aN+aS+aP0,(6)aP0=ρcPΔVΔt,b=S˙ΔV+aP0TPn,ΔV=0.5rN+rSΔrΔz.Figure 1
Schematic of mesh generation for fuel pellets or absorbers.


## 3.3. Radiation Heat Transfer

The radiative heat transfer is calculated at the fuel pin level. Due to the shielding effect, central rodP only radiates heat directly to the surrounding 16 rods as shown in Figure 2. According to the geometric position, the surrounding rods can be divided into three types: 1, 2, and 3. The radiative view factors are X1, X2, and X3, respectively. Then, the effective radiation for central rod P is calculated as follows:(7)JPn+1=εσTP4+1−ε∑k=1,4J1,kn+1X1+1−ε∑k=1,4J2,kn+1X2+1−ε∑k=1,8J3,kn+1X3,where J1, J2, and J3 are the effective radiation at the outer surface of rods 1, 2, and 3, respectively. The radiative heat flux at the outer surface of rod P is given by(8)qP=ε1−εσTP4−JPn+1.Figure 2
Schematic of the radiation model between rods.


## 3.4. Cladding Oxidation

The oxygen diffusion model proposed by Olander [16] is adopted to calculate the oxidation of cladding, in which the cladding is divided into the oxide layer (ZrO2) and the metal layer (including a-Zr(O) and β-Zr). The oxidation rate depends on the amount of steam available, the thickness of the metal layer as well as the diffusion rate of oxygen in the cladding (namely temperature). According to the amount of available steam, the steam state in the flow channel can be divided into two cases: steam-rich and steam-starved. The amount of steam available is calculated by the molar steam balance:(9)VchMgdyidt=Fgyi−1−yi−Aclw˙,in which Mg=ptot/RT is the molar density of gas. yi is the steam mole fraction in the bulk gas in node i. The total gas molar flow rate is denoted by Fg and Acl is the external surface area of the cladding. The steam absorption rate is given by(10)w˙a=kgMgyi−ycl,where ycl is the steam molar fraction at the cladding outer surface. The oxygen diffusion rate in the cladding is given by(11)w˙d=ρZrDoCs−Caδ,in which the molar density of Zr in ZrO2 is denoted as ρZr and Do is the diffusion coefficient of oxygen in the oxide layer. Cs, Ca, and Cb are the O/Zr ratio at the outer surface of the oxide layer, at the interface in the oxide layer side, and at the interface in the metal layer side, respectively, which are determined by Zr-O phase diagram. The oxidation rate w˙ is calculated as the minimum of w˙a and w˙d. In the case of oxidation when both oxide and metal phases are present, the thickness of the oxide layer δ is calculated as follows:(12)δn+1=δn+ΔtGCa−Cb−DMCb4E2ξnc¯nCb+GDoδnCs−Ca,where DM is the diffusion coefficient of oxygen in α-Zr(O). The average oxygen content in the metal layer c¯ is determined by the following equation:(13)c¯n+1=c¯n+CbΔtDM2Eξn2c¯Cb−Gξn1−c¯nCbdδdtn+1.For the case of oxidation in all-metal cladding or in all-oxide cladding, other controlling equations are used to calculate the oxide thickness and the average oxygen content in the metal layer.


## 3.5. Dissolution of UO2 and ZrO2

Hofmann [7], Hayward and George [8], and Kim and Olander [9] established dissolution models based on their experiments, respectively. However, the results of each experiment in the precipitate period are significantly different. In the calculation, all the three models as well as the IBRAE model are used to simulate the dissolution of UO2 and ZrO2 and will be verified by the RIAR test in the next section.


### 3.5.1. Kim UO2 Dissolution Model

When the Ra number in the (U, Zr, O) melt is greater than 1000, the natural convection model is adopted [9]:(14)CUn+1=CUS−CUS−CUn−1/4+0.54S4VgDU3νd1/4Δt−4.Otherwise, the diffusion model is used [9]:(15)CUn+1=CUn2+4DUπe2CUS2Δt0.5,in which CU is the mass fraction of U in the bulk melt. S, V, and ν are the UO2, (U, Zr, O) contact surface, the volume, and viscosity of (U, Zr, O), respectively. d is the diameter corresponding to S and e is the (U, Zr, O) mixture thickness.


### 3.5.2. Hayward UO2 Dissolution Model

Based on crucible experiments, Hayward found that the mass fraction of U in the melt can be calculated as follows [8, 17]:(16)CUn+1=CUS−CUS−CUn·exp−kUΔt,in which kU is a constant related to the contact area and the initial oxygen content.


### 3.5.3. Hofmann UO2 and ZrO2 Dissolution Model

The crucible experiments of Hofmann show that UO2 is rapidly dissolved in the initial period, and the mass fraction of UO2 in the melt rapidly increases to 35.8%. Subsequently, UO2 is continuously dissolved, and the UO2 mass fraction in the melt follows the parabolic law [7]:(17)CUO2n+1=0.358+CUO2n−0.3582+Kp1Δt0.5,(18)Kp1=1.02×1011exp−6.77×105RT.The ZrO2 dissolution model is calculated in a similar form of equation (7):(19)CZrO2n+1=0.233+CZrO2n−0.2332+Kp2Δt0.5,(20)Kp2=1.47×1010exp−6.66×105RT.


### 3.5.4. Simultaneous Dissolution of UO2 and ZrO2

For the calculation of the simultaneous dissolution of UO2 and ZrO2, the first method is to assume it as two separate processes. Within one time step, the dissolution of UO2 is calculated firstly according to the models mentioned above, and then the dissolution of ZrO2 with Hofmann correlation is calculated until the mixture composition reaches the liquidus in the (U, Zr, O) ternary phase diagram. The second method is to adopt the IBRAE model proposed by Veshchunov et al. [18, 19], which can be used to calculate the simultaneous dissolution of UO2 and ZrO2 as well as the oxidation of (U, Zr, O) mixture. In the IBRAE model, both intact fuel rods and the molten pool can be considered as UO2, (U, Zr, O) mixture, and ZrO2 three-layer structure. The interface movement rate can be derived from the flux match. For the interface of UO2 and (U, Zr, O),(21)dr1dt=−DOUO2∂ρO/∂rI1−kOI1cOmax−cOBρO−g1+g2cMρU/cM.For the interface of (U, Zr, O) and ZrO2,(22)r3dr3dt−r2dr2dt=r2DOZrO2∂φO/∂rI2−kOI2cOmax−cOBφO−g1ρU/cM,where r1, r2, and r3 are positions of boundaries of the solid UO2 and ZrO2 layers. ρi, ci, and φi are the molar densities of different components (i=O, U, and Zr) in UO2, (U, Zr, O) mixture, and ZrO2, respectively. The molar density of the melt is denoted by cM. g1 and g2 are temperature-dependent parameters which determine the position of liquidus line in the ternary phase diagram.The oxygen concentration in the melt in (21) and (22) can be calculated by the mass conservation equations, while the oxygen concentration gradient at interface I1 is calculated by the diffusion equation of oxygen in solid UO2 with moving boundary conditions. For the oxygen concentration gradient at the interface I2,(23)if steam−rich,DOZrO2∂φ0∂rI2≈DOZrO2φ0B−φ0I1r3−r2,if steam−starved,DOZrO2∂φ0∂rI2≈kgMgyi.


## 3.5.1. Kim UO2 Dissolution Model

When the Ra number in the (U, Zr, O) melt is greater than 1000, the natural convection model is adopted [9]:(14)CUn+1=CUS−CUS−CUn−1/4+0.54S4VgDU3νd1/4Δt−4.Otherwise, the diffusion model is used [9]:(15)CUn+1=CUn2+4DUπe2CUS2Δt0.5,in which CU is the mass fraction of U in the bulk melt. S, V, and ν are the UO2, (U, Zr, O) contact surface, the volume, and viscosity of (U, Zr, O), respectively. d is the diameter corresponding to S and e is the (U, Zr, O) mixture thickness.


## 3.5.2. Hayward UO2 Dissolution Model

Based on crucible experiments, Hayward found that the mass fraction of U in the melt can be calculated as follows [8, 17]:(16)CUn+1=CUS−CUS−CUn·exp−kUΔt,in which kU is a constant related to the contact area and the initial oxygen content.


## 3.5.3. Hofmann UO2 and ZrO2 Dissolution Model

The crucible experiments of Hofmann show that UO2 is rapidly dissolved in the initial period, and the mass fraction of UO2 in the melt rapidly increases to 35.8%. Subsequently, UO2 is continuously dissolved, and the UO2 mass fraction in the melt follows the parabolic law [7]:(17)CUO2n+1=0.358+CUO2n−0.3582+Kp1Δt0.5,(18)Kp1=1.02×1011exp−6.77×105RT.The ZrO2 dissolution model is calculated in a similar form of equation (7):(19)CZrO2n+1=0.233+CZrO2n−0.2332+Kp2Δt0.5,(20)Kp2=1.47×1010exp−6.66×105RT.


## 3.5.4. Simultaneous Dissolution of UO2 and ZrO2

For the calculation of the simultaneous dissolution of UO2 and ZrO2, the first method is to assume it as two separate processes. Within one time step, the dissolution of UO2 is calculated firstly according to the models mentioned above, and then the dissolution of ZrO2 with Hofmann correlation is calculated until the mixture composition reaches the liquidus in the (U, Zr, O) ternary phase diagram. The second method is to adopt the IBRAE model proposed by Veshchunov et al. [18, 19], which can be used to calculate the simultaneous dissolution of UO2 and ZrO2 as well as the oxidation of (U, Zr, O) mixture. In the IBRAE model, both intact fuel rods and the molten pool can be considered as UO2, (U, Zr, O) mixture, and ZrO2 three-layer structure. The interface movement rate can be derived from the flux match. For the interface of UO2 and (U, Zr, O),(21)dr1dt=−DOUO2∂ρO/∂rI1−kOI1cOmax−cOBρO−g1+g2cMρU/cM.For the interface of (U, Zr, O) and ZrO2,(22)r3dr3dt−r2dr2dt=r2DOZrO2∂φO/∂rI2−kOI2cOmax−cOBφO−g1ρU/cM,where r1, r2, and r3 are positions of boundaries of the solid UO2 and ZrO2 layers. ρi, ci, and φi are the molar densities of different components (i=O, U, and Zr) in UO2, (U, Zr, O) mixture, and ZrO2, respectively. The molar density of the melt is denoted by cM. g1 and g2 are temperature-dependent parameters which determine the position of liquidus line in the ternary phase diagram.The oxygen concentration in the melt in (21) and (22) can be calculated by the mass conservation equations, while the oxygen concentration gradient at interface I1 is calculated by the diffusion equation of oxygen in solid UO2 with moving boundary conditions. For the oxygen concentration gradient at the interface I2,(23)if steam−rich,DOZrO2∂φ0∂rI2≈DOZrO2φ0B−φ0I1r3−r2,if steam−starved,DOZrO2∂φ0∂rI2≈kgMgyi.


## 3.6. Criteria of Cladding Oxidation Layer Failure

In the steam-rich case, the oxidation reaction in the cladding dominates. Oxidation of cladding continues until it is completely oxidized. In the calculation, it is assumed that the fuel pellets and cladding will melt together once the temperature reaches the eutectic point (2673 K).In the steam-starved condition, the dissolution reaction in the cladding dominates. In the calculation, it is assumed that the cladding oxide layer will fail once the (U, Zr, O) mixture temperature exceeds theα-Zr(O) melting point (2248 K) and the oxide layer thickness is less than 300 μm.


## 3.7. Corium Relocation

Based on the explanation of Veshchunov and Palagin on CORA W1 and W2 experimental results [15], the corium relocation downward in the bundle is considered as the slug flow with fixed velocity 1 × 10−3m/s. During the relocation, convective heat exchange with the steam, radiant heat exchange with the surrounding nonadjacent rods, and the heat generated by corium oxidation are considered. The energy equation during relocation is given by(24)ρcpdTmpVmpdt=Qox−AmphTmp−Tg−Qrad.


## 4. Separate Effect Validation

In this section, the physical models for both the oxidation reaction and dissolution reaction which are the two of the most important phenomena in the early phase of core degradation will be verified.


### 4.1. Oxidation of Zirconium Cladding

The oxidation model is validated with theoretical cases of a single fuel rod oxidation proposed by Olander [16]. Both steam-rich and steam-starved conditions are simulated. In order to simplify the model and constrain other variables, the effect of heat generation due to oxidation is not considered. Instead, the fuel rod temperature of each node is determined by the following equations:(25)T=T0+RRLt,T<Ttrans,(26)T=Ttrans+RRHt−ttrans,T≥Ttrans,in which T0 is the initial temperature. The temperature at which the microstructure of the oxide layer begins to change from tetragonal to cubic is denoted by Ttrans with a value of 1773 K. RRL and RRH are the temperature ramp rates below and above the transition temperature.In the steam-starved case, the inlet steam flow rate is 2× 10−2 m/s. The temperature ramp rates RRL and RRH are 1.0 K/s and 2.5 K/s, respectively. The steam mole fraction in the flow channel at different positions is illustrated in Figure 3. The cladding oxygen absorption rate w˙ increases as the temperature of the fuel rod rises. Therefore, the steam mole fraction in the flow channel gradually decreases. When the steam-starved condition occurs, the oxygen absorption from the flow channel is limited by the amount of steam available. Once the oxygen absorption rate of the cladding from the flow channel is less than the diffusion rate of oxygen from the oxide layer to the metal layer, the oxide layer will become thinner, as shown in Figure 4.Figure 3
Comparison of steam mole fractions in the flow channel.Figure 4
Comparison of fractional oxide layer thickness in cladding.In the steam-rich case, the inlet steam flow rate is 1 m/s. The temperature ramp ratesRRL and RRH are 0.7 K/s and 1.0 K/s, respectively. The steam mole fraction in the flow channel at different positions is presented in Figure 5. It can be seen that the steam mole fraction in the flow channel decreases steeply at around 600 s. It is explained by the fact that the temperature of the rod increases to around 1773 K at which the microstructure of the oxide layer changes from tetragonal to cubic, causing the diffusion of oxygen in the oxide layer increasing rapidly and large amount of steam consumed. After 750 s, the steam mole fraction in the 210 cm and 390 cm flow channels increase since the cladding below 210 cm has been sufficiently oxidized. Then, the consumption of steam is reduced so that the steam mole fraction arriving above 210 cm increases. After 880 s, the steam mole fraction in the flow channel drops again due to the increase in the temperature of the fuel rod causing the increase of the oxygen diffusion rate in the oxide layer and the metal layer. The variation of fractional oxide layer thickness in cladding is shown in Figure 6.Figure 5
Comparison of steam mole fractions in the flow channel.Figure 6
Comparison of fractional oxide layer thickness in cladding.It can be seen from Figures3–6 that the calculation results in this paper are consistent with the results of Olander, whether it is steam-rich or steam-starved. Then the rationality of the oxidation model in this paper can be verified.


### 4.2. Dissolution of UO2 and ZrO2

RIAR test [20] is used to verify physical models for the simultaneous dissolution of UO2 and ZrO2 and perform the uncertainty analysis.Figures7 and 8 show the U mass fraction in the (U, Zr, O) mixture calculated by different models compared with experimental values and Muller’s simulation results [20]. Since there is no limit on the solubility, the solid UO2 is continuously dissolved by the molten mixture, which causes the U mass fraction in the melt to increase approximately linearly with time in the Hofmann model. But in Kim and Hayward models, the dissolution of UO2 is limited by its solubility in the (U, Zr, O) mixture. It can be inferred from the ternary phase diagram that there is a competition between UO2 dissolution and ZrO2 dissolution. It means that the solubility of UO2 is reduced due to the dissolution of ZrO2. In the simulation with Kim and Hayward models at 2200°C, the dissolution rate of ZrO2 calculated by Hofmann parabolic correlation is overpredicted. Excessive and rapid dissolution of ZrO2 makes the U in the mixture quickly reach saturation, thereby reducing the volume of solid UO2 dissolved. In the IBRAE model, the simultaneous dissolution of ZrO2 and UO2 is characterized by a complete mathematical model instead of the experimental correlations. The results calculated by the IBRAE model in this paper are in agreement with Muller’s results and experimental values.Figure 7
U mass fraction in the melt changes with time at 2100°C.Figure 8
U mass fraction in the melt changes with time at 2200°C.Muller’s results are considered as a reference to perform uncertainty analysis. Dissolved volumes of UO2 and ZrO2 at 2200°C are presented in Figures 9 and 10, respectively. Within 100 s after the occurrence of UO2 dissolution, dissolved volumes of UO2 calculated by each model are approximately consistent, during which the maximum uncertainty range does not exceed 20%. Then, due to the premature saturation of U in the melt, the dissolution of UO2 in Hayward and Kim models is ceased with 35% uncertainty at 200 s. However, because of no limit on solubility and high dissolution rate, the volume of UO2 dissolved using the Hofmann model is maximum and with the uncertainty of 60% at 200 s. As shown in Figure 10, the uncertainty of ZrO2 dissolved volume is close to 300%. The calculated volumes of ZrO2 dissolved by Kim, Hayward, and Hofmann models are considerably different from Muller’s results, which can be explained by the overpredicted dissolution rate of ZrO2 caused by the Hofmann correlation. Therefore, it can be inferred that since the dissolution rate of ZrO2 is high but the dissolution rate of UO2 is low for Hayward and Kim models, it may cause the failure time of the cladding oxide layer to be advanced while the mass of UO2 dissolved at this stage is underestimated. On the contrary, it will overestimate the volume of UO2 dissolved by using the Hofmann model, which is suitable for conservative calculations.Figure 9
Comparison of dissolved UO2 volume for different models at 2200°C.Figure 10
Comparison of dissolved ZrO2 volume for different models at 2200°C.In severe accidents, the dissolution of UO2 by molten zirconium is usually accompanied by cladding oxidation. In the steam-rich case, the zirconium cladding may be completely oxidized to ZrO2. The increase of oxygen content in the (U, Zr, O) mixture will greatly reduce the solubility of U. Therefore, only when the temperatures of UO2 and ZrO2 exceed the melting point or reach the eutectic temperature, UO2 and ZrO2 can relocate to the bottom of the core as corium. But in the steam-starved condition, the zirconium water reaction occurs but is not intense, forming a thin ZrO2 oxide layer on the outside of cladding. In this case, the molten zirconium dissolves UO2 inward and dissolves ZrO2 outward. As a result, some of the fuel pellets melt at the temperature far below its melting point. The dissolved UO2 is relocated downward in the form of (U, Zr, O) mixture, while the undissolved UO2 in the same pellet collapses in the form of debris. Therefore, the temperature corresponding to the meltdown of core is greatly reduced.


## 4.1. Oxidation of Zirconium Cladding

The oxidation model is validated with theoretical cases of a single fuel rod oxidation proposed by Olander [16]. Both steam-rich and steam-starved conditions are simulated. In order to simplify the model and constrain other variables, the effect of heat generation due to oxidation is not considered. Instead, the fuel rod temperature of each node is determined by the following equations:(25)T=T0+RRLt,T<Ttrans,(26)T=Ttrans+RRHt−ttrans,T≥Ttrans,in which T0 is the initial temperature. The temperature at which the microstructure of the oxide layer begins to change from tetragonal to cubic is denoted by Ttrans with a value of 1773 K. RRL and RRH are the temperature ramp rates below and above the transition temperature.In the steam-starved case, the inlet steam flow rate is 2× 10−2 m/s. The temperature ramp rates RRL and RRH are 1.0 K/s and 2.5 K/s, respectively. The steam mole fraction in the flow channel at different positions is illustrated in Figure 3. The cladding oxygen absorption rate w˙ increases as the temperature of the fuel rod rises. Therefore, the steam mole fraction in the flow channel gradually decreases. When the steam-starved condition occurs, the oxygen absorption from the flow channel is limited by the amount of steam available. Once the oxygen absorption rate of the cladding from the flow channel is less than the diffusion rate of oxygen from the oxide layer to the metal layer, the oxide layer will become thinner, as shown in Figure 4.Figure 3
Comparison of steam mole fractions in the flow channel.Figure 4
Comparison of fractional oxide layer thickness in cladding.In the steam-rich case, the inlet steam flow rate is 1 m/s. The temperature ramp ratesRRL and RRH are 0.7 K/s and 1.0 K/s, respectively. The steam mole fraction in the flow channel at different positions is presented in Figure 5. It can be seen that the steam mole fraction in the flow channel decreases steeply at around 600 s. It is explained by the fact that the temperature of the rod increases to around 1773 K at which the microstructure of the oxide layer changes from tetragonal to cubic, causing the diffusion of oxygen in the oxide layer increasing rapidly and large amount of steam consumed. After 750 s, the steam mole fraction in the 210 cm and 390 cm flow channels increase since the cladding below 210 cm has been sufficiently oxidized. Then, the consumption of steam is reduced so that the steam mole fraction arriving above 210 cm increases. After 880 s, the steam mole fraction in the flow channel drops again due to the increase in the temperature of the fuel rod causing the increase of the oxygen diffusion rate in the oxide layer and the metal layer. The variation of fractional oxide layer thickness in cladding is shown in Figure 6.Figure 5
Comparison of steam mole fractions in the flow channel.Figure 6
Comparison of fractional oxide layer thickness in cladding.It can be seen from Figures3–6 that the calculation results in this paper are consistent with the results of Olander, whether it is steam-rich or steam-starved. Then the rationality of the oxidation model in this paper can be verified.


## 4.2. Dissolution of UO2 and ZrO2

RIAR test [20] is used to verify physical models for the simultaneous dissolution of UO2 and ZrO2 and perform the uncertainty analysis.Figures7 and 8 show the U mass fraction in the (U, Zr, O) mixture calculated by different models compared with experimental values and Muller’s simulation results [20]. Since there is no limit on the solubility, the solid UO2 is continuously dissolved by the molten mixture, which causes the U mass fraction in the melt to increase approximately linearly with time in the Hofmann model. But in Kim and Hayward models, the dissolution of UO2 is limited by its solubility in the (U, Zr, O) mixture. It can be inferred from the ternary phase diagram that there is a competition between UO2 dissolution and ZrO2 dissolution. It means that the solubility of UO2 is reduced due to the dissolution of ZrO2. In the simulation with Kim and Hayward models at 2200°C, the dissolution rate of ZrO2 calculated by Hofmann parabolic correlation is overpredicted. Excessive and rapid dissolution of ZrO2 makes the U in the mixture quickly reach saturation, thereby reducing the volume of solid UO2 dissolved. In the IBRAE model, the simultaneous dissolution of ZrO2 and UO2 is characterized by a complete mathematical model instead of the experimental correlations. The results calculated by the IBRAE model in this paper are in agreement with Muller’s results and experimental values.Figure 7
U mass fraction in the melt changes with time at 2100°C.Figure 8
U mass fraction in the melt changes with time at 2200°C.Muller’s results are considered as a reference to perform uncertainty analysis. Dissolved volumes of UO2 and ZrO2 at 2200°C are presented in Figures 9 and 10, respectively. Within 100 s after the occurrence of UO2 dissolution, dissolved volumes of UO2 calculated by each model are approximately consistent, during which the maximum uncertainty range does not exceed 20%. Then, due to the premature saturation of U in the melt, the dissolution of UO2 in Hayward and Kim models is ceased with 35% uncertainty at 200 s. However, because of no limit on solubility and high dissolution rate, the volume of UO2 dissolved using the Hofmann model is maximum and with the uncertainty of 60% at 200 s. As shown in Figure 10, the uncertainty of ZrO2 dissolved volume is close to 300%. The calculated volumes of ZrO2 dissolved by Kim, Hayward, and Hofmann models are considerably different from Muller’s results, which can be explained by the overpredicted dissolution rate of ZrO2 caused by the Hofmann correlation. Therefore, it can be inferred that since the dissolution rate of ZrO2 is high but the dissolution rate of UO2 is low for Hayward and Kim models, it may cause the failure time of the cladding oxide layer to be advanced while the mass of UO2 dissolved at this stage is underestimated. On the contrary, it will overestimate the volume of UO2 dissolved by using the Hofmann model, which is suitable for conservative calculations.Figure 9
Comparison of dissolved UO2 volume for different models at 2200°C.Figure 10
Comparison of dissolved ZrO2 volume for different models at 2200°C.In severe accidents, the dissolution of UO2 by molten zirconium is usually accompanied by cladding oxidation. In the steam-rich case, the zirconium cladding may be completely oxidized to ZrO2. The increase of oxygen content in the (U, Zr, O) mixture will greatly reduce the solubility of U. Therefore, only when the temperatures of UO2 and ZrO2 exceed the melting point or reach the eutectic temperature, UO2 and ZrO2 can relocate to the bottom of the core as corium. But in the steam-starved condition, the zirconium water reaction occurs but is not intense, forming a thin ZrO2 oxide layer on the outside of cladding. In this case, the molten zirconium dissolves UO2 inward and dissolves ZrO2 outward. As a result, some of the fuel pellets melt at the temperature far below its melting point. The dissolved UO2 is relocated downward in the form of (U, Zr, O) mixture, while the undissolved UO2 in the same pellet collapses in the form of debris. Therefore, the temperature corresponding to the meltdown of core is greatly reduced.


## 5. Integral Effect Validation

In-pile experiments Phebus FPT provide a reliable benchmark for numerical simulation of the core degradation process. In this paper, the integral effect of the core early degradation phase is verified by Phebus FPT0 experiment [21, 22]. The input bundle power and steam flow rate are shown in Figure 11. Considering the pressure tube is immersed in circulating cooling water with temperature 438 K, the temperature of the outer surface of the pressure tube is approximately supposed as 438 K.Figure 11
Bundle power and steam flow rate in Phebus FPT0.Figure12 shows the comparison of temperature at 0.6 m in the calibration phase (0∼10000 s) between the simulation results of this paper, ATHLET CD, and ICARE2 as well as experimental values. In this period, due to the increase in power and the decrease in steam flow, the temperature variation of fuel rods can be divided into three rising stages and three platform stages. It can be seen from Figure 13 that the temperature of the fuel rod at 0.7 m increases sharply at around 12000 s. This period is so-called oxidation escalation phase (10000∼13000 s), during which the power of the bundle increases causing the temperature of the bundle rises again. After the cladding temperature exceeds 1500 K, cladding oxidation becomes significant, generating a large amount of heat and hydrogen. Figures 14 and 15 show the hydrogen generation rate and the cumulative hydrogen generation mass, respectively. The hydrogen generated mass calculated in this paper during this period is 90.83 g, which is close to the experimental value 85 g. Figure 16 presents the variation of temperature at 0.6 m in the corium relocation phase (13000 s∼18000 s). The distribution of temperature along the height throughout the simulation is shown in Figure 17.Figure 12
Variation of outer ring fuel rod temperature at 0.6 m.Figure 13
Variation of inner ring fuel rod temperature at 0.7 m.Figure 14
Hydrogen generation rate.Figure 15
Cumulative hydrogen generation mass.Figure 16
Variation of outer ring fuel rod temperature at 0.6 m.Figure 17
Temperature distribution along the height.Figures13 and 16 show that the maximum temperatures at 0.6 m and 0.7 m exceed 2500 K in the oxidation escalation phase, while the maximum temperature at 0.8 m is only 2250 K, as shown in Figure 18. It can be explained that in this phase, massive steam consumed in the middle and lower regions of the bundle causes a short period of steam starved in the upper region, as shown in Figure 19. In the steam-starved case, it is assumed in the calculation that once the temperature of fuel rods exceeds the α-Zr(O) melting point (2248 K) and oxide layer thickness is less than 300 μm, the cladding oxide layer will fail and the previous hold-up (U, Zr, O) melt will relocate downwards. As presented in Figure 17, some dissolved fuel pellets begin to relocate downward after 11970 s, which is quite in agreement with the experimental value 12050 s. Based on the simulation results, it can be inferred that the reason for rapid accident process of LBLOCA may be not resulted from the fuel rapidly rising to its melting point due to lack of cooling, but a large scale of cladding melts, forming a molten metal pool in the flow channel, and dissolving the adjacent fuel pellets, causing the fuel to melt down at much lower temperatures compared with its melting point 3113 K.Figure 18
Variation of inner ring fuel rod temperature at 0.8 m.Figure 19
Variation of steam mole fraction at different heights.


## 6. Conclusions

In this paper, based on the detailed analysis of the mechanism of controlling phenomena in the early phase of core degradation, fine models at the fuel pin level are established in 3D and programed by Fortran to calculate (1) core heating up, (2) zirconium cladding oxidation, (3) dissolution reaction between molten zirconium and fuel pellets, and (4) formation of molten pool in the core active section.Analysis of dissolution models shows that there is a considerable uncertainty in the Hofmann ZrO2 dissolution correlation, which may cause the failure time of the cladding oxide layer to be advanced while the mass of UO2 dissolved is underestimated by using Kim or Hayward models. According to the simulation results of Phebus FPT0, it can be inferred that the dissolution reaction between the molten zirconium and fuel pellets is the main reason for the melting of UO2 at low temperatures. In the steam-starved case, part of the fuel pellets can melt down even at 2248 K and relocate to the bottom of the core, which is much lower than the melting point of UO2.



---

*Source: 6798738-2020-08-03.xml*