Implementation of the Jäger contact model
for discrete element simulations
Abstract
In threedimensional discrete element (DEM) simulations, the particle motions within a granular assembly can produce bewildering sequences of movements at the contacts between particle pairs. With frictional contacts, the relationship between contact movement and force is nonlinear and pathdependent, requiring an efficient means of computing the forces and storing their histories. By cleverly applying the principles of Cattaneo, Mindlin, and Deresiewicz, Jäger (2005) developed an efficient approach for computing the full threedimensional force between identical elastic spheres that have undergone difficult movement sequences (J. Jäger, New Solutions in Contact Mechanics, WIT Press, Southampton, U.K.). The paper presents a complete Jäger algorithm that can be incorporated into DEM codes. The paper also describes three special provisions for DEM simulations: (1) a method for handling particle pairs that undergo complex tumbling and twirling motions in threedimensions; (2) a compact data structure for storing the loading history of the many contacts in a large assembly; and (3) an approximation of the Jäger algorithm that reduces memory demand. The algorithm addresses contact translations between elastic spheres having identical properties, but it does not resolve the tractions produced by twisting or rolling motions. A performance test demonstrates that the algorithm can be applied in a DEM code with modest increases in computation time but with more substantial increases in required storage.
16000010 \runningheadsM. R. Kuhn Jäger Contact for DEM Simulations \corraddrUniversity of Portland, 5000 N. Willamette Blvd., Portland, OR, 97203 USA, Tel. (1)5039437361, Fax (1)5039437316 \cgsnNational Science FoundationNEESR936408 \noreceived \norevised \noaccepted
1 Introduction
The discrete element (DEM) method is a computational approach for simulating granular materials by tracking the interactions of discrete, individual grains at their contacts. Realistic simulations require a contact displacementforce model that is faithful to the physical grains that are being modeled. In early implementations, simple linear springs were used to model the normal and tangential forces between particles, with the tangential force being limited by an abrupt frictional threshold [1]. Simulations that employ such simple models can yield qualitative agreement with some granular phenomena, but they are inadequate for the quantitative analysis of many situations. For example, linear contact models are incapable of accounting for the observed dependence of bulk material stiffness and wave speeds upon the confining pressure, which is an essential element in geotechnical problems. Linear models also cannot account for the progressive and hysteretic loss of stiffness in materials that are cyclically sheared with small strain amplitudes, as in liquefaction problems.
To account for these and other effects, Hertz theory is now available in many DEM codes to model the normal forces that develop between smooth elastic spheres that are pressed together [2]. The corresponding tangential shearing tractions that arise when particles shift in transverse, tangential directions were first solved by Cattaneo and Mindlin, who reported the progressive nature of frictional slip within the small circular contact area between two spheres [3, 4]. Such shearing tractions are historydependent, and the systematic analysis of loading histories usually begins with the work of Mindlin and Deresiewicz [5] who cataloged a canon of eleven basic histories and derived analytic expressions for the resulting tractions. This significant work was limited in two respects: (1) tangential loading and unloading were assumed to occur in a single tangential direction with no provision for transverse tangential movements within the contact plane, and (2) the loading histories were limited to a few stages, with each additional stage bringing further analytical difficulties.
DEM simulations present far more complex situations. Figure 1 shows the movement of a typical contact in a DEM simulation of slow (quasistatic) monotonic biaxial compression of a dense material. This particular contact, one of several thousand, became engaged after the peak stress state had been reached, and it remained engaged through an assembly strain of about 5%, during which the normal movement and the two components of tangential movement, and , were tracked. The contact displays a bewildering sequence of loading and unloading in all directions, and this for a monotonic loading. Although rapid collisional flows produce briefer and simpler particle interactions, the behavior in Figure 1 is typical of slow flows of dense materials in which the particles are in persistent contact across extended periods of bulk deformation. Cyclic loading of dense assemblies will produce even more varied motions at the contacts.
Seridi and Dobry [6] approached the first limitation by analyzing an additional loading history in which tangential motion in one direction was abruptly followed by an infinitesimal tangential movement in the perpendicular direction. Placed in a continuum elastoplasticity context, they found that the result resembled a loading probe tangent to a yield surface. With this insight, they and their colleagues then resolved both of the MindlinDeresiewicz limitations by developing a general algorithm for the threedimensional contact displacementforce relation that was analogous to incremental elastoplasticity with kinematic hardening [7]. The algorithm approximates the force by using small increments of movement and by tracking the loading history as a sequence of yield cones. Because of its computational demands, the method has not been widely adopted.
More recently, VuQuoc and Zhang [8, 9] developed another tangential forcedisplacement algorithm that uses an incrementalstiffness form of the MindlinDeresiewicz model. They derived incremental stiffnesses for four possible loading cases: combinations of tangential and normal force increments that are each either increasing or decreasing. The hysteretic behavior is captured by storing the tangential forces at which reversals in the loading direction occur (termed “turning points”). They also employed an improved form of a partiallylatching spring system for the normal forces [10, 8]. This system models plastic deformation in the normal direction. Zhang and VuQuoc [11] demonstrated their contact by modeling rapid, collisional flows of soybeans. Most DEM codes now use a Hertzian normal force combined with some incrementalstiffness form of the MindlinDeresiewicz model (for example, [12, 13, 14]).
In a remarkable series of papers and book, Jäger developed an elegant approach to the CattaneoMindlinDeresiewicz problem, one that permits large and arbitrary movements to be analyzed exactly and in whole (notably [15, 16, 17, 18, 19, 20]). The method expresses an otherwise complex distribution of shearing tractions as a superposition of simple CattaneoMindlin functions which provide a compact means of chronicling the essential elements of quite convoluted loading histories. The author suggests that this method be named the Jäger contact or the Jäger algorithm. The paper presents a complete coding of the algorithm that is outlined in the Jäger text ([20], pp. 129–130). The coding is intended for displacementdriven DEM simulations and includes additions that assure its computational dependability for arbitrary displacement histories.
The Jäger algorithm computes the contact force for each step in an arbitrary sequence of normal and tangential movements, such as the movements in a series of DEM time steps. The algorithm has several distinguishing characteristics:

Each normaltangential movement can be of arbitrary size without the need for incremental stiffnesses and without dividing a large movement into smaller, incremental submovements. In this sense, the Jäger algorithm reproduces the closedform solutions of Mindlin and Deresiewicz (such as equations 7 and 14 in [5]) that give the accumulated contact displacement in terms of the full contact force. (MindlinDeresiewicz also present instantaneous, incremental stiffnesses and compliances, but these are not part of the Jäger method.) The method only requires that the total normal movement is small when compared with the particle radius and that each movement is monotonic, such that motion advances in a continuous and proportional manner within each movement.

The algorithm includes a rigorous means of identifying load reversals (turningpoints) in a full threedimensional setting, and it incorporates the effect of such reversals upon subsequent loading steps.

As shown by Mindlin and Deresiewicz [5], certain combinations of tangential and normal movements will suppress slip within the circular contact area between elastic spheres. Jäger refers to such movements as producing a stick condition, but such movements are also called elastic or nonsimple [8]. Even though slip is suppressed, such movements will alter the distribution of contact tractions and, hence, will affect the onset and extent of slip in future movements. Mindlin and Deresiewicz [5] do not address nonsimple loading sequences in which a general elastic movement is followed by a nonelastic movement. Jäger derived a solution to this problem by following the tractions in the current and past slip zones of a contact area and the manner in which the CattaneoMindlin tractions must be superposed (see [17, 18] for succinct derivations).
Because both elastic and nonelastic movements can affect future slip, the Jäger algorithm thoroughly chronicles any changes in the directions of either elastic or nonelastic movements, including those at turningpoints. The algorithm then accounts for the effect that past elastic and nonelastic movements will have upon the current and future forces. In short, the Jäger algorithm exactly computes the contact force for nonsimple histories, and it does so without the use of incremental stiffnesses or compliances.
Besides presenting the coding for the Jäger algorithm, the paper also describes three special provisions that apply to DEM simulations: (1) a method that is necessary for applying the Jäger algorithm to particle pairs that undergo the twirling and tumbling motions that can be expected in threedimensional DEM simulations, i.e. motions that will also twirl and rotate the entire contact force; (2) an approximation of the Jäger contact that reduces memory demand; and (3) a compact data structure for storing the loading histories of the many contacts in a large assembly.
The algorithm (and the paper) is not without limitations, as it only concerns contact translations and does not account for tractions produced by the twisting (torsional) rotations of two particles about their contact normal direction, nor does it account for tractions produced by the rolling between particles (rolling friction). The paper’s underlying material model is one of isotropic linear elasticity for two spheres having equal elastic properties, an assumption carried over from the original works of Hertz, Cattaneo, Mindlin, and Deresiewicz. The method does not account for contact adhesion or for plasticity, viscoelasticity, or fracture and breaking of the spheres. Investigators have recently modeled the contact problem as an elastoplastic process [10, 21]. VuQuoc and Zhang modeled the normal forcedisplacement relationship of elastoplastic spheres [22, 23] and, more recently, developed an algorithm for the tangential and normal interactions of elastoplastic spheres [24, 25]. Experimental validation of this advanced normal force model is provided in [26] for polymer spheres. The paper assumes that the particles are spherical at their contact, although Jäger’s theory encompasses other shapes [20].
The next section describes the Jäger contact and provides detailed pseudocode of its algorithm. Section 3 supplies three additional provisions (listed above) that were not part of Jäger’s original method but are required for an efficient DEM implementation. The algorithm’s performance in a large DEM simulation is described in Section 4, which also includes a simple example that demonstrates the exact correspondence of the Jäger and MindlinDeresiewicz results and compares the Jäger solution with an incrementalstiffness solution.
2 Jäger Algorithm
We consider two identical isotropicelastic spheres that undergo an arbitrary sequence of finite, perhaps large, translations. Contact forces and are the compressive normal force and the 2vector tangential force. Displacements and are the cumulative normal approach (overlap) and the 2vector tangential shift of the particles’ centers. Figures 2 and 3 present the Jäger algorithm that is outlined in [16] and [20], §7.2.
The two figures also include the author’s additions that improve computational dependability and efficiency. The figures present a procedure (function) that would be called within each time step and for each contact within a DEM assembly.
For the given input (line 1) the algorithm returns the contact’s normal and 2vector tangential forces, and (line 2). The input includes the contact halfoverlap (indentation) ; the 2vector tangential halfshift movement ; the friction coefficient ; the modulus , which depends on the particles’ shear modulus and Poisson ratio ; the ratio of elastic normal and tangential stiffnesses, ; the particles’ radii of curvature, and ; the previous tangential force ; a small tolerance parameter ; a pointer to the top of an equivalent load history stack, described in Section 3.3; and a parameter to reduce memory demand (Section 3.1). Lists that store the equivalent load history are included as both input and output (line 3). The contents of these lists (i.e., , , , and ) are described later.
The Hertz normal force only depends on the indentation [2]:
(1) 
where the modulus (line 13). For spheres of different radii, an approximate average radius is adopted on line 8 [27].
In the work of Mindlin and Deresiewicz [5], most equations for tractions, compliances, and displacements refer to the radii of various contact and stick areas. For example, the full contact area has radius
(2) 
and they gave the radii of stick areas the symbols , , etc. Although these radii could be used within the Jäger algorithm, it is more convenient to use and as simpler proxies. In this regard, equations (1) and (2) can be used to shift between the results in the paper and those of Mindlin and Deresiewicz.
Once the particles touch, the tangential force will depend upon the DEM history of the finite movements (translations) between time steps, the and , that lead to the contact’s current position . The current force is produced by a force increment relative to the previous force (i.e., from the previous time step):
(3) 
(line 71). The tangential increment is entirely elastic, producing no slip, when the current movement satisfies the following two conditions: and (line 15, as demonstrated in [5], §14, 15, and 18). When either condition is violated, the current movement will produce frictional slip within an outer annular ring of the circular contact area between the spheres [5].
The current movement will not only affect the current increment of tangential force but can affect future increments as well. A journal of the contact’s past loading, in the form of an equivalent load history is maintained, so that the current and future loads will be reconciled with this history. The contact’s equivalent (star) load history is a compact recording of an equivalent sequence of loading steps that would lead to the previous contact force and the cumulative displacement . The equivalent load history is stored as lists of the normal indentations and the corresponding normal forces — the lists , and , — along with a list of 2vector directions of tangential force, , . Each  pair designates the apex of a yield cone in displacementspace and forcespace. In another sense, the and are proxies for the radii of past slip areas, as in equation (2) and in reference [5]. This “” data gives information about the points of load reversals (turningpoints) as well as information from past elastic movements that will affect the onset and the extent of future frictional slip. Upon entering the procedure, the history’s index will range from to . The zeroth point is permanently initialized to the unloaded condition: and . Within the procedure, the lists’ lengths can be increased by 1, remain the same, or be reduced to as small as 2 (lines 52 and 70). Before leaving the procedure, the most recently computed values of , , and are appended to the equivalent history (lines 70–71). Because these lists are more compactly stored as linked lists, the indices 0,1, are, in reality, pointers to locations within such linked lists (Section 3.3).
The concept of an equivalent load history is illustrated in Figure 4. In this simple twodimensional example, the tangential movements are colinear, so only single components of the twovectors and are relevant. The first part, Figure 4a, shows an example sequence of nine movements between two spheres having friction coefficient and , where .
The particles first touch at stage “0” and then undergo a sequence of nine proportional (straight) normaltangential movements: , ,, , These successive movements produce the corresponding sequence of forces in the second part, Figure 4b, labeled as the points 1, 2, , 9. The solid line in Figure 4b gives the final equivalent load history at the end of step 9. This equivalent history (path 0–1–1a––7a–9) would produce the same final displacement (, ), the same contact force (, ), and the same tractions as the original nine steps (0–1–2––8–9). The dashed lines in Figure4b trace the intermediate equivalent histories that would lead to the cumulative intermediate forces , , etc. These intermediate (dashed) segments were eventually eliminated within the equivalent history that leads to the final force . For example, the sequence 0–1–1a–5 will produce the same cumulative displacement (, ) as the sequence 0–1–2–3–4–5, so only the truncated (solid) sequence is stored as the equivalent load history at step 5. Although five of the nine movements are elastic (1, 2, 3, 6, and 7), four have slopes that violate the elastic condition and thus produce annular slips (4, 5, 8, and 9, shown as the darker lines in Figure 4a). The evolution of the equivalent load history is listed in Figure 4c, with all segments constrained to a slope no steeper than . The four slip movements result in equivalent path segments that have a slope magnitude : the segments 2a–4, 1a–5, 6a–8, and 7a–9. Along the final (solid) equivalent load history in Figure 4b, three of the segments are shown with heavier solid lines, as they are nonelastic (slip) loadings that lead to the final force . The remaining (thinner) solid lines are elastic segments, yet they must be stored since they can also affect the future onset and extent of frictional slip (i.e., the loading sequence is nonsimple). Figure 4d shows the final distribution of shear traction across the circular contact area of radius , revealing remnants of the seven equivalent force steps as seven superposed CattaneoMindlin functions. If needed, the traction distribution can be constructed with the Jäger algorithm using principles given in Section 4.1.
Although the equivalent load sequence of the twodimensional ninestep example in Figure 4 is somewhat complex, the full threedimensional sequences in DEM simulations are far more tortuous: instead of a length of 7, the lengths are often greater than 100 (see Figure 1).
We now return to the Jäger algorithm of Figures 2 and 3. A finite movement advances the normal indentation, , and the new normal force is computed with equation (1) (line 13). Jäger showed that the tangential movement can be expressed in terms of three indentations: the current indentation , the previous indentation (stored as ), and an effective (equivalent) indentation :
(4) 
Each indentation corresponds to a normal force (, , and ) and a contact radius through equations (1) and (2). The first case in equation (4) is for an elastic movement, and (4) can be rearranged to give the direction of tangential force, (line 17).
The second case (4) applies when the movement produces annular slip within the contact area. In this case, the norm is capped at the friction coefficient , and direction is aligned with the tangential movement, (lines 62 and 65). When used in a displacementdriven DEM algorithm, the second case serves as the consistency condition for establishing the most recent point of the equivalent history, which is used later to find . (This situation is similar to using a consistency principle to locate the backstress in conventional elastoplasticity with kinematic hardening.) When the previous movement is elastic (with but ), equation (4) is rearranged as
(5) 
where
(6) 
as in lines 24 and 37 (see [20], p. 128). On the other hand, when the previous movement also produces slip (with ), we must distinguish between current movements that produce continued loading — thus expanding the current yield cone with apex at — and movements that produce unloading within the cone (in Figure 4b, the former is represented by steps 4, 5, and 8; the latter by step 9). Loading occurs when the conditions on lines 29–30 are satisfied. When the movement produces unloading, equation (4) is rearranged as
(7) 
as on line 34.
The equivalent force that corresponds to is computed with equation (1) (lines 40–42). Once and are determined, Jäger’s forcespace complement of (4) is used to compute the tangential force increment :
(8) 
as on lines 18 and 67.
The sequence must be monotonically increasing (in an elastoplasticity setting, this requirement obviates the intersection of yield cones). When slip occurs, equations (5) and (7) will give a value , and we must amend the load history: (1) replacing and with and , (2) shifting and to and , and (3) altering and so that they originate from and (lines 49–51 and 57). In Figure 4, these operations would replace the sequence 0–1–2–3–4 with its equivalent sequence 0–1–2–2a–4, and they replace 0–1a–5–6–7–8 with 0–1a–5–6–6a–8. In some cases, , so that the intermediate episode must be eliminated in a repetitive manner (lines 22 and 50–53): in this way, the sequence 0–1–2–2a–4–5 is replaced with 0–1–1a–5.
The algorithm in Figures 2 and 3 is based on the one outlined by Jäger [20]. The lines 28–35 are added for the special case of two successive time steps that produce slip, as in equation (7) (since equation 5 would otherwise involve division by zero). The small tolerance is added in line 28 to avoid computational problems when comparing slopes of force (the author has used an in his trials). Lines 25 and 41 reduce the computational demands of applying equation (1) by retrieving previously computed values of from memory. Three additional provisions are described in Section 3.
The algorithm has been implemented by the author in the OVAL DEM code to simulate the quasistatic loading of large assemblies of particles [28]. At the time of writing, the algorithm has been robust through about fifty billion procedure calls. The algorithm’s performance is also described in Section 4.
3 Implementation Details
3.1 An approximate load history
The algorithm requires sufficient memory to store the three lists , , and . The demand for memory is particularly acute during gradual or quiescent DEM simulations in which contacts remain elastic across many time steps, causing these lists to grow with each passing step. In Figure 5, the list for this contact grows to length 6 during the gradual sequence of elastic steps 0 to 5; whereas the single inelastic step 6 suddenly reduces the list to length 4 (0–1–1a–6) and releases memory that can be used elsewhere.
As an approximation, the lists can be reduced artificially by systematically combining several similar elastic segments of the equivalent loading history. This approximation is achieved by inserting the code of Figure 6 into that of Figure 2, lines 72–81.
When the slopes of two successive elastic segments differ by less than a prescribed small amount (an input parameter with range 0 to ), the approximation combines the two segments into one, thus releasing memory (lines 76–78, in which an average slope is computed and the intermediate point, , is removed). This process is applied repeatedly to combine successive groups of segments (line 74). The approximation algorithm is illustrated in Figure 5 in which the four segments 1–2, 2–3, 3–4, and 4–5 have a similar slope, with , and are combined into the single (stippled) segment 1–5. Although this modification introduces a small error into the history that can affect future inelastic loading steps, the memory savings are considerable. The error in approximating is about 3% for the example in Figure 5, whereas the historylists can be reduced by over 30 percent when the ratio is as small as 0.10.
3.2 Contact frame rotations
In the previous sections, the two particles were assumed only to translate, thus producing the contact movements and . In a more general setting, two contacting particles can both move and rotate in several modes, and the DEM code must extract the single mode of contact translation before the Jäger algorithm can be applied. The four modes of movement between two particles are a contact translation; a rotational contact twisting; a rolling of the particles at their contact; and a rigid rotation of the pair (see [29]). Although twisting and rolling are objective motions that can alter the traction within the circular contact area, they are not considered in the paper — only contact translations are assumed to produce the deformations that give rise to and . Rigid rotations must be considered, however, as they will rotate a contact’s coordinate frame. In an extreme case, two particles can twirl and tumble as a rigid (glued) pair while producing no contact translations at all, even as their contact force is twirled and reoriented within the global (assembly) frame.
The contact translation of particle relative to particle is
(9) 
where , , , and are the particle translations and rotations, and and are the vectors from the centers of particles and to the contact. In this equation, all overbar quantities (for example, ) designate vectors viewed within the global (assembly) coordinate frame. This view will differ from that within the local contact frame of the particle pair (no overbar), and it is within this local frame that the Jäger algorithm operates. The translation vector can be split into normal and tangential parts,
(10)  
(11) 
where is the outward unit normal vector of at the contact.
In a DEM code, we compute the movement in the global frame (equation 11); transform (rotate) the global and vectors into the local contact frame (as and ); use the Jäger algorithm to find in the local contact frame; and then transform the local back into the global so that this force can be used to advance the particles’ positions and orientations. The global frame (overbar vectors with three components) is used throughout the DEM code, except within the Jäger algorithm, where the local contact frame is in effect (i.e., the bare vectors with two possibly nonzero components:  and ). During these coordinate transformations, the list can remain within the contact frame, since it is not needed outside of the Jäger algorithm.
The incremental rotation of the local frame, , is the sum of two types of rigid rotation:
(12) 
a rigid twirling of the two particles about axis ,
(13) 
and a rigid tilting of the normal vector,
(14) 
The latter expression involves the curvature tensors and as described in [30]. For the simple case of two spherical surfaces,
(15) 
where and are the surface radii of curvature.
Quaternions offer a compact and efficient means of shifting between coordinate frames [31]. They are commonly used in DEM codes to store particle orientations and effect particle rotations [14]. The orientation of the contact plane, however, is distinct from the orientations of the two particles, and we will use quaternions to rotate the contact plane and transform between the – and – pairs.
The frame orientation of a contact is expressed as the fourcomponent unit quaternion having a unit norm, . When two particles first touch, their contact’s quaternion is initialized as, say,
(16) 
which is formed from two Euler angles of the contact frame: and . With every time step and before entering the Jäger algorithm, the of each contact is updated to account for the contact’s incremental rotation during the previous step (computed with equation 12) which is embedded in a quaternion :
(17) 
using the quaternion product that is defined in the Appendix. Equation (17) is an approximation that can cause to drift slightly from the unit condition . If desired, can be renormalized, either exactly as or with the Katz approximation [32]. Before entering the Jäger algorithm, the contact’s tangential movement and previous force are rotated into the contact frame,
(18)  
(19) 
using the quaternion conjugate and the double product defined in the Appendix. After leaving the Jäger algorithm, the returned tangential force is rotated back into the global frame:
(20) 
where the force then would be included in Newton’s equations to advance the particles’ positions and orientations.
3.3 Data structures
The equivalent load history of each contact is stored in three lists , , and . These lists are stack structures with lastin/firstout access and with no need for random access into the stacks. Linked lists are compact means of storing and accessing data of this form (see [33], §2.2.3). Rather than creating separate arrays for each contact, three masterlists , , and can store the data of all contacts. Figure 7 gives code for accessing and modifying these lists at various points within Figures 2, 3, and 6.
The input argument of the Jäger procedure (line 1 of Figure 2) is a pointer to the top position of a contact’s history. A fourth list stores descending pointers to the next (lower) position of each contact’s stack, and it also stores ascending pointers to the next available position for placing new data. In this sense, or (i) is the top of the stack of a particular contact; (j(i)) is the next lower position, , in its stack; etc. The j(Zero) value points to the next available position; j(j(Zero)) points to the following available position; etc. Pointer is associated with a single contact, and it can also point to other data (lists) for the contact: in particular, its and data. These lists have a length as long as the number of contacts. The pointer list , however, is associated with the Jäger equivalent load histories of all contacts, and the four lists , , , and are as long as the number of contacts times the average length of the load history per contact: the capacity parameter size_j in Figure 7.
4 Algorithm Performance
4.1 Performance in a simple loadingunloading sequence
In section 17 of [5], Mindlin and Deresiewicz give a closedform solution of movement and tangential force for a simple sequence of proportional loading and unloading. This sequence will serve to illustrate the accuracy of the Jäger algorithm and to compare it with an incrementalstiffness solution — the displacementdriven algorithm of Zhang and VuQuoc [25], employed as an elasticfrictional formulation. In the first stage of loading, two elastic spheres are pressed together with force (see inset of Figure 8a).
The tangential and normal forces are then increased proportionally, with , until nears the frictional limit, . The tangential and normal forces are then reversed in the same proportion until . The solid lines in Figure 8 are the MindlinDeresiewicz solution ([5], equations 80, 83, and 84). This closedform solution is faithfully reproduced by both the Jäger method and by the incrementalstiffness method, provided that the displacement steps (,) are infinitesimal. A more telling test is when the loading sequence is broken into four large steps, with two displacement steps between and , and two more steps between and . Although these steps are much larger than would be expected in DEM simulations, the Jäger method gives the exact MindlinDeresiewicz solution for each large loading step (Figure 8a); indeed, when the Jäger solution is expanded, it coincides with the closedform MindlinDeresiewicz solution. The incrementalstiffness solution is seen to yield an approximation of the MindlinDeresiewicz solution (Figure 8b), although this approximation can be greatly improved by subdividing each large step into smaller subincrements.
Even though shearing tractions within the contact area are rarely computed in DEM simulations, a simple change to the Jäger algorithm yields the shearing traction in addition to the tangential force . We now consider such tractions in order to illustrate the exact correspondence of the Jäger and MindlinDeresiewicz solutions for this simple example of proportional loading and unloading ([5], §14, 15, and 16). Whenever a tangential force (, , , etc.) appears within a line in Figures 2, 3, and 6, the corresponding tangential traction can be computed by simply replacing this force with a traction distribution that is computed by replacing, in turn, that line’s normal forces (, , , etc.) with CatteneoMindlin distributions,
(21) 
In this distribution, is the radial distance from the center of the circular contact area, and “” is a radius computed with equation (2) using the corresponding normal force (, , , etc.).
4.2 Performance with a large assembly of spheres
The algorithm’s computational speed and memory demand are examined with a DEM simulation of biaxial planestrain compression of an assembly of 4096 spheres. The densely packed assembly was compressed in the vertical direction while maintaining constant stress in one horizontal direction and zero strain in the other horizontal direction. As is typical in such tests, the initial deformation is primarily elastic but is followed by inelastic behavior and intense dilation, with the assembly reaching a peak compressive stress at a strain of about [34]. The simulation produced vertical strains of 0–4% (40,000 time steps) during which the assembly contained an average of about 8600 contacts. The two simulations were run on a singlecore singlethread Pentium4 2.66GHz processor with 333MHz memory.
Standard DEM simulations use a central difference explicit time integration scheme in which numerical stability can usually be assured when the time step is less than some limiting value. Tavarez and Plesha [35] give the maximum time step as , where is the damping expressed as a fraction of critical damping, and is the highest vibrational frequency of the entire assembly. Because is difficult to assess, most simulations use a simplified approach based on the contact stiffnesses and particle masses, in which the time step is less than some fraction of (e.g., [36]). With Hertz contacts, stiffness depends on the contact force, and the time step must be adjusted accordingly. The tangential stiffness of the CattaneoMindlinDeresiewicz contact is no greater than the normal (Hertz) stiffness, and the tangential behavior softens with increasing obliquity of the contact force. For these reasons, an upper bound of the time step can be based upon the normal forces and the corresponding Hertz stiffnesses (see [14, 37, 38] for instructions on adaptive time steps with Hertz contacts).
Table 1 compares the computation time for a simulation that employed the Jäger algorithm with one that used a simpler incrementalstiffness modifiedMindlin model, a model that can capture only a single force reversal (e.g., [13]). A runtime profiler (gprof) was used to determine the CPUtime spent computing the contact forces within each simulation. The CPUtime spent within the Jäger algorithm was considerably greater that within the simpler procedure (by a factor of ). A DEM simulation requires other calculations besides force computation, and the total time to run the simulations was only modestly longer when using the Jäger algorithm (by a factor of ).
The maximum memory required by the Jäger algorithm occurred at a strain of %, when the assembly contained about 8700 contacts — during a period of substantial inelastic loading but prior to the peak stress. The corresponding average length of the equivalent load history per contact is also shown in Table 1. The lists had average lengths as great as 260 per contact with the Jäger algorithm. Using four 8byte floatingpoint numbers for , , and and a 4byte integer for , an average history of length 260 requires 9630 bytes of memory per contact — about 83 MBytes for the 8700 contacts. Without this demand, the entire DEM code would require only 10 MBytes of memory, so the Jäger algorithm does impose a substantial memory demand on the simulations. By using the approximation described in Section 3.1, however, the required memory can be reduced considerably. With a parameter of 0.10, which would produce only small errors in the force calculations, the equivalent history is reduced in length from 260 to 171; a larger parameter of 0.20 reduces the average length to 107 (about 40% of the length 260 required for exact force calculations).
5 Conclusion
The Jäger algorithm is an efficient approach to computing the threedimensional CattaneoMindlinDeresiewicz contact force for arbitrary contact translations, regardless of size or sequence. Each movement episode is solved exactly and in whole, without the need to divide the movement into smaller increments. The algorithm is particularly useful for DEM simulations in which particles are persistently in nonsimple contact across many time steps, such as simulations of large inelastic deformation or of cyclic loading. The paper details the necessities for incorporating the algorithm into DEM codes. The computation time is only about 50% greater than when a much simpler contact algorithm is employed, a substantial increase but certainly not an obstacle to its adoption. More problematic, however, is the considerable storage required by the algorithm. Based on the results of a fairly severe performance test, 1 GByte of memory can accommodate DEM simulations of about 40,000 particles, and perhaps 100,000 particles when an approximation is used within the algorithm. Because it authentically captures the contact mechanics of ideal elastic particles, however, the algorithm should be considered for simulations that aim for fidelity to the true behavior of such granular materials.
This material is based upon work supported by the National Science Foundation under Grant No. NEESR936408.
Appendix
References
 [1] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique 1979; 29(1):47–65.
 [2] Timoshenko S, Goodier JN. Theory of Elasticity. 3rd edn., McGrawHill Book Co., Inc.: New York, NY, 1970.
 [3] Cattaneo C. Sul contatto di due corpi elasticiti. Accademia dei Lincei, Rendiconti, Series 6 1938; 27:342–348.
 [4] Mindlin RD. Compliance of elastic bodies in contact. J. Appl. Mech. 1949; 16:259–268.
 [5] Mindlin R, Deresiewicz H. Elastic spheres in contact under varying oblique forces. J. Appl. Mech. 1953; 19(1):327–344.
 [6] Seridi A, Dobry R. An incremental elasticplastic model for the forcedisplacement relation at the contact between two elastic spheres. Research report, Dept. of Civil Engrg., Rensselaer Polytechnic Institute, Troy, NY 1984.
 [7] Dobry R, Ng TT, Petrakis E, Seridi A. General model for contact law between two rough spheres. J. Engrg. Mech. 1991; 117(6):1365–1381.
 [8] VuQuoc L, Zhang X. An accurate and efficient tangential forcedisplacement model of elastic frictional contact in particleflow simulations. Mech. of Materials 1999; 31(4):235–269.
 [9] VuQuoc L, Zhang X. Erratum. Mech. of Materials 1999; 31(11):761–762.
 [10] Walton OR, Braun RL. Viscosity, granulartemperature, and stress calculations for shearing assemblies of inelastic, frictional disks. J. Rheology 1986; 20(5):949–980.
 [11] Zhang X, VuQuoc L. Simulation of chute flow of soybeans using an improved tangential forcedisplacement algorithm. Mech. of Materials 2000; 32(2):115–129.
 [12] Thornton C, Randall CW. Applications of theoretical contact mechanics to solid particle system simulation. Micromechanics of Granular Materials, Satake M, Jenkins J (eds.). Elsevier Science Pub. B.V.: Amsterdam, The Netherlands, 1988; 133–142.
 [13] Lin X, Ng TT. Contact detection algorithms for threedimensional ellipsoids in discrete element modelling. Int. J. Numer. and Analytical Methods in Geomech. 1995; 19(9):653–659.
 [14] VuQuoc L, Zhang X, Walton OR. A 3D discreteelement method for dry granular flows of ellipsoidal particles. Comput. Methods Appl. Mech. Engrg. 2000; 187(3):483–528.
 [15] Jäger J. Elastic contact of equal spheres under oblique forces. Arch. Appl. Mech. 1993; 63(6):402–412.
 [16] Jäger J. Stepwise loading of halfspaces in elliptical contact. J. Appl. Mech. 1996; 63(3):766–773.
 [17] Jäger J. A new principle in contact mechanics. J. Tribology 1998; 120(4):677–684.
 [18] Jäger J. Uniaxial deformation of a random packing of particles. Arch. Appl. Mech. 1999; 69(3):181–203.
 [19] Jäger J. Properties of equal bodies in contact with friction. Int. J. Solids Structures 2003; 40(19):5051–5061.
 [20] Jäger J. New Solutions in Contact Mechanics. WIT Press: Southampton, UK, 2005.
 [21] Thornton C. Coefficient of restitution for collinear collisions of elasticperfectly plastic spheres. J. Appl. Mech. 1997; 64(2):383–386.
 [22] VuQuoc L, Zhang X. An elastoplastic contact forcedisplacement model in the normal direction: displacementdriven version. Proc. R. Soc. Lond. A 1999; 455:4013–4044.
 [23] VuQuoc L, Zhang X, Lesburg L. A normal forcedisplacement model for contacting spheres accounting for plastic deformation: forcedriven formulation. J. Appl. Mech. 2000; 67(2):363–371.
 [24] VuQuoc L, Lesburg L, Zhang X. An accurate tangential forcedisplacement model for granularflow simulations: contacting spheres with plastic deformation, forcedriven formulation. J. Comp. Phys. 2004; 196(1):298–326.
 [25] Zhang X, VuQuoc L. An accurate elastoplastic frictional tangential forcedisplacement model for granularflow simulations: displacementdriven formulation. J. Comp. Phys. 2007; 205(1):730–752.
 [26] Plantard G, Papini M. Mechanical and electrical behaviors of polymer particles. experimental study of the contact area between two particles. experimental validation of a numerical model. Granular Matter 2005; 7(1):1–12.
 [27] Johnson KL. Contact Mechanics. Cambridge Univ. Press, 1985.
 [28] Kuhn MR. OVAL and OVALPLOT: programs for analyzing dense particle assemblies with the Discrete Element Method. http://faculty.up.edu/kuhn/oval/oval.html 2002; .
 [29] Kuhn MR, Bagi K. Contact rolling and deformation in granular media. Int. J. Solids Structures 2004; 41(21):5793–5820.
 [30] Kuhn MR, Bagi K. Alternative definition of particle rolling in a granular assembly. J. Engrg. Mech. 2004; 130(7):826–835.
 [31] Altmann SL. Rotations, Quaternions, and Double Groups. Oxford University Press: New York, 1986.
 [32] Katz A. Computational Rigid Vehicle Dynamics. Krieger Publ. Co.: Malabar, Florida, 1997.
 [33] Knuth DE. The Art of Computer Programming: Fundamental Algorithms, vol. 1. AddisonWesley Pub.: Reading, Mass., 1973.
 [34] Kuhn MR. Micromechanics of fabric and failure in granular materials. Mech. of Materials 2010; 42(9):827–840.
 [35] Tavarez FA, Plesha ME. Discrete element method for modelling solid and particulate materials. Int. J. Numer. Methods Engrg. 2007; 70(4):379–404.
 [36] O’Sullivan C, Bray JD. Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme. Engrg. Computations 2004; 21(2/3/4):278–303.
 [37] Ashmawy AK, Sukumaran B, Hoang VV. Evaluating the influence of particle shape on liquefaction behavior using discrete element modeling. Proc. of the 13th Intl. Offshore and Polar Engrg. Conf., vol. 2, ISOPE, Honolulu, 2003; 542–549.
 [38] Uthus L, Hopkins MA, Horvli I. Discrete element modelling of the resilient behaviour of unbound granular aggregates. Int. J. Pavement Engrg. 2008; 9(6):387–395.
5.00in
\toprule  Algorithm  
modifiedMindlin  Jäger  
\midruleForce calculation time (minutes)  3.4  8.2 
Total simulation time (minutes)  17.6  26.2 
\midruleAvg. length of equiv. history, exact Jäger forces  –  260 
Avg. length of equiv. history,  –  171 
Avg. length of equiv. history,  –  107 
\bottomrule 