1 Introduction
The discrete element method (DEM) is a numerical technique for simulating the behavior of granular systems and investigating the grainscale mechanics of these systems ^{1}. The method uses an explicit time integration to update the position and rotation for each particle and at each time in a series of timesteps, requiring the calculation of the graintograin contact forces of every contact and at every timestep. A welldefined, precise, and robust relationship between contact movement and contact force is essential for DEM codes, and the most common movement–force relationship, by far, is the linear–frictional contact. With this model, the force components that are normal and tangential to a contact surface are separately computed. The normal (compressive) contact force between two particles, , at the time is simply the accumulated overlap of the particles’ idealized profiles multiplied by a normal contact stiffness . The change in tangential force, , that occurs during the timestep
is equal to the two particles’ relative tangential movement, vector
, during the timestep multiplied by the tangential stiffness , but the magnitude of the accumulated tangential force, , is limited to a friction coefficient multiplied by the normal force. These two rules are conventionally written as(1)  
(2)  
(3) 
where and are the tangential forces at the start and end of the timestep, and is the unit vector normal to the contact plane at the end of the timestep. The two rules may seem similar, but they are fundamentally different, with the tangential rule being the more complex. The normal force is simply a function of the accumulated overlap of the particles; whereas, the tangential rule is an incremental relationship. The normal force is readily expressed as a scalar; whereas, the accumulated tangential force and the tangential increment are vectors that lie in the contact’s tangential plane but may have different directions. These difficulties are particularly complex for threedimensional (3D) particles, and the complexity of the tangential relationship is the subject of this brief note, which is focused on the “” changes that occur within a single DEM timestep. The drawbacks of conventional frictional models have been noticed by other authors ^{2}, and although subtle, the evolution of within a timestep affects the full increment , an effect that is rarely included in DEM codes.
In this note, we consider four refinements of the tangential force in Eq. (2):

Fresh contact between two particles can occur within the span of a single timestep , and we correct for any small displacement that occurs before contact is established.

For an established contact that is not initially slipping but is slipping at the end of , slip does not occur throughout the full timestep, and we correct for the elastic displacement that precedes slip.

As a contact undergoes slip, the direction of the tangential force can differ from the direction of tangential movement , causing the direction of the tangential force to change within the timestep. We correct for this directional change, which will usually produce a force increment that is not aligned with the movement .

The two particles can roll and twirl during a timestep, and we correct for any rigid rotation of the particles and of the contact force between them.
In this short communication, we develop an algorithm, detailed in Fig. 1, to resolve these refinements. We also provide means of calculating the mechanical work done within a contact, both elastic and dissipated. Two examples are presented that illustrate the benefits of applying the four refinements.
2 Algorithm
Because the refinements listed above apply within the course of a timestep, we assume that the contact’s movements, and , measured from the start of are proportional and uniform within , so that they are parameterized as
(4) 
where proceeds from 0 to 1 during the full step . Note that the normal movement is a scalar; whereas, tangential movement is a vector, although its bold font might be indistinct. With the assumption of Eq. (4), the paper describes an algorithm that incorporates the four refinements listed above. The algorithm, summarized in Fig. 1, requires the following input information (line 1): the particles’ overlaps at times and , the two particles’ movements during , the contact’s stiffness and its friction coefficient, and other secondary information, including a small that gives the machine precision.
With DEM, one computes the overlap of the particles from the geometric descriptions of their shapes, and this overlap at is used to find the normal force , using Eq. (1) (lines 7 and 20). If the contact overlap is negative, the contact force is, of course, zero (lines 4 and 5). If the particles are touching, then we find the normal and tangential movements in Eq. (4) from the full contact movement of particle relative to particle at their contact (lines 8 and 9),
(5)  
(6) 
where the , , , and are the particles’ translations and rotations, and are the contact vectors that connect the centers of and with the contact point, and is the contact’s unit normal vector at the end of the time step, directed outward from . Note that a positive incremental overlap produces a compressional normal force.
If the particles are touching at the end of but were not touching at the start of the timestep, then we compute the value at which the particles first touch (line 11),
(7) 
and we replace and with the movements that occur after the particles have touched,
(8) 
where the arrow “” represents substitution (assignment) when coding the algorithm for contact force (line 12). If the contact is touching at both the start and end of , then in Eq. (8). The values of and are then used with Eq. (2) to compute a tangential force increment and the normal force (lines 13 and 14). If the particles were not touching at the start of , then the full tangential force is equal to the increment (line 14). The tangential force must also be checked to satisfy the friction limit of Eq. (2_{3}), as in line 16.
If the particles were already touching at the start of , then the increment in Eq. (2_{2}
) is merely an estimate of the actual increment, as described now. Having established the instant of contact, we then compute the tangential contact force. If the magnitude
at exceeds the friction limit , then the contact will slide and we must find the instant at which sliding commences. The situation is illustrated in Fig. 2, which represents a contact’s tangent plane, such that the vectors of tangential force and movement lie within this plane.At the start of (i.e., when ), the limiting magnitude of the tangential force is , but this limit will change during the course of , as proceeds from 0 to 1, due to a concurrent change in the normal force (see Eq. 4_{1}). The changing friction limit is represented by circular yield surfaces in the figure. The tangential force also changes during , and at some point , the tangential force will reach the friction limit and sliding will commence (represented by the dashed circle in Fig. 2). Within the interval , the behavior is elastic, and the tangential force at time is given by , as in Eq. (2). The corresponding friction limit at is , as in Eqs. (1) and (4_{1}). To solve for at the instant the tangential force reaches the friction limit, we equate the scalar magnitude and the scalar limit by applying the law of cosines,
(9) 
where the squared quantity on the left is the radius of the dashed circle in Fig. 2. We solve for (lines 22–27), as
(10)  
(11)  
(12)  
(13) 
Note that if the contact has already reached the friction limit at the start of , then and is the larger of and (line 25). When , we apply L’Hôspital’s to find (line 29):
(14) 
In the algorithm of Fig. 1, we use a small input parameter to test the proximity of and to zero (lines 25 and 26).
Sliding commences at time , when the normal and tangential forces, and , are
(15) 
and for the next set of calculations, we reset to zero, and replace with , with , with , and with (line 30).
Having established the start of sliding, we proceed to find the tangential force at the end of . At the start of sliding, the tangential force has the unit direction . This direction will likely differ from the direction of the tangential movement during . If so, the final direction of the tangential force at time will not coincide with either the initial direction , the direction of the movement , or the direction of the estimated tangential force , as given in Eq. (2). Rather, sliding will begin in the initial direction of the tangential force, , but this sliding will occur concurrently with the elastic movement that is orthogonal to the slip direction, and this elastic movement will alter the force direction during the increment . The final direction of the tangential force, , will be intermediate between the directions and .
This situation is illustrated in Fig. 3, which represents the tangential plane of the contact.
Recall that we have reset to 0 and have reset , , and to their values at the start of sliding, . The tangential plane of the figure contains the tangential force at the start of sliding (the in Eq. 15); the tangential movement ; and the final tangential force . The unit vectors and lie within this plane, with aligned with the initial force , and is aligned with . We also define a unit vector that is perpendicular to the direction of tangential movement (the vector can be computed with the crossproduct , where is the contact’s normal direction at the end of the timestep, line 31). At the start of sliding, when , direction makes angle with the initial force direction , such that , as in line 32. Although the directions and do not change during the timestep (recall Eq. 4), the force orientation and its corresponding angle will change, with starting at and ending at , and we now describe the method for finding this final angle. The circle in Fig. 3a is a yield surface that defines the frictional limit of the tangential force. Any tangential contact movement that is tangent to this circle (that is, movement within the tangent plane that is perpendicular to the current tangential force direction, ) is elastic. A change in the angle by the increment corresponds to an increment of elastic movement perpendicular to , equal to the tangential force increment divided by the elastic stiffness (Fig. 3b). This elastic force increment also equals the radius of the yield surface, , multiplied by . Noting that direction makes angle with (Fig. 3b, where the are shown as negative in this figure), the magnitude of the full increment of movement during increment is given by
(16) 
where proceeds from 0 to 1. Note that angle will be reduced during the timestep, as becomes more aligned with , which is consistent with the negative in the equation. In this equation, the normal force will change during timestep and, as such, is a function of (see Eqs. 1 and 4),
(17) 
Equation (16) is a firstorder separable differential equation on the domain , in which angle starts at and ends at . The solution is (lines 36–40)
(18)  
(19) 
where the parameters and are as follows: and . Note that the final term on the right, , is replaced with when (line 38). When , the direction of movement is aligned with the initial force , and (line 34). The that is computed with Eq. (18) will always be positive, so we must intervene so that its sign conforms with that of , as in Eq. (19).
Solving with the nonlinear Eq. (18
) can be computationally taxing, and rather than using Newton’s method (or another iterative approach), we can simply create a lookup table and interpolate to determine
to find the term on the right, and we can create an inverse lookup table and interpolate for to determine . These master tables are efficiently reused with every contact and at every timestep. Once is computed, we find the direction of the final tangential force, and its full value is the vector sum of its and components (lines 42–43),(20)  
(21) 
where we have applied the friction limit on magnitude , as in Eq. (2_{3}).
Figure 4 shows the effect of the movement magnitude on the final direction of the tangential force.
In this example, the movement vector is applied orthogonal to the initial tangential force , with , and the normal force is assumed to remain constant during the timestep . Under these conditions, the tangential force will rotate as the movement increases, to become more closely aligned with the direction of . When the movement vector is very small, the direction of is barely altered, and the final angle remains near . Conversely, with large movements, the force rotates to conform with the movement direction, with , approaching . The figure also shows the commonly used approximation of , in which the increment is imply added to the initial force and then projected onto the final yield surface to find the final force, ignoring the progressive change in the direction of slip during , as in Eq. (16). Although the approximation works well when the increment is small compared with the current normal force (i.e., small compared with the quotient ), large errors occur when the force increment is larger than 20% of the initial tangential force.
Having found the new tangential force in the manner of Eqs. (20) and (21), we must remember that the two particles are moving and rotating during timestep , thus rotating the tangent plane. As such, the tangential force can drift from its tangent plane after a series of timesteps. The fourth (and final) refinement of the force calculation accounts for this directional drift and for twirling of the tangential force within the tangent plane. In DEM codes, the drift is usually prevented by projecting the new force onto the new tangent plane ^{3, 4} (line 48). Less common, however, is an adjustment that must be applied within the tangent plane. If the two particles rotate (twirl) about their contact normal as a rigid pair, then the tangential force will also rotate, even in the absence of a relative tangential movement between the particles. This rotation of the tangential force is consistent with the principle of objectivity, as no rotation would be seen by an observer who rotates in unison with the twirling pair; whereas, a stationary observer would see an equal rotation of the two particles and of the tangential force. The projected adjustment and the the rotated (twirled) adjustment are computed successively as (lines 48 and 49)
(22)  
(23) 
where and are the particles’ incremental rotation vectors (as in Eq. 5), and the twirling is simply the average of the two rotations about the contact normal .
As two particles interact, elastic energy is stored or released from their contact, and energy can also be dissipated in friction. We are often interested in the microscale destination of the mechanical work that is delivered to a contact. In Fig. 1, we compute the total, accumulated work done by the normal force and the work increment done by the tangential force. The normal work is entirely elastic, and its total stored energy can be computed directly from the final force, as , in line 51. The tangential work is computed as an increment during and is the sum of elastic (reversible) and dissipated (irreversible) parts, and . The reversible tangential increment in line 52 simply depends on the initial and final forces, but the increment of total tangential work depends upon the tangential force–movement path, which is quite complex for a sliding contact (this complexity is partly expressed in Eq. 16, which can be used to find the relationship between movement and the evolving tangential force ). For tangential movements that precede slip, the work increment is simply equal to the elastic increment (lines 18 and 28). Once sliding commences, we approximate the tangential work as the inner product of the movement and the average of the tangential forces at the start of sliding and at the end of (line 45). The increment of dissipation is equal to the total tangential work minus the elastic increment (lines 52–53).
Viscous dissipation is typically used as a means of stabilizing DEM simulations, and Cundall and Strack ^{1} describe two forms of such damping: massdamping of a particle body, which is proportional to the particle’s velocity, and contactdamping, which is proportional to the contact velocity (see Eq. 5). If contactdamping is used, it is effected by simply adding a damping force to the contact force in line 50 of Fig. 1, where is the damping constant.
3 Examples
Two examples are presented: the first illustrates the effect of making the first three adjustments to the conventional linear–frictional algorithm; whereas, the second example illustrates the effects of including the rolling and twirling of tangential forces (see the four refinements outlined in the Introduction). We should emphasize, however, that these improvements will not appreciably affect certain types of DEM results, but can become significant when DEM simulations are used for studying other elements of granular behavior. The macroscale behavior over large spans of strain is fairly insensitive to the contact details: for example, similar strength results are obtained from both DEM and Contact Dynamics (CD) simulations at large strains. The two simulation methods are quite different, and the latter altogether neglects contact stiffness and treats the contacts as rigidfrictional. As another illustration, the macroscale behavior of 2D disks and 3D spheres are qualitatively similar, and one might not even guess a simulation’s dimensionality by looking at the stressstrain response across large spans of strain. The effect of the first three linear–frictional refinements will become more pronounced for small spans of strain, as with stressprobe studies, and under conditions in which many contacts are either sliding or are undergoing combined sliding and elastic movement. These conditions are explored in the next paragraph, and the effect of the fourth refinement is addressed further below.
A straightforward simulation of triaxial compression was conducted on a 3D assembly of 10,648 smooth nonconvex sphereclusters contained within periodic boundaries. The dense initial particle arrangement was isotropic with an initial porosity 0.363 (void ratio of 0.570), an interparticle friction coefficient , and equal normal and tangential contact stiffnesses, . The preliminary stage of loading was drained isobaric (constant) triaxial compression in the direction, using the refined model of this paper. At peak strength, the ratio of deviator stress and mean stress was about 1.7, but our example was conducted at a smaller strain, when was about 0.8, where we suspended the loading and conducted two series of stressprobes: one series with the conventional unimproved linearfrictional model and the other series with the model in the paper. At this smaller strain, the material was in the early stage of strain hardening and about 11% of the contacts were either sliding or with a mobilized friction within 0.001% of .
For each of the two linear–frictional algorithms, over 70 axisymmetric probes were conducted with a technique previously described by the authors ^{5}, such that both elastic and plastic strain increments were measured for a suite of incremental loading directions (more precisely, reversible and irreversible increments). All probes were axisymmetric and conducted within the  Rendulic plane, with strainprobes of magnitude . We focused upon the plastic increments for each probe, which we processed and plotted in a manner that permits direct evaluation of the yield surface, flow direction, and plastic modulus as well as conformance with conventional elastoplasticity. The plotting technique is more extensively described elsewhere by the authors ^{5} and is illustrated in Fig. 5. The result of a single probe is represented by a single point in the Rendulic plane of volumetric and deviatoric strains and stresses, with the axes normalized to a common orthonormal basis (volumetric strain , and deviatoric strain , with ). Each probepoint is located along a direction that is aligned with the stress increment, and the radial distance from the origin is equal to the magnitude of the resulting plastic strain increment divided by the magnitude of the stress increment (Fig. 5a). Although this form of plotting is unusual, it presents information about both the strain and stress increments, and it provides the following information in a compact form (Fig. 5b): (i) the locus of probepoints will form a single circle in the Rendulic plane that passes through the origin (and a sphere in the full –– space), provided that the material conforms with conventional elastoplasticity theory; (ii) the diameter of the circle is the inverse of the plastic modulus; and (iii) the orientation of the circle equals the direction of the yield surface. The plot can also include the flowdirection for each stressdirection, although we did not use this feature with the simplified plots of the paper. Here we use the plot to compare the hardening moduli from the series of DEM simulations for the two linear–frictional contact algorithms.
Figure 5c shows the results for both the conventional linear–frictional algorithm and that of the paper. Each locus includes over 70 points, and each point is the result of a DEM stressprobe in which the strain increment had a magnitude of that was achieved with 30 DEM timesteps. Each locus is shown to be a circle within the Rendulic plane, although this is not the case when plotted in the deviatoric plane (see the authors’ work^{5}, indicating a failure of elastoplastic principles when applied to granular materials). The tilt of the larger circle indicates the yield direction, and the presence of a small complementary circle simply means that, contrary to elastoplasticity principles, plastic strains occur in the “unloading” direction (within the presumed elastic region), a result that has been demonstrated by the authors ^{5} for strains as small as .
The results in the figure show that the hardening modulus (inverse of the larger circle’s diameter) for the conventional linear–frictional algorithm has an error of about 5%, when compared with the paper’s algorithm. Although small, the error will increase if probes are done with fewer than 30 timesteps, and the error will also increase as loading proceeds toward the peak stress, bringing more contacts to the friction limit.
We now consider the effect of the fourth adjustment that is listed in the Introduction and accomplished with Eqs. (22–(23). After computing a new tangential force that results from the contact movements and , the tangential force must be adjusted in two ways: (a) by projecting each contact’s force onto its tangent plane and (b) by twirling the force about the contact normal as a result of any corotations of the two particles. Only the projected force is typically computed ^{3}. The importance of both adjustments, however, is illustrated in a simple example. As with the previous example, the assembly of 10,648 particles was loaded in constant triaxial compression to a of about 0.8 by compressing the assembly in the direction. The simulation was then stopped at this stress , and the entire assembly was rotated as a rigid body about the axis in a sequence of 1000 rotation increments to an cumulative rotation of 90. Because this rigid rotation produces no relative movements among the contacts, neither the nor movements, the tangential contact forces do not change, even though the forces must rotate with the entire assembly (stated differently, an observer who rotates with the assembly would see no change in the forces, although a stationary observer would watch the forces rotating). This rotation is not considered in Eqs. (1)–(21), but is realized with the adjustments of Eqs. (22) and (23). The stress after the 90 rotation, , was computed from the contact data with the usual Love–Weber (Piola) summation of dyads of branch vectors and contact forces. To illustrate the effects of the projection and twirling adjustments, we restored the stress tensor to the original (prerotated) frame, designated as stress , by applying the tensor transformation
(24) 
where the rotation tensor effects a full counterrotation of 90. Standard DEM codes will correctly account for rotations of the branch vectors as well as rotations of the normal forces, but a proper algorithm must also rotate the tangential forces. If correct, the rotated (and then counterrotated) stress will coincide with the original stress tensor . Figure 6 shows four stress tensors: (a) the original stress before the 90 rotation, (b) the stress when neither the projection nor twirling adjustments are made during the 90 rotation, (c) the stress when only the projection adjustments are made during the 90 rotation, and (d) the stress when both the projection and twirling adjustments are made.
The results show that unless both adjustments are made, the original stress is not recovered, and a substantial error results. Indeed, borrowing a phrase of Wolfgang Pauli, the stress tensors of Figs. 6b and 6c, which do not include the proper adjustments, “are not even wrong”: they are not even tensors.
4 Discussion
In the paper, we document an algorithm for computing the contact force of a linearfrictional contact, which is the simplest and most common model used in DEM simulations. The algorithm resolves several subtle difficulties that arise when using a finite timestep but are typically not addressed in conventional codes. The first three difficulties involve incremental changes in the tangential force and lead to differences between the refined approach of the paper and the more conventional approach of simply projecting an approximate final force onto the slip (yield) surface. A stressprobe example shows that modest errors result in an assembly’s incremental stiffness when the three adjustments are not made. The error can certainly be reduced by using smaller (and more) timesteps, but at the expense of using larger probes and additional computation time. Unless it is properly resolved, the fourth difficulty produces significant accumulating errors in the tangential forces, errors that are unaffected by the size of the timestep. Without applying the fourth adjustment, the resulting stress tensor is not only incorrect, it is no longer a tensor. Although these are practical and compelling reasons for using an improved algorithm, a more thoughtful argument is that a numerical simulation, which is intended to model or understand a physical process, should be faithful to the underlying physics of the process.
References
 1 Cundall P. A., Strack O. D. L.. A discrete numerical model for granular assemblies. Géotechnique. 1979;29(1):47–65.
 2 Matuttis H.G., Chen J.. Understanding the discrete element method: simulation of nonspherical particles for granular and multibody systems. Singapore: John Wiley & Sons; 2014.
 3 Lin X., Ng T.T.. A threedimensional discrete element model using arrays of ellipsoids. Géotechnique. 1997;47(2):319–329.
 4 VuQuoc L., Zhang X., Walton O. R.. A 3D discreteelement method for dry granular flows of ellipsoidal particles. Comput. Methods Appl. Mech. Eng.. 2000;187(3):483–528.
 5 Kuhn M. R., Daouadji A.. Multidiirectional behavior of granular materials and its relation to incremental elastoplasticity. Int. J. Solids Struct.. 2018;152153:305–323.