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 elastic linear normal force for contact, normal viscosity force, constant normal cohesion force at contact, Coulomb friction tangent force, and rolling resistance moment

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\). When \(d_n > 0\), this force is zero. The value of \(k_n\) can be set with the keyword knContact

2. The viscuous 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 viscuous damping. There are two ways to set the value of \(\alpha_n\):

1. The keyword en2Contact that is the energy normal restitution rate (\(e_n^2\)). In this case, the viscuous damping rate will be set associated

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

2. The keyword en2ContactFromViscRate that will

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

Note

Because of the viscuous part of the normal force, the total normal force \(f_n = f_n^\text{elas} + f_n^\text{visc}\) can be negative. Although the physical meaning of this negative value is arguable, the normal force is by default restricted to remain positive or zero.

2. Tangential component

The tangential force represents the frictional resistance between particles when they slide against each other. It is calculated based on the relative tangential velocity (\(v_t\)) and a tangential stiffness parameter (\(k_t\), keyword ktContact). The Coulomb friction model is used to limit the tangential force, ensuring that it does not exceed the product of the friction coefficient (\(\mu\), keyword muContact) and the normal elastic force (\(f_n^\text{elas}\)). The tangential force is incrementally updated with the following relation, at each time step \(\Delta t\), since the onset of contact:

\[\Delta f_t = \left [ k_t v_t \Delta t \right ]_{\pm \mu f_n^\text{elas}}\]

The friction force \(f_t\) is cancelled as soon as contact is lost.

3. Moment vector

The resistant moment (\(\underline{M}\)) accounts for rotational resistance between particles. It is calculated based on the rotational velocity difference (\(\underline{\omega}_i - \underline{\omega}_j\)) and the rotational stiffness parameter (\(k_r\), keyword krContact). Additionally, the moment is limited by the product of the moment length coefficient (\(\mu_r\), keyword murContact), the normal force (\(f_n\)), and the length of the branch vector (\(\ell\)). The resistant moment is updated, similarly to the friction force, as follows:

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

Note

The moment resistance is applied only when \(k_r > 0\).

Law for rock avalanches (keyword Avalanche)

This is historically the first law that has been implemented in Rockable (actually in DEMbox, its ancestor).

The force-law implemented is specifically designed for simulating rock avalanches at Laboratoire 3SR. It governs the interactions between particles or bodies in a granular material, considering both normal and tangential forces as well as resistant moments due to rotational motion.

1. Normal component

   The normal force (I.fn) between particles is calculated to represent the contact forces that resist compression and prevent interpenetration. The calculation involves considering the change in the normal displacement (I.dn) and the normal stiffness parameter (kn). The normal force may also undergo unloading or loading depending on the change in normal displacement.

2. Tangential component

   The tangential force (I.ft) represents the frictional resistance between particles when they slide against each other. The calculation considers the tangential velocity (vt) and a tangential stiffness parameter (kt). Frictional forces are limited by the friction coefficient (mu) and the normal force (I.fn).

3. Resistant moment

   The resistant moment (I.mom) accounts for the rotational resistance between particles. The calculation considers the rotational velocity difference (box->Particles[I.j].vrot - box->Particles[I.i].vrot) and the rotational stiffness parameter (kr). Additionally, there is a correction applied to the branch vector (branch) to account for different scenarios when the free body is either particle i or particle j.

Note

Weighted Interaction Parameters. The force-law allows for adjusting the interaction parameters (kn, kt, kr) based on a weight factor (w). This weighting is determined by the function box->ctcPartnership.getWeight, which can influence the forces and moments between the particles.

Note

The code includes conditional blocks for certain options (FT_CORR) which may affect the actual computations. Additionally, the equations use specific parameters (en2 and mur), which should be defined and obtained from relevant data sources in the simulation setup.

Breakable elastic solid bonds (keyword StickedLinks)

This law has been introduced since the PhD work of Marta Stasiak.

Two states can exist: cohesive bond and contact.

1. Cohesive bond

   When a cohesive bond exists between two particles, the force-law employs the following equations:

   1. Normal component

      The normal force (\(f_n = f_n^\text{elas} + f_n^\text{visc}\)) between the particles is calculated using Hooke’s elastic contact, combined with viscous damping:

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

      where \(d_n^0\) is the initial normal distance. This time, the parameter \(k_n\) is set with the keyword knInnerBond or knOuterBond. The viscous normal force is:

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

      where \(\alpha_n\) is indirectly set with the keyword en2InnerBond or en2OuterBond.

   2. Tangential component

      The tangential elastic force (\(f_t\)) is incrementally calculated based on the relative tangential velocity (\(v_t\)) and a tangential stiffness parameter (\(k_t\), keyword ktInnerBond or ktOuterBond):

      \[\Delta f_t = k_t (v_t \Delta t)\]

   Optionally, an additional term for tangential viscosity (\(f_t^\text{visc} = 2 \alpha_n \sqrt{k_n m_\text{eff}} \, v_t\)) can be added to the elastic part of the tangential force. However, this approach accumulates viscosity over time and may not be accurate:

   3. Moment vector

      The moment carried by the bond is accumulated incrementally from the relative rotation velocity, with the rolling stiffness \(k_r\) (keyword krInnerBond or krOuterBond):

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

      Unlike the resistant moment of a contact, this moment is not capped: as long as the bond exists it grows freely, and it is the rupture criterion below that limits it through \(M_0\).

   4. Rupture criterion

      The rupture criterion determines if the cohesive bond should break. If the yield value is greater than 0, the bond is irreversibly broken (actually, the breakage of the whole interface is scheduled as soon as possible). The yield function reads:

      \[\varphi = \left( \frac{\vert f_t \vert}{f_t^0} \right)^p + \left(\frac{\Vert\underline{M}\Vert}{M_0}\right)^p - \frac{f_n}{f_n^0} - 1\]

2. Contact

   When particles are in contact (no cohesive bond), the force-law employs the following equations:

   1. Normal component

      The normal force (I.fn) between the particles is calculated similar to the cohesive bond case, but with different interaction parameters (kn, kt, damp, etc.):

      \[ \begin{align}\begin{aligned}\text{vn} = \mathbf{I.vel} \cdot \mathbf{I.n}\\\text{fne} = -\text{kn} \cdot \text{I.dn}\\\text{fnv} = \text{damp} \cdot \text{vn}\\\text{I.fn} = \text{fne} + \text{fnv}\end{aligned}\end{align} \]

   2. Tangential component (Friction)

      The tangential force (I.ft) due to friction between the particles is calculated in a similar manner to the cohesive bond case but with different interaction parameters (kt, mu, etc.).

   3. Resistant moment vector

      If the contact allows for rotational motion (kr parameter), a resistant moment (I.mom) is calculated based on the relative rotational velocity of the particles.

Note: The code includes conditional blocks for certain options (FT_CORR and box->ctcPartnership.getWeight != nullptr) which may affect the actual computations.

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 stiffnesses of an interface are shared between the sub-bonds it holds. If the interface carries \(N_b\) bonds, each of them uses

\[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 values of \(k_n\) and \(k_t\) come from knInnerBond and ktInnerBond, or from the interface itself when ParamsInInterfaces is 1.

The normal force combines the elastic response, measured from the distance \(d_n^0\) recorded at gluing, with the viscous damping:

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

and the tangential force is accumulated incrementally, then completed by a viscous term:

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

Warning

As in StickedLinks, the tangential viscous term is added to the accumulated force rather than to its elastic part alone, so the viscosity piles up over time. Keep this in mind when interpreting the tangential response of a bond.

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} \frac{1}{2} k_n \left( d_n^b - d_n^0 \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.

Note

The stiffnesses used in this energy balance are the ones read from the group-pair table, not the per-bond values divided by \(N_b\) that the force computation uses.

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.

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/kt/krInnerBond, kn/kt/krOuterBond	no	no	yes	inner only	no
fn0/ft0/mom0/pow*Bond, en2*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.