Time reversibility of the discrete element method

An algorithm is presented for discrete element method simulations of energy‐conserving systems of frictionless, spherical particles in a reversed‐time frame. This algorithm is verified, within the limits of round‐off error, through implementation in the LAMMPS code. Mechanisms for energy dissipation such as interparticle friction, damping, rotational resistance, particle crushing, or bond breakage cannot be incorporated into this algorithm without causing time irreversibility. This theoretical development is applied to critical‐state soil mechanics as an exemplar. It is shown that the convergence of soil samples, which differ only in terms of their initial void ratio, to the same critical state requires the presence of shear forces and frictional dissipation within the soil system.

explicit time integrator for DEM as these integrators strike a good balance between computational cost, stability, and accuracy. 9 Velocity-Verlet and similar integrators are symplectic and time-reversible, neglecting round-off errors in typical computer implementations. 10 The rationale for exploring the time reversibility of DEM is different from that of molecular dynamics. Time reversibility is of significance in molecular dynamics because it guarantees the conservation of energy and angular momentum. Most DEM simulations, by contrast, dissipate a large proportion of the energy input to the system through mechanisms such as frictional sliding, damping, particle crushing or bond breakage. Hence, the primary motivation for investigating the time reversibility of DEM is its implications for existing theories of granular behavior. In this paper, the implications of time reversibility on critical-state soil mechanics is considered as an exemplar.
In Section 2, a typical soft-sphere DEM algorithm is summarised. A reversed-time algorithm for an energy-conserving system of frictionless spherical particles is presented in Section 3 along with a description of its implementation in the LAMMPS code 11 and some results of verification cases. Sections 4 and 5 consider the time reversibility of more complicated particle systems including shear forces, friction or damping. The implications of these findings for critical-state soil mechanics are discussed in Section 6 before concluding in Section 7.

CONVENTIONAL FORWARD-TIME DEM ALGORITHM
The second-order velocity-Verlet integration scheme is considered here due to its widespread use in DEM codes. As with all explicit numerical integration schemes, its stability is contingent on choosing a time step, Δt, which is smaller than some critical value. 12 For simplicity, we assume a fixed Δt and spherical particles. However, it is noted that the assumption of a fixed time step is not essential to ensure the time reversibility of a geometric numerical integration scheme. In particular, the Störmer-Verlet method can be formulated with an adaptive time step size while retaining its time reversibility. 13,14 Consider one particle of radius r, mass m and moment of inertia I = 2 5 mr 2 within a system of particles at a time t during a simulation. At the start of the time step t → t + Δt, the particle's translational velocity v t and rotational velocity t are partially updated based on the net (resultant) force and moment, f t and M t , according to Equations (1) and (2): The displacement x of the particle is updated using this "half-step" velocity: The rotation of the particle may be updated similarly using t+ Δt 2 though the accumulated rotation is not routinely stored for spheres. The interparticle overlap and increment of relative shear displacement at each contact c, respectively, denoted as U n c,t+Δt and U s c,t+Δt , may then be updated. These are used to calculate the contact forces f c,t+Δt . The net force acting on the particle is an accumulation of j interparticle contact forces, as well as body forces such as gravity or forces imposed by moving fluids on the particle such as drag, buoyancy or lift, where present. Damping also affects the net force acting on the particle. Similarly, the net moment (or torque) is a summation of contributions from j interparticle contacts along with optional contributions from fluids or rotational (rolling and/or twisting) resistance.
Considering only the interparticle contacts: where n and s superscripts are used to distinguish normal and shear/tangential components of contact force. The contribution of each interparticle contact to moment, calculated as the product of shear force and the displacement of the interparticle contact from the particle's centroid x r c,t+Δt , does not require incremental calculation. Neither does each normal contact force (Equation (6)), generally calculated using either a linear contact model where k n is a constant normal spring stiffness or a Hertzian model where H n c = 4G is a function of the radii of the contacting particles (r a , r b ), particle shear modulus G and Poisson's ratio .
Each tangential contact force requires incremental calculation, that is, f s c,t+Δt is computed from the value at the preceding time step, f s c,t . Before updating f s c,t , it must be mapped onto the contact plane at time step t + Δt as the contact normal may have changed and/or the particle may have rotated about the contact normal between time steps t and t + Δt. There are several approaches to perform this "shear force remapping" in DEM codes; 15 for this study, the following three-step approach has been adopted: f s,A c,t is the projection of the existing shear force at time step t onto the new contact plane defined by the normal unit vector n. 16 Equation (8) where k s is a constant shear spring stiffness and H s c = H n c √ U n c,t+Δt The accumulated shear force is set to zero when an interparticle contact is lost. Finally, the magnitude of f s c,t+Δt is limited to |f n c,t+Δt | by a Coulomb friction criterion where is the interparticle friction coefficient.
At the end of the time step, the velocity update is completed using Equations (11) and (12):

REVERSIBILITY OF UNDAMPED FRICTIONLESS SYSTEMS
Consider in the first instance the idealised case of frictionless particles without damping or other complications of the basic DEM algorithm. This significantly simplifies the algorithm presented in Section 2 by eliminating shear forces, and hence the net moment, Equation (5), becomes zero. Since the angular velocities are unchanged throughout the simulation, Equations (2) and (12) do not require calculation. Thus, the algorithm in the forward-time frame simplifies to: (14) 3 2 if Hertzian.
Even though the velocity-Verlet integration scheme is symplectic, at first appearance it seems that Equations (13), (14), (15), (16), taken as a whole, cannot be reversed in time. In the forward-time frame, the net force is obtained from a summation of contact forces; in the reversed-time frame, it might seem unavoidable that the net force would have to be apportioned in some way to the contacts which cannot be done uniquely. In fact, that is not the case: the contact forces can be obtained from the particle displacements, which can be updated prior to the force calculation.
There is therefore a change in the order in which the equations are computed in the reversed-time frame compared to the forward-time frame as shown in Equations (17), (18), (19), (20): (18) if Hertzian.
Gravitational force, as a body force of constant magnitude, may be readily incorporated into this reversed-time algorithm by addition to f t .

Numerical implementation
The open-source LAMMPS code 11 was adapted to operate in a reversed-time frame. Implementing Equations (17), (18), (19), (20) in LAMMPS was straightforward. The signs of Equations (13), (14), and (16) were switched in the initial_integrate and final_integrate functions of FixNVESphere to decrement rather than increment time. Time steps still counted forward to avoid more substantial modifications of the code. The order of function calls in Verlet::run was altered so that end_of_step functions, ordinarily invoked at the end of a time step as their name implies, were invoked at the beginning of a time step. Furthermore, data output was moved to the very start of a time step to ensure that data files are consistent with the conventional forward-time LAMMPS implementation. In addition to the fundamental calculation cycle described in Section 2, every practical MD/DEM code contains algorithms to detect contacting particles efficiently. LAMMPS uses link-cell methods to create and update Verlet neighbor lists as required. 11 Furthermore, most commonly used DEM codes are parallelised to take advantage of modern hardware and enable very large-scale simulations. LAMMPS employs a domain decomposition approach, using MPI to transfer particle data between adjoining processors on each time step. All of this "bookkeeping" works equally as well in a reversed-time frame as in the conventional forward-time frame and did not require modification for this study.
One addition made to the main LAMMPS code was local damping. This was used to prepare a sample for oedometric compression in Section 3.2.2. Local damping acts on a particle. 19 Different variants exist; the simple one implemented in LAMMPS employs the magnitude of the net force |f t+Δt |, the sign of particle velocity sgnv t+ Δt 2 and a scalar local damping coefficient l : where the d superscript denotes damping.
The other important addition made to LAMMPS was the computation of stress for a system of spheres, both the representative 3 × 3 stress tensor for each particle and the average stress in an assembly. The former is obtainable for a particle of volume V p at time t as while the latter is calculated as for N particles within a region of volume V. A detailed derivation is provided in literature. 19,20 In the forward-time frame, the representative particle stresses are computed after calculating the contact forces using Equations (6),(7),(8),(9),(10). These stresses remain valid until the same point in the following time step, that is, p t and t are active during the time step spanning t → t + Δt. In the reversed-time frame, the stresses are out of sync: p t and t are active during the time step t → t − Δt. In practice, the error introduced is small as the average stresses on the assembly, which are required for the quasi-static, stress-controlled verification simulation in Section 3.2.3, change little during successive time steps. That simulation made use of periodic boundary conditions. Whenever a periodic cell changes in size or shape, particles within the cell are remapped to the deforming space. 21 In the absence of friction, this remapping is unproblematic in the reversed-time frame.

Verification cases
Three verification cases were run after implementing the frictionless reversed-time algorithm as described in Section 3.1.
In each case, a conventional forward-time DEM simulation was run until the desired end state was reached. That state was saved in a restart file which was imported as the starting point for the analogous reversed-time simulation. The reversed-time simulation was run for the same number of time steps before comparing with the forward-time simulation. It is noted that these verification simulations are not perfectly reversible due to round-off errors. In principle these round-off errors could be eliminated by using integer arithmetic, 4 but that was not considered to be necessary for this demonstration.

Two-particle collision
Two identical particles, A and B, of diameter 1 cm and density 1000 kg m −3 were placed in a large cell. A was centred at the origin; B was placed a distance of 0.0101 m along the x-axis to give an initial separation between A and B of 100 m. A and B were assigned initial translational velocities of [0.9; 0.1; 0.0] m s −1 and [0.0; 0.0; 0.1] m s −1 , respectively. A Hertzian contact model was adopted with a shear modulus G = 2 MPa and Poisson's ratio = 0.2. Δt was set at 100 ns: far lower than the critical value of 325 s appropriate for a normal incident velocity around 1 m s −1 . 22 15 000 time steps, that is, 1.5 ms, were run in both the forward-and reversed-time frames: sufficient time for the particles to collide and fully separate. Figure 1 shows that much of the momentum of particle A is transferred to particle B by the collision. The plots for the forward-and reversed-time frames are overlapping, showing an excellent correspondence between the algorithms.

Strain-controlled oedometer test
Nine hundred spherical particles were poured into a vertical cylinder of radius 0.1 m, closed off at the bottom (z = 0 m) by a flat platen as illustrated in Figure 2A. The particle diameters in mm were randomly selected from the uniform distribution U (9,11). The particles were allowed to settle under gravity in the presence of local damping (Equation (21)) until all of the kinetic energy had dissipated. Then a second flat platen was introduced at z = 0.2 m to fully confine the The upper platen was moved downwards at a constant velocity of −0.5 m s −1 to compress the sample before pausing and allowing the kinetic energy to dissipate through damping. The forces on the lower and upper platens after this process were −323 N and 318 N, respectively (5 N difference in magnitudes due to particle weight). A solid fraction of 0.606 was attained; the sample at this stage is sketched in Figure 2C. Damping was then disabled and the sample was further compressed by moving the upper platen at a constant low velocity of −0.01 m s −1 for a total of 1.5× 10 6 time steps, that is, 75 ms. During this oedometric compression, the solid fraction increased to 0.625 while the forces on the platens increased to approximately 700 N in magnitude, creating the final state represented in Figure 2D. This final state was imported to the reversed-time frame code, the velocity of the upper platen was reversed to 0.01 m s −1 and a further 1.5× 10 6 time steps were run to recover the original solid fraction of 0.606 and forces of −323/318 N. Although excellent agreement was achieved between the forward-and reversed-time simulations as shown in Figure 3, the oscillations, which are particularly marked for the upper platen, make it difficult to say definitively whether this is truly a reversible process. These oscillations arise due to the lack of an energy dissipation mechanism during oedometric compression. The reversibility of energy dissipation mechanisms is discussed in Section 5.

3.2.3
Face-centered cubic packing Thornton 23 developed expressions for the peak stress ratios in a face-centered cubic (fcc) assembly of uniform rigid spheres subjected to plane strain and triaxial conditions. For the axisymmetric triaxial case with 2 = 3 , 23 where 1∕2∕3 are the major, intermediate, and minor principal stresses, respectively. 1 3 = 2 in the frictionless case ( = 0).
Each sphere has 12 contacts with neighbouring spheres in a fcc assembly. The failure mechanism described by Thornton 23 results in the loss of four contacts under triaxial conditions, that is, a drop in coordination number from 12 to 8. One hundred and twenty-eight monosized spheres of diameter 40 m were created in a fcc packing as shown in the inset on Figure 4. A linear contact model was used with k n = 1 × 10 11 N m: a very large value to approximate the rigidity required by the analytical solution. The assembly was compressed isotropically to 150 kPa before "drained" triaxial compression at a fixed strain rate of −1 × 10 −6 s −1 in the z direction, with fixed 2 = 3 = 150 kPa in the x and y directions. Figure 4 shows that the expected analytical solutions were obtained: 1 3 = 2 and the coordination number drops from 12 to 8. The final state was imported to the reversed-time frame code, with a reversal of shearing direction, in order to recover the isotropic starting state.
Although the reversibility of this simulation is extremely good, it is not perfectly reversible even if round-off errors are disregarded. The stress calculations in the forward-and reversed-time frames differ slightly, as discussed in Section 3.1. In addition, a proportional controller adjusts the positions of the periodic boundaries based on the average stresses on the assembly and a user-defined gain. The disparity between the forward-and reversed-time simulations depends on the proportional gain: a larger gain reduces the disparity but instability would occur if too large a value were chosen.

SHEAR FORCES
Consider now the somewhat more realistic scenario in which shear forces are computed in the forward-time frame using Equations (7),(8),(9), (10). Energy conservation is assumed, necessitating an effectively infinite interparticle friction coefficient to allow the shear forces to be investigated independently of Coulomb friction. The three-step shear force remapping given by Equations (7),(8),(9) is time-reversible. So too is Equation (10) as long as particles remain in contact, that is, as long as interparticle contacts are not lost during the forward-time simulation. If f s c,t+Δt is known at the start of a time step, f s c,t can be recovered. However, this is never the case for all contacts in a practical DEM simulation: particles are in constant motion and the contact network evolves continually, even in very dense systems. 24 The irreversibility arises because once particles separate in the forward-time frame, the shear force is set to zero irrespective of its magnitude prior to the loss of contact. Upon reversal, the magnitude of that shear force is unknowable. As a simple demonstration, consider the two-particle collision shown in Section 3.2.1 with the addition of an infinite interparticle friction coefficient. Figure 5 shows the magnitudes of normal and shear force in the forward-time frame. The shear force reaches a maximum magnitude of 66 mN over the course of the collision, decreasing to 35 mN at contact separation: a non-negligible shear force that could not be recovered in a reversed-time frame.
Another source of irreversibility arises when periodic boundary conditions are adopted. Particles are remapped to the deforming space whenever a periodic cell changes in size or shape. The strain rate of the deformation is considered when computing the relative velocities used in the shear force calculation. 21 The strain rates in the forward-and reversed-time frames can differ by one time step, for example, if set by a servo controller rather than being fixed. This is a similar error to the one discussed for the calculation of representative particle stresses (Equation (22)).

Friction
In a DEM simulation, friction requires the computation of shear forces which are, in general, irreversible for the reasons described in Section 4. However, frictional dissipation introduces an additional source of irreversibility. If at some time t + Δt we have |f s c,t+Δt | = |f n c,t+Δt |, a Coulomb friction criterion may have been applied. While f n c,t is easily obtained in a reversed-time frame, the first step in recovering f s c,t is reversing the imposition of the Coulomb friction limit. The unscaled shear force would have a magnitude of |f s c,t+Δt | where 1 ≤ ≤ max and max depends on U s c,t+Δt , the contact stiffness and the shear force remapping that has taken place. The irreversibility arises because can take a range of possible values; exactly which value is required to recover f s c,t cannot be known.

Damping
Local and viscous damping are widely used and available in many commercial and open-source DEM codes. The expression implemented in LAMMPS for local damping of the net force is given in Section 3.1 as Equation (21). A similar expression may also be applied to moment instead of force. Equation (21) is not time-reversible as can be seen from Equations (17), (18), (19), (20). One would need to know f d t+Δt and sgnv t+ Δt Viscous damping acts at a contact, multiplying the relative velocity between particles by a scalar viscous damping coefficient v . As for local damping, there is a mismatch between forces and velocities. In the forward-time frame, f c,t+Δt is concurrent with particle velocities at t + Δt 2 ; in the reversed-time frame, f c,t+Δt is concurrent with particle velocities at t + 3Δt 2 .
Neither local nor viscous damping is time-reversible using the conventional velocity-Verlet integration scheme shown in Section 2. One practical solution may lie in the augmentation of the Verlet integrator with a time-reversible velocity predictor. 25 This approach has been applied successfully to extended Lagrangian methods that include a velocity-dependent term in the equations of motion; the errors obtained for a Nosé-Hoover thermostat were acceptable. 25 Assessing these time-reversible velocity predictors for the reversibility of various damping approaches in DEM is a subject for future investigation.
The imposition of a drag force is a time-reversible means of dissipating energy from a DEM simulation that does not require any modifications of the integration scheme. However, it is only reversible under specific and highly restrictive conditions. A drag force can be incorporated into the forward-time algorithm presented in Section 3 by multiplying the particle velocity v t+ Δt 2 by d before Equation (14) is computed to yield v d t+ Δt 2 . d < 1 so that particle velocities reduce over time and kinetic energy is lost from the simulation. Equations (13), (14), (15) are unchanged by the imposition of damping during the first time step and precede Equation (25).
The reversed-time equivalent of this is the following: The computation of Equations (18), (19), (20) follow Equation (28) as in Section 3.
The results of this form of damping, and its limitations, are shown in Figure 6. In each of these three cases, a single sphere of diameter 1 cm and density 1000 kg m −3 was dropped under gravity (g= −9.81 m s −2 ) onto a horizontal plane. The particle was allocated a shear modulus of 2 MPa and Poisson's ratio of 0.2 while the plane was assigned values of 1 GPa and 0.2, respectively. The time step was 10 ns and d was fixed at 0.99999. The three cases shown in Figure 6 differ in two respects: the initial height of the particle above the plane, and the duration of the simulation in the forward/reversed-time frame. Case A began with the particle 20 m from the plane and the total simulation duration was 50 ms. For case B, the equivalent numbers were 20 m and 80 ms; for case C, 300 m and 80 ms.
Case A is clearly time-reversible whereas cases B and C are not. Case B is run for longer than case A, allowing the particle to come fully to rest on the plane; some small oscillations are still present at 50 ms for case A although these are imperceptible on Figure 6. Once the particle's velocity reaches exactly zero, there is no way to recover a nonzero velocity or net force upon time reversal. Case C demonstrates a related limitation of this damping approach. The particle reaches its terminal velocity of gΔt 1− d = −0.00981 m s −1 before colliding with the plane. The expected velocity is recovered in the reversed-time frame until a constant terminal velocity has been reached. Thereafter it is impossible to know when the velocity becomes less than the terminal velocity in magnitude. The premature deviation of the forward-time and reversed-time plots on Figure 6C are due to the propagation of tiny numerical errors in the reversed-time frame.

Rotational resistance
Simulating nonspherical particles in DEM is significantly more computationally expensive than using spheres, regardless of the approach chosen to model the nonspherical particles. A cheaper alternative which has sometimes been adopted is to add an artificial moment term at each interparticle contact. These moments restrict rolling and/or twisting between particles and hence capture some of the effect of incorporating nonspherical shapes on the bulk behavior of the particle system. Many models have been proposed for rotational resistance. 26 The more physically justifiable models, categorised as Type C by Ai et al 26 limit the moments that can be applied, for example,. [27][28][29] The imposition of a limit introduces an irreversibility similar to the one described for Coulomb friction in Section 5.1.

Particle crushing/bond breakage
There are two broad approaches by which particle crushing is simulated in DEM: 30 either the fundamental particles (spheres, polyhedra, etc.) are replaced with smaller fundamental particles when some failure criterion has been met, or an agglomerate composed of unbreakable fundamental particles is assembled using bonds with finite strengths. Neither is time-reversible for obvious reasons. In the first case, one would need to replace small particles at specific instants during the reversed-time simulation with larger particles, regardless of the forces or stresses experienced by the small particles.
In the second case, fundamental particles would need to be connected by bonds as the reversed-time simulation proceeds. There are infinitely many ways of accomplishing both.

Other complications
There are numerous other complications of the basic DEM scheme which introduce time irreversibility. Some normal contact force models exist which cannot be reversed, for example, elastic-plastic models with distinct loading and unloading curves. 31 The coupling of DEM with another solver, for example, CFD, SPH, LBM, MPM, MBD, or FEM, is likely to introduce some time irreversibility; elucidating the nature of these irreversibilities is outside the scope of this paper.

EXEMPLAR APPLICATION: CRITICAL-STATE SOIL MECHANICS
The general irreversibility of the energy dissipation mechanisms presented in Section 5 severely restrict the practical applications of reversed-time DEM. For example, it would be interesting to introduce perturbations into the final state of a frictionless system and then reverse time to explore the achievable range of initial states. It would be useful to reverse time using a different loading protocol, resulting in two simulations which follow different paths but reach identical final states. However, the large oscillations that would arise in the absence of energy dissipation make these potential applications infeasible. Hence the main application of reversed-time DEM is its implications for existing theories of granular behavior. One potential application area is the packing and jamming of frictionless grains. This has been the subject of significant research attention among the physics community in the last 20 years. 32,33 Another application area, chosen here as an exemplar, is a less idealized one: critical-state soil mechanics (CSSM). This is a framework for describing the behavior of soil which was developed during the 1950s and 1960s. 34,35 One of the key principles of CSSM is that once a soil has been sheared to a sufficiently large strain, continued shearing no longer changes either the soil's volume or stress state: the soil attains a "critical state." This is shown schematically in Figure 7 for two samples of the same soil, differing only in terms of their initial void ratio (the volume of voids divided by the volume of solids). In soil mechanics terminology, the black line is characteristic of a dense soil which dilates prior to critical state, while the grey line is characteristic of a loose soil which contracts prior to critical state. Both samples become indistinguishable at critical state: not only do they have the same void ratio and stress state but DEM simulations have demonstrated that the samples will also have the same micro-scale measures such as coordination number and deviatoric fabric.
It is clear from Figure 7 that samples with very different initial states can reach the same critical state when subjected to the same loading protocol. This has been proved by countless laboratory experiments and DEM simulations. Consider one of the numerical samples on Figure 7 which has been sheared to a large strain, that is, is at the right-hand side of Figure 7. A reversed-time simulation must be able to recover both the loose and dense paths shown on the figure (and indeed the infinitely many other possible paths). The only way to recover a diverse range of initial states from a reversed-time DEM simulation is the inclusion of time irreversibilities such as friction in the simulation. In the absence of friction or other irreversibilities, there is a one-to-one relationship between the input and output states. In a frictionless, energy-conserving system, it would be impossible for samples differing only in terms of their initial void ratios to converge to the same critical state when following the same shearing protocol. The introduction of friction and shear forces vastly expand the range of attainable states and enable the dense/loose behavior seen in Figure 7. This once again highlights the importance of interparticle friction to CSSM. 36

CONCLUSIONS
For the first time, an algorithm has been presented for DEM simulations in a reversed-time frame. This algorithm avoids the requirement to apportion the net force on a particle to contacts, which cannot be done uniquely, by changing the order in which the equations are computed compared to the forward-time frame. This algorithm is restricted to an energy-conserving system of frictionless, spherical particles. Following its implementation in LAMMPS, the algorithm was verified using three cases: a two-particle collision, a strain-controlled oedometric compression, and a stress-controlled triaxial compression of a fcc assembly of monosized spheres. Time irreversibility is inevitable once shear forces have been incorporated into a simulation: whenever particles separate in the forward-time frame, the interparticle shear force is set to zero, so its magnitude before being zeroed cannot be known upon time reversal. In addition, energy dissipation mechanisms such as frictional sliding, local/viscous damping, rotational resistance, particle crushing or bond breakage are shown not to be time-reversible in general using a standard velocity-Verlet integration scheme.
The absence of a practical, time-reversible mechanism for energy dissipation limits the applications of reversed-time DEM. Its main value is likely to lie in extending theories of granular behaviour. The exemplar application considered in this paper is CSSM. In the absence of friction or other irreversibilities, there is a one-to-one relationship between the input and output states. Hence the convergence of initially dense and loose samples, sheared in the same manner, to the same critical state is possible only because of the presence of friction and shear forces within the soil system. This emphasises the importance of interparticle friction to CSSM.