Establishing stable time-steps for DEM simulations of non-collinear planar collisions with linear contact laws

SUMMARY The discrete element method, developed by Cundall and Strack, typically uses some variations of the central difference numerical integration scheme. However, like all explicit schemes, the scheme is only conditionally stable, with the stability determined by the size of the time-step. The current methods for estimating appropriate discrete element method time-steps are based on many assumptions; therefore, large factors of safety are usually applied to the time-step to ensure stability, which substantially increases the computa- tional cost of a simulation. This work introduces a general framework for estimating critical time-steps for any planar rigid body subject to linear damping and forcing. A numerical investigation of how system damp- ing, coupled with non-collinear impact, affects the critical time-step is also presented. It is shown that the critical time-step is proportional to q mk if a linear contact model is adopted, where m and k represent mass and stiffness, respectively. The term which multiplies this factor is a function of known physical parame- ters of the system. The stability of a system is independent of the initial conditions. © 2016 The Authors.


INTRODUCTION
Contact mechanics primarily involves the analysis of the impulses and reaction forces that develop as a result of two or more bodies colliding or being in contact, and the governing dynamics which ensues [1][2][3][4]. The main objectives of contact mechanics are to develop and understand methods for calculating the changes in motion, for the case of a collision between bodies [5], and to calculate interaction forces for bodies already in contact [6]. The understanding of the mechanics of contact has far-reaching applications in areas such as civil engineering, mechanical engineering, vehicle collision analysis and sport science [7].
One such method for modelling contact problems is the discrete element method (DEM) [8]. DEM is typically used to describe the behaviour of bulk granular materials in which there is a large number of particles and multiple simultaneous contacts [9]. The bodies are assumed to be rigid, but a small amount of overlap at the contact region is permitted. The impact process is modelled using a combination of springs, dampers and sliders. Using Newton's second law together with Euler's second law, the contact process can be described mathematically as a system of second-order differential equations (ODEs). 188 S.J. BURNS AND K.J. HANLEY planar cases in Section 4. In Section 5, we carry out a numerical investigation of the changes in the critical time-step under parameter variation. The paper concludes with a discussion in Section 6 that provides an insight for engineers and other researchers working with DEM.

STABILITY AS A FUNCTION OF MASS AND STIFFNESS
In this section, we will introduce and discuss one of the most frequently used techniques for selecting a stable time-step.

Stability using Rayleigh's theorem
In this paper, when we refer to current methodologies, we are referring to the approach whereby the critical time-step is calculated as a function of the mass and stiffness of the system particles. This commonly used approach makes use of a corollary of Rayleigh's theorem, the essential details of which will be given here.
Consider the following second-order system of differential equations where M is the mass matrix, C is the damping matrix, K is the stiffness matrix, F is the column vector of external forces and q, P q and R q are the displacement, velocity and acceleration vectors, respectively. Belytschko [14] derives the stability criterion for the discretised form of (1) by making use of modal decomposition to reduce the system to a single degree of freedom system. The maximum stable time-step t crit is a function of the eigenvalues of the amplification matrix and is given by for a linear, undamped system [11,15]. Belytschko [14] gives the following equation where max is the maximum eigenvalue of Equation (3) can be derived using an extension of Rayleigh's bounding theorem [14], which relates the eigenvalues of any two systems which are equivalent except for linear constraints. It is important to note, however, that this condition was derived for a second-order system of the form Recasting this system as a system of first-order differential equations gives the following equivalent relation A major shortcoming of this analysis is that it requires that the modal equations of motion will decouple. In order for this to be possible, it is necessary to impose certain restrictions on C , which are often unphysical. The method employed in [14] assumes Rayleigh damping in that the damping matrix C is directly defined as a linear combination of the mass and stiffness matrices. The main consequence of this is that it may not always be correct to use (2) to calculate the critical time-step for a system in which damping is present. This method has also not been developed to allow for collisions of non-symmetric bodies. 189

Natural frequency of the system
The principle of natural frequency of a system can be used to theorise numerical stability for the discretised form of the system [21]. The time at which the compression phase ends, t c , can be calculated in terms of the natural frequency of oscillation of the system [1] and is given by where ! is the natural frequency of oscillation with respect to the normal or tangential direction. Equation (7) is valid only when the tangential and normal dynamics are decoupled, i.e., for collinear collisions. If it is assumed that all forces and deformations occur around the contact point, it is also intuitive to consider the natural frequency of oscillation at a frame of reference located at the contact point [1]. Choosing a time-step greater than (7) could lead to particles passing through each other without an impact being detected or could lead to an excessively large overlap between particles, which, as a result, would give an unphysical reaction force leading to instability in the system. The time at which the compression period ends, given by (7), is almost the same as the critical time-step predicted by (2). It makes sense physically that one should not choose a time-step greater in value than the theoretical time for compression. As discussed in Section 2.1, a number of assumptions are necessary in order to derive (2). We propose using (7) as an alternative way of calculating critical time-steps such that Equation (8) does not require any restrictions on the system damping other than linearity with respect to the generalised coordinates. For non-collinear collisions, (8) will not be exact. However, it serves as an adequate bound for estimating the critical time-step. It is also important to highlight the physical meaning of the critical time-step, and the link between natural frequency of oscillation and theoretical time for compression facilitates this. In the next section, we will present a general framework for deriving the natural frequency of oscillation for any planar rigid body.

GENERAL FRAMEWORK
Consider the two bodies H and H 0 in free flight as shown in Figure 1(a) and in contact at point P as shown in Figure 1(b). The position and rotation of the centre of mass G of body H can be described in the n 1 n 2 plane by the coordinates q 1 and q 2 and the angle Â, and similarly, q 0 1 and q 0 2 are the coordinates and Â 0 the rotation of the centre of mass G 0 of body H 0 (Figure 1(b)).
We let be the displacement, velocity and acceleration vectors of the centre of mass of H and H 0 , respectively. In order to translate the contact forces to the centre-of-mass reference frame, we need to consider @q P @q D where L is the distance from G to P and where L 0 is the distance from G 0 to P . Further, we can define e and e 0 as the eccentricity of collision such that e D L cos.Â/ and e 0 D L 0 cos.Â 0 /: Then, using Newton's second law, the equations of motion for the two bodies at contact can be written as and where M and M 0 are, respectively, the mass matrices for H and H 0 given by and where the external forces, F and F 0 , which respectively act on H and H 0 , are given by The subscripts 1 and 2 respectively represent the components of the external forces acting in the tangential and normal directions, and R and R 0 are the external torques acting on the bodies as shown in Figure 1(a). Similarly, we define and 0 as the forces generated at impact of each body given by It is important to note that this set-up is general and does not specify the mechanism that generates the normal and tangential forces, N and T , respectively. In the context of the analysis presented 191 in this article, when we refer to external forces, we are referring to all forces which do not feature in the contact force .
As discussed in Section 2.2, it is intuitive to use a reference frame located at the contact point between the contacting rigid bodies. For this purpose, we consider (9) translated to the contact point at coordinates q P so that (12) This is the rate of change of the contact point velocities as a function of time, in which´ 1 and .´0/ 1 are the symmetric matrices given bý For the analysis of the system at contact, we do not need to consider terms which do not change throughout the impact phase, and therefore, we can neglect the functions f and f 0 ( [22]). This framework is completely general and can be used with any planar geometry, provided the moment of inertia is known, together with any contact law which is linear with respect to the generalised coordinates.

Linear contact interaction
We model the contact process using linear springs and dampers acting in the normal and tangential directions (see Figure 2). Let k 1 and c 1 be the spring and damping constants for the components acting in the tangential direction and let k 2 and c 2 be the spring and damping constants for the components acting in the normal direction. The contact forces and 0 are then given by where the particle overlap is given by the relative displacement of the contact points and is given by Defining the relative acceleration of the contact points of H and H 0 as R Q q P WD R q 0 P R q P , the equations of motion now take the form where This general framework can be used for any planar rigid-body provided the moment of inertia, together with the orientation terms given in A, B and D, are known. In the following section, we will derive the general solution to (13) in order to determine the numerical stability for this system.

TIME-STEP CALCULATIONS
In this section, we will discuss the general solution to the ODE system derived in Section 3.1 and discuss what features of the solution affect the stability of the discretised system. We will then calculate analytical expressions for the critical time-step for a number of cases.
In order to solve (13), it is first necessary to recast the second-order system as a system of firstorder ODEs. We let (13) can now be recast as the following system 0 where the system matrix 193 subject to the initial conditions The general solution to (14) is determined by calculating the eigenvalues and eigenvectors of (15) and using the initial condition vector (16) to give where 1 , 2 , 3 , 4 , and V 1 , V 2 , V 3 , V 4 are the eigenvalues and corresponding eigenvectors of (15), respectively. It is important to note that the eigenvalues and eigenvectors are not a function of the initial conditions (16). The initial conditions form part of the constants˛1,˛2,˛3,˛4, which scale the solution. As discussed in Section 2, the stability condition is only a function of the system eigenvalues. This means that for linear contact interactions, the stability of the discretised system is entirely independent of the initial conditions, and further that the natural frequency does not depend on the initial velocity at contact, and furthermore, the compression time is independent of the initial velocity.

Collinear case
For the special case of a collinear collision O B D 0; Â D Â 0 D 2 Á , the normal and tangential dynamics decouple, and (8) becomes exact as a result. This represents the case in which the centre of mass of each body can be joined by the line normal to the tangent plane passing through the contact point P , for example, in the case of two spheres colliding. The eigenvalues of (15) are then given by Using (18), (20), (8) and (6) gives Note that, given the generality of this formulation, it is possible to have contact with initial relative tangential velocity as well as relative normal velocity.

Zero damping
Consider the special case of zero damping (c 1 D c 2 D 0). The eigenvalues of (15) are now given by Using (23), (8) and (6) gives

Two-disc collision
We will consider an example system of a two-disc rigid collision to allow for a direct comparison with the current state of the art. For the case of a collision between discs of radii r and r 0 and masses m and m 0 , we have that O B D 0; Â D Â 0 D 2 . Considering the specific case analysed in [15] of an undamped pair of discs (c 1 D c 2 D 0), (22) reduces to For a two-disc collision of discs of equal mass m and equal spring constants k, (28) further reduces to This value is slightly more conservative than the q m k estimate predicted in [15]. For the case of close-packed configurations, as also considered in [15] and [16], the 0:64 q m k estimate reported here is less conservative. This is due to the fact that a close-packed configuration will tend to increase the natural frequency of the system. The estimate predicted by (28) shows that when the spring constants are equal, the tangential frequency dictates the system stability. In practice, however, the normal spring stiffness is typically three times greater than the tangential spring stiffness.

Full system
The eigenvalues of (15) can be calculated analytically; however, the expression for each eigenvalue contains a large number of terms, making it non-intuitive to give here. We instead consider a numerical investigation to examine the calculated stable time-step for a variety of system parameters as shown in Section 5.

NUMERICAL RESULTS
Here, we will focus on the numerical analysis of how the critical time-steps calculated in Section 4 are affected by the system parameters. The purpose of this is threefold. Firstly, we want to show how system damping affects the critical time-step. Secondly, we want to investigate the effect the eccentricity of collision has on the critical time-step, and thirdly, we will exploit the fact that our analysis allows for time-step calculations for any planar rigid body to investigate the effect of shape on the critical time-step.

Undamped
In this section, we will consider the critical time-steps for a number of undamped cases. The purpose of this is to enable a comparison with the damped case in order to correctly interpret the effect of damping on the system. To determine the effect of eccentricity on the critical time-step, we will use the eccentricity parameters e and e 0 defined in Section 3. For this purpose, we define a value for the moment of inertia of each body and vary the eccentricity as shown in Figure 3. This acts as a rotation about the contact point P of each body away from the collinear configuration. In Figure 3(a), we show different contours of the spring constant k 1 for zero damping in the system. For each value of k 1 , the critical time-step decreases with increasing eccentricity, showing that the system  becomes less stable the further we move away from a collinear configuration. In Figure 3(b), we show different contours of the spring constant k 2 for zero damping in the system. As for Figure 3(a), we see the same effect of decreasing stability with increasing eccentricity. It is interesting to note how the critical time-step changes with increasing k 1 . When we increase k 1 from 1 10 6 kg s 1 to 2 10 6 kg s 1 in Figure 3(a), we do not obtain a p 2 decrease in t crit as expected. However, when we increase k 1 from k 1 D 2 10 6 kg s 1 to k 1 D 4 10 6 kg s 1 , we observe the p 2 decrease in t crit . The reason for this 'nonlinearity' is due to the fact that the contour k 1 D 1 10 6 kg s 1 corresponds to the eigenvalue 1 , but the contours k 1 D 2 10 6 kg s 1 , k 1 D 4 10 6 kg s 1 correspond to the eigenvalue 3 given by (23) and (25), respectively.
In Figure 5, we investigate the effect of particle size on the critical time-step. We consider a rigid disc-rigid disc collision for systems with constant mass and investigate how the stability is affected by decreasing the ratio m m 0 , where m 6 m 0 . In Figure 4(a), we consider different contours of the total system mass, and in Figure 4(b), we consider varying contours of the spring constant k 2 . In both cases, it is clear that as the ratio m m 0 7 ! 0, t crit 7 ! 0. In Figure 4(a), we see that as the system mass Figure 5. Schematic illustration of the rigid-body orientations used for the analysis shown in Figure 6 and 10. In (a), a rigid rectangular plate-rigid rectangular plate collision is shown, in which an isolated point contact as opposed to a surface contact is assumed. In (b), a rigid rectangular plate-rigid disc collision is considered, and in (c), a rigid disc-rigid disc collision is shown.

197
is increased, the stability also increases, consistent with previous findings. Similarly, in Figure 4(b), we see that as k 2 is increased, the stability decreases, as expected.
To examine the effect that particle shape and inertia has on the critical time-step, we consider two types of laminar, namely, a rigid rectangular plate and a rigid disc (see Figure 5). In Figure 6(a), we consider a rigid rectangular plate-rigid rectangular plate collision as a function of tangential to normal spring constant ratio k 1 k 2 for varying contours of eccentricity of collision. We see that the collinear case is the most stable, as expected, and that the stability decreases as the ratio k 1 k 2 increases. In Figure 6(b), we consider three cases, namely, rigid rectangular plate-rigid rectangular plate collision, rigid rectangular plate-rigid disc collision and a rigid disc-rigid disc collision. We see that the non-collinear rigid rectangular plate-rigid rectangular plate collision is the least stable and that the collinear rigid rectangular plate-rigid disc collision is the most stable. This is perhaps a surprising result in that it is a more stable configuration than the disc-disc collision. Similar to 6(a), we see that the stability decreases as the ratio k 1 k 2 increases. In both cases, we see that for certain conditions, t crit remains constant with respect to k 1 k 2 . This corresponds to the transition between the eigenvalues 1 and 3 given by (23) and (25). A similar transition was observed in the different contours in Figure 3

Damped
In Figure 7, we present the same analysis as in Figure 3 but with system damping included. For each value, the critical time-step decreases with increasing eccentricity and increasing damping, showing that the system becomes less stable with increased damping and eccentricity. In Figure 7(b), we see the same trend. In both cases, we see that the critical time-step is almost four times smaller than in the undamped case, showing the significant effect system damping has on the stability. In Figure 8, we investigate the effect that damping in the normal direction has on the critical time-step. Figure 8(a) shows different contours for the tangential damping constant, and Figure 8(b) shows different contours for the normal spring constant. The critical time-step decreases with increased damping and increased normal force. In Figure 8(b), the critical time-step is constant until c 2 Ñ 390 kg s 1 for the contour k 2 D 1 10 6 kg s 1 . This corresponds to a transition between the eigenvalues in Section 4.4. Once a large enough damping force has been reached, one of the system eigenvalues dominates and dictates the system stability. The result that the system becomes less stable with increased damping is consistent with the findings in [23].
To examine the effect that particle shape and inertia have on the critical time-step for the damped case, we consider two different numerical investigations as shown in Figures 9 and 10. In Figure 9,  we consider a similar set-up to that considered in Figure 5 in Section 5.1. Figure 9(a) considers the critical time-step as function of the mass ratio for varying contours of system damping. Figure 9(b) considers the critical time-step as function of the mass ratio for varying contours of k 2 for constant non-zero damping. In both cases, the same trends as in Figure 5 are observed with the heightened effect because of the presence of damping. This shows that when there is a collision between two discs with one significantly larger than the other, then it is necessary to adopt a very small time-step to ensure numerical stability, particularly when damping is included. Figure 10 considers the effect that the shape of the rigid body has on the critical time-step for the damped case. In this analysis, three different cases are considered, namely, a rigid rectangular plate-rigid rectangular plate collision, rigid rectangular plate-rigid disc collision and a rigid discrigid disc collision. In Figure 10(a), the critical time-step is calculated as a function of the damping constant ratio c 1 c 2 while in Figure 10(b), the critical time-step is calculated as a function of the forcing constant ratio k 1 k 2 . The critical time-step decreases with increasing ratios of damping or forcing. Both Figure 10(a) and (b) shows the change in the critical time-step with increasing shape complexity, with the rigid rectangular plate-rigid disc collision being the most stable, as in the undamped case in Figure 6(b). In Figure 10(a), we see that for certain conditions t crit remain constant with respect to c 1 c 2 . This corresponds to the transition between the eigenvalues 1 and 3 given by (23) and (25), similar to the trend observed in Figure 6(a) and (b). However, we do not see this behaviour in Figure 10(b). The presence of damping, for this configuration, has the effect that 3 is always greater than 1 in this region of parameter space.

DISCUSSION AND CONCLUSIONS
In the framework of using system eigenvalues to determine numerical simulation stability, an extension to current methodologies has been proposed for taking into account system damping together with complexity in particle geometry. This method is based on a general framework, which analyses the contact phase in terms of contact point equations of motion. This general framework is compatible with any planar rigid body subject to a linear contact law.
A comparison between the classical critical time-step calculation, based on a single degree of freedom, undamped system, and the theoretical expression for the compression time of the contact phase is made and suggested as an alternative way for calculating the critical time-step. The framework presented in this work however can be used with the current methodology at the cost of possible inaccuracy because of the damping assumptions.
For all of the cases considered here, the critical time-step calculations were found to be proportional to q m k , consistent with current methodologies. In each instance, this factor was scaled by some combination of the symmetric matrix terms O A, O B and O D. These terms represent the mass distribution of the system. It is well known that in practice, it is necessary to multiply q m k by a 'factor of safety' in order to achieve a stable numerical simulation. This work has given a physical meaning to this safety factor in terms of physical parameters of the system, albeit just for the general, linear planar cases presented here. We have also shown that, for the linear cases analysed here, the stability of the system is independent of the initial conditions.
We performed a numerical investigation of how the critical time-step varies as a function of the system damping together with the eccentricity of collision. This analysis is completely novel and provides useful insights for engineers and other researchers who want to perform stable numerical simulations for non-symmetric systems with damping. The analytical results presented here have also confirmed previous observations in the literature showing how simulations become less stable with increased system damping. We have also shown that there are regions in parameter space in which the critical time-step remains constant with respect to increasing spring constant and damping constant. This is a novel result and is something which is of great use to researchers considering DEM simulations.
This work, in its current form, is only applicable for two contacting particles. The framework presented here, however, is general and can be extended to incorporate a close packing of DEM particles: a configuration which is often encountered in many DEM applications. As mentioned previously, the effect of increasing the number of particles will be to increase the natural frequency of the system and thus reduce the critical time-step. We expect that there may also be some coupling between the eccentricity of each body, which may yield some interesting results.
This work is the first step in an attempt to calculate critical time-steps for a number of general configurations. We anticipate that when this analysis is extended to nonlinear interaction laws, the factor of safety will be a function of the initial conditions, particularly the initial contact relative velocity.