Critical time step for discrete element method simulations of convex particles with central symmetry

A time step must be selected for any explicit discrete element method (DEM) simulation. This time step must be small enough to ensure a stable simulation but should not be overly conservative for computational efficiency. There are established methods to estimate critical time steps for simulations of spherical particles. However, there is a comparable lack of guidance on choosing time steps for DEM simulations involving nonspherical particles: a fact which is increasingly problematic as simulations of nonspherical particles become more commonplace. In this article, the eigenvalues of the amplification matrix are used to develop an explicit formula for the critical time step for a range of shapes including ellipsoids, convex superquadrics and convex, central symmetric polyhedra. This derivation is based on a linear analysis and applies to both underdamped and overdamped systems. The dependence on the particle mass and contact stiffness expected for a system of spheres is recovered. For a fixed particle mass, as particle shape becomes increasingly nonspherical, the critical time step decreases nonlinearly. Thus, estimating a critical time step by assuming a sphere of equivalent volume may not always be conservative.


INTRODUCTION
The discrete element method (DEM) was proposed by Cundall and Strack 1 for modeling systems of particles and their interactions. The particles in DEM are modeled as rigid bodies which are allowed to overlap when they come into contact, from which interparticle forces are calculated using a contact model. Particle motions can be described as a system of second-order differential equations which require solution using a suitable numerical integration scheme. Typically an explicit, conditionally stable integration scheme of second-order accuracy is chosen. 2 Stability is contingent on choosing a sufficiently small time step. However, as reducing the time step increases the computational cost of a simulation, it is desirable that the chosen time step is not unnecessarily small. This is an important consideration since DEM simulations often require considerable computational effort: a fact which has motivated the adoption of high-performance computing and the development of parallelized codes. DEM often uses disks (2D) or spheres (3D) as the shape of element, again to reduce the computational cost. When a linear contact model is adopted, the critical time step has a √ m k dependency where m and k, respectively, represent some measure of particle mass and contact stiffness. Cundall and Strack 1 estimated a time step with this dependency for a single particle of mass m connected to the ground by a spring of stiffness k. Belytschko 3 proposed modal decomposition of the equations of motion to obtain a simple system with a single degree of freedom for explicit finite element analyses. The eigenvalues of the amplification matrix for this system cannot exceed 1 in magnitude 4 which gives the maximum stable time step. O'Sullivan and Bray 5 drew an analogy between particles in a DEM simulation and nodes of a finite element mesh. Based on this, they found critical time steps for regular packings of uniform disks and spheres. Even though these approaches are not strictly applicable to nonlinear systems, they have nonetheless been extended to systems of disks/spheres with nonlinear contact models [6][7][8] and rotational resistance. 9 Otsubo et al. 10 estimated the critical time step using a nonlinear contact model based on eigenvalue decomposition of the global mass and stiffness matrices for idealized and randomly packed assemblies of spheres. Tavarez and Plesha 11 constructed mass and stiffness matrices for a DEM unit cell of monosized spheres to obtain an upper-bound estimate of the critical time step. Burns et al. 12 developed a framework for selecting a stable time step for both linear and nonlinear interactions of spheres by analyzing the equations of motion as a nonlinear map. The most common approach to estimate a critical time step for nonlinear contact between spheres is fundamentally based on the propagation of Rayleigh waves through a particle assembly. 13,14 These prior studies of the critical time step considered disks or spheres. However, advances in computational power have made it increasingly feasible to run DEM simulations using particle shapes other than disks or spheres. Lu et al. 15 discuss the increased scientific interest in modeling systems of nonspherical particles and the various methods being used to achieve this, for example, ellipsoids, superquadrics (a generalization of ellipsoids), polyhedra or rigid clusters composed of multiple spheres. At present, time steps for nonspherical particle systems are often found by trial and error though some limited progress has been made to rigorously establish critical time steps. Hart et al. 16 conservatively recommended a critical time step proportional to √ m k for the block-based 3DEC code, taking m as the smallest particle mass and k as the largest contact stiffness in the system. The commercial DEM code PFC 17 automatically estimates the critical time step based on two constraints: one using contact stiffness matrices under the assumption that the degrees of freedom decouple from each other; the other a kinematic constraint which depends on the maximum components of translational velocity and acceleration among all particles in the system. The critical time step is taken as the minimum value. 17 This methodology is applied for all particle shapes including spheres, multisphere clusters, and polyhedra. However, as the maximum components of translational velocity and acceleration cannot be known before analysis, the critical time step is recomputed during the simulation based on the continually changing particle velocities and accelerations which "may be relatively expensive". 17 Burns et al. 18 obtained critical time steps proportional to √ m k for any planar rigid body subject to linear damping and forcing. A limitation of this analysis is that particle rotations have one degree of freedom, within the contact plane, which does not represent all possible particle kinematics.
In this article, we build upon this previous work to establish critical time steps for a subset of nonspherical particles including ellipsoids, convex superquadrics, and polyhedra which are both convex and central symmetric. Two analyses are compared: one based on the application of Belytschko's criterion to the characteristic equation at the interparticle contact and the other based on the eigenvalues of the amplification matrix. Both assume a linear contact model and damping.

RELATIVE ACCELERATION BETWEEN TWO CONTACTING PARTICLES
Consider two arbitrarily shaped particles with centers of mass x 1 and x 2 , contacting at x c as shown in Figure 1. Force F and moment M act on particle 2. The global moment on each particle can be expressed for any arbitrary global coordinate system as F I G U R E 1 Illustration of two nonspherical particles in contact Moments on each particle may be converted from a global to a local coordinate system using rotation matrices: where The rotation matrices can be written in terms of unit quaternions: Suppose the moments of inertia of these two particles are: These can be computed knowing each particle's geometry. We define a, b, c to be the half-lengths in the x, y, z directions of an axis-aligned bounding box around any particle in its local coordinate system. The local rotational accelerations of the particles arėL where L1 and L2 are the rotational velocities of the two particles in their local coordinate systems. These can be returned to the global coordinate system using the rotation matrices: The translational accelerations of the two particles at the contact point arë Finally, the relative translational and rotational accelerations at the contact may be computed:

MAXIMUM RELATIVE TRANSLATIONAL ACCELERATION AT THE CONTACT
In the commercial PFC 17 code, the kinematic time step constraint, which in addition to a stiffness constraint determines the simulation time step, is based on the maximum components of translational velocity and acceleration among all particles in the system. For this two-particle case, we have made a similar assumption that the maximum relative translational acceleration at the contact, that is, the maximum magnitude ofẍ c , determines the critical time step.ẍ c appears in the differential equation describing the system's dynamics at the contact: where M, C, and K refer to the particle mass, damping, and contact stiffness matrices, respectively. 12 Imposing contact models and integrators, an amplification matrix can be obtained. The magnitude ofẍ c depends on up to 18 independent parameters for moments of inertia, as shown by Equations (9) and (10). In order to find a practical solution, several additional assumptions are made to constrain the analysis: 1. Only a single contact is possible between two particles; 2. Each particle is central symmetric with a local coordinate system coinciding with the principal axes of moment of inertia, that is, I ij = 0 when i ≠ j; 3. The semiaxis lengths a ≥ b ≥ c and I zz ≥ I yy ≥ I xx (the former ensuring the latter for most shapes).
Therefore, each particle has three nonzero parameters defining its moment of inertia: I xx , I yy , I zz . These assumptions restrict the applicability of this analysis to a subset of shapes including ellipsoids, convex superquadrics, and polyhedra which are both convex and central symmetric. It is also noted that this analysis is limited to a single contact per particle; multiple simultaneous contacts require the critical time step to be reduced from the two-particle estimate derived in this article. 10 The location of the contact point x c on the surface of particle 1 and the orientation of particle 2 relative to particle 1 both affectẍ c . Figure 2 shows three of the infinitely many configurations of particles 1 and 2 in center-to-center and tip-to-tip contact, in which the two ellipsoids both have semiaxis lengths of a = 3 m, b = 2 m, and c = 1 m, respectively, from major to minor. It can be shown that the maximumẍ c arises for orthogonal "tip-to-tip" contacts. A graphical demonstration of this is given in Appendix A1. The expressions for the relative translational acceleration at the contact are shown as Equations (22)-(24) for Figure 2(A-C), respectively, obtained by combining Equations (1)- (19):

F I G U R E 2 Three examples
showing particle 2 changing its position and orientation relative to particle 1 We can write a generic expression forẍ c as where A x , A y , A z are comprehensive inverse inertia coefficients and Ω x , Ω y , Ω z are terms related to the particles' rotational velocities. With the assumption that a ≥ b ≥ c and I zz ≥ I yy ≥ I xx , we can state that Equation (25) shows that, for nonspherical particles, the relative translational acceleration at the contact is influenced by two factors: the contact force, which is related to the interparticle overlap, and particle rotation. We denote these two terms as the "translation-dominant" and "rotation-dominant" terms which we analyse separately, that is, critical time steps are obtained by neglecting either the "rotation-dominant" or "translation-dominant" terms in Sections 4 and 5, respectively. In practice, the translation-dominant time step is several orders of magnitude smaller than the rotation-dominant time step in most scenarios and hence determines the critical time step.

Belytschko's criterion
For an explicit finite element method analysis, Belytschko 3 expressed the critical time step as for a linear system with Rayleigh damping where max is the largest natural angular frequency of the system. The linear contact model is given by Equation (29) in which k j and c j are the contact stiffness and viscous damping coefficient, respectively, in the jth direction: Combining Equations (27) and (29), the ordinary differential equation of the contact point is which has the characteristic equation The roots of the characteristic equation are The overdamped case, for which c 2 j A 2 max − 4k j A max > 0, does not have a natural angular frequency. Hence Δt cri cannot be derived from Belytschko's criterion. The underdamped case, for which c 2 j A 2 max − 4k j A max < 0, has a maximum natural angular frequency 19 According to Equation (28), the critical time step is In the undamped case, c j = 0 so Equation (34) has some similarity to the collision time between particles reported by Silbert et al., 20 although that has no dependence on the particle shape: A fraction of this collision time has been adopted in other studies as a suitable simulation time step. The critical damping is obtained by setting c 2 j A 2 max − 4k j A max = 0: Since Belytschko's criterion (Equation (28)) requires the system to have a natural angular frequency, it is applicable only to the underdamped case. The normal and shear damping coefficients need to be maintained below the critical value (Equation (37)) to make this approach valid: a significant limitation. In subsection 4.2, we propose an alternative which does not have this limitation, that is, is applicable irrespective of the choice of damping coefficients.

Amplification matrix method
Linearity is a requirement of the amplification matrix method 18 by which the critical time step can be derived from the eigenvalues of the amplification matrix. The linear contact model adopted is shown as Equation (29). It has previously been shown that the Euler and velocity Verlet integrators give identical critical time steps for undamped systems. 21 An Euler integrator is adopted which gives two equations: n and n + 1 denote successive time steps. Combining this integrator with the linear contact model results in Equation (40) where K j = −k j A max and C j = −c j A max , j = 1, 2, 3: Thus, the amplification matrix is 3 × 3. Its eigenvalues can be obtained analytically using the characteristic polynomial given by Equation (41). It is noted that the commonly used velocity Verlet integrator and its variants would result in a 5 × 5 amplification matrix.
The three eigenvalues are | | < 1 is required for stability, with the critical time step being found from | | = 1. Considering only the real , which has two solutions for D: . Equating these solutions to Equation (43) gives two meaningful values for the time step: A third solution, Δt = − 1 3C j , is disregarded since this gives D = 0 on substitution into Equation (43).
Since C j , K j < 0, < 0. Therefore, the critical time step may be expressed as In the undamped case, C j = 0 so which matches the critical time step from Belytschko's criterion in Equation (35). For physical realism, 22 the normal contact stiffness k n ≥ k s , the shear contact stiffness. Conservatively taking the larger, for the undamped case. Even though Equation (48) contains k n , the shear direction is usually the one which limits the critical time step as implied by Equations (22)-(24). This is consistent with the findings of Tu and Andrade 7 and Burns et al., 12 among others.

Special case of undamped spheres
For two identical spheres of radius r and mass m, a 1 = a 2 = r, I xx1 = I xx2 = 2 5 mr 2 so A max = 7 m . Substituting into Equation (48), we obtain Thus, the expected result of is considerably more conservative than = √ 2 found by Burns et al. 12 for the same two-sphere configuration. A large factor of safety must be applied for multicontact situations comprising more than two spheres 10 ; for multiple contacts, a wide range of values have been recommended, for example, 0.17, 5 0.5, 11 or 1.27. 7

CRITICAL TIME STEP FOR ROTATION-DOMINANT SCENARIOS
One could envisage a scenario in which particles have large initial rotational velocities, or particles acquire high rotational velocities via moments induced by fluid motion, for example, through coupling with computational fluid dynamics. If the translation-dominant term in Equation (25) is neglected, using Equations (1)- (18) and (20) the following expressions are obtained for Figure 2(A-C), respectively: Defining the maximum magnitude of particle rotational velocity max as max(| G |),  . (54) In addition, note that Combining Equations (53) and (55), . (56) For a system containing particles of different types, B max can be the maximum obtainable, leading to the following rotation-dominant critical time step: This expression requires estimation or evaluation of the maximum particle angular velocity during a simulation.

NUMERICS
As the translation-dominant time step is generally orders of magnitude smaller than the rotation-dominant time step, that is the focus of these numerics. However, as an illustration, subsection 6.4 compares the translation-dominant and rotation-dominant time steps for a study of particle aspect ratio.

Damping coefficient
Three equations are presented in subsections 4.1 and 4.2 to determine a critical time step when viscous damping is valid: 1. Equation (34) derived based on Belytschko's criterion. Applicable only to underdamped systems; considers particle shape in its formulation 2. Equation (36) from Silbert et al. 20 Applicable only to underdamped systems; no consideration of particle shape 3. Equation (46) derived based on the amplification matrix method. Applicable to both underdamped and overdamped systems; considers particle shape in its formulation To enable a fair comparison across a range of damping coefficients, consider the spherical particle case in subsection 4.3 where A max = 7 m . Suppose k n = k s = k and c n = c s = c. We define a normalized dimensionless number = c √ mk to quantify the amount of damping present. The three equations listed above can be rewritten for identical spherical particles in terms of as follows: Δt cri for Equations (34), (36), and (46), respectively. Δt cri √ k m is dimensionless. Figure 3 shows the evolution of Δt cri √ k m with . The curves for Equation (58), derived based on Belytschko's criterion, 3 and Equation (59) (60) are similar-and are identical in the undamped case as already seen from Equations (35) and (48). Since the critical time step based on the amplification matrix method is the most generally applicable, that is adopted for the remaining numerics.

Contact stiffness
Consider the two ellipsoids shown in Figure 2 as representative particles which satisfy the assumptions in Section 3. Each ellipsoid is identical, with density of 1000 kg m −3 and contact stiffness k n = k s ranging from 1 × 10 10 to 10 × 10 10 N m −1 . We introduce a new dimensionless number, N, to quantify the damping for nonspherical particles in terms of A max : c n = c s = N √ k n ∕A max . The evolution of Δt cri with k n is shown in Figure 4 for N = 0, 1, 2, 3, 4, or 5. Decreasing the particle stiffness increases the critical time step irrespective of damping.

Particle density
The same two ellipsoids are adopted as in subsection 6.2. Now the contact stiffness is fixed at k n = k s = 1 × 10 10 N m −1 while the density is varied from 500 to 5000 kg m −3 . Figure 5 shows that increasing the particle density, and hence mass, increases the critical time step. This is the motivation for "density scaling" in which the particle densities are artificially increased by orders of magnitude in quasi-static simulations, enabling an increase of time step and a reduction of computation time.

Particle aspect ratio
For a nonspherical particle, one would expect a change of aspect ratio to change a particle's moment of inertia, and thereby change the critical time step in accordance with Section 4. Choosing the two standard ellipsoids with density of 1000 kg m −3 , contact stiffness of k n = k s = 1 × 10 10 N m −1 and no damping, we fix the volume V = 4 3 abc = 8 m 3 of each ellipsoid and systematically vary two aspect ratios, defined as e 1 = a b , e 2 = b c . Since a ≥ b ≥ c, e 1 , e 2 ≥ 1. The maximum rotational velocity has been set at max = 1 rad s −1 . Figure 6 shows the evolution of the critical time step with aspect ratios e 1 or e 2 , considering both the translation-and rotation-dominant expressions for the critical time step. Using the F I G U R E 6 Critical time step for two contacting ellipsoids as the particle aspect ratios e 1 and e 2 are systematically varied from 1 to 5 at a fixed normal contact stiffness k n = 1 × 10 10 N m −1 , particle density = 1000 kg m −3 , and particle volume V = 8 m 3 translation-dominant expression, the critical time step decreases as the aspect ratios e 1 or e 2 increase, that is, a nonspherical particle necessitates a smaller critical time step to ensure stability than a spherical particle of equivalent volume. The same is true using the rotation-dominant expression. It is therefore not advisable to assume equivalent diameters in order to calculate a simulation time step if particle aspect ratios differ significantly from one. The translation-dominant time step is more than four orders of magnitude smaller than the rotation-dominant time step for this representative situation.

CONCLUSIONS
This article has proposed a practical criterion for determining the critical time step for nonspherical particles which are convex and possess central symmetry, for example, ellipsoids, convex superquadrics and certain polyhedra. The derivation using the amplification matrix method is based on linearity so is limited to a linear contact model with optional viscous damping. The derived expression, which applies across the whole range of damping from undamped to overdamped, has the expected √ m k dependence on particle mass and contact stiffness. Hence, as the mass increases or the stiffness decreases, the critical time step increases. A second methodology based on Belytschko's criterion yielded the same critical time step as the amplification matrix method in the absence of damping. As particle shapes become increasingly nonspherical, quantified in this article by the deviation from unity of the aspect ratios of ellipsoids, the critical time step reduces nonlinearly. This indicates that the assumption of an "equivalent sphere" for estimating a critical time step is potentially nonconservative.

APPENDIX A. MAXIMUM RELATIVE ACCELERATION ASSOCIATED WITH ORTHOGONAL CONTACT
To demonstrate that orthogonal contact leads to the maximum relative translational acceleration at the contact, that is, the maximum magnitude ofẍ c , consider the same two ellipsoids as in Section 6, both with unit density. The center of ellipsoid 1 is fixed at the origin (0, 0, 0) with an axis-aligned orientation (q x = q y = q z = 0, q w = 1). The contact force is fixed at (1, 1, 1) for this analysis and the particle rotational velocities are zero for both particles in this analysis. Overlaps are restricted to negligible values. The configuration of the contacting ellipsoid 2 is altered in two different ways: Change center of ellipsoid 2 Move the center of ellipsoid 2 without changing its fixed orientation, giving a range of contact points on the surface of ellipsoid 1 Change orientation of ellipsoid 2 Change the orientation of ellipsoid 2 without changing the contact point between the ellipsoids Change center of ellipsoid 2 The orientation of ellipsoid 2 is fixed, with a rotation of 45 • around the minus y axis (q x = q z = 0, q y = − sin( 8 ), q w = cos( 8 )). 23 Its center is moved within Octant I (x, y, z > 0) as shown in Figure A1(A) to generate contact points covering one-eighth of the surface of ellipsoid 1. Two scenarios are considered as depicted in Figure A1(A,C). The only difference between these scenarios is a 90 • rotation of ellipsoid 2 around its local x axis. Figure A2 shows the partial surface of ellipsoid 1, colored by the magnitude of relative translational acceleration at the contact. For both scenarios, the maximum acceleration occurs at the longitudinal tip of ellipsoid 1. The relative positions of the ellipsoids whenẍ c is a maximum is shown in Figure A1(B,D) for scenarios 1 and 2, respectively.

Change orientation of ellipsoid 2
The contact point is fixed, within a small tolerance, at (3, 0, 0): the tip of ellipsoid 1 at whichẍ c is a maximum. The orientation of ellipsoid 2 is changed by moving its center point within Octant I as indicated in Figure A3(A,C). The same two scenarios are considered as before, differing only by a 90 • rotation of ellipsoid 2 around its local x axis. Figure A4 shows the surface formed from the center points of ellipsoid 2, colored as before by the magnitude ofẍ c . For both scenarios, the magnitude ofẍ c approaches a similar maximum value at three distinct orientations. Two of these orientations, shown in red at the left and bottom of Figure A4, correspond to nonorthogonal contact configurations. The

F I G U R E A1
Illustration of the two scenarios for the center of ellipsoid 2 moving relative to ellipsoid 1 (A and C), and the configuration of the two ellipsoids when the magnitude of relative translational acceleration at the contact is a maximum (B and D)

F I G U R E A3
Illustration of the two scenarios in which the orientation of ellipsoid 2, and hence the position of its center, are changed relative to ellipsoid 1 (A and C), and the configuration of the two ellipsoids when the magnitude of relative translational acceleration at the contact is a maximum (B and D)

F I G U R E A4
Surface formed by the center points of ellipsoid 2 as its orientation is changed, colored by the magnitude of relative translational acceleration at the contact third, at z = 3 m, is an orthogonal configuration which is more convenient for analysis than the nonorthogonal configurations. Therefore, the orthogonal configurations shown in Figure A3(B,D), for scenarios 1 and 2, respectively, are identified as those giving the maximum magnitude of relative translational acceleration at the contact.