Implementation of numerically exact Coulomb friction for discrete element simulations in three dimensions

As a follow-up of an earlier work on the numerically exact Coulomb friction in two-dimensional simulations, we present here the relations and implementation for three-dimensional discrete element particles.


Introduction
The Coulomb friction law (see [1] and references therein for an outline of the historical development) consists of two rules: For finite sliding velocities, dynamic friction (with the absolute value of the normal force |F N | and friction coe cient µ) has a simple functional dependence on the sign of the velocity v and is a dissipative force. The static friction is only determined by the inequality and is a constraint force which opposes any external forces which try to change the relative velocity away from zero. Its treatment as non-smooth force started with [2] and lead to Filippov's use of the convex hull for the solution for static friction [3]. In this tradition, we elaborated the many-body friction for the two-dimensional case [4] and will here treat the formalism for three dimensions. As in [4], we employ only one coe cient of friction µ. Deviations between static and dynamic friction coe cients are usually due to changes in the surface chemistry and proportional to the resting time [5] and will be neglected here, though the formalism permits also the implementation of time dependent coe cients of friction.

Discrete Element Method
The MATLAB-implementation for this research is based on the FORTRAN-code for the discrete element method (DEM) with soft polyhedra (motion according to the corresponding rigid body equations, forces proportional to the geometrical overlap) used in [6,7] which is described in ⇤ e-mail: krengel-dominik-ds@ynu.ac.jp ⇤⇤ e-mail: jchen@jamstec.go.jp ⇤⇤⇤ e-mail: hg@mce.uec.ac.jp Figure 1. Geometry of two particles P 1 , P 2 contacting with the overlap P 1 \ P 2 with the dotted seam @P 1 \ @P 2 in (a), enlarged overlap region, including the triangular contact segments and their normal vectorsñ i (length proportional to the triangle area) in (b) under a di↵erent angle, the averaged contact Plane C and the normal vectorn in (c). more detail in [8]. The overlap-geometry is illustrated in Fig. 1: For two polyhedra P 1 , P 2 , the forcepoint (used also for the computation of the torques) is given by the centroidal point P of the overlap polyhedron P 1 \ P 2 . The connections of the corners of the "seam" @P 1 \ @P 2 (between @P 1 and @P 2 , dashed line in Fig. 1) to the centroid P define triangles: Their normalsñ i , see Fig. 1 (b) are averaged (weighted with the triangle area) ton, which defines the normal of what will later be calld the "contact plane" C. Elastic force and normal damping act parallel tô n through point P. The elastic inter-particle force is computed from the volume V O of the overlap P 1 \P 2 and the Young's modulus Y, which fixes the sound velocity of a space filling particle packing to the continuum material by characteristic length The dissipative normal force is proportional to and the change of the overlap, where V O /⌧ is the change of the overlap volume during a timestep ⌧ in analogy to the harmonic oscillator. The usual prefactor p m red Y with reduced mass of the particles with mass m 1 , m 2 renders the force law scale invariant. The damping F damp is truncated so that it cannot reverse the sign of the (repulsive) normal force, The Coulomb friction acts in the contact plane C which is orthogonal toñ.

Exact implementation of friction
We derive the conditions for static and dynamic friction, as well as the equation for the Lagrange parameter for the one-dimensional case, following [9]. From this, the three-dimensional formalism follows directly for the (onedimensional) projections of three-dimensional quantities into the two-dimensional contact plane C.

One-dimensional case
Let us start with an ordinary di↵erential equation in x for a frictional oscillator with damping, with friction coe cient µ, normal force F N and viscous damping constant 2 . We can now separate the solution according to v > 0 and v < 0 intȯ which means that theẋ-dependency drops out for v = 0. Next, we look at the solution in the phase space of positions and velocities (x(t),ẋ(t)) in Fig. 2. For the velocitydependent force f = mẍ, in eq. (1), we discriminate the cases v > 0 and v < 0, while at the boundary between region I and region II, for v = 0 the inequality eq. (2) is taken into account as the convex hull of f I and f II with a Lagrange-parameter , With the auxiliary parameters the four possible cases at the boundary v = 0 are: 1. a II < 0, a I > 0 signifies dynamic friciton, likewise 2. a I < 0a II > 0, where in both cases there is a jump of the friction force of µF N . For 3. a I < 0, a II < 0: flow is pulled into the constraintẋ = 0 from above and below, which is the case of static friction. 4. a I > 0, a II > 0 is the Painlevé paradox, see discussions in [4]. The Lagrange-parameter 0   1 for static friction is obtained in the convex hull of the dynamic friction for the "indicator function" (constraint function) which is not a specific contact velocity, but the whole field of (v = 0, x) in the dynamic system. Its time derivative gives f (ẋ, ) = 0 for eq. (13) and therefore This yields the value for the Lagrange parameter = a I a I + a II .
Accordingly, for the evaluation of the static friction, only the (relative tangential) accelerations used in a I , a II are necessary, never the (in practice numerically contaminated) contact velocities.

Friction in three dimensions
For three-dimensional particles i = 1, 2 with centroid at r (i) C , velocityṽ (i) C and accelerationã (i) C , the coordinates, velocities and accelerations in the contact point P arẽ with the angular velocity! (i) of particle i. The relative velocities and accelerations at the contact point are accordinglyṽ The one-dimensional case carries over to three dimensions when we consider the vector projection of a vectorw (ṽ rel P orã rel P ) into the contact plane C with normal (but not necessarily normalized) vectorñ defined in Fig. 1 proj plane,n (w) =w w ·ñ ||ñ || 2ñ (23) As the formalism deals only with the tangential motion, we have to use the tangential reduced mass m red,t with the masses m 1 , m 2 of the contacting particles, the inverses of the inertia tensors I 1 , I 2 and the vectors from the centroids to the contact pointr (1) CP ,r (2) CP . That this is the correct form follows from the obvious invariance under rotational coordinate transformations.
If the time integration is not stopped and restarted as in [9], so that the velocity is close to zero, various conditions must be evaluated whether the velocity is small enough or reversed compared to the next or previous timestep. The corresponding conditions, eqs. (35-44) in [4], with sg = 1, 0, 1, are scalars in two dimensions and become two-dimensional vectors for three-dimensional particles, with the tangential reduced mass in eq. (24) for m. The corresponding scalar condition for velocity reversal in two dimensions becomes an inner product of twodimensional vectors for interparticle friction in three dimensions. For finite relative velocities v at the contact, an additional force m red,tang v must be added to bring the particles into relative tangential rest.

Test cases
We test our approach for several example cases to corroborate the validity of di↵erent aspects of our implementation: In sec. 4.2, we investigate the numerical noise, in sec. 4.2 we show why the tangential reduced mass must be employed and in sec. 4.3 we demonstrate the feasibility of our approach for multiple particles and contacts.

Block on a plane
A block on a slope allows to verify that the friction coecient works qualitatively right. with friction coe cient µ = 0.6, we test the motion for slope inclinations below, at and above the critical angle ↵ 1 = arctan(0.8µ), ↵ 2 = arctan(1.0µ) and ↵ 3 = arctan(1.2µ). The height over time is drawn in Fig. 3. As can be seen, the block slides down with a parabolic time dependence of the velocity for ↵ 3 ⇡ 37 , with a constant velocity at for ↵ 2 ⇡ 31 and for ↵ 1 ⇡ 25 , it is practically stuck on the slope. We can see some creep downhill, with very low velocity (⇡2% of the particle length over 2 seconds). This is due to the numerical errors (discretisationand rounding errors), as in the two-dimensional case [4], but also due to some residual rocking motion of the particle, which drives the particle slowly downhill. Reducing the timestep ⌧ reduces the creep marginally, as can be seen in the insert of Fig. 3.

Egg standing up
Typical dynamical problems which involve rotation as well as solid friction are tippe-tops and cooked eggs with rising figure axes, because the center of mass moves faster than the contact point; for viscous friction, such rising is absent [10]. The rising is not limited to "perfectly symmetric" bodies or those with smooth surfaces: Plastic models of lemons and potatoes show the same behaviour as egg-shapes and ellipsoids. Therefore With suitable initial conditions as in Fig. 4 the figure axis does indeed rise, which is a corroboration that the static friction is implemented correctly. Next we compare the behavior for m red,tang and the reduced mass for central forces in eq. (6) for the computation of the friction force. As we can see from Fig. 5, the velocity of the contact point (the "sharp" end of the ellipsoid) on the ground for the tangential reduced mass m red,tang eq. (24) is used for the friction, which corroborates this choice, while with the reduced mass m red eq. (6) for central forces, there is an erratic zig-zagging as the resulting forces are always slightly too large, even for this small friction coe cient. Nevertheless, in spite of the edges and corners of the polyon, the character of solid friction is essentially preserved.

A small heap
From the point of granular materials, the ultimate aim of having an exact friction law is to use it for many particle systems. In particular, a convincing implementation of dry friction allows the formation of heaps on surfaces without roughness but realistic coe cients of friction. Our implementation fulfils this criterion, as can be seen from the heap in Fig. 6 with 91 particles, Y=10 6 N/m 2 , ⇢ = 1 g/cm 3 and a friction coe cient of µ = 0.6.

Summary and conclusions
We have explained the formalism for the numerically exact implementation of Coulomb friction in three dimensions in the framework of analytical mechanics of constraint systems: The concepts for one and two dimensions carry over directly, only the tangential reduced mass is more complex.
In the current, direct implementation of the constraint equations, the accuracy is limited by drift due to the discretisation error and with a spurious coupling of linear forces to the angular degrees of freedom. As in the twodimensional case, the drift is due to the fact that the noise in the angular degree of freedom cannot be fixed by the correct static friction alone. Currently, we are investigating how we can get rid of the residual rotation in a physically justifiable (frame-invariant) way.