Overlap Computation
5.3 Modeling of the tangential Force
The need to model the tangential force arise from the difficulty to implement Coulomb fric-tion numerically exact for static fricfric-tion, a problem not only limited to DEM simulafric-tions of granular material, but a general problem for numerical approaches for the dynamics of many-body systems. While many-body systems of rigid particles with Coulomb fric-tion forces have been dealt with some time ago via “contact dynamics” [60] (see secfric-tion 1.4.3), for the case where the Coulomb friction forces are reactions to elastic normal forces with velocity constraints of many-body systems, a general solution is still outstanding. In
112 FORCE MODELING
granular research, friction is usually “approximated” with “breaking spring” model dating back to Cundall and Strack [24] who first introduced the discrete element model (DEM) for granular material research, or worse, by “fudging” a velocity-dependent viscous-like damping (e.g. Chong et al. [110]), or, finally it is totally neglected [111]. So far, for nu-merical exact solutions for static friction, for one- and two-dimensional problems, solving the static friction as non-holonomic constraint has been discussed [74, 68, 112]. For three dimensions, there are still open questions, due to foreseeable difficulties, among them the vanishing of the uniqueness of the tangential direction. While in one or two dimensions, the external forces between particles can be compensated uniquely by the static friction force at the contact, in three dimensions the force directions may be skewed and thus no unique direction for the friction exists a priori. Instead of seeking the exact solution for Coulomb friction, we generalize the Cundall-Strack friction model, which worked well in two dimensions [106, 108], to three dimensions. In this section, we explain how the Cundall and Strack’s friction model works in two dimensions first and then how we adapt and implement it in three dimensions for our current DEM code. The caveats for using such models are also discussed.
5.3.1 Cundall-Strack friction in two dimensions
t
t stopped
2
5
contact
t sliding t approach
t sliding
0 1
3 4
t sliding
Fig. 5.15 Deformation of the “tangential spring” in Cundall-Strack friction model in a contact process, the friction in the model is approximated as a spring in the tangential direction.
In two dimensions, the tangential contact force introduced by Cundall and Strack is modeled incrementally [24]: at timet, the magnitude of the tangential forceftis given by
||ft(t)||=||ft(t−dt)||+kt||vt|| ·dt, (5.10)
MODELING OF THE TANGENTIAL FORCE 113
in which kt is the tangential stiffness, vt the tangential velocity; the direction of ft is in the direction opposite tovt; additionally, the magnitude of the tangential force is checked against the maximal possible valueµ||fn||to apply the rule
ft(t) =sgn(ft(t))·µ||fn(t)|| if||ft(t)||> µ||fn||, (5.11) where sgn is the sign function, µ the friction coefficient and fn the normal contact force given by Eq. (5.7). Since the direction of the tangential force is obtained from the tangen-tial velocity, the Cundall-Strack model uses scalar increment of magnitude. As kt must have the dimension of a spring constant, the Cundall-Strack model is sometimes referred to as a model of breaking tangential springs.
µ f n
t
t t
t t
t
static friction
0 1
2 3
4 5
t
Cundall−Strack friction force dynamic friction
µ f n
Fig. 5.16 Intended behavior for the Cundall-Strack friction (thick black dashed line) and actual oscillatory behavior (thin black full line).
The behavior we want to obtain with this model is the one in Fig. 5.16, a monotonous increase of the tangential force, until saturation occurs at the theoretical value. Unfortu-nately, the behavior will not be monotonous: If we divide Eq. (5.10) bydt, we obtain (in scalar form)
df/dt=−ktvt. (5.12)
If we integrate this equation over dt, we see that it is essentially the equation for the harmonic oscillator. That means the tangential force does not always act strictly opposite to the actual velocity, but due to the inertia of the “harmonic oscillator”, it may due to delay, even act in the direction of the actual velocity. Only if the tangential friction reaches the value for sliding friction (condition before Eq. (5.11)), energy is dissipated. Without additional damping, the tangential force will lead to oscillatory behavior (see thin full line in Fig. 5.16), only additional damping will lead to a purely monotonous curve (black
114 FORCE MODELING
dashed line in Fig. 5.16). While the damping may reduce the oscillations, one also sees that the Cundall-Strack model leads to a tangential friction whose action is delayed in comparison to the “exact” friction.
v (t−dt) t
|f (t−dt)| t
t
|f (t−dt)|
t
v (t)
t
f (t)
−kv (t) dt
t
Fig. 5.17 In the two-dimensional case, the tangential friction from the previous time-step ft(t−dt) is used irrespective of a possible shift in the direction at the new time-step.
As shown in Fig. 5.17, we can interpret the directions via vector projection: the tangen-tial force from the previous time-step ft(t−dt) is projected onto the tangential direction of the current time-step along the current contact velocity vt(t) while its magnitude is kept by rescaling the projection of ft(t−dt); Then the new tangential forceft is obtained by the vector addition of the incremental force −ktvt during the time intervaldt and the magnified projection of previous tangential forceft(t−dt).
5.3.2 Cundall-Strack friction in three dimensions
To generalize the Cundall-Strack model to three dimensions, several issues should be ad-dressed. In three dimensions, the manifold of the tangential movement is two-dimensional (r1(t), r2(t)), while the actual trajectory of the contact point is three dimensional (x(t), y(t), z(t)). Therefore, we need a way to relate the contact manifold(r1(t), r2(t))with the con-tact trajectory (x(t), y(t), z(t))in a unique way. At the same time, we want to retain the incremental feature of the Cundall-Strack model in two dimensions that the magnitude of the contact force at the previous timestep is used irrespective of a shift in the direction.
The three concepts, projection, rescaling and vector addition, are crucial to adapt the Cundall-Strack model for three dimensional contacts:
1. Projection: During the advance from time-stept−dt to time-stept, where we have the new contact normal n(t)ˆ and the new tangential velocity vt(t), we project the old
MODELING OF THE TANGENTIAL FORCE 115
tangential forceft(t−dt) onto the new tangential plane,
ft(t−dt)p=ft(t−dt)−(ft(t−dt)·n(t)) ˆˆ n(t). (5.13) 2. Rescaling: We then rescale the projection ft(t−dt)p by the magnitude of the previous tangential force||ft(t−dt)||,
ft(t−dt)r=||ft(t−dt)|| · ft(t−dt)p
||ft(t−dt)p||. (5.14) 3. Vector addition: the rescaled projectionft(t−dt)r is then added by the incremental vector for the intervaldt,
ft(t) =ft(t−dt)r−ktvt(t)dt. (5.15) A cut-off is applied, if the results from the vector addition exceeds the maximal friction allowed (the dynamic friction)
ft(t) =sgn(ft(t))·µ||fn(t)|| if||ft(t)||> µ||fn||, (5.16) to obtain the tangential force as dynamic friction.
Applying the above scheme, we can obtain a three-dimensional vector ft(dt), which is in the tangential plane of the contact and of incremental property, which will be used as an approximation of friction in our DEM simulations. A test case of a block on an inclined plane will be used in section 5.4 as a validation for the generalization of Cundall-Strack friction in three dimensions.
In practice, the incremental term−ktvt(t)dtin Eq. (5.15) deserves more consideration:
similar to the elastic force model, where the one-dimensional “spring constant“ k (with unit [kg/s2] or equivalent to “force/length”) of the harmonic oscillator is replaced by a
“two-dimensional” modulusY (with unit [kg/(m·s2)] or equivalent to “force/area”) in our force model, for the tangential stiffness constant in a three dimensional problem, we need the characteristic length Lc (Eq. ( 5.3)) to fix the unit for the tangential force increment
∆ft=Yt·Lcvt(t)dt, (5.17) when using a “tangential Young’s modulus” Yt [106] instead of a spring constant kt. The tangential Young’s modulus Yt is often chosen asYt=νY (0< ν <1), in which ν is func-tionally an analogue to Poisson’s ratio, as a measurement of a particle in force equilibrium (which leads to static friction) to resist the change of position along the tangential direc-tion (since Y is used as a measure to resist the change along normal direction). In our DEM simulations,Yt= 2/7Y is used as in our two dimensional simulations for polygonal particles [25, 106, 108]. IfYt is used instead of kt in Eq. (5.15) and the length dimension
116 FORCE MODELING
Lc is missing, for cases that the characteristic length are small (Lc<1, which are usually the cases in our DEM simulations with small particles), we will have a large increment term
Ytvtdt >(Ytvtdt)Lc, (5.18) which may lead to the “incremental” mechanism (Eq. (5.15)), which is the essence of the model, dysfunctional. For example, for a simulation with Lc is of the order of 10−3, the increment tangential force from Eq. (5.15) would be 100 to 1000 times larger than the normal increment ∆ft in Eq. (5.17). Due to the cut-off in Eq. (5.16), the tangential force would be maximal almost all the time in simulation, which means the “incremental”
mechanism is not working and the model is equivalent to dynamic friction, neglecting the static friction all the time.
5.3.3 Dissipative force in tangential direction
The damping force in tangential direction for the polygonal particles in two dimensions [106] takes the following form,
ftd =−γt
√
YtMredt vt, (5.19)
in which theMredt is the reduced “tangential mass”
Mredt = 1
1/mi+ 1/mj+r2i/Ii+rj2/Ij, (5.20) which takes into account the moments of inertia Ii and Ij of the particles. To adapt it to three dimensions, we need to consider two things, first the dimension of the equation and second the form of the reduced “tangential mass”. Since Ytis with unit [kg·s−2·m−1] compare to [kg·s−2] in two dimensions, the characteristic length Lc enters the tangential damping force in three dimensions
ftd=−γt
√
YtMredt Lcvt. (5.21)
Mredt is the reduced “tangential mass” in three dimensions, as the moments of inertia of the particles in three dimensions would be a matrix rather than scalar quantity in Eq. (5.20).
We consider the tangential velocity at the contact point for two particles iand j
vt= [v1+ω1×r1−(v2+ω2×r2)]t. (5.22) Its time derivative gives the tangential acceleration
v˙ = ˙v1+ ˙ω1×r1+ω1×r˙1−( ˙v2+ ˙ω2×r2+ω2×r˙2), (5.23)
MODELING OF THE TANGENTIAL FORCE 117
where the subscript “t” which indicates the tangential direction is dropped. Substituting r˙i=ωi×ri into Eq. (5.23) and rearranging the terms gives
v˙ = ˙v1−v˙2+ ˙ω1×r1−ω˙2×r2+ω1×(ω1×r1)−ω2×(ω2×r2), (5.24) whereωi×(ωi×ri)would be vanished for round particles since it is along theri direction and for round particles it has no component in the tangential direction. We denoteA = ω1×(ω1×r1)−ω2 ×(ω2 ×r2) since it does not contain the force terms explicitly and substitute
˙
ωi =I−i 1((Iiωi)×ωi+ri×fi) into Eq. (5.24) yields
v˙ = ˙v1−v˙2+(
I−11(r1×f1))
×r1−(
I−21(r2×f2))
×r2+A+B, (5.25) where
B=(
I−11(I1ω1)×ω1
)×r1−(
I−21(I2ω2)×ω2
)×r2
is another term which does not contain the force explicitly. Since the two particles are in contact and we ignore the influences from the other particles as well as the gravity, then we have f2 =−f1 and the Eq. (5.25) becomes
v˙ = ( 1
m1 + 1 m2
) ft+ (I−11(r1×ft))
×r1+(
I−21(r2×ft))
×r2+ (I−11(r1×fn))
×r1+(
I−21(r2×fn))
×r2+A+B,
(5.26)
whereft and fn are the tangential and the normal components of f1. For round particles, the item(
I−i 1(ri×fn))
×rican be dropped sinceriis parallel tofn, while for polyhedral (or polygonal) particles usually it can not. The reduced tangential mass for two dimensional problems, Eq. (5.20), comes from the simplification of Eq. (5.26) for two dimensional discs
v˙ = ( 1
m1
+ 1 m2
+ r21 I1
+r22 I2
) ft+B
= 1
Mredt ft+B,
(5.27)
which could be treated as a physical harmonic oscillator with massMredt for whose viscous damping force would be
fd=−√
Mredt ktv.
Since it is difficult to decouple the tangential and the normal component of the contact force in Eq. (5.26) and write in the form like a harmonic oscillator along the tangential direction alone as we can do for two dimensional discs in Eq. (5.27), we did not proceed to derive for the reduced “tangential mass” for polyhedral particles. Moreover, the derivation is only exact for two-particle problems, while most of the time, our particles will have more
118 FORCE MODELING
contacts. Instead, we use the same value as for the dissipative force in normal direction Mredt = 1
1/mi+ 1/mj
. (5.28)
for polyhedral particles Pi andPj in our DEM code. Compared with its two dimensional counterpart in Eq. (5.28) it may overdamped since the its effect on the rotary degrees of freedom have not been taken into account. Cundall and Strack suggested to choose the same value for the normal damping constant γn and the tangential damping constant γt [24]. As an compensation for our overdamping, we could take the damping constant for tangential direction half of that for the normal direction. In practice, this is not necessary, since we don’t want the oscillation along the tangential direction anyway, so an overdamped system would be preferable. What really should be taken into consideration is whether the choice of “mass” and “stiffness” in the damping model would not increase the damping, so that the corresponding system of equations would become stiff. This would make the timesteps in explicit integrators prohibitively small, but it will not a problem with our BDF-integrators. The total force in tangential direction would then be
ft = fet+fd, (5.29a)
ft = sgn(ft)·µ||fn|| if||ft||> µ||fn||, (5.29b) where fet is the tangential force from Eq. (5.15) or its cut-off, Eq. (5.16).
5.3.4 Caveats
In the Cundall-Strack model [24] and its generalization [113, 114], the tangential force is incremented from the previous time step, with some additional modifications for the direction, and there is a cut-off due to the normal force and the friction coefficient. Such models are widely used and are believed to yield reasonable results when the velocity variations in the system are not large, e.g. in shearing a granular system or in dynamic process like heap formation etc. There are several drawbacks of a friction model compared to the (currently not accessible) analytic solution. While actual solid friction is practically instantaneous, there is a delay in the Cundall-Strack model which depends actually on the spring constant. This delay is smaller, the largerktin eq. (5.12) is chosen. Nevertheless, a large choice ofktnecessitates a reduction of the timestepdtfor the time integration, which else depends only on Young’s modulus and the particle mass, i.e. the necessary amount of computer time increases for “better” modeling. Another problem is that the “tangential spring” is an unphysical degree of freedom which may randomly store and release energy and therefore act as a noise term. Especially for strongly oscillatory problems with high amplitudes, the consequences for the verisimilitude of the simulation are difficult to fathom.
A SIMPLE TEST: A CUBE ON AN INCLINED PLANE 119