Assessment of the potentials of implicit integration method in discrete element modelling of granular matter

Computers and Chemical Engineering 49 (2013) 183–193

Contents lists available at SciVerse ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Assessment of the potentials of implicit integration method in discrete element modelling of granular matter K. Samiei a , B. Peters a,∗ , M. Bolten b , A. Frommer b a b

Research Group in Engineering Science, Faculty of Science, Technology and Communication, University of Luxembourg, L-1359 Luxembourg City, Luxembourg Department of Mathematics, Faculty C-Mathematics and Natural Sciences, University of Wuppertal, D-42097 Wuppertal, Germany

a r t i c l e

i n f o

Article history: Received 22 May 2012 Received in revised form 25 August 2012 Accepted 20 October 2012 Available online 3 November 2012 Keywords: Implicit integration DEM Numerical analysis Granular material

a b s t r a c t Discrete element method (DEM) is increasingly used to simulate the motion of granular matter in engineering devices. DEM relies on numerical integration to compute the positions and velocities of particles in the next time step. Typically, explicit integration methods are utilized in DEM. This paper presents a systematic assessment of the potentials of implicit integration in DEM. The results show that though the implicit integration enables larger time steps to be used compared to the common explicit methods, the overall speed up is overruled by higher computational costs of the implicit method. © 2012 Elsevier Ltd. All rights reserved.

1. Introduction In order to predict and optimize the behaviour and motion of granular matter in engineering devices, numerical simulation tools are increasingly employed (Cleary, 2004). To date the discrete element method (also called distinct element method) is the leading approach to simulate the dynamics of granular media. The DEM is a numerical approach where statistical measures of the global behaviour of a phenomenon are computed from the individual motion and mutual interactions of a large population of elements (Cundall & Strack, 1979). Modelling is straightforward: the grains are the elements, they interact through local, pairwise contacts, yet are also subject to external factors such as gravitation or contacts with surrounding objects, and they otherwise obey Newtons laws of motion (Kozicki & Donzé, 2009). In contrast to the continuum approach, DEM analysis accounts for inter-particle contacts. However, run-time efficiency is still a limiting factor in large scale applications and therefore research and investigation of alternative methods and algorithms and evaluations of their costs and benefits are of grave importance. The Lagrangian time driven method is applied to the discrete particles of a moving ensemble which is regarded as a system of a finite number of visco-elastic particles with a given shape and material properties. The state of particles is obtained by time

∗ Corresponding author. Tel.: +352 466644 5496; fax: +352 466644 5200. E-mail address: [email protected] (B. Peters). URL: http://www.xdem.de (B. Peters). 0098-1354/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compchemeng.2012.10.009

integration of the dynamics equations derived from the classical Newtonian mechanics approach based on the Newton’s second law for translation and rotation of each particle in the particle ensemble. All the forces and moments acting on each particle are evaluated at every time step. The state of particles at the next time step can be calculated from the state of particles at the current and/or previous time steps. This will simplify solving the equation because the new positions could be expressed as explicit functions of the already known values. If the state of particles at the next time step is calculated not only from the current and previous time steps but also from the next time step, the equation of motion will be implicit in new positions. Generally the implicit method is computationally more expensive because it requires a system of equations to be solved at each time step. On the other hand, relatively larger time steps could be used in implicit methods due to higher numerical stability. The common practice in DEM simulations is the explicit updating and the use of implicit methods has been very limited. A method called discontinuous deformation analysis (DDA) (Ke & Bray, 1995), was one of the first works presenting an implicit method for two dimensional simulation of particulate media. This method is claimed to solve systems of few thousands of particles in “reasonable times” though no comparison with explicit methods was presented. More recently, Schäfer and Negrut (2010) evaluated the potential of implicit integration methods in molecular dynamics simulation of biological molecules. Although they report good energy conservation response by the implicit methods, the increase in the time step was limited due to loss of convergence of the iterative method. Tuley, Danby, Shrimpton, and

184

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

Palmer (2010) investigated the performance of different integration schemes including implicit schemes in solving the motion of a particle during a normal collision with a wall. This work was later advanced by Jasion, Shrimpton, Danby, and Takeda (2011) including the response of the same integration schemes to a collision in the tangential direction of motion. Their results show that more simple integrators such as a first order symplectic Euler scheme are more suitable than costly and conditionally more accurate integrators for simulations containing thousands of particles. Overall, the experience with implicit methods in the context of DEM has been limited and reporting mixed, i.e. both positive and negative, results. Despite several investigatory studies discussing the advantages and disadvantages of the implicit methods, very few have drawn their conclusions out of the results of an implemented implicit method. This could be, among other reasons, due to the fact that implicit methods are in general harder to implement than the explicit integrators. Therefore, the objective of this research is to implement an implicit integration algorithm for numerical simulation of granular matter and to asses its costs and benefits compared to the conventional explicit methods.

Fig. 1. Forward, backward and central difference approximation applied to the velocity curve during the loading phase of a linear spring.

2. Integration methods

2.2. Taylor

The objective of integrating the Newton’s equations of motion is to derive the velocity and position of the particles at the new time step based on the forces acting on them. The contact force in DEM is in general a function of the positions and velocities of particles in contact. On the other hand, it determines the particle’s acceleration according to the Newton’s second law of motion. Therefore the differential equation describing the motion of a particle involves the first and second derivatives of the position as expressed in Eq. (1).

In this method, the Taylor polynomial of degree one is used to calculate the velocity:

2

∂ ri 1  = F mi i ∂t 2

 ri ,

∂ri ∂t

 (1)

Various integration methods exist to solve such an Ordinary Differential Equation (Ferziger & Peric, 2002; Negrut, 1998). The integration formulas corresponding to the schemes used in this paper are briefly described in this section. All the integration schemes discussed in this paper can be derived by approximating the derivatives of a function by Taylor series: f (t + t) = f (t) +

f  (t) f  (t) f (3) (t) t + t 2 + t 3 + · · · 1! 2! 3!

(2)

2.1. Symplectic Euler

(3)

The symplectic Euler method is different from the standard Euler in calculation of the position because it uses the backward difference approximation. The new position is then can be written as: rt+1 = rt + vt+1 t

(5)

The position is derived from the Taylor polynomial of degree two: rt+1 = rt + vt t +

1  t t 2 a 2

(6)

This method can be derived from constant acceleration motion, too. In other words, this integration method assumes the acceleration during the time step to be equal to the acceleration at the beginning of that time step. This is an explicit integration scheme. 2.3. Position Verlet The position Verlet method uses the central difference approximation for the second derivative of a function to compute the position:  t t 2 rt+1 = 2rt − rt−1 + a

(7)

The velocity can be derived from central difference approximation of the first derivative:

The standard Euler method is based on the forward difference approximation for the first derivative of a function. So the velocity at the new time step can be written as follows:

vt+1 = vt + a t t

vt+1 = vt + a t t

(4)

This is a modification to standard Euler formula where only forward difference approximation is used for both velocity and position (Hairer, Lubich, & Wanner, 2003). The symplectic Euler method is also called the semi-implicit Euler method because it utilizes the new velocity in prediction of the new position though velocity itself is solved explicitly in terms of old values.

vt =

rt+1 − rt−1 2t

(8)

However, the velocity in Eq. (8) is one step behind the position. The acceleration is assumed constant during a time step in this method but the value at the middle of the time step is used to approximate the constant acceleration rather than the beginning or end of the time step in forward or backward approximations. Fig. 1 illustrates the forward, backward and central difference approximations of the new velocity vt+t during loading of a linear spring. As shown in the figure, the central difference approximation is generally more accurate than the forward and backward approximations. The backward difference approximation is more accurate than the forward difference in the loading phase as shown in Fig. 1. On the contrary, the forward difference is a better approximation than the backward difference in the unloading phase when the acceleration is declining. The position verlet method is an explicit second order integrator.

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

185

2.4. Velocity Verlet As noted above, the basic position Verlet method does not provide velocities at the next time step. To overcome this problem, a modification to the basic Verlet method is widely used (Hairer et al., 2003). The position in the velocity Verlet method is computed based on the central difference approximation similarly to the basic Verlet method. Replacing rt−1 from Eq. (8) into Eq. (7), the position is derived as: rt+1 = rt + vt t +

1  t t 2 a 2

(9)

The velocity is approximated as follows: 1 2

vt+1 = vt + (a t + a t+1 )t

(10)

To implement this, a predictor-corrector scheme is developed. The predictor, calculates the position from previous values according to Eq. (9) and predicts the velocity as: 1 2

 t t vpr = vt + a t+1

(11) Fig. 2. Geometry of contact between spherical particles.

At the next step, the force calculation function is called and accordingly the new acceleration is calculated. The corrector, then, corrects the new velocity according to Eq. (10). This implementation does not require more memory storage than the basic Verlet method but the calculation of the new acceleration is based on the predicted velocity value instead of the corrected value. This method can be considered semi-implicit in calculation of the velocity while it is explicit in calculation of position.

therefore, categorised as a semi-implicit method. It requires only one force evaluation per time step. 2.6. Numerov The Numerov’s method, also called Cowel’s method or sometimes the implicit Störmer method (Skeel, Zhang, & Schlick, 1997), can be written as follows:

2.5. Gear The Gear algorithm (Gear, 1966, 1971) is known for its numerical stability and is implemented in two steps: predictor and corrector. At the first step the Gear-predictor predicts the new position, velocity, acceleration by extrapolating the current values using the Taylor expansion. The degree of the Taylor polynomial depends on the order of the Gear algorithm. For the forth order Gear scheme, 3 the time derivatives up to ∂ rt /∂t 3 are required: pr rt+1

 t t + vpr = vt + a t+1 t +  pr =a a t+1

(12)

3

1 ∂ rt t 2 2 ∂t 3

(13)

3

∂ rt t ∂t 3

(14)

At the next step, the force calculation function is called and accordingly the new acceleration based on the predicted values for positions and velocities is estimated. Then the Gear-corrector corrects the positions, velocities and accelerations based on the difference between the acceleration estimate and the acceleration =a  est − a  pr . The correction equations for the forth prediction: a order Gear are shown in Eqs. (15)–(18). pr corr = rt+1 + rt+1

1  t 2 a 12

(15)

pr + vcorr t+1 = v t+1

5  ta 12

(16)

 corr  pr  a t+1 = at+1 + a 3

t 2 t + a  t−1 )  t+1 + 10a (a 12

(19)

The Numerov’s method does not provide the velocity by itself. The backward difference approximation can be used to predict the velocity:

vt+1 =

rt+1 − rt t

(20)

The Numerov method is an implicit integrator.

3

1 1 ∂ rt  t t 2 + = rt + vt t + a t 3 2 6 ∂t 3

∂ rt+1 ∂t 3

rt+1 = 2rt − rt−1 +

corr

3

=

∂ rt+1 ∂t 3

(17) pr

+

 a t

(18)

By using predicted values instead of actual values, the Gear scheme converts the implicit equations to explicit ones. It can be,

3. Linear spring dashpot force model The damped linear spring contact force model, based on the work by Cundall and Strack (1979), is chosen in this study. This is a simple model where the analytical solution is quite easily available which makes it suitable for validation of the results. The normal force in this model is composed of a linear spring force and a dashpot force where the spring produces an elastic repulsive force and the dashpot contributes to the damping. The magnitude of the normal force can be written as ˙ Fn = mij ı¨ = −(kn ı + cn ı)

(21)

where ı is the overlap depth between the contacting pair, mij , kn and cn are the reduced mass, the normal spring stiffness and the normal damping (dissipation) coefficients, respectively. The overlap depth can be written as: ıij = Ri + Rj − |ri − rj |

(22)

where ıij stands for the overlap depth between the two spheres and R and r represent the radius and the position vector of the particles, respectively. The normal direction of contact is simply the unit vector connecting the centre of the two spheres as illustrated in Fig. 2. Referring to Eq. (22) which defines the overlap depth between ˙ two contacting spheres, the first derivative of the overlap (ı)

186

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

corresponds to the relative normal velocity in physical sense. Sub¨ corresponds to sequently, the second derivative of the overlap (ı) the relative normal acceleration. The reduced mass can be derived from Newton’s third law of motion affirming the forces of each particle on the other are of equal magnitude and in the opposite direction: mij ı¨ = −mij (¨ri − r¨ j ) = −mi r¨ i = mj r¨ j

Read Input Data

Detect Contacts

(23)

Calculate Forces

where m and r represent the mass and position of the contacting particle. The reduced mass can then easily be derived as: 1 1 1 = + mij mi mj

and Torques Integrate for Posi-

(24)

tion and Orientation

The force differential equation presented in Eq. (21) has an analytical solution which with boundary conditions ı(0) = 0 and ˙ ı(0) = v0 , could be written as follows: ı(t) = e

c t

n − 2m

ij

2v0 · sin · 

 t  2



with  =

4mij kn − cn 2 mij

Increment Time

(25)

No

The relative normal velocity ı˙ can be derived from the analytical solution as follows: ˙ =e ı(t)

−cn t/2mij

·

v0





cn − sin mij

 t  2

+  cos

 t 

initial

Terminate



4mij kn

(28)

2 + ln e2

The total duration of contact can also be derived from the analytical solution, Eq. (25): Tcontact =

mij (2 + ln e2 )

(29)

kn

The total duration of contact is independent of the impact velocity in linear force models and it is therefore expressed entirely in terms of material properties. The maximum overlap of the contact could also be calculated from the analytical solution:





ımax = e

computations for all time steps simulates the motion of the system in the specified duration of time. As already stated, the Newton’s second law is employed to calculate the motion of particles. 4.1. Explicit integration





Fig. 3. The basic flowchart for explicit time-driven method.

(27)

where e denotes the coefficient of restitution. A coefficient of restitution of 1 corresponds to no damping at all and a value of 0 means 100% dissipation of the impact energy. Noting that ı˙ final is when ı(t) in Eq. (25) re-approaches zero, the normal dissipation coefficient can then be derived from Eqs. (26) and (27) as follows: cn = ln e

Yes

(26)

2

Employing the analytical solution, the normal dissipation coefficient can be derived in terms of the coefficient of restitution. Coefficient of restitution is an important characteristic of the impact which determines the effect of energy dissipation. It is defined as the ratio between the normal component of the relative velocity of the contacting pair after and before the impact:

ı˙ final e= ı˙

End Time?

−c/2

4mij kn

−c 2



2mij v0 4mij kn − c 2



(30)



which in case of no damping will reduce to ımax = v0 mij /kn . The frictional forces in the tangential direction of contact are not part of this study. 4. Algorithm The principle behind the time-driven DEM is to discretize time into small time steps where the motion of the system of particles and boundaries is computed for each time step. Summing up the

The basic flowchart of the method is sketched in Fig. 3 and is briefly described below. After an initialization step, the time loop starts. First, the collision partners for each particle need to be detected based on Eq. (22). Once two particles are identified as collision partners, their contact force is calculated. The torque applied on each particle is then calculated based on its contact force and the position of impact. The contact forces and torques are then used to calculate the acceleration and angular acceleration of each particle based on Newton’s second law for translation and rotation. Explicit integration leads to the acceleration’s first and second integrals: velocities and positions. Each time step concludes once the positions, orientations, velocities and angular velocities of all particles in the domain space are updated. The time loop continues until the defined end time of the simulation has reached at which point the program safely terminates. 4.2. Implicit integration Choosing an implicit integration scheme will fundamentally change the DEM algorithm. Before presenting the flowchart for the implicit DEM, the implicit approach needs to be explained. Normal repulsive force between spherical particles is commonly modelled as elastic spring force with linear or non-linear forcedisplacement relations depending on the force model. The general differential equation describing the motion of particle i due to contact with any other rigid body j can then be written as:   ∂ ri kij ˛ (ri − rj ) = (R + Rj − |ri − rj |) · + g mi i |ri − rj | ∂t 2 2

N

j=1

(31)

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

187

where k, m, R and r represent the stiffness, mass, radius and position, respectively, N the total number of bodies including particles and boundary shapes in the system, g the gravity vector and ˛ is the coefficient dependent on the force model. Of course not all other shapes have contacts with particle i necessarily for which the overlap and consequently the force would be zero. With i = 1, 2, . . ., n in Eq. (31), a system of non-linear differential equations is formed with the same number of unknowns as the number of equations. There is a unique answer for such a system which is to be found. Applying the Numerov formula, Eq. (19), to the system of equations given in Eq. (31) would lead to:

Read Input Data

Reset Jacobian and Right Hand Side Detect Contacts Assemble Jacobian and Right Hand Side

ri,t+1 =

t 2  i,t + a  i,t−1 ) + 2ri,t − ri,t−1 , i = 1, 2, 3, . . . , n  + 10a (a 12 i,t+1 (32)

Solve Linear System

where  i,t+1 = a

N ˛

kij (Ri + Rj − |ri,t+1 − rj,t+1 |) (ri,t+1 − rj,t+1 ) j=1

 i,t = a

mi |ri,t+1 − rj,t+1 |

N ˛

kij (Ri + Rj − |ri,t − rj,t |) (ri,t − rj,t ) j=1

 i,t−1 = a

mi |ri,t − rj,t |

+ g

N ˛

kij (Ri + Rj − |ri,t−1 − rj,t−1 |) (ri,t−1 − rj,t−1 ) j=1

mi |ri,t−1 − rj,t−1 |

+ g

(33)

kn (ıt+1 + 10ıt + ıt−1 ) 12

No

Convergence?

Yes

(34)

Integrate for Orientation

+ g

(35)

The non-linear system of equations can be solved by applying the Newton method (Polyak, 2007). The unknowns are the new positions, ri,t+1 while the current and old positions ri,t , ri,t−1 are known values. The Newton method adopts an iterative strategy where a linear equation is solved at every iteration assuming an initial guess for the positions. The equation is solved and the guesses are improved until convergence is achieved. The solution method briefly explained above forms the basic algorithm for the implicit method which is summarized in the flowchart in Fig. 4. Contact detection in the implicit method is performed in the same way as the explicit method but it is carried out at every iteration. This can be explained by recalling that the force in an implicit scheme depends on the next overlap in addition to the current and old overlaps as shown in Eq. (36) for the implicit Numerov formula: Fn = −

Increment Iteration

(36)

Therefore, at every iteration the new solved positions need to be used for calculating new overlaps based on Eq. (22). Identifying contacts at every iteration has its advantages and disadvantages. On one hand contact detection is a costly part of the DEM algorithm and performing this at every iteration adds to the computational burden. On the other hand contacts can be detected even during a time step in the implicit method which is advantageous over explicit method where contacts are detected only at the beginning of the time step. The difference between implicit and explicit methods in contact detection is demonstrated in depth in Section 5.2. To set up the linear system inherent in the Newton method, the right hand side matrix and the Jacobian need to be assembled. The entries of the Jacobian matrix are basically the partial derivatives of the right hand side matrix. The matrices need to be updated at every iteration by new positions. Once the linear system is set up, a linear solver can be employed to solve it. The Gauss–Seidel iterative method is chosen to solve the system (Barrett et al., 1994). Convergence in the Newton method is enforced by a tolerance

Increment Time

No

End Time?

Yes Terminate

Fig. 4. The basic flowchart for implicit method.

value. Convergence is defined strictly in the algorithm in the sense that all particles in the system need to converge before the program can proceed to the next time step. The equation of rotational motion of the particles are not solved implicitly but the torque acting on shapes is calculated as the average of the two torques at the beginning and the end of that time step. In this way, the positions at the next time step are also used to compute the torque acting on particles. This is in contrast to the explicit method where only the force at the beginning of the time step is used to calculate the torque on the particles. The implicit algorithm is implemented in C++programming language and is incorporated into the extended discrete element method (XDEM) software which is an advanced numerical simulation tool for particulate applications. 5. Validation results and discussion The motion of particles can be studied in three different situations: • during the free motion when gravity is the only constant force applied on the particles, • during a collision with another object, • during the interface between the two situations above when the commencement or termination of the contact is first detected.

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

The performance of different integration schemes is not much different during the free motion of particles. The challenging cases are the two last situations. Section 5.1 compares the response of the discussed integration schemes to solve the motion of a particle during collision. Another case is presented in Section 5.2 which explores the difference between the explicit and implicit methods in detecting an upcoming contact. While Section 5.1 focuses on the accuracy of the integration schemes, the cases presented in Section 5.3 focuses on the numerical stability of implicit and explicit methods. As the value of the time step is dependent on the specific studied case, it has no general meaning when analysing different cases. Therefore it would be more meaningful to use the contact resolution defined in Eq. (37) for comparisons. Tcontact t

1.4

1 0.8 0.6 0.4 0.2

(37)

0

Since a linear force model is used in the following verification tests, the contact duration is not dependent on the impact velocity and can be precisely calculated from the size of the particles and their mechanical properties. Whenever the problem includes polydisperse particles, the minimum contact duration is considered. A survey of previous DEM work shows that the size of time step is often 50 times smaller than the duration of contact (Babic, Shen, & Shen, 1990; Campbell, 2002; Ketterhagen, Curtis, & Wassgren, 2005; Thompson & Grest, 1991; Walton & Braun, 1986).

-0.2

CR =

5.1. During collision The objective in this validation test is to compare how accurate different integration methods predict the position of the particle during collision with a wall. The time step chosen here is four times smaller than the total duration of contact. In other words, the collision is resolved in four discrete time steps. The position Verlet method needs the position at two previous time steps in order to predict the next position. Therefore and again for a fair comparison, the initial positions for the first two time steps are provided to all the integration schemes. The graphs presented in this section are based on normalized values in order to represent general trends independent of any specific case. The values are normalized over the maximum point of the analytical solution in the absence of energy dissipation.

Analycal Symplecc Euler Taylor Posion Verlet Velocity Verlet Gear-order4 Implicit Numerov

1.2

Normalized Posion [-]

188

5.1.2. Normal energy dissipation The energy dissipation is modelled as a linear dashpot force dependent on the relative normal velocity during contact. There are different ways to integrate the normal energy dissipation into the Numerov’s method. Before arriving at comparisons between

0.5

0.75

1

Normalized Time [-] Fig. 5. Prediction of the position of a particle during collision with a wall, comparison of different integration methods, contact-resolution = 4, e = 1.

implicit and explicit methods, the methodology to integrate the energy dissipation in the implicit method needs to be determined. The normal velocity, being the first derivative of the position, can be approximated by the backward difference formula. The backward difference or backward Euler method incorporated into the Numerov’s method leads to Eq. (38): Fn = −

kn ıt − ıt−1 (ıt+1 + 10ıt + ıt−1 ) − cn · 12 t

(38)

Using this approach, the dissipation term would not contain the new position and can be formulated explicitly in terms of the known values. Alternatively, the forward difference approximation will introduce the new position in the dissipation term as shown in Eq. (39): Fn = −

kn ıt+1 − ıt (ıt+1 + 10ıt + ıt−1 ) − cn · 12 t

(39)

Finally, the Numerov’s method in conjunction with the central difference approximation of the normal velocity leads to Eq. (40): Fn = −

kn ıt+1 − ıt−1 (ıt+1 + 10ıt + ıt−1 ) − cn · 12 2t

(40)

Fig. 6 shows how these different combinations predict the particle’s position assuming a contact resolution of 4 and a coefficient 1.2

Analycal Numerov-forward Numerov-backward Numerov-central

1

Normalized Posion [-]

5.1.1. Normal repulsive force Initially energy dissipation is excluded and the contact resolution, defined in Eq. (37), is set to 4. Fig. 5 shows the results of the prediction of the particle’s position. The results show that the Taylor method’s prediction is less accurate than the other methods. This is expected as the Taylor method is a first order integrator. The symplectic Euler method is also a first order integrator but it is semi-implicit in calculation of the new position and is therefore quite accurate. In the absence of velocity-dependent dissipative term, the position Verlet produces the same prediction for the particle’s position as the velocity Verlet does, as expected. Considering both velocity and position, the velocity Verlet and the Gear integration schemes provide quite accurate predictions in the absence of energy dissipation. Best predictions are provided by Gear and implicit Numerov schemes. Gear’s method predicts the maximum overlap better than Numerov’s method while on the other hand, the position at the end of the collision is better predicted by Numerov’s method.

0.25

0

0.8

0.6

0.4

0.2

0

0

0.25

0.5

0.75

1

Normalized Time [-] Fig. 6. Prediction of the position of a particle during collision with a wall, comparison of different methods to integrate energy dissipation, contact-resolution = 4, e = 0.5.

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193 1.6

Analycal Symplecc Euler Taylor Posion Verlet Velocity Verlet Gear-order4 Implicit Numerov

1.4 1.2

Normalized Posion [-]

189

1

Fig. 8. Detection of contacts in the explicit solution.

0.8 0.6 0.4 0.2

Fig. 9. Detection of contacts in the implicit solution.

0 - 0.2

0

0.25

0.5

0.75

1

Normalized Time [-] Fig. 7. Prediction of the position of a particle during collision with a wall, comparison of different integration methods, contact-resolution = 4, e = 0.5.

of restitution of 0.5. It is clear that the backward difference approximation of the velocity does not provide an accurate prediction comparing to the two other methods. The forward difference approximation is a better prediction during the loading phase while the central difference approximation is a better match with the analytical solution during the unloading phase. Fig. 7 shows the comparison of the implicit Numerov method combined with the central difference approximation with the explicit methods. The results show that by adding dissipation, the integration methods deviate from the analytical solution. This is because the dissipation term in Eq. (21) increases the loading force while the change to the unloading force becomes more sharp. It is observed that in all the explicit integration schemes, the first prediction of the position after impact and consequently the maximum overlap are over-estimated. This is because an explicit integration scheme detects a collision at most one time step after it happens. Contact detection with explicit and implicit methods will be further discussed in Section 5.2. The prediction of the Numerov implicit scheme is conservative and underestimates the maximum overlap but altogether is a better prediction than the ones produced by the explicit methods. Admittedly, the deviations will diminish if smaller time steps are utilized but the relatively large time step chosen here is to magnify and analyse the differences. Such deviations will still exist with smaller time steps but to a less degree. 5.2. Contact detection The objective of this validation test is to illustrate the fundamental difference between the explicit and implicit methods in detecting a contact at the first place. The explicit method employed for comparison with the implicit results in this test is the Position Verlet scheme presented in Eq. (7). The Verlet method to predict the position is a very common explicit method (Leach, 2001). Energy dissipation is excluded in this test and therefore both position Verlet and velocity Verlet methods lead to the same prediction of the particle’s position. Both Verlet and Numerov methods require two previous values in order to predict the next position which makes them an appropriate pair for conducting a fair comparison. Figs. 8 and 9 show how the implicit and explicit methods, respectively, respond to a particle approaching another object for example a wall. Snapshots 1–4 in both figures correspond to consecutive time steps. A relatively large time step is selected in order to illustrate the difference in response more clearly. Snapshot 1 is the

initial set up when the particle is moving towards the wall with an initial velocity. At the next time step, snapshot 2, the particle is close to the wall but with no contact yet. Here comes the difference between the implicit and explicit solutions. To calculate the new positions from here, the explicit method uses the current values which indicate no contact and decide the new position. The implicit solution uses the Newton–Raphson iterative method to find the new position and by each iteration it gets closer to the new position and therefore it discovers there will be a contact. At snapshot 2, the implicit solution already knows that there will be a contact in the next time step and it does consider this contact in calculating the new position. In other words, the explicit solution discovers a contact only after it happened while the implicit discovers the contact as it is happening. Due to this early detection of contact, the first overlap calculated by the implicit scheme is smaller than the first overlap in the explicit method as seen in snapshot 3 in Figs. 8 and 9. Obviously with a larger overlap, the explicit scheme pushes the particle more forcefully than the implicit scheme which results in a higher rise in snapshot 4. Apart from the size of the time step, the numerical error in this case also depends on how synchronized the time step and the collision start. To clarify this, Fig. 10 shows the fraction of a time step before the commencement of the collision. If the collision starts between the time step i and i + 1, the pre-collision fraction of the time step can be written as: tpre−collision =

tcollision−start − ti t

(41)

In order to discover the influence of tpre−collision on the prediction results, the simple case of a particle approaching a wall, shown in Figs. 8 and 9, is carried out with different values of tpre−collision . Prior to impact, the particle is moving towards the wall with a constant velocity of 0.5 m/s while the gravitational force is excluded. The value of the initial velocity does not have any effect on the duration of contact in the linear spring force model and consequently the resolution of collision and the results are therefore qualitatively independent of the impact velocity. In order to magnify the trend, a large time step corresponding to contact resolution of 1.2 is initially chosen. Fig. 11 shows how the predictions of the explicit and implicit methods depend on the pre-collision fraction of the time step. It shows that the exact time when a collision starts during a time step has less effect on the implicit method than the explicit method. This is explained by the fact that the implicit method can detect a contact during a time step in contrast to the explicit method which only checks for contacts at the start of the time step. Also deduced from Fig. 11 is that the worst prediction in both explicit and implicit methods is when the collision starts at the same time as the time step starts, i.e. when tpre−collision approaches zero or one. This is when the integration method is most “shocked” at the next

190

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

Fig. 10. Fraction of a time step before the commencement of the collision.

0.6

worst case scenario of collision commencement is attained, i.e. tpre−collision = 0. The results in Fig. 12 show that the differences between the explicit and implicit methods in detecting and solving a contact is mainly in very large time steps. With contact resolutions bigger than 5, the differences become negligible and both methods would give a very good prediction of the maximum overlap compared with the analytical solution.

Analycal Explicit Verlet Implicit Numerov

0.5

δmax /R [-]

0.4

0.3

5.3. Numerical stability 0.2

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

ΔTpre-collision [-] Fig. 11. Explicit and Implicit methods’ prediction of the maximum overlap occurring when a particle approaches a wall in relation to the pre-collision fraction of the time step, contact resolution = 1.2.

time step because considerable changes in the force are happening during the time step. Except for small test cases, there is no control, in general, over the start of the collisions in DEM. Therefore, the performance of the integration scheme in detecting the contact should be evaluated for the worst case scenario. Both explicit and implicit methods’ predictions would improve, of course, with increasing time step. Fig. 12 shows the effect of the time step on the maximum overlap predicted by explicit and implicit methods. The test is configured in such a way that the

As also discussed by Dziugys and Peters (2001), in addition to satisfy the required accuracy, a proper integration scheme is expected to be stable. This section explores numerical stability of implicit and explicit methods. Similar to the previous case, the Position Verlet and Numerov methods are the chosen explicit and implicit schemes and energy dissipation is excluded. 5.3.1. Particle sandwich A simple verification test is presented in this section where the analytical solution is easy to determine. The case chosen is a spherical particle in continuous contact with two parallel walls, as shown in Fig. 13. The distance between the two walls is equal to the sphere’s diameter. The sphere is initially positioned in contact with one of the walls which creates an initial overlap. In this type of “sandwich configuration”, the particle has no room for free motion and is continuously pushed back by the repulsive force from either of the walls. Such a continuous exposure to contact forces from different objects resembles the reality of many industrial applications involving dense granular matter. There are several features in the case described above which makes it easy to solve. These are:

0.6

• The motion of the particle is one-dimensional. It reverses its direction at maximum overlap points but remains on the same line of movement. • There is only one moving object. • The particle has only one contact at a time. It is either in contact with the lower wall or the upper one and not with both at the same instance. • The force model is linear.

Analycal Explicit Verlet Implicit Numerov

0.5

δmax/R [-]

0.4

0.3

0.2

0.1

0

1

2

3

4

5

6

7

8

Collision Resoluon [-] Fig. 12. Explicit and Implicit methods’ prediction of the maximum overlap occurring when a particle approaches a wall in relation to the time step, tpre−collision = 0.

Fig. 13. A particle sandwiched between two fixed walls.

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193 Table 1 Physical properties of the spherical glass bead used in the test and the corresponding duration of contact, linear spring force model. Value

Unit

Density,  Stiffness constant, k Diameter, D Coefficient of restitution, e

2520 250 8 1

kg/m3 kg/s2 mm

Duration of contact, Tcollision

0.005

s

Considering these simplifications, the equation of motion in Eq. (31) can be rewritten for this specific case, as shown in Eq. (42).

1

Analycal Explicit Verlet Implicit Numerov

0.9

Normalized Posion of Parcle [-]

Parameter

191

0.8 0.7 0.6 0.5 0.4 0.3 0.2

2

kij ∂ ri kij + r = (R + rj ) mi i mi i ∂t 2

(42)



ri = C1 sin



kij mi



t



+ C2 cos

kij mi



t

+ Ri + rj

(43)

The boundary conditions are assumed to be: ri (0) = Ri + rj + 0.15ri (44)

r˙ i (0) = 0

These initial conditions mean the particle is positioned, with zero velocity, experiencing an overlap with the upper wall equal to 15% of its radius. The value of initial overlap being 15% of the particle’s radius is chosen arbitrary and does not affect the generality of the case. This is due to the fact that the duration of collision in the linear spring force model is independent of the maximum overlap, as explained in Section 3. Applying the boundary conditions to specify the constants C1 and C2 , would lead to the particular solution to the differential equation:



ri = 0.15ri cos

kij mi



t

+ Ri + rj

(45)

The material chosen for both the spherical particle and the walls is glass with a density of 2520 kg/m3 and the stiffness constant is chosen to be 250 kg/s2 . Radius of 4 mm is assumed for the sphere and since no damping is assumed, the coefficient of restitution is 1. Based on these properties and using Eq. (29), the duration of contact is calculated to be approximately 0.005 s. These parameters are summarized in Table 1. The results are calculated for each time step by exact analytical solution and numerical solutions, both explicit and implicit, for the total duration of 0.45 s. Initially, two different time steps are used for the numerical solutions, one relatively small time step and one relatively big time step. The small time step is 0.0001 s, which is less than 50 times smaller than the duration of contact given in Table 1 and the big time step is 0.001 s which is more than 5 times smaller than the duration of contact. In other words each contact is solved in approximately 50 and 5 time steps respectively for the small and the big time steps. The comparison of the results are performed in the arbitrary chosen range of 0.4–0.45 s. The very start of the simulation is not chosen to be analysed because as the time advances, the numerical error of the integration is being added up and the errors become more representative. The results presented

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Normalized Time [-] Fig. 14. Prediction of the particle’s position in the sandwich setting by analytical, explicit and implicit methods with a relatively small time step, contact resolution = 50.

are normalized over the maximum value for the position of the particle and time since the absolute values do not make sense in the comparison. Fig. 14 shows the results of the analytical, explicit and implicit solutions for the small time step. It is observed in Fig. 14 that at such a small time step both explicit and implicit methods provide a good prediction of the particle’s position compared with the analytical solution. As a next step, the comparison is carried out for the large time step and the results are shown in Fig. 15. It is observed that the results of the explicit prediction are shifted away from the analytical solution at such a large time step while the implicit results are still consistent with the analytical solution. What is seen in Fig. 15 is that the prediction of the explicit method of how the particle moves is quite accurate but there is a delay in the response of the explicit method. In other words, the explicit prediction is slightly behind the analytical and implicit solutions in time. Such responses from the explicit and implicit methods are consistent with the theory, too, recalling that the explicit method remains insensitive to the changes happening during a time step and is in fact one step behind the reality. The

1

Analycal Explicit Verlet Implicit Numerov

0.9

Normalized Posion of Parcle [-]

where subscript i represents the moving particle and subscript j represents the wall in contact with the particle. The right hand side of the equation is constant. The general solution to such a second order linear ordinary differential equation, Eq. (42), is shown in Eq. (43) (Boyce & DiPrima, 2000; Kamke, 1977; Polyanin & Zaitsev, 2003).

0.1

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Normalized Time [-] Fig. 15. Prediction of the particle’s position in the sandwich setting by analytical, explicit and implicit methods with a relatively large time step, contact resolution = 5.

192

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

implicit method being numerically stable even with large time steps.

Fig. 16. Comparison of the relative error of the implicit and explicit methods in prediction of the particle’s position in the sandwich setting with different time steps.

implicit method does not “freeze” as such since the position is updated at each iteration during a time step. While the explicit method’s prediction of the maximum overlap is quite accurate with the contact resolution of 5 as shown in Fig. 15, it will deviate from the analytical maximum overlap as the time step increases. At the relatively very large time step of 0.0035 s, which is in fact close to the duration of contact, the explicit method totally failed with the particle moving outside the walls. The implicit method’s error with this time step was still bounded. There is a limit, however, for the implicit method as well in the size of the time step as it cannot exceed the duration of collision. The test was performed with other time steps too. Fig. 16 shows the biggest relative errors in prediction of the particle’s position by explicit and implicit methods for different time steps. Understandably as the time step goes up, the error in both implicit and explicit methods increase but the error in the implicit method remains bounded. 5.3.2. Hopper discharge In order to verify such a difference between implicit and explicit methods in error propagation, a hopper discharge sample case is run. The problem includes 153 spherical glass beads in a circular hopper. Two different particle diameters of 3 and 4 mm are used which leads to a poly-disperse problem. Small and big particles are visualized using different colours in the animation of results for better illustration of the discharge flow. The mechanical and contact properties of the particles are summarized in Table 2. The case is run using a relatively large time step of 0.0005 s which corresponds to the contact resolution of 2.4. The video components (1) and (2) in Appendix A show the animation of hopper discharge solved by explicit and implicit methods, respectively. The exact same boundary conditions are used in both explicit and implicit runs. The video component (1) in Appendix A demonstrates the failure of explicit integration in predicting the flow of particles with big time steps. Particles start to unrealistically blow up and pass the walls due to formation of large overlaps and consequently high velocities. On the other hand, the prediction of the implicit method demonstrated in video component (2) in Appendix A is good. No particle passes the wall and the discharge flow is normal and physically realistic. Such a response confirms the theory proving the Table 2 Physical properties of the spherical glass particles in the hopper discharge simulation and the corresponding duration of contact, linear spring force model. Parameter

Value

Unit

Density,  Stiffness constant, k Diameter, D Coefficient of restitution, e

2520 250 3 and 4 1

kg/m3 kg/s2 mm

Duration of contact, Tcollision

0.0012

s

6. Conclusions Several conclusions can be drawn from the results of the comparisons and verification tests presented in this study which are summarized below. • Implicit integration schemes are more robust than the explicit schemes in large time steps, i.e. contact resolutions smaller than 5. Earlier detection of contacts and smoother transition between different phases of collision results in errors being bounded in implicit methods. However, the size of the time step is also limited in implicit methods as it cannot be bigger than the duration of collision. • In smaller time steps, i.e. contact resolutions bigger than 5, the gain in accuracy in implicit methods is too small to have a decisive impact on the global behaviour of the granular system. • Though implicit integration allows bigger time steps than explicit methods, it still remains a risky approach to utilize big time steps. This is because the duration of collision in non-linear impact models is not absolutely predictable from the material properties. Therefore, the time step selected should always consider safety margins for uncertainties in impact velocities. Besides, the implicit solution risks non-convergence of the Newton method in large time steps. This issue, though, can be managed by adapting the tolerance of the Newton method according to the convergence history. This approach, however, does sacrifice some degree of accuracy. • Effective detection of contacts in both explicit and implicit methods does not only depend on the size of the time step but also on how late a collision starts relative to the start of the time step. Since there is no reasonably efficient way in DEM to synchronize the start of the time step with the start of each collision, the worst case scenario should always be considered when evaluating a method. The implicit method is less sensitive to pre-collision fraction of the time step and is therefore more robust than explicit methods in contact detection. • Implicit methods are computationally more intensive than explicit methods. A system of non-linear equations needs to be solved at each time step iteratively. This makes implicit methods

K. Samiei et al. / Computers and Chemical Engineering 49 (2013) 183–193

too expensive for the benefits. There are methods to improve the efficiency of the Newton method such as less frequent Jacobian assembly or less frequent collision detection. However, such methods should be employed with caution because they also affect the accuracy of the solution while reducing the computations. • The examples presented in this study confirms the promising potential of the implicit integration methods in the context of DEM. However, with the current state of the art where the computational costs of force computation is the main drawback of the DEM, explicit method will yet remain the convenient method of integration within DEM community. Future research should focus on improving the implicit algorithms to efficiently reduce the computational overhead. Acknowledgement Experiments presented in this paper were carried out using the HPC facility of the University of Luxembourg. Appendix A. Supplementary Data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compchemeng. 2012.10.009. References Babic, M., Shen, H. H., & Shen, H. T. (1990). The stress tensor in granular shear flows of uniform, deformable disks at high solids concentrations. Fluid Mechanics, 219 Barrett, R., Berry, M., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., et al. (1994). Templates for the solution of linear systems: Building blocks for iterative methods (2nd ed.). Philadelphia, PA: SIAM. Boyce, W., & DiPrima, R. (2000). Elementary differential equations (7th ed.). New York: Wiley. Campbell, C. S. (2002). Granular shear flows at the elastic limit. Fluid Mechanics, 465, 261–291. Cleary, P. W. (2004). Large scale industrial DEM modelling. Engineering Computations, 21, 169–204.

193

Cundall, P., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Geotechnique, 29, 47–65. Dziugys, A., & Peters, B. (2001). An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers. Granular Matter, 3, 231–265. Ferziger, J. H., & Peric, M. (2002). Computational methods for fluid dynamics. SpringerVerlag. Gear, C. W. (1966). The numerical integration of ordinary differential equations of various orders. Technical Report Argonne National Laboratory. Gear, C. W. (1971). Numerical initial value problems in ordinary differential equations. Englewood Cliffs: Prentice-Hall. Hairer, E., Lubich, C., & Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer/Verlet method. Acta Numerica, 12, 399–450. Jasion, G., Shrimpton, J., Danby, M., & Takeda, K. (2011). Performance of numerical integrators on tangential motion of DEM within implicit flow solvers. Computers and Chemical Engineering, 35, 2218–2226. Kamke, E. (1977). Differentialgleichungen: Lösungsmethoden und Lösungen, I, Gewöhnliche Differentialgleichungen. Leipzig: B.G. Teubner. Ke, T.-C., & Bray, J. (1995). Modeling of particulate media using discontinuous deformation analysis. Journal of Engineering Mechanics, 121, 1234–1243. Ketterhagen, W. R., Curtis, J. S., & Wassgren, C. R. (2005). Stress results from twodimensional granular shear flow simulations using various collision models. Physical Review, E, 71. Kozicki, J., & Donzé, F. V. (2009). YADE-OPEN DEM: an open-source software using a discrete element method to simulate granular material. Engineering Computations, 26(7), 786–805. Leach, A. R. (2001). Molecular modelling: Principles and applications. Prentice Hall. Negrut, D. (1998). On the implicit integration of differential algebraic equations of multibody dynamics. PhD thesis University of Iowa. Polyak, B. (2007). Newton’s method and its use in optimization. European Journal of Operational Research, 181, 1086–1096. Polyanin, A., & Zaitsev, V. (2003). Handbook of exact solutions for ordinary differential equations (2nd ed.). Boca Raton: Chapman and Hall/CRC. Schäfer, N., & Negrut, D. (2010). A quantitative assessment of the potential of implicit integration methods for molecular dynamics simulation. Computational and Nonlinear Dynamics, 5. Skeel, R. D., Zhang, G., & Schlick, T. (1997). A family of symplectic integrators: stability, accuracy, and molecular dynamics applications. SIAM Journal on Scientific Computing, 18, 203–222. Thompson, P., & Grest, G. (1991). Granular flow-friction and the dilatancy transition. Physical Review letters, 67, 1751–1754. Tuley, R., Danby, M., Shrimpton, J., & Palmer, M. (2010). 2010. Computers and Chemical Engineering, 34(6), 886–899. Walton, O. R., & Braun, R. L. (1986). Viscosity, granular-temperatures, and stress calculations for shearing assemblies of inelastic, frictional disks. Rheology, 60, 949–980.