Rigid Body Collisions


Collision handling for rigid bodies:This example shows a simple method handling collisions between rigid bodies and a floor:
A rigid body has the following properties:
1. Detect collisionsIn 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 ycoordinate 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 impulses2.1. Relative velocity in contact pointsFirst, 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 velocityIn 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 computationNow we can compute an impulse that should stop the motion in normal direction by using the impulsebased 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 impulsebased 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 impulseFinally, 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 frictionA 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 forcesSo 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 ycoordinate 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 integrationThe rigid bodies are advected by numerical time integration. In our case we use a symplectic Euler method: $$\begin{align*} \mathbf v(t + \Delta t) &= \mathbf v(t) + \frac{\Delta t}{m} \mathbf f^{\text{ext}} \\ \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t) \\ \omega(t + \Delta t) &= \omega(t) + \Delta t I^{1} \tau \\ \alpha(t + \Delta t) &= \alpha(t) + \Delta t \omega(t + \Delta t), \end{align*}$$ where $\mathbf f^{\text{ext}}$ are the external forces, $\omega$ is the angular velocity and $\alpha$ is the rotation angle. 