 Coupled normal-contact penalty modes

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$, soThe contact-Jacobian force/velocity pairing and stiff compliant-contact model come from the primary multibody formulation; the local modal reduction here follows by substitution.^[ref]

Although $AK$ need not be symmetric, it is similar to the symmetric positive-semidefinite matrix $K^{1/2}A{K}^{1/2}$. Set $z=K^{1/2}\xi$ to obtainIts 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$. ThenWith equal contact stiffness $k$, the equal-gap mode is vertical translation and the opposite-gap mode is rocking. Their angular frequencies areTake $m=2\mathrm{kg}$, $I=0.005\mathrm{kg}{\mathrm{m}}^2$, $a=0.10\mathrm{m}$ and $k=10^{5}N/m$. Then $A$ has diagonal entries $2.5\mathrm{k}{\mathrm{g}}^{-1}$ and off-diagonal entries $-1.5\mathrm{k}{\mathrm{g}}^{-1}$. Its eigenvalues are $1$ and $4\mathrm{k}{\mathrm{g}}^{-1}$, giving $316.2$ and $632.5rad/s$. Treating the contacts independently would assign each $500rad/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 equationThe 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\approx 3.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.