 Contact-mode timestep stability

Contact-mode timestep stability concerns whether a numerical integrator introduces growing errors when advancing the fast dynamics of an active penalty contact. It depends on the integration rule and the local modal frequencies, not just the size of the modeled overlap. The primary compliant-contact treatment identifies stiff contact springs as a source of numerically stiff differential equations.^[ref]

Freeze the active contacts and linearize about a compressed equilibrium. An undamped mode has $\ddot{x}+{\omega }^{2}x=0$. For step $h$, put $z=h\omega$. The Forward Euler method updates displacement and velocity with matrixwhose eigenvalues are $1\pm iz$. Their magnitudes are $\sqrt{1+z^2}>1$ whenever $z>0$. Thus the method has no nonzero stable step for indefinite integration of this undamped oscillator, although smaller steps reduce its error on a fixed finite interval.

The central-difference recurrence instead givesFor $0<z<2$ its two distinct amplification roots lie on the unit circle. At $z=2$ the repeated root is $-1$, and solutions can include $n(-1{)}^n$, so the strict bound matters. A collection of frozen undamped modes therefore requires $h{\omega }_{\max }<2$ for this oscillatory stability condition. In coupled normal-contact penalty modes, ${\omega }_{\max }$ comes from the coupled system rather than an arbitrary total-mass estimate.

Stability still allows substantial contact-mode numerical phase error. For central differences the phase per step is $\theta =2\arcsin (z/2)$, giving numerical angular frequency $\theta /h$. At $\omega =1000rad/s$ and $h=1\mathrm{ms}$, $z=1$ is stable but the numerical frequency is about $1047.2rad/s$, $4.72\%$ high. A timing-sensitive contact event can consequently drift even when oscillation amplitudes do not grow.

Implicit schemes illustrate another tradeoff. For the Backward Euler method, the amplification factors are $(1\mp iz{)}^{-1}$ and have magnitude $1/\sqrt{1+z^2}$. At $z=1$, each step reduces the oscillatory amplitude by about $29.3\%$, despite the continuous undamped system conserving energy. The trapezoidal rule instead has factorswith phase $2\arctan (z/2)$. The same example gives $927.3rad/s$, $7.27\%$ low. Unconditional stability for this linear model does not imply correct timing at a large step.

A practical error check therefore compares both amplitude and phase while reducing $h$. For instance, a small target dimensionless step $h{\omega }_{\max }=0.1$ gives much smaller phase error than one, but the acceptable value depends on the observable and time span; it is not a universal accuracy certificate.

The frozen-mode argument has a defined scope. Contact opening changes the active stiffness, nonlinear force laws change tangent stiffness, and damping introduces additional eigenvalue scales. An integrator can satisfy the local oscillator bound while missing an intervening collision or contact release. A proposed test starts with a small compressed-contact perturbation, checks frequencies and energy behavior, then repeats across a known contact-opening event with decreasing steps. Record force histories and event times as well as penetration, so numerical damping or phase drift is not mistaken for a material response.