Linear-frictional contact model for 3D discrete element simulations of granular systems

by   Matthew R. Kuhn, et al.

The linear-frictional contact model is the most commonly used contact mechanism for discrete element (DEM) simulations of granular materials. Linear springs with a frictional slider are used for modeling interactions in directions normal and tangential to the contact surface. Although the model is simple in two dimensions, its implementation in 3D faces certain subtle challenges, and the particle interactions that occur within a single time-step require careful modeling with a robust algorithm. The paper details a 3D algorithm that accounts for the changing direction of the tangential force within a time-step, the transition from elastic to slip behavior within a time-step, possible contact sliding during only part of a time-step, and twirling and rotation of the tangential force during a time-step. Without three of these adjustments, errors are introduced in the incremental stiffness of an assembly. Without the fourth adjustment, the resulting stress tensor is not only incorrect, it is no longer a tensor. The algorithm also computes the work increments during a time-step, both elastic and dissipative.


A Fully Implicit Method for Robust Frictional Contact Handling in Elastic Rods

Accurate frictional contact is critical in simulating the assembly of ro...

Dynamic Models of Planar Sliding

In this paper, we present a principled method to model gen- eral planar ...

An efficient diffusion generated motion method for wetting dynamics

By using the Onsager variational principle as an approximation tool, we ...

Traction chain networks: Insights beyond force chain networks for non-spherical particle systems

Force chain networks are generally applied in granular materials to gain...

An Efficient Model Order Reduction Scheme for Dynamic Contact in Linear Elasticity

The paper proposes an approach for the efficient model order reduction o...

Realistic Haptic Rendering of Interacting Deformable Objects in Virtual Environments

A new computer haptics algorithm to be used in general interactive manip...

Approximating intractable short ratemodel distribution with neural network

We propose an algorithm which predicts each subsequent time step relativ...

1 Introduction

The discrete element method (DEM) is a numerical technique for simulating the behavior of granular systems and investigating the grain-scale 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 time-steps, requiring the calculation of the grain-to-grain contact forces of every contact and at every time-step. A well-defined, 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 time-step

is equal to the two particles’ relative tangential movement, vector

, during the time-step 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


where and are the tangential forces at the start and end of the time-step, and is the unit vector normal to the contact plane at the end of the time-step. 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 three-dimensional (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 time-step. The drawbacks of conventional frictional models have been noticed by other authors 2, and although subtle, the evolution of within a time-step 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):

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

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

  3. 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 time-step. We correct for this directional change, which will usually produce a force increment that is not aligned with the movement .

  4. The two particles can roll and twirl during a time-step, 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 time-step, we assume that the contact’s movements, and , measured from the start of are proportional and uniform within , so that they are parameterized as


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.

Figure 1:  Algorithm for computing the contact force and mechanical work.

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),


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 time-step, then we compute the value at which the particles first touch (line 11),


and we replace and with the movements that occur after the particles have touched,


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. (23), as in line 16.

If the particles were already touching at the start of , then the increment in Eq. (22

) 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.

Figure 2:  Contact tangent plane when the initial movement is elastic, followed by slip.

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. 41). 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 (41). 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,


where the squared quantity on the left is the radius of the dashed circle in Fig. 2. We solve for (lines 22–27), as


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):


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


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.

Figure 3:  Evolution of the tangential force and sliding direction: (a) yield surfaces and the initial and final directions of the tangential force, and (b) increment in movement tangent to the yield surface. Note that the angles are negative in the directions of this figure.

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 cross-product , where is the contact’s normal direction at the end of the time-step, 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 time-step (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


where proceeds from 0 to 1. Note that angle will be reduced during the time-step, as becomes more aligned with , which is consistent with the negative in the equation. In this equation, the normal force will change during time-step and, as such, is a function of (see Eqs. 1 and 4),


Equation (16) is a first-order separable differential equation on the domain , in which angle starts at and ends at . The solution is (lines 36–40)


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 look-up table and interpolate to determine

to find the term on the right, and we can create an inverse look-up table and interpolate for to determine . These master tables are efficiently reused with every contact and at every time-step. 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),


where we have applied the friction limit on magnitude , as in Eq. (23).

Figure 4 shows the effect of the movement magnitude on the final direction of the tangential force.

Figure 4:  Example showing the effect of the magnitude of the tangential movement on the final direction of the tangential force. The movement is applied orthogonal to the initial 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 time-step . 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 time-step , thus rotating the tangent plane. As such, the tangential force can drift from its tangent plane after a series of time-steps. 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)


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 micro-scale 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: mass-damping of a particle body, which is proportional to the particle’s velocity, and contact-damping, which is proportional to the contact velocity (see Eq. 5). If contact-damping 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 macro-scale 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 rigid-frictional. As another illustration, the macro-scale behavior of 2D disks and 3D spheres are qualitatively similar, and one might not even guess a simulation’s dimensionality by looking at the stress-strain 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 stress-probe 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 non-convex sphere-clusters contained within periodic boundaries. The dense initial particle arrangement was isotropic with an initial porosity 0.363 (void ratio of 0.570), an inter-particle 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 stress-probes: one series with the conventional unimproved linear-frictional 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 strain-probes 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 elasto-plasticity. 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 probe-point 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 probe-points 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 elasto-plasticity 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 flow-direction for each stress-direction, 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 5:  Plots of incremental plastic strains from stress-probes in the Rendulic plane of volumetric and deviatoric stress and strain: (a) technique for plotting the plastic increment of a single probe; (b) locus of probe points, as predicted by elasto-plastic theory; and (c) locus of probe points for two the linear–frictional 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 stress-probe in which the strain increment had a magnitude of that was achieved with 30 DEM time-steps. 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’ work5, indicating a failure of elasto-plastic 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 elasto-plasticity 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 time-steps, 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 co-rotations 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 (pre-rotated) frame, designated as stress , by applying the tensor transformation


where the rotation tensor effects a full counter-rotation 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 counter-rotated) 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.

Figure 6:  Stress tensors computed for assembly of 10,648 particles that has been loaded with constant- triaxial compression in the direction. The units are kPa. Subfigure (a) gives the original tensor . The three other tensors have been restored to the original frame after a 90 rotation about the axis: (b) stress when neither the projection nor twirling adjustments are made during the 90 rotation, (c) stress when only the projection adjustments are made during the 90 rotation, and (d) stress when both the projection nor 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 linear-frictional 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 time-step 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 stress-probe 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) time-steps, 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 time-step. 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.


  • 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 non-spherical particles for granular and multi-body systems. Singapore: John Wiley & Sons; 2014.
  • 3 Lin X., Ng T.-T.. A three-dimensional discrete element model using arrays of ellipsoids. Géotechnique. 1997;47(2):319–329.
  • 4 Vu-Quoc L., Zhang X., Walton O. R.. A 3-D discrete-element method for dry granular flows of ellipsoidal particles. Comput. Methods Appl. Mech. Eng.. 2000;187(3):483–528.
  • 5 Kuhn M. R., Daouadji A.. Multi-diirectional behavior of granular materials and its relation to incremental elastoplasticity. Int. J. Solids Struct.. 2018;152-153:305–323.