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.

Contact-mode numerical phase error is the difference between the oscillation phase produced by a time integrator and the phase of the modeled contact dynamics. It can grow while the numerical solution remains bounded. Contact-mode timestep stability controls a different property; stable integration alone does not ensure accurate force timing or contact-event timing. The physical motivation is the fast oscillation introduced by a stiff compliant contact model.^[ref]

For a frozen undamped contact mode $\ddot{x}+{\omega }^{2}x=0$, central differences givePut $x_n=e^{in\theta }$ and use $e^{i\theta }+e^{-i\theta }=2\cos \theta$. On the stable branch $0<h\omega <2$,The numerical frequency exceeds the modeled frequency. A Taylor series givesFor the trapezoidal rule, its amplification-factor phase instead givesThat scheme preserves modal amplitude in this undamped linear setting while lagging in phase. These formulas isolate a single mode with fixed frequency and consistent initialization; general nonlinear contacts require additional analysis.

Accumulated contact-mode phase drift after a duration $T$ is $\Delta \varphi =T(\tilde{\omega }-\omega )$ at time-grid points. For an original example, let $\omega =1000rad/s$ and $T=0.10\mathrm{s}$. With $h=1\mathrm{ms}$, the central-difference frequency is about $1047.20rad/s$, producing $4.720\mathrm{rad}$ of phase lead. The step satisfies the strict stability inequality, yet the phase lead exceeds half a cycle.

Reducing $h$ to $50\mu \mathrm{s}$ gives about $0.01042\mathrm{rad}$ of lead over the same duration. At $25\mu \mathrm{s}$ it is about $0.002604\mathrm{rad}$, approximately one quarter as large, matching the leading second-order error. For a small desired phase budget ${\varphi }_{\max }$, the central-difference expansion suggestsWith ${\varphi }_{\max }=0.01\mathrm{rad}$, this estimate is about $49.0\mu \mathrm{s}$, much smaller than the $2\mathrm{ms}$ stability boundary. Check the exact frequency formula when the small-step expansion is not sufficiently accurate.

In a contact experiment or simulation, a shifted oscillation phase changes when force maxima and gap crossings are predicted. Damping, changing tangent stiffness and a contact opening can invalidate the fixed-frequency model, so the phase estimate is not by itself an event-time error certificate. A proposed convergence check compares force histories, dominant frequency and event times at successively halved steps, while keeping contact geometry and constitutive parameters fixed. Agreement in maximum penetration alone can miss a phase error that accumulates over many cycles.

Accumulated contact-mode phase drift is the phase difference that builds over many timesteps when a numerical contact mode oscillates at a slightly incorrect frequency. A small error per step can matter after many cycles even when the amplitude remains bounded. The fast oscillations arise naturally in stiff compliant-contact models, whose numerical integration difficulty is described in the primary multibody notes.^[ref]

For a fixed undamped mode of angular frequency $\omega$, step $h$, and numerical phase advance $\theta (h\omega )$ per step, after $N$ steps the unwrapped phase error isUnwrapped phase counts complete cycles; reducing it modulo $2\pi$ would hide a missed or added cycle. This expression assumes the same fixed mode and a uniform step throughout the comparison.

For central differences, $\theta (z)=2\arcsin (z/2)$ on $0<z<2$. Expanding givesFor the trapezoidal rule, $\theta (z)=2\arctan (z/2)$ instead givesBoth are second-order phase errors at fixed duration, but their signs and leading coefficients differ. Halving $h$ approximately quarters the accumulated phase error; doubling the simulated duration approximately doubles it when the mode remains unchanged.

For an original example, choose $\omega =600rad/s$, $h=0.5\mathrm{ms}$ and $T=0.20\mathrm{s}$, so $N=400$ and $h\omega =0.3$. Central differences give $\tilde{\omega }\approx 602.273rad/s$ and a phase lead about $0.4546\mathrm{rad}$, or $26.05^{\circ }$. The relative frequency error is only about $0.379\%$, yet the phase lead has become readily visible. The leading-order estimate is $0.4500\mathrm{rad}$.

An approximate phase budget $|\Delta {\varphi }_N|\le {\varphi }_{\max }$ for central differences suggestsWith the preceding mode and duration and ${\varphi }_{\max }=0.020\mathrm{rad}$, it gives $h\lesssim 105.4\mu \mathrm{s}$. Taking $100\mu \mathrm{s}$ gives an exact phase lead about $0.01801\mathrm{rad}$. This is far below the $3.333\mathrm{ms}$ oscillator stability boundary; contact-mode timestep stability and timing accuracy impose different requirements.

The associated timing shift is approximately $-\Delta \varphi /\omega$ for a fixed phase feature of a nearly sinusoidal mode with small drift. It is not a general contact-event error bound: a changing amplitude, a grazing zero, changing tangent stiffness or an active-contact transition can invalidate that approximation. A variable-step method also needs its actual update analysis; summing uniform-step formulas without checking how successive maps compose is not automatically exact.

A proposed convergence test records a small compressed-contact oscillation over a stated duration, compares unwrapped numerical and reference phases, and repeats with half the step and twice the duration. Compare the force waveform and event times as well as the maximum penetration. Material damping changes the physical response, while numerical phase drift changes the timing of that modeled response; they should be tested separately.