Penalty approximation of unilateral contact (source code)

= Penalty approximation of unilateral contact
{title2=$N=k\max(-d,0)$}

A <penalty approximation of unilateral contact> allows a small modeled overlap and uses a stiff restoring <force> instead of enforcing an exact rigid nonpenetration constraint. With signed gap $d$, positive in separation, a simple nonadhesive law is
$$
N=k\max(-d,0),\qquad k>0.
$$
An open or just-touching contact has zero <force> in this spring law; a sustained compressive load requires negative $d$. Thus it approximates <unilateral contact> rather than satisfying its rigid gap condition exactly. Tedrake's primary compliant-contact treatment uses stiff springs and damping to resist penetration and explains the resulting numerical <stiffness>.<Underactuated Robotics by Russ Tedrake>[Multi-Body Dynamics, Compliant Contact Models]{ref}

For an original gravity-supported example, a $2\,\mathrm{kg}$ <mass> with no other vertical load has $N=mg=19.62\,\mathrm N$. With $k=2\times10^5\,\mathrm{N\,m^{-1}}$, its static overlap is $mg/k=98.1\,\mathrm{\mu m}$. Increasing $k$ by a factor of 100 reduces that overlap to $0.981\,\mathrm{\mu m}$, but also makes the contact dynamics faster.

While the contact stays compressed, linearization about equilibrium gives
$$
m\ddot\xi+k\xi=0,\qquad\omega_c=\sqrt{k/m}.
$$
The <normal contact penalty time scale> is therefore proportional to $\sqrt{m/k}$. For the first <stiffness> above, $\omega_c\approx316.2\,\mathrm{rad\,s^{-1}}$; after increasing <stiffness> it is ten times larger. At fixed numerical step size, a smaller overlap error need not mean a more accurate trajectory.

For example, the <Forward Euler method> applied to this undamped oscillator has amplification factors $1\pm i\omega_c\Delta t$, whose magnitudes exceed one for every positive step. Merely shrinking the step does not make that method asymptotically stable for indefinite undamped oscillation, although it can reduce error on a finite interval. A suitable integration method, damping model and convergence check must be chosen. The source notes specifically identify stiff differential equations as a challenge for compliant contact simulation.<Underactuated Robotics by Russ Tedrake>[Multi-Body Dynamics, Compliant Contact Models, numerical integration discussion]{ref}

A penalty <stiffness> is a chosen model parameter, not automatically the actual material <stiffness>. Impact dissipation and tangential friction still need their own constitutive choices; a purely elastic normal spring does not imply a perfectly inelastic collision.