Physics-Based Simulation

Rigid Body Dynamics

Before jumping into ABD, let's review the traditional approach to simulate rigid body dynamics.

Let’s begin by reviewing Newtonian mechanics — specifically, Newton’s second law for a particle at position $\mathbf{x}(t)$ with mass $m$:

$$\begin{array}{rl}\frac{d\mathbf{x}}{dt}& =\mathbf{v},\\ m\frac{d\mathbf{v}}{dt}& =\mathbf{f}.\end{array}$$

Now, compute $\frac{d(m\mathbf{x}\times \mathbf{v})}{dt}$:

$$\begin{array}{rl}\frac{d(m\mathbf{x}\times \mathbf{v})}{dt}& =m\frac{d\mathbf{x}}{dt}\times \mathbf{v}+m\mathbf{x}\times \frac{d\mathbf{v}}{dt}\\ & =m\mathbf{v}\times \mathbf{v}+\mathbf{x}\times \mathbf{f}\\ & =0+\mathbf{x}\times \mathbf{f}\\ & =\mathbf{x}\times \mathbf{f}.\end{array}$$

Here, $m\mathbf{x}\times \mathbf{v}$ is called the angular momentum, and $\mathbf{x}\times \mathbf{f}$ is the torque.

For a rigid body $\mathcal{B}$, we can integrate over its volume:

$$\frac{d}{dt}\left({\int }_{\mathcal{B}}\rho \mathbf{x}\times \mathbf{v}d\mathbf{x}\right)={\int }_{\mathcal{B}}\rho \mathbf{x}\times \mathbf{f}(\mathbf{x})d\mathbf{x},$$ where $\rho =\rho (\mathbf{x})$ is the mass density of $\mathcal{B}$ at position $\mathbf{x}$.

Now let $\mathbf{c}$ denote the center of mass (COM) of $\mathcal{B}$. Then:

$$\begin{array}{rl}\frac{d}{dt}\left({\int }_{\mathcal{B}}\rho \mathbf{c}\times \mathbf{v}d\mathbf{x}\right)& ={\int }_{\mathcal{B}}\rho \mathbf{c}\times \frac{d\mathbf{v}}{dt}d\mathbf{x}\\ & ={\int }_{\mathcal{B}}\rho \mathbf{c}\times \mathbf{f}(\mathbf{x})d\mathbf{x}.\end{array}$$

Combining both:

$$\frac{d}{dt}\left({\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})\times \mathbf{v}d\mathbf{x}\right)={\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})\times \mathbf{f}(\mathbf{x})d\mathbf{x}.$$

Since ${\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})\times {\mathbf{v}}_{\mathbf{c}}d\mathbf{x}=({\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})d\mathbf{x})\times {\mathbf{v}}_{\mathbf{c}}=0\times {\mathbf{v}}_{\mathbf{c}}=0$, we can further simplify:

$$\frac{d}{dt}\left({\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})\times (\mathbf{v}-{\mathbf{v}}_{\mathbf{c}})d\mathbf{x}\right)={\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})\times \mathbf{f}(\mathbf{x})d\mathbf{x},$$ where the right-hand side is the torque $\mathbf{\tau }$ about the center of mass.

Denoting the right-hand side as the torque ${\mathbf{\tau }}_{\mathbf{c}}$, compute the left-hand side:

$$\begin{array}{rl}\frac{d}{dt}({\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})\times (\mathbf{v}-{\mathbf{v}}_{\mathbf{C}})d\mathbf{x})& =\frac{d}{dt}({\int }_{\mathcal{B}}\rho (\mathbf{x}-\mathbf{c})\times (\mathbf{\omega }\times (\mathbf{x}-\mathbf{c}))d\mathbf{x})\\ & =\frac{d}{dt}({\int }_{\mathcal{B}}\rho (⟨\mathbf{x}-\mathbf{c},\mathbf{x}-\mathbf{c}⟩\mathbf{\omega }-⟨\mathbf{x}-\mathbf{c},\mathbf{\omega }⟩(\mathbf{x}-\mathbf{c}))d\mathbf{x})\\ & ={\mathbf{I}}_{\mathcal{B}}^{\mathbf{c}}\mathbf{\omega }.\end{array}$$

Here, the inertia tensor ${\mathbf{I}}_{\mathcal{B}}^{\mathbf{c}}$ is defined as:

$${\mathbf{I}}_{\mathcal{B}}^{\mathbf{c}}:={\int }_{\mathcal{B}}\rho \left(\Vert \mathbf{x}-\mathbf{c}{\Vert }_2^2\mathbb{I}-(\mathbf{x}-\mathbf{c})(\mathbf{x}-\mathbf{c}{)}^T\right)d\mathbf{x}.$$

Since $\mathbf{x}$ depends on $t$, so does ${\mathbf{I}}_{\mathcal{B}}^{\mathbf{c}}$. Letting $\mathbf{R}(t)$ be the rotation matrix (the rotation of $\mathcal{B}$ from time $0$ to $t$). It's easy to verify that:

$${\mathbf{I}}_{\mathcal{B}}^{\mathbf{c}}(t)=\mathbf{R}(t){\mathbf{I}}_{\mathcal{B}}^{\mathbf{c}}(0)\mathbf{R}(t{)}^T.$$

Differentiating yields:

$${\mathbf{\tau }}_{\mathbf{c}}=\mathbf{\omega }\times ({\mathbf{I}}^{\mathbf{c}}\mathbf{\omega })+{\mathbf{I}}^{\mathbf{c}}\frac{d\mathbf{\omega }}{dt}.$$

Thus, the rigid body dynamics law becomes:

$$\begin{array}{rl}\frac{d}{dt}(m{\mathbf{v}}_{\mathbf{c}})& =\mathbf{f},\\ \mathbf{\omega }\times ({\mathbf{I}}^{\mathbf{c}}\mathbf{\omega })+{\mathbf{I}}^{\mathbf{c}}\frac{d\mathbf{\omega }}{dt}& ={\mathbf{\tau }}_{\mathbf{c}}.\end{array}\text{\hspace{1em}(25.1.1)}$$

These equations resemble Newton's second law. The first describes how forces affect the linear momentum of the rigid body, while the second describes how torques influence its angular momentum.

With a simple explicit integrator, we can simulate rigid body motion via:

$$\begin{array}{rl}{\mathbf{v}}_{n+1}& ={\mathbf{v}}_n+\Delta t\frac{{\mathbf{f}}_n}{m},\\ {\mathbf{x}}_{n+1}& ={\mathbf{x}}_n+\Delta t{\mathbf{v}}_n,\\ {\mathbf{\omega }}_{n+1}& ={\mathbf{\omega }}_n+\Delta t({\mathbf{I}}^{\mathbf{c}}{)}^{-1}({\mathbf{\tau }}_{\mathbf{c}}-\mathbf{\omega }\times ({\mathbf{I}}^{\mathbf{c}}\mathbf{\omega })),\\ {\hat{\mathbf{q}}}_{n+1}& ={\mathbf{q}}_n+\frac{\Delta t}{2}\mathbf{\omega }{\mathbf{q}}_n,\\ {\mathbf{q}}_{n+1}& =\frac{{\hat{\mathbf{q}}}_{n+1}}{|{\hat{\mathbf{q}}}_{n+1}|}.\end{array}\text{\hspace{1em}(25.1.2)}$$

Here, $\mathbf{q}$ is the unit quaternion representing the body’s rotation. Its update is based on a first-order approximation:

$$\begin{array}{rl}\tilde{\mathbf{q}}-\mathbf{q}& =e^{\frac{\Delta t}{2}\mathbf{\omega }}\mathbf{q}-\mathbf{q}\\ & =(\cos \frac{\Vert \mathbf{\omega }\Vert \Delta t}{2}-1+\frac{\mathbf{\omega }}{\Vert \mathbf{\omega }\Vert }\sin \frac{\Vert \mathbf{\omega }\Vert \Delta t}{2})\mathbf{q}\\ & =(\Theta (\Delta t^2)+\frac{\Delta t}{2}\mathbf{\omega })\mathbf{q},\end{array}$$ which leads to: $$\underset{\Delta t\to 0}{\lim }\frac{\tilde{\mathbf{q}}-\mathbf{q}}{\Delta t}=\frac{1}{2}\mathbf{w}\mathbf{q},$$

where $\mathbf{w}=({\mathbf{\omega }}_x,{\mathbf{\omega }}_y,{\mathbf{\omega }}_z,0)$; $\tilde{\mathbf{q}}$ is the composition of the rotation represented by the rotation vector $\Delta t\mathbf{\omega }$ and the rotation represented by $\mathbf{q}$.

Remark 25.1.1 (Rotation Vector to Quaternion). A rotation of angle $\theta$ around the axis defined by the unit vector $\mathbf{u}=({\mathbf{u}}_x,{\mathbf{u}}_y,{\mathbf{u}}_z)={\mathbf{u}}_x\mathbf{i}+{\mathbf{u}}_y\mathbf{j}+{\mathbf{u}}_z\mathbf{k}$ can be represented by the quaternion $$\mathbf{q}=e^{\frac{\theta }{2}({\mathbf{u}}_x\mathbf{i}+{\mathbf{u}}_y\mathbf{j}+{\mathbf{u}}_z\mathbf{k})}=\cos \frac{\theta }{2}+({\mathbf{u}}_x\mathbf{i}+{\mathbf{u}}_y\mathbf{j}+{\mathbf{u}}_z\mathbf{k})\sin \frac{\theta }{2}=\cos \frac{\theta }{2}+\mathbf{u}\sin \frac{\theta }{2}$$