Force laws

We provide here the interaction models existing in Rockable. The equations are explained with the associated parameters. The parameters [PARAMETER] are to be defined in the input conf-file in the following way:

  • [PARAMETER] (int) firstGroupNumber (int) secondGroupNumber (double) value

For example, to set the value \(10^8\) to the normal stiffness \(k_n\) between the elements that belong to the group number 0 and the elements that belong to the group number 2, we will write:

input.txt
      ...
knContact 0 2 1e8
      ...

Default model (keyword Default)

This is the force law to be used in most cases. In short, the model includes a linear elastic normal force, a normal viscous damping, a Coulomb friction, and a rolling resistance moment. There is no cohesion: the interaction exists only while the bodies overlap (\(d_n \leq 0\)), and all its forces are reset to zero as soon as the contact is lost.

  1. Normal component

  1. The elastic part of the normal contact force is

\[f_n^\text{elas} = -k_n d_n\]

where \(d_n \leq 0\) is the normal distance, so \(f_n^\text{elas} \geq 0\). The value of \(k_n\) is set with the keyword knContact.

  1. The viscous part of the normal contact force is

\[f_n^\text{visc} = 2 \alpha_n \sqrt{k_n m_\text{eff}} \, v_n\]

where \(v_n\) is the relative normal velocity, \(m_\text{eff}=(m_i m_j)/(m_i+m_j)\) is the effective mass, and \(\alpha_n \in [0, 1[\) is the rate of normal viscous damping. There are two ways to set the value of \(\alpha_n\):

  1. The keyword en2Contact, which is the squared normal restitution coefficient (\(e_n^2\)). The viscous damping rate is then

\[\alpha_n = \frac{- \ln e_n}{\sqrt{\ln^2 e_n + \pi^2}}\]

  1. The keyword en2ContactFromViscRate, which gives \(\alpha_n\) directly; the corresponding restitution is

\[e_n^2 = \exp \left(-\frac{2 \alpha_n \pi}{\sqrt{1 - \alpha_n^2}}\right)\]

Note

Because of the viscous part, the total normal force \(f_n = f_n^\text{elas} + f_n^\text{visc}\) could be negative. It is restricted to remain positive or zero.

Note

With preventCrossingLength \(L > 0\), the elastic part is multiplied by \(L / (L + d_n)\), which stiffens the repulsion at large overlap and prevents the bodies from crossing each other. The factor is capped at 100. This holds for Default and GeoVisc.

  1. Tangential component

The tangential force is accumulated incrementally from the relative tangential velocity \(\underline{v}_t\), with the tangential stiffness \(k_t\) (keyword ktContact), and projected onto the Coulomb cone defined by the friction coefficient \(\mu\) (keyword muContact) and the elastic part of the normal force:

\[\underline{f}_t \leftarrow \underline{f}_t + k_t \underline{v}_t \Delta t, \qquad \left\Vert \underline{f}_t \right\Vert \leq \mu f_n^\text{elas}\]

Compiling with ALLOW_NEGATIVE_THRESHOLDS uses the total normal force \(f_n\) instead of \(f_n^\text{elas}\) in the threshold.

  1. Moment vector

When the rolling stiffness \(k_r\) (keyword krContact) is strictly positive, a resistant moment is accumulated from the relative rotation velocity, and capped by the rolling resistance \(\mu_r\) (keyword murContact), which is a length, times the normal force:

\[\underline{M} \leftarrow \underline{M} + k_r \left( \underline{\omega}_j - \underline{\omega}_i \right) \Delta t, \qquad \left\Vert \underline{M} \right\Vert \leq \mu_r f_n\]

Note

When \(k_r = 0\), no moment is computed.

Note

When a contact partnership is active (keyword ContactPartnership), the stiffnesses \(k_n\), \(k_t\), \(k_r\) of each sub-contact are multiplied by its weight \(w\), and the damping coefficient by \(\sqrt{w}\).

Law for rock avalanches (keyword Avalanche)

This is historically the first law that has been implemented in Rockable (actually in DEMbox, its ancestor), to simulate rock avalanches at Laboratoire 3SR. It differs from Default by the way the energy is dissipated at the contacts: there is no viscous damping, but a hysteretic normal force.

  1. Normal component

    The normal force is computed incrementally while the bodies get closer (loading), and follows a softer unloading branch when they separate:

    \[\begin{split}f_n \leftarrow \begin{cases} f_n - k_n \, \Delta d_n & \text{if } \Delta d_n < 0 \text{ (loading)} \\ - e_n^2 \, k_n \, d_n & \text{if } \Delta d_n > 0 \text{ (unloading)} \end{cases} \qquad f_n \geq 0\end{split}\]

    where \(\Delta d_n\) is the change of the normal distance since the previous step. The unloading stiffness is \(e_n^2 k_n\), so that the ratio of the energy given back to the energy stored is \(e_n^2\) (keyword en2Contact). The force is kept positive or zero.

  2. Tangential component

    As in Default, accumulated incrementally with \(k_t\) and capped by the Coulomb cone, here with the total normal force: \(\Vert \underline{f}_t \Vert \leq \mu f_n\).

  3. Resistant moment

    Accumulated as in Default, \(\underline{M} \leftarrow \underline{M} + k_r (\underline{\omega}_j - \underline{\omega}_i) \Delta t\) (even when \(k_r = 0\), which then gives no moment), and capped by

    \[\left\Vert \underline{M} \right\Vert \leq \mu_r \, \ell \, f_n\]

    where \(\mu_r\) (keyword murContact) is here dimensionless, and \(\ell\) is the distance from the contact point to the rotation axis of the free body (a rock block), that is, the length of the branch vector from its centre to the contact point once its component along the rotation vector is removed. When body \(i\) is driven, the free body is \(j\); otherwise it is \(i\). When the free body does not rotate at all, the whole branch vector is used.

Note

The contact partnership weights \(k_n\), \(k_t\) and \(k_r\) as in Default. FT_CORR applies the objectivity correction to the tangential force, as in the other laws.

Bonded Cell Method (keyword BCM)

This law was introduced to model a cohesive solid discretised into cells that share faces, rather than a packing of grains glued at a few points. The interface between two cells carries an area, and the rupture is governed by a fracture energy rather than by a force threshold.

Important

BCM only makes sense once the interfaces have been created by the pre-processing command stickBCM, which is the command that computes the area of each interface (see Pre-processing commands). All the interfaces built that way are flagged as inner, so BCM reads only the *InnerBond parameters.

Cohesive bond

The values of \(k_n\) and \(k_t\) come from knInnerBond and ktInnerBond, or from the interface itself when ParamsInInterfaces is 1. By default, they are the stiffnesses of each bond. With bondStiffnessPerInterface 1, they are the stiffnesses of the whole interface, shared between its \(N_b\) bonds:

\[k_n \leftarrow \frac{k_n}{N_b} \qquad k_t \leftarrow \frac{k_t}{N_b}\]

so that refining the discretisation of an interface does not stiffen it. The damping coefficient of the bonds and the breakage energy below use the same weighted stiffnesses. For the cubic cells of Deformable particles (DCEM), which have 4 bonds per face, knInnerBond 0 0 1e5 per bond and knInnerBond 0 0 4e5 per interface give the same result.

The normal force combines the elastic response, measured from the distance \(d_n^0\) of the bond when it was created by stickBCM (each bond has its own, so that a bond carries no force in its initial position), with the viscous damping:

\[f_n = -k_n \left( d_n - d_n^0 \right) + c\, v_n\]

and the tangential force is the sum of an elastic part, accumulated incrementally, and of a viscous part:

\[\underline{f}_t^e \leftarrow \underline{f}_t^e + k_t\, \underline{v}_t \Delta t \qquad \underline{f}_t = \underline{f}_t^e + c\, \underline{v}_t\]

The viscous part of the previous step is removed before the elastic increment is added, so that the viscosity is not cumulated over time.

Neither the tangential force nor the moment is capped while the bond holds: the bond is elastic until it breaks.

Rupture criterion

The criterion is energetic and applies to the interface as a whole, not to its individual bonds. At each step, the elastic energy stored in the interface is summed over its bonds. Only the bonds in tension contribute to the normal part:

\[E_n = \sum_{b\,\mid\, d_n^b > d_n^{0,b}} \frac{1}{2} k_n \left( d_n^b - d_n^{0,b} \right)^2 \qquad E_t = \sum_{b} \frac{1}{2} k_t \left\Vert \underline{s}_b \right\Vert^2\]

where \(\underline{s}_b\) is the accumulated tangential sliding of the bond. The whole interface breaks as soon as

\[E_n + E_t > 2\, A\, G_c\]

with \(A\) the area of the interface and \(G_c\) the fracture energy (keyword gcInnerBond, or gcOuterBond for an outer interface). The factor 2 accounts for the two surfaces created by the crack.

The stiffnesses and \(G_c\) are those of the force computation: per bond or per interface (bondStiffnessPerInterface), and taken from the interface itself when ParamsInInterfaces is 1.

The area, and the \(d_n^0\) and \(\underline{s}_b\) of each bond, are saved in the conf-files (see Format of configuration files (conf-file)), so that a computation restarted from a conf-file breaks its interfaces at the same time as the uninterrupted one.

When the criterion is met the interface is not broken immediately: it is added to a list and all its bonds are released once every force has been computed. This postponement avoids introducing an asymmetry in the incrementally computed tangential forces.

Contact

Once a bond is broken, or between two cells that were never glued, the interaction falls back to a law identical to Default: linear elastic normal repulsion with viscous damping and a positivity clamp, Coulomb friction capped at \(\mu f_n\), and a resistant moment capped at \(\mu_r f_n\). The parameters are then knContact, ktContact, muContact, krContact and murContact, weighted by the contact partnership when one is active.

The interface damage accumulated over a simulation can be followed with the TrackDamage data extractor (see Data extractors).

Slow landslides (keyword GeoVisc)

GeoVisc is a variant of Default designed for slow landslides, where the relevant resistance is viscous rather than elastic-frictional. Its specific feature is to make the friction force viscous while still yielding on the Coulomb cone.

Normal component

Identical to Default: linear elastic repulsion plus viscous damping, with the total clamped to remain positive.

\[f_n = \left[ -k_n d_n + c\, v_n \right]_{\geq 0}\]

The keywords are knContact and en2Contact.

Tangential component

This is where the law departs from Default. The tangential force is not accumulated incrementally from a tangential stiffness; it is proportional to the current sliding velocity, through a tangential viscosity \(\eta_t\) (keyword viscTContact):

\[\underline{f}_t = \eta_t\, \underline{v}_t\]

It is then projected back onto the Coulomb cone, using the elastic part of the normal force as the reference:

\[\left\Vert \underline{f}_t \right\Vert \leq \mu \left| f_n^\text{elas} \right|\]

Note

Using \(f_n^\text{elas}\) rather than the total \(f_n\) keeps the threshold free of the viscous fluctuations of the normal force. Compiling with ALLOW_NEGATIVE_THRESHOLDS switches the reference to the total \(f_n\).

Because the tangential force has no memory, ktContact is not used by this law: a body at rest carries no tangential force, and the model cannot sustain a static shear. It is meant for creeping motion, not for equilibrium. The exercise 4 of the tutorial examples/columnCollapse shows it: a deposit that would stop with Default keeps creeping with GeoVisc.

Moment vector

As in Default, and only when \(k_r > 0\):

\[\Delta \underline{M} = \left[ k_r \left( \underline{\omega}_j - \underline{\omega}_i \right) \Delta t \right]_{\pm \mu_r f_n}\]

The keywords are krContact and murContact.

Summary of the parameters used by each law

Parameter

Default

Avalanche

StickedLinks

BCM

GeoVisc

knContact, en2Contact

yes

yes

yes

yes

yes

ktContact

yes

yes

yes

yes

no

muContact, krContact, murContact

yes

yes

yes

yes

yes

viscTContact

no

no

no

no

yes

kn/ktInnerBond, kn/ktOuterBond

no

no

yes

inner only

no

krInnerBond, krOuterBond

no

no

yes

no

no

en2InnerBond, en2OuterBond

no

no

yes

inner only (damping set by stickBCM)

no

fn0/ft0/mom0/pow*Bond

no

no

yes

no

no

gcInnerBond, gcOuterBond

no

no

no

yes

no

A parameter that a law does not read may be left undefined in the conf-file. Conversely, every pair of groups that can come into contact must have the parameters of the selected law defined, otherwise they default to zero.