Physics-Based Simulation

Collision Constraints

In any physical simulation, preventing the interpenetration of objects is paramount for achieving plausible results. Position-Based Dynamics provides a unified approach to this challenge. Collisions are not treated as a separate post-processing step involving impulses or penalty forces; instead, they are formulated as unilateral inequality constraints, just like IPC [Li et al. 2020] and integrated directly into the core PBD solver loop.

Triangle Self-Collisions

For thin-shell objects like cloth or other deformable surfaces, a primary challenge is handling self-collision, where the object folds and interacts with itself. The most common scenario is a vertex penetrating a triangle from another part of the mesh. This interaction is modeled with a unilateral constraint involving all four participating particles.

Consider a vertex with position $\mathbf{q}$ and a triangle with vertices ${\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3$. A collision constraint can be formulated by defining a minimum separation distance, or thickness $h$, from the triangle plane. The unilateral constraint function is: $$C(\mathbf{q},{\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3)=(\mathbf{q}-{\mathbf{p}}_1)\cdot \mathbf{n}-h\ge 0\text{\hspace{1em}(33.2.1)}$$ where $\mathbf{n}=({\mathbf{p}}_{2,1}\times {\mathbf{p}}_{3,1})/\Vert {\mathbf{p}}_{2,1}\times {\mathbf{p}}_{3,1}\Vert$ is the unit normal of the triangle. In case of $C<0$, all four particles ($\mathbf{q},{\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3$) are involved, and their positions are corrected according to their respective inverse masses to resolve the penetration while conserving momentum. It is also essential to check the barycentric coordinates of the vertex's projection onto the triangle plane to ensure the contact point lies within the triangle's boundaries before applying the correction.

Particle-Environment Collisions

The simplest collision scenario involves a dynamic particle interacting with a static, immovable piece of geometry, such as a floor plane or a convex container. This geometry acts as a boundary for the simulation domain. For a particle at position $\mathbf{p}$ interacting with a static plane, the non-penetration condition can be formulated as a constraint on the particle's signed distance from the plane.

Let the plane be defined such that ${\mathbf{n}}^T\mathbf{p}-d_{\text{rest}}=0$ for any point $\mathbf{p}$ on it, where $\mathbf{n}$ is the unit normal and $d_{\text{rest}}$ is an offset. A non-penetration constraint for a particle $\mathbf{p}$ is then written as an inequality: $$C(\mathbf{p})={\mathbf{n}}^T\mathbf{p}-d_{\text{rest}}\ge 0\text{\hspace{1em}(33.2.2)}$$ When a particle violates this constraint (i.e., $C(\mathbf{p})<0$), the PBD solver projects its position back to the surface. Since only one particle is dynamic, its inverse mass is effectively infinite compared to the static geometry, so it receives the full position correction required to satisfy $C(\mathbf{p})=0$. The correction simply moves the particle along the plane normal $\mathbf{n}$ to resolve the penetration. You can imagine how this can be simply extended to more complex shapes.

Particle-Particle Collisions

For simulating systems composed of discrete elements, such as granular materials, we have to consider direct particle-particle collision. This constraint prevents any two particles from overlapping.

For two particles at positions ${\mathbf{p}}_i$ and ${\mathbf{p}}_j$ with corresponding radii $r_i$ and $r_j$, the non-penetration constraint is a simple inequality based on their center-to-center distance: $$C({\mathbf{p}}_i,{\mathbf{p}}_j)=\Vert {\mathbf{p}}_i-{\mathbf{p}}_j\Vert -(r_i+r_j)\ge 0\text{\hspace{1em}(33.2.3)}$$ Unlike the linear plane constraint, this function is non-linear. This can be solved similarly to the stretching constraints considered in the previous section. The solver calculates a correction that pushes the two particles apart along the vector connecting their centers, distributing the correction based on their inverse masses to conserve linear momentum.

Frictional Effects at the Position Level

Friction is a dissipative contact force that opposes relative tangential motion between surfaces. To be more robust we can incorporate friction directly into the position-level constraint solve. This method is applied after an interpenetration constraint between two particles, $i$ and $j$, has been resolved.

Let ${\mathbf{p}}_i$ and ${\mathbf{p}}_j$ be the particle positions at the start of the time step. Let $\mathbf{p}{\ast }_i$ and $\mathbf{p}{\ast }_j$ be the current candidate positions, which have already been corrected to resolve penetration. The core idea is to compute a frictional position correction that opposes the tangential component of the particles' relative displacement during the timestep.

First, we determine the relative displacement vector over the timestep, $\Delta {\mathbf{p}}_{\text{rel}}$: $$\Delta {\mathbf{p}}_{\text{rel}}=(\mathbf{p}{\ast }_i-{\mathbf{p}}_i)-(\mathbf{p}{\ast }_j-{\mathbf{p}}_j)$$ Next, we find the tangential component of this displacement, $\Delta {\mathbf{p}}_t$, by projecting it onto the contact plane defined by the contact normal $\mathbf{n}$. $$\Delta {\mathbf{p}}_t=\Delta {\mathbf{p}}_{\text{rel}}-(\Delta {\mathbf{p}}_{\text{rel}}\cdot \mathbf{n})\mathbf{n}\text{\hspace{1em}(33.2.4)}$$ The friction model determines the magnitude of the correction based on a comparison between the tangential displacement $\Vert \Delta {\mathbf{p}}_t\Vert$ and the static friction threshold, which is proportional to the penetration depth $d$ and the coefficient of static friction ${\mu }_s$.

A position correction vector, $\Delta {\mathbf{p}}_{\text{friction}}$, is calculated to oppose $\Delta {\mathbf{p}}_t$. This correction is then distributed between the two particles. The correction for particle $i$ is given by: $$\Delta {\mathbf{p}}_i=-\frac{w_i}{w_i+w_j}\times \{\begin{array}{llc}\Delta {\mathbf{p}}_t& \text{if }\Vert \Delta {\mathbf{p}}_t\Vert \le {\mu }_{s}d& \text{(Static Friction)}\\ \Delta {\mathbf{p}}_t\min \left(\frac{{\mu }_{k}d}{\Vert \Delta {\mathbf{p}}_t\Vert },1\right)& \text{otherwise}& \text{(Kinetic Friction)}\end{array}\text{\hspace{1em}(33.2.5)}$$ where $w_i=1/m_i$ is the inverse mass and ${\mu }_k$ is the coefficient of kinetic friction. The negative sign is crucial, as friction must oppose the tangential displacement. In the static case, the correction fully cancels out the relative tangential movement. In the kinetic case, the correction is limited by the kinetic friction force (Coulomb's law).