A Fast Algorithm for the Discrete Element Method by Contact Force Prediction

The discrete element method (DEM) takes enormous calculation time because it requires a very small time step, one small enough to represent the large frequency in the contact dynamic model. In general, the equations of motion of particles are solved using the second-order Adams-Bashforth method, which estimates the values of contact force in the following calculation time by linear extrapolation, or by multi-step methods such as the predictor-corrector method. Inspired by these two conventional methods, we propose a Contact Force Prediction Method that makes a larger time step possible. Our method uses the predicted values of contact force at every contact point, which are exact solutions or numerical solutions of dif ferential equations that represent two particle contacts. It has been confirmed experimentally that the proposed method gives reasonable results of packing and discharge simulations, and accelerates DEM calculation 3 8 times.


Introduction
Discrete Element Method (DEM) simulations on a number of granular systems have reported positive results since this method can estimate many effects at particle level 1) .However, it requires enormous computation time because it needs a sufficiently small time step to follow high fluctuations in the contact dynamic model between particles.To overcome the deficit of conventional DEM, speeding up DEM is proposed in this study.There are several approaches to speeding up DEM, such as promoting efficiency in the algorithm to detect contact between particles, increase in the time step and limitation of the number of particles to calculate.In this study, we report a method to increase the time step.
In DEM calculation, the time step must be set smaller than the particle size and the density becomes smaller or the spring constant becomes larger, in order to obtain a stable numerical solution.Several methods to increase the time step have been reported, such as methods to set the spring constant smaller 2)-4) and the particle size larger 5) than the real value.This method is an effective method of increasing the time step, but it must be used after careful examination of the influence on the phenomenon to which DEM calculation is applied, because delays in contact detection cannot be avoided.
The purpose of this study is to develop an original DEM in which a large time step is possible by improving the algorithm needed to obtain the numerical solution to the contact dynamic equation between particles.We propose a Contact Force Prediction Method in which predicted values of contact force at every contact point are calculated by a more stable method than the conventional DEM and used to solve particle velocity and location.The packing (Fig. 1) and discharge (Fig. 2) system is chosen here as an example of a phenomenon in which contact between many particles is important, and it is shown that the proposed method is effective for calculating the packing ratio or the discharge rate.

Delays in contact search and stability of calculation
When the time step is increased, DEM calculation may become incapable of reproducing a real phenomenon for the following reasons: one is delays in the contact detection and another is the error of the differential method.In the following, the effects of these factors with the conventional method are examined.

Outline of DEM
Contact force f cij in DEM is calculated by projecting on the normal direction (f nij ) and the tangential direction (f sij ) as shown below: where (t) means a function of time t.The normal direction force is calculated from the following equations: where k is the spring constant, h is the damping constant, r is the radius of the particle, d is the overlap between particles, n is the normal unit vector, v is the velocity vector, and x is the location vector.Subscript i denotes the physical values of particle i; j denotes particle j; ij denotes between particles i and j, n denotes the normal direction value, and r denotes the relative value.a 1 and a 2 are constants that depend on the type of contact dynamic model: a 1 ҃1 and a 2 ҃0 for the linear spring model and a 1 ҃1.5 and a 2 ҃0.25 for the Hertz model 6) .The tangential direction force is calculated from the following equations: v sij (t)҃v ij (t)Ҁv nij (t)ѿr i i (t)ѿr j j (t)҂n ij (t), ( 9) where s is the tangential unit vector and is the angular vector velocity.a 3 and a 4 are constants that depend on the type of the contact dynamic model: a 3 ҃a 4 ҃0 for the linear spring model and a 3 ҃0.5 and a 4 ҃0.25 for the Mindlin model 7), 8) .The total tangential force is limited by the Coulomb frictional limit.
where m is the friction constant.
In the linear spring model, it is difficult to theoretically decide the spring constant value except when the particle shape is a disk, but it is usually used because of its convenience in the calculation 9) .

An error from delays in contact detection
Contact detection between particles may produce a maximum delay that equals the time step.This delay can be ignored when the time step is small, but the following two problems occur when the time step is large.
One is the large contact force between particles due to the large strain, which is not practical, because two particles that are in contact with each other continue to come close during a delay in contact detection without receiving contact force from each other.A method whereby the calculation time is returned to the starting point when two particles make contact with each other and where particle locations are recalculated, which overcomes this problem, has already been reported 5), 10) .In this study, this problem is solved using the contact force prediction method described below.
Another problem is that several particles make contact with the same particle at same time.Two contacts with a small time lag may be misjudged as occurring at almost the same time because of the large time step.This phenomenon changes the behavior of particles after they make contact because the consumption of the collision energy is different.
As an example, we consider a difference between a real phenomenon and calculation result using the three particles below.As shown in Fig. 3, there are three particles, A, B and C, without rotation in a onedimensional line, and the particles of both ends, A and C, collide at the same velocity with the central particle, B, whose initial velocity is zero.Fig. 4 shows the calculated velocity of each particle after contact when the beginning of the contact between particle C and B is after ∆t * from the beginning of contact between particle A and B. Although natural gas hydrate pellets 11), 12) have been chosen in this simulation, parameter values used in the simulation can be substituted for the physical properties of ice pellets 13) , as shown in Table 1.The particle diameter is 10 mm, and a linear model is used.The contact time between particles in this simulation is 3.7҂10 Ҁ5 seconds, and the time step is set up in 1/100 of the contact time, in which the calculation error becomes sufficiently small.∆t * in Fig. 4  value of ∆t, the difference between the beginning of the contact between particles A and B and the beginning of the contact between particles C and B, divided by t c , the contact time between two particles: The velocity is also normalized by the velocity of particle A and C before they make contact.Fig. 5 is the simulation result of the normalized velocity error of particle A calculated from Fig. 4, when the contact detection between particle A and B is shifted by the time step in the system in which the contact between A and B and the contact between C and B take place at the same time.It can be concluded that DEM with a large time step is the method that allows the error shown in Fig. 5 compared with DEM with a small time step so that the effect of the contact between the other particles can be ignored.Therefore, the appropriate numerical solution for DEM with a large time step is not a highly precise method, but its calculation is sufficiently fast and stable for a large time step.

Stability of calculation
Fig. 6, except its shadowed area, is the general calculation procedure for DEM.The explicit method is used in DEM and particle i's velocity v i (tѿ∆ t), angular velocity i (tѿ∆ t), location x i (tѿ∆ t) and angular location i (tѿ∆ t) after time step ∆ t are calculated using contact force f i (t) and torque T i (t).In the following, several solutions with the conventional DEM are compared using the linear spring model in the ∆t t c system in which a particle without rotation makes contact with a wall in one dimension.
Although several numerical methods of solving the differential equation that represents the contact between two particles can be considered for the conventional DEM, the easiest with first-order accuracy is the explicit Euler method 14) .On the other hand, the second-order Adams-Bashforth 15) , Leapfrog 16) and predictor-corrector methods 10) are used with the conventional DEM.Although the predicted value of the contact force in the second-order Adams-Bashforth method is calculated by linear extrapolation, there is the problem of second-order accuracy not being strictly guaranteed unless appropriate processing is performed at the discontinuous points, such as the beginning and end points of the contact.Neither is second-order accuracy guaranteed with the Leapfrog method, unless the value of the velocity during the contact is recalculated using an appropriate method, because the particle velocity is defined by only the middle point of a time step.Moreover, although the correction process is usually repeated until an error decreases enough in the predictor-corrector method, it is common to limit the correction process to 1 or 2 times to shorten the computational time in conventional DEM.It is necessary to store all contact states until the calculation of predicted values for the velocity and location are finished because predicted values are calculated after all contact forces are calculated with this method.This processing increases computational time.Table 2 compares the 4 above-mentioned methods.One-dimensional collision of one particle is assumed, and subscript i and n are omitted.In addition, subscript k is defined with k҃t/∆ t and ∆ t; time step, is simplified with h.The contact force is calculated assuming a linear spring model as follows: The analysis contents are explained as follows.In the case of the explicit Euler method, the explicit Euler method is expressed in the matrix description as follows: The determinant of A is as follows: where m is the mass of the particle.From Equation (15), it is found that the eigenvalue of matrix A may be larger than 1.This shows the possibility of deviating from a stability domain of a calculation with the explicit Euler method.The calculation stability in Table 2 is calculated in the same way as when h҃0. ∆t * in Table 2 is defined as follows: where t c is the contact time when h҃0.
The conservation of energy in the explicit Euler method is examined when h҃0 as follows: Equation (18) shows that energy is not strictly conserved, but there is a conservative quality for a conversion when the eigenvalue of matrix A is 1.More concretely, this means that the calculated energy value shows the minute fluctuations up and down around the true energy value and does not greatly decrease or increase.The energy conservation in Table 2 is examined in the same way as when h҃0. "Good" shows the case where the energy is not conserved but the eigenvalue is 1.
Numerical equations in terms of the energy conservation in Table 2 show the error of energy conservation, which is defined as follows: n k•kҀ1 and d k•kҀ1 are defined as follows: The calculation time in Table 2 shows the time needed to calculate one-dimensional collision of a particle when the calculation time in the explicit Euler method is 1.Adscript ѿ in the second-order Adams-Bashforth method shows the calculation time increases when appropriate processing is performed at discontinuous points, such as the beginning and end points of the contact to guarantee second-order accuracy.With the Leapfrog method, the following formulae for the beginning and end points are used to guarantee second-order accuracy: However, Equation ( 22) cannot be used for the explicit method when h 0, so more appropriate processing such as the multi-step method is needed to guarantee second-order accuracy.Adscript ѿ with the Leapfrog method shows the calculation time increases for the reasons given above.However, it is necessary to note that the contact detection holds most of the calculation time rather than the calculation of contact force in a state of crowded particles.
Although the Leapfrog method has the above-mentioned problem, Table 2 shows that it is the superior of the three methods with second-order accuracy in calculation stability and energy conservation.
With the above-mentioned consideration, a onedimensional collision of one particle with a wall is assumed.It is necessary for the time step to set up a smaller value in the case of collision of many particles, because the frequency of the contact dynamic model becomes larger when the number of contact points increases.In the following section, we propose a method whereby a stable numerical solution can be found for a large time step in a system of multiple particle collisions.

Outline of contact force prediction method
A predicted value of the contact force at the stage where two particles collide is solved by this method.
where " ˆ" denotes a predicted value.We refer to this method in which a predicted value of the contact force is used in this way as the "contact force prediction method" in this study.Fig. 6 shows the calculation procedure for this method.As compared to the conventional calculation procedure, the procedure for predicting the contact force at every contact point, the shadowed area in Fig. 6, is added.Thus, it is necessary for the calculation time needed to obtain the predicted value to be small.The meaning of using a predicted value of a twoparticle collision for a multi-particle collision is examined here from the viewpoint of numerical analysis.As an example of a multi-particle collision, we assumed a one-dimensional collision of many particles without rotation.The dynamic equations of particle i and j in a one-dimensional collision system are as follows: where the second term of the right-hand side in each equation is a summation of forces acting on particles i and j except for the contact force between particles i and j.For example, achieving second-order accuracy for the velocity of particle i depends on the following dn j (t) dt dn i (t) dt equation: where the terms above third-order accuracy for ∆ t are omitted.Equation ( 24) is substituted for Equation (26) as follows: On the other hand, the following equation relates to Equation (28).
From Equation (28) and Equation (30), it is found that second-order accuracy is not achieved unless the inf luence of the other contact forces, ∑ f i (t) and ∑ f j (t), is considered in the predicted value.However, it is not necessary to achieve second-order accuracy with this method because DEM with a large time step dn i (t) dt permits an error of a certain order for the calculation of the contact force, as mentioned in Section 2. Equation (27) shows that first-order accuracy is achieved even if the inf luence of the other contact forces, ∑ f i (t) and ∑ f j (t), is ignored.In the following, a method whereby the inf luence of other contact forces is ignored for the predicted value is considered.
Although various methods can be used to obtain the predicted value, four methods whereby stable solutions can be achieved in a large time step are selected and shown in Table 3.The manner of consideration is similar to that used in Table 2, and the stability of calculation and energy conservation are examined when h҃0. "Exactly Good" in terms of energy conser vation means that the particle energy is conserved strictly according to the examination similar to Equation ( 18).The modified Runge-Kutta method with second-order accuracy is as follows: This differs from the conventional Runge-Kutta method with second-order accuracy in that Equation 12 No Good  (34) is the trapezoid method.Although the calculation when h҃0 is unstable with the conventional Runge-Kutta method with second-order accuracy, it is interesting that this method is stable even with such little modification."Exact Solution" in Table 3 means a method whereby an exact solution to the differential equation of the contact dynamic model is used for the predicted value."Stable" in terms of the stability of calculation in Table 3 means that a stable solution for the predicted value of the contact force can be achieved with a large time step.However, there is another limit for the time step in DEM calculation, as mentioned in Section 2.2.The time step has to be set to less than the contact time t c , i.e., ∆t * ͨ1.
From Table 3, the implicit trapezoid method and the exact solution method are superior to the other two methods in terms of the large time step and energy conservation, but these methods cannot be applied to the non-linear spring model.The method that can be applied to the non-linear spring model and permits a large time step is the fourth-order Runge-Kutta method.The calculation time is smallest with the implicit trapezoid method and the exact solution method.However, the effect of the time step that can be set larger is bigger for the total speedup of DEM calculation than the calculation time of the predicted value is small.This is because the calculation time of the contact detection is larger than that of the calculation of the contact force, as mentioned above.
A predicted value is calculated by projecting on the normal direction and the tangential direction similar to the contact force.dnij (tѿ∆ t) and dsij (tѿ∆ t) define the predicted value of the particle location in the normal and tangential direction, respectively; nnij (tѿ∆ t) and nsij (tѿ∆ t) define the predicted particle velocity in the normal and tangential direction, respectively.These values can be obtained by either method, as mentioned above, and the predicted value of the contact force in the normal direction, fnij (tѿ∆ t), is obtained as follows: where a 1 and a 2 are similar to the constants in Equation (2).With the conventional DEM, tension in the normal direction is not permitted by the no-tension joint, so the predicted value of the contact force in the normal direction is also not permitted as follows: The predicted value in the tangential direction, fsij (tѿ∆ t), is obtained as follows: where a 3 and a 4 are similar to the constants in Equation ( 6).∆ tcij is the time when the two-particle collision finishes, i.e., f * nij ҃0, and is obtained by the following linear approximation equation, for example: Moreover, the predicted value of the contact force in the tangential direction is limited by the Coulomb frictional limit: (40)

Application examples of the contact force
prediction method Some examples of DEM calculation using the abovementioned original method are introduced here.In this study, the packing and discharge system are used as an example of multi-particle collision.
The properties of the particles have already been shown in Table 1.The particle size is 50 mm, the number of particles is 200, and the calculation space is a quadratic column whose bottom is 0.3 m҂0.3 m.The contact time between two particles in this calculation is 2.0҂10 Ҁ4 seconds.The periodic boundary condition is used in all side walls in order to moderate the wall effect.The linear spring model is used for the contact force model, and the spring constant is determined from the contact time by the Hertz spring model.
In all calculations, the exact solution method is used to obtain the predicted value of the contact force with the contact force prediction method.

Packing simulation
The packing simulation results produced by the free fall of particles are shown here.Particles fall freely from a height of 1 m.Fig. 7 shows the time change of the energy summation of all particles in the case of the contact force prediction method or the conventional method, respectively.n t is the number of partitions in the contact time as follows: where t t is the contact time needed for one particle to contact a wall.Although the calculation is dispersed in n t 3 with the conventional method, it is not dispersed in n t ͨ3 with the proposed method.Fig. 8 shows the comparison between the conventional method and the proposed method for the calculation results of the packing ratio.The calculated packing ratio becomes larger as the time step becomes t c ∆ t larger.It is assumed that this is because the overlap at the beginning of contact between two particles is maintained.How fine the time step must be set with the proposed method depends on how much calculation error is permitted.In Fig. 8, although the packing ratio does not depend on the time step in n t 300 with the conventional method, n t 30 with the proposed method.In other words, the proposed method can set the time step 10 times larger than that of the conventional method.The calculation time of 1 step with the proposed method is 1.2 times larger than with the conventional method.Thus, the proposed method accelerates DEM calculation 8 times when the time step is set 10 times larger than conventionally.
Incidentally, there have been various considerations regarding the random packing of mono-dispersed particles in the long term 17), 18) .The random packing ratio should be between the closest packing ratio, 0.74, and the packing ratio in a simple cube lattice, theoretically 0.52.Moreover, experimental results of mono-dispersed spherical particles by Westman and White 17) show that their packing ratio is between 0.553 and 0.63.Therefore, the packing ratio 0.57 obtained in this study is the appropriate value.

Discharge simulation
In order to simulate a system in which particle movement is more dynamic than in the packing system, a 0.2-m diameter circular outlet is opened at the bottom of the packed particles bed, and the particles are allowed to fall freely from the outlet to the f loor, a distance of under 1 m.
Similar to Fig. 7, Fig. 9 shows the time change of the energy summation of all particles.Although the calculation is dispersed in n t ҃3 with the conventional method, it is not dispersed in n t ͨ3 with the proposed method.Fig. 10 shows a comparison between the conventional method and the proposed method for the calculation results of the number of discharged particles from the outlet.As for the discharge behavior, a big difference between the conventional method and the proposed method is not found, even with the large time step.It is believed that this is because the calculation error is canceled out due to the movement of all particles.In Fig. 9, the kinetic energy is not dispersed in n t ҃6 with the conventional method and n t ҃1.5 with the proposed method, so the proposed method can set the time step 4 times larger than that of the conventional method.When the calculation time of 1 step is considered, the proposed method accelerates DEM calculation 3 times.

Conclusion
DEM in which a large time step is possible is developed by improving the algorithm to solve the contact force between two particles.In general, the equations of motion of particles are solved by the second-order Adams-Bashforth method, which estimates the values of the contact force in the following calculation time by linear extrapolation, or by multi-step methods such as the predictor-corrector method.Inspired by these two conventional methods, we propose an original method called the contact force prediction method, which uses the predicted value of the contact force at each contact point of two particles.As for the method of obtaining the predicted value of the contact force with the proposed method, four methods in which a large time step can be set are examined, and their characteristics are shown.Moreover, the packing and discharge system is chosen as an example of the multi-particle collision, and the proposed method accelerates DEM calculation 3Ҁ8 times.
In this paper, some calculation results of mono-dispersed particles are introduced as an imminent example, but we have already confirmed that the contact force prediction method is useful for two sizes of particle at a number of approximately 10 million.Under the present parameters, it is expected that the number of particles calculated in DEM will be limited to 100 million, even if several techniques to accelerate DEM calculation, such as the method proposed here, are applied.Thus, an additional acceleration, such as parallel computation or combination with a continuous method, is needed for application to the general scale system.However, it is confirmed that the contact force prediction method proposed in this study is useful for DEM acceleration.

Fig. 5 Fig. 6
Fig. 5 Error of velocity of particle A

Fig. 7 Fig. 8
Fig. 7 Kinetic energy of total particles in packing simulation

Table 3
Comparison of 4 methods of calculating the estimated value of force (k҃t/∆t, h҃∆t)