Numerical Study on the Corium Pool Heat Transfer With OpenFOAM

In order to ensure the successful implementation of in-vessel retention strategy to terminate molten corium pool evolution, it is necessary to evaluate the heat flux distribution on the lower head wall of reactor pressure vessel. In the present study, the models of internal heating and natural convention buoyancy, as well as the models of WMLES turbulence and phase changing, were applied in the open source CFD software OpenFOAM to perform numerical simulations for the COPRA single-layer molten pool experiment. The distributions of temperature, heat flux, and crust thickness were obtained. The simulation results were in good comparison with COPRA experimental data, proving the validity of the developed model for corium pool heat transfer characteristics. The simulation method and results could be applied to further in-depth study of thermal behavior in the corium pool.


INTRODUCTION
The nuclear reactor core may melt into high-temperature corium and then relocate into the reactor pressure vessel's lower plenum, where the molten material may form a corium pool after a severe reactor accident occurred. In-Vessel Retention (IVR) has been proposed to prevent the deterioration of the severe nuclear accident. The decay heat generated in the corium pool in the lower plenum can be removed by cooling water outside the reactor vessel. The thermal load distributions along the lower head wall were determined by the natural convection heat transfer characteristics inside the corium melt pool, which imposed great significance for ensuring the successful implementation of IVR strategy (Heofanous et al., 1996).
The main components of the relocated corium were non-eutectic mixtures in compositions of metal and metal oxides from fuel rods and brackets. In the early stages of corium pool formation, the lower oxidation fraction may lead to the stratification of melt layers. However, in the late stages, with increasing dissolved oxygen in the corium, the melt pool was in a homogeneous configuration of ceramic mixture according to the Zr-U phase diagram (Asmolov et al., 2004). Therefore, the simulation of homogenous pool remained an important issue with the solidification character of non-eutectic melt .
Heat transfer characteristics of melt in the reactor's lower head were studied experimentally and numerically in past years. The main purpose of these studies was to understand the fundamental behavior of the melt pool inside the reactor vessel . These experiments have been simplified due to the complexity of the actual molten pool and the experimental cost. The facilities can be divided into three categories based on the geometry: quarter-circular slice, semi-circular slice, and hemispherical pool . Famous experiments such as BALI (Bonnet and Seiler, 1999), COPO (Kymäläinen et al., 1992), SIMECO (Sehgal et al., 1999), and COPRA  have provided important data of heat transfer characteristics based on different boundary conditions.
The key points of experiment simulation were the modeling of appropriate turbulence models and solidification models. Horvat and Mavko (2004) performed the simulation of UCLA's experiment (Asfia and Dhir, 1996) using the SST turbulence model and good comparisons were obtained. Dinh and Nourgaliev (1997) pointed out that the k-ε turbulence model could not represent the thermal distributions accurately in the cases with high Rayleigh number turbulence. Fukasawa et al. (2008) compared the turbulence models of the k-ε and large eddy simulation based on the BALI experiment (Bonnet and Seiler, 1999). The results showed that the LES method was capable of describing flow physics and heat transfer characteristics. Tran et al. (Tran et al., 2010) indicated that the implicit LES worked quite well in predicting natural convection heat transfer for strong turbulent pools. Tran and Dinh (2009) developed a phasechange effective convection model (PECM) to simulate the heat transfer characteristics inside molten pool. However, this method failed to analyze the inner temperature and flow field. Zhang et al. (2014) developed the 2D numerical model based on the SIMPLE algorithm and modified k-ε model for the simulation of the partial solidification process with diffusive convection.
In this paper, the simulation work based on LES turbulence model and solidification model was performed for the COPRA experiments in the platform of OpenFOAM. The distributions of thermal parameters were obtained and compared to validate the simulation model.

SOLIDIFICATION MODEL
There were two numerical methods for dealing with phasechange problems. The dynamic mesh technique was able to track the evolution of interfaces but was also computation consuming. Therefore, a traditional grid technology was chosen with theoretical and empirical formula to analyze the process of phase changing. For the melt pool with natural convection and solidification, the governing equations can be established based on the assumption of incompressible fluid.
Mass equation: Momentum equation: where µ eff is the effective dynamic viscosity; S b is the buoyancy force related to density ρ; β is the thermal expansion coefficient; T s is the solidus temperature; S m is the Darcy source term in the mushy region with similar flow resistance of porous medium. Enthalpy energy equation: where λ eff is the effective thermal conductivity and Q is the inner heating density. The solidification will occur if the temperature is lower than the solidus temperature. The enthalpy-porosity model proposed by Voller and Prakash (1987) was used to describe the solidification process inside the molten pool. The Darcy term of the momentum equation was derived from Darcy's law (Darcy, 1856) in equation as follows: where K is the permeability and µ is the fluid viscosity. Due to the non-eutectic characteristics for corium simulant, there existed a mushy zone in front of the curst interface with both solid and liquid phases. It was assumed that the mushy zone was generally characterized with dendrites, which can be treated as porous medium (Flemings, 1974). Then the Darcy equation was obtained as: The Darcy source term in momentum equation was formulated as: where Φ is the liquid fraction; C is the mushy zone constant depending on the morphology of the porous media and recommended here as 1.6 × 10 3 kg/(s·m −3 ) (Voller and Prakash, 1987); ε is a minimum value of 0.001 to prevent the error of division by zero. The porosity in mushy zone was considered as a function of temperature and was expressed in form of error function as (Rösler and Brüggemann, 2011): where T l and T m are the liquidus temperature and average temperature of the mushy zone respectively. Both sensible heat and latent heat of fusion were included in the enthalpy energy equation: Frontiers in Energy Research | www.frontiersin.org where L h is the latent heat of fusion. Then the enthalpy energy equation was rewritten into: Incorporating the models of internal heating and natural convention buoyancy, as well as the models of WMLES turbulence and phase changing into solvers in the open source platform of OpenFOAM, the numerical work was performed for the COPRA single-layer corium pool experiments with strong turbulence.

SIMULATION METHOD
The COPRA single-layer corium pool facility was a twodimensional 1/4 circular slice vessel to simulate the Chinese advanced PWR reactor lower head at full scale. The inner radius of the vessel was 2.2 m and the height was 1.9 m. The curved vessel (30 mm thickness) was enclosed from outside with the cooling path and the vessel's top surface was adiabatic. The molten salt of NaNO 3 -KNO 3 (in mole fraction) compositions was applied as test material. The selected mixture of 20% to 80% has the maximum temperature difference of 60 K between solidus and liquidus line. Dinh et al. (2000) pointed out that the binary mixture of NaNO 3 -KNO 3 has similar characteristics in phase diagrams compared to the real melt. The internal Rayleigh number could reach up to 10 14 -10 16 in COPRA melt pool . Before performing the numerical simulation, some basic assumptions need to be done. The melt was an incompressible Newtonian fluid and the volumetric heat source in the pool was homogenously distributed. The Boussinesq hypothesis was employed, except for the density in the buoyancy term in momentum equation. Other properties are all treated as constants. The geometry of the mesh was the same as the COPRA melt pool with pool height of 1.9 m. The mesh setup was shown in Figure 1. The grid independence was verified by comparisons with 1.2, 1.35, and 1.4 million mesh. The simulation results from 1.35 million were selected with good enough results considering calculation effectiveness. The mesh encryption was performed in the regions near the inner wall with a minimum of 0.6 mm to better interpret the solidification process.
The internal heating density was set to be 10,500 W/m 3 estimated from test data. The initial temperature of the molten pool was set as 600 K and the isothermal boundary of 303.15 K was adopted on the outside curved wall surface. The upper surface of the molten pool was set as radiative with emissivity of 0.44. The other vertical boundary conditions were all set as adiabatic. The transient simulation was required for the large   (Zhang L. T. et al., 2016). It should be noted that the thermal conductivity of the curved wall was modified to consider the additional thermal resistance introduced by the 0.3 mm air gap between the crust and the curved wall.

RESULTS AND DISCUSSION
The temperatures along pool height in steady state were extracted to compare with experimental data in Figure 2. It can be seen that the temperature increasing rate was relatively rapid within a distance of 200 mm above the pool bottom. The weak turbulence flow and thick crust formation was predicted in the bottom region, leading to larger temperature gradient. With pool height increasing, the rate of temperature increase slowed down and temperature distribution tended to be flat. Overall, the trend  of temperature distribution agreed well with the experimental results in the full height regions.
The temperature and velocity distribution in the corium pool from the simulation were presented in Figures 3, 4 respectively. The results clearly illustrated the temperature stratification across the pool region with increasing temperature in the top region. The lower temperatures at the bottom represented the solidified crust along the curved wall. Similarly, the velocity reflecting the turbulence intensity was also higher near the top and side boundaries due to radiation loss and direct cooling. The natural convection inside the corium pool driven by the internal heating gradually transferred heat and energy toward the top region and cooling boundary to achieve its final thermal balance.  The heat flux distribution of the curved wall with the polar angle was compared in Figure 5. The results showed that the heat fluxes from the simulation and experiments were in good comparison at all regions except for the top part. The heat flux decreased near the top surface from experiment because of obvious thermal loss. The calculated heat flux didn't represent this thermal reduction.
When the molten pool solidifies, the crust was formed on the inner wall surface to effectively reduce the heat flux, which was beneficial to the realization of IVR. The calculated crust thickness and experimental data were compared in Figure 6 below 60 • . It was shown that the crust thickness was decreasing with increasing angle, which was in opposite trend of heat flux. It indicated that the simulated crust thickness was obviously lower than the experimental results at polar angles greater than 40 • . In addition, the crust thickness from simulation nearly disappeared over 60 • , resulting in the heat flux continuously increasing. This phenomenon may result from the arrangement of bottom heating rods without enough heating generated near the curved wall in COPRA experiment. The completely homogenous heating from the simulation resulted in little existing crust in the upper region with strong turbulence and higher heat fluxes.

CONCLUSIONS
In order to understand the natural convection heat transfer characteristics of the molten pool under the high Rayleigh number, the models of internal heating and natural convention buoyancy, as well as the models of WMLES turbulence and solidification, were applied in the open source CFD software OpenFOAM to perform numerical simulations for the COPRA single-layer molten pool experiment.
The temperature distributions along pool height and the heat flux distributions along curved wall were all in good comparisons with experimental results. The temperature and velocity distributions across the pool region were obtained to clearly demonstrate the thermal stratification due to natural convection. The calculated thickness distribution of the crust differs over polar angle of 60 • because the internal heating in experiment was not homogenous near wall. Different from the previous simple model of solidification, the present solidification model is specially treated with Darcy source term in this paper. The overall simulation results reasonably reflected the heat transfer characteristics of the molten pool, proving the validity of the developed model for corium pool thermal behavior. The simulation method and results could be applied to further in-depth study of thermal behavior in the corium pool.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
For the enclosed manuscript entitled Numerical study on the corium pool heat transfer with OpenFOAM. The author contributions are listed as follows: ZX contributed the simulation work and also contributed the original writing and revision. HG and YH contributed conceptions. LZ provided experimental data and contributed constructive discussions. ZM, WS, and SB contributed review and suggestions. LP contributed the supervision of the study.