**Sedimentation and Bed Formation Behaviors Using a Hybrid Method**

**Md Abdur Rob Sheikh 1, Xiaoxing Liu 1, Tatsuya Matsumoto 1, Koji Morita 1,\*, Liancheng Guo 2, Tohru Suzuki <sup>3</sup> and Kenji Kamiyama <sup>3</sup>**


Received: 1 August 2020; Accepted: 21 September 2020; Published: 24 September 2020

**Abstract:** 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.

**Keywords:** hybrid simulation method; multi-fluid model; discrete element method, sedimentation, bed formation

#### **1. 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–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–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.

#### **2. 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–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]:

$$\frac{\partial \overline{\rho}\_m}{\partial t} + \nabla \cdot \left(\overline{\rho}\_m \, \overline{\boldsymbol{\upsilon}}\_q\right) = -\Gamma\_{m\prime} \tag{1}$$

$$\begin{split} \frac{\partial \overline{\rho}\_{q} \overrightarrow{\boldsymbol{v}}\_{q}}{\partial t} + \sum\_{m \in q} \nabla \cdot \left( \overline{\rho}\_{m} \overrightarrow{\boldsymbol{v}}\_{q} \overrightarrow{\boldsymbol{v}}\_{q} \right) + \alpha\_{q} \nabla p - \overline{\rho}\_{q} \overrightarrow{\boldsymbol{g}} + \mathcal{K}\_{q\mathcal{S}} \overrightarrow{\boldsymbol{v}}\_{q} - \sum\_{q'} \mathcal{K}\_{qq'} \left( \overrightarrow{\boldsymbol{v}}\_{q'} - \overrightarrow{\boldsymbol{v}}\_{q} \right) - \overrightarrow{\mathcal{V}\mathcal{M}}\_{q} \\ = -\sum\_{q'} \Gamma\_{qq'} \left[ \mathcal{H} \{ \Gamma\_{qq'} \} \overrightarrow{\boldsymbol{v}}\_{q} - \mathcal{H} \{ -\Gamma\_{qq'} \} \overrightarrow{\boldsymbol{v}}\_{q'} \right], \end{split} \tag{2}$$

$$\begin{split} \frac{\partial \overline{\rho}\_{M} c\_{M}}{\partial t} + \sum\_{m \in M} \nabla \cdot \left( \overline{\rho}\_{m} c\_{m} \overline{\boldsymbol{\upsilon}}\_{q} \right) + p \left[ \frac{\partial \boldsymbol{a}\_{M}}{\partial t} + \nabla \cdot \left( \boldsymbol{a}\_{M} \overline{\boldsymbol{\upsilon}}\_{q} \right) \right] - \frac{\overline{\rho}\_{M}}{\overline{\rho}\_{q}} \left[ \sum\_{q'} K\_{qq'} \left( \overline{\boldsymbol{\upsilon}}\_{q} - \overline{\boldsymbol{\upsilon}}\_{q'} \right) \left( \overline{\boldsymbol{\upsilon}}\_{q} - \overline{\boldsymbol{\upsilon}}\_{q'} \right) \right] \\ + K\_{qS} \overline{\boldsymbol{\upsilon}}\_{q'} \left( \overline{\boldsymbol{\upsilon}}\_{q} - \overline{\boldsymbol{\upsilon}}\_{qS} \right) - \overline{V} \overline{M}\_{q'} \left( \overline{\boldsymbol{\upsilon}}\_{q} - \overline{\boldsymbol{\upsilon}}\_{GL} \right) = Q\_{N,M} + Q\_{M,M} + Q\_{H,M} \end{split} \tag{3}$$

where subscripts *m*, *q* and *M* denote the density, velocity and energy component, respectively; *t* is the time; ρ*<sup>m</sup>* and Γ*<sup>m</sup>* are the macroscopic density and total mass transfer rate per unit volume from the density component *<sup>m</sup>*, respectively; *<sup>v</sup> <sup>q</sup>* is the vector of the velocity field *<sup>q</sup>*; *<sup>p</sup>* denotes the pressure, *g* the gravitational acceleration, *KqS* the momentum exchange function between the velocity field *q* and structure component, *Kqq* the momentum exchange function between the velocity field *q* and *q* , *VMq* 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 *eM* and α*<sup>M</sup>* are the specific internal energy and volume fractions of component *M*, respectively; *QN*, *QM* and *QH* 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

$$m\_i \frac{d^2 \overrightarrow{r\_i}}{dt} = \sum\_j \overrightarrow{F\_{ij}} + \overrightarrow{F\_i} + \overrightarrow{F\_i} \,\prime \tag{4}$$

$$
\overrightarrow{\dot{v}\_i} = \frac{d\overrightarrow{r}\_i}{dt}\prime\tag{5}
$$

$$
\overrightarrow{I}\_i \frac{d^2 \overrightarrow{\theta\_i}}{dt^2} = \sum\_j \overrightarrow{A}\_{ij\prime} \tag{6}
$$

$$
\overrightarrow{\alpha\_i} = \frac{d\overrightarrow{\theta\_i}}{dt},
\tag{7}
$$

where *mi* is the mass of particle *i*; *ri* and *<sup>v</sup> <sup>i</sup>* are the vectors of its position and velocity, respectively; <sup>θ</sup>*<sup>i</sup>* and <sup>ω</sup>*<sup>i</sup>* are the vectors of its angular displacement and angular velocity, respectively; *F c ij* is the contact force of particle *<sup>i</sup>* with neighboring particle *<sup>j</sup>* or wall; *F f <sup>i</sup>* the total solid–fluid interaction force on particle *i*; *F g <sup>i</sup>* is the gravitational force; *<sup>I</sup> <sup>i</sup>* and *Aij* and ω*<sup>i</sup>* 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

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.

#### **3. Verification Experiment and Numerical Simulations**

#### *3.1. 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 *Dc* 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 *dn* of 40 mm. The initial water level is 480 mm, and the exit of the nozzle is located at a height *Nh* 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.

**Figure 1.** Schematic illustration of the setup for the particle sedimentation experiments.

#### *3.2. 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.



**Figure 2.** Schematic illustration of the setups for the simulation (color bar: volume fraction of water phase): (**a**) initial arrangement of solid particles; (**b**) arrangement of dummy particles.

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−<sup>5</sup> 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.** Values of model parameters used in the DEM calculation.

#### **4. Results and Comparisons**

#### *4.1. Homogeneous Particles*

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 Al2O3 particles during their fall, the corresponding lateral dispersion of the Al2O3 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 Al2O3 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 Al2O3 particles with *dp* = 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 *dp* = 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.

**Figure 3.** Falling and sedimentation behavior of homogeneous particles with *dp* = 6 mm: (**a**) Case 1; (**b**) Case 2.

**Figure 4.** Side-view height of the homogeneous particle bed in the radial direction: (**a**) Case 1; (**b**) Case 2.

#### *4.2. 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 Al2O3 particles are dispersed widely during their fall in the pool. Similar falling behavior is observed in Case M2, for which larger falling Al2O3 particles with *dp* = 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.

**Figure 5.** Falling and sedimentation behaviors of mixed particles with different properties: (**a**) Case M1; (**b**) Case M2; (**c**) Case M3; (**d**) Case M4; (**e**) Case M5.

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.

#### *4.3. 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, photographs are shown of the cross-sections at four axial positions for binary mixtures of Al2O3 and SS particles with different particle sizes in equal and different volume proportions. Here, the mixtures are (a) Al2O3 and SS particles with *dp* = 6 mm in equal volume proportions, (b) Al2O3 particles with *dp* = 4 and 6 mm in equal volume proportions, (c) Al2O3 and SS particles with *dp* = 6 mm in a volume proportion of 1:3 and (d) Al2O3 and SS particles with *dp* = 6 mm in a volume proportion of 3:1. The axial positions indexed as 1, 2, 3 and 4 are the top, middle and bottom of the conical bed mound, and the middle of the cylindrical bed basement, respectively. For Case M1 (Figure 7a), the heavier SS particles concentrate around the bed center at the axial four positions, whereas the light Al2O3 particles disperse largely during their fall and after their impact on the bed mound. This behavior was also observed in Cases M4 and M5 with mixtures of Al2O3 and SS particles (Figure 7c,d). A similar particles distribution behavior occurs in Case M2 (Figure 7b), where larger or heavier Al2O3 particles are seen more in the periphery of the bed, although the behaviors of the dispersion and distribution do not substantially contrast those of Cases M1, M4 and M5. Nevertheless, the observed particles distribution behaviors were reproduced reasonably well by the present simulation in all cases.

The apparent areas of Al2O3 and SS particles mixed with different sizes and densities on the cross-sections at the four axial positions are shown in Figure 8, specifically (a) Case M1: Al2O3 and SS particles with *dp* = 6 mm, (b) Case M2: Al2O3 particles with *dp* = 4 and 6 mm, (c) Case M4: volume mixing ratio 1:3, Al2O3 and SS particles with *dp* = 6 mm and (d) Case M5: volume mixing ratio 3:1, Al2O3 and SS particles with *dp* = 6 mm. The axial positions of the cross-sections are the same as those in Figure 7. In this figure, the proportions of Al2O3 and SS particles were defined as the ratios of apparent areas of Al2O3 and SS particles to the cross-sectional area of the bed. In the experiments, the apparent areas of the particles in different colors were measured from the cross-sectional shot images of the bed using an image analysis tool. The ratios of the two different particles obtained in the simulations and the experiments are indicated in the upper and lower sides, respectively, at each axial position. A strong agreement between the simulated results of the particles distributed area at different axial positions with the experimental observations was established (Figure 8a–d). In Case M1, the central and lower parts of the bed are mainly occupied by SS particles and this behavior is seen also in Case M4. However, in Case M2, the ratio of the smaller particles in the bed increases in the height direction. Moreover, the lighter Al2O3 particles are seen in all parts of the bed in Case M5. Overall, the particles distribution behavior of the lighter and larger Al2O3 particles observed in Cases M1, M2, M4 and M5 is reproduced quantitatively in the present simulations.

**Figure 6.** Radial variation of bed height for mixed particles with different properties in the radial direction: (**a**) Case M1; (**b**) Case M2; (**c**) Case M3; (**d**) Case M4; (**e**) Case M5.

**Figure 7.** Top-view of mixed particle bed at four cross-sectional positions: (**a**) Case M1; (**b**) Case M2; (**c**) Case M4; (**d**) Case M5.

**Figure 8.** Particle distribution at four axial positions in the mixed particle bed: (**a**) Case M1; (**b**) Case M2; (**c**) Case M4; (**d**) Case M5.

#### **5. 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.

**Author Contributions:** Conceptualization, K.M., T.S. and K.K.; methodology, K.M., T.S. and K.K.; software, X.L. and L.G.; validation, T.M.; formal analysis, M.A.R.S.; investigation, M.A.R.S. and T.M.; resources, K.M. and T.M.; data curation, M.A.R.S.; writing—original draft preparation, M.A.R.S.; writing—review and editing, K.M.; visualization, M.A.R.S.; supervision, K.M. and T.M.; project administration, K.M., T.S. and K.K.; funding acquisition, K.M., T.S. and K.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Japan Society for the Promotion of Science (JSPS) KAKENHI [Grant No. 16K06960] and Japan Atomic Energy Agency.

**Acknowledgments:** The authors would like to thank R. Kawata and Y. Ohara for assistance with the numerical simulations. We also thank E. Son for assistance with the experiments. Finally, we are grateful to Richard Haase, from Edanz Group (www.edanzediting.com/ac) for editing a draft of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
