Implementation of the Jager contact model for discrete element simulations

TL;DR

This paper introduces a comprehensive implementation of the Jager contact model for 3D discrete element method simulations, enabling efficient force calculations and history storage for complex particle contact sequences.

Contribution

It presents a complete Jager algorithm adapted for DEM, including methods for handling complex motions, optimized data storage, and memory-efficient approximations.

Findings

Algorithm can be integrated into DEM with modest computational overhead.

Supports complex tumbling and twirling contact motions.

Requires increased storage but maintains efficiency.

Abstract

In three-dimensional discrete element method (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 non-linear and path-dependent, requiring an efficient means of computing the forces and storing their histories. By cleverly applying the principles of Cattaneo, Mindlin, and Deresiewicz, Jurgen Jager developed an efficient approach for computing the full three-dimensional force between identical elastic spheres that have undergone difficult movement sequences (J. Jager, New Solutions in Contact Mechanics. WIT Press: Southampton, U.K.). This paper presents a complete Jager algorithm that can be incorporated into DEM codes and also describes three special provisions for DEM simulations: (1) a method for handling particle pairs that undergo complex tumbling and twirling motions in three-dimensions; (2) a compact data structure for storing the loading history of the many contacts in a large assembly; and (3) an approximation of the Jager 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.

Equations

$$F_{\text{n}}=\overline{E}\sqrt{{\zeta }^{3}R}$$

$$a={\left(\frac{F_{\text{n}}R}{\overline{E}}\right)}^{1/3}$$

$${\mathbf{F}}_{\text{t}}={\mathbf{F}}_{\text{t\_old}}+\Delta {\mathbf{F}}_{\text{t}}$$

$$\Delta \mathbf{\xi }=\{\begin{array}{ll}\kappa \mathbf{p}(\zeta -{\zeta }_{\text{old}}),& \Delta \zeta >0\text{ and }|\Delta \mathbf{\xi }/(\kappa \Delta \zeta )|<f\\ \kappa \mathbf{p}(\zeta -{\zeta }^{\ast })-\kappa {\mathbf{p}}_i^{\ast }({\zeta }_{\text{old}}-{\zeta }^{\ast })\text{ and }|\mathbf{p}|=f,& \text{otherwise}\end{array}$$

$${\zeta }^{\ast }=\zeta -{\left[\sqrt{{\left(\frac{\Delta {\mathbf{e}}^{\prime }\cdot {\mathbf{p}}_i^{\ast }}{|\Delta {\mathbf{e}}^{\prime }{|}^2}\right)}^2+\frac{f^2-|{\mathbf{p}}_i^{\ast }{|}^2}{|\Delta {\mathbf{e}}^{\prime }{|}^2}}-\frac{\Delta {\mathbf{e}}^{\prime }\cdot {\mathbf{p}}_i^{\ast }}{|\Delta {\mathbf{e}}^{\prime }{|}^2}\right]}^{-1}$$

$$\Delta {\mathbf{e}}^{\prime }=\frac{1}{\kappa }\Delta \mathbf{\xi }-{\mathbf{p}}_i^{\ast }(\zeta -{\zeta }_{\text{old}})$$

$${\zeta }^{\ast }=\zeta +\frac{1}{2}\frac{|\Delta {\mathbf{e}}^{\prime }{|}^2}{\Delta {\mathbf{e}}^{\prime }\cdot {\mathbf{p}}_i^{\ast }}$$

$$\Delta {\mathbf{F}}_{\text{t}}=\{\begin{array}{l}\mathbf{p}(F_{\text{n}}-F_{\text{n\_old}})\\ \mathbf{p}(F_{\text{n}}-F_{\text{n}}^{\ast })-{\mathbf{p}}_i^{\ast }(F_{\text{n\_old}}-F_{\text{n}}^{\ast })\end{array}$$

$$\Delta {\overline{\mathbf{u}}}^{\text{def}}=\Delta {\overline{\mathbf{u}}}^q-\Delta {\overline{\mathbf{u}}}^p+\left(\Delta {\overline{\mathbf{\theta }}}^q\times {\overline{\mathbf{r}}}^q-\Delta {\overline{\mathbf{\theta }}}^p\times {\overline{\mathbf{r}}}^p\right)$$

$$\Delta \zeta$$

$$\Delta \overline{\mathbf{\xi }}$$

$$\Delta \overline{\mathbf{\omega }}=\Delta {\overline{\mathbf{\omega }}}^{\text{twirl}}+{\overline{\mathbf{n}}}^p\times \Delta {\overline{\mathbf{n}}}^p$$

$$\Delta {\overline{\mathbf{\omega }}}^{\text{twirl}}=\frac{1}{2}\left[\left(\Delta {\overline{\mathbf{\theta }}}^p+\Delta {\overline{\mathbf{\theta }}}^q\right)\cdot {\overline{\mathbf{n}}}^p\right]{\overline{\mathbf{n}}}^p$$

$$\Delta {\overline{\mathbf{n}}}^p=\Delta {\overline{\mathbf{\theta }}}^p\times {\overline{\mathbf{n}}}^p+{\overline{\mathbf{K}}}^p\cdot {\left({\overline{\mathbf{K}}}^p+{\overline{\mathbf{K}}}^q\right)}^{-1}\cdot \left[(\Delta {\overline{\mathbf{\theta }}}^q-\Delta {\overline{\mathbf{\theta }}}^p)\times {\overline{\mathbf{n}}}^p-2{\overline{\mathbf{K}}}^q\cdot \Delta \overline{\mathbf{\xi }}\right]$$

$$\Delta {\overline{\mathbf{n}}}^p=\Delta {\overline{\mathbf{\theta }}}^p\times {\overline{\mathbf{n}}}^p+\frac{R^q}{R^p+R^q}\left[(\Delta {\overline{\mathbf{\theta }}}^q-\Delta {\overline{\mathbf{\theta }}}^p)\times {\overline{\mathbf{n}}}^p+\frac{2}{R^q}\Delta \overline{\mathbf{\xi }}\right]$$

$$\stackrel{\circ }{\mathbf{q}}=⟨\cos (\beta /2),\sin (\beta /2)[\sin (\eta ),-\cos (\eta ),0]⟩$$

$$\stackrel{\circ }{\mathbf{q}}\approx {\stackrel{\circ }{\mathbf{q}}}_{\text{old}}+\frac{1}{2}\Delta \stackrel{\circ }{\omega }\circ {\stackrel{\circ }{\mathbf{q}}}_{\text{old}}$$

$$\Delta \mathbf{\xi }$$

$${\mathbf{F}}_{\text{t\_old}}$$

$${\overline{\mathbf{F}}}_{\text{t\_old}}=\stackrel{\circ }{\mathbf{q}}{}^{\ast }\circ {\mathbf{F}}_{\text{t\_old}}\circ \stackrel{\circ }{\mathbf{q}}$$

$${\sigma }_{\text{C-M}}(\rho )=\{\begin{array}{ll}\frac{3\overline{E}}{2\pi R}{\left(a^2-{\rho }^2\right)}^{1/2}& \rho <a\\ 0& \rho \ge a\end{array}$$

$$\frac{1}{2}\Delta \stackrel{\circ }{\omega }\circ {\stackrel{\circ }{\mathbf{q}}}_{\text{old}}=\frac{1}{2}\left[\begin{array}{cccc}0& -\Delta {\overline{\omega }}_1& -\Delta {\overline{\omega }}_2& -\Delta {\overline{\omega }}_3\\ \Delta {\overline{\omega }}_1& 0& -\Delta {\overline{\omega }}_3& \Delta {\overline{\omega }}_2\\ \Delta {\overline{\omega }}_2& \Delta {\overline{\omega }}_3& 0& -\Delta {\overline{\omega }}_1\\ \Delta {\overline{\omega }}_3& -\Delta {\overline{\omega }}_2& \Delta {\overline{\omega }}_1& 0\end{array}\right]\left[\begin{array}{c}{q}_{\text{old},0}\\ q_{\text{old},1}\\ q_{\text{old},2}\\ q_{\text{old},3}\end{array}\right]$$

$$\stackrel{\circ }{\mathbf{q}}\circ \mathbf{a}\circ \stackrel{\circ }{\mathbf{q}}{}^{\ast }$$

$$\stackrel{\circ }{\mathbf{q}}{}^{\ast }\circ \mathbf{a}\circ \stackrel{\circ }{\mathbf{q}}$$

$$\left[\mathbf{Q}\right]=\left[\begin{array}{ccc}\frac{1}{2}-q_2^2-q_3^2& q_0{q}_3+q_1{q}_2& q_1{q}_3-q_0{q}_2\\ q_1{q}_2-q_0{q}_3& \frac{1}{2}-q_1^2-q_3^2& q_0{q}_1+q_2{q}_3\\ q_1{q}_3+q_0{q}_2& q_2{q}_3-q_0{q}_1& \frac{1}{2}-q_1^2-q_2^2\end{array}\right]$$

Full text

\NME

16000010 \runningheadsM. R. Kuhn Jäger Contact for DEM Simulations

\corraddrUniversity of Portland, 5000 N. Willamette Blvd., Portland, OR, 97203 USA, Tel. (1)-503-943-7361, Fax (1)-503-943-7316

\cgsnNational Science FoundationNEESR-936408 \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 displacement-force 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 history-dependent, 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 (quasi-static) 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 $\zeta$ and the two components of tangential movement, $\xi_{1}$ and $\xi_{2}$, 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 elasto-plasticity 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 Mindlin-Deresiewicz limitations by developing a general algorithm for the three-dimensional contact displacement-force relation that was analogous to incremental elasto-plasticity 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, Vu-Quoc and Zhang [8, 9] developed another tangential force-displacement algorithm that uses an incremental-stiffness form of the Mindlin-Deresiewicz 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 partially-latching spring system for the normal forces [10, 8]. This system models plastic deformation in the normal direction. Zhang and Vu-Quoc [11] demonstrated their contact by modeling rapid, collisional flows of soybeans. Most DEM codes now use a Hertzian normal force combined with some incremental-stiffness form of the Mindlin-Deresiewicz model (for example, [12, 13, 14]).

In a remarkable series of papers and book, Jäger developed an elegant approach to the Cattaneo-Mindlin-Deresiewicz 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 Cattaneo-Mindlin 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 displacement-driven 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 normal-tangential movement can be of arbitrary size without the need for incremental stiffnesses and without dividing a large movement into smaller, incremental sub-movements. In this sense, the Jäger algorithm reproduces the closed-form 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. (Mindlin-Deresiewicz 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 is inherently three-dimensional, delivering the final three-dimensional force for an arbitrary sequence of normal and 2-vector tangential movements [19, 20].

• •

The algorithm includes a rigorous means of identifying load reversals (turning-points) in a full three-dimensional 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 non-simple [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 non-simple loading sequences in which a general elastic movement is followed by a non-elastic 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 Cattaneo-Mindlin tractions must be superposed (see [17, 18] for succinct derivations).

Because both elastic and non-elastic movements can affect future slip, the Jäger algorithm thoroughly chronicles any changes in the directions of either elastic or non-elastic movements, including those at turning-points. The algorithm then accounts for the effect that past elastic and non-elastic movements will have upon the current and future forces. In short, the Jäger algorithm exactly computes the contact force for non-simple 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 three-dimensional 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, visco-elasticity, or fracture and breaking of the spheres. Investigators have recently modeled the contact problem as an elasto-plastic process [10, 21]. Vu-Quoc and Zhang modeled the normal force-displacement relationship of elasto-plastic spheres [22, 23] and, more recently, developed an algorithm for the tangential and normal interactions of elasto-plastic 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 pseudo-code 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 Mindlin-Deresiewicz results and compares the Jäger solution with an incremental-stiffness solution.

2 JÄGER ALGORITHM

We consider two identical isotropic-elastic spheres that undergo an arbitrary sequence of finite, perhaps large, translations. Contact forces $F_{\text{n}}$ and $\mathbf{F}_{\text{t}}$ are the compressive normal force and the 2-vector tangential force. Displacements $2\zeta$ and $2\boldsymbol{\xi}$ are the cumulative normal approach (overlap) and the 2-vector 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 2-vector tangential forces, $F_{\text{n}}$ and $\mathbf{F}_{\text{t}}$ (line 2). The input includes the contact half-overlap (indentation) $\zeta$; the 2-vector tangential half-shift movement $\Delta\boldsymbol{\xi}$; the friction coefficient $f$; the modulus $\overline{E}=8G/(3(1-\nu))$, which depends on the particles’ shear modulus $G$ and Poisson ratio $\nu$; the ratio $\kappa$ of elastic normal and tangential stiffnesses, $\kappa=(2-\nu)/(2(1-\nu))$; the particles’ radii of curvature, $R_{A}$ and $R_{B}$; the previous tangential force $\mathbf{F}_{\text{t\_old}}$; a small tolerance parameter $\epsilon$; a pointer $i$ to the top of an equivalent load history stack, described in Section 3.3; and a parameter $|\Delta\mathbf{p}_{\text{max}}|^{2}$ 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., $\zeta^{\ast}_{\ldots}$, $F^{\ast}_{\text{n,}\ldots}$, $\mathbf{p}^{\ast}_{\ldots}$, and $j_{\ldots}$) are described later.

The Hertz normal force only depends on the indentation $\zeta$ [2]:

[TABLE]

where the modulus $\overline{E}=8G/(3(1-\nu))$ (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

[TABLE]

and they gave the radii of stick areas the symbols $b$, $c^{\ast}$, etc. Although these radii could be used within the Jäger algorithm, it is more convenient to use $\zeta$ and $F_{\text{n}}$ 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 $\mathbf{F}_{\text{t}}$ will depend upon the DEM history of the finite movements (translations) between time steps, the $\Delta\zeta$ and $\Delta\boldsymbol{\xi}$, that lead to the contact’s current position $(\zeta,\boldsymbol{\xi})$. The current force $\mathbf{F}_{\text{t}}$ is produced by a force increment $\Delta\mathbf{F}_{\text{t}}$ relative to the previous force (i.e., from the previous time step):

[TABLE]

(line 71). The tangential increment $\Delta\mathbf{F}_{\text{t}}$ is entirely elastic, producing no slip, when the current movement $(\Delta\zeta,\Delta\boldsymbol{\xi})$ satisfies the following two conditions: $\Delta\zeta>0$ and $|\Delta\boldsymbol{\xi}/(\kappa\Delta\zeta)|<f$ (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 ($\ast$-star) load history is a compact recording of an equivalent sequence of loading steps that would lead to the previous contact force $\mathbf{F}_{\text{t,old}}$ and the cumulative displacement $\boldsymbol{\xi}_{\text{old}}$. The equivalent load history is stored as lists of the normal indentations and the corresponding normal forces — the lists $\zeta^{\ast}_{0}$, $\zeta^{\ast}_{1},\;\ldots$ $\zeta^{\ast}_{i}$ and $F_{\text{n,}0}^{\ast}$, $F_{\text{n,}1}^{\ast},\;\ldots$ $F_{\text{n,}i}^{\ast}$ — along with a list of 2-vector directions of tangential force, $\mathbf{p}^{\ast}_{0}$, $\mathbf{p}^{\ast}_{1},\;\ldots$ $\mathbf{p}^{\ast}_{i}$. Each $\zeta^{\ast}$-$F_{\text{n}}^{\ast}$ pair designates the apex of a yield cone in displacement-space and force-space. In another sense, the $\zeta^{\ast}$ and $F_{\text{n}}^{\ast}$ are proxies for the radii of past slip areas, as in equation (2) and in reference [5]. This “$\ast$” data gives information about the points of load reversals (turning-points) 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 [math] to $i$. The zeroth point is permanently initialized to the unloaded condition: $\zeta^{\ast}_{0}=F_{\text{n,}0}^{\ast}=0$ and $\mathbf{p}^{\ast}_{0}=\mathbf{0}$. 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 $\zeta$, $F_{\text{n}}$, and $\mathbf{p}$ are appended to the equivalent history (lines 70–71). Because these lists are more compactly stored as linked lists, the indices 0,1,$\ldots i$ 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 two-dimensional example, the tangential movements are colinear, so only single components of the two-vectors $\boldsymbol{\xi}$ and $\mathbf{F}_{\text{t}}$ are relevant. The first part, Figure 4a, shows an example sequence of nine movements between two spheres having friction coefficient $f=0.5$ and $\kappa f=0.75$, where $\kappa=(2-\nu)/(2(1-\nu))$.

The particles first touch at stage “0” and then undergo a sequence of nine proportional (straight) normal-tangential movements: $(\Delta\zeta_{1},\Delta\xi_{1})$, $(\Delta\zeta_{2},\Delta\xi_{2})$,$\ldots$, $(\Delta\zeta_{9},\Delta\xi_{9})$, These successive movements produce the corresponding sequence of forces in the second part, Figure 4b, labeled as the points 1, 2, $\ldots$, 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–$\ldots$–7a–9) would produce the same final displacement ($\zeta_{9}$, $\boldsymbol{\xi}_{9}$), the same contact force ($F_{\text{n,9}}$, $\mathbf{F}_{\text{t,9}}$), and the same tractions as the original nine steps (0–1–2–$\ldots$–8–9). The dashed lines in Figure4b trace the intermediate equivalent histories that would lead to the cumulative intermediate forces $\mathbf{F}_{\text{t,1}}$, $\mathbf{F}_{\text{t,2}}$, etc. These intermediate (dashed) segments were eventually eliminated within the equivalent history that leads to the final force $\mathbf{F}_{\text{t,9}}$. For example, the sequence 0–1–1a–5 will produce the same cumulative displacement ($\zeta_{5}$, $\boldsymbol{\xi}_{5}$) 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 $\Delta\xi/\Delta\zeta$ 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 $\mathbf{p}$ no steeper than $\pm f$. The four slip movements result in equivalent path segments that have a slope magnitude $|\Delta\mathbf{F}_{\text{t}}/\Delta F_{\text{n}}|=|\mathbf{p}|=f$: 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 non-elastic (slip) loadings that lead to the final force $\mathbf{F}_{\text{t,9}}$. 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 non-simple). Figure 4d shows the final distribution of shear traction across the circular contact area of radius $a$, revealing remnants of the seven equivalent force steps as seven superposed Cattaneo-Mindlin 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 two-dimensional nine-step example in Figure 4 is somewhat complex, the full three-dimensional 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 $(\Delta\zeta,\Delta\boldsymbol{\xi})$ advances the normal indentation, $\zeta=\zeta_{\text{old}}+\Delta\zeta$, and the new normal force $F_{\text{n}}$ is computed with equation (1) (line 13). Jäger showed that the tangential movement $\Delta\boldsymbol{\xi}$ can be expressed in terms of three indentations: the current indentation $\zeta$, the previous indentation $\zeta_{\text{old}}$ (stored as $\zeta^{\ast}_{i}$), and an effective (equivalent) indentation $\zeta^{\ast}$:

[TABLE]

Each indentation corresponds to a normal force ($F_{\text{n}}$, $F_{\text{n,old}}$, and $F_{\text{n}}^{\ast}$) and a contact radius through equations (1) and (2). The first case in equation (4) is for an elastic movement, and (41) can be rearranged to give the direction of tangential force, $\mathbf{p}=\Delta\mathbf{F}_{\text{t}}/\Delta F_{\text{n}}=\Delta\boldsymbol{\xi}/(\kappa\Delta\zeta)$ (line 17).

The second case (42) applies when the movement produces annular slip within the contact area. In this case, the norm $|\mathbf{p}|$ is capped at the friction coefficient $f$, and direction $\mathbf{p}$ is aligned with the tangential movement, $\mathbf{p}=f\,\Delta\boldsymbol{\xi}/|\Delta\boldsymbol{\xi}|$ (lines 62 and 65). When used in a displacement-driven DEM algorithm, the second case serves as the consistency condition for establishing the most recent point $\zeta^{\ast}$ of the equivalent history, which is used later to find $\Delta\mathbf{F}_{\text{t}}$. (This situation is similar to using a consistency principle to locate the back-stress in conventional elasto-plasticity with kinematic hardening.) When the previous movement is elastic (with $|\mathbf{p}^{\ast}_{i}|<f$ but $|\mathbf{p}|=f$), equation (42) is rearranged as

[TABLE]

where

[TABLE]

as in lines 24 and 37 (see [20], p. 128). On the other hand, when the previous movement also produces slip (with $|\mathbf{p}^{\ast}_{i}|=f$), we must distinguish between current movements $\Delta\boldsymbol{\xi}$ that produce continued loading — thus expanding the current yield cone with apex at $\zeta^{\ast}_{i-1}$ — 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 (42) is rearranged as

[TABLE]

as on line 34.

The equivalent force $F_{\text{n}}^{\ast}$ that corresponds to $\zeta^{\ast}$ is computed with equation (1) (lines 40–42). Once $\zeta^{\ast}$ and $F_{\text{n}}^{\ast}$ are determined, Jäger’s force-space complement of (4) is used to compute the tangential force increment $\Delta\mathbf{F}_{\text{t}}$:

[TABLE]

as on lines 18 and 67.

The sequence $\zeta^{\ast}_{0},\;\zeta^{\ast}_{1},\;\ldots\;\zeta^{\ast}_{i}$ must be monotonically increasing (in an elasto-plasticity setting, this requirement obviates the intersection of yield cones). When slip occurs, equations (5) and (7) will give a value $\zeta^{\ast}\leq\zeta^{\ast}_{i}(=\zeta_{\text{old}})$, and we must amend the load history: (1) replacing $\zeta^{\ast}_{i}$ and $F_{\text{n,}i}^{\ast}$ with $\zeta^{\ast}$ and $F_{\text{n}}^{\ast}$, (2) shifting $\zeta_{\text{old}}$ and $F_{\text{n\_old}}$ to $\zeta^{\ast}$ and $F_{\text{n}}^{\ast}$, and (3) altering $\mathbf{F}_{\text{t\_old}}$ and $\Delta\boldsymbol{\xi}$ so that they originate from $F_{\text{n\_old}}$ and $\zeta_{\text{old}}$ (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, $\zeta^{\ast}\leq\zeta^{\ast}_{i-1}$, so that the intermediate episode $\zeta^{\ast}_{i}$ 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 $\epsilon$ is added in line 28 to avoid computational problems when comparing slopes of force (the author has used an $\epsilon=1\times 10^{-10}$ in his trials). Lines 25 and 41 reduce the computational demands of applying equation (1) by retrieving previously computed values of $F_{\text{n}}^{\ast}$ 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 quasi-static 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 $\zeta^{\ast}_{0,1,\ldots,i}$, $F^{\ast}_{\text{n},0,1,\ldots,i}$, and $\mathbf{p}^{\ast}_{0,1,\ldots,i}$. 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 $\mathbf{p}^{\ast}$ slopes of two successive elastic segments differ by less than a prescribed small amount $|\Delta\mathbf{p}_{\text{max}}|$ (an input parameter with range 0 to $2f$), the approximation combines the two segments into one, thus releasing memory (lines 76–78, in which an average slope $\mathbf{p}_{i}^{\ast}$ is computed and the intermediate point, $i-1$, 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 $|\mathbf{p}^{\ast}_{i}-\mathbf{p}^{\ast}_{i-1}|\leq|\Delta\mathbf{p}_{\text{max}}|$, 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 $\mathbf{F}_{\text{t,}6}$ is about 3% for the example in Figure 5, whereas the history-lists can be reduced by over 30 percent when the ratio $|\Delta\mathbf{p}_{\text{max}}|/f$ 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 $\Delta\zeta$ and $\Delta\boldsymbol{\xi}$. 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 $F_{\text{n}}$ and $\mathbf{F}_{\text{t}}$. 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 $\Delta\overline{\mathbf{u}}^{\text{def}}$ of particle $q$ relative to particle $p$ is

[TABLE]

where $\Delta\overline{\mathbf{u}}^{p}$, $\Delta\overline{\mathbf{u}}^{q}$, $\Delta\overline{\boldsymbol{\theta}}^{p}$, and $\Delta\overline{\boldsymbol{\theta}}^{p}$ are the particle translations and rotations, and $\overline{\mathbf{r}}^{p}$ and $\overline{\mathbf{r}}^{q}$ are the vectors from the centers of particles $p$ and $q$ to the contact. In this equation, all overbar quantities (for example, $\Delta\overline{\mathbf{u}}^{p}$) 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,

[TABLE]

where $\overline{\mathbf{n}}^{p}$ is the outward unit normal vector of $p$ at the contact.

In a DEM code, we compute the $\Delta\overline{\boldsymbol{\xi}}$ movement in the global frame (equation 11); transform (rotate) the global $\Delta\overline{\boldsymbol{\xi}}$ and $\overline{\mathbf{F}}_{\text{t\_old}}$ vectors into the local contact frame (as $\Delta\boldsymbol{\xi}$ and $\mathbf{F}_{\text{t\_old}}$); use the Jäger algorithm to find $\mathbf{F}_{\text{t}}$ in the local contact frame; and then transform the local $\mathbf{F}_{\text{t}}$ back into the global $\overline{\mathbf{F}}_{\text{t}}$ 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 non-zero components: $\Delta\xi_{1}$-$\Delta\xi_{2}$ and $F_{\text{t},1}$-$F_{\text{t},2}$). During these coordinate transformations, the $\mathbf{p}^{\ast}$ 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, $\Delta\overline{\boldsymbol{\omega}}$, is the sum of two types of rigid rotation:

[TABLE]

a rigid twirling of the two particles about axis $\overline{\mathbf{n}}^{p}$,

[TABLE]

and a rigid tilting $\Delta\overline{\mathbf{n}}^{p}$ of the normal vector,

[TABLE]

The latter expression involves the curvature tensors $\overline{\mathbf{K}}^{p}$ and $\overline{\mathbf{K}}^{q}$ as described in [30]. For the simple case of two spherical surfaces,

[TABLE]

where $R^{p}$ and $R^{q}$ 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 $\Delta\overline{\boldsymbol{\xi}}$–$\overline{\mathbf{F}}_{\text{t}}$ and $\Delta\boldsymbol{\xi}$–$\mathbf{F}_{\text{t}}$ pairs.

The frame orientation of a contact is expressed as the four-component unit quaternion $\overset{\circ}{\mathbf{q}}=\langle q_{0},\mathbf{q}\rangle=\langle q_{0},q_{1},q_{2},q_{3}\rangle$ having a unit norm, $|\overset{\circ}{\mathbf{q}}|^{2}=\sum{q_{i}^{2}}=1$. When two particles first touch, their contact’s quaternion is initialized as, say,

[TABLE]

which is formed from two Euler angles of the contact frame: $\beta=\cos^{-1}(\overline{n}^{p}_{3})$ and $\eta=\tan^{-1}(\overline{n}^{p}_{2}/\overline{n}^{p}_{1})$. With every time step and before entering the Jäger algorithm, the $\overset{\circ}{\mathbf{q}}$ of each contact is updated to account for the contact’s incremental rotation $\Delta\overline{\boldsymbol{\omega}}$ during the previous step (computed with equation 12) which is embedded in a quaternion $\Delta\overset{\circ}{\boldsymbol{\omega}}=\langle 0,\Delta\overline{\boldsymbol{\omega}}\rangle$:

[TABLE]

using the quaternion product $\overset{\circ}{\mathbf{a}}\circ\overset{\circ}{\mathbf{b}}$ that is defined in the Appendix. Equation (17) is an approximation that can cause $\overset{\circ}{\mathbf{q}}$ to drift slightly from the unit condition $|\overset{\circ}{\mathbf{q}}|=1$. If desired, $\overset{\circ}{\mathbf{q}}$ can be renormalized, either exactly as $\overset{\circ}{\mathbf{q}}\leftarrow\overset{\circ}{\mathbf{q}}/|\overset{\circ}{\mathbf{q}}|$ or with the Katz approximation [32]. Before entering the Jäger algorithm, the contact’s tangential movement $\Delta\overline{\boldsymbol{\xi}}$ and previous force $\overline{\mathbf{F}}_{\text{t\_old}}$ are rotated into the contact frame,

[TABLE]

using the quaternion conjugate $\overset{\circ}{\mathbf{q}}\!\,^{\ast}$ 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:

[TABLE]

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 $\zeta^{\ast}_{0,1,\ldots,i}$, $F^{\ast}_{\text{n},0,1,\ldots,i}$, and $\mathbf{p}^{\ast}_{0,1,\ldots,i}$. These lists are stack structures with last-in/first-out 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 master-lists $\zeta^{\ast}$, $F_{\text{n}}^{\ast}$, and $\mathbf{p}^{\ast}$ 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 $i$ 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 $j$ 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, $\zeta^{\ast}_{i}$ or $\zeta^{\ast}$(i) is the top of the $\zeta^{\ast}$-stack of a particular contact; $\zeta^{\ast}$(j(i)) is the next lower position, $\zeta^{\ast}_{i-1}$, 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 $i$ is associated with a single contact, and it can also point to other data (lists) for the contact: in particular, its $\overset{\circ}{\mathbf{q}}$ and $\mathbf{F}_{\text{t\_old}}$ data. These lists have a length as long as the number of contacts. The pointer list $j$, however, is associated with the Jäger equivalent load histories of all contacts, and the four lists $j$, $\zeta^{\ast}$, $F_{\text{n}}^{\ast}$, and $\mathbf{p}^{\ast}$ 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 loading-unloading

sequence

In section 17 of [5], Mindlin and Deresiewicz give a closed-form 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 incremental-stiffness solution — the displacement-driven algorithm of Zhang and Vu-Quoc [25], employed as an elastic-frictional formulation. In the first stage of loading, two elastic spheres are pressed together with force $F_{\text{n},0}$ (see inset of Figure 8a).

The tangential and normal forces are then increased proportionally, with $dF_{\text{t}}/dF_{\text{n}}=4>f$, until $F_{\text{t}}$ nears the frictional limit, $F_{\text{t,B}}=0.9fF_{\text{n,B}}$. The tangential and normal forces are then reversed in the same proportion until $F_{\text{t,D}}=-0.9fF_{\text{n,D}}$. The solid lines in Figure 8 are the Mindlin-Deresiewicz solution ([5], equations 80, 83, and 84). This closed-form solution is faithfully reproduced by both the Jäger method and by the incremental-stiffness method, provided that the displacement steps ($\Delta\zeta$,$\Delta\xi$) are infinitesimal. A more telling test is when the loading sequence is broken into four large steps, with two displacement steps between $F_{\text{t,0}}$ and $F_{\text{t,B}}$, and two more steps between $F_{\text{t,B}}$ and $F_{\text{t,D}}$. Although these steps are much larger than would be expected in DEM simulations, the Jäger method gives the exact Mindlin-Deresiewicz solution for each large loading step (Figure 8a); indeed, when the Jäger solution is expanded, it coincides with the closed-form Mindlin-Deresiewicz solution. The incremental-stiffness solution is seen to yield an approximation of the Mindlin-Deresiewicz solution (Figure 8b), although this approximation can be greatly improved by sub-dividing each large step into smaller sub-increments.

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 $\boldsymbol{\tau}$ in addition to the tangential force $\mathbf{F}_{\text{t}}$. We now consider such tractions in order to illustrate the exact correspondence of the Jäger and Mindlin-Deresiewicz solutions for this simple example of proportional loading and unloading ([5], §14, 15, and 16). Whenever a tangential force ($\mathbf{F}_{\text{t}}$, $\mathbf{F}_{\text{t\_old}}$, $\Delta\mathbf{F}_{\text{t}}$, 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 $\boldsymbol{\tau}(\rho)$ that is computed by replacing, in turn, that line’s normal forces ($\mathbf{F}_{\text{n}}$, $\mathbf{F}_{\text{n\_old}}$, $\mathbf{F}_{\text{n}}^{\ast}$, etc.) with Catteneo-Mindlin distributions,

[TABLE]

In this distribution, $\rho$ is the radial distance from the center of the circular contact area, and “$a$” is a radius computed with equation (2) using the corresponding normal force ($\mathbf{F}_{\text{n}}$, $\mathbf{F}_{\text{n\_old}}$, $\mathbf{F}_{\text{n}}^{\ast}$, etc.).

Figure 9 shows the Jäger and Mindlin-Deresiewicz solutions at three stages in the loading and unloading sequence of Figure 8a. The two solutions coincide exactly.

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 plane-strain 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 $2.5\%$ [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 single-core single-thread Pentium-4 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 $2(\sqrt{1-\xi}-\xi)/\omega_{\text{max}}$, where $\xi$ is the damping expressed as a fraction of critical damping, and $\omega_{\text{max}}$ is the highest vibrational frequency of the entire assembly. Because $\omega_{\text{max}}$ 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 $\sqrt{m/k}$ (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 Cattaneo-Mindlin-Deresiewicz 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 incremental-stiffness modified-Mindlin model, a model that can capture only a single force reversal (e.g., [13]). A run-time profiler (gprof) was used to determine the CPU-time spent computing the contact forces within each simulation. The CPU-time spent within the Jäger algorithm was considerably greater that within the simpler procedure (by a factor of $2.4$). 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 $1.5$).

The maximum memory required by the Jäger algorithm occurred at a strain of $0.5$%, 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 8-byte floating-point numbers for $\zeta^{\ast}$, $F_{\text{n}}^{\ast}$, and $\mathbf{p}^{\ast}$ and a 4-byte integer for $j$, 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 $|\Delta\mathbf{p}_{\text{max}}|$ of 0.10$f$, 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$f$ 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 three-dimensional Cattaneo-Mindlin-Deresiewicz 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 non-simple 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.

\ack

This material is based upon work supported by the National Science Foundation under Grant No. NEESR-936408.

APPENDIX

The quaternion product in (17) is the matrix product

[TABLE]

The vector rotations on the right of Eqs. (18), (19), and (20) are the matrix products

[TABLE]

where $\overset{\circ}{\mathbf{q}}$ is a unit quaternion, $\overset{\circ}{\mathbf{q}}\!\,^{\ast}$ is its conjugate, $\mathbf{a}$ is a 3-vector, and

[TABLE]