= Coupled normal-contact penalty modes
{title2=$\ddot\xi+(JM^{-1}J^{\mathsf T})K\xi=0$}
<Coupled normal-contact penalty modes> are small oscillations generated by several simultaneously compressed penalty contacts acting through shared body inertia. A displacement at one contact can accelerate another because both load the same translating or rotating bodies. Separate scalar <effective mass in a normal-contact penalty mode> estimates do not generally give the coupled <frequencies>.
Freeze a contact configuration, use independent generalized displacement perturbations $q$, positive-definite <mass matrix> $M$, and gap perturbations $\xi=Jq$. Let $K$ be a positive diagonal matrix of normal penalty stiffnesses. After subtracting equilibrium loads, the restoring generalized <force> is $-J^{\mathsf T}KJq$, so
$$
M\ddot q+J^{\mathsf T}KJq=0,\qquad
\ddot\xi+AK\xi=0,\quad A=JM^{-1}J^{\mathsf T}.
$$
The contact-Jacobian <force>/velocity pairing and stiff compliant-contact model come from the primary multibody formulation; the local modal reduction here follows by substitution.<Underactuated Robotics by Russ Tedrake>[Multi-Body Dynamics, The Dynamics of Contact, contact Jacobian and generalized force equation; Compliant Contact Models]{ref}
Although $AK$ need not be symmetric, it is similar to the symmetric positive-semidefinite matrix $K^{1/2}AK^{1/2}$. Set $z=K^{1/2}\xi$ to obtain
$$
\ddot z+(K^{1/2}AK^{1/2})z=0.
$$
Its positive eigenvalues are squared angular <frequencies>. A zero eigenvalue may describe redundant gap coordinates outside the realizable displacement subspace; it should not automatically be interpreted as an extra physical oscillation. Free generalized motions satisfying $Jq=0$ are a separate possibility.
For an original planar example, a rigid bar has vertical translation $y$, small rotation $\theta$, <mass> $m$, center-of-<mass> moment of inertia $I$, and two vertical contacts at horizontal offsets $-a,+a$. Then
$$
J=\begin{pmatrix}1&-a\\1&a\end{pmatrix},\qquad
M=\begin{pmatrix}m&0\\0&I\end{pmatrix},\qquad
A=\begin{pmatrix}1/m+a^2/I&1/m-a^2/I\\1/m-a^2/I&1/m+a^2/I\end{pmatrix}.
$$
With equal contact <stiffness> $k$, the equal-gap mode is vertical translation and the opposite-gap mode is rocking. Their angular <frequencies> are
$$
\omega_{\rm translate}=\sqrt{2k/m},\qquad
\omega_{\rm rock}=\sqrt{2ka^2/I}.
$$
Take $m=2\,\mathrm{kg}$, $I=0.005\,\mathrm{kg\,m^2}$, $a=0.10\,\mathrm m$ and $k=10^5\,\mathrm{N/m}$. Then $A$ has diagonal entries $2.5\,\mathrm{kg^{-1}}$ and off-diagonal entries $-1.5\,\mathrm{kg^{-1}}$. Its eigenvalues are $1$ and $4\,\mathrm{kg^{-1}}$, giving $316.2$ and $632.5\,\mathrm{rad/s}$. Treating the contacts independently would assign each $500\,\mathrm{rad/s}$ and miss the faster rocking mode.
<Contact-mode timestep stability> depends on the largest active <frequency> and on the integrator. For the undamped central-difference recurrence, one mode has amplification equation
$$
r^2-2(1-h^2\omega^2/2)r+1=0.
$$
The distinct oscillatory roots lie on the unit circle for $0<h\omega<2$; the endpoints need separate treatment. The example therefore requires $h<2/632.5\approx3.162\,\mathrm{ms}$ for this linear stability condition. Resolving the waveform accurately normally requires a smaller step. This is not a stability rule for every integration method.
The derivation applies while both contacts remain compressed and geometry is locally constant. A sufficiently large rocking amplitude can unload one contact, changing the active model. A proposed check preloads both supports, applies small symmetric and antisymmetric perturbations, and compares measured or simulated <frequencies> with the two predictions. Repeat with smaller amplitudes and smaller timesteps to separate linearization error, contact opening and numerical error.
Back to article page
