On the optimal numerical time integration for DEM using Hertzian force models

Computers and Chemical Engineering 58 (2013) 211–222

Contents lists available at ScienceDirect

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

On the optimal numerical time integration for DEM using Hertzian force models Matthew Danby a,1 , John Shrimpton a,∗ , Mark Palmer b a b

Energy Technology Research Group, School of Engineering Sciences, Highfield Campus, University of Southampton, Southampton SO17 1BJ, UK Inhaled Product Development, GlaxoSmithKline R&D, Park Road, Ware, Hertfordshire SG12 0DP, UK

a r t i c l e

i n f o

Article history: Received 25 May 2011 Received in revised form 13 June 2013 Accepted 27 June 2013 Available online 29 July 2013 Keywords: Discrete element method Numerical integration Integration schemes

a b s t r a c t The numerical accuracy of a selection of different time integration techniques used to solve particle motion is investigated using a normal collision employing the non-linear Hertzian contact force. The findings are compared against the linear force model where it has been found that the expected order of accuracy of higher-order integration schemes is not realised (Tuley et al., 2010). The proposed mechanism for this limitation has been cited as the errors in integration which occur across the force profile discontinuity. By investigating the characteristics of both the non-linear elastic and the non-linear damped Hertzian contact models, it has been found that higher orders of accuracy are recoverable and depends on the degree of the governing non-linear equation. The numerical errors of the linear and non-linear force models are however markedly different in character. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction This contribution focuses on the relative scheme performance of a range of algorithms to integrate the equations of motion as applied to a binary particle collision. By examining the behaviour of ten different algorithms, both in terms of the differences in performance of the integrator characteristics (e.g. multi-point versus single-point), nominal scheme order (e.g. first-, second- and fourthorder) and how these characteristics are affected by the underlying collision model (e.g. linear versus non-linear) a greater insight as to which numerical integration scheme(s) to use in larger-scale Discrete Element Method (DEM) simulations may be gained. The schemes chosen represent a cross-section of implicit and explicit schemes, all of which have their main characteristics summarised in Table 2 with the full numerical procedure being listed elsewhere (Tuley, Danby, Shrimpton, & Palmer, 2010). A common feature of all DEM approaches is the discretisation of time into distinct time-points, which may be spaced regularly in time, ‘fixed timestep’ or non-uniformly, ‘adaptive timestep’. In order to advance time from the current time point to the next, by considering the resultant DEM force only, the particle trajectories are calculated then integrated using the equations of motion over the current timestep. The scheme used to integrate time has

∗ Corresponding author. E-mail addresses: [email protected], [email protected] (J. Shrimpton). 1 Partially supported by GlaxoSmithKline plc. 0098-1354/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compchemeng.2013.06.018

been shown to greatly influence the accuracy of the resulting trajectory (Tuley et al., 2010), with factors such as the magnitude of the force profile discontinuity at impact, the placement of the time points about the instant of particle collision and the degree of collision inelasticity all affecting the resultant performance. In addition, numerical schemes which have an expected lower order of accuracy may outperform those which have a higher expected order of accuracy. However, the degree by which the Hertzian DEM collision force model influences the integrator scheme performance has as yet been unaddressed. The remainder of the paper is organised as follows: Section 2 provides details of the problem configuration and a general consideration of collision force models. Section 3 gives details of how the numerical performance of these models is parameterized. Here also a brief discussion of the analytical solution is given, along with further details in Appendix A. Section 4 gives an overview of the numerical time integration schemes tested. Section 5 discusses the results and Section 6 lists the main conclusions of the study. 2. Test configuration 2.1. Definition A normal collision between two particles is used to asses the numerical scheme accuracy. The test comprises of a stationary particle q and a moving particle p as given in Fig. 1. Particle q, which is initially located at the co-ordinate system origin, is given an initial velocity towards q such as to impact normally with q which remains at rest throughout the test. The position of the stationary

212

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

test. The physical parameters employed for the numerical tests are presented in Table 1.

List of symbols a cˇ t e Fe , F k m NFM r t tcollision ti u x xpq

collision resolution damping coefficient time step collision coefficient of restitution contact force, elastic and viscous contact stiffness particle mass free motion timestep fraction particle radius time collision duration time at impact velocity position particle centre separation

2.2. The collision force model Interparticulate force models can be grouped into the following four categories (Kruggel-Emden, Simsek, Rickelt, Wirtz, & Scherer, 2007): continuous potential, linear, non-linear and hysteretic models. The Hertzian contact force is a non-linear model. For the numerical testing described here, the approach given by Tuley et al. (2010) is followed, where the linear force model has been investigated. The non-linear force model describes collision behaviour using an elastic spring and viscous dissipation force. According to Newton’s second law, the total force applied to particle p is equal to its rate of rate of change of linear momentum. The equation of motion for particle p can be expressed in terms of its mass m, time t, velocity u and the two force components acting during the colliding state Fe and F (elastic and viscous damping force components, respectively).

Greek letters ˛ degree of contact stiffness nonlinearity ˇ degree of contact damping nonlinearity particle overlap ıpq εxamp position amplitude error particle density 

dxp dt m

Subscripts n integration time point p, q particles p and q

xpq = xq − xp



ıpq =





rp + rq − xpq  , 0,

Fe = −k˛ ıpq

stationary particle q

Fig. 1. Diagram of test case particles p and q.

(2)

  xpq  < rp + rq   xpq  ≥rp + rq

(3)

The elastic repulsive force arises from the magnitude of the overlap between the particle surfaces, ıpq and the normal contact stiffness, generally k˛ . The exact dependance of this force on particle overlap depends on which type of model is employed, the general case is:

particle q is chosen to give a specific number of timesteps prior to the moment of impact, hence controlling the value of NFM (see Section 3.3). The test is run until particle p rebounds from q and returns back to its starting point at the origin. Using the DEM algorithm to resolve the particle collision, the position and velocity of particle p is solved by numerical integration using one of the schemes under

x

(1)

The DEM collision forces only operate when the two surfaces of the impacting particles are overlapping. Given the radii of particles p and q, rp and  rq , respectively, and the relative position of their centres, xpq , the magnitude of overlap is given as:

Abbreviations AB2 second order Adams–Bashforth fourth order Adams–Bashforth AB4 ABM fourth order Adams–Bashforth–Moulton first order back-differenced BD1 BD2 second order back-differenced fourth order back-differenced BD4 DEM discrete element method FE forward Euler RK2 second order Runge–Kutta fourth order Runge–Kutta RK4 SE symplectic Euler

moving particle p

⎧   ⎨ Fe + F , xpq  < rp + rq =   ⎩0, xpq  ≥rp + rq

The relative position between particle centres p and q is:

Superscripts non-dimensionalised *

velocity, u

dup dt

= up

˛

(4)

In the case of applying Hooke’s law, ˛ is unity, thus being a linear elastic force. The linear elastic force has proven to be a popular choice in the literature since the seminal publication of Cundall and Strack (1979). However, by considering the theoretical pressure distribution acting upon two conforming curved surfaces, Hertz was the first to provide a description of the forces acting upon two ideal elastic non-cohesive spheres (Hertz, Jones, & Schott, 1896). If the centres of two conforming spheres are brought together, integration of the resulting normal pressure distribution from the surface deformation results in an elastic repulsive force which is related to the surface overlap non-linearly, in this case ˛ = 3/2 is well-established and was first employed by Tsuji, Tanaka, and Ishida (1992) for practical DEM simulation. For ‘real’ materials the degree of non-linearity may be physically interpreted as the response of the material properties (which may soften or harden with compression) of the contact geometry itself. Although DEM methods are capable of simulating the interaction of particles of complex geometries (Schinner, 1999), the computational complexities of simulating large-scale particle interaction is greatly simplified by using spherical particles since the task of collision detection and interparticle force calculation is

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

simpler. Therefore, it is of practical relevance to focus the investigation towards Hertzian contacts given its wide-ranging applications to the simulation of spherical granular media. Despite being well established for spherical contacts, Hertzian theory does not account for the energy losses which occur during real particle impacts. A viscous damping term is required to account for the energy dissipation occurring during an inelastic collision, hence an additional term is required. The most common general form of the viscous damping being:



F = −cˇ ıpq

ˇ

dxpq dt



(5)

where cˇ is the damping coefficient. For the case where ˇ = 0 and  = 1, the linear energy dissipation model, or ‘dashpot damping’ is recovered. Linear damping has been widely noted for its nonphysical behaviour in the literature (Zhang & Whiten, 1996) since it produces a non-zero viscous force despite a vanishing overlap and a spurious attractive force prior to separation. To overcome the problem of a damping term which should increase from zero, it is made a function of displacement (which does increase from zero) by choosing a positive ˇ term. The literature reveals that the choice of ˇ depends on the type of DEM simulation, with some authors suggesting ˇ = 1/2 being theoretically and experimentally justified for sphere–sphere contacts (Kuwabara & Kono, 1987), and ˇ = 1/4 being better suited for sphere–wall collisions (Falcon, Laroche, Fauve, & Coste, 1998). Despite including this term however, the spurious reversal of the viscous force component during particle separation remains. Of the literature reviewed, no authors have proposed choosing a viscous damping force using anything other than as a linear function of impact velocity, i.e.  = 1. This parameter could be a topic of investigation, since if it were made negative, then the damping term could be attenuated at the point of impact and separation, reducing the non-physical effects of force reversal nearing separation. To summarise, two force-model cases are are examined here: the purely Hertzian contact (i.e. ˛ = 3/2, cˇ = 0) and the non-linear / 0, ˇ = 1/2). For brevity, theses are referred damped (i.e. ˛ = 3/2, cˇ = to as elastic and damped cases, respectively. 3. The test model For the discussion which follows, the subscripts on xp and upq are dropped to simplify the notation, since particle q remains stationary throughout. The terms x and u hereinafter refer to the position and velocity of particle p. The force applied to particle p is equal to its rate of change of linear momentum. This can be expressed as an equation of motion in terms of time t, mass m, velocity up and the elastic (Fe ) and viscous damping (F ) force components, respectively.

213

Fig. 2. NFM is defined as the fraction of timestep t occurring before the nearest integration time point to the moment of impact.

Since the error is defined relative to the maximum analytical particle overlap, the initial separation of particles p and q does not affect the error measurement therefore results from different NFM are directly comparable. Details of how the analytical maximum overlap is calculated can be found in Appendix A. As previously reported by Tuley et al. (2010), the error analysis is conducted on the position profile since it has a clearly defined maxima. Since the numerical discretisation is not guaranteed to coincide with the peak displacement, a second-degree polynomial is constructed using the three integration time points closest to the peak to obtain a better estimate of the maxima at the point of inflection. 3.2. Analytical solution Initially, particle p is in free motion with velocity u = 4 m s−1 directed towards the stationary particle q; no forces act and therefore the position linearly increases. Upon contact with q, the elastic and viscous damping force terms increase from zero to repel particle p. The collision ends when the particles velocity has completely reversed and p returns to a state of free motion. The solution to the equation of motion can therefore be decomposed into two parts: free motion and colliding states. The trivial solution to the free motion part holds regardless of the collision force model. The complete analytic solution to the colliding state is dependent on the force model; for the linear collision model there exists an exact analytical solution for the particle motion by application of the under-damped spring-mass-damper harmonic oscillator (Kreyszig, 2006), however this is not so for the non-linear force models. Instead, for the non-linear case with damping the approach given by Ramirez, Poschel, Brilliantov, and Schwager (1999) and Schwager and Poschel (1998) is followed, which consists of a series expansion solution. For the undamped non-linear elastic case (i.e. a purely Hertzian collision) an analytical solution to determine the maximum interparticle overlap and collision duration can be used to determine the scheme error (Hertz et al., 1896). The full details of these solutions are given in Appendix A. During free motion, the time to impact ti is given by: ti =

xpq,0 − (rp + rq ) u0

(7)

and the position prior to impact is simply:

3.1. Integration scheme error measurement

x (t) = u(0)t

The approach taken to quantify the error of the numerical integration schemes is conducted by considering a single normal particle collision between particles p and q as shown in Fig. 1. After the trajectory is allowed to evolve such that the mobile particle has returned back to its origin, the resultant position profile is used to assess the error. By comparison of the absolute maxima in the numerical position profile ( xnumerical ) to that from the analytical solution ( xanalytical ) and normalising by the analytical maximum

The timestep t is non-dimensionalised (t ) as a fraction of the analytical collision duration, tcollision . Here this is defined as the time between contact and separation. The details of how tcollision is determined can also be found in Appendix A.

x particle overlap (ı

pq analytical ) the error εamp is thus defined as:    xnumerical − xanalytical 

εxamp =

ı

pq analytical

(6)

(t < ti )

(8)

3.3. Collision boundary discretisation During free motion, there are no net forces acting on the mobile particle and all of the numerical integration methods give exact solutions. The free motion is represented by NFM and is the fraction of the timestep existing immediately prior to the initial point of interparticle contact (see Fig. 2). The placement of the integration time points around the start of a collision affects the accuracy of the

214

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

Table 1 Default physical parameters of the standard (Fig. 1) test cases. Property

Particle p

Particle q

Radius, r (m) Density,  (kg m−3 ) Initial velocity, u(0) (m s−1 ) Contact stiffness, k (N m−3/2 ) Coefficient of restitution, e

0.01 1000 4.0

0.1 ∞ 0.0 105 1 (elastic) 0.8 (damped)

numerical solution. In the test setup here, the non-dimensional NFM represents the fraction of timestep which occurs before the start of the collision, i.e.: NFM =

ti − tnearest t

(9)

where tnearest is the integration time point closest in time to the moment of impact and ti is exact time of impact found analytically from Eq. (7), hence −0.5 ≤ NFM ≤ 0.5. This is illustrated in Fig. 2, in this case tnearest would be defined as t2 . For NFM = 0 the force, velocity and position are evaluated precisely at the start of the collision, NFM = +0.5 and NFM = −0.5 are equivalent, and time of impact occurs exactly halfway between evaluation time points. 4. Numerical integration 4.1. The schemes evaluated Ten time integration algorithms are examined in the context of a single inter-particle collision. The rationale for their choice is based upon comparison of the performance of: nominal scheme order (i.e. nominally first, second and fourth), stability (i.e. implicit and explicit) and expense (i.e. single-point and multi-point). A summary of the schemes investigated together with these properties is listed in Table 2. The schemes tested here are not an exhaustive list; there are other time integration schemes which have been used in the literature too (see Dziugys & Peters, 2001; Kruggel-Emden, Sturm, Wirtz, & Scherer, 2008). For example, the multi-point third-order Verlet method (Aoki & Akiyama, 1995; Kopf, Paul, & Dunweg, 1997; Satoh, 1995a, 1995b) is widely used in molecular dynamics (Haile, 1992). A limitation of the Verlet algorithm is that it is not well suited to collisions involving a velocity dependant force, since velocities and positions are treated differently. However, modifications such as the second-order half-step leap-frog Verlet has proven a popular choice by DEM researchers (Langston & Tuzun, 1994; Langston,

Tuzun, & Heyes, 1995a, 1995b; Lan & Rosato, 1997) which staggers the evaluation of the position and velocity by half a timestep and can be used with velocity dependant force models. A discussion of integration methods for DEM would not be complete without mention of the popular Gear’s predictor–corrector method (Haile, 1992). This predictor stage comprises of a Taylor series expansion of the position and higher order time derivatives. By evaluating the discrepancy in particle acceleration (as calculated from the force derived from the predicted position and velocity) against the acceleration calculated in the prediction step itself, a correction stage is applied which revises the position and higher time derivatives using appropriate Gear coefficients (Rougier, Munjiza, & John, 2004) depending on the scheme order. A study of the literature reveals that explicit time integration schemes are favoured. The use of implicit methods such as those we have listed in Table 2 have not been widely reported in the literature, two examples from the literature are cited in the discussion. A comprehensive discussion of the characteristics of the ten schemes employed here has been reported previously (Tuley et al., 2010). 5. Results The numerical results presented here allows a direct comparison between the integration schemes and an assessment of the effect of varying the contact force parameters, in this case the coefficient of restitution. Dziugys and Peters (2001) present the main requirements by which time integration methods should comply. These include scheme stability, low computational time and storage and conservation of energy and momentum. In addition to these prerequisites, the integration scheme should also be robust enough to give repeatable DEM results. In order to determine how robust the schemes are, the sensitivity to the fraction of free motion timestep, NFM , prior to the instance of collision is examined. This test allows investigating of how the initial conditions and subsequent interparticle collision events in a practical DEM simulation affect the resultant numerical accuracy. The realised actual scheme order is compared to the nominal scheme order, which has been found in some cases to be lower. Since we find that the integration scheme accuracy is a function of the pre-collision free motion (to varying degrees, depending on the scheme), results for the two extreme cases, i.e. NFM = 0 and NFM = 0.5 are presented. The effect of collision inelasticity by means of increasing the interparticulate coefficient of restitution is also investigated. It is well recognised that the linear dashpot model for collision energy

Table 2 Summary of integration schemes tested. Key: Forward Euler (FE), Symplectic Euler (SE), Second order Runge–Kutta (RK2), Fourth order Runge–Kutta (RK4), Second order Adams–Bashforth (AB2), Fourth order Adams–Bashforth (AB4), Fourth order Adams–Bashforth–Moulton (ABM), First order back-differenced (BD1), Second order backdifferenced (BD2), Fourth order back-differenced (BD4), Explicit (E), Implicit (I), Number of iterations (Nh ). Scheme

Nominal order

Type

DEM force evaluation (per timestep)

Memory required (preceding time-point data required)

Existing literature usage

FE

1

E

1

1

Taguchi (1992) Sadd, Tai, and Shukla (1993) Tsuji, Kawaguchi, and Tanaka (1993) Tsuji et al. (1992)

SE RK2 RK4

1 2 4

E E E

1 2 4

1 1 1

AB2 AB4 ABM BD1 BD2 BD4

2 4 4 1 2 4

E E E I I I

1 1 2 Nh Nh Nh

2 4 3 1 2 4

Ovesen, Petersen, and Perram (1996) Shida, Suzuki, and Kawai (1997) Sundaram and Collins (1996)

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

loss produces non-physical forces at the point of interparticulate contact and release. This side-effect should be improved with the non-linear damping term (see Section 2.2). It is also recognised that the coefficient of restitution is not solely an intrinsic physical property, but also has a dependency on the impact velocity in the non-linear case (Schwager & Poschel, 2007). However, by reducing the coefficient of restitution for the collision, the sensitivity of the schemes accuracy as compared to the linear damping model can be determined. Unlike the linear elastic force model examined by Tuley et al. (2010), for the purely Hertzian force model (i.e. the undamped non-linear force model with ˛ = 3/2) the collision duration (tcollision ) does depend on the relative particle impact velocity (for ˛ = / 1, see Eq. (17)). Therefore the timestep resolution of the collision is affected by the collision velocity. For the purely Hertzian contact, if the impact velocity is increased by factor m, the collision duration decreases by m−1/5 . So rather than investigating a range of impact velocities, its effect on numerical accuracy is the same as adjusting the contact resolution by the appropriate factor. In addition to reducing the collision duration for Hertzian contacts, increasing impact velocities also increase the maximum interparticle overlap distance. In the linear elastic case, doubling the impact velocity would double the maximum particle overlap however, in the purely Hertzian contact, this consequence is less severe; by increasing the impact velocity by a factor m, the maximum interparticle overlap is increased by m4/5 . This is an advantageous result, since large particle overlaps are generally regarded as undesirable and unrealistic: for example Cleary and Sawley (2002) state that ideal overlap distances are in the range of 0.1–1.0% of particle diameter, and Matuttis, Luding, and Herrmann (2000) state that particle overlap should be “much smaller than the particle radius”. However, some authors allow much greater particle–particle overlap distances – Hogue and Newland (1994) permit overlap up to the smallest particle radius. Although ‘realistic’ or ‘ideal’ maximum particle–particle overlap may be in the range of 0.1–1.0%, there is no numerical reason to constrain the overlap to such limits. The only numerical constraint is that problems will be encountered if the inter-particle overlap exceeds the sum of the two particle radii, due to the instantaneous switch of normal direction between particle centre-points. In the results to follow, Section 5.1 examines the error in the damped and undamped non-linear force models with respect to the proximity of the time step boundary to the collision time point. Section 5.2 examines the actual versus the expected order of the scheme, and in terms of the degree of non-linearity of the force model. Section 5.3 examines this point further, in terms of the elastic force model order, and Section 5.4 examines the effect of damping. 5.1. Absolute error analysis: collision boundary discretisation In this set of tests, the undamped non-linear and damped / 0, non-linear force models (i.e. ˛ = 3/2, cˇ = 0 and ˛ = 3/2, cˇ = ˇ = 1/2, respectively), are examined, herein referred to as elastic and damped for brevity. The initial separation distance between the two spheres p and q is chosen to produce a specific amount of ‘free motion’ prior to impact. This is to compare the elastic and damped force models sensitivity to the timestep placement about the commencement of collision. Figs. 4–6 show that at collision resolution of about 300 timesteps per collision, all the integration schemes tested have uniform order of accuracy for diminishing timestep, hence the non-dimensional timestep, t for this set of tests is 3 × 10−3 . The results presented in Fig. 3 show that this temporal discretisation about the instant of impact may yield considerably different accuracy in the numerical solution, depending on the integration scheme used and whether or not damping is present.

215

For the nominally first-order integration schemes, it is clear that there is no notable dependance on the NFM on the resultant error for the FE, BD1 or SE schemes regardless of whether the collision is damped or elastic (see Section 3.1). The only notable observation is that the SE scheme has consistently better accuracy (the error is lower in both undamped and damped cases), and while the FE and BD1 schemes are insensitive to whether or not damping is present, the SE scheme loses two orders of magnitude in accuracy if the collision is inelastic. This contrasts startlingly with the findings in Tuley et al. (2010), where the SE scheme was the only nominally first-order scheme found to sensitive to the choice of NFM , spanning six orders of magnitude. Therefore the SE scheme behaves more robustly with the non-linear force model, so for repeatable and consistent DEM results, this reliability is likely to be beneficial for practical calculations including multiple particle collisions under the full spectrum of NFM . From Fig. 3(c), just as in the case of the nominally first-order schemes, there is no marked influence of NFM on the accuracy of the AB2, BD2 or RK2 integration schemes in the undamped case; AB2 and BD2 schemes perform similarly with single-point RK2 achieving greater accuracy by an order of magnitude. Contrast this finding with the damped case in Fig. 3(d) and it is found that all the nominally second-order schemes investigated are sensitive to the fraction of ‘free-motion’ prior to collision. The dependance on error of these schemes to NFM appears disorderly. Whereas the BD2 scheme appears to favour NFM ≈ 0 for best accuracy, the AB2 scheme favours NFM ≈ 0.5. With damping, the RK2 scheme loses its accuracy advantage over the other two schemes and the resultant error has a complex dependance on NFM also. Again, these results are in contrast to the findings of the linear force model (Tuley et al., 2010), where the scheme accuracy showed symmetry about NFM = 0 and maxima and minima were located at either NFM = 0 or NFM = 0.5 for all schemes investigated. The sensitivity in error of the nominally fourth-order schemes with NFM is greatest of all the different order schemes investigated. For the elastic case, where lower-order schemes showed insensitivity to NFM , the nominally fourth-order schemes display a more complex dependance. The nominally fourth order scheme RK4 demonstrates the greatest variation in accuracy with NFM spanning three orders of magnitude (see Fig. 3(e)), although this is less than the equivalent case studied in Tuley et al. (2010). For all NFM , the damped case has lower accuracy for all nominally fourth-order schemes than the undamped case. A trend found with the linear force model is that the AB4, ABM and BD4 schemes perform almost identically with varying NFM ; the same trend is true here also, with only slight discrepancies in the elastic case. It can be concluded that the non-linear force model affects the scheme error with NFM quite differently than in the linear case; for nominally first-order schemes there is almost no dependance and likewise for undamped nominally second-order schemes. The remaining schemes have an unpredictable behaviour, with maxima and minima in error located at fractions of NFM rather than in the linear force model case, where best- and worst-case errors were generally located at NFM = 0 or NFM = 0.5 and all errors being symmetric about NFM = 0. The variation in accuracy of the schemes with the discretisation around the start of the collision implies care must be taken when comparing the accuracy against other parameters (i.e. the NFM should remain constant). The effects of introducing damping are generally either to make the error more sensitive to NFM (in the nominally second-order schemes) or to reduce the numerical accuracy across all NFM (in the nominally fourth-order schemes). For geometrically simple particles such as the spherical ones employed here, at the moment a collision event is detected, it is possible to find the exact moment which they made contact (Xu & Yu, 1997) by tracking back in time. It is therefore possible to ensure the case that NFM = 0 by using a variable timestep integration method.

216

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222 0

0

10

10

−1

−1

εxamp (%)

10

−2

10

x

εamp (%)

10

−3

−3

10

10

−4

10 −0.5

−2

10

−4

0

10 −0.5

0.5

N

(b) Nominally 1st order schemes (damped)

−2

−2

εxamp (%)

10

−3

10

x

εamp (%)

10

−4

−3

10

−4

0

10 −0.5

0.5

NFM

−2

−2

10

−4

−4

10

εamp (%) x

10

x

0.5

(d) Nominally 2nd order schemes (damped)

10

εamp (%)

0

NFM

(c) Nominally 2nd order schemes (undamped)

−6

10

−8

−6

10

−8

10

−0.5

0.5

FM

(a) Nominally 1st order schemes (undamped)

10 −0.5

0

N

FM

10

0

0.5

NFM

(e) Nominally 4th order schemes (undamped)

−0.5

0

0.5

NFM

(f) Nominally 4th order schemes (damped)

Fig. 3. (a–f) Integration schemes: FE ( ), SE ( ), BD1 ( ), RK2 ( ), AB2 ( ), BD2 ( ), RK4 ( ), AB4 ( ), ABM ( ), BD4 ( ), position amplitude error dependence on the duration (in timesteps) of particle free motion before the collision starts. The undamped case is a purely Hertzian collision, the damped case has a coefficient of restitution of 0.8.

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

Resolution, a 3

2

10

2

10

10

Resolution, a 1

10

3

2

10

2

10

10

217

Resolution, a 1

10

3

2

10

2

10

10

Resolution, a 1

10

1

0

2

10

1

10

1

10

0

10

3

10

2

10 10

10

0

0

10 −2

10

−2

10

10

−1

−1

10 −4

−4

10

10

−2

10

−2

10

10

−3

−6

−3

10

−6

10

10

10

−4

−4

10

−3

10

−2

10

Δt∗

−1

−4

10

10

−3

10

−2

10

Δt∗

−1

10

10 −4 10

−4

−3

10

−2

10

Δt∗

−1

10

10 −4 10

(a) Elastic, NF M = 0 (b) Elastic, NF M = 0.5 (c) Damped, NF M = 0 Fig. 4. (a–d) Nominally 1st order schemes: FE (

), SE (

Resolution, a 3

2

10

2

10

10

) and (

3

0

2

10

2

10

10

0

10

10

−2

3

2

10

2

10

−1

10

(d) Damped, NF M = 0.5

Resolution, a 1

10

−2

10

Δt∗

) to guide the eye. The ordinate for these figures is εxamp (%).

) with first order gradient (

Resolution, a 1

10

−3

10

10

Resolution, a 1

10

3

0

2

10

2

10

10

1

10

0

10

10

−2

10

10

−2

−2

10 −4

10

−4

10

10

−4

−4

10 −6

10

−6

10

10

−6

−4

10

−3

10

−2

10

Δt∗

−1

10

−4

10

−3

10

−2

10

Δt∗

−1

10

10 −4 10

−6

−3

10

−2

10

Δt∗

−1

10

10 −4 10

(a) Elastic, NF M = 0 (b) Elastic, NF M = 0.5 (c) Damped, NF M = 0 Fig. 5. (a–d) Nominally 2nd order schemes: RK2 ( εxamp (%).

), AB2 (

Resolution, a 3

2

10

10

) and (

) with second order gradient (

Resolution, a 1

3

10

2

10

10

3

2

10

10

−2

Δt∗

10

−1

10

(d) Damped, NF M = 0.5

) to guide the eye. The ordinate for these figures is

Resolution, a 1

10

−3

10

Resolution, a 1

3

10

2

10

10

1

10

0

0

10

0

10

10

0

10

−5

−5

10

−5

10

10

−5

10

−10

−10

−10

10

10

10

−10

10

−4

10

−3

10

−2

10

Δt∗

−1

10

−4

10

−3

10

−2

10

Δt∗

−1

10

−4

10

−3

10

−2

10

Δt∗

−1

10

(a) Elastic, NF M = 0 (b) Elastic, NF M = 0.5 (c) Damped, NF M = 0 Fig. 6. (a–d) Nominally 4th order schemes: RK4 ( these figures is εxamp (%).

), AB4 (

), ABM (

), BD4 (

) with fourth order gradient (

−4

10

−3

10

−2

10

Δt∗

−1

10

(d) Damped, NF M = 0.5 ) to guide the eye. The ordinate for

218

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

However, we have seen here that, unlike the linear force model (where some schemes perform best with NFM = 0), numerical integration of the non-linear force model yields results which are either insensitive to NFM , or somewhat chaotic. Therefore it would not be beneficial to ensure that NFM = 0 for each collision, since greatest accuracy is not located at NFM = 0. 5.2. Scheme order analysis: elastic and damped collision The dependence of integration accuracy on timestep size is one of the most important measures of scheme performance since it dictates the rate at which the accuracy can be improved by reducing timestep size. The nominal order of the integration scheme (see Section 4.1) indicates what the expected actual order would theoretically be, but other factors such as discontinuities in the force and velocity profiles and derivatives may affect the quality of the numerical solution. Here we compare the actual order of the schemes against the nominal order. A physical constraint on the length of the timestep is that it must be “considerably smaller” than the duration of the particle collision, tcollision . A maximum timestep can be evaluated as a fraction of the collision time, t|max =

t collision a

(10)

Various authors propose different values for a: Thompson and Grest (1991) use a = 50, Dury and Ristow (1997) use a = 15, Langston et al. (1995a) used a = 30. The value of a can be regarded as a critical collision ‘resolution’, the minimum number of discrete points that need to be evaluated over the course of a single collision. For explicit schemes, stability needs to be considered when choosing the resolution a. Using the stability of criterion of O’Sullivan and Bray (2004) for particles with no rotation allowed, a minimum value of a = 9 is required for a stable explicit solution. To investigate the effects of poorly resolved collision events, our results have a minimum resolution of a = 5. Figs. 4–6 show the relationship between the amplitude error in position, εxamp and non-dimensional timestep, t for both an elastic and damped particle collision with different pre-collision free motion. Since Fig. 3 indicates that different scheme errors may be produced by different NFM , the two extremes (NFM = 0 and NFM = 0.5) are chosen to investigate the realised scheme order. For the nominally first order schemes, Fig. 4 confirms that FE and BD1 schemes achieve first order, regardless of whether the collision is elastic or damped, or whether the value NFM coincides with the collision start or not. The SE scheme exhibits second order accuracy in the elastic collision test, again regardless of the choice of NFM . However, with the addition of damping to the interparticulate collision, the SE integration scheme is reduced to first order accuracy, but achieves an order of magnitude less error compared to the other nominally first-order schemes. This trend might be expected from Fig. 3(a) and (b), since it is observed that the nominally first-order schemes examined here are insensitive to NFM and the SE scheme is shown to give greater accuracy. The behaviour of the nominally first-order schemes is similar using both the linear and non-linear force models, except for the fact that the SE scheme is unaffected by NFM and maintain second-order accuracy for both NFM = 0 and NFM = 0.5. The nominally second order schemes require greater computational expense or greater storage requirement than first order schemes, but the error is expected to decrease more rapidly as the timestep length is reduced. Fig. 5 illustrates that without damping, all the nominally second-order schemes examined here clearly achieve second order accuracy. When damping is added to the non-linear force model, the schemes again achieve second-order

accuracy albeit with some minor deviations from the trend. This finding differs from the corresponding linear force model (Tuley et al., 2010) where it was found that the addition of damping can reduce some of the numerical schemes observed order of accuracy. This finding was anticipated, since the non-linear damping term avoids producing a non-physical, non-zero viscous force at the onset of collision. Since the magnitude of the damping force component increases from zero, there is a smooth transition from non-colliding to colliding states; the force profile is no longer zeroorder discontinuous with the non-linear damping term. The nominally fourth-order schemes are expected to have a fourth-order rate of error reduction as timestep size decreases. Fig. 6 shows the actual trend, which for all schemes and configuration (elastic or damped and regardless of NFM ) does not consistently achieve this higher error reduction rate. Under elastic conditions, regardless of NFM (see Fig. 6(a) and (b)), it appears that the multi-point methods (AB4, ABM and BD4 schemes) do achieve fourth-order accuracy for collision resolutions less than about 50. For better resolved collisions, the order of the nominally fourth-order schemes is reduced to order 2.5. With collision damping enabled, there is little evidence of a fourth-order accuracy regime as previously, and all schemes are reduced to an accuracy order of just 1.5. These results contrast with the linear force model results (Tuley et al., 2010) whereby the multi-point and singlepoint schemes behaved dissimilarly. Here, the non-linear scheme does not produce much difference in the schemes performance for each configuration examined. These results are significant because it would suggest that if an inelastic non-linear collision force model is employed, it would be more advantageous to use any of the nominally second-order scheme since they would achieve a higher rate of error reduction than a fourth order scheme. Although the same cannot be said of the non-linear elastic case, it could be questioned whether using a nominally fourth-order scheme in well resolved collisions would offer significant advantages over a second-order scheme. It has been reported that the ‘switching effect’ of the interparticle repulsive forces does not provide a physically realistic DEM model or accurate numerical integration (Fraige & Langston, 2004). It has been speculated that the transition between free-motion and colliding states is the cause for the difference in the nominal (expected) scheme order and the actual order (Tuley et al., 2010). The results here question this suggestion. That is, If the undamped linear force model can produce fourth-order accurate results using the single-point RK4 method at NFM = 0 (avoiding the problem of integrating across a discontinuity), why is the same not true for the non-linear force model under the same conditions? This finding, along with the discrepancies in the measured order of accuracies between the linear (Tuley et al., 2010) and non-linear force models, suggests that not only is the measured scheme order affected by integration across the force discontinuity (at least in the linear case) but also by the degree of the underlying equation. 5.3. Recovery of fourth order behaviour To examine the suggestion that the degree of the underlying force model accounts for the differences which exist between the linear and non-linear scheme accuracy, the same error analysis can be conducted on the nominally fourth-order schemes assigning different values of ˛ in Eq. (4). The numerical tests in Fig. 7 show that observed order of accuracy of the numerical schemes both a function of the degree of the force model, ˛ (see Fig. 7(a)) and to a lesser extent NFM (see Fig. 7(b)). For increasing ˛, the observed order of accuracy becomes closer to the schemes nominal order. For ˛ ≥ 2 all of the the nominally fourth-order schemes exhibit fourth-order accuracy for NFM = ±0.5 and NFM = ±0, with small deviations from this trend for other values

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

219

Resolution, a 2

10

0

10

3

2

1

3

2

1

3

2

1

10 10 10

10 10 10

10 10 10

α=

α=1

α=

1 2

3

2

1

10 10 10

α=2

3 2

3

2

1

−3

−2

−1

3

2

1

10 10 10

α=

5 2

−2

10

εx

amp

(%)

−4

10

−6

10

−8

10

−10

10

−12

10

−14

10

−3

−2

−1

10 10 10

−3

−2

−1

10 10 10

−3

−2

−1

10 10 10

Δt∗

−3

−2

−1

3

2

1

10 10 10

10 10 10

(a) Elastic, NF M = 0

Resolution, a 2

10

0

10

3

2

1

10 10 10 NF M = ±0.5

3

2

1

10 10 10 NF M = −0.3

3

2

1

10 10 10 NF M = −0.1

10 10 10 NF M = 0.1

10 10 10 NF M = 0.3

−2

10

εx

amp

(%)

−4

10

−6

10

−8

10

−10

10

−12

10

−14

10

−3

−2

−1

10 10 10

−3

−2

−1

10 10 10

−3

−2

−1

10 10 10

Δt∗

−3

−2

−1

10 10 10

−3

−2

−1

10 10 10

(b) Elastic, α = 2

Fig. 7. (a and b) Nominally 4th order schemes: RK4 ( ), AB4 ( ), ABM ( show the dependance of the schemes actual order of accuracy with ˛ and NFM .

), BD4 (

of NFM . This finding confirms that the multi-point methods (AB4, ABM and BD4 schemes) may integrate across a discontinuity while preserving the nominal scheme order.

The nominally second-order schemes examined here exhibit the most complex trends with damping of all the different schemes. In the case of NFM = 0, there is general trend of increasing accuracy with increasing elasticity, with the AB2 and RK2 schemes exhibiting a peak in accuracy before the restitution coefficient becomes unity. The AB2 scheme has the smallest variation in error with restitution coefficient when NFM = 0, but one of the largest variations when NFM = 0.5. In terms of a numerical accuracy comparison between the linear and non-linear force models, when the coefficient of restitution is 0.5, the nominally second-order multi-point methods have approximately 1.5 orders of magnitude greater accuracy when NFM = 0 using the non-linear approach. However, there is no significant difference when NFM = 0.5. The single-point RK2 method on the other hand has the opposite characteristic; for NFM = 0 the non-linear model slightly lower accuracy at a coefficient of restitution of 0.5, but when NFM = 0.5, the non-linear scheme produces two orders of accuracy less error. It would appear then that the damping terms of the non-linear force model is capable of producing better accuracy results than the linear force model but only under certain conditions (such as heavily damped collisions). The nominally fourth-order schemes presented in Fig. 8(e) and (f) show that as the damping is increased, the accuracy of the

5.4. Damping Unlike the linear force model which had a damping term independent of the interparticle overlap (i.e. ˇ = 0 in Eq. (5)), the non-linear force model has a damping term is also a function of the particle overlap (in this case, ˇ = 1/2 in Eq. (5)). The advantage of this modification is that it eliminates the step-change in the force profile both at the start and end of the collision event. It is therefore of interest to see how the removal of this discontinuity affects the scheme accuracy, as compared against the linear case. All the results in Fig. 8 are presented with a non-dimensional timestep, t of 3 × 10−3 to facilitate a direct comparison with the results in Tuley et al. (2010). The results from the nominally first-order schemes in Fig. 8(a) and (b) show a trend similar to that in the linear case; BD1 and FE schemes have a constant error with coefficient of damping (regardless of NFM ) with the SE schemes error reducing as collision inelasticity is reduced (albeit without reaching the same high accuracies as in the linear case).

) with fourth order gradient (

) to guide the eye. These figures

220

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222 0

0

10

10

−1

−1

(%) εamp x

10

−2

10

εx

amp

(%)

10

−3

−3

10

10

−4

10

0.5

−2

10

−4

0.6

0.7

0.8

0.9

10

1

0.5

Restitution Coefficient, e

(a) Nominally 1st order schemes, NF M = 0 −2

0.8

0.9

1

−2

10

−3

−3

εx

amp

εxamp (%)

10

(%)

10

−4

−4

10

10

−5

10

0.5

−5

0.6

0.7

0.8

0.9

10

1

0.5

Restitution Coefficient, e

−2

0.8

0.9

1

−2

10

−3

−3

10

(%) εamp x

10

−4

εx

10

−5

−4

10

−5

10

10

−6

0.5

0.7

(d) Nominally 2nd order schemes, NF M = 0.5

10

10

0.6

Restitution Coefficient, e

(c) Nominally 2nd order schemes, NF M = 0

(%)

0.7

(b) Nominally 1st order schemes, NF M = 0.5

10

amp

0.6

Restitution Coefficient, e

−6

0.6

0.7

0.8

0.9

1

Restitution Coefficient, e

(e) Nominally 4th order schemes, NF M = 0

10

0.5

0.6

0.7

0.8

0.9

1

Restitution Coefficient, e

(f) Nominally 4th order schemes, NF M = 0.5

Fig. 8. (a–f) Integration schemes: FE ( ), SE ( ), ( ), RK2 ( ), AB2 ( ), BD2 ( ), RK4 ( position amplitude error dependence on the coefficient of restitution and particle free motion before the collision starts.

), AB4 (

), ABM (

) and BD4 (

)

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

scheme reduces non-linearly; this is true regardless of NFM . Unlike the nominally second-order schemes, here the multi-point nominally fourth-order schemes behave almost identically in each test. Again, by comparison to the linear force model in Tuley et al. (2010), it is found that across all comparable damped collision, with NFM = 0, the non-linear force model yields superior accuracy by approximately two orders of magnitude. Although the same finding is not true for NFM = 0.5, for moderate damping levels (coefficient of restitution >0.8) the non-linear scheme provides superior accuracy to the linear model.

6. Conclusions A number of conclusions can be drawn from this work comparing numerical time integration algorithms for the non-linear DEM force model, they are outlined here. In addition to the timestep size, inter-particle contact stiffness and damping, integrator accuracy also depends on the pattern of discretisation around collision boundaries. More significantly, the actual order of accuracy produced by the collision force model is a function of the degree of the underlying equation. With second degree (or greater) elastic repulsion terms, all of the nominally fourth-order schemes tested here exhibit fourth-order behaviour. The best performing nominally first-order scheme is, as has been found with the linear force model, the SE scheme. Despite not achieving the higher orders of accuracy as in the linear case, this scheme exhibits robustness (accuracy independent of NFM ) and is consistently more accurate than the other nominally first-order schemes tested here. Unlike previous tests have shown with the linear force model, there is no significant differences of the observed scheme order of the nominally second-order schemes; the multi-point and singlepoint schemes behave similarly, with the RK2 scheme achieving marginally better accuracy in the elastic non-linear case. Unlike the results from the linear force model however, all the nominally second-order schemes exhibit second-order accuracy regardless of NFM with the non-linear collision model. Although fourth-order accuracy is not achieved using any of the nominally fourth-order integrators, the Hertzian force model does generally produce higher orders of accuracy than the equivalent linear force model test (the only exception being the RK4 scheme under linear elastic conditions with NFM = 0). As with the linear force model, the nominally second-order schemes perform with a greater order of accuracy than the nominally fourth-order schemes with damping present – however, the fourth-order schemes have greater absolute accuracy at small collusion resolutions used for practical DEM simulations. The non-linear damping term yields more similar results between multi-point and single-point schemes, unlike the linear damping term which may produce greater differences in the orders of accuracy.

Appendix A. In order to assess the accuracy of the numerical schemes tested here, an analytical solution is used as the benchmark to compare against. Depending on which force model is employed, a different analytical solution is sought. Rather than explicitly requiring a complete equation of motion, only the peak interparticle overlap,

ı

pq analytical and collision duration, tcollision is required. The following sections contain the derivations for calculating the analytical parameters for the non-linear elastic and non-linear damped collision force models.

221

A.1. Non-linear, undamped analytical solution By neglecting the damping term (cˇ = 0) and applying Newton’s second law, the governing collision equation becomes: d2 x k˛ ˛ + x =0 m dt 2

(11)

multiplication by dx/dt and application of the chain-rule gives: d dt

 dx 2 dt

+

 dx 2 k

˛ ˛

m

dt

x =0

(12)

integration w.r.t. time yields:

 dx 2 dt

+

2 k˛ x˛+1 +C =0 m ˛+1

(13)

where the constant of integration, C can be found by applying the initial conditions: x|t=0 = 0

(14)

u|t=0 = u0

(15)

∴ C = −u20

(16)

to find the maximum interparticle overlap, the condition that dx/dt = 0 is applied, hence: ı

pq =

(˛ + 1) m 2 u0 2k˛

1/(˛+1) (17)

To find evaluate the collision duration, Eq. (13) is rearranged and integrated w.r.t. x, i.e.:

 ı pq 

u20

tcollision = 2 0

2 k˛ x˛+1 − m ˛+1

−1/2 dx

(18)

for which the general solution to which is:



tcollision = 2

ı

pq u0





F1 2

˛+1

1 2 k˛ ı

1 1 pq , ;1 + ; 2 ˛+1 ˛ + 1 u2 m ˛ + 1 0

 (19)

where 2 F1 (.) is Gauss’ hypergeometric function (see, for example, Gradshteyn & Ryzhik, 1994), which may be simplified further depending on the value taken by ˛. A.2. Non-linear, damped analytical solution Here we include the damping term, with ˇ = 1/2 and  = 1. In this case ˛ is chosen to be 3/2 which corresponds to the Hertzian elastic repulsion between two spheres. By applying Newton’s second law, the governing collision equation becomes: cˇ √  dx k˛ 3 d2 x 2 + + x =0 x m m dt dt 2

(20)

Unlike the non-linear elastic case, the non-linear damping term prevents a simple analytical expression for the maximum overlap and the time at which it occurs. The approach given by Ramirez et al. (1999) is followed, whereby a series solution is used to provide a semi-analytic solution. The initial conditions remain unchanged: x|t=0 = 0

(21)

u|t=0 = u0

(22)

After adopting the same procedure, it is noted that there is a discrepancy in the resultant numerical coefficients, yn between those previously published by Ramirez et al. (1999) and the current work. The comparison presented in Table 3 shows inconsistencies for y2

222

M. Danby et al. / Computers and Chemical Engineering 58 (2013) 211–222

Table 3 Numerical comparison of the coefficients yn . Note the discrepancies in values y2 and y3 between published data and those found in this work. Numerical comparison, yn n

Ramirez et al. (1999)

Current work

0 1 2 3 4 5 6 7 8 9 10 11 12 13

1.093362 −0.504455 0.260542 −0.136769

1.0933620739 −0.5044548944 0.2840430188 −0.1702207784 0.1055007091 −0.0668487139 0.0430394924 −0.0280510843 0.0184608507 −0.0122461857 0.0081775892 −0.0054914650 0.0037055047 −0.0025108978

and y3 . It is assumed that the resulting numerical coefficients used in this work are correct. These values are derived from the preliminary bn coefficients, of which b120 is the maximum coefficient calculated in this work. These have been made available through the publisher. Note that the key finding of higher-order accuracy recovery in this contribution is based an exact analytical solution of the non-linear force model. References Aoki, K. M., & Akiyama, T. (1995). Simulation studies of pressure and density wave propagations in vertically vibrated beds of granules. Physical Review E, 52(3), 3288–3291. Cleary, P. W., & Sawley, M. L. (2002). Dem modelling of industrial granular flows: 3D case studies and the effect of particle shape on hopper discharge. Applied Mathematical Modelling, 26(2), 89–111. Cundall, P. A., & Strack, O. D. L. (1979). Discrete numerical-model for granular assemblies. Geotechnique, 29(1), 47–65. Dury, C. M., & Ristow, G. H. (1997). Radial segregation in a two-dimensional rotating drum. Journal de Physique I, 7(5), 737–745. Dziugys, A., & Peters, B. (2001). An approach to simulate the motion of spherical and non-spherical fuel particles in combustion chambers. Granular Matter, 3(4), 231–265. Falcon, E., Laroche, C., Fauve, S., & Coste, C. (1998). Behavior of one inelastic ball bouncing repeatedly off the ground. European Physical Journal B, 3(1), 45–57. Fraige, F. Y., & Langston, P. A. (2004). Integration schemes and damping algorithms in distinct element models. Advanced Powder Technology, 15(2), 227–245. Gradshteyn, I. S., & Ryzhik, I. M. (1994). Table of integrals, series, and products (5th ed.). Boston: Academic Press. Haile, J. M. (1992). Molecular dynamics simulation: Elementary methods. New York: John Wiley & Sons Inc. Hertz, H., Jones, D. E., & Schott, G. A. (1896). Miscellaneous papers. London: Macmillan and Co. Hogue, C., & Newland, D. (1994). Efficient computer-simulation of moving granular particles. Powder Technology, 78(1), 51–66. Kopf, A., Paul, W., & Dunweg, B. (1997). Multiple time step integrators and momentum conservation. Computer Physics Communications, 101(1-2), 1–8. Kreyszig, E. (2006). Advanced engineering mathematics (9th ed.). London: John Wiley & Sons. Kruggel-Emden, H., Simsek, E., Rickelt, S., Wirtz, S., & Scherer, V. (2007). Review and extension of normal force models for the discrete element method. Powder Technology, 171(3), 157–173.

Kruggel-Emden, H., Sturm, M., Wirtz, S., & Scherer, V. (2008). Selection of an appropriate time integration scheme for the discrete element method (DEM). Computers & Chemical Engineering, 32(10), 2263–2279. Kuwabara, G., & Kono, K. (1987). Restitution coefficient in a collision between 2 spheres. Japanese Journal of Applied Physics Part 1: Regular Papers Short Notes & Review Papers, 26(8), 1230–1233. Lan, Y. D., & Rosato, A. D. (1997). Convection related phenomena in granular dynamics simulations of vibrated beds. Physics of Fluids, 9(12), 3615–3624. Langston, P. A., & Tuzun, U. (1994). Continuous potential discrete particle simulations of stress and velocity-fields in hoppers – Transition from fluid to granular flow. Chemical Engineering Science, 49(8), 1259–1275. Langston, P. A., Tuzun, U., & Heyes, D. M. (1995a]). Discrete element simulation of granular flow in 2D and 3D hoppers – Dependence of discharge rate and wall stress on particle interactions. Chemical Engineering Science, 50(6), 967–987. Langston, P. A., Tuzun, U., & Heyes, D. M. (1995b]). Discrete element simulation of internal-stress and flow-fields in funnel flow hoppers. Powder Technology, 85(2), 153–169. Matuttis, H. G., Luding, S., & Herrmann, H. J. (2000). Discrete element simulations of dense packings and heaps made of spherical and non-spherical particles. Powder Technology, 109(1-3), 278–292. O’Sullivan, C., & Bray, J. D. (2004). Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme. Engineering Computations, 21(2-4), 278–303. Ovesen, J. H., Petersen, H. G., & Perram, J. W. (1996). Comparison of two methods for solving linear equations occurring in molecular dynamics applications. Computer Physics Communications, 94(1), 1–18. Ramirez, R., Poschel, T., Brilliantov, N. V., & Schwager, T. (1999). Coefficient of restitution of colliding viscoelastic spheres. Physical Review E, 60(4), 4465–4472. Rougier, E., Munjiza, A., & John, N. W. A. (2004). Numerical comparison of some explicit time integration schemes used in DEM FEM/DEM and molecular dynamics. International Journal for Numerical Methods in Engineering, 61(6), 856–879. Sadd, M. H., Tai, Q. M., & Shukla, A. (1993). Contact law effects on wave-propagation in particulate materials using distinct element modeling. International Journal of Non-Linear Mechanics, 28(2), 251–265. Satoh, A. (1995a]). Molecular-dynamics simulations on internal structures of normal shock-waves in Lennard–Jones liquids. Journal of Fluids Engineering: Transactions of the ASME, 117(1), 97–103. Satoh, A. (1995b]). Stability of computational algorithms used in molecular dynamics simulations. Journal of Fluids Engineering: Transactions of the ASME, 117(3), 531–534. Schinner, A. (1999). Fast algorithms for the simulation of polygonal particles. Granular Matter, 2(1), 35–43. Schwager, T., & Poschel, T. (1998). Coefficient of normal restitution of viscous particles and cooling rate of granular gases. Physical Review E, 57(1), 650–654. Schwager, T., & Poschel, T. (2007). Coefficient of restitution and linear-dashpot model revisited. Granular Matter, 9(6), 465–469. Shida, K., Suzuki, R., & Kawai, T. (1997). Numerical error of total energy. dependence on timestep. Computer Physics Communications, 102(1-3), 59–67. Sundaram, S., & Collins, L. R. (1996). Numerical considerations in simulating a turbulent suspension of finite-volume particles. Journal of Computational Physics, 124(2), 337–350. Taguchi, Y. H. (1992). New origin of a convective motion – Elastically induced convection in granular-materials. Physical Review Letters, 69(9), 1367–1370. Thompson, P. A., & Grest, G. S. (1991). Granular flow – Friction and the dilatancy transition. Physical Review Letters, 67(13), 1751–1754. Tsuji, Y., Kawaguchi, T., & Tanaka, T. (1993). Discrete particle simulation of 2dimensional fluidized-bed. Powder Technology, 77(1), 79–87. Tsuji, Y., Tanaka, T., & Ishida, T. (1992). Lagrangian numerical-simulation of plug flow of cohesionless particles in a horizontal pipe. Powder Technology, 71(3), 239–250. Tuley, R., Danby, M., Shrimpton, J., & Palmer, M. (2010). On the optimal numerical time integration for Lagrangian DEM within implicit flow solvers. Computers & Chemical Engineering, 34(6), 886–899. Xu, B. H., & Yu, A. B. (1997). Numerical simulation of the gas–solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chemical Engineering Science, 52(16), 2785–2809. Zhang, D., & Whiten, W. J. (1996). The calculation of contact forces between particles using spring and damping models. Powder Technology, 88(1), 59–64.