The Development of a Multiphysics Coupled Solver for Studying the Effect of Dynamic Heterogeneous Configuration on Particulate Debris Bed Criticality and Cooling Characteristics

Abstract

For a sodium-cooled fast reactor, the capability for stable cooling and avoiding re-criticality on the debris bed is essential for achieving in-vessel retention when severe accidents occur. However, an unexploited uncertainty still existed regarding the compound effect of the heterogeneous configuration and dynamic particle redistribution for the debris bed’s criticality and cooling safety assessment. Therefore, this research aims to develop a numerical tool for investigating the effects of the different transformations of the heterogeneous configurations on the debris bed’s criticality/cooling assessment. Based on the newly proposed methodology in this research, via integrating the Discrete Element Method (DEM) with Computational Fluid Dynamics (CFD) and Monte-Carlo-based Neutronics (MCN), the coupled CFD–DEM–MCN solver was constructed with the originally created interface to integrate two existing codes. The effects of the different bed configurations’ transformations on the bed safety assessments were also quantitively confirmed, indicating that the effect of the particle-centralized fissile material had the dominant negative effect on the safety margin of avoiding re-criticality and particle re-melting accidents and had a more evident impact than the net bed-centralized effect. This coupled solver can serve to further assess the debris bed’s safety via a multi-physics simulation approach, leading to safer SFR design concepts.

1. Introduction

It is essential to ensure plant safety in a sodium-cooled fast reactor (SFR), a next-generation nuclear energy system for the future proposed by the Generation IV International Forum. When a Severe Accident (SA) occurs, the Core Disruptive Accident (CDA) is one of the significant safety concerns of SFRs. According to the progression status of the postulated CDA, the developed sequence can be categorized into four phases: the Initiating phase, the Transition phase, the Material relocation phase, and the Heat removal phase [1]. During the material relocation phase, the molten core material may be discharged downward to the lower sodium plenum through the potential path (e.g., Control Rod Guide Tube, CRGT [2,3]), and the molten material can be quenched and broken into particulate fragments due to the hydrodynamic instability and the thermal interaction with the surrounding sodium coolant under Fuel–Coolant Interaction (FCI) [4,5]. Eventually, these fragmented particles (debris particles) accumulate and are deposited on the in-vessel debris catcher in the lower plenum to form the particulate debris bed. Based on previous studies, it can be recognized that a significant power excursion due to a re-criticality accident of the degraded molten materials would damage the integrity of the Reactor Vessel (RV), and the insufficient cooling of the generated decay heat from the debris bed to the coolant would also cause the coolant to dry out and particle re-melting following further melting through the RV boundary [1,2,3,6]. In order to secure the safety margin of In-Vessel Retention (IVR), ensuring the avoidance of sudden power excursions and insufficient cooling of the debris bed are the prerequisite concerns. Therefore, investigating the factors that may affect the criticality potential and the cooling characteristics of the heat-generating particulate debris bed are crucial for realizing the goal of IVR against CDA to escalate the plant safety of SFRs.

Some studies regarding the debris bed’s safety analysis have been conducted on the aspects of the criticality and cooling of the generated-decay heat. As for the issues of the debris bed criticality, the Japan Atomic Energy Agency (JAEA) conducted fuel debris criticality experiments by addressing the scenario of Molten-Core-Concrete-Interaction (MCCI) in the STAtic experiment Critical facility (STACY) [7]. Two configurations of simulated fuel debris (under-moderated and over-moderated) were prepared by varying the interval between fuels. The results showed that the influence of the concrete composition on the effective multiplication factor ($k_{eff}$ [-]) was manifested more in under-moderated configuration than over-moderated ones. Freiría López et al. [8] conducted numerical criticality studies on a virtual spherical debris bed under an oxide fuel (UO₂)/water system. They indicated that the factors of the burnup of the fuel material, particle size, porosity, and debris bed size affect the necessary boron concentration for securing the sub-critical state. Li et al. [9] numerically investigated the criticality safety of the conical debris bed under a Mixed Oxide (MOX) fuel/liquid system. Except for the effect of the debris bed’s porosity, this also demonstrated that, based on the assumption of the same volume, the lower slope angle of the conical debris bed contributed a negative effect on the $k_{eff}$ value of the debris bed. As for the issue of the cooling characteristics of the heat-generating debris bed, a criterion of the Dryout Heat Flux (DHF), which is the maximum heat flux prior to the incipient of the possible dryout, was selected to conduct cooling capability investigations on a particular debris bed’s configuration in previous studies. For a one-dimensional debris bed, early researchers attempted to numerically predict the DHF value for a specific configuration of the heat-generating particulate bed, and several models were proposed based on the analytical or empirical corrections [10,11,12]. These models can describe how the debris bed’s properties, such as particle size, porosity, and bed thickness, affected the DHF value of the heat-generating porous bed. For the more plausible scenario of the multi-dimensional cases, the cooling characteristics between the top-flooded cylindrical debris bed and the heap-like conical one with multi-dimensional flooding were investigated in the COOLOCE experimental facility [13]. The results demonstrated that the multi-dimensional cooling of the conical one resulted in a higher DHF value than the top-flooded cylindrical one of the same bed heights, but the top-flooded cylindrical with lower debris heights could have a higher DHF value if two debris beds are assumed to have equal volume and bottom radius.

Thus far, most previous studies on the proposed empirical models for investigating the effects of the debris bed properties on the criticality potential and the cooling characteristics are still restricted to the assumption that the debris bed is composed of homogeneous solid particles. However, based on previous studies, the debris bed in the actual scenario would be composed of a mixture of heterogeneous particles with a non-uniformity of the material composition (such as fuel particles and structure material particles [14,15,16]), accompanying a shape and size distribution [4,5,17,18,19,20] due to the stochastic process of the molten material’s discharging and fragmentation. Hence, there would exist a debris bed with a heterogeneous configuration created by the mixture of heterogeneous solid particles in the realistic SA scenario, which would address the uncertainty source of the bed safety assessments from previous studies. In addition, during the bed formation, the heterogeneous configuration of the debris bed would process further dynamic transformation due to the particle redistribution, such as particle avalanching (when the slope angle of the debris bed exceeded the critical repose angle during bed formation [21,22,23]) or possible self-leveling (due to the mechanical energy of the evaporated coolant within the particulate debris bed [24,25,26,27,28]), causing the partial fluidization of the top layer and the subsequent flattening of the bed via lateral particle redistribution. For example, considering the heterogeneous debris bed composed of the fuel and structure particles (as shown in Figure 1) under the process of particle redistribution, except for the thickness of the debris bed being suppressed, the coherent phenomenon of gradual local accumulation of the heavier fuel particles due to the density difference will also occur [14,15,16,29], with the accompaniment of the local higher neutron/heat flux in the debris bed. Although previous studies have implied that the suppressed bed thickness could escalate the cooling capability and be beneficial to avoid criticality accidents [9,13], whether the dynamic transformation of the heterogeneous configuration would increase or decreases the $k_{eff}$ value or whether it promotes or inhibits the flow within the porous bed has not been quantitatively elucidated as yet. To the authors’ best knowledge, no published document that registers the compound effect of the heterogeneous configuration and dynamic particle redistribution on the bed safety assessment has been made available. Therefore, this research aims to develop a tool for investigating the effects of the different dynamic transformations of the heterogeneous configurations on both the criticality potential and cooling characteristics of the debris bed, which is beneficial for further securing the sufficient safety margin of IVR through eliminating the source of uncertainties of this compound effect from the bed safety assessment.

Even for the fixed-particle condition, there are still existing complex multi-mechanisms when considering bed safety assessment. In addition to the mechanism of neutronics for evaluating the criticality safety, the thermal hydraulics in the fluid phase, the heat transfer between the particles, and the momentum/energy transfer between particle-fluid phases also should also be taken into account for conducting the cooling characteristics assessment of a decay heat-generating porous bed, as shown Figure 1. Moreover, the importance of the multi-physics interplay is further increased when considering dynamic particle redistribution into the debris bed’s safety assessment. For example, dynamic particle redistribution can address the extra effect on the flow velocity distribution and following heat convection capability, which can further affect the temperature distribution of the solid particles and coolant fluid. On the other hand, the altered flow velocity distribution can also address the effects of the subsequent particle movement by the solid-fluid interfacial force, subsequently contributing to the heterogeneous configuration transformation, whose effect will be further reflected in the neutronics calculation.

The transience with complex multiphysics interplay phenomena makes it hard to measure all the required timely data for criticality/cooling analysis simultaneously via experimental approaches; therefore, numerical analysis could be a suitable alternative to capture the main physical features along the transient time domain and yield the assessments in a broad spectrum, which was selected as the main direction of this study. However, the typical numerical method used to conduct thermal-hydraulic analysis in the fluid phase, Computational Fluid Dynamics (CFD), which assigns the area of the particulate bed as the homogeneous porous media by applying the empirical pressure drop model [30], has the inherent limitation that the particles are assumed to have stationary status, and Monte-Carlo-based Neutronics (MCN, to conduct neutronic analysis) are under the same limitation. This can also cause the effect of the multi-physics interplay accompanied by particle redistribution being unable to be integrated into the bed safety assessment. Hence, in order to realize the purpose of this research when considering the interaction of multi-physics, a newly coupled methodology is proposed here by incorporating the Discrete Element Method (DEM) with CFD and MCN methods, which is free from the constraint of the static particle conditions and even accessible for the scenario of the debris bed’s heterogenous configuration transformation under particle redistribution when considering the multi-physics interaction between the criticality/cooling analysis. The characteristics of the coupled elements are explained below:

• DEM, a Lagrange-based analytical method developed by Cundall and Strack for considering multi-body collision and movement [31], is suitable for simulating particle redistribution. In addition, since DEM can introduce individually heterogeneous solid DEM particles into the computational domain, it can reproduce various heterogeneous configurations, and should be allowed to consider the heat conduction transfer between solid particles;

• CFD is the Eulerian-based method for conducting thermal-hydraulic analysis in the fluid phase. Additionally, the momentum exchange (e.g., drag force, buoyance) and energy exchange (e.g., heat convection transfer) between fluid and solid particle can also be estimated using sub-models based on the calculated information from both CFD the DEM sides, which also contributes to the further effect on the subsequent heterogeneous configuration’s transformation, predicted in the DEM side;

• Because of its flexible geometry description capability, MCN can construct highly sophisticated and heterogeneous geometries in detail, and it is possible to reflect information from both the CFD and DEM sides to perform neutronic evaluations for the debris bed’s criticality analysis.

In summary, through the newly proposed coupled CFD–DEM–MCN methodology described in this research, the transient data of the fluid phase information, particle movements with temperature, and the $k_{eff}$ eigenvalue of the bed system can be acquired during the dynamic process of the heterogeneous bed configuration’s transformation, considering the multiphysics interplay.

For automatically conducting the criticality/cooling bed safety assessment, the tool of the coupled CFD–DEM–MCN solver was also newly built in this research by integrating existing codes. The CFD solver, DEM solvers, and the sub-models for evaluating momentum/energy exchange between fluid-particle phases were executed via the commercial code of STAR-CCM+ (Version 2020.3) [32], a general-purpose tool allowing the processing of three-dimensional calculation among multi-phases with high fidelity and flexibility. The MCN solver can be conducted on the in-house code of MVP3, developed by the JAEA [33] and validated via conducting benchmarking on the debris particles system with two other Monte-Carlo-based codes (Serpent and MCU), showing that the three modern codes are capable of performing neutronic analysis for the fuel debris particles system without significant discrepancy and with high accuracy [34]. The interface, which is liable for exchanging information and the workflow control between STAR-CCM+ and MVP3, was originally created in this research to realize the automated multi-physics-concerning function in the coupled CFD–DEM–MCN solver. In order to examine the capability of the coupled CFD–DEM–MCN solver, numerical investigations of the different heterogeneous bed configurations’ transformations under the SFR’s postulated SA scenarios were carried out, and the effects of these alternative processes on the criticality/cooling bed safety assessments were also quantitatively confirmed in this research. The following content will further explain the detailed structure and the examined results of the CFD–DEM–MCN coupled solver.

2. Analytical Methodology

This chapter introduces the detailed structure of the CFD–DEM–MCN coupled solver, aiming to apply to the criticality/cooling bed safety assessment with the availability of considering multiphysics interplay under the heterogeneous bed configuration’s transformation due to particle redistribution in the SFR’s postulated SA scenarios. CFD–DEM (introduced in Section 2.1), executed on the STAR-CCM+ code, was in charge of conducting the cooling characteristics analysis of the bed safety assessment and evaluating the transformation process of the heterogeneous bed configuration, in which the analysis in the continuous fluid phase and discrete particle phase was through CFD and DEM solvers, respectively. As for the momentum and energy exchanges between phases, the sub-models provided in the STAR-CCM+ were utilized in this research to realize the physical interplay between fluid and moving heterogeneous particles. Following the explanation of the MCN solver (executed on the code of MVP3) in Section 2.2, interface for integrating STAR-CCM+ and MVP3 created initially is introduced in Section 2.3, with the explanations of the functions of workflow control and the data exchange between two codes for accomplishing the physics interplay of the heterogeneous bed configuration’s transformation and neutronic calculation.

2.1. CFD–DEM

2.1.1. Continuous Fluid Phase

The thermal-hydraulics analysis in the continuous fluid phase of this research was conducted by the CFD solver provided from the code of STAR-CCM+, which stemmed from the fundamental core set of the governing equations. For the fluid-particle multi-phase system, the fluid was not fully occupied in a CFD cell due to the discrete particles existing in the fluid region. Hence, the CFD governing equations applied in this research consider the fluid volume occupied in a CFD cell in order to introduce the fluid volume fraction (${\epsilon }_f$ [-]) of a CFD cell, defined as below:

(1)${\epsilon }_f=\frac{V_f}{V_{cell}},$

where $V_{cell}$ [m³] and $V_f$ [m³] are the CFD cell and fluid volumes in a CFD cell, respectively. The governing equations of the CFD solver, considering ${\epsilon }_f$, are listed below, including the conservation equations for mass (continuity equation, Equation (2)), momentum (Navier–Stokes equation, Equation (3)), and energy (Equation (4)), respectively [32]:

(2)$\frac{\mathsf{\partial }{\epsilon }_f{\rho }_f}{\mathsf{\partial }t}+\nabla \cdot ({\epsilon }_f{\rho }_f{\overrightarrow{V}}_f)=0,$

(3)$\frac{\mathsf{\partial }({\epsilon }_f{\rho }_f{\overrightarrow{V}}_f)}{\mathsf{\partial }t}+\nabla \cdot ({\epsilon }_f{\rho }_f{\overrightarrow{V}}_f{\overrightarrow{V}}_f)=-{\epsilon }_f(\nabla p+{\rho }_f\overrightarrow{g})+\nabla \cdot ({\mu }_f{\epsilon }_f\nabla {\overrightarrow{V}}_f)+{\overrightarrow{S}}_m,$

(4)$\frac{\mathsf{\partial }{\epsilon }_f{\rho }_f{C}_{p,f}{T}_f}{\mathsf{\partial }t}+\nabla \cdot ({\epsilon }_f{\rho }_f{C}_{p,f}{\overrightarrow{V}}_f{T}_f)=\nabla \cdot (k_f{\epsilon }_f\nabla T_f)+S_E,$

where ${\rho }_f$ [kg/m³], ${\overrightarrow{V}}_f$ [m/s], and ${\mu }_f$ [kg/m·s] are the fluid density, velocity, and dynamic viscosity, respectively. $p$ [kg/m·s²] and $\overrightarrow{g}$ [m/s²] are the pressure and gravity. $C_{p,f}$ [J/kg·K], $k_f$ [W/m·K], and $T_f$ [K] stand for the heat capacity, thermal conductivity, and the temperature of the fluid. ${\overrightarrow{S}}_m$ [N/m³] and $S_E$ [W/m³] are presented as the momentum and energy source terms contributed by the feedback from the discrete particle phase. The details of ${\overrightarrow{S}}_m$ and $S_E$ will be further explained in the latter Section 2.1.3.

2.1.2. Discrete Particle Phase

The motion of the discrete particle phase was evaluated by the DEM solver provided by the STAR-CCM+ code in this research. In the DEM solver, particle translation and rotation are considered based on Newton’s law of motion, with no mesh requirement.

The governing equation for simulating the particle translation is subjected to the particle linear momentum conservation, given below [32]:

(5)$m_i\frac{d{\overrightarrow{V}}_i}{dt}=m_i\overrightarrow{g}+{\displaystyle \sum _{j=1,j\ne i}({\overrightarrow{F}}_{cont,ij})}+{\overrightarrow{F}}_{f,i},$

where $m_i$ [kg] and ${\overrightarrow{V}}_i$ [m/s] are the mass and the velocity of the particle i. ${\overrightarrow{F}}_{f,i}$ [kg·m/s²] is the fluid-particle force on the particle i, and is explained in Section 2.1.3. ${\overrightarrow{F}}_{cont,ij}$ is the contact force from particle j to particle i, evaluated by the contact model (spring-dashpot model), referring to the particle properties (e.g., Young’s modulus, Poisson’s ratio, static friction coefficient, rolling resistance coefficient, and particle size). The detailed explanations of the spring-dashpot model applied in this research can be found in our previous research (regarding the CFD–DEM method’s validation, which is also conducted based on the code of STAR-CCM+ [35]).

The equation of the energy conservation for the discrete particle phase can be expressed as below [32]:

(6)$m_i{C}_{p,i}\frac{d{T}_i}{dt}={\displaystyle \sum _{j=1,j\ne i}^{}(Q_{cond,ij})+Q_{s,i}+Q_{f,i}},$

where $C_{p,i}$ [J/kg·K] and $T_i$ [K] are the particle i’s specific heat and temperature. $Q_{f,i}$ [W] is the heat transfer term from the fluid phase to the particle i, which will be explained in Section 2.1.3. $Q_{s,i}$ [W] is the heat source term of the particle i, such as the generated decay heat in the debris particle. $Q_{cond,ij}$ [W] represents the heat conduction from the particle j to i, and it can be obtained as below [32]:

(7)$Q_{cond,ij}=4{\gamma }_c{k}_{eq}(T_j-T_i),$

where ${\gamma }_c$ [m] is the contact area radius. $k_{eq}$ [W/m·K] is the equivalent thermal conductivity between the particle i and j in the form of:

(8)$\frac{1}{k_{eq}}=\frac{1}{k_i}+\frac{1}{k_j},$

where $k_i$ [W/m·K] and $k_j$ [W/m·K] are the thermal conductivities of particle i and j, respectively.

2.1.3. Momentum and Energy Exchange between Phases

This section aims to explain the momentum/energy exchange source terms in the continuous fluid phase (${\overrightarrow{S}}_m$ and $S_E$) and discrete particle phase (${\overrightarrow{F}}_{f,i}$ and $Q_{f,i}$). It should be noted in advance that the momentum/energy exchange between fluid and particle phases in this research is based on the conception of the two-way coupling, which means that the momentum/energy transferred to a specific particle from the fluid phase of the CFD mesh where the particle centroid is located will provide the equal and opposite feedback to the fluid phase of this CFD mesh [32]. In addition, the CFD mesh coupling with DEM simulation is always unresolved in STAR-CCM+, which means that the momentum/energy exchanges of a CFD mesh are calculated by collecting the feedback of all the particles whose centroids are located at this CFD mesh, without solving the detail flow redirection around the solid particles. In other words, the shape of the solid particle is unresolved in the CFD mesh [32,36].

Firstly, the forces exerted from the fluid phase to the particle i considered can be expressed as below:

(9)${\overrightarrow{F}}_{f,i}={\overrightarrow{F}}_{p,i}+{\overrightarrow{F}}_{d,i},$

where ${\overrightarrow{F}}_{p,i}$ [kg·m/s²] and ${\overrightarrow{F}}_{d,i}$ [kg·m/s²] are pressure gradient and viscous drag forces, respectively. ${\overrightarrow{F}}_{p,i}$ can be counted as the buoyance, and formulated as below [32]:

(10)${\overrightarrow{F}}_{p,i}=-V_i\nabla p,$

where $V_i$ [m³] is the volume of the particle i. ${\overrightarrow{F}}_{d,i}$ is due to the relative velocity between particle and viscous fluid, as follows [32]:

(11)${\overrightarrow{F}}_{d,i}=\frac{1}{2}{C}_{d,i}{\rho }_f{A}_i\left|{\overrightarrow{V}}_f-{\overrightarrow{V}}_i\right|({\overrightarrow{V}}_f-{\overrightarrow{V}}_i),$

where $A_i$ [m²] is the projected area of the particle i. $C_{d,i}$ [-] is the drag coefficient of particle i.

In our previous study, the different drag coefficient models provided by the STAR-CCM+ were verified and validated with experiments which aimed to simulate the binary-density particle redistribution under the consideration of the postulated SA scenario of the prototype SFR [35]. In the experiment, a binary-density particle bed, composed of lighter-density Fluorinated Ethylene Propylene (FEP, density = 2000 kg/m³) and a heavier-density stainless steel sphere (density = 7870 kg/m³), was settled in a transparent container, as in Figure 2a. The initial position of the stainless steel particle was located on the center top of the particle bed. Then, the water was injected from the bottom to observe the movement of the binary-density particle bed. Since the refractive index of the FEP particle was close to the water at visible wavelengths, it appears almost transparent in the water, and the historical data of the stainless steel particle’s movement can be tracked even in the middle of the FEP particles. Subsequently, the CFD–DEM model was created. Two types of simulating particles were placed into the flow region, as shown in Figure 2b. In both simulations and experiments, the historical data of the heavier steel particle’s movement were transformed via Fast Fourier Transform (FFT) to observe the particle movement pattern on the frequency domain, as shown in Figure 2c. By comparing the data between the STAR-CCM+’s simulation and the experiments, the validated results reveal that the Gidaspow drag coefficient model yields a modest error deviation, which is applied in this research. Further details of the model validation and the Gidaspow drag coefficient model were explained in our previous study [35].

The heat transferred from the fluid phase to the particle i can be expressed as the heat convection term ($Q_{conv}{}_{,i}$ [W]),

(12)$Q_{f,i}=Q_{conv}{}_{,i},$

and $Q_{conv}{}_{,i}$ can be formulated as [32]:

(13)$Q_{conv}{}_{,i}=h_i{A}_{s,i}(T_f-T_i),$

where $h_i$ [W/m²·K] and $A_{s,i}$ [m²] are the heat convection coefficient and particle surface area of particle i. In the numerical simulation, the value of $h_i$ is granted in terms of the particle i’s Nusselt number (${\mathrm{Nu}}_{p,i}$ [-]) via the form below:

(14)${\mathrm{Nu}}_{p,i}\equiv \frac{h_i{D}_i}{k_f},$

where $D_i$ [m] is the particle i’s diameter. For estimating ${\mathrm{Nu}}_{p,i}$, the Ranz–Marshall correlation provided by STAR-CCM+ was utilized here, which was suitable for the spherical particle with an applicable ${\mathrm{Re}}_{p,i}$ range to apply to the hypothetical SA scenarios in the prototype SFR [26,32,37]. The Ranz–Marshall correlation is formulated as below [32,37]:

(15)${\mathrm{Nu}}_{p,i}=2(1+0.3{\mathrm{Re}}_{p,i}^{1/2}{\mathrm{Pr}}_f^{1/3}),$

where ${\mathrm{Pr}}_f$ [-] and ${\mathrm{Re}}_{p,i}$ [-] are the Prandtl number of the fluid phase and the particle Reynolds number of the particle i. The definitions of ${\mathrm{Pr}}_f$ and ${\mathrm{Re}}_{p,i}$ are as below:

(16)${\mathrm{Pr}}_f\equiv \frac{C_{p,f}{\mu }_f}{k_f},$

(17)${\mathrm{Re}}_{p,i}=\frac{{\rho }_f\left|{\overrightarrow{V}}_f-{\overrightarrow{V}}_i\right|D_i}{{\mu }_f},$

Finally, after solving ${\overrightarrow{F}}_{f,i}$ and $Q_{f,i}$, the momentum/energy exchanges from the particles to the fluid phase of a CFD mesh are based on the feedback from the total particles whose centroids are located at the corresponding CFD mesh, as below [32]:

(18)${\overrightarrow{S}}_m\equiv \frac{-{\displaystyle \sum _{i=1}^{}{\overrightarrow{F}}_{f,i}}}{V_{cell}},$

(19)$S_E\equiv \frac{-{\displaystyle \sum _{i=1}^{}{Q}_{f,i}}}{V_{cell}},$

2.2. MCN

The MCN method is utilized for simulating the neutron generations approximately using the discrete batches. In each batch, the computationally pragmatic numbers of neutrons are initialized. Each neutron’s history from birth (because of the external source or a fission event) to death (due to absorption or leakage out of the calculation domain) is simulated independently. The user can state the number of batches and the number of histories in a batch. Based on the Monte Carlo spirit, sampling via employing the proper probability density functions (cross sections) with generated random numbers can determine the history of each neutron in a batch with a nature of probabilistic occurrences, such as whether the neutron scattering or absorption events occur when the incident neutron collides with the medium, whether the fission event is induced if the neutron absorption event happens, how many numbers of neutrons are generated, and what the new directions will be and how far they would travel before their subsequent collision events happened in a fission event. The positions where the fission events happen at the end of a batch will serve as the new neutron sources’ distribution for the next batch. Meanwhile, the number size of the fission neutrons will be normalized as the neutron weight during the calculation to preserve each batch’s initiating fission neutrons number.

The specific contributions of the histories in a batch can be tallied further to reckon the $k_{eff}$ eigenvalue of that batch. MVP3 provides two methods to evaluate $k_{eff}$ to solve the eigenvalue problems of a system [33]. The first method is based on the definition of the neutron multiplication, as follows:

(20)$k_{eff}=\frac{F}{N},$

where $F$ [-] is the contribution from the weight of neutrons created by the fission events, and $N$ [-] is from the weight of source neutrons. The other method is based on the equation of neutron balance, as below:

(21)$k_{eff}=\frac{F}{L+NA},$

where $L$ [-] is the contribution from the weight of neutrons leaking out of the calculation domain, and $NA$ [-] means the weight of the neutrons under the net absorption events. In MVP3, there are three alternative estimators (Collision estimator, Analog estimator, and Track length estimator [33]) to evaluate the contributions related to the incident neutrons colliding with the medium of a batch ($F$ and $NA$). The ways to evaluate $F$ and $NA$ via the Collision estimator are as below:

(22)$F={\displaystyle \sum _n{W}_n\frac{({\overline{\nu }}_n{\Sigma }_{f,n})}{{\Sigma }_{t,n}}},$

(23)$NA={\displaystyle \sum _m{W}_m\frac{{\Sigma }_{NA,m}}{{\Sigma }_{t,m}}},$

where $n$ and $m$ express the indexes for the neutrons inducing the fission and the net absorption events in a batch, respectively. $W_n$ [-] and $W_m$ [-] are the weight of the neutron n and m. ${\Sigma }_{f,n}$ [1/m] and ${\overline{\nu }}_n$ [-] are the macroscopic fission cross section and the average number of generated neutrons in the neutron $n$’s fission event. ${\Sigma }_{t,n}$ [1/m] and ${\Sigma }_{t,m}$ [1/m] are the macroscopic total cross sections for the neutron $n$ and $m$, respectively. ${\Sigma }_{NA,m}$ [1/m] stands for the macroscopic net absorption cross section in the net absorption event of neutron $m$, which can be expressed as below:

(24)${\sum }_{NA,m}={({\sum }_a-{\sum }_{n,2n}-2{\sum }_{n,3n}-3{\sum }_{n,4n})}_m,$

where ${\sum }_a$ [1/m] is the macroscopic absorption cross section. ${\sum }_{n,2n}$ [1/m], ${\sum }_{n,3n}$ [1/m], and ${\sum }_{n,4n}$ [1/m] are the macroscopic cross sections for (n, 2n), (n, 3n), and (n, 4n) reactions, respectively. The forms used to estimate F and NA by the Analog estimator are based on the Collision estimator, just substituting the microscopic types for the macroscopic cross sections [33]. As for the Track length estimator, F and NA are evaluated based on the corresponding neutron length traveled in the medium, as follows:

(25)$F={\displaystyle \sum _n{W}_n{\delta }_n({\overline{\nu }}_n{\Sigma }_{f,n})},$

(26)$NA={\displaystyle \sum _m{W}_m{\delta }_m{\Sigma }_{NA,m}},$

where ${\delta }_n$ [m] and ${\delta }_m$ [m] mean the distance travelling by the neutron since the last event. Therefore, with the combination of the two $k_{eff}$ evaluation methods and three different estimators, six estimated $k_{eff}$ eigenvalues can be provided for each batch.

With the batch marching forward iteratively, the $k_{eff}$ values estimated by a particular combination of an estimator and the $k_{eff}$ evaluation method from each batch are accumulated until the end of the user-stated batch, and averaged; hence, six averaged $k_{eff}$ from the different estimated methods can be obtained at the end of all simulated batches separately. Nevertheless, the $k_{eff}$ value estimated from a batch can be meaningful only after the spatial distribution of the fission neutron source converges. Therefore, the estimated $k_{eff}$ values from the various initial batches (called the inactive batches) would be discarded until the fission neutron source has been converged, and the $k_{eff}$ values after inactive batches (called the active batches) have been collected. Eventually, the most confident $k_{eff}$ is obtained using the maximum likelihood method based on the final six averaged $k_{eff}$ values from different estimated methods [33]:

(27)$k_{eff}=\frac{{\sum }_a{\sum }_b{T}_{a,b}{X}_a}{{\sum }_a{\sum }_b{T}_{a,b}},$

where $X_a$ stands for each final average $k_{eff}$ value from a specific combination of an estimator and the $k_{eff}$ evaluation method, and $T_{a,b}$ is the element of the inverse matrix of the covariance estimated by these six final average $k_{eff}$ values.

2.3. Coupling Methodology and Process of the CFD–DEM–MCN Solver

The CFD–DEM–MCN coupled solver built in this research was executed automatically in a cyclical manner, composed of two separately executing codes (STAR-CCM+ code for CFD–DEM solver, MVP3 code for MCN solver) and coupled externally using the originally created interface (MVP3 and STAR-CCM+ interface). The MVP3 and STAR-CCM+ interface, based on the Python language with the assistance of Java script, is in charge of data exchange and workflow control between two codes. The flowchart of the CFD–DEM MCN coupled solver is also presented in Figure 3.

After initializing the MVP3 and STAR-CCM+ models, the MVP3 input file is called by the interface and transferred to the MVP3 code to solve the $k_{eff}$ eigenvalue problem. The solved $k_{eff}$ value saved in the MVP3 output file will be extracted and transformed to the CSV file type to be collected by the interface. Subsequently, if the debris bed is under a sub-critical state ($k_{eff}$ < 1), the calculation will keep processing forwards; otherwise, the power excursion due to re-criticality can be treated as an occurrence, and the calculation of the CFD–DEM–MCN coupled solver will be terminated here. For the case of the sub-critical state, the Java macro file, which is the Javascript-based file executable for the STAR-CCM+ workflow being controlled externally, will be automatically created and called in the MVP3 and STAR-CCM+ interface to activate the STAR-CCM+ calculation of the current CFD timestep ($t_{CFD}^n$) with the parameters updating, such as particle heat source. In this research, based on the conservative concerns in the cooling characteristic analysis, the particle heat source transferred from the interface to STAR-CCM+ is assigned as the constant value without utilizing the decreasing decay heat with time passing.

After the STAR-CCM+ calculation has been activated, the DEM solver in the STAR-CCM+ will start to calculate the force and the heat transfer from the fluid phase to particles by referring to the fluid phase information in the CFD solver (referring to the initial condition if at the first timestep), subsequently estimating the force and heat transfer between particles. This process in one elapse of the DEM timestep ($\Delta t_{DEM}$) will be repeated continuously until the accumulated $\Delta t_{DEM}$ is equal to one elapse of the CFD timestep ($\Delta t_{CFD}$), eventually transferring the solved results of particle information and the momentum/energy feedback to the CFD solver. After updating the fluid volume fraction and momentum/energy source terms in the CFD cells based on the transferred data from the DEM solver, the fluid information will be evaluated until the result has converged in the CFD solver, transferring the converged fluid phase information to the DEM solver as the reference for the next CFD timestep ($t_{CFD}^{n+1}$). In addition, the solved particle information (e.g., particle index, particle centroid position, velocity, temperature) and the fluid phase information will also be extracted in the CSV file type and transferred to the MVP3 and STAR-CCM+ interface, temporarily stopping the STAR-CCM+ calculation here.

Before stepping into the next timestep, if the end time ($t_{end}^{}$) has been reached, the CFD–DEM–MCN coupled solver will be stopped here; otherwise, the calculation will proceed to the next timestep with the updating of the heterogeneous bed configuration. Based on each particle’s index and the corresponding particle centroid position obtained from the STAR-CCM+, the MVP3 and STAR-CCM+ interface will automatically generate the new MVP3 input file to reflect the updated heterogenous bed configuration in the MCN solver. Therefore, the CFD–DEM–MCN coupled solver built in this research will be available to investigate the effect of the alternative heterogeneous bed configuration’s transformation on the criticality potential (by collected history data of the $k_{eff}$ values) and cooling characteristics (by collected transient information of particle/fluid phase) simultaneously, considering the momentum/energy exchange between the dynamic heterogeneous particle and fluid phase (two-way coupling in STAR-CCM+) linking to the neutronic calculation (one-way coupling from STAR-CCM+ to MVP3, for the conservation concerned in this research).

3. Dynamic Heterogeneity Investigation into Debris Bed Safety Analysis

3.1. Explanation of the Investigated Debris Bed Models

Because the primary purpose of this research is to examine the capability of the CFD–DEM–MCN coupled solver among different heterogeneous bed configuration’s transformation processes during particle redistribution, the hypothetically intermediate states of the bed formation in prototype SFRs were prepared as the initial conditions of the investigated bed models for conducting bed safety assessments. The concept of the model settings is explained in this section.

Firstly, all the initial states of the investigated debris bed models in this research are located at the center of the in-vessel debris catcher in SFR’s lower plenum, assuming that all the model geometries are the typical conical shape [24,25] with the same initial bed size and slope angle. This research was not concerned with the simulating process of FCI and the falling of the fragmented particles; therefore, in order to initialize the transformation process of the heterogeneous bed configuration, the initial slope angle in this study was placed higher than the reposed angle to induce particle redistribution (based on the gravity-driven particle avalanching).

Subsequently, the sources used to provide different heterogeneous bed configuration transformations among the investigated debris bed models were from (1) the non-uniformity of the debris particles and (2) the initial distribution of these non-uniform particles. As in the primitive examination of this coupled solver on the bed safety assessment, the particle size, shape, and surface roughness were kept the same, and only the non-uniformity factor of the material between particles was selected in this research, realized by introducing binary mixture particles, including the Fuel particle (composed of MOX fuel) and the Structure particle (composed of stainless steel) shown in Figure 4. Because the composite materials differ between Fuel and Structure particles, other accompanying non-uniformity factors, such as decay heat generating rate, particle density, and other relative physical properties, were also considered in this bed safety assessment.

This research’s first investigated heterogeneous model was prepared by randomly positioning Fuel particles and Structure particles within the debris bed as the initial heterogeneous configuration, named the “Mixed model”. On the other hand, considering the potential scenario in which the molten fuel, when not accompanying the molten structure, creates different discharge orders from the core area to the lower plenum [2], the other investigated heterogeneous debris bed model, named the “Stratified model”, was also prepared in this research by settling Fuel particles on the bottom layer of the debris bed to create an initial stratified configuration, as shown in Figure 4. In addition to the reference material used to compare the heterogeneous models, the “Homogeneous model”, composed of homogeneous particles, was also prepared in this study. It should be noted that all the investigated debris bed models in this research were assumed from the same discharged molten core materials with the same volume; hence, the homogeneous particles’ material composition and other related physical properties were derived from the volume-weighted average values of the total Fuel particles and Structure particles in either the Mixed or Stratified models. For example, the particle density of the Homogeneous particle (${\rho }_{i,H}$) can be derived as below:

(28)${\rho }_{i,H}=\frac{{\rho }_{i,F}{V}_{T,F}}{V_{T,F}+V_{T,S}}+\frac{{\rho }_{i,S}{V}_{T,S}}{V_{T,F}+V_{T,S}},$

where ${\rho }_{i,F}$ and ${\rho }_{i,S}$ are the particle density of the Fuel particle and the Structure particle. $V_{T,F}$ and $V_{T,S}$ are the total particle volumes of Fuel particles and Structure particles in either the Mixed model or the Stratified model.

To sum up, the initial geometry of the conical particulate debris bed and also the summation of volume and the summation of the material composition in a bed were the same among the Homogeneous, Mixed, and Stratified models. In addition, the particle shape, size, and surface roughness were also appointed as fixed values among all the particles (including Homogeneous particles, Fuel particles, and Structure particles). The heterogeneity sources among the three debris models were from the non-uniform fissile material distributions between particles or within the bed, whose configuration transformations are induced via assigning an initial slope angle higher than the repose one. Through the setting described above, the effects of different heterogeneous configuration transformation processes on the debris bed safety assessments can be investigated, and further classified as the (1) particle-centralized fissile material effect, by comparing the results from the Mixed model with the Homogeneous model, and (2) the bed-centralized fissile material effect by comparing the results between the Mixed model and the Stratified model, respectively.

3.2. Simulation Settings

This section discusses the input parameters of three debris bed models (Homogeneous, Mixed, and Stratified models) used in the bed safety assessment via the CFD–DEM–MCN solver in this research. These models rely on the assumptions discussed in the previous section and the current understanding of stochastically occurring processes during severe accidents.

Following the explanation in Section 3.1, the calculating domain in the CFD–DEM solver was set as a simplified three-dimensional cylindrical container referring to the prototype SFR’s lower plenum [38]. The calculating domain contains the working fluid of the liquid sodium (continuous liquid phase in the CFD solver) for removing the decay heat generated from the particulate debris bed (discrete solid phase in DEM solver). The input parameters used to describe the liquid sodium’s physical properties in the CFD solver are summarized in Table 1, and the polynomial function of the temperature is applied to the equation of the state in liquid sodium [39]. The geometry parameters are presented in Figure 5, and the origin of the coordinates was set at the bottom center of the calculating domain. On the top surface of the domain, four cold legs and an outlet were pertained. The inlet flowrates of these cold legs were settled to zero in the CFD solver; therefore, the mechanism of decay heat removal from the debris bed to the outlet (set as a constant pressure boundary of 1 atm) was mainly based on the natural circulation. In addition, the adiabatic condition was also assigned to the other surfaces of the calculating domain.

As for the preparation of the debris bed’s initial configurations among the three models, the particle injector in the DEM solver (provided by STAR-CCM+) was utilized to fill discrete solid particles in an assigned three-dimensional conical space within the calculating domain; thus, the parameters of the debris bed, such as the total mass of the debris bed, discrete particles’ shapes and size, and the geometry size of the assigned conical space, were the required information, and are summarized in Table 2. Firstly, under the postulated severe accident in the prototype SFR, the total mass of the debris bed in this research was estimated based on the discharged molten core materials composed of the binary mixture materials of MOX fuel (as the molten fuel) and stainless steel (to simulate the molten structure). From the previous numerical analysis of the hypothetical severe accident (Unprotected Loss Of Flow, ULOF) in SFR, the results show that around 15% of the total molten core material will be discharged to the lower plenum [3], and one third of the total molten core material leaked to the lower plenum in the prototype SFR is assumed here based on the conservative concern. Based on this conservative assumption, considering the referred parameters of the total MOX fuel inventory [40] and the volume ratio of fuel-to-structure (0.64: 0.36 [41]) in the prototype SFR, the MOX fuel and stainless steel mass of the debris bed in this research were estimated as 27,045.45 kg and 10,856.60 kg, and the total mass of the debris bed could also be derived from the summation of the previous two values. Subsequently, the total materials of the debris bed were discrete because of the simplified spherical particles used to fill in the assigned conical space in the calculating domain, with a postulated radius of 0.065 m referred from the DEFOR-A experiments [17], which consider the potential phenomenon of agglomeration during the fragmentation of the molten core materials. Finally, the geometry size of the assigned conical space can be determined by the total volume of the MOX fuel and stainless steel in the debris bed, the average porosity of the debris bed, and the initial slope angle. The total volume of the MOX fuel and the stainless steel can be estimated by their densities (11,000 kg/m³ for MOX fuel [26], 7850 kg/m³ for stainless steel [42]) with their mass discussed above. An average porosity value of 0.59 is assumed here, referred from the FARO/THERMOS experimental results of spatial porosity distribution via injecting molten UO₂ into the liquid sodium [20]. To the best of our knowledge, there are currently no data available on the slope angle needed to cause particle avalanching in the debris bed of the prototype SFR. Instead, research for Light Water Reactor (LWR) has shown that the slope angle to induce particle avalanching for the debris bed has a lower boundary of 22° [22,27,28], and an initial slope angle of 30° was selected in this study. Therefore, the height (2.485 m) and the bottom radius (1.434 m) of the assigned conical space could be evaluated and utilized to set the initial bed configurations for this research.

Except for the initial geometry parameters, the input parameters for discerning different materials between particles (e.g., Fuel particles made of MOX fuel; Structure particles made of stainless steel) in the calculation are also necessary information for creating the different initial heterogeneous configurations of the debris bed models, as shown in Figure 4. The input parameters of the physical properties of Fuel and Structure particles in the DEM solver are summarized in Table 3. The Young’s modulus and Poisson’s ratio, relating to the particle–particle momentum exchange due to collision, are referenced to the physical properties of MOX fuel [43] and stainless steel [42], respectively. The parameters of either static friction coefficient or rolling resistance coefficient (also relating to the particle–particle collision calculation) are settled as the same number between Fuel and Structure particles based on the same surface roughness assumption between particles (discussed in Section 3.1), and the values are the default numbers provided in STAR-CCM+. As for the heat capacity and thermal conductivity relating to the particle–particle energy exchange, the settled parameters in the Fuel particle and the Structure particle are also referred from the material properties of the MOX fuel [44] and stainless steel [42]. The decay heat generation rate settled for the Fuel particle in this research was referred from the hypothetical scenario of ULOF in an SFR, which implemented a fuel assembly with an inner duct structure (Fuel Assembly with Inner Duct Structure: FAIDUS). The previous studies suggest that the time from the start of core melting to the onset of FCI is about a few minutes, after which the debris bed is estimated to be piled up on the debris catcher of the lower plenum in around one minute, considering the debris falling rate [1,45]. Therefore, the decay heat generation rate of the debris bed in such a short time can be considered equivalent to the decay heat soon after the normal reactor shutdown (6–7% of operating power) [46]. Based on conservative concern, the total decay heat generation rate of the debris bed in this research was estimated according to the assumption of one third of the total core molten fuel of the prototype SFR, maintaining at 7% of one third of its thermal operating power [40]. The value of the heat generation rate for each Fuel particle is derived by dividing the total volume of the MOX fuel in the debris bed (without setting the heat source in Structure particles). Finally, the Homogeneous particle physical properties are based on the conception that the MOX fuel and stainless steel in the debris bed are evenly distributed to each particle. Hence, the input parameters of each Homogeneous particle can be derived using the volume-weighted average value of the total Fuel particles and Structures from either the Mixed model or the Stratified model (via the same approach as Equation (28), applying the volume ratio of MOX fuel to stainless steel in Table 2). The results of the derived Homogeneous particle’s input parameters are also listed in Table 3. Here, based on the information in Table 2 and Table 3, the initial configuration of the Homogeneous model can be created by randomly filling Homogeneous particles into the assigned conical space in the calculating domain, and the Mixed model’s initial configuration can be created by randomly filling all the Fuel and Structure particles into the exact geometry of the assigned conical space. For preparing the initial configuration of the Stratified model, the Fuel particles in the Mixed model are filled into the bottom layer of the same assigned conical space, and the Structure particles left in the Mixed model are also filled into the upper layer of the conical space.

A mesh-independent test for the CFD–DEM solver was also conducted via different structures of the computational meshes in the fluid region based on the Homogeneous model, as shown in Figure 6.

The polyhedral meshes were utilized to discretize the fluid region. Based on the time discretization scheme of the implicit unsteadiness (timestep of 0.01 s in the CFD solver, five sub-steps in the DEM solver), at the time after 100 s of initializing the CFD–DEM solver in the Homogenous model, the temperature distribution of the liquid sodium due to the decay heat generated from the Homogeneous particles became stable. The different average temperatures of the liquid sodium in the calculating domain ($T_{f,Avg.}$ [K]) corresponding to different CFD mesh structures are also listed in Table 4. In this study, it was noticed that refining the CFD meshes excessively to a certain number (53,149) did not have a significant impact on the deviation of $T_{f,Avg.}$ from the coarser mesh structure. As a result, the Mixed and Stratified models will be further investigated, employing the same mesh structure of the 53,149 mesh number.

In the MCN solver, the geometry setting of the calculating domain was the same as the one applied in the CFD–DEM solver, and the particle positions in the calculating domain were updated following the results calculated from the CFD–DEM solver. In addition, as well as the conception in the CFD–DEM solver, the Homogeneous particles, Fuel particles, and Structure particles were also discerned by assigning corresponding material compositions in the MCN solver, and the input parameters of the material compositions for different particles are summarized in Table 5. The material composition of the Structure particle was referred from the stainless steel [47]. In terms of the Fuel particle, the fresh MOX fuel composition, without accounting for the depletion process in SFRs, was utilized in this study for conservative concerns in the criticality safety assessment. Considering the transition stage from LWRs to SFRs, both types of reactors would be in service concurrently. This drives the possible scenario of supplying SFRs MOX fuel with transuranium (TRU) acquired from recycled LWR spent fuel [40]. Therefore, the fresh MOX fuel composition in this study was referred from the TRU composition of recycled spent fuel from Advanced LWR (ALWR) based on the scenario that the waiting time for decay heat cooling was 40 years from the state of an average discharge burnup of 60 GWd/t [40]. The contribution of the discharge boron materials from the backup CRGTs in prototype SFR is also considered in this research [2]. As the same approach was used to derive input parameters of physical properties for Homogeneous particles in the CFD–DEM solver, the composition materials of each Homogeneous particle in the MCN solver were generated using the volume-weighted average value of the total Fuel particles and Structures in either the Mixed model or the Stratified model. This study employed the continuous cross-section data library of JENDL-4.0 [48], setting the number of neutron histories and active batches as 10,000 and 220, respectively. In addition, the number of inactive batches excluded from statistical processing was set as 20. Based on these calculation conditions, the uncertainties of all the computed values of $k_{eff}$ in this research were below 0.035%.

3.3. Results and Discussion

3.3.1. Transformation of Heterogeneous Configuration

Figure 7 displays the transformation processes of the bed configuration during the particle redistribution on these three debris bed models, including the Homogeneous, Mixed, and Stratified models. The particle in the grey in Figure 7 stands for the Homogeneous particle as the element of the Homogenous model. On the other hand, the black and white particles in the Mixed and Stratified model represent the Fuel particle and Structure particle, respectively.

In the beginning, since the slope angle of the conical debris bed model was higher than the stable repose angle, particle avalanching was induced, letting particles slide downwards and being redistributed in all three models. During the particle redistribution, all three debris bed models exhibited a similar trend of being flattened gradually, suppressing their slope angles and the heights of the debris beds. In the end, particle redistribution in all three models ceased after a new stable bed configuration with lower slope angles and bed heights within 10 s had been achieved, as shown in Figure 7.

Meanwhile, the heterogeneous debris bed models (Mixed model and Stratified model) revealed different characteristics of the bed configuration’s transformation during particle redistribution from the Homogeneous model. As for the Mixed model, both light Structure particles and weighty Fuel particles moved downwards during particle redistribution, and the light Structure particles tended to gradually disperse outwards in the debris bed after a series of particle–particle collisions with density differences, as shown in Figure 7b. In addition, this indicated that the weighty Fuel particles tended to stay and be concentrated within the debris bed. The characteristics of the Mixed model’s configuration transformation also reflect the previous experimental results, which demonstrated that heavier particles take over the bottom layer’s larger ratio after the sedimentation of binary-density particles [14,15]. Finally, because of the effect of the dispersed light Structure particles, the height of the debris ($H_{bed}$ [m]) was lower than the other two models after particle redistribution, as shown in Figure 8. For the Stratified model, although both Fuel particles and Structure particles also moved downwards, the initial stratified configuration was not altered during the particle redistribution since the Structure particles with lower density on the upper layer found it tough to intrude into the bottom high-density layer of the Fuel particles, as in Figure 7c. In the end, the $H_{bed}$ of the Stratified model was higher than the other two models after particle redistribution, as demonstrated in Figure 8. The values of the $H_{bed}$ after particle redistribution among the three models are also listed in Table 6.

3.3.2. Criticality Analysis

This section examines how changes in the bed’s heterogeneous configurations impact the assessment of its criticality, including analyzing the $k_{eff}$ alteration during particle redistribution and the safety margin for criticality safety after particles have ceased moving, using re-criticality accidents as an indicator. The effects of various transformations in heterogeneous configurations can be categorized into (1) the Particle-Centralized fissile material effect (PC effect), which compares the Mixed model to the Homogeneous model, and (2) the effect of Bed-Centralized fissile material at the bottom of the debris bed (BC effect), which compares the Stratified model to the Mixed model.

Regarding the PC effect, the Mixed and Homogeneous models display a decrease in $k_{eff}$ eigenvalues during particle redistribution, as demonstrated in Figure 9. This decline continues until the redistribution of particles stops, indicating that the decreased height of the debris bed during particle redistribution effectively suppresses the non-leakage probability ($P_{NL}$ [-]). To delve further into the PC effect on the $k_{eff}$ change, Figure 10 illustrates the relationship between the change in $k_{eff}$ ($\Delta k_{eff}(t)$ [-]) and the change in debris bed height ($\Delta H_{bed}(t)$ [m]) for the Mixed and Homogeneous models.

The definitions of $\Delta k_{eff}(t)$ and $\Delta H_{bed}(t)$ are as below:

(29)$\Delta k_{eff}(t)=k_{eff}(t)-k_{eff}(t=0),$

(30)$\Delta H_b{}_{ed}(t)=H_{bed}(t)-H_{bed}(t=0),$

where $t$ [s] stands for the time after particle redistribution starting. Figure 10a shows that the sensitivity of $\Delta k_{eff}(t)$ to $\Delta H_{bed}(t)$ in the Mixed model is higher than in the Homogeneous model, showing that the accumulated fissile materials’ movements would strengthen the effect of the declined $P_{NL}$ when the debris bed height is suppressed during particle redistribution. In addition, the Mixed model also exhibited a greater decreased range in the debris bed height than the Homogeneous model during the particle redistribution (Figure 10b); therefore, the PC effect in the Mixed model caused a higher decreased range of $k_{eff}$ than the Homogeneous model after particle redistribution, as in Figure 10c. Nevertheless, since the PC effect in the Mixed model also contributes to the higher thermal utilization factor  [-] (due to the higher concentrated fissile material in Fuel particles) and the greater resonance escape probability $P_p$ [-] (due to the self-shielding effect from the Structure–Fuel–Structure spatial arrangement [49]), the value of $k_{eff}$ in the Mixed model was still higher than the Homogeneous model after particle redistribution, in which the results are listed in Table 7.

As for the BC effect, the Stratified model also presented a declining trend during particle redistribution until particles had ceased, as shown in Figure 9. Following this, in the initial stage of the particle redistribution, the Stratified model displayed a higher sensitivity of $\Delta k_{eff}(t)$ to $\Delta H_{bed}(t)$ than the Mixed model (Figure 10a), showing that the downward movement of the accumulated fissile material layer in the Stratified model could further strengthen the phenomenon of the $P_{NL}$’s declination when debris bed height was suppressed. However, in the latter stage of the particle redistribution, the BC effect in the Stratified model for $P_{NL}$’s declination was weakened when the bottom fissile material layer became more flattened, and the Stratified model demonstrated a lower sensitivity of $\Delta k_{eff}(t)$ to $\Delta H_{bed}(t)$ than the Mixed model, in the end. Meanwhile, in Figure 10b, the Stratified model also showed a minor decreased range of the debris bed height compared to the Mixed model; hence, the BC effect in the Stratified model brought a lower decreased range of $k_{eff}$ than the Mixed model, as shown in Figure 10c. Nonetheless, since the BC effect in the Stratified model also provided a negative effect on $P_p$, due to the loss of the spatial heterogeneity on the bottom fissile material layer, the Stratified model was grounded on the lower $k_{eff}$ value compared to the Mixed model after particle redistribution in the end, and the results are also listed in Table 7.

Based on the results above, the PC and BC effects on the $k_{eff}$’s change ($\Delta k_{eff}$) and safety margin of criticality safety ($1-k_{eff}$) after particle redistribution can be evaluated in the form used below:

(31)$\mathrm{PC} \mathrm{effect}=\frac{X_{\mathrm{Mixed} \mathrm{model}}-X_{\mathrm{Homogeneous} \mathrm{model}}}{X_{\mathrm{Homogeneous} \mathrm{model}}},$

(32)$\mathrm{BC} \mathrm{effect}=\frac{X_{\mathrm{Stratified} \mathrm{model}}-X_{\mathrm{Mixed} \mathrm{model}}}{X_{\mathrm{Mixed} \mathrm{model}}},$

where $X$ means the target result (such as $\Delta k_{eff}$ or $1-k_{eff}$), and the subscripts of $X$ stand for the corresponding debris bed models. The results of PC and BC effects in the criticality assessment, applying Equations (31) and (32), are listed in Table 8.

In summary, PC and BC effects can contribute to the positive and negative effects on $\Delta k_{eff}$ (in Table 8) due to the strengthening (in PC effect) and weakening (in BC effect) of $P_{NL}$’s declination during particle redistribution, respectively. For the criticality safety margin of $1-k_{eff}$, however, PC and BC effects were the negative and positive dominant effects, respectively (in Table 8), since the $P_p$ was strengthened by the heterogeneous spatial arrangement of the Fuel–Structure particles in the PC effect, and vice versa in BC effect. Meanwhile, all the models were kept at sub-critical states during the particle redistribution process (as shown in Figure 9). Thus far, it can be proven that the CFD–DEM–MCN coupled solver is capable of assessing the impact of various transformations of the heterogeneous configuration of the bed on the criticality safety assessment of hypothetical SAs in SFRs.

3.3.3. Cooling Characteristics Analysis

The purpose of this section is to investigate the effects of the alternative heterogeneous configurations’ transformations (divided into PC and BC effects in this research) on the debris bed’s cooling safety assessment, including the heat removal efficiency of the debris bed (taking the average ${\mathrm{Nu}}_p$ of the debris bed as the criteria in this research) and the safety margin of stable particle cooling (measured by the occurrence of particle re-melting accidents).

Since the debris bed’s cooling was still under development after 10 s, the following transient behaviors of the bed’s decay heat removal are shown in Figure 11 and Figure 12.

Figure 11 demonstrates the flow pattern of the liquid sodium from 10 s to 100 s in three debris bed models. The decay heat generated from the debris particles induced the buoyance-driven convective flow upwards to remove the debris bed’s heat. In the meantime, the liquid sodium adjacent to the debris bed was drawn into the porous area due to the pressure drop from the left convective flow, facing the flow resistance offered by the drag force as the obstacle to the flow circulation inside the debris bed. Subsequently, the upward convective flow brought the decay heat to the outlet of the calculating domain (center of the top surface), and it sank to the lower part afterward. Ultimately, the natural circulation in the clockwise flow pattern was achieved at 100 s for all debris bed models. Figure 12 presents the temperature distribution of the liquid sodium from 10 s to 100 s among three debris bed models. Corresponding to the flow pattern, the liquid sodium’s temperature gradient was gradually developed from the debris bed to the outlet, and, in the end, the maximum temperatures of the liquid sodium ($T_{f,Max}$ [K]) in Figure 13 also approached stability without liquid sodium boiling in all the three models (Sodium’s boiling temperature is 1153 K at 1 atm [39]). Therefore, the subsequent investigation into the PC and BC effects on the cooling assessment was based on the information collected at 100 s.

As for the PC effect, firstly, it could affect the flow around the heat-generating particles (Fuel particles in the Mixed model; Homogeneous particles in the Homogeneous model), and the resulting flow velocity distributions collected from where the heat-generating particles were located are shown in Figure 14, including the radial direction and the vertical direction of the debris bed. $\left|{\overrightarrow{V}}_{f(p)}\right|$ [m/s] stands for the flow velocity where the heat-generating particles’ centroids are located, and $R_p$ [m] and $H_p$ [m] mean the heat-generating particles’ radial and vertical locations, defined as below:

(33)$R_p=\sqrt{x^2+y^2},$

(34)$H_p=z,$

where $x$ [m], $y$ [m], and $z$ [m] stand for the particle’s coordinates in the calculating domain. Figure 14a shows that $\left|{\overrightarrow{V}}_{f(p)}\right|$ presented the declined trend of the inward flow velocity in both Mixed and Homogeneous models, since $\left|{\overrightarrow{V}}_{f(p)}\right|$ along the radial direction was suppressed by the flow resistance (drag force) in the porous area. The average $\left|{\overrightarrow{V}}_{f(p)}\right|$ (${\left|{\overrightarrow{V}}_{f(p)}\right|}_{Avg.}$ [m/s]) in the radial direction also showed similar values between the two models (listed in Table 9). By contrast, $\left|{\overrightarrow{V}}_{f(p)}\right|$ presented the increased trend of the upward vertical flow velocity in both Mixed and Homogeneous models (as shown in Figure 14b), which reflected that $\left|{\overrightarrow{V}}_{f(p)}\right|$ along the vertical direction was mainly affected by the dominant buoyance force relating to the heat-generating particles. Furthermore, compared with the Homogenous model, although the PC effect in the Mixed model caused the higher increased $\left|{\overrightarrow{V}}_{f(p)}\right|$ gradient along the vertical direction (in Figure 14b) because of the higher value of the accumulated decay heat source along the vertical direction in the Mixed model, the ${\left|{\overrightarrow{V}}_{f(p)}\right|}_{Avg.}$ on the vertical direction in the Mixed model was still lower than the Homogeneous model (listed in Table 9) since the heat-generating particles in the Mixed model were located at the relatively lower place of the debris bed after particle redistribution. Meanwhile, the average ${\mathrm{Nu}}_p$ of the heat-generating particles in a debris bed (${\mathrm{Nu}}_p{}_{,Avg.}$ [-]), which can represent the decay heat removal efficiency of the debris bed, also reflected the influence of the affected $\left|{\overrightarrow{V}}_{f(p)}\right|$. ${\mathrm{Nu}}_p$’s accumulated frequency curve of Mixed/Homogeneous models and the corresponding ${\mathrm{Nu}}_p{}_{,Avg.}$ are presented in Figure 15 and Table 9, respectively. The lower ${\mathrm{Nu}}_p{}_{,Avg.}$ in the Mixed model compared to the Homogeneous model (as shown in Table 9) presented the PC effect, suppressing the ${\left|{\overrightarrow{V}}_{f(p)}\right|}_{Avg.}$ along the vertical direction in the Mixed model, and playing the dominant factor in deteriorating the performance of the Mixed model‘s ${\mathrm{Nu}}_p{}_{,Avg.}$. Finally, the particle temperature distributions of Mixed and Homogeneous models after stable cooling are also exhibited in Figure 16. Although the relatively random distributions of the particle temperature existed in the Mixed model because of the uneven distribution of the heat-generating Fuel particles within the debris bed, both Mixed and Homogenous models still showed a similar trend in the spatial particle temperature distribution, which was related to the built flow pattern. Along the radial direction, the particle temperatures increased towards the center of the debris bed due to the declined inward flow velocity in both models. On the other hand, both models displayed the decreased gradient of the particle temperature along the upward vertical direction due to the enhanced convective flow velocity. Hence, for both models, the particle with maximum temperature after stable cooling was located at the bottom center area of the debris bed, and the values of the maximum particle temperature ($T_{p,Max}$ [K]) for the corresponding models are listed in Table 9. The $T_{p,Max}$ in the Mixed model (from the Fuel particle) was ultimately higher than in the Homogeneous model, which reflects the overall impact of the PC effect in the Mixed model, including the weakened ${\mathrm{Nu}}_p{}_{,Avg.}$ and the heightened heat source of each heat-generating particle due to the particle-centralized fissile material (Fuel particle) present in the Mixed model.

The BC effect in the Stratified model is subsequently discussed by comparing it with the Mixed model. As shown in the Mixed model, the Stratified model also exhibits a decreased trend of inward flow velocity along the radial direction of the debris bed (Figure 14a), and an increased trend of upward flow velocity along the vertical direction (Figure 14b). Moreover, the BC effect also causes a lower value of ${\left|{\overrightarrow{V}}_{f(p)}\right|}_{Avg.}$ in the vertical direction in the Stratified model compared to the Mixed model (as shown in Table 9) due to the comparative lower position distribution of the heat-generating particles than in the Mixed model (as shown in Figure 14b) after particle redistribution. In the meantime, the BC effect on the ${\left|{\overrightarrow{V}}_{f(p)}\right|}_{Avg.}$, along the vertical direction also contributed to the influence to weaken ${\mathrm{Nu}}_p{}_{,Avg.}$, resulting in the value of ${\mathrm{Nu}}_p{}_{,Avg.}$ in the Stratified model being lower than in the Mixed model, eventually (shown in Table 9). Regarding the particle temperature distribution in the Stratified model in Figure 16, it also displayed a similar trend with other models, in that the particle temperature is increased with suppressed inward flow velocity along the radial direction, and the particle temperature is decreased with the enhanced upward convective flow along the vertical direction. Finally, the $T_{p,Max}$ of the Stratified model (from the Fuel particle), located at the bottom central area of the debris bed, as well, was higher than in the Mixed model (shown in Table 9). This can reflect the total impact of the BC effect in the Stratified model on $T_{p,Max}$, including the weakened ${\mathrm{Nu}}_p{}_{,Avg.}$ of the debris bed and the deteriorated particle–particle heat conductive transfer for the Fuel particles due to the evident bed-centralized fissile material area on the bottom of the Stratified model without Structure particles inside.

Following the same procedure in the criticality assessment, the PC and BC effects on the ${\mathrm{Nu}}_p{}_{,Avg.}$ of the debris bed and the cooling safety margin for avoiding the particle re-melting accident (employing $T_{p,melting}-T_{p,Max}$ as the criteria in this research, where $T_{p,melting}$ is the MOX’s melting point of 2800 K [44]) can be evaluated using Equations (31) and (32), and the results are listed in Table 10.

To sum up, both PC and BC effects provided a dominant negative effect on the ${\mathrm{Nu}}_p{}_{,Avg.}$ (in Table 10), mainly due to the relatively lower flow velocity around the heat-generating particles, relating to the lower located position after the particle redistribution. As for the cooling safety margin of $T_{p,melting}-T_{p,Max}$, both PC and BC effects also exhibited total negative impacts (as shown in Table 10). Except for the effects from the suppressed ${\mathrm{Nu}}_p{}_{,Avg.}$ of the debris bed, the higher heat source of a heat-generating particle due to the particle-centralized fissile material (in the PC effect), and the deteriorated particle–particle heat conductive transfer for a heat-generating particle due to being located at the bed-centralized fissile material area without Structure particles inside (in BC effect), would also contribute to the influence on their own final negative impact on $T_{p,melting}-T_{p,Max}$, eventually delivering a more conspicuous negative impact from the PC effect than the BC effect (Table 10). Here, it can be proved that the coupled CFD–DEM–MCN solver has the capability to evaluate the effects of the different transformations of the heterogeneous bed configurations on the cooling safety assessment for the postulated SAs in SFRs.

4. Conclusions

This research aims to develop a numerical tool for investigating the effects of different dynamic transformations of the heterogeneous bed configurations on the debris bed’s criticality/cooling safety assessment under hypothetical SAs in SFRs. To realize this goal, the coupled CFD–DEM–MCN methodology was newly proposed in this research, which can relieve the bed’s safety assessment from the constraint of the static particle conditions, considering the multiphysics interplay between the thermal-hydraulics and neutronics. In addition, the coupled CFD–DEM–MCN solver used in an automatically cyclical manner was developed successfully using the originally created interface (workflow control and data exchange), integrating two other existing codes (the STAR-CCM+ code for the CFD–DEM solver, and the MVP3 code for the MCN solver). The effects of different transformation processes of the heterogeneous bed configuration induced by the particle redistribution, including (1) the effect of particle-centralized fissile material (PC effect, by comparing the Mixed model with the Homogeneous model) and (2) the effect of the bed-centralized fissile material at the bottom of the debris bed (BC effect, by comparing the Stratified model with Mixed model), were quantitatively confirmed in this research. The analytical results showed that PC and BC effects had negative and positive dominant effects on securing the criticality safety margin ($1-k_{eff}$) for the criticality assessment after particle redistribution had ceased. For the cooling assessment, when attaining the stable cooling state, the results indicated that both PC and BC effects had dominant negative effects ensuring the safety margin of avoiding particle re-melting accidents ($T_{p,melting}-T_{p,Max}$), showing a more evident impact from the PC effect than the BC effect. Therefore, this can prove that the newly developed tool, the coupled CFD–DEM–MCN solver, is capable of investigating the effects of different dynamic heterogeneous bed configurations on the criticality and cooling safety assessment of debris beds, considering multiphysics interaction.

As the preliminary examination of the newly developed numerical tool, the current study focused on the capability of clarifying the uncertainty from the factor of the dynamic heterogeneous configuration on the debris bed’s safety assessment by assuming the hypothetically intermediate configurations of the debris bed in the process of the debris sedimentation as the initial state in this research. Hence, the whole debris bed’s sedimentation process, including the fragmented particles’ falling in the lower plenum of SFRs after FCI, will be integrated into future bed safety assessments by conducting further model validations with more referable materials for the SAs analysis. In addition, along with considering the impacts of different molten fuel compositions and decay heat transients corresponding to various conceivable core burnups in future studies, the effects from other additional heterogeneities of the debris bed, such as the particle size, shape, and the ratio of the fuel to structure material, will also be further investigated with the applicable range, in order to cover potential scenarios existing in the stochastic SAs’ development processes. Therefore, the newly developed coupled CFD–DEM–MCN solver can still be expected to be an essential cornerstone in the area of SA safety analysis.

Author Contributions

Conceptualization, K.O., T.T. and N.E.; methodology, C.-Y.L., N.E. and M.P.; software, C.-Y.L. and M.P.; formal analysis, C.-Y.L.; investigation, C.-Y.L. and K.W.; resources, K.O. and Y.O.; data curation, C.-Y.L.; writing—original draft preparation, C.-Y.L.; writing—review and editing, A.U. and K.W.; visualization, C.-Y.L. and A.U.; supervision, K.O. and Y.O. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Acknowledgments

The authors acknowledge the support of Kenji YOKOYAMA at the Japan Atomic Energy Agency, Japan, for the assistance of simulation analysis via MVP3.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The schematic of the mechanisms regarding the debris bed safety assessment under the scenario of heterogeneous particle redistribution.

Figure 2. CFD–DEM model validation: (a) experimental setup, (b) CFD–DEM model, and (c) FFT analysis on the stainless steel particle’s movement.

Figure 3. The flowchart of the CFD–DEM–MCN coupled solver.

Figure 4. Schematics of the three debris bed models’ initial configurations in the CFD–DEM–MCN coupled solver (vertical cross-section of the debris bed).

Figure 5. Schematics and geometry parameters of the calculating domain in the CFD–DEM–MCN coupled solver.

Figure 6. CFD mesh configurations in the mesh-independent examinations.

Figure 7. Transient process of the particle redistribution among three debris bed models: (a) Homogeneous model, (b) Mixed model, and (c) Stratified model.

Figure 8. Historical data of the debris bed’s height during particle redistribution among three models.

Figure 9. History data of the [Forumla omitted. See PDF.] eigenvalues during particle redistribution among three debris bed models.

Figure 10. Neutronic analysis during particle redistribution among three models: (a) ratio of the ’s change to the bed height’s change, (b) the bed height’s change, and (c) [Forumla omitted. See PDF.]’s change.

Figure 11. Liquid sodium’s flow pattern for 10–100 s (vertical cross-section of the calculation domain): (a) Homogeneous model, (b) Mixed model, and (c) Stratified model.

Figure 12. Liquid sodium’s temperature distribution for 10–100 s (vertical cross-section of the calculation domain): (a) Homogeneous model, (b) Mixed model, and (c) Stratified model.

Figure 13. Historical data of the liquid sodium’s maximum temperature among three debris bed models.

Figure 14. The flow velocity distribution collected from where heat-generating particles are located after stable cooling (at 100 s) among three models: (a) along the radial direction (0.13 m ≤ Hp ≤ 0.26 m), and (b) along the vertical direction (0 m ≤ Rp ≤ 0.13 m).

Figure 15. The accumulated frequency of the heat-generating particle’s Nup after stable cooling (at 100 s), among three models.

Figure 16. Particle temperature distribution after stable cooling (at 100 s) (half region of the cylindrical calculation domain): (a) Homogeneous model, (b) Mixed model, and (c) Stratified model.

Table 1

Input parameters of the liquid sodium’s physical properties in the CFD solver.

Physical Properties	Value
Dynamic viscosity [kg/m·s]	2.01 × 10−4
Heat capacity [J/kg·K]	1252
Thermal conductivity [W/m·K]	58.34

Table 2

Geometry parameters of the three debris bed models’ initial configurations.

Parameter	Value
Debris bed shape	Conical shape
Debris bed height [m]	1.434
Debris bed bottom radius [m]	2.485
Slope angle [°]	30
Total debris bed mass (MOX fuel + stainless steel) [kg]	37,902.05
MOX fuel mass [kg]	27,045.45
Stainless steel mass [kg]	10,856.60
Average porosity [-]	0.59
Volume ratio of MOX fuel to stainless steel	0.64: 0.36
Debris shape	Spherical shape
Debris particle radius [m]	0.065

Table 3

Input parameters of the physical properties among different particles in the DEM solver.

Physical Properties	Fuel Particle	Structure Particle	Homogeneous Particle
Density, [kg/m3]	1.10 × 104	7.85 × 103	9.84 × 103
Young’s modulus [GPa]	220	200	212.8
Poisson’s ratio [-]	0.33	0.29	0.32
Static friction coef. [-]	0.61	0.61	0.61
Rolling resistance coef. [-]	0.001	0.001	0.001
Heat capacity [J/kg·K]	330	500	319.2
Thermal conductivity [W/m·K]	3	16.2	7.75
Decay heat generation rate [W/m3]	3.388 × 107	0	2.168 × 107

Table 4

Results of mesh-independent examinations on the Homogeneous model.

Number of CFD Meshes [-]	Tf, Avg. [K]	Relative Deviation ofTf, Avg. [%]
34,516	960.34	-
53,149	961.36	0.106
88,088	961.83	0.049

Table 5

Input parameters of the material composition among different particles in the MCN solver.

Material	Fuel Particle *	Structure Particle *	Homogeneous Particle *
U235	4.91 × 10−5	-	3.11 × 10−5
U238	1.93 × 10−2	-	1.23 × 10−2
Pu238	5.65 × 10−5	-	3.58 × 10−5
Pu239	2.77 × 10−3	-	1.76 × 10−3
Pu240	1.64 × 10−3	-	1.04 × 10−3
Pu241	2.18 × 10−4	-	1.38 × 10−4
Pu242	1.97 × 10−4	-	1.25 × 10−4
Np237	2.74 × 10−5	-	1.74 × 10−5
Am241	1.08 × 10−4	-	6.83 × 10−5
Am243	5.34 × 10−5	-	3.39 × 10−5
Cm244	5.32 × 10−5	-	3.37 × 10−5
O16	4.89 × 10−2	-	3.10 × 10−2
B10	6.28 × 10−4	6.28 × 10−4	6.28 × 10−4
B11	2.28 × 10−4	2.28 × 10−4	2.28 × 10−4
Fe	-	7.01 × 10−2	2.57 × 10−2
Cr	-	1.20 × 10−2	4.38 × 10−3
W	-	7.19 × 10−4	2.63 × 10−4
Ti	-	9.87 × 10−4	3.61 × 10−4

* Unit: Atomic number density [10²⁴ atoms/cm³].

Table 6

Height of the debris bed after particle redistribution stopped among three models.

	Homogeneous Model	Mixed Model	Stratified Model
$H_{bed}$ [m]	0.589	0.536	0.640

Table 7

Neutronic analysis results after particle redistribution had been stopped (at 10 s) among three models.

	Homogeneous Model	Mixed Model	Stratified Model
$\Delta k_{eff}/\Delta H_{bed}$ [1/m]	0.100	0.130	0.119
$\Delta H_{bed}$ [m]	−0.845	−0.898	−0.794
$\Delta k_{eff}$ [-]	−0.084484	−0.116748	−0.094490
$k_{eff}$ [-]	0.865216	0.879252	0.841981

Table 8

PC and BC effects on the criticality assessment results after particle redistribution had been stopped (at 10 s).

	PC Effect [%]	BC Effect [%]
$\Delta k_{eff}$	38.19	−19.06
$1-k_{eff}$	−10.41	30.87

Table 9

Data collected from the heat-generating particles after stable cooling (at 100 s) among three models.

	Homogeneous Model	Mixed Model	Stratified Model
${\left|{\overrightarrow{V}}_{f(p)}\right|}_{Avg.}$$\mathrm{along} R_p$ [m/s] *	0.109	0.110	0.099
${\left|{\overrightarrow{V}}_{f(p)}\right|}_{Avg.}$$\mathrm{along} H_p$ [m/s] **	0.087	0.081	0.064
${\mathrm{Nu}}_{p,}{}_{Avg.}$[-]	21.14	21.05	20.65
$T_{p,Max}$ [K]	1118.88	1161.04	1170.37

* Calculated from the radial direction (0.13 m ≤ H_p ≤ 0.26 m), as in Figure 14a. ** Calculated from the vertical direction (0 m ≤ R_p ≤ 0.13 m), as in Figure 14b.

Table 10

The PC and BC effects on the cooling assessment of the heat-generating particles after stable cooling (at 100 s).

	PC Effect [%]	BC Effect [%]
${\mathrm{Nu}}_p{}_{,Avg.}$	−0.43	−1.90
$T_{p,melting}-T_{p,Max.}$	−2.51	−0.57

References

1. Suzuki, T.; Kamiyama, K.; Yamano, H.; Kubo, S.; Tobita, Y.; Nakai, R.; Koyama, K. A scenario of core disruptive accident for Japan sodium-cooled fast reactor to achieve in-vessel retention. J. Nucl. Sci. Technol.; 2014; 51, pp. 493-513. [DOI: https://dx.doi.org/10.1080/00223131.2013.877405]

2. Suzuki, T.; Tobita, Y.; Nakai, R. Evaluation of recriticality behavior in the material-relocation phase for Japan sodium-cooled fast reactor. J. Nucl. Sci. Technol.; 2015; 52, pp. 1448-1459. [DOI: https://dx.doi.org/10.1080/00223131.2015.1005719]

3. Suzuki, T.; Tobita, Y.; Kawada, K.; Tagami, H.; Sogabe, J.; Matsuba, K.; Ito, K.; Ohshima, H. A preliminary evaluation of unprotected loss-of-flow accident for a prototype fast-breeder reactor. Nucl. Eng. Technol.; 2015; 47, pp. 240-252. [DOI: https://dx.doi.org/10.1016/j.net.2015.03.001]

4. Iwasawa, Y.; Abe, Y. Melt jet-breakup and fragmentation phenomena in nuclear reactors: A review of experimental works and solidification effects. Prog. Nucl. Energy; 2018; 108, pp. 188-203. [DOI: https://dx.doi.org/10.1016/j.pnucene.2018.05.009]

5. Matsuba, K.I.; Kamiyama, K.; Toyooka, J.I.; Tobita, Y.; Zuyev, V.A.; Kolodeshnikov, A.A.; Vassiliev, Y.S. Experimental discussion on fragmentation mechanism of molten oxide discharged into a sodium pool. Mech. Eng. J.; 2016; 3, pp. 15-00595. [DOI: https://dx.doi.org/10.1299/mej.15-00595]

6. Tentner, A.M.; Parma, E.; Wei, T.; Wigeland, R. Severe Accident Approach-Final Report Evaluation of Design Measures for Severe Accident Prevention and Consequence Mitigation; ANL-GENIV-128 Argonne National Laboratory: Argonne, IL, USA, 2010; [DOI: https://dx.doi.org/10.2172/973483]

7. Gunji, S.; Tonoike, K.; Izawa, K.; Sono, H. Study of experimental core configuration of the modified STACY for measurement of criticality characteristics of fuel debris. Prog. Nucl. Energy; 2017; 101, pp. 321-328. [DOI: https://dx.doi.org/10.1016/j.pnucene.2017.03.002]

8. Freiría López, M.; Buck, M.; Starflinger, J. A Criticality Evaluation of Fukushima Daiichi Unit 1 Fuel Debris. Proceedings of the 2018 26th International Conference on Nuclear Engineering (ICONE26); London, UK, 22–26 July 2018; ICONE26-81148 [DOI: https://dx.doi.org/10.1115/ICONE26-81148]

9. Li, C.Y.; Takata, T.; Okamoto, K. Numerical Investigation and Analysis on Debris Bed Coolability and Re-criticality for Japan Sodium Cooled Fast Reactor (JSFR). Proceedings of the 2019 27th International Conference on Nuclear Engineering (ICONE27); Ibaraki, Japan, 19–24 May 2019; ICONE27-1319 [DOI: https://dx.doi.org/10.1299/jsmeicone.2019.27.1319]

10. Lipinski, R.J. Model for Boiling and Dryout in Particle Beds; NUREG/CR-2646 Sandia National Laboratories: Albuquerque, NM, USA, 1982.

11. Hu, K.; Theofanous, T.G. On the measurement and mechanism of dryout in volumetrically heated coarse particle beds. Int. J. Multiph. Flow.; 1991; 17, pp. 519-532. [DOI: https://dx.doi.org/10.1016/0301-9322(91)90047-7]

12. Schulenberg, T.; Müller, U. An improved model for two-phase flow through beds of coarse particles. Int. J. Multiph. Flow.; 1987; 13, pp. 87-97. [DOI: https://dx.doi.org/10.1016/0301-9322(87)90009-7]

13. Takasuo, E.; Holmström, S.; Kinnunen, T.; Pankakoski, P.H. The COOLOCE experiments investigating the dryout power in debris beds of heap-like and cylindrical geometries. Nucl. Eng. Des.; 2012; 250, pp. 687-700. [DOI: https://dx.doi.org/10.1016/j.nucengdes.2012.06.015]

14. Sheikh, M.A.R.; Son, E.; Kamiyama, M.; Morioka, T.; Matsumoto, T.; Morita, K.; Matsuba, K.; Kamiyama, K.; Suzuki, T. Sedimentation behavior of mixed solid particles. J. Nucl. Sci. Technol.; 2018; 55, pp. 623-633. [DOI: https://dx.doi.org/10.1080/00223131.2017.1419888]

15. Sheikh, M.A.R.; Liu, X.; Matsumoto, T.; Morita, K.; Guo, L.; Suzuki, T.; Kamiyama, K. Numerical simulation of the solid particle sedimentation and bed formation behaviors using a hybrid method. Energies; 2020; 13, 5018. [DOI: https://dx.doi.org/10.3390/en13195018]

16. Phan, L.H.S.; Ngo, P.M.; Miura, R.; Tasaki, Y.; Matsumoto, T.; Liu, W.; Morita, K. Self-leveling behavior of mixed solid particles in cylindrical bed using gas-injection method. J. Nucl. Sci. Technol.; 2019; 56, pp. 111-122. [DOI: https://dx.doi.org/10.1080/00223131.2018.1531077]

17. Kudinov, P.; Karbojian, A.; Tran, C.T.; Villanueva, W. Agglomeration and size distribution of debris in DEFOR-A experiments with Bi2O3–WO3 corium simulant melt. Nucl. Eng. Des.; 2013; 263, pp. 284-295. [DOI: https://dx.doi.org/10.1016/j.nucengdes.2013.06.011]

18. Johnson, T.R.; Pavlik, J.R.; Baker, L., Jr. Postaccident Heat Removal: Large-Scale Molten-Fuel-Sodium Interaction Experiments; ANL-75-12 Argonne National Laboratory: Argonne, IL, USA, 1975; [DOI: https://dx.doi.org/10.2172/4228270]

19. Mitchell, G.W.; Ottinger, C.A.; Meister, H. D10 Experiment: Coolability of UO₂ Debris in Sodium with Downward Heat Removal; NUREG/CR-4055 Sandia National Laboratories: Albuquerque, NM, USA, 1984; [DOI: https://dx.doi.org/10.2172/6112398]

20. Magallon, D.; Hohmann, H.; Schins, H. Pouring of 100-kg-scale molten UO₂ into sodium. Nucl. Technol.; 1992; 98, pp. 79-90. [DOI: https://dx.doi.org/10.13182/NT92-A34652]

21. Yakush, S.E.; Konovalenko, A.; Basso, S.; Kudinov, P. Effect of particle spreading on coolability of ex-vessel debris bed. Proceedings of the 16th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-16); Chicago, IL, USA, 30 August–4 September 2015; pp. 1210-1222.

22. Basso, S.; Konovalenko, A.; Kudinov, P. Empirical closures for particulate debris bed spreading induced by gas–liquid flow. Nucl. Eng. Des.; 2016; 297, pp. 19-25. [DOI: https://dx.doi.org/10.1016/j.nucengdes.2015.10.016]

23. Robinson, D.A.; Friedman, S.P. Observations of the effects of particle shape and particle size distribution on avalanching of granular media. Phys. A Stat. Mech. Appl.; 2002; 311, pp. 97-110. [DOI: https://dx.doi.org/10.1016/S0378-4371(02)00815-4]

24. Zhang, B.; Harada, T.; Hirahara, D.; Matsumoto, T.; Morita, K.; Fukuda, K.; Yamano, H.; Suzuki, T.; Tobita, Y. Self-leveling onset criteria in debris beds. J. Nucl. Sci. Technol.; 2010; 47, pp. 384-395. [DOI: https://dx.doi.org/10.1080/18811248.2010.9711969]

25. Cheng, S.; Tanaka, Y.; Gondai, Y.; Kai, T.; Zhang, B.; Matsumoto, T.; Morita, K.; Fukuda, K.; Yamano, H.; Suzuki, T. et al. Experimental studies and empirical models for the transient self-leveling behavior in debris bed. J. Nucl. Sci. Technol.; 2011; 48, pp. 1327-1336. [DOI: https://dx.doi.org/10.1080/18811248.2011.9711823]

26. Morita, K.; Matsumoto, T.; Nishi, S.; Nishikido, T.; Tagami, H.; Suzuki, T.; Tobita, Y. A new empirical model for self-leveling behavior of cylindrical particle beds. J. Nucl. Sci. Technol.; 2016; 53, pp. 713-725. [DOI: https://dx.doi.org/10.1080/00223131.2015.1126204]

27. Basso, S.; Konovalenko, A.; Yakush, S.E.; Kudinov, P. The effect of self-leveling on debris bed coolability under severe accident conditions. Nucl. Eng. Des.; 2016; 305, pp. 246-259. [DOI: https://dx.doi.org/10.1016/j.nucengdes.2016.05.020]

28. Basso, S.; Konovalenko, A.; Kudinov, P. Effectiveness of the debris bed self-leveling under severe accident conditions. Ann. Nucl. Energy; 2016; 95, pp. 75-85. [DOI: https://dx.doi.org/10.1016/j.anucene.2016.04.048]

29. Tatemoto, Y.; Mawatari, Y.; Okazaki, S.; Noda, K. Separation of solid particles by density difference in a liquid–solid fluidized bed. J. Chem. Eng. Jpn.; 2005; 38, pp. 264-270. [DOI: https://dx.doi.org/10.1252/jcej.38.264]

30. Saha, K.; Agarwal, A.K.; Ghosh, K.; Som, S. Two-Phase Flow for Automotive and Power Generation Sectors; Springer: Singapore, 2019; [DOI: https://dx.doi.org/10.1007/978-981-13-3256-2]

31. Cundall, P.A.; Strack, O.D.L. A discrete numerical model for granular assemblies. Geotechnique; 1979; 29, pp. 47-65. [DOI: https://dx.doi.org/10.1680/geot.1979.29.1.47]

32. Siemens Digital Industries Software. Simcenter STAR-CCM+ User Guide, Version 2020.3; Siemens Digital Industries Software: Plano, TX, USA, 2020.

33. Nagaya, Y.; Okumura, K.; Sakurai, T.; Mori, T. MVP/GMVP; Version 3; General purpose Monte Carlo codes for neutron and photon transport calculations based on continuous energy and multigroup methods; JAEA-DATA/CODE-2016-018 Japan Atomic Energy Agency: Ibaraki, Japan, 2017; [DOI: https://dx.doi.org/10.11484/jaea-data-code-2016-018]

34. Smirnov, A.; Bogdanova, E.; Pugachev, P.; Ternovykh, M.; Saldikov, I.; Tikhomirov, G.; Takezawa, H.; Nishiyama, J.; Obara, T. Monte Carlo codes benchmarking on sub-critical fuel debris particles system for neutronic analysis. J. Nucl. Sci. Technol.; 2022; 59, pp. 34-43. [DOI: https://dx.doi.org/10.1080/00223131.2021.1950068]

35. Li, C.Y.; Wang, K.; Pellegrini, M.; Erkan, N.; Okamoto, K. Numerical simulation and validation of debris bed self-leveling behavior with mixed-density particles using CFD-DEM coupling algorithm. Nucl. Technol.; 2021; 208, pp. 843-859. [DOI: https://dx.doi.org/10.1080/00295450.2021.1973181]

36. Razavi, F.; Komrakova, A.; Lange, C.F. CFD–DEM simulation of sand-retention mechanisms in slurry flow. Energies; 2021; 14, 3797. [DOI: https://dx.doi.org/10.3390/en14133797]

37. Ranz, W.E.; Marshall, W.R. Evaporation from drops—Parts I. Chem. Eng. Prog.; 1952; 48, pp. 141-146.

38. Etoh, M.; Kamishima, Y.; Okamura, S.; Watanabe, O.; Ohyama, K.; Negishi, K.; Kotake, S.; Sakamoto, Y.; Kamide, H. Conceptual Design Study of JSFR (2)—Reactor System. Proceedings of the International Conference on Fast Reactors and Related Fuel Cycles: Challenges and Opportunities (FR09); Kyoto, Japan, 7–11 December 2009; IAEA-CN-176/08-11FP

39. Fink, J.K.; Leibowitz, L. Thermodynamic and Transport Properties of Sodium Liquid and Vapor; ANL/RE-95/2 Argonne National Laboratory: Argonne, IL, USA, 1995; [DOI: https://dx.doi.org/10.2172/94649]

40. Naganuma, M.; Ogawa, T.; Ohki, S.; Mizuno, T.; Kotake, S. Minor actinide-bearing oxide fuel core design study for the JSFR. Nucl. Technol.; 2010; 170, pp. 170-180. [DOI: https://dx.doi.org/10.13182/NT10-A9455]

41. Kasahara, N. Fast Reactor System Design; Springer: Singapore, 2017; [DOI: https://dx.doi.org/10.1007/978-981-10-2821-2]

42. Peckner, D.; Bernstein, I.M. Handbook of Stainless Steels; McGrawHill: New York, NY, USA, 1977.

43. Marchetti, M.; Laux, D.; Seibert, A.; Bottomley, P.D.; Wiss, T.; Rondinella, V.V.; Despaux, G. Elastic properties of severely degraded fuels. J. Nucl. Mater.; 2020; 529, 151918. [DOI: https://dx.doi.org/10.1016/j.jnucmat.2019.151918]

44. Carbajo, J.J.; Yoder, G.L.; Popov, S.G.; Ivanov, V.K. A review of the thermophysical properties of MOX and UO₂ fuels. J. Nucl. Mater.; 2001; 299, pp. 181-198. [DOI: https://dx.doi.org/10.1016/S0022-3115(01)00692-4]

45. Bachrata, A.; Bertrand, F.; Marie, N.; Serre, F. A Comparative Study on Severe Accident Phenomena Related to Melt Progression in Sodium Fast Reactors and Pressurized Water Reactors. J. Nucl. Eng. Radiat. Sci.; 2021; 7, 030801. [DOI: https://dx.doi.org/10.1115/1.4047921]

46. Lamarsh, J.R.; Baratta, A.J. Introduction to Nuclear Engineering; 3rd ed. Prentice Hall: Upper Saddle River, NJ, USA, 2001.

47. Hoelzer, D.T.; Pint, B.A.; Wright, I.G. A microstructural study of the oxide scale formation on ODS Fe–13Cr steel. J. Nucl. Mater.; 2000; 283, pp. 1306-1310. [DOI: https://dx.doi.org/10.1016/S0022-3115(00)00248-8]

48. Shibata, K.; Iwamoto, O.; Nakagawa, T.; Iwamoto, N.; Ichihara, A.; Kunieda, S.; Chiba, S.; Furutaka, K.; Otuka, N.; Ohsawa, T. et al. JENDL-4.0: A new library for nuclear science and engineering. J. Nucl. Sci. Technol.; 2011; 48, pp. 1-30. [DOI: https://dx.doi.org/10.1080/18811248.2011.9711675]

49. Freiría López, M.; Buck, M.; Starflinger, J. Neutronic modeling of debris beds for a criticality evaluation. Ann. Nucl. Energy; 2019; 130, pp. 164-172. [DOI: https://dx.doi.org/10.1016/j.anucene.2019.02.034]