Rigid Body Collisions

Collision handling for rigid bodies:

1. detect collisions
2. compute impulses to resolve collisions and to simulate repulsion
3. compute penalty forces to resolve penetrations due to numerical drift
4. time integration to get new rigid body positions, rotations and velocities

A rigid body has the following properties:

• a mass $m$,
• a position $\mathbf{x}$ (position of its center of mass),
• a velocity $\mathbf v$,
• an inertia tensor $I$ (note: in 3D this is not a scalar but a 3x3 tensor),
• a rotation angle $\alpha$ (note: in 3D we need a rotation matrix or a quaternion),
• an angular velocity $\omega$.

1. Detect collisions

In our example we use a very simple collision detection. Since all rigid bodies have a box shape, we just test all four vertices of the box if they are below the floor. This is done by checking if the y-coordinate is less than 0 (which defines the floor plane y=0).

Typically the four vertices of the box geometry are stored in local coordinates. We denote this vertices as $\mathbf b^\text{local}_i$. So before the collision test we have to transform the vertices to global coordinates. To transform the four vertices of the box of a rigid body to global coordinates, we need the angle $\alpha$ and the position $\mathbf x$ of the rigid body: $$\begin{align*} \mathbf R &= \begin{pmatrix} \cos(\alpha) & -\sin(\alpha) \\ \sin(\alpha) & \cos(\alpha) \end{pmatrix} \\ \mathbf b^\text{global} &= \mathbf R \mathbf b^\text{local} + \mathbf x. \end{align*} $$

2. Compute collision impulses

2.1. Relative velocity in contact points

First, we compute the relative point velocities for the two contact points of a collision: $$\mathbf u_\text{rel} = \mathbf u_1 - \mathbf u_2,$$ where the point velocity for a point in a rigid body is defined as $$\mathbf u = \mathbf v + \omega \times \mathbf r$$ and $\mathbf r = \mathbf c - \mathbf x$ is the world space vector from center of mass to the contact point $\mathbf c$.

Note that in our 2D case the cross product between the scalar $\omega$ and the vector $\mathbf r$ is defined as $\omega \times \mathbf r = [-\omega * \mathbf r_y,\ \omega * \mathbf r_x]^T$

2.2. Decomposition in normal and tangential velocity

In the next step we decompose the velocity vector in a normal part (in direction of the contact normal $\mathbf n$) and a tangential part: $$\begin{align*} \mathbf u^\text{n} &= (\mathbf u_\text{rel} \cdot \mathbf n) \mathbf n \\ \mathbf u^\text{t} &= \mathbf u_\text{rel} - \mathbf u^\text{n}. \end{align*}$$

2.3. Impulse computation

Now we can compute an impulse that should stop the motion in normal direction by using the impulse-based dynamic simulation approach (see also the example ibds.html). The velocity constraint is defined as: $$C(\mathbf u_1, \mathbf u_2) = (\mathbf u_1 - \mathbf u_2) \cdot \mathbf n = 0.$$ So the relative velocity in normal direction should be zero. For the impulse-based simulation we need the first derivatives with respect to $\mathbf v$ and $\omega$ and finally we get the required impulse as $$\mathbf p' = -\frac{1}{\frac 1 m + (\mathbf r \times \mathbf n)^2 I^{-1}} \mathbf u^\text{n}, $$ where the vector cross product in 2D is defined as $\mathbf r \times \mathbf n = \mathbf r_x \mathbf n_y - \mathbf r_y \mathbf n_x$.

Since we do not only want to stop the motion in normal direction but also want to simulate repulsion, we compute the final impulse using the restitution coefficient $e$ as $$\mathbf p = (1+e) \mathbf p'.$$

2.4. Apply impulse

Finally, we check if the resulting impulse points in direction of the contact normal ($\mathbf p \cdot \mathbf n > 0$) to avoid sticking. If this is the case, we apply the impulse and adapt the velocities of the corresponding rigid body as: $$\begin{align*} \mathbf v &:= \mathbf v + \frac{1}{m} \mathbf p \\ \omega &:= \omega + I^{-1} (\mathbf r \times \mathbf p). \end{align*}$$

2.5. Dynamic friction

A simple dynamic friction impulse can be determined according to Coulomb's law as $$\mathbf p^\text{friction} = -\mu \| \mathbf p \| \frac{\mathbf u^\text{t}}{\| \mathbf u^\text{t}\|}, $$ where $\mu$ is the friction coefficient.

3. Compute penalty forces

So far the system was solved on velocity level by computing impulses. Therefore, position errors (penetration) due to numerical errors cannot be considered. These errors sum up over the simulation and the bodies will penetrate the floor. This is also known as numerical drift.

To solve this problem we add a penalty force for each penetration that counteracts the violation on position level. The violation is determined by the penetration depth $d$ of the contact point. In our example this is simply the negative y-coordinate of the box vertex (since the floor was defined by y=0). Then the penalty force and a corresponding friction force are determined as $$\begin{align*} \mathbf f^\text{n} &= k d \mathbf n, \\ \mathbf f^\text{friction} &= -\mu \| \mathbf f^\text{n} \| \frac{\mathbf u^\text{t}}{\| \mathbf u^\text{t}\|}. \end{align*}$$ Finally, we also have to update the torque accordingly: $$\tau := \tau + (\mathbf r \times (\mathbf f^\text{n} + \mathbf f^\text{friction})). $$

4. Time integration