# Luding contact model¶

This contact model, published in [Lud08], is suitable for flexible handling and testing of various adhesion/plasticity/viscosity parameters.

This model separately treats 4 contact senses, each including mostly-decoupled elastic-plastic-viscous behavior.

*Normal*sense, where compressive plasticity during loading increases contact stiffness (up to value defined by plasticity limit), unloading/reloading following elastic path. Adhesion effects are proportional to historically maximum compression, thus suitable for modeling stickiness induced by prior compression (e.g. during impact). Further energy is dissipated by viscosity.*Tangential*sense, with Coulomb-style friction with different limits for static/dynamic friction (related via a dimensionless parameter) and viscosity. This sense is 2-dimensional, perpendicular to the normal sense.*Rolling*sense, describing rotation along the same 2 axes as tangential sliding/slipping, with Coulomb-style friction and viscosity.*Twisting*sense, describing rotation along the normal axis, with Coulomb-style friciton and viscosity.

The model is designed in such a way that all tangential, rolling and twisting senses can use the same routine, only with different parameters and different dimensionality: 2d for tangential and rolling, 1d for twisting).

## Formulation¶

The contact exists if there is geometrical overlap (\(\delta\geq0\)) and is dissolved for \(\delta<0\).

### Normal¶

Elastic stiffness \(k_{n2}\) varies based on the historically maximum normal overlap \(\delta_{\rm max}\) and is interpolated between \(k_{n1}\) (for \(\delta_{\rm max}=0\)) and \(\hat{k}_{n2}\) (for \(\delta_{\rm max}\geq\delta_{\rm lim}\)). This is shown in Fig. 29. Algebraically,

with variable stiffness

and the plastic limit overlap (constat for each contact) being defined as

where \(a_1\), \(a_2\) are equivalent radii of contacting particles and

is “reduced radius” of the contact.

Dissipation (work) has two components:

viscous contribution, computed incrementally as

\[\Delta W_{nv}=\gamma_n v_n^2 \Dt,\]plastic contribution, which is the easiest to compute in a total manner when the contact breaks, using \(\delta_{\rm max}\) as history variable and current elastic stiffness \(k_{n2}\equiv k_{n2}(\delta_{\rm max})\); we compute \(\delta_0\) (elastically unloaded zero-force overlap) as auxiliary, and \(\delta_{\rm min}\) to evaluate adhesive contribution (see Fig. 29 which shows the work shaded):

\[ \begin{align}\begin{aligned}\delta_0&=\delta_{\rm max}\left(1-\frac{k_{n1}}{k_{n2}}\right)\\-k_{na}\delta_{\rm min}&=-(\delta_0-\delta_{\rm min})k_{n2}\\\delta_{\min}&=\delta_0\frac{k_{n2}}{k_{n2}+k_{na}}\\W_{np}&=(k_{n1}\delta_{\rm max}+k_{na}\delta_{\rm min})\delta_0\end{aligned}\end{align} \]

Rescaled damping coefficient, typical contact duration and restitution coefficient can be respectively computed as

with reduced mass \(m_{12}=m_1 m_2 /(m_1+m_2)\).

### Tangential¶

Sliding (yield) force is defined in dependence on the normal force (including adhesion) as

trial force (2d vector in tangential plane) being computed from tangential stiffness \(k_t\), elastic displacement \(\vec{\xi}_t\) (2d vector), viscous coefficient \(\gamma_t\) and relative velocity at contact \(\vec{v}_t\) (2d vector) as

where tangential stiffness is defined as \(k_{t}=\hat{k}_{n2}\frac{k_t}{k_n}\).

If \(|\vec{f}_t|\leq f_y^y\), static friction is active (zero plastic dissipation) and

\[ \begin{align}\begin{aligned}\vec{f}_t&=\vec{f}_t^t,\\\nnext{\vec{\xi}}_t&=\pprev{\vec{\xi}}_t+\vec{v}_t\Dt.\end{aligned}\end{align} \]Otherwise, dynamic friction is activated as

\[ \begin{align}\begin{aligned}\vec{f}_t&=\frac{\vec{f}_t^t}{|\vec{f}_t^t|}f_t^y\phi_d,\\\nnext{\vec{\xi}}_t&=-\frac{1}{k_t}\left(\vec{f}_t+\gamma_t\vec{v}_t\right).\end{aligned}\end{align} \]where \(\phi_d\) is the dynamic/static friction coefficients ratio. Plastic dissipation (work) is equal to

\[\Delta W_{tp}=\frac{1}{2}(|\vec{f}_t^t|+f_t^y\phi_d)\frac{1}{k_t}(|\vec{f}_t^t|-f_t^y\phi_d).\]

Viscous dissipation is equal to

### Rolling¶

Rolling and twisting are computed in terms of torques (rather than forces) but the same equations apply; reduced radius \(a_{12}\) is used as length to convert dimensionalities of forces/torques.

Rolling stiffness is \(k_r=\hat{k}_{n2}\frac{k_r}{k_n}a_{12}\), yield torque \(m_r^y=\tan\phi_r(f_n+k_{nc}\delta)a_{12}\). Elastic rotation \(\vec{\xi}_r\) is initially zero, the contact is damped using viscous coefficient \(\gamma_r\). Plastic and viscous work (\(\Delta W_{rp}\), \(\Delta W_{rv}\)) are computed the same as in the tangential sense, with indices replaced.

### Twisting¶

Twisting is identical to rolling, only by demoting 2d-vectors to 1d-vectors (scalars) with indices replaced, i.e. twisting stiffness \(k_w=\hat{k}_{n2}\frac{k_w}{k_n}a_{12}\), yield torque \(m_w^y=\tan\phi_w(f_n+k_{nc}\delta)a_{12}\), elastic twist \(\xi_w\), viscous coefficient \(\gamma_w\). Work contributions are \(\Delta W_{wp}\), \(\Delta W_{wv}\).

## Parameters¶

The following table summarizes all model parameters:

normal | tangent | rolling | twisting | |
---|---|---|---|---|

stiffness | \(\hat{k}_{n2}\) | \(k_t/k_n\) | \((k_r/k_n)a_{12}\) | \((k_w/k_n)a_{12}\) |

normal plasticity | \(\frac{k_{n1}}{\hat{k}_{n2}}\) \(\delta_{\rm lim}\) | – | – | – |

adhesion | \(k_{nc}/k_n\) | – | – | – |

static friction | – | \(\tan\phi_t\) | \(\tan\phi_r\) | \(\tan\phi_w\) |

dynamic friction | – | \(\phi_d\) | \(\phi_d\) | \(\phi_d\) |

viscous coefficient | \(\gamma_n\) | \(\gamma_t\) | \(\gamma_r\) | \(\gamma_w\) |

This table summarizes corresponding attributes in the WooDEM code. Since properties of every comtact are computed from contacting materials, these values are marked as * for averaging, † for minimum and # for averaging with scaling by \(a_{12}\):

normal | tangent | rolling | twisting | |
---|---|---|---|---|

stiffness | `k2hat` |
`ktDivKn` * |
`krDivKn` # |
`kwDivKn` # |

normal plasticity | `k1DivKn` * `deltaLimRel` * |
– | – | – |

adhesion | `kaDivKn` * |
– | – | – |

static friction | – | `tanPhi` † |
`statR` † |
`statW` † |

dynamic friciton | – | `dynDivStat` † |
`dynDivStat` † |
`dynDivStat` † |

viscous coefficient | `viscN` † |
`viscT` † |
`viscR` † |
`viscW` † |

Tip

Ask questions or report issues at github.