Drake
The Details of Computing Contact Forces

Drake includes a compliant contact model.

Compliant models determine contact forces by assuming that all bodies are, to some extent, pliable. As two bodies collide, they deform and that deformation produces the contact normal force. In practice, this is achieved by allowing undeformed geometry to penetrate and then infer the deformation which would eliminate the penetration and compute the force from that inferred deformation. More particularly, Drake's compliant model is adapted from Simbody's Hertz/Hunt & Crossley/Stribeck model described in [Sherman 2011] with further details drawn directly from [Hunt 1975].

The model naturally decomposes into the normal and tangential components of the contact force and the discussion follows this decomposition: normal force tangent force.

## Computing the Normal Component of the Contact Force

### Overview of Compliant Normal Force

Generally, the component of the contact force due to pressure (i.e., not frictional forces) is simply the pressure on the contact surface integrated over that surface's area. If we assume the contact surface is planar then we have a well defined normal direction, and the force due to pressure is parallel with that normal and can be defined precisely as:

fₙ = p(q, v)·A(q)

where A(q) is the contact patch area, and p(q, v) is the average contact pressure on that patch, the details of which depend on how the contact is characterized (e.g., penetration depth and approaching speed), and the geometry and material parameters of the contacting bodies. Notice there is no approximation in this equation; it is always true since this equation is the definition of the average contact pressure p(q, v). Thus, determining the contact normal force consists of determining the geometry of the contact patch, and the average pressure on that patch.

### Hunt-Crossley Model

Drake uses the Hunt-Crossley model [Hunt 1975] for computing a normal force fₙ that accounts for both elasticity and dissipation effects. This is a continuous model based on Hertz elastic contact theory. We'll examine the underlying Hertz contact model and then show its extension.

Hertz Contact Model

The Hertz contact model applies to scenarios where the contacting surfaces can be locally represented by two principal axes of curvature (e.g., spheres, ellipsoids, cylinders, planes, etc.) For the sake of simplicity, this discussion focuses on spheres where a single radius of curvature (per surface) is sufficient to describe the contacting geometry. The primary advantage of this model is that the magnitude of the normal force can be expressed in closed form.

Consider two spheres made of the same material and of the same size (with Young's modulus E and radius R). If the amount that the two spheres are penetrating is x, then the resultant normal force (according to the Hertz contact model) would be:

fₙ = ⁴/₃⋅E⋅√R⋅√x³.

The two spheres compress such that they are touching along a disk with radius √(R⋅x). Based on this, we can equate it to the earlier function:

• fₙ = p(q, v)·A(q) = ⁴/₃⋅E⋅√R⋅√x³
• p(q, v) = 4/(3π)·E·√(x/R)
• A(q) = π⋅R⋅x

This can be generalized to spheres of different sizes and different materials by creating effective radii of curvature and Young's modulus. Assuming body i has radius Rᵢ and Young's modulus Eᵢ, we can define an effective radius of curvature and Young's modules as functions of the contacting geometries: R(R₁, R₂) and E(E₁, E₂), respectively. Generally, the details of those functions depend on the geometry and scenario. See the implementation details below for more details.

This pattern can be extended to contact between other shapes with characteristic radii of curvature:

• Sphere and plane: same as sphere and sphere with the plane treated as a sphere with infinite radius of curvature.
• Two perpendicularly crossed cylinders of equal radius R: same as sphere and sphere of radius R.
• Vertical cylinder of radius R and a plane: fₙ = 2⋅E⋅R⋅x.

Hunt-Crossley Extension

The Hunt-Crossley model extends the Hertz contact by adding a dissipation factor which depends on the rate of change of penetration. The resultant force is:

fₙ = H(x)(1 + ³/₂⋅d⋅ẋ),

where H(x) is the Hertz factor appropriate for the contact geometry, x is the penetration depth, ẋ is penetration rate (positive for increasing penetration and negative during rebound), and d is a dissipation term that captures the empirically observed velocity dependence of the coefficient of restitution, e = (1 - d⋅v), for (small) impact velocity v. d has units of 1/velocity and, in theory at least, can be measured right off the coefficient of restitution-vs.-impact velocity curves; it is the negated slope at low impact velocities.

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. For steel, bronze or ivory, [Hunt 1975] reports values of d between 0.08-0.32 sec/m.

By definition fₙ should always be positive, so that the contact force is a repulsive force. Mathematically, for arbitrary x and ẋ, it is possible for fₙ to become negative, creating an attractive or "sucking" force. This case will be achieved if ẋ < -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 (ẋ < 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:

• s ∈ [0, 1): the coefficient of friction rises smoothly from zero to the static coefficient of friction, μs.
• s ∈ [1, 3): The coefficient of friction falls smoothly from μs to the dynamic (kinetic) coefficient of friction, μd.
• s ∈ [3, ∞): Coefficient of friction is held constant at μd.

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: Drake Contact Implementation