Comparing numerically exact and modelled static friction

. Currently there exists no mechanically consistent "numerically exact" implementation of static and dynamic Coulomb friction for general soft particle simulations with arbitrary contact situations in two or three dimension, but only along one dimension. We outline a di ﬀ erential-algebraic equation approach for a "numerically exact" computation of friction in two dimensions and compare its application to the Cundall-Strack model in some test cases.


Introduction
Coulomb sliding friction is a dissipative force, given by the equation where μ is the coefficient of friction, F normal the normal force and v the velocity. Its implementation is computationally as simple as for viscous friction. Static Coulomb friction, on the other hand, is a constraint of motion −μF normal ≤ F Fr, dyn ≤ μF normal , which fixes two bodies tangentially to each other. Its implementation is algorithmically much more demanding: For zero tangential velocity, the tangential inter-particle force must compensate the external forces and fullfill inequality (2). Currently, there exists no mechanically consistent theory of Coulomb friction for two-or three dimensional particle aggregates for arbitrary contact situations. Existing exact algorithms are only suitable for one dimension, see Sec. 2. In this research, we outline the formalism for the numerically exact implementation of Coulomb friction in generalisation of the one-dimensional formalism and show the implementations for some test cases. Rolling friction is neglected, it is generally one or two orders of magnitude lower than sliding friction. The same friction coefficient is used for dynamic and static friction, measurable differences are usually due to time-dependant surface changes (e.g. accumulation of humidity in interstices) which cannot be absorbed into mere constants anyway, see Ref. [9].
relative sliding distance of the contacting particles. The drawback is that this models a degree of freedom, which may release energy into the system when there is no preferred direction of motion (e.g. in vibrated systems), not a constraint force. Stewart [11] proposed a mathematical formulation where the contact situation is absorbed into a linear system which is solved to obtain the friction values. Nevertheless, only examples for contacts in one dimension were treated and it was left open how the inequality eq. (2) can be satisfied mathematically. While formalisms in convex analysis allow the description of tangential interparticle forces which fullfill the formal conditions, there is no guarantee that the variation during time evolution is physical or even smooth. Moreau et al. [8,10] have treated the problem in their "contact mechanics" in a generalization of Newtonian kinematics which we would like to avoid (also for the sake of applicability of numerical integrators). Further, in contact mechanics the normal and tangential forces are evaluated simultaneously, while we would like to obtain friction forces as reactions to the existing normal forces. Finally, the approach used rigid particles, which makes the investigation of some problems unfeasible, be it because it is the limit of small strains on the one hand or because the sound velocity is infinity so that issues of shock propagation cannot be dealt with.

Linear oscillator with Coulomb friction
The ordinary differential equation (ODE) of a simple linear oscillator with dry friction is not continuous for sign changes ofẋ, which require a special treatment. The mathematical theory goes back to A.F. Filippov [3], while in the numerical treatment, we will follow Hairer and Wanner [4]. The jump in eq. (3) is dealt with by separating the equations formally into with region I for positive and region II for negative velocity. Forẋ = 0, additional formalism is necessary. The are the sums of all forces in region I, respectively II, where for eq. (3) we have used as external force f ext = −kx. The velocityẋ is used as "indicator function" g(x,ẋ) =ẋ to determine whether the system is in region I or II (Fig. 1). The solution atẋ = 0 depends on the direction of the flow of the ODE, which is given by the (signed) scalar product between the velocity and the corresponding acceleration ∇ = ∇ẋ is the directional derivative with respect to the velocity. There are three scenarios for the sign of the a i 's (compare Fig. 1): 1. a I · a II < 0: the friction is dynamic, the flow crosses the x-axis at (v = 0, |x|, x > 0.4 in Fig. 1) 2. a I < 0, a II < 0: the flow pushes into the constraint at Fig. 1) 3. a I > 0, a II > 0: the flow pulls away from the constraint, which usually only occurs for an unfortunate choice of initial conditions (Painlevé paradox, see [7], section 3.3.1) . For case 2, to maintain the constraint, the right hand side of eqs. (7,8) is neither f I nor f II , but a solution in the convex hull, where 0 ≤ λ ≤ 1 is a Lagrange parameter. This means solving Differentiating eq. (11) with respect to time allows to express λ in terms of a I and a II , Using eqs. (7,8) and (in one dimension) ∇g = ±1, the equation further simplifies to This is a remarkable result, as the highly non-linear nonsmooth static friction can be computed by solving a mere linear scalar equation. The linearity in the f ... of eqs. (7,8) allows to write the force as That means that atẋ = 0, static friction compensates the external forces numerically exact, f (ẋ = 0, λ) = 0 (within the error of the solver). The scalar products eqs. (7,8) can be numerically approximated as The δ I , δ II have to be chosen small enough, so that the flow cannot cross the manifold, but large enough, so that the f I , f II are still taken into account. The δ I , δ II are given by

Particle on an inclined slope
The most simple case for a two-dimensional particle system is a square block sliding down a slope. It is a single contact system, with gravity g = 9.81[m/sec 2 ] as external force. We initialise the square particle orientated parallel to the slope without overlap. The slope has an inclination angle of α = π/6 = 30 • , slightly lower than the critical angle of ≈ 31 • for our μ = 0.6, so the block should come to rest. The particle and the slope are modelled as convex polygons. For the interaction we use the centroid of overlap region as the force-point, the normal direction of the contact is the weighted average of the connecting lines between the centroid and the intersection points. The elastic force is proportional to the particle overlap A and the Young's modulus Y and inverse proportional to the inverse characteristic length l. (see Ref. [7] for more details). The damping force is proportional to the change in overlap area In the Cundall-Strack model for Coulomb friction, the tangential force is incremented by and cut off at μF normal . Without cutoff, it integrates out to an harmonic oscillator F = −k t x. Accordingly, as long as there is no dissipation, a particle on a slope arctan(μ) will oscillate around the equilibrium position. Adding as damping term to eq. (19) can lead to permanent slipping close to the critical angle, as the dv is opposite to v and therefore the grip of the tangential force is retarded. Instead, an instantaneous damping ∝ −v is conventionally introduced, which is nevertheless out of phase with eq. (19), so very small prefactors must be used. Physical parameters are given in Tab. 1, and remain the same unless specified otherwise. For CS, one can speed up the grip by using a larger prefactor in the incrementation: Nevertheless, this may make the use of a smaller timestep necessary. Our exact friction implementation (DAE) shows stronger drift far from the critical angle compared to the Cundall-Strack (CS) model (Fig. 2). At the critical angle, the CS model shows a later onset of static friction, i.e. a later "grip", compared to the DAE implementation (Fig. 3), in particular the damped DAE implementation has a significantly earlier onset than any of the other configurations. As can be seen in Fig. 4, the drift decreases with the timestep, a typical accuracy effect in the implementation of DAEs. The difficulty of implementing "numerically exact" static friction for discrete element simulations lies in the fact, that stationary configurations must be obtained by dynamic simulations in the presence of noise, which is a superposition of discretisation error, defective position information, as well as zero-point oscillations of particles around the equilibrium position. Obtaining static equilibria is a much greater challenge than computing mere dynamics: In dynamic computer simulations, bodies may move due to the noise amplitudes alone, while static equilibria have to be obtained within the numerical accuracy of the integrator. The drift increases with the timestep. The implicit Backward-Difference-Formulae (BDF, Gear Predictor-Corrector) show less drift than the explicit Heun-integrator ("2.nd order Runge-Kutta", implemented in predictor-corrector fashion, see [7], p. 80),  with the second order BDF integrator being more accurate than the fifth order integrator. On the other hand, particle size has a much lower impact on the drift for implicit integrators. Fig. 5 shows that only larger particle sizes (>0.1m) cause an increase in the drift for the BDFintegrators. For the Heun-integrator, the drift shows significant increase both for very small and for very large particle sizes. This is consistent with the common practice in numerical analysis that implicit integrators are used for the solution of DAE's due to better stability compared to explicit integrators. If the drift is too large for the time scale considerations, it would be possible to "fix" the particle by velocity projection (see [5], p. 471). This would correspond to setting the velocity zero for a non-moving floor, and to project onto zero relative velocities for moving contacting particles.   Figure 3. Horizontal position for a 1 cm particle on a slope with the critical angle (30 • ) for CS and DAE friction implementations with and without damping.

Many particle systems
We adapt the approach to many-body systems in twodimension in the following way: The constraints are evaluated for each contact. Instead of the single particle mass, the reduced (tangential) mass m eff,tan = (1.0/m 1 + 1.0/m 2 + r 2 1,⊥ /I 1 + r 2 2,⊥ /I 2 ) −1 (21) is computed from the masses m ... and the moments of inertia I ... with the projection r ...,⊥ of the contact vectors (from the centres of mass to the contact point) onto the normal direction. The external forces f ext in eqs. (5,6) are replaced by the relative accelerations, scaled by the effective tangential mass,   Figure 5. Dependence of the numerical drift on the particle size l for dt = 5 · 10 −6 s.
Like solids form in molecular dynamics through equilibration of forces for the inter-particle interactions, the equilibration of the tangential force for each contact leads to many-particle agglomerates in static equilibrium. As no matrix formulation is necessary, no condition number considerations are necessary and no instantaneous propagation of forces occurs. Mathematical convergence proofs are not possible, as the stability depends on the physical stability: even for experimental heaps, vibration can destroy the stability, so it is the entirety of parameters (initial velocities, friction coefficients, damping, state of motion of boundaries, . . . ) which determines whether e.g. a heap is stable or not. For the many particle system, we compare both friction approaches with via a granular heap of 1625 particles formed from a hopper, see Fig. 6). The particles have between 5 and 10 corners, distributed randomly on an ellipse with 2.5 mm mean radius and eccentricity ≈ 0.8 (for more details, see [7], chapter 7), and the same parameters as in section 4. For slopes of ≈ 30 • , the drift is shown in Fig. 7: For the CS-model (equivalent to Baumgartestabilization a DAE-problem [1], i.e. a parametrization of an oscillation around the constraint), the drift is considerably smaller, as it is , while the drift of the DAE-approach depends on the timestep.
1dm Figure 6. Heap with particles dropped from a hopper. Some particles are still falling and appear airborne in the graphic.

Conclusion
We have implemented "numerically exact" Coulomb friction based on a DAE-formalism where the constraints are computed based on the interparticle forces for many body problems in two dimensions, which does not need a matrix formulation, by evaluating the frictional constraints. Compared to the Cundall-Strack model, the "grip" of the numerically exact friction is significantly faster and there is no internal degree of freedom which can release energy in the system. This makes the approach superior for granular systems where realistic damping is necessary, for flu-idisation under vibration (as for earthquakes) or the onset of landslides. Depending on the particle size, the Cundall-Strack model is either close to optimal, as it is equivalent to a Baumgarte-stabilisation of a DAE-constraint, or it fails as the grip acts so slowly that no mechanical equilibrium is reached at all. While the weak point of the Cundall-Strack model is the parameter dependence, a weak point of the DAE-approach in many body configurations is the numerical drift due to the finite time discretiztion accuracy. This is a problem shared by all approaches which compute the friction as equilibrium force for DAE-solutions. Reducing the timestep considerably is not an option, as the simulation times will become too long. If necessary, the drift can be eliminated by stabilisation by projection (in particular velocity projection) [5,6], which for single-particle problems corresponds to setting the velocity, while for many-body problems, the algorithm is more complicated.