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 (Gaus-Liu and Miassoedov, 2014).

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 (Zhang Y. et al., 2015). 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 (Zhang L. et al., 2015). Famous experiments such as BALI (Bonnet and Seiler, 1999), COPO (Kymäläinen et al., 1992), SIMECO (Sehgal et al., 1999), and COPRA (Zhang Y. et al., 2016) 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 phase-change 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 phase-change 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:

$\begin{array}{ll}\frac{\partial \rho }{\partial t}+\nabla \blacksquare (\rho u)=0& (1)\end{array}$

Momentum equation:

$\begin{array}{ll}\frac{\partial \rho u}{\partial t}+\nabla \blacksquare (\rho uu)={\nabla }^2({\mu }_{eff}u)-\nabla p+S_b+S_m& (2)\end{array}$

$\begin{array}{ll}{S}_b=\rho g\left\{1-max\left[\beta (T-T_s),0\right]\right\}& (3)\end{array}$

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:

$\begin{array}{ll}\frac{\partial (\rho h)}{\partial t}+\blacksquare (\rho uh)=\nabla \blacksquare ({\lambda }_{eff}\nabla T)+Q& (4)\end{array}$

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:

$\begin{array}{ll}u=-\frac{K}{\mu }\nabla p& (5)\end{array}$

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:

$\begin{array}{ll}K=\frac{{\varphi }^3}{C\left(1-{\varphi }^2\right)}& (6)\end{array}$

$\begin{array}{ll}\nabla p=-C{\frac{\mu (1-\varphi )}{{\varphi }^3}}^{2}u& (7)\end{array}$

The Darcy source term in momentum equation was formulated as:

$\begin{array}{ll}{S}_m=-C\frac{{(1-\varphi )}^2}{{\varphi }^3+\epsilon }u& (8)\end{array}$

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³ kg/(s·m⁻³) (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):

$\begin{array}{ll}\varphi =0.5erf\left[\frac{4(T-T_m)}{T_l-T_s}\right]+0.5& (9)\end{array}$

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:

$\begin{array}{ll}h=c_{p}T+\Delta h& (10)\end{array}$

$\begin{array}{ll}\Delta h=\{\begin{array}{cc}{L}_h& T\ge T_l\\ \varphi \blacksquare L_h& T_l\ge T\ge T_s\\ 0& T<T_s\end{array}& (11)\end{array}$

where L_h is the latent heat of fusion. Then the enthalpy energy equation was rewritten into:

$\begin{array}{ll}\frac{\partial \left(\rho c_{p}T\right)}{\partial t}+\nabla \blacksquare (\rho c_{p}uT)=\nabla \blacksquare ({\lambda }_{eff}\nabla T)+S_h+Q& (12)\end{array}$

$\begin{array}{ll}{S}_h=-L_h\left[\frac{\partial (\rho \varphi )}{\partial t}+div(\rho u\varphi )\right]& (13)\end{array}$

$\begin{array}{ll}{S}_h=-\rho L_h\frac{4exp\left[\frac{4(T-T_m)}{T_l-T_s}\right]}{(T_l-T_s)\sqrt{\pi }}\left(\frac{\partial T}{\partial t}+u\nabla t\right)& (14)\end{array}$

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 two-dimensional 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₃-KNO₃ (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₃-KNO₃ has similar characteristics in phase diagrams compared to the real melt. The internal Rayleigh number could reach up to 10¹⁴-10¹⁶ in COPRA melt pool (Zhang Y. et al., 2016).

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.

FIGURE 1

Figure 1. Mesh setup of the COPRA.

The internal heating density was set to be 10,500 W/m³ 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 eddy simulation method but only steady results were compared with experimental data.

Some physical property parameters for the COPRA salt simulant were listed in Table 1 (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.

TABLE 1

Table 1. Main properties of COPRA experiment.

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.

FIGURE 2

Figure 2. Comparison of experimental and simulated temperature distribution.

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.

FIGURE 3

Figure 3. Temperature distribution in the corium pool.

FIGURE 4

Figure 4. Velocity distribution in the corium pool.

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.

FIGURE 5

Figure 5. Comparison of experimental and simulated heat flux distribution.

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.

FIGURE 6

Figure 6. Comparison of experimental and simulated thickness.

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.

Funding

This work is financially supported by the National Natural Science Foundation of China (Grant Nos. 11705188 and 11805026) and the China Postdoctoral Science Foundation funded project (Grant No. 2019T120811).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Abbreviations

CFD, Computational Fluid Dynamics; LES, Large Eddy Simulation; IVR, In-Vessel Retention; OpenFOAM, Open Field Operation and Manipulation; WMLES, Wall-Modeled LES model.

References

Asfia, F. J., and Dhir, V. K. (1996). An experimental study of natural convection in a volumetrically heated spherical pool bounded on top with a rigid wall. Nucl. Eng. Des. 163, 333–348. doi: 10.1016/0029-5493(96)01215-0

CrossRef Full Text | Google Scholar

Asmolov, V. G., Bechta, S. V., Khabensky, V. B., Gusarov, V. V., Vishnevsky, V. Y., Kotova, S. Y., et al. (2004). “Partitioning of U, Zr and Fe between molten oxidic and metallic corium,” in Proceeding of MASCA Seminar (Aix-en-Provence), 1–13.

Google Scholar

Bonnet, J. M., and Seiler, J. M. (1999). “Thermal hydraulic phenomena in corium pools: the bali experiment,” in the 7th International Conference on Nuclear Engineering (ICONE-7) (Tokyo: ASME).

Google Scholar

Darcy, H. (1856). Les Fontaines Publiques De la ville De Dijon (Paris: Victor Dalmont).

PubMed Abstract | Google Scholar

Dinh, T. N., Konovalikhin, M. J., and Sehgal, B. R. (2000). Core melt spreading on a reactor containment floor. Prog. Nucl. Energy 36, 405–468. doi: 10.1016/S0149-1970(00)00088-3

CrossRef Full Text | Google Scholar

Dinh, T. N., and Nourgaliev, R. R. (1997). Turbulence modelling for large volumetrically heated liquid pools. Nucl. Eng. Des. 169, 131–150 doi: 10.1016/S0029-5493(96)01281-2

CrossRef Full Text | Google Scholar

Flemings, M. C. (1974). Solidification Processing. Metall. Trans. 5, 2121–2134. doi: 10.1007/BF02643923

CrossRef Full Text | Google Scholar

Fukasawa, M., Hayakawa, S., and Saito, M. (2008). Thermal-hydraulic analysis for inversely stratified molten corium in lower vessel. J. Nucl. Sci. Technol. 45, 873–888. doi: 10.1080/18811248.2008.9711489

CrossRef Full Text | Google Scholar

Gaus-Liu, X., and Miassoedov, A. (2014). “Live experimental results of melt pool behaviour in the pwr lower head with insulated upper lid and external cooling,” in Proceedings of the 21st International Conference on Nuclear Engineering 2013 (Chengdu: ASME). doi: 10.1115/ICONE21-15204

CrossRef Full Text | Google Scholar

Heofanous, T. G., Liu, C., Additon, S., Angelini, S., Kymäläinen, O., and Salmassi, T. (1996). In-vessel Coolability and Retention of Core Melt. 1:DOE/ID-10460. doi: 10.2172/491623

CrossRef Full Text | Google Scholar

Horvat, A., and Mavko, B. (2004). “Numerical investigation of natural convection heat transfer in volumetrically heated spherical segments,” in Conference of Nuclear Energy for New Europe (Portorož).

Google Scholar

Kymäläinen, O., Tuomisto, H., and Theofanous, T. G. (1992). “Critical heat flux on thick walls of large, naturally convecting loops,” in ANS Proc 1992 National Heat Transfer Conference (San Diego, CA: HTD) 6, 44–50.

Google Scholar

Rösler, F., and Brüggemann, D. (2011). Shell- and- tube type latent heat thermal energy storage: Numerical analysis and comparison with experiments. Heat Mass Transfer 47, 1 027–1 033. doi: 10.1007/s00231-011-0866-9

CrossRef Full Text | Google Scholar

Sehgal, B. R., Nourgaliev, R. R., Dinh, T. N., Karbojian, A., Green, J. A., and Bui, V. A. (1999). “FOREVER experiments on thermal and mechanical behavior of a reactor pressure vessel during a severe accident,” in Conference: Workshop on In-Vessel Core Debris Retention and Coolability (Garching).

Google Scholar

Tran, C. T., and Dinh, T.-N. (2009). The effective convectivity model for simulation of melt pool heat transfer in a light water reactor pressure vessel lower head. Part II: model assessment and application. Prog. Nucl. Energy 51, 860–871. doi: 10.1016/j.pnucene.2009.06.001

CrossRef Full Text | Google Scholar

Tran, C. T., Kudinov, P., and Dinh, T. N. (2010). An approach to numerical simulation and analysis of molten corium coolability in a boiling water reactor lower head. Nucl. Eng. Des. 240, 2148–2159. doi: 10.1016/j.nucengdes.2009.11.029

CrossRef Full Text | Google Scholar

Voller, V. R., and Prakash, C. (1987). A fixed-grid numerical modeling methodology for convection-diffusion mushy region phase- change problems. Heat Mass Transfer 30, l709–1720. doi: 10.1016/0017-9310(87)90317-6

CrossRef Full Text | Google Scholar

Zhang, L., Zhou, Y., Zhang, Y., Tian, W., Qiu, S., and Su, G. (2015). Natural convection heat transfer in corium pools: a review work of experimental studies. Prog. Nucl. Energy 79, 167–181. doi: 10.1016/j.pnucene.2014.11.021

CrossRef Full Text | Google Scholar

Zhang, L. T., Zhang, Y. P., Zhou, Y. K., Su, G. H., Tian, W. X., and Qiu, S. Z. (2016). COPRA experiments on natural convection heat transfer in a volumetrically heated slice pool with high rayleigh numbers. Ann. of Nucl. Energy 87, 81–88. doi: 10.1016/j.anucene.2015.08.021

CrossRef Full Text | Google Scholar

Zhang, Y., Niu, S., Zhang, L., Qiu, S., Su, G., and Tian, W. (2015). A review on analysis of LWR severe accident. J. Nucl. Eng. Radiat. Sci. 1:041018. doi: 10.1115/1.4030364

CrossRef Full Text | Google Scholar

Zhang, Y., Zhang, L., Zhou, Y., Tian, W., Qiu, S., Su, G., Zhao, B., Yuan, Y., and Ma, R. (2016). The COPRA experiments on the in-vessel melt pool behavior in the RPV lower head. Ann. Nucl. Energy 89, 19–27. doi: 10.1016/j.anucene.2015.11.013

CrossRef Full Text | Google Scholar

Zhang, Y. P., Su, G. H., Qiu, S. Z., Tian, W. X., Gaus-Liu, X., Kretzschmar, F., and Miassoedov, A. (2014). Numerical study on the heat transfer characteristics of LIVE-L4 melt pool with a partial solidification process. Prog. Nucl. Energy 74, 213–221. doi: 10.1016/j.pnucene.2014.03.011

CrossRef Full Text | Google Scholar

Nomenclature