Effective mass in a normal-contact penalty mode ($m_{\mathrm{eff}}=(J{M}^{-1}{J}^{\mathsf{T}}{)}^{-1}$)

The effective mass in a normal-contact penalty mode is the inertia seen by one specified relative normal displacement, rather than automatically the total mass of either body. A contact away from a center of mass can excite rotation as well as translation, changing that inertia and the normal contact penalty time scale.

Use independent generalized velocities $v$, a positive-definite mass matrix $M$, and a normal gap $g$ whose local velocity is $\dot{g}=Jv$. A scalar normal force $\lambda$ produces generalized force $J^{\mathsf{T}}\lambda$. The primary multibody formulation gives this contact-Jacobian velocity/force pairing explicitly.^[ref] At a frozen configuration, the part of gap acceleration caused by that force isIf the scalar $a=J{M}^{-1}{J}^{\mathsf{T}}$ is positive, defineThe definition uses a physical normal gap measured in metres and a matching normal force measured in newtons. Arbitrarily rescaling the gap coordinate changes the associated inertia. Constraints must first be accounted for by independent coordinates or the appropriate constrained response.

A short normal impulse $P$ gives $\Delta v=M^{-1}{J}^{\mathsf{T}}P$, hence $\Delta \dot{g}=P/m_{\mathrm{eff}}$. This supplies a direct interpretation: less effective mass means a larger relative normal-velocity change for the same impulse. At a fixed local configuration, a compressed linear penalty spring of stiffness $k$ gives the small-motion estimateIt omits changing contact geometry, other driving forces and damping; the source separately explains why stiff compliant contact challenges numerical integration.^[ref]

For two freely translating particles, the relative gap sees $a=1/m_1+1/m_2$, so $m_{\mathrm{eff}}=m_1{m}_2/(m_1+m_2)$. This is the reduced mass. Equal $2\mathrm{kg}$ particles give $m_{\mathrm{eff}}=1\mathrm{kg}$; replacing one by a fixed wall gives $2\mathrm{kg}$ instead.

For a planar rigid body above a fixed plane, let its mass be $m$, rotational moment of inertia about its center be $I$, and the contact normal be vertical. A contact a horizontal distance $r$ from the center hasFor $m=2\mathrm{kg}$, $I=0.020\mathrm{kg}{\mathrm{m}}^2$ and $r=0.10\mathrm{m}$, $m_{\mathrm{eff}}=1\mathrm{kg}$. At the center it is $2\mathrm{kg}$. With $k=10^{5}N/m$, the off-center angular frequency is $316.2rad/s$, compared with $223.6rad/s$ at the center. Using the full mass in both cases misses this rotational response.

Several simultaneous contacts require more than separate scalar replacements. Their matrix $A=J{M}^{-1}{J}^{\mathsf{T}}$ generally has off-diagonal terms: a force at one contact changes another contact’s acceleration. With frozen geometry and positive diagonal stiffness matrix $K$, compressed gap perturbations satisfy $\ddot{\xi }+AK\xi =0$; its nonzero modal squared frequencies are those of the symmetric matrix $K^{1/2}A{K}^{1/2}$. Coupled normal-contact penalty modes can therefore be faster than estimates treating each contact independently. A proposed verification applies a small known impulse or displacement perturbation at a stated configuration and checks relative contact motion, followed by a timestep convergence study.

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.

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.