An Improved Distinct Element Method for High Packing Fraction Stochastic Media Modeling

Due to the generality and flexibility of Monte Carlo method in geometric modeling, Monte C arlo method plays an important role in accurate simulation of random media. At present, rand om sequential addition method (RSA) and distinct element method (DEM) are more accurate and mature explicit modeling methods. The former approach has the problem of upper limit of packing fraction, which is suitable for stochastic geometry with lower filling rate. DEM m ethod can fill random medium model with packing fraction higher than 60%, but DEM is not suitable for non-contact dispersed particles based on the interaction between particles. There fore, an improved DEM method is proposed to solve the problem of modeling non-contact p articles dispersed in the stochastic media with high packing fraction. The virtual surfaces are constructed outside of the outer layer of particles to make them in contact with each other. T hus, the particle system is suitable for DEM method. The construction of virtual surface does not affect the neutron transport process. The correctness of the improved DEM is verified by comparing the total filling particle number and calculation results of keff with RSAmethod. At the same time, according to the distribution of filling particles, the improved DEM metho d fills the particles more uniformly.


INTRODUCTION
After the Fukushima Daiichi nuclear disaster, the need for secure nuclear reactors seems to be more important than ever. One of the possible candidates of these fourth generation reactors is the so-called pebble bed reactor (PBR) [1]. However, the double heterogeneity of PBR with stochastic mixture poses challenges for traditional deterministic transport methods. Therefore, Taking the advantages of the flexible geometry modeling and the use of continuous-energy nuclear cross sections, Monte Carlo (MC) codes play an important role in the accurate simulation of stochastic medium [2]. Several MC codes are capable of dealing with random geometries, such as MCNP, Serpent and TRIPOLI-4 [3][4][5].
To take into account the stochastic effect, several methods have been developed, such as implicit method and explicit method. Implicit method does not build stochastic geometry before neutron tracking process. The location of a spherical fuel is sampled probabilistically along the neutron flight path from the spatial probability distribution of fuel particles. Implicit method does not have packing fraction limitation and the time of stochastic medium geometry building process can be saved compared with explicit methods. But, implicit methods cannot obtain right packing fraction. This phenomenon has been observed in stochastic medium packing fraction tests using Reactor Monte Carlo code (RMC). Meanwhile, Serpent developer also suggests not using implicit modeling method when modeling HTR geometries [3].
As for the explicit modeling methods, two main accurate and mature approaches have been applied in PBR simulation, random sequential addition method (RSA) and distinct element method (DEM). RSA method is firstly developed for chemistry field to solved thermodynamic problems [6]. RSA method samples the location of particles one by one with no overlap, therefore, the filling process does not need any dynamic calculation and thus, much quicker than DEM method. RSA method is used for generating TRISO particles distributed in fuel balls with relatively low packing fraction (generally around 30%). DEM is a kind of dynamic approach to create the sphere packages with particles interacting with each other [7]. Li and Wei applied DEM method to pebble bed reactor simulations with CFD-DEM coupled model. The pebble-coolant interaction effects on pebble motions and coolant dynamics are investigated [8]. Rintala et al. also use DEM method to produce randomly packed pebble beds and DEM pebble packing simulation technique will be used to provide detailed data also for the thermal-hydraulic model [3]. However, Renzo also mentioned the drawback of DEM method. In the field of granular motion simulations, it is impossible to incorporate long-range inter-particle forces in the module [9] [10]. Therefore, an improved DEM method is proposed to solve the problem of modeling non-contact particles dispersed in the stochastic media with high packing fraction.
The remainder of this paper is organized as follows. Section 2 describes the two main explicit methods and improved DEM approach. This section also presents calculation tests and results. General conclusions are discussed in the final section.

Explicit modeling methodology and limitation
The theory and application limitation of random sequential addition method and distinct element method will be clarified in this part.

The random sequential addition method(RSA)
The study of the structure and macroscopic properties of hard spheres in physical dimensions (d= 2 or 3) has a rich history, dating back to at least the work of Boltzmann [11]. Widom developed random sequential addition method for a system consisting of many hard spheres and compared the configurations with thermodynamic equilibrium one [6]. The procedure is simple, choose at random a point in the volume, and place a sphere there with its center at the chosen point. From the diminished volume now accessible to the center of next sphere, choose another center randomly. Torquato analyzed the packing fraction of RSA of hard spheres in high euclidean dimensions [12]. For d-dimensional hard spheres, we call φs ≡ φ(∞) the saturation (infinite-time) limit of the density. Feder postulated that the infinite-time limit for d-dimensional hard spheres follows the algebraic behavior: [12] Where τ represents a dimensionless time. Base on the previous investigators work, Zhang has fitted their data of saturation densities as a function of d for 2≤d≤8 using a variety of different functions and found best fit form: Where a1 = 1.0801, a2 = 0.32565, and a3 = 0.11056 are parameters. Therefore, the approximate threedimensional saturation density value is 0.3815.
Meanwhile, we also calculated the packing time with different packing fractions. The fitting curve is shown in Fig. 1.

Figure 1. Fitting of the packing fraction by Eq. (1)
And the saturation density of three-dimension configuration is obtained through formula (1). The saturation density is 0.356, a little bit lower than the value calculated by Eq. (2).

The distinct element method (DEM)
The distinct element method was firstly proposed by Cundall and Strack to characterize the behavior of particles in granular flow simulations [14]. The method is based on the use of an explicit numerical scheme in which the interaction of the particles is monitored contact by contact and the motion of the particles modeled particle by particle. The calculations performed in the distinct element method alternate between the application of Newton's second law to the discs and a force-displacement law at the contacts [14].
As mentioned above, the calculation cycle including two main steps: 1) Newton's second law gives the displacement of particles resulting from the total contact forces; 2) The force-displacement law is applied to obtain the new contact forces from motion received in step (1).
To accelerate the calculation, virtual meshes technique also applied to the DEM method. Therefore, when calculating the contact forces, the global searching will be simplified to local mesh searching. Normal overlap and tangential motion is considered when calculate the contact forces. Therefore, based on the linear spring-dashpot model, the contact forces are written as: Where, Fnij is normal force and Ftij is tangential force. The deformation α is the overlapping distance between the particles. The tangential overlap displacement is δ. The spring is characterized by the stiffness in normal and tangential directions (kn, kt), while damping factors (cn, ct) determine the property of dashpot. The direction of forces is related to the slip velocity vsij and relative velocity vrij.
The Hertz-Mindlin model is used for the stiffness and damping factors calculation: Then, the net contact force Fi and contact torque Ti on particle i could be obtained by: Using the net forces and torques, the translational and rotational speeds are determined by The vector of gravitational acceleration is g, the moment of inertia for particle i is Ii. Meanwhile, it is also important to determine the critical time step given by: [14] 0.163 0.877 Where, γ, G represent Poisson's ratio and shear modulus respectively.

The improved DEM method
According to the analysis above, an improved DEM method should be developed for the non-contact particles dispersed in the stochastic media with high packing fraction. First, we need to let the particles contact with each other by adding virtual surfaces outside the outer layer of particles. These surfaces will disappear when the neutron tracking process begin so that the virtual surfaces won't affect real geometry conditions. With the expanded particles, the packing fraction will be adjusted to a larger value and DEM method could be used to construct the particles distribution module. The radius of virtual sphere surfaces can be written as: where, PFreal is actual packing fraction of TRISO particles, and PFsystem is a fixed packing fraction value to make sure the expanded particles contact with each other.

Test case and results analysis
The typical TRISO fuel particle is used, of which the dimensions are listed in table 1. Single rod filled with TRISO particles is built for method validation and comparison. The reflective boundaries are used for the outer boundaries. The TRISO particles are distributed in the rod and don't contact with each other. With the use of improved DEM method, the 45% packing fraction configuration is obtained and is shown in Fig. 2. The particle distribution is shown in Fig. 2, the gap between each two particles is small and difficult to see. This is due to the system packing fraction value setting to 50%. Therefore, the radius ratio of expanded ball and original one is 1.0357, very close to 1.
To further study the effect of particle distribution to the keff, the packing fraction is set to a lower value, 30%. The calculation results are shown in table 2. The keff is compared with two different explicit method. Even the packing configuration is different, the keff value is almost the same. Therefore, the calculation results imply that both two methods are suitable for the low packing fraction circumstance.

CONCLUSIONS
Random sequential addition method is a kind of efficient filling method compared to the DEM method. This approach is often used for generating TRISO particles distributed in fuel balls with relatively low packing fraction. However, RSA method has up limitation of packing fraction, and this upper limit value has been studied in this work. Meanwhile, DEM method is based on the contact force analysis. Therefore, due to the application limitation of RSA method and DEM method, a new kind of explicit modeling method is proposed in this work based on the DEM approach. The configuration with 45% packing fraction has been successfully obtained. The correctness of the improved DEM is verified by comparing the total filling particle number and calculation results of keff with RSA method. And the similar results implies that particle distributions generated by those two methods do not affect the keff value. This might because both two methods can generate particles randomly, even if the configuration of improved DEM method might be more uniformly.