Computing Contact Forces
Collaboration diagram for Computing Contact Forces:

Consider the definition of a contact with interpenetrating collision elements A and B.

The single, computed contact point serves as the origin of a contact frame C, so is designated Co; we'll shorten that to just C when it is clear we mean the point rather than the frame. We define the normal as pointing outward from the deformed surface of B towards the deformed interior of A. The C frame's z-axis is aligned along this normal (with arbitrary x- and y-axes). We are also interested in the points of bodies A and B that are coincident with Co, which we call Ac and Bc, respectively. Because the two forces are equal and opposite, we limit our discussion to the force f acting on A at Ac (such that -f acts on B at Bc).

Figure 1: Illustration of contact between two spheres.

The computation of the contact force is most naturally discussed in the contact frame C (shown in Figure 1).

The contact force, f, can be decomposed into two components: normal, fₙ, and tangential, fₜ such that f=fₙ+fₜ. The normal force lies in the direction of the contact frame's z-axis. The tangential component lies in the contact frame's x-y plane. In Drake's compliant contact model, although these components are orthogonal, they are not independent; the tangential force is a function of the normal force.

The model described here is adapted from Simbody's Hertz/Hunt & Crossley/Stribeck model described in [Sherman 2011]. We will summarize the elements of this model below.

Hunt-Crossley Normal Force

Drake uses the Hunt-Crossley model [Hunt 1975] for computing a normal force fₙ that accounts for both stiffness and dissipation effects. This is a continuous model based on Hertz elastic contact theory, which correctly reproduces the empirically observed velocity dependence of coefficient of restitution, where e=(1-dv) for (small) impact velocity v and a material property d with units of 1/velocity. In theory, at least, d can be measured right off the coefficient of restitution-vs.-impact velocity curves: it is the negated slope at low impact velocities.

Given a collision between two spheres, or a sphere and a plane, we can generate a contact force from this equation fₙ = kxᵐ(1 + mdẋ) where k is a stiffness constant incorporating material properties and geometry (to be defined below), x is penetration depth and is penetration rate (positive during penetration and negative during rebound). Exponent m depends on the surface geometry and captures the change in contact patch area with penetration. For Hertz contact where the geometry can be approximated by sphere (or sphere-plane) interactions, m=3/2. For contact geometry where the patch area is independent of penetration depth, m=1.

For Drake, we liberally extend the application of this model to point contact and (somewhat implausibly) take m=1. So, for Drake, the normal component of the contact force is currently:

fₙ = kx(1 + dẋ),

where k > 0 and d > 0. Please note, d is not a damping coefficient (which would have units of force/velocity), but is a dissipation factor, with units of 1/velocity, modeling power loss as a rate-dependent fraction of the conservative deformation-dependent force.

By definition fₙ should always be positive, so that the contact force is a repulsive force. Mathematically, for arbitrary x and ẋ, it is possible for fₙ to become negative, creating an attractive or "sucking" force. This case will be achieved if ẋ < -1 / d. To prevent sucking forces, the normal component is clamped to zero. In this regime, it is still possible for there to be a repulsive force for bodies that are drawing apart (ẋ < 0), as long as the relative velocities are small. This approximately models recovery of energy in the deformed material with observed hysteresis effects as described in [Hunt 1975]. However, a surface can't return to its undeformed state arbitrarily quickly. If the bodies are pulled apart faster than the surface can recover, the bodies will separate before the potential energy of deformation can be converted to kinetic energy, resulting in energy loss (to heat, vibration, or other unmodeled effects).

Stribeck Friction Tangential Force

Static friction (or stiction) arises due to surface characteristics at the microscopic level (e.g., mechanical interference of surface imperfections, electrostatic, and/or Van der Waals forces). Two objects in static contact need to have a force fₚ applied parallel to the surface of contact sufficient to break stiction. Once the objects are moving, dynamic (kinetic) friction takes over. It is possible to accelerate one body sliding across another body with a force that would have been too small to break stiction. In essence, stiction will create a contrary force canceling out any force too small to break stiction (see Figure 2).

Figure 2: Idealized Stiction/Sliding Friction Model

In idealized stiction, tangent force fₜ is equal and opposite to the pushing force fₚ up to the point where that force is sufficient to break stiction (the red dot in Figure 2). At that point, the tangent force immediately becomes a constant force based on the dynamic coefficient of friction. There are obvious discontinuities in this function which do not occur in reality, but this can be a useful approximation and can be implemented in this form using constraints that can be enabled and disabled discontinuously. However, here we are looking for a continuous model that can produce reasonable behavior using only unconstrained differential equations. With this model we can also capture the empirically-observed Stribeck effect where the friction force declines with increasing slip velocity until it reaches the dynamic friction level.

Figure 3: Stribeck function for stiction

The Stribeck model is a variation of Coulomb friction, where the frictional (aka tangential) force is proportional to the normal force as:

fₜ = μ⋅fₙ

In the Stribeck model, the coefficient of friction, μ, is replaced with a slip speed-dependent function:

fₜ = μ(s)⋅fₙ,

where s is a unitless multiple of a new parameter: slip tolerance (vₛ). Rather than modeling perfect stiction, it makes use of an allowable amount of relative motion to approximate stiction. When we refer to "relative motion", we refer specifically to the relative motion of the two points Ac and Bc on the corresponding bodies that are coincident in space with the contact point C.

The function, as illustrated in Figure 3, is a function of the unitless multiple of vₛ. The domain is divided into three intervals:

Other than the residual "creep" velocity limited by vₛ, which can be arbitrarily small (in theory; see next section for practical considerations), this model produces a reasonably good approximation of Coulomb friction. Its primary drawback is that the model is numerically very stiff in the stiction region, which requires either small step sizes with an explicit integrator, or use of a more-stable implicit integrator.

Next topic: Working with Contacts in Drake