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.

1. 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.

2. 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.

3. Rolling sense, describing rotation along the same 2 axes as tangential sliding/slipping, with Coulomb-style friction and viscosity.

4. 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. 28. Algebraically,

\[\begin{split}f_n=\underbrace{\gamma_n v_n}_{\hbox{viscosity}}+\begin{cases}k_{n1}\delta & k_{n2}(\delta-\delta_0)\geq k_{n1}\delta\quad\text{(plastic loading)} \\ -k_{nc}\delta & -k_{nc}\delta\geq k_{n2}(\delta-\delta_0)\quad\text{(adhesive)} \\ k2(\delta-\delta_0)&\text{otherwise (in-between unloading/reloading)} \end{cases}\end{split}\]

with variable stiffness

(27)\[\begin{split}k_{n2}(\delta_{\rm max})&=\begin{cases} \hat{k}_{n2} & \delta_{\rm max}\geq\delta_{\rm lim} \\ k_{n1}+(\hat{k}_{2}-k_{n1})\frac{\delta_{\rm max}}{\delta_{\rm lim}} & \delta_{\rm max}<\delta_{\rm lim} \end{cases}\end{split}\]

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

\[\delta_{\rm lim}=\frac{\hat{k}_{n2}}{\hat{k}_{n2}-k_{n1}}\phi_f a_{12}.\]

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

\[a_{12}=\frac{2a_1 a_2}{a_1+a_2}\]

is “reduced radius” of the contact.

Fig. 28 Diagram showing contact behavior in the normal sense: plastic loading with stiffness \(k_{n1}\), elastic unloading/reloading with stiffness \(k_{n2}(\delta_{\rm max})\), adhesion with \(k_{na}\) stiffness.

Dissipation (work) has two components:

1. viscous contribution, computed incrementally as

   \[\Delta W_{nv}=\gamma_n v_n^2 \Dt,\]

2. 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. 28 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

\[ \begin{align}\begin{aligned}\eta_0&=\frac{\gamma_n}{2m_{12}},\\t_c&=\frac{\pi}{\omega}=\frac{\pi}{\sqrt{k_n/m_{12}-\eta_0^2}},\\r&=v_n'/v_n=\exp(-\pi \eta_0/\omega)=\exp(-\eta_0 t_c)\end{aligned}\end{align} \]

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

\[f_t^y=\tan\phi_t(f_n+k_{nc}\delta),\]

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

\[\vec{f}_t^t=-k_t\pprev{\vec{\xi}}_t-\gamma_t\vec{v}_t\]

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

1. 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} \]

2. 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

\[\Delta W_{tv}=\gamma_v|\vec{v}_t|^2\Dt.\]

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

Report issues or inclarities to github.