Numerical Simulation of the Solid Particle Sedimentation and Bed Formation Behaviors Using a Hybrid Method

In the safety analysis of sodium-cooled fast reactors, numerical simulations of various thermal-hydraulic phenomena with multicomponent and multiphase flows in core disruptive accidents (CDAs) are regarded as particularly difficult. In the material relocation phase of CDAs, core debris settle down on a core support structure and/or an in-vessel retention device and form a debris bed. The bed’s shape is crucial for the subsequent relocation of the molten core and heat removal capability as well as re-criticality. In this study, a hybrid numerical simulation method, coupling the multi-fluid model of the three-dimensional fast reactor safety analysis code SIMMER-IV with the discrete element method (DEM), was applied to analyze the sedimentation and bed formation behaviors of core debris. Three-dimensional simulations were performed and compared with results obtained in a series of particle sedimentation experiments. The present simulation predicts the sedimentation behavior of mixed particles with different properties as well as homogeneous particles. The simulation results on bed shapes and particle distribution in the bed agree well with experimental measurements. They demonstrate the practicality of the present hybrid method to solid particle sedimentation and bed formation behaviors of mixed as well as homogeneous particles.


Introduction
In the latter stage of an unprotected loss-of-flow (ULOF) event, which is a representative initiator of a core disruptive accident (CDA) in sodium-cooled fast reactors (SFRs), molten core materials may relocate into a lower sodium plenum due to a gravity-driven discharge through potential paths such as control rod guide tubes. During this material relocation phase, the molten core materials are rapidly quenched and fragmented in the subcooled sodium and become solidified into smaller-sized particles [1]. The discharged particles, called debris, which are roughly 300 µm in Sauter mean diameter in SFR [2], settle down on a core support structure and/or an in-vessel retention (IVR) device such as a multi-layered debris tray installed in the bottom region of the reactor vessel and form debris beds [3].
In the safety analysis of CDAs, the sedimentation behavior of these particles is a crucial issue as it has a significant influence on heat removal from the fuel debris with decay heat. Numerous studies have focused on debris bed characteristics [4,5], but the shape of the debris bed is chosen arbitrarily, e.g., cylindrical [6], conical or Gaussian-shaped [7], heap-like [8,9], conical [10] or homothetically piled [11]. Moreover, the debris bed was assumed homogeneous and uniformly distributed above the lower-most level of the reactor pool as supposed in a preceding study of the cooling capability of a debris bed in a reactor accident [4].
Experimental investigations with simulant materials composed of homogeneous solid particles were conducted recently [12,13] and, in connection with this study, an exploration was performed in our previous studies for mixtures of particles with different properties [14,15]. These experimental analyses discovered the nature of the shape characteristics of the particle bed and its dependence on some important factors, such as particle density and diameter as well as the diameter and length of the nozzle, which is used to discharge solid particles into the water pool [12][13][14][15]. As the experimental process consumes a high cost and long duration, the evaluation of such sedimentation behaviors is performed knowing the limited scope of reproducibility. Although the high cost and long duration of such experiments can be reduced, it is difficult to imitate the reactor conditions rigorously. However, without experimental data, a justifiable extrapolation to reactor conditions would be tough. This is because numerical analysis, which reflects findings obtained by experiments with the addition of possible reproducibility of the reactor conditions, enables us some extrapolation, even though the experimental scope is limited. Consequently, we require the numerical study on sedimentation behavior using solid particles as simulated debris.
The first practical tool, the SIMMER-II code, was developed at the Los Alamos National Laboratory (LANL). This code was used in many experimental and reactor analyses [16,17]. Subsequently, SIMMER-III was introduced by the Japan Nuclear Cycle Development Institute (JNC), presently called the Japan Atomic Energy Agency (JAEA). The SIMMER-III code is a two-dimensional, multi-velocity-field, multiphase, multicomponent, Eulerian, fluid dynamics combined with a spaceand energy-dependent neutron kinetics model [18]. The major drawback of SIMMER-III is its two-dimensional treatment, although the computational technology of SIMMER-III has improved the safety evaluation of liquid-metal-cooled fast reactors (LMFRs). The lack of dimensionality in SIMMER-III requires certain conservatism because of large uncertainties in accident analysis. Consideration of the three-dimensional distribution of core materials including neutron absorbers would reasonably simulate the reactivity effect, which could lead to severe re-criticality, caused by the relocation of disrupted core materials [19,20]. Although we need large computational resources in large-scale three-dimensional simulations of safety analyses, the parallelization technology will reduce the computational load.
Reasonable numerical simulations of the transient debris behavior involving the thermal-hydraulic phenomena of the surrounding fluid phases require a comprehensive computational tool, considering a complicated multiphase mechanism for the debris bed cooling capability. A reliable and trustworthy tool, the reactor safety analysis code, SIMMER-IV [19,20], has been applied to key phenomena such as fuel discharge and relocation [21] and pool sloshing [22], as well as safety analyses of the transition phase [19,20] and the post-disassembly expansion phase [23], in CDAs of an LMFR successfully. It is a three-dimensional multi-velocity-field, multiphase, multicomponent, Eulerian, fluid dynamics code coupled with a neutron kinetics model and a fuel pin model. This code has recently been productively applied to simulations of critical thermal-hydraulic phenomena in CDAs as well as to reactor safety analysis. However, mechanical interactions among solid particles and discrete phase characteristics of solid particles are not modeled directly in the associated fluid dynamics calculations in the simulations of multiphase flows with rich solid particles. Appropriate calculations of the interaction between the fluid phases and particles, as well as the interactions between particles themselves, are necessary for simulating such phenomena reasonably. A discrete method can predict the solid particle interactions and motions directly based on the solid-phase continuity at a microscopic level. In principle, with proper initial and boundary conditions, Newton's equations of motions are solved together. Unlike conventional grid methods, it is not necessary to assume uniform constituency and constitutive relations for discrete solid particles [24]. Cundall and Strack [25] advanced and introduced the discrete element method (DEM), which is the most successful numerical approach among this category of numerical methods. Multi-body collisions can be calculated accurately with an explicit force model. Besides, DEM gives transient trajectories and velocities of particles as local information. Moreover, the fluid-particle interactions are calculated by coupling the DEM with computational fluid dynamics methods using drag force terms straightforwardly [26].
A theoretical model was established to explore the characteristics of the dense particle bed in simulations of the decay heat removal from the fuel debris formed in the lower plenum of SFRs [27]. The inter-particle interaction is anticipated to be a substantial phenomenon in the dense particle bed. Inter-particle collisions and contacts were demonstrated to be influential. In addition, this model is able to explain that characteristics of the debris movement initiated by sodium vapor flow in debris beds.
Guo et al. [28][29][30][31] combined DEM with the multi-fluid model and developed a hybrid computational method based on the theoretical background introduced above. It is expected that this numerical method reproduces the particle behavior involved in the thermal-hydraulic phenomena reasonably. A semi-implicit time factorization approach solves the governing equations of the fluid phases, whereas DEM calculates particle movements through the coupling algorithm. In their governing equations, the multiphases are then coupled via drag terms explicitly. The fundamental applicability of the developed hybrid method has been validated by simulating typical gas-solid fluidized beds [28], the self-leveling of particle beds in a rectangular pool [29] and gas-liquid particle three-phase flows [30] in two-dimensional systems. Recently, a three-dimensional numerical simulation was performed for the self-leveling of the particle beds in a cylindrical pool, and a fundamental validation of the developed hybrid method was demonstrated successfully [31]. In addition, a DEM-based numerical study on the sedimentation behavior of solid particles in two-dimensional systems was accomplished [32]. The primitive numerical simulation study was not coupled with the fluid flow and particle dynamics. Therefore, a numerical analysis in a three-dimensional system that covers both a mixture of particles as well homogeneous ones is essential.
The objective of the present numerical study is to demonstrate the practicality of the hybrid method to particle sedimentation and bed formation behaviors of mixed particles as well as homogeneous ones. In the present study, three-dimensional simulations were performed and compared with results obtained in a series of particle sedimentation experiments in which gravity-driven solid particles [12,15] are discharged into a quiescent cylindrical water pool. Although the experiments using simulant materials did not reproduce the particle sedimentation and bed formation behaviors that are expected to occur under reactor accident conditions, the hybrid method will be validated for fundamental characteristics of the behaviors, which were measured in the experiments performed under controlled conditions.

Mathematical Treatment
In the present study, we performed three-dimensional calculations for solid particle sedimentation and bed formation behaviors of mixed as well as homogeneous particles using the hybrid computational method, which combines DEM with the multi-fluid model of the SIMMER-IV code. A detained explanation of the method is given by Guo et al. [28][29][30][31].
In the multi-fluid model of the SIMMER-IV code, the governing equations of multi-fluid phases are the conservation equations of mass, momentum and energy, in terms of the local mean variables over a computational cell, in abbreviated form [19]: ∂ρ M e M ∂t where subscripts m, q and M denote the density, velocity and energy component, respectively; t is the time; ρ m and Γ m are the macroscopic density and total mass transfer rate per unit volume from the density component m, respectively; v q is the vector of the velocity field q; p denotes the pressure, g the gravitational acceleration, K qS the momentum exchange function between the velocity field q and structure component, K qq the momentum exchange function between the velocity field q and q , VM q the virtual mass term of the gas phase, H(x) the Heaviside unit function, Γ qq the mass transfer rate from component q and q , and e M and α M are the specific internal energy and volume fractions of component M, respectively; Q N , Q M and Q H are, respectively, the nuclear heating rate, the rate of energy interchange due to the mass transfer and other heat transfer rates. The subscript GL indicates terms such as the averaged velocity relevant at interfaces between the vapor and liquid.
In the DEM calculations, the solid phase is represented by a discrete phase. Newton's law [24] is used to describe the translational motion of a particle as where m i is the mass of particle i; r i and v i are the vectors of its position and velocity, respectively; θ i and ω i are the vectors of its angular displacement and angular velocity, respectively; F c ij is the contact force of particle i with neighboring particle j or wall; F f i the total solid-fluid interaction force on particle i; F g i is the gravitational force; I i and A ij and ω i are the moment of inertia, torque and angular velocity, respectively.
In the present DEM, which assumes that particles are spheres, a viscoelastic contact model [33] is applied to calculate the contact forces between the particles, as well as between particles and the wall. The DEM calculation is coupled with the fluid dynamics one of the SIMMER-IV code, which is based on a time factorization time-splitting approach, through the terms of solid-fluid interactions in the momentum conservation Equation (2). Time-step sizes for the fluid dynamics and DEM calculations are controlled independently. In the SIMMER-IV code, the Courant condition, the optimum pressure iteration condition and the excessive vaporization/condensation iteration condition are used mainly Energies 2020, 13, 5018 5 of 15 to optimize the time-step size of fluid dynamics calculations. For the DEM calculation, in which velocities and positions of particles are updated based on an explicit scheme, a large time-step size is limited so as to prevent unphysical disturbances of the particle in the particle-fluid interactions. In addition, the smallest time period, by which two particles reach their maximum overlap from the initial contact in the collision, is considered. In general, the latter time-step size is much smaller in the DEM calculations.

Particle Sedimentation Experiments
A series of particle sedimentation experiments was performed for gravity-driven discharges of solid particles into a quiescent cylindrical water pool [12,15]. For the experimental setup ( Figure 1), a transparent cylindrical tank filled with water with an inner diameter D c of 375 mm and a height of 1040 mm is the main apparatus. Solid particles, which are kept in an overhead hopper initially, were poured into the water tank though a nozzle having an inner diameter d n of 40 mm. The initial water level is 480 mm, and the exit of the nozzle is located at a height N h of 473 mm from the bottom of the tank. In the experiments, the solid particles discharged from the nozzle fall into the water pool, and are deposited on the bottom of the tank finally forming a particle bed with a conical shape. The particle materials include Al 2 O 3 and SS. For homogeneous spherical particles, their bulk volume is 5 L, whereas the two types of spherical particles with a bulk volume of 2.5 L are mixed up as a binary mixture of particles with different properties.

Particle Sedimentation Experiments
A series of particle sedimentation experiments was performed for gravity-driven discharges of solid particles into a quiescent cylindrical water pool [12,15]. For the experimental setup (Figure 1), a transparent cylindrical tank filled with water with an inner diameter of 375 mm and a height of 1040 mm is the main apparatus. Solid particles, which are kept in an overhead hopper initially, were poured into the water tank though a nozzle having an inner diameter of 40 mm. The initial water level is 480 mm, and the exit of the nozzle is located at a height of 473 mm from the bottom of the tank. In the experiments, the solid particles discharged from the nozzle fall into the water pool, and are deposited on the bottom of the tank finally forming a particle bed with a conical shape. The particle materials include Al2O3 and SS. For homogeneous spherical particles, their bulk volume is 5 L, whereas the two types of spherical particles with a bulk volume of 2.5 L are mixed up as a binary mixture of particles with different properties.

Simulation Conditions
We calculated several cases of particle sedimentation experiments using binary mixtures of spherical Al2O3 and SS particles with different densities or different sizes as well as homogeneous particles. In these arrangements, we considered two cases with homogeneous particles and five cases with binary mixed particles. In the binary mixture, we chose three cases with different densities and sizes but set an equal volume mixing ratio. In addition, we took the setup as two cases with binary mixed particles of different volume mixing ratios. However, in this simulation process, we selected particles of larger sizes to avoid prolonged calculations times. The simulation conditions and particle properties of the selected cases are listed in Table 1. Two cases, labeled Cases 1 and 2, correspond to homogeneous solid particles and five types of binary mixtures, labeled Cases M1-M5, were chosen for the simulations presented. Figure 2 shows the initial position of the solid particles and fluid phases in the computational domain of Case M1. The number of computational cells for fluid phases is 25 × 25 × 25 in a X-Y-Z geometry, and the width of each cell is 24 mm. For solid particles, which filled the hopper initially, 31,568 DEM particles were used in simulations in the same bulk volume as the experiment.

Simulation Conditions
We calculated several cases of particle sedimentation experiments using binary mixtures of spherical Al 2 O 3 and SS particles with different densities or different sizes as well as homogeneous particles. In these arrangements, we considered two cases with homogeneous particles and five cases with binary mixed particles. In the binary mixture, we chose three cases with different densities and sizes but set an equal volume mixing ratio. In addition, we took the setup as two cases with binary mixed particles of different volume mixing ratios. However, in this simulation process, we selected particles of larger sizes to avoid prolonged calculations times. The simulation conditions and particle properties of the selected cases are listed in Table 1. Two cases, labeled Cases 1 and 2, correspond to homogeneous solid particles and five types of binary mixtures, labeled Cases M1-M5, were chosen for the simulations presented. Figure 2 [34]. Smaller artificial values (MPa scale or smaller) are applied to maintain the time-step sizes of order 10 −5 s with sufficient accuracy. The particle hopper, which initially retains the solid particles, and the cylindrical tank are shaped by dummy DEM particles fixed in space to consider interactions with moving solid particles. The simulated wall surfaces of the hopper and the tank are not smooth and do not represent the real ones exactly. Although particle motions are affected in varying degrees by these artificially rough surfaces, especially of the inner wall of the nozzle, it is expected that the particle sedimentation and bed formation behaviors are not affected largely by the rough surface of the nozzle because the bulk particle flow in the nozzle is dense during the particle injection.  Table 2 lists the model parameter values used in the DEM calculation. The time-step size was very small, as explained in Chapter 2, if real values (GPa scale) of Young's modulus E and the modulus of rigidity G are used in the calculation. To solve this problem, we refer the reader to [34]. Smaller artificial values (MPa scale or smaller) are applied to maintain the time-step sizes of order 10 −5 s with sufficient accuracy. The particle hopper, which initially retains the solid particles, and the cylindrical tank are shaped by dummy DEM particles fixed in space to consider interactions with moving solid particles. The simulated wall surfaces of the hopper and the tank are not smooth and do not represent the real ones exactly. Although particle motions are affected in varying degrees by these artificially rough surfaces, especially of the inner wall of the nozzle, it is expected that the particle sedimentation and bed formation behaviors are not affected largely by the rough surface of the nozzle because the bulk particle flow in the nozzle is dense during the particle injection.   Figure 3 shows visual snapshots of homogeneous particles falling in the water pool and their sedimentation behaviors at 1.0, 9.0 and 17 s after initiating particle injection for Case 1 and similarly at 1.0, 8.0 and 17 s for Case 2. In Case 1, the SS particles form a jet-like particle flow in the pool. From the figures, the lateral dispersion of the SS particles is most likely to be underestimated in the present simulations. This may be caused by the lack of a turbulence model in the present simulations. Nevertheless, although a wide-spreading behavior is observed in Case 2 for the lighter Al 2 O 3 particles during their fall, the corresponding lateral dispersion of the Al 2 O 3 particles in the water pool is underestimated in the simulation. Moreover, for Case 2, the turbulence flow induced by the falling particles may influence the particle motion and enhance the dispersal of the lighter Al 2 O 3 particles, however the multi-fluid model used in the present simulation does not consider turbulence effects in the water, which will have an impact on the motion of the solid particles. In addition, in another experiment using Al 2 O 3 particles with d p = 6 mm in air, less lateral spreading of the particles is observed because turbulence effects on the solid particles were relatively small. Therefore, the difference in observations of the particle falling behavior may be caused by the lack of a turbulence model in the present simulation. Nevertheless, for both Cases 1 and 2, the shape of the bed formed after the completion of particle sedimentation is comparable in both simulations and experiments. The radial variation of the bed height of the homogeneous particles with d p = 6 mm is presented in Figure 4; specifically, the experiment results, which did not show strong axial asymmetry in bed shapes, are compared with (a) the calculated radial bed height of SS particles, Case 1, and (b) the calculated radial bed height of Al 2 O 3 particles, Case 2. In this figure, the side-view height was measured from the side-view shot image of the particle bed. A strong agreement is seen between the two sets of results for the side-view profile of both SS and Al 2 O 3 particles. This is because even if the lateral spreading of the particles is Energies 2020, 13, 5018 8 of 15 smaller, particles settling higher up on the bed will tumble down the mound, and hence the bed will finally have a similar repose angle. after the completion of particle sedimentation is comparable in both simulations and experiments. The radial variation of the bed height of the homogeneous particles with = 6 mm is presented in Figure 4; specifically, the experiment results, which did not show strong axial asymmetry in bed shapes, are compared with (a) the calculated radial bed height of SS particles, Case 1, and (b) the calculated radial bed height of Al2O3 particles, Case 2. In this figure, the side-view height was measured from the side-view shot image of the particle bed. A strong agreement is seen between the two sets of results for the side-view profile of both SS and Al2O3 particles. This is because even if the lateral spreading of the particles is smaller, particles settling higher up on the bed will tumble down the mound, and hence the bed will finally have a similar repose angle.   after the completion of particle sedimentation is comparable in both simulations and experiments. The radial variation of the bed height of the homogeneous particles with = 6 mm is presented in Figure 4; specifically, the experiment results, which did not show strong axial asymmetry in bed shapes, are compared with (a) the calculated radial bed height of SS particles, Case 1, and (b) the calculated radial bed height of Al2O3 particles, Case 2. In this figure, the side-view height was measured from the side-view shot image of the particle bed. A strong agreement is seen between the two sets of results for the side-view profile of both SS and Al2O3 particles. This is because even if the lateral spreading of the particles is smaller, particles settling higher up on the bed will tumble down the mound, and hence the bed will finally have a similar repose angle.

Mixed Particles
Visual snapshots of the falling and sedimentation behaviors of mixed particles in the water pool are similarly presented ( Figure 5). A comparison of snapshots between simulation and experimental results is given at three instants in time after particle injection has begun for Cases M1-M5. In Case M1, the bulk particle jet flow is formed mainly by heavy SS particles, and lighter Al 2 O 3 particles are dispersed widely during their fall in the pool. Similar falling behavior is observed in Case M2, for which larger falling Al 2 O 3 particles with d p = 6 mm spread in the pool, whereas smaller ones form a particle jet flow. In Case M3, particles dispersion from the main jet flow is not observed. However, for the particle mixtures with different volume mixing ratios, Cases M4 and M5, the dispersion of the lighter Al 2 O 3 particles still occurs. Although in simulations the lateral dispersion of particles is likely underestimated during their fall compared with observations, the observed falling and sedimentation behaviors are reasonably well reproduced in these five cases.
Energies 2020, 13, 5018 9 of 15 dispersed widely during their fall in the pool. Similar falling behavior is observed in Case M2, for which larger falling Al2O3 particles with = 6 mm spread in the pool, whereas smaller ones form a particle jet flow. In Case M3, particles dispersion from the main jet flow is not observed. However, for the particle mixtures with different volume mixing ratios, Cases M4 and M5, the dispersion of the lighter Al2O3 particles still occurs. Although in simulations the lateral dispersion of particles is likely underestimated during their fall compared with observations, the observed falling and sedimentation behaviors are reasonably well reproduced in these five cases. In Figure 6, comparisons of the radial bed height variation after the completion of particle sedimentation are shown for the five mixed particle cases, in which strong axial asymmetry was not observed in the bed shapes. Although a general agreement between the simulation and experimental results is obtained, the slight discrepancies in bed height may be due to differences in falling and sedimentation behavior. In Figure 6, comparisons of the radial bed height variation after the completion of particle sedimentation are shown for the five mixed particle cases, in which strong axial asymmetry was not observed in the bed shapes. Although a general agreement between the simulation and experimental results is obtained, the slight discrepancies in bed height may be due to differences in falling and sedimentation behavior.

Analysis of Particles Distributions
In Cases M1 and M2, M4 and M5, after the particle sedimentation was completed, top-view pictures of particle beds were taken at four cross-sectional positions by removing particles from the bed using a vacuum suction device. In Figure 7 In Figure 6, comparisons of the radial bed height variation after the completion of particle sedimentation are shown for the five mixed particle cases, in which strong axial asymmetry was not observed in the bed shapes. Although a general agreement between the simulation and experimental results is obtained, the slight discrepancies in bed height may be due to differences in falling and sedimentation behavior.

Analysis of Particles Distributions
In Cases M1 and M2, M4 and M5, after the particle sedimentation was completed, top-view pictures of particle beds were taken at four cross-sectional positions by removing particles from the bed using a vacuum suction device. In Figure 7

Conclusions
In the present validation study of the hybrid numerical simulation method, three-dimensional simulations were performed and compared with results obtained in a series of particle sedimentation experiments of gravity-driven discharges of solid particles, the simulants of core debris, into a quiescent water pool. The governing equations of the multi-fluid phases were calculated using the conventional multi-fluid model, whereas for solid particles, the equations for the translational motion were calculated considering contacts among particles based on DEM. The simulation results on the particle falling behavior, bed height and particle distribution in the bed agree well with the experimental results for binary mixtures of particles with different densities or sizes as well as homogeneous particles. The particles distribution behavior of the lighter or larger particles in the binary mixtures was successfully reproduced using the present hybrid method. Although the lack of a turbulence model may underestimate the dispersion behavior of lighter and larger Al2O3 particles observed during their fall in the pool, the results of the present simulation demonstrate the practicality of the present hybrid method to the solid particle sedimentation and bed formation behaviors of both mixed and homogeneous particles. It might be of interest to consider quantifying the effect of the rough surfaces of the hopper and the tank walls, which were modeled by rather large dummy DEM particles, on the motion of solid particles in future work. Further validation is needed of the present method for a wide range of experimental conditions such as particle size and density, that have a strong influence on particle sedimentation and bed formation behaviors.

Conclusions
In the present validation study of the hybrid numerical simulation method, three-dimensional simulations were performed and compared with results obtained in a series of particle sedimentation experiments of gravity-driven discharges of solid particles, the simulants of core debris, into a quiescent water pool. The governing equations of the multi-fluid phases were calculated using the conventional multi-fluid model, whereas for solid particles, the equations for the translational motion were calculated considering contacts among particles based on DEM. The simulation results on the particle falling behavior, bed height and particle distribution in the bed agree well with the experimental results for binary mixtures of particles with different densities or sizes as well as homogeneous particles. The particles distribution behavior of the lighter or larger particles in the binary mixtures was successfully reproduced using the present hybrid method. Although the lack of a turbulence model may underestimate the dispersion behavior of lighter and larger Al 2 O 3 particles observed during their fall in the pool, the results of the present simulation demonstrate the practicality of the present hybrid method to the solid particle sedimentation and bed formation behaviors of both mixed and homogeneous particles. It might be of interest to consider quantifying the effect of the rough surfaces of the hopper and the tank walls, which were modeled by rather large dummy DEM particles, on the motion of solid particles in future work. Further validation is needed of the present method for a wide range of experimental conditions such as particle size and density, that have a strong influence on particle sedimentation and bed formation behaviors.