Physics-Based Simulation

Preliminaries

To fully appreciate the PBD approach, it is essential to first understand the classical formulation of Lagrangian dynamics it gets inspired. The state of a dynamic system composed of $N$ particles is described by their individual mass $m_i$, position ${\mathbf{p}}_i\in {\mathbb{R}}^3$, and velocity ${\mathbf{v}}_i\in {\mathbb{R}}^3$. The evolution of this system is governed by a pair of first-order ordinary differential equations derived from Newton's second law: $${\dot{\mathbf{v}}}_i=\frac{1}{m_i}{\mathbf{f}}_i\text{\hspace{1em}(32.1.1)}$$ $${\dot{\mathbf{p}}}_i={\mathbf{v}}_i\text{\hspace{1em}(32.1.2)}$$ where ${\mathbf{f}}_i$ is the sum of all forces on particle $i$.

Rigid bodies require additional attributes to describe their rotational state: an orientation quaternion ${\mathbf{q}}_i\in \mathbb{H}$, an angular velocity ${\mathbf{\omega }}_i\in {\mathbb{R}}^3$, and an inertia tensor ${\mathbf{I}}_i\in {\mathbb{R}}^{3\times 3}$. The rotational motion is then described by the Newton-Euler equations: $${\dot{\mathbf{\omega }}}_i={\mathbf{I}}_i^{-1}\left({\mathbf{\tau }}_i-\left({\mathbf{\omega }}_i\times ({\mathbf{I}}_i{\mathbf{\omega }}_i)\right)\right)\text{\hspace{1em}(32.1.3)}$$ $${\dot{\mathbf{q}}}_i=\frac{1}{2}{\hat{\mathbf{\omega }}}_i{\mathbf{q}}_i\text{\hspace{1em}(32.1.4)}$$ where ${\mathbf{\tau }}_i$ is the net torque and ${\hat{\mathbf{\omega }}}_i$ is the quaternion $[0,{\mathbf{\omega }}_i]$.

To simulate this evolution, these continuous equations are discretized using a numerical integrator. The Symplectic Euler method updates the velocity first, then uses this new velocity to update the position, improving stability over standard explicit Euler. For a particle, with time step $\Delta t$, the update is: $${\mathbf{v}}_i(t_0+\Delta t)={\mathbf{v}}_i(t_0)+\Delta t\frac{1}{m_i}{\mathbf{f}}_i(t_0)\text{\hspace{1em}(32.1.5)}$$ $${\mathbf{p}}_i(t_0+\Delta t)={\mathbf{x}}_i(t_0)+\Delta t{\mathbf{v}}_i(t_0+\Delta t)\text{\hspace{1em}(32.1.6)}$$ This procedure is applied analogously for rigid body states.

Remark 32.1.1 (Quaternion Normalization). Due to numerical integration error, the quaternion ${\mathbf{q}}_i$ may drift from unit length. It is essential to re-normalize the quaternion after each integration step to maintain a valid rotational state.

Finally, interactions and physical limits are modeled using holonomic constraints, which depend only on positions and orientations, but not velocities! Constraints are kinematic restrictions in the form of equations and inequalities that constrain the relative motion of bodies. An equality (bilateral) constraint takes the form $C_j(\mathbf{p},\mathbf{q})=0$, while an inequality (unilateral) constraint is $C_j(\mathbf{p},\mathbf{q})\ge 0$. In classical dynamics, these are satisfied by computing constraint forces and adding them to ${\mathbf{f}}_i$ in Equation (32.1.1). It is this specific mechanism for handling constraints that PBD fundamentally changes.