An efficient calculation method for particle motion in discrete element simulations

: ..... ~ .-~

..

POWD TECHNO ELSEVIER

P o w d e r T e c h n o l o g y 98 ( ! 9 9 8 ) 2 2 3 - 2 3 0

An efficient calculation method for particle motion in discrete element simulations D. Zhang, W.J. Whiten * ,lulius Kruttschmitt Mineral Research Centre, University of Queensland, Indooroopilly. Queenshmd 4068. Austraha

Received I December 1997: received in revised fornl 30 January 1998: accepted 18 March 1998

Abstract

Discrete element models can be widely applied in simttlation of many types of physical systems including mineral processing and mining. The simulation process requires the calculation of the motion of each discrete element. The existing calculation methods are time consuming and often have been simplilied to give unrealistic behaviour. A fast calculation method using a realistic force model has been developed. This method separates the solution of the differential equation of motion from the discrete element simulation process, it can re:;ult in improved accuracy and speed in certain types of simulations. © 1998 Elsevier Science S.A. All rights reserved. Kevwords: Mineral processing: Discrete element simulation

1. I n t r o d u c t i o n

Discrete element models can be widely applied in simulation of mineral processing and m!ning. Some applications in ball mill modelling and other grinding mill modelling have been investigated by other researchers [ 1.2 ]. Discrete element models simulate the physical process by dividing the system to be simulated into a number of discrete particles and by dividing the simulation time into small time steps. Given the initial velocity and position of the particles, the particles' velocity and position are updated at each time step. The procedure can be repeated until the requited inlbrmation is obtained from the simulation. The motion of impacting particles can be calculated in two ways: the coefficient of restitution method and the force summation method. When two bodies collide, the coefficient of restitution e is defined as the negative ratio of the relative speed alter the collision to the relative speed before collision:

ever, experimental results show that the coefficient of restitution for normal impacts is a function of the impact velocity [3]. Also, the coefficient of restitution is a function of temperature for some materials. By using the coefficient of restitution, the particle velocity after imp:let can be easily calculated. Velocity calculations using the coefficient of restitution assurne an instantaneous collision. When the collision time is long compared with other events, such as a second collision or changes in other forces, calculation using the coefficient of restitution is not apprc 9riatt. The alternative method is to integrate the sum of the contact forces to determine the motion of the particles [4-6]. The lbrce of each contact point is calculated by a contact lbrce model. The velocity and displacement are then updated by the sum of the total forces acting on the particle as fi)liows' •f'""'

(2 )

m

V,,+ ~= V,, +a,,At e=-

I'il

--

(!)

I" b

The coefficient of restitution has some value between zero and unity. According to Newton's theory, as described in many physics text books, the coefficient of resti!ution is a constant value depending on the impact body material. How* Corresponding author.

X,,+ ~=X,,+

V,,+V,,, ~At 2

(3)

(4)

The total three is usually assumed constant during the updating time step. A low order integration method as given above is generally used as the contact Ibrces make both discrete changes and rapid changes, which makes the problem unsuited to higher order methods. It has been found that most

0032-5910/98/$ - see front matter ~) 1998ElsevierScience S.A. All rights reserved. PHS0032-5910(98)00055-2

........ ~+~J~

~!!~~~~J~i~!~!~!~!i ¸~~:¸¸i~i!i'!~! ~ili~i!!~!i,¸!!i:!i~i!i!i!i~ii!i!i!i!~ii~ii~!ilil

D. Zhang, W.J. W h i t e n / P o w d e r Technology 98 (1998) 223-230

of the existing models using this method use an incorrect model for the contact forces [71. Furthermore, this method requires a huge amount of calculation during the simulation process to obtain accurate results. If the amount of calculation could be reduced the time required for the simulation process would be significantly smaller or alternatively larger simulations could be attempted. In this paper, a third method using velocity and position summation to calculate the particle motion in normal direction is introduced. For an accurate calculation, the particle motion calculation is divided into normal motion calculation and tangential motion calculation. This paper only discusses the normal motion calculation, and the tangential motion calculation will be described in a following paper. The contact equations are integrated numerically once and dimensional analysis is used to construct a two dimensional table which stores the pre-integrated solution. The two dimensional table method separates the integration of the contact force equation and the simulation process. Therefi~re, the time required for the simulation process is reduced.

It is noted that the contact forces are both significantly larger and change more rapidly than the other forces and hence they are treated in a separate manner. The effect of the gravitational force can be determined as in the force summation method, i.e.:

aV~=gAt 1

(lO) ,

AXe= ~gAt-

( !1 )

Similarly the velocity and displacement caused by other external forces can be calculated as: AV,- = ~ A t

(12)

AX,,= ]'~'~'At" 2m

(13)

m

Therefore, combining the external and gravitational as:

a,,=g+t~hn

(14)

The total velocity and displacement (Eqs. (6) and (7)) are now given as: 2. Using velocity and displacement summation simulation

V,,+ i=V,+.%V,.+a,At

(15) i

The velocity and displacement summation method takes the force equation:

d'X(t) m dt " =.L +f. +l'~

(5)

where f,. is the contact force,.t~, is the external force, andre, is the gravitational Ibrce. Integrating the components gives: t~ A t

,f

V,, , , = V,, + - tll

At

t,

,f

f dt + -ttt

!

,f

l~,dt + - " tit

.1"~dt

X.'l-X,,+v,,m+l f,,,f ,/~.

dt dt

'ff

'

17)

t

' If t

,,,

.I~.dt dt

ttt

t

f

(17)

The second term in Eq. (6) and the third tern1 in Eq. ( 7 ) art_, the velocity change and displacement change caused by the contact forces. They are denoted as follows: r ¢ Ar

av,= l,. j" f,d,

To calculate the forces between particles in the normal direction is typically done by assuming spring and dashpot components at the contact points. Relative motion of the particles using the spring and dashpot model can be described by the following differential equation based on Hertzian contact theory 18 !:

t ~ . ~ t t + Xt

,l~,dtdt + --

lit

(7). 3. The spring and dashpot contact model

r~xtt,~t

+ --

The velocity and displacement can be calculated by Eqs. (14)-(16) when the velocity and displacement changes caused by the contact force are known. The tbllowing sections discuss a calculation method for the effect of the contact force using dimensional analysis and the integrals in Eqs. (6) and

(6 )

!

and integrating again gives the change in position as the sum of double integrated of the threes.

'

(16)

t~At

I

t ÷At t~.~t

X,+ ,=X,,+V,,At+AX,.+ 2 a,,At 2

~8)

The plots of velocity t. vs. displacement .r using Eq. (17) with different initial velo,:ity is plotted in Fig. I. It can be seen that the curves corresponding to different initial conditions all have a similar shape.

t t + ..%t t + A t

I

I

4. Dimensional analysis Eq. (17) can be normalised to the following form with the initial conditions x(0) = 0 and t,(O) = 1"

D. Zhang. W.J. Whiten/Powder Technology 98 t ! 998) 223-230

225

0.8 0.6 0.4" 0.2"

° .2 ° 4 ° .6 °.8

-

0.15

l'

Fig. 1. Velocity vs. displacement for different initial conditions (m = 10, k = 2000, a = 0.4).

(18)

where x( t ) = (mvJVk) 2/~x*( t* ), t = ( (mt',,'-/k)'-~s/v.)t*,

Fig. 3. z Against c* using Eq. (19) (a =0.4). to the normalised values. Thus z(t) is defined as a function of t* and parametric plots of z(t*) vs. v*(t*) and z(t*) vs. x*(t*) can be obtained from the solution of Eq. (18). v*(t*)

=

/r(t)2 + .k.r(t)sn

c( t ) = t',~t'( t'* ) The curves shown in Fig. I therefore all can be normalised

into one curve with the v axis divided by Vothe initial velocity and the x axis divided by (mt,o2/k) ~'/~, shown in Fig. 2. When using the discrete element method, it is required to discretize time into small time intervals to calculate and update the particle position at next time step by the current velocity and displacement. Using the normalised equation will simplify the calculation process. As shown in Fig. 2, with a given a, Eq. (18) gives the normalised v-x curve. This curve can be denormalised to the original curve with the given initial velocity. Therefore, it is necessary to know the initial impact velocity of each impact. However, the initial impact velocity is unknown because the particle velocity is the result of the total impacts on the particle. Using the particle current velocity and position (t', x) which are always known, the initial impact velocity for each contact can be derived. Therefore, the particle velocity in next time step will be determined by the initial velocity, current velocity and current position. To calculate the initial velocity, a new parameter z is chosen. The new parameter z is independent of initial velocity and can be used to convert from actual velocity and position

r(t)

z(t) =

(19)

~/v*(t*)2 + x * ( t * ) "~/~"

m

For a given a, the trajectory (v*, x*) can be calculated by Eq. (18). Fig. 3 shows the relationship between z(t*) and t,*(t*) while Fig. 4 shows the relationship bet~'een z( t* ) and x* ( t* ) . Therefore, at the beginning of each simulation, the normalised differential equation (Eq. (18)) with given a is solved numerically. The value of z(~*) for each point is calculated. The values z(t*), v*(t*) and t* are stored and used for each step calculation. The velocity and displacement change caused by the contact force can be calculated by the dimensional analysis method using the following steps: 1. Calculate the current z,, from the current relative velocity t',, and relative displacement x, using Eq. t 19). 2. Find the v,,* from the z,, value using a table for the curve in Fig. 3. 3. Calculate the initial velocity v,, as t'o t',/t,,,* or if v, is too small from co= ~/(x,/x,,* )5/2k/m. 4. Calculate t,,* from v,,*. 5. Calculate the normalised time step size At* from the initial velocity and hence t,, + ,*. =

1' 0.8

0

0.6 o.¢

o

0.2

!

.

.

.

.

"

.

.

.

.

"

-o-°iI --0.

.

.

.

.

'

.

.

.

.

.

'

"

I

' --

Fig. 2. Nonnalised velocity vs. normalised displacement (a = 0.4).

-o . ~ -o 4" -o 6" -0

0.2

!

!

0.

.

4

~

8"

Fig. 4. z Againstx* using Eq. (19) ta=0.4).

x

226

D, Zhang, W,J. Whiten/ Powder Technology 98 (1998) 223-230

appropriate grid points. The particle velocity and displacement at next time step are then simply calculated by incrementing the current values. The time step size used in the. simulation is preferably equal to the time. step size used in the table. For the time steps that are a multiple of the table step size, the table can be used multiple times to give the total effect of the time step and when the time step size used in the simulation is not an exact multiple of the time step size used in the table, the remaining part is calculated by linear interpolation.

6. Obtain the normalised velocity c,, + ~* and displacement .~,,+ ~* at t, + ~* from the stored solution of Eq. ( i 8) with given a. 7. Calculate the velocity change At,,, = v , + ~ - t , , and displacement change A~,,=x,,+ ~ - x , at current time step using the normalised velocity c, + ~* and displacement .~,,+ ~*, and the initial velocity co. 8. Calculate the velocity change AV, and displacement change AX, of each particle from At,,, and Ax,,.

5. Velocity and position increment table

6. Examples To implement a simulation using the relations of the previous sections, a two dimensional new motion table is created to save the change in velocity and displacement lbr the next time step. Instead of solving the motion equation to update the velocity and displacement, the velocity and displacement are updated by reading the velocity and displacement changes for the current time step from the new motion table. This method separates the generation of an accurate solution of the differential equation from the simulation process. As reading from a table is simpler than calculating the contact forces and a much larger time step can be used, the simulation speed is increased. The new motion table has displacement and velocity axes. It is generated by the dimensional analysis calculation described in Section 4. Each type of contact has its own table. This table contains velocity change and displacement change at each grid point, as shown below:

Two examples are used to demonstrate the differences between the new velocity summation method and the force summation method. These examples include free motion, single contacts and multiple contacts. An accurate solution to these problems was obtained by integrating using Numerical Recipes routines [9l to solve the equations tbr each section of the motion. As the two simulation methods behave differently with regard to step size they were both run with a range of step sizes. Step sizes that give roughly similar accuracy were then chosen so that the simulation times could be compared. 6. I. Exanqde h one ball bmmces on a moface with gravity

These objects have the lollowing properties ( input data to the program ):

I

At=

dt = c,,,, - t',,

(21))

III

-.L,-c,,At

Object

X,

X,

V,

V,

Mass

Radius

A B

20.9 t)-3().0

9.5 1)

0 0

- 1.0 0

4(1 le20 (inlinity)

I N/A

( 21 I

Contact paralr, eters: stiffness ( k ): 2000, damping ( a ) : 0.4,

"It) update the particle motion, the velocity and displacement change for the next time step are interpolated using the

gravity (g): 9.8. Pantmeter used in solution table: dt: 0.01, dr: I, dx: 0.05. The simulation results from the two simulation methods are shown in Fig. 5. in Fig. 6, the step sizes lbr the two

Ax=

dtdt=x,,.,, Ill III

I() -

7~

'

velocitysummation method At =0.01

5'

E

-2.5 0

! 2

I 4

Time Fig. 5. Displacementvs. time (example I ).

D. Zhang, W.J. Whiten/Powder Technology 98 (1998) 223-230

0.02

force summation method

227

At =0.0001

velocity summation method At =0.01

0.015

0.01 0.005 Ua 0 -0.005

-0.01 0

I

I

I

2

4

6

Time Fig. 6. The error in displacement lbr the two methods in Fig. 5.

,,s

Correct

20

•

At =0,001

o

At =0.01

+

At =0.1

++ +++ 4"

4"÷+~,¢,N~+4"+

15

•

4"+

,o'p~ I "k

s-! 0

++

++

4

-5" 0

++

4"

%

4-

++

4"

4"

4"

+

4" 4"

4"

++

4"+

4"

+

+

•

= 4":--

-

I

I

I

I

2

4

6

8

Time Fig. 7. Displacementvs. time using force summationmethod tbr different time step sizes (example I }. methods have been chosen to give similar accuracy. The velocity summation method required 3.3 s to make up the new motion table and 5. l s to finish the calculation, and the force summation method required 65 ! s for the calculation to produce a slightly less accurate result. For comparison integration using Numerical Recipes routines [9J the calculation for one collision took 3. ! s. The force summation method assumes the tbrces acting on the objects are constant within the time step. However, the contact force increases rapidly until it reaches the maximum and then reduces to zero during the contact process. It is a continuous process. Thus, the force summation method eliminates the smooth change of the force within the updating time step. As shown in Fig. 7, the smaller the time step used, the smaller the error that occurred in the calculations. Furthermore, the assumption of a constant force during the time step is less accurate when objects change from non-contact to contact and contact to non-contact processes since the changes in force magnitude in these two steps are much larger than in the other steps when the objects are in contact. At the step in which objects change from non-contact to contact, the force summation method assumes that there is

no contact force acting. The object's motion is updated to B instead of the real position B' shown in Fig. 8. It can be seen in Fig. 8 that the calculated initial velocity then is higher than the actual initial velocity. Therefore, the calculation trajectory lies on another curve. The particle gains energy from each

A

m

v

X

Fig. 8. An example to show the calculation of the step from non-contact to contact in force summation method.

,¸!¸¸¸!'!i !i!!!!!!!!iiii !!i!!i!!i!i!i!!:ii!li!i!i!iii!!i!i!!ii!!ili

=: D. Zhang, W.Z Whiten I Powder Technology 98 (1998) 223-230

228

impact. ~ i s can be seen in Fig. 7 for the simulation with time step size 0.1. In the contact to non-contact step, the contact force changes to zero and contact objects separate. Again, the force summation method adds an extra energy to the calculation in this step. The error in the force summation method can be reduced by using a very small time step. As shown in Fig. 7, the simulation error reduced dramatically with smaller time step size. However, the simulatio~l time increased dramatically as the time step size decreased; as the following results show: Time step size Simulation time (s)

0,1

0.01 7.7

2.0

0.001 65.9

for a range of time step sizes. To use a smaller time step size might require an additional table with a finer v - x grid which would result in a slightly longer simulation time. 6.2. Example 2: two balls bounce on a surface with gravi~.

These objects have the following properties (input data to the program)-

0,0001 651

The velocity summation method considers all the changes in forces during the calculation time step as the integrals of the three are used. For this method, the new motion table must be constructed initially which takes some time. However, this table can be used lbr the whole simulation process and only needs to be constructed once. In complicated simulation cases, a huge number of simulation steps are required. The time saved in the calculation of each step not only compensates for the time used in constructing the solution table but also can very significantly shorten the time tbr large simulations. The calculation times for three step sizes are given below: Time step size Time to construct the table (s) Calculation time I s )

O. I 3.8 0.7

0,01 3.3 5. I

Object

X,

X,.

V,

V,,

Mass

Radius

A B C

20.9 20.9 0-30.0

! 3.5 11.5 0

0 0 0

- 1.0 - 25.0 0

4 40 le20 (inlinity)

0. ! 0.5 N/A

Contact parameters:

A-B B-C

Stiffness (k)

Damping ( a )

50 000 50 000

0.9 0.4

Parameters used in the solution table:

A-B B-C

0.001 3.5 50.6

dt

dr

da

0.001 0.001

0.5 0.5

0.1105 0.005

Again, the time steps have been chosen so that the two methods give similar accuracy, as shown in Figs. 10 and 11. However, the velocity summation method shows an obvious advantage in simulation speed in this example. The force summation method required 10 394 s to finish the calculation while the velocity summation method spent 16.4 s to construct the solution table and 76.4 s Ibr the calculation. A multiple contact happened at the end part of example 2. Multiple collisions require exact timing for accurate simulation and are thus very sensitive to small variations. With the chosen time steps, both methods gave similar accuracy in this region.

The velocity summation method gives reasonable solutions regardless of the time step size. However, as seen in Fig. 9, the simulation with time step size 0.001 gives less accuracy. It is because the time step size used in the calculation should be larger than or equal to the time step size defined in the new motion table (time step size used in the solution table is 0,01 in this case). This is caused by the error in interpolating the solution from the table to obtain the required time step length. Each solution table is onl~ suitable

he tab!e step size is 0.01 ?5"

correct

.~" 8 -23"

-2.5

I

0

2

!

Time

4

+

At = 0.001

•

At =0.01

o

At - 0 . 1

I

•

6

s

Fig. 9. Displacements vs. time using velocity summation method with varied time step size (example 1 ).

D. Zhang, W.£ Whiten / Powder Technology 98 ( 19981223-230

~

.

¢"

correct

'~. \

/ _,

•

229

force summation method At =O.0000l

15

/ t

~

'

velocitysummation method

....,

5-

O-

-5 0

U 5

275

"' 715

Time Fig. 10. Displacementvs. time (example 2 ).

r.

force sun|nlalion method

(}.()0001

At =

qiAIh

............o ...........

velocity sumnlalion

method

AI = 0 , 0 0 1

O.O4

0.02

I

-

-I).(12

-0.04

I

0

2

.

.

I

.

.

.

.

.

.

.

4

!

6

Time Fig. I !. Error in displacement liarFig. 10.

7. Conclusions The velocity summation method is robust in the sense that changes in the size of time step have little effect. The solution table is constructed using a very small time step size and a higher order method in integration, which gives an accurate calculation result. Higher order and adaptive integration methods can be used for this integration as it does not involve time steps that have a change from non-contact to contact or contact to hen-contact. Different time step sizes normally only affect the ~imulation accuracy slightly. This allows the user have the flexibility to choose the time step unlike in the force summation method where time steps that are too large result in an energy gain. The velocity summation method is much faster than the force summation method in both examples for similar accuracy. The main reason for the reduced computation time is that a much larger time step can be used to obtain the same accuracy. Even in the small examples given, the time taken to prepare the update table is much smaller than that saved

by its use. It can be deduced that the velocity summation method will significantly decrease the time needed for a simulation in complicated simulation cases.

8. List of symbols

api

dt

dy e

f g k m t

Damping coefficient ( In the non-linear equation, it is an empirical constant that can be related to the coefficient of restitution. ) Acceleration at step n Time step size used in the new motion table Width of cells used in the new motion table Height of cells used in the new motion table The coefficient of restitution Force Acceleration due to gravity Stiffness coefficient Particle mass Time

D. Zhang, W.J. Whiten/Powder Technology 98 (1998) 223-230

230 At

V l" l'a l' b UO

At,

X x

~x

Time step Absolute velocity of particle Relative velocity Relative velocity after collision Relative velocity before collision Initial relative velocity Velocity change Absolute displacement of particle The amount of normal contact overlap The change of the amount of normal contact overlap The parameter used to convert to normalised values.

badicators and subscripts * c e g n n+ I x y

Normalised value Value related to the contact three Value related to the external force Value related to the gravity Value at step n Value at step n + I x Coordinate y Coordinate

Acknowledgements Authors would like to thank the JKMRC management for support of this project.

References [ ! ] B.K. Mishra, R.K. Rajamani, The di~cre',e element method for the simulation of bali mills, Appl. Math. Modelling 16 (1992) 598-604. [ 2 ] T. Inoue, K. Okaya, Designing the centrifugal mill by discrete element method simulation, Proceedingsof the XX International Mineral Processing Congress, 1997, pp. 41 !-421. [ 3 ] C.V. Raman, The photographi~ study of impact at minimal ve|ocities, Phys. Rev, 12 (1918) 442--447. [4l B.K. Mishra, Study of media mechanics in tumbling mills by the discrete element method, PhD Thesis, The University of Utah, 1991. [ 51 P.A. Cundall, O.D.L. Strack, A discrete numerical model tbr granular assemblies, Geotechnique 29 ( 1) (1979) 47-65. 16] R.E. Barbosa-Carrillo, Discrete element models for granular materials and rock masses, PhD thesis, The University of Illinois, 1985. 171 D. Zhang, W.J. Whiten, The calculation of contact forces between particles with spring and damping models, Powder Technol. 88 (1996) 59-64. [81 Y. Tsuji, T. Tanaka, T. ishida, Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe, Powder Technol. 71 (1992) 239-250. [91 W.H. Press et al., Numerical Recipes ir. C: The Art of Scientilic Computing, Cambridge Univ. Press, 1992.