- Email: [email protected]
Step size control for efficient discrete element simulation
0892-6875(01)00148-0
Minerais Engineering, Vol. 14, No. 10, pp. 1341-1346,2001 © 2001 Elsevier Science Ltd Ali rights reserved 0892-6875/01/$ - see front matter
STEP SIZE CONTROL FOR EFFICIENT DIS CRETE ELEMENT SIMULATION*
D. ZHANG and W.J. WHITEN Julius Kruttschnitt Mineral Research Centre, University of Queensland, 4072 Australia Email: [email protected] (Received 9 April 2001; accepted 31 May 2001)
ABSTRACT
The step size determines the accuracy of a discrete element simulation. The position and velocity updating calculation uses a pre-calculated table and hence the control of step size can not use the integration formulas for step size control. A step size control scheme for use with the table driven velocity and position calculation uses the difference between the calculation result from one big step and that from two smalt steps. This variable time step size method chooses the suitable time step size for each particle at each step automaticalty according to the conditions. Simulation using fixed time step method is compared with that of using variable time step method. The difference in computation time for the same accuracy using a variable step size (compared to the fixed step) depends on the particular problem. For a simple test case the times are roughly similar. However, the variable step size gives the required accuracy on the first run. A fixed step size may require several runs to check the simulation accuracy or a conservative step size that results in longer run times. © 2001 Elsevier Science Ltd. Alt rights reserved. Keywords Modelling; simulation
INTRODUCTION The discrete element simulations model the behaviour of fluid and solids by using an assembly of discrete elements. This method can be widely applied in mining and mineraI processing areas. The discrete element models divides the simulation time into small time steps. For each time step the velocity and position of each particle are updated using the equation of motion and the forces acting on the particle. The size of the time step is critical to the calculations since it determines the stability and accuracy of the calculation in the discrete element model. A time step size which is too small results in an unnecessarily long simulation time. However, a time step size which is too large produces incorrect simulation results.
• Presented at Comminution '01, Brisbane, Australia, March 2001
1341
1342
D. Zhang and W. J. Whiten
According to Cundall and Strack (1979), the following relation must be satisfied to give stable calculation:
I1t < 2.Jm/ k,
(1)
k m
where
11 t
stiffness coefficient, particle mass, is size of time step.
Tsuji (1993) found that the calculation became unstable when the time step near the limit value given by this equation. Therefore, Tsuji proposed that the time step should be proportional to the oscillation period of the spring-mass system. The oscillation period of the system is given by
I1t < 27r.Jm / k
(2)
Tsuji suggested a suitable time step can be obtained by dividing one half of the natural oscillation period by a small integer n, eg 5. Both Tsuji's and Cundall's time step equation have a similar form with Tsuji recommending a slightly smaller value. However, the particle separation time is less than one half of the oscillation period, because the impacting particles will separate when impact force equals to zero (Zhang and Whiten, 1996), and a second collision or changes in other forces could happen within this time step to affect the particles impact. A smaller time step may be needed to obtain accurate calculation of the particle motion. Often it is necessary to check the simulation results using different step sizes. Zhang & Whiten (1998) introduced a method of using accurately integrated (Press et al., 1992) tables of updates to the velocity and position of each discrete element. A step doubling scheme to control the accuracy of this velocity and position updating by varying the step size has been developed. This method allows the time step size to be adjusted to a suitable size automatically according to the calculation error in each part of the simulation. The step doubling scheme calculates each simulation step twice: firstly, it calculates the simulation results for the whole step, and secondly, it calculates the results by using two half steps. The calculation error for each step is defined as the difference of the simulation results from one full step and the results from two half steps. If the calculation error is within a predefined error tolerance, the simulation moves on to the next step. The step size used for the next step is related to the calculation error of the previous step. If the calculation error is larger than the error tolerance, the step size is reduced and the calculations repeated until the error is within the error tolerance. The step size needed to give the required accuracy typically varies in different parts of a discrete element simulation and at different times in the simulation. The extra computation needed to check the accuracy of a given step size can often be justified by the reduction in computation obtained by being able to use larger steps in parts of the calculation. In addition, the controlled accuracy avoids the need to choose a conservative step size and to make multiple runs to check the accuracy. Both of these require additional calculation.
STEP DOUBLING SCHEME USED IN DISCRETE ELEMENT SIMULATIONS The step doubling scheme used in discrete element simulations is as follows: 1.
Predefine the maximum acceptable error (error tolerance), maximum time step size and minimum time step size. 2.· Calculate the velocity and displacement of every object using the time step determined from the previous step. 3. Calculate the velocity and displacement of every object using two half time steps. 4. Calculate the difference of the two results and compare with the predefined maximum error; for the particles which have the error larger than the predefined error, repeat steps 2 and 3 using half of the old time step as the full time step.
Step size control for efficient discrete element simulation
5. 6.
1343
Calculate a more accurate value for the motion by combining the simulation results from two half time steps. Determine the time step size for the next time step from the calculation error. If the error is much less than the predefined error, a larger time step size will be used for next time step.
This scheme can be applied to any method of updating the velocity and displacement of the simulation objects but is particularly suited to the table driven method developed in Zhang and Whiten (1998), where the updates are not made using a numerical formula for the solution of differential equations. The repeat of calculation steps 2 & 3 from step 4 needs to be done only for the discrete elements that do not satisfy the accuracy test. For these particles multiple steps are calculated until all particles in the simulation are updated to the same time. The time step size is constructed by a time step base times a time step factor which is determined by the overall conditions of the contact. The time step base remains the same for the whole simulation. Local time step factors determined by the calculation error are applied to determine the time step size of each element. For the particles that are not in contact, the difference between one full step and two half steps is very small. So a larger time step size can be applied. Therefore the step size control scheme can avoid unnecessary use of small time step sizes. For the particles that are in contact, a smaller time step size must be used as the calculation error becomes larger especially when the particles change from non-contact to contact and contact to non-contact. For practical usage, the time step sizes used in the step doubling are limited between predefined maximum time step size and minimum time step size. The error term for control of the step size needs to involve both position and velocity errors and the optimum combination of these can vary with the type of simulation. The following formula has been found suitable for a range of simulations: (3)
where
At t
-
(At hl + Ath2 )
IdXA +dXmin
This formula gives additional weight to the velocity errors as the se tend to have a cumulative effect as the simulation proceeds.
A COMPARISON OF THE CALCULATION RESULTS USING FIXED TIME STEP SIZE VERSUS VARIABLE TIME STEP SIZE
This section compares the calculation results using fixed time step size and variable time step size controlled by the step doubling scheme. The efficient calculation methods (Zhang and Whiten, 1998, 1999) which generally give more accurate and efficient simulations are used for both fixed time step and variable time step simulations. This example includes calculations in normal and tangential direction and includes: sliding, sliding changing to rolling, rolling, and rolling changing to a static contact. During the initial contact sliding occurs and then this changes to rolling. The rolling velocity gradually reduces and eventually stops. The baIl stops on the slope as the rolling friction is larger than the gravity component in the tangential direction.
D. Zhang and W. J. Whiten
1344
The two objects have the following properties (input data to the program):
Object
Xx
Vx
A
2 0-60.0
10.0
B
o
Mass
radius
4
1 NIA
1e20 (infinity)
Contact parameters: stiffness (k): 2000, damping (a): 0.4, gravity (g): 9.8, sliding coefficient (f.1.): 0.7, rolling coefficients (q ): 0.4
(% ): 0.3, static friction coefficient:
10.
The simulation was run using a range of fixed time step sizes. It was found that the time step size I:1t 0.01 is suitable and gives a similar accuracy to the variable step size calculation. Larger time steps give incorrect results, smaller time steps give an unnecessary long simulation time.
=
The calculation time for the simulations are given as below: time step size time to contruct table (s) calculation time (s)
I:1t = 0.1
I:1t=O.Ol
1.8 1.5
1.8 14.3
I:1t
= 0.001 1.8 153.2
variable 2.0 11.9
=
As shown in Figures 1 and 2, the simulation with fixed time step size I:1t 0.01 and variable time step size give similar accuracy in the normal direction and use a similar time for simulation. As is described in Zhang and Whiten (1999), the tangentia1 calcu1ation is not very sensitive to the time step size.
correct 8
~,,
7
"-.
.,.~
y
x
variable
0
I:1t =0.01
•
I:1t=O.l
12
14
'\
6
5
2
4
6
8
10
16
x Fig.1 Trajectory of the simulation with variable step size and fixed step size. Both simulations correctly calculate the sliding friction and the rolling friction during the contact. When the ve10city in the tangential direction progressively reduced to zero, the static friction is applied. The particle then remained in the static state since gravity is not large enough to overcome the maximum static friction. Ali the switching points in this example such as switching from sliding mode to rolling mode and a1so from rolling mode to static mode are correctly detected and calculated properly by both methods.
Step size control for efficient discrete element simulation
1345
0.06
,1t=O.Ol
....................
0.05
variable tirne step
0.04 1-1
0
0.03
1-1 1-1
~
0.02
0.01
0 0
2
3
4
5
6
7
8
Tirne Fig.2 The error in trajectory in Figure 1 with two simulations.
CONCLUSIONS A variable time step method has been developed for the table driven velocity and position updating method of discrete element simulation and this method has been compared with the use of a fixed time step. The simulation results show: 1.
The variable time step size method al ways gives a reliable solution while the fixed time step method with too large a time step size gives inaccurate simulation results. The fixed time step method with reasonable size steps gives stable simulation results.
2.
To obtain a given accuracy with a fixed time step method, several runs with different time steps may be required.
3.
The variable time step size controlled by the step doubling scheme can find a suitable step sizes for ail parts of simulation process.
4.
The simulation using the variable time step method is faster th an the simulation with fixed step size when most of the particIes are free particIes, since many fewer steps are required.
5.
The simulation using the variable time step method required more calculation steps th an that with fixed step size in general during the multiple contact process and particIe oscillation process.
Depending on the mix of conditions, the simulations using a fixed time step size with known suitable time step size gives roughly similar simulation times. However, the suitable time step size is al ways unknown. The step doubling simulation method automatically chooses the time step according to the conditions avoiding the need for multiple runs.
ACKNOWLEDGMENTS Authors would like to thank the JKMRC management for support of this project.
1346
D. Zhang and W. J. Whiten
REFERENCES Cundall, P.A. and Strack, O.D.L, A discrete numerical model for granular assemblies, Geotechnique, 1979, 29,47-65. Press, W.H., Teulolsky, S.A., Vetterling, W.T. and Flannery, B.P., Numerical Recipes in C: The Art of Scientific Computing, 2nd edn. 1992, Cambridge University Press. Tsuji, Y., Kawaguchi, T. and Tanaka, T., Discrete particle simulation of two-dimensional fluidized bed, Powder Technology, 1993,77,79-87. Zhang, D. and Whiten, W.J., The Calculation of Contact Forces between Particles with Spring and Damping Models, Powder Technology, 1996,88,59-64. Zhang, D. and Whiten, W.J., An efficient calculation method for particle motion in discrete element simulations, Powder Technology, 1998,98,223-230. Zhang, D. and Whiten, W.J., A New Calculation Method for Particle Motion in Tangential Direction in Discrete Element Simu.1ations, Powder Technology, 1999,102,235-243.
Correspondence on papers published in Minerais Engineering is invited bye-mail to [email protected]
Minerais Engineering, Vol. 14, No. 10, pp. 1341-1346,2001 © 2001 Elsevier Science Ltd Ali rights reserved 0892-6875/01/$ - see front matter
STEP SIZE CONTROL FOR EFFICIENT DIS CRETE ELEMENT SIMULATION*
D. ZHANG and W.J. WHITEN Julius Kruttschnitt Mineral Research Centre, University of Queensland, 4072 Australia Email: [email protected] (Received 9 April 2001; accepted 31 May 2001)
ABSTRACT
The step size determines the accuracy of a discrete element simulation. The position and velocity updating calculation uses a pre-calculated table and hence the control of step size can not use the integration formulas for step size control. A step size control scheme for use with the table driven velocity and position calculation uses the difference between the calculation result from one big step and that from two smalt steps. This variable time step size method chooses the suitable time step size for each particle at each step automaticalty according to the conditions. Simulation using fixed time step method is compared with that of using variable time step method. The difference in computation time for the same accuracy using a variable step size (compared to the fixed step) depends on the particular problem. For a simple test case the times are roughly similar. However, the variable step size gives the required accuracy on the first run. A fixed step size may require several runs to check the simulation accuracy or a conservative step size that results in longer run times. © 2001 Elsevier Science Ltd. Alt rights reserved. Keywords Modelling; simulation
INTRODUCTION The discrete element simulations model the behaviour of fluid and solids by using an assembly of discrete elements. This method can be widely applied in mining and mineraI processing areas. The discrete element models divides the simulation time into small time steps. For each time step the velocity and position of each particle are updated using the equation of motion and the forces acting on the particle. The size of the time step is critical to the calculations since it determines the stability and accuracy of the calculation in the discrete element model. A time step size which is too small results in an unnecessarily long simulation time. However, a time step size which is too large produces incorrect simulation results.
• Presented at Comminution '01, Brisbane, Australia, March 2001
1341
1342
D. Zhang and W. J. Whiten
According to Cundall and Strack (1979), the following relation must be satisfied to give stable calculation:
I1t < 2.Jm/ k,
(1)
k m
where
11 t
stiffness coefficient, particle mass, is size of time step.
Tsuji (1993) found that the calculation became unstable when the time step near the limit value given by this equation. Therefore, Tsuji proposed that the time step should be proportional to the oscillation period of the spring-mass system. The oscillation period of the system is given by
I1t < 27r.Jm / k
(2)
Tsuji suggested a suitable time step can be obtained by dividing one half of the natural oscillation period by a small integer n, eg 5. Both Tsuji's and Cundall's time step equation have a similar form with Tsuji recommending a slightly smaller value. However, the particle separation time is less than one half of the oscillation period, because the impacting particles will separate when impact force equals to zero (Zhang and Whiten, 1996), and a second collision or changes in other forces could happen within this time step to affect the particles impact. A smaller time step may be needed to obtain accurate calculation of the particle motion. Often it is necessary to check the simulation results using different step sizes. Zhang & Whiten (1998) introduced a method of using accurately integrated (Press et al., 1992) tables of updates to the velocity and position of each discrete element. A step doubling scheme to control the accuracy of this velocity and position updating by varying the step size has been developed. This method allows the time step size to be adjusted to a suitable size automatically according to the calculation error in each part of the simulation. The step doubling scheme calculates each simulation step twice: firstly, it calculates the simulation results for the whole step, and secondly, it calculates the results by using two half steps. The calculation error for each step is defined as the difference of the simulation results from one full step and the results from two half steps. If the calculation error is within a predefined error tolerance, the simulation moves on to the next step. The step size used for the next step is related to the calculation error of the previous step. If the calculation error is larger than the error tolerance, the step size is reduced and the calculations repeated until the error is within the error tolerance. The step size needed to give the required accuracy typically varies in different parts of a discrete element simulation and at different times in the simulation. The extra computation needed to check the accuracy of a given step size can often be justified by the reduction in computation obtained by being able to use larger steps in parts of the calculation. In addition, the controlled accuracy avoids the need to choose a conservative step size and to make multiple runs to check the accuracy. Both of these require additional calculation.
STEP DOUBLING SCHEME USED IN DISCRETE ELEMENT SIMULATIONS The step doubling scheme used in discrete element simulations is as follows: 1.
Predefine the maximum acceptable error (error tolerance), maximum time step size and minimum time step size. 2.· Calculate the velocity and displacement of every object using the time step determined from the previous step. 3. Calculate the velocity and displacement of every object using two half time steps. 4. Calculate the difference of the two results and compare with the predefined maximum error; for the particles which have the error larger than the predefined error, repeat steps 2 and 3 using half of the old time step as the full time step.
Step size control for efficient discrete element simulation
5. 6.
1343
Calculate a more accurate value for the motion by combining the simulation results from two half time steps. Determine the time step size for the next time step from the calculation error. If the error is much less than the predefined error, a larger time step size will be used for next time step.
This scheme can be applied to any method of updating the velocity and displacement of the simulation objects but is particularly suited to the table driven method developed in Zhang and Whiten (1998), where the updates are not made using a numerical formula for the solution of differential equations. The repeat of calculation steps 2 & 3 from step 4 needs to be done only for the discrete elements that do not satisfy the accuracy test. For these particles multiple steps are calculated until all particles in the simulation are updated to the same time. The time step size is constructed by a time step base times a time step factor which is determined by the overall conditions of the contact. The time step base remains the same for the whole simulation. Local time step factors determined by the calculation error are applied to determine the time step size of each element. For the particles that are not in contact, the difference between one full step and two half steps is very small. So a larger time step size can be applied. Therefore the step size control scheme can avoid unnecessary use of small time step sizes. For the particles that are in contact, a smaller time step size must be used as the calculation error becomes larger especially when the particles change from non-contact to contact and contact to non-contact. For practical usage, the time step sizes used in the step doubling are limited between predefined maximum time step size and minimum time step size. The error term for control of the step size needs to involve both position and velocity errors and the optimum combination of these can vary with the type of simulation. The following formula has been found suitable for a range of simulations: (3)
where
At t
-
(At hl + Ath2 )
IdXA +dXmin
This formula gives additional weight to the velocity errors as the se tend to have a cumulative effect as the simulation proceeds.
A COMPARISON OF THE CALCULATION RESULTS USING FIXED TIME STEP SIZE VERSUS VARIABLE TIME STEP SIZE
This section compares the calculation results using fixed time step size and variable time step size controlled by the step doubling scheme. The efficient calculation methods (Zhang and Whiten, 1998, 1999) which generally give more accurate and efficient simulations are used for both fixed time step and variable time step simulations. This example includes calculations in normal and tangential direction and includes: sliding, sliding changing to rolling, rolling, and rolling changing to a static contact. During the initial contact sliding occurs and then this changes to rolling. The rolling velocity gradually reduces and eventually stops. The baIl stops on the slope as the rolling friction is larger than the gravity component in the tangential direction.
D. Zhang and W. J. Whiten
1344
The two objects have the following properties (input data to the program):
Object
Xx
Vx
A
2 0-60.0
10.0
B
o
Mass
radius
4
1 NIA
1e20 (infinity)
Contact parameters: stiffness (k): 2000, damping (a): 0.4, gravity (g): 9.8, sliding coefficient (f.1.): 0.7, rolling coefficients (q ): 0.4
(% ): 0.3, static friction coefficient:
10.
The simulation was run using a range of fixed time step sizes. It was found that the time step size I:1t 0.01 is suitable and gives a similar accuracy to the variable step size calculation. Larger time steps give incorrect results, smaller time steps give an unnecessary long simulation time.
=
The calculation time for the simulations are given as below: time step size time to contruct table (s) calculation time (s)
I:1t = 0.1
I:1t=O.Ol
1.8 1.5
1.8 14.3
I:1t
= 0.001 1.8 153.2
variable 2.0 11.9
=
As shown in Figures 1 and 2, the simulation with fixed time step size I:1t 0.01 and variable time step size give similar accuracy in the normal direction and use a similar time for simulation. As is described in Zhang and Whiten (1999), the tangentia1 calcu1ation is not very sensitive to the time step size.
correct 8
~,,
7
"-.
.,.~
y
x
variable
0
I:1t =0.01
•
I:1t=O.l
12
14
'\
6
5
2
4
6
8
10
16
x Fig.1 Trajectory of the simulation with variable step size and fixed step size. Both simulations correctly calculate the sliding friction and the rolling friction during the contact. When the ve10city in the tangential direction progressively reduced to zero, the static friction is applied. The particle then remained in the static state since gravity is not large enough to overcome the maximum static friction. Ali the switching points in this example such as switching from sliding mode to rolling mode and a1so from rolling mode to static mode are correctly detected and calculated properly by both methods.
Step size control for efficient discrete element simulation
1345
0.06
,1t=O.Ol
....................
0.05
variable tirne step
0.04 1-1
0
0.03
1-1 1-1
~
0.02
0.01
0 0
2
3
4
5
6
7
8
Tirne Fig.2 The error in trajectory in Figure 1 with two simulations.
CONCLUSIONS A variable time step method has been developed for the table driven velocity and position updating method of discrete element simulation and this method has been compared with the use of a fixed time step. The simulation results show: 1.
The variable time step size method al ways gives a reliable solution while the fixed time step method with too large a time step size gives inaccurate simulation results. The fixed time step method with reasonable size steps gives stable simulation results.
2.
To obtain a given accuracy with a fixed time step method, several runs with different time steps may be required.
3.
The variable time step size controlled by the step doubling scheme can find a suitable step sizes for ail parts of simulation process.
4.
The simulation using the variable time step method is faster th an the simulation with fixed step size when most of the particIes are free particIes, since many fewer steps are required.
5.
The simulation using the variable time step method required more calculation steps th an that with fixed step size in general during the multiple contact process and particIe oscillation process.
Depending on the mix of conditions, the simulations using a fixed time step size with known suitable time step size gives roughly similar simulation times. However, the suitable time step size is al ways unknown. The step doubling simulation method automatically chooses the time step according to the conditions avoiding the need for multiple runs.
ACKNOWLEDGMENTS Authors would like to thank the JKMRC management for support of this project.
1346
D. Zhang and W. J. Whiten
REFERENCES Cundall, P.A. and Strack, O.D.L, A discrete numerical model for granular assemblies, Geotechnique, 1979, 29,47-65. Press, W.H., Teulolsky, S.A., Vetterling, W.T. and Flannery, B.P., Numerical Recipes in C: The Art of Scientific Computing, 2nd edn. 1992, Cambridge University Press. Tsuji, Y., Kawaguchi, T. and Tanaka, T., Discrete particle simulation of two-dimensional fluidized bed, Powder Technology, 1993,77,79-87. Zhang, D. and Whiten, W.J., The Calculation of Contact Forces between Particles with Spring and Damping Models, Powder Technology, 1996,88,59-64. Zhang, D. and Whiten, W.J., An efficient calculation method for particle motion in discrete element simulations, Powder Technology, 1998,98,223-230. Zhang, D. and Whiten, W.J., A New Calculation Method for Particle Motion in Tangential Direction in Discrete Element Simu.1ations, Powder Technology, 1999,102,235-243.
Correspondence on papers published in Minerais Engineering is invited bye-mail to [email protected]