Physics-Based Simulation

Position-Based Fluids: Density and Surface Constraints

The Position-Based Dynamics framework can be elegantly extended from solids to simulate incompressible fluids. The resulting method, known as Position-Based Fluids (PBF) [Macklin & Muller 2013], replaces the complex pressure solves of traditional Smoothed Particle Hydrodynamics (SPH) [Koschier et al. 2020] with a set of geometric constraints. This approach inherits the stability and efficiency of PBD, allowing for large time steps suitable for real-time applications, while effectively enforcing the constant-density condition that characterizes incompressible flow.

The Per-Particle Density Constraint

The fundamental principle of PBF is to ensure that the density around each fluid particle remains constant. The density at a given particle's location is estimated using a kernel-based summation over its neighbors.

The density ${\rho }_i$ for a particle $i$ at position ${\mathbf{p}}_i$ is estimated as: $${\rho }_i=\sum _j{m}_{j}W({\mathbf{p}}_i-{\mathbf{p}}_j,h)\text{\hspace{1em}(33.4.1)}$$ where the sum is over all neighboring particles $j$, $m_j$ is the mass of particle $j$, $h$ is the smoothing kernel radius, and $W$ is a radially symmetric smoothing kernel function. For simplicity, we can assume all particles have equal mass and absorb it into the density calculation.

The goal is to enforce that this estimated density ${\rho }_i$ matches a user-defined rest density ${\rho }_0$. This is formulated as a per-particle constraint function $C_i$.

Definition 33.4.1 (PBF Density Constraint). For each particle $i$, the density constraint $C_i$ is defined as the scaled deviation from the rest density ${\rho }_0$: $$C_i({\mathbf{p}}_1,...,{\mathbf{p}}_N)=\frac{{\rho }_i}{{\rho }_0}-1\text{\hspace{1em}(33.4.2)}$$ The constraint is satisfied when $C_i=0$.

To solve this system of constraints, PBF computes a position correction $\Delta {\mathbf{p}}_i$ for each particle. Following the PBD methodology, we first compute the gradient of $C_i$ with respect to the position of a particle $k$: $${\nabla }_{{\mathbf{p}}_k}{C}_i=\frac{1}{{\rho }_0}\sum _j{m}_j{\nabla }_{{\mathbf{p}}_k}W({\mathbf{p}}_i-{\mathbf{p}}_j,h)=\{\begin{array}{ll}\frac{1}{{\rho }_0}{\sum }_j{m}_j\nabla W({\mathbf{p}}_i-{\mathbf{p}}_j,h)& \text{if }k=i\\ -\frac{m_k}{{\rho }_0}\nabla W({\mathbf{p}}_i-{\mathbf{p}}_k,h)& \text{if }k\ne i\end{array}\text{\hspace{1em}(33.4.3)}$$ The standard PBD solver computes a single Lagrange multiplier ${\lambda }_i$ for each constraint $C_i$ using the formula ${\lambda }_i=-C_i/({\sum }_k\Vert {\nabla }_{{\mathbf{p}}_k}{C}_i{\Vert }^2)$. The position correction for a particle $\Delta {\mathbf{p}}_i$ is then derived from the influence of its own constraint and the constraints of its neighbors. Due to the symmetry of the kernel gradient (${\nabla }_{{\mathbf{p}}_i}W({\mathbf{p}}_i-{\mathbf{p}}_j,h)=-{\nabla }_{{\mathbf{p}}_j}W({\mathbf{p}}_i-{\mathbf{p}}_j,h)$), this results in a simple and efficient final update rule for the position correction of particle $i$: $$\Delta {\mathbf{p}}_i=\frac{1}{{\rho }_0}\sum _j({\lambda }_i+{\lambda }_j)\nabla W({\mathbf{p}}_i-{\mathbf{p}}_j,h)\text{\hspace{1em}(33.4.4)}$$ Here, a different kernel, the Spiky kernel, is typically used for its non-vanishing gradient, which prevents particle clustering.

Remark 33.4.1 (Robustness and Constraint Softening). A practical issue arises when a particle has few neighbors, as the denominator ${\sum }_k\Vert {\nabla }_{{\mathbf{p}}_k}{C}_i{\Vert }^2$ can approach zero, leading to large, unstable position corrections. To prevent this, a small relaxation parameter $\epsilon$ is added to the denominator, softening the constraint. This is known as Constraint Force Mixing (CFM) [Smith & others 2005]. Then: $${\lambda }_i=-\frac{C_i({\mathbf{p}}_1,...,{\mathbf{p}}_N)}{{\sum }_k\Vert {\nabla }_{{\mathbf{p}}_k}{C}_i{\Vert }^2+\epsilon }\text{\hspace{1em}(33.4.5)}$$

Correcting for Tensile Instability

A well-known artifact in SPH-based fluid simulations is the "tensile instability," where particles near a free surface clump together due to an inability to satisfy the rest density. PBF offers a direct solution by incorporating an artificial pressure term.

Specifically, an artificial pressure term, $s_{\text{corr}}$, is introduced to create a short-range repulsive force between nearby particles. This term is a function of $W$: $$s_{\text{corr}}=-k{\left(\frac{W({\mathbf{p}}_i-{\mathbf{p}}_j,h)}{W(\Delta \mathbf{q},h)}\right)}^n\text{\hspace{1em}(33.4.6)}$$ where $k$ and $n$ are small positive constants (e.g., $k=0.1,n=4$), and $\Delta \mathbf{q}$ is a vector representing a small distance within the kernel radius (e.g., $\Vert \Delta \mathbf{q}\Vert =0.3h$). This term is only applied to push particles apart and is incorporated directly into the position correction calculation from Equation (33.4.4): $$\Delta {\mathbf{p}}_i=\frac{1}{{\rho }_0}\sum _j({\lambda }_i+{\lambda }_j+s_{\text{corr}})\nabla W({\mathbf{p}}_i-{\mathbf{p}}_j,h)\text{\hspace{1em}(33.4.7)}$$ This prevents clumping, and implicitly creates a surface tension-like effect at the fluid's boundary.

Velocity-Level Corrections: Viscosity and Vorticity

While the primary constraints in PBF operate on positions, velocity-level updates are still necessary to model phenomena like viscosity and to reintroduce lost energy. These are applied as a post-process after the main PBD solver loop has updated positions and velocities.

XSPH Viscosity: To impart a more coherent, liquid-like motion and reduce particle intermixing, XSPH viscosity is applied. It adjusts each particle's velocity to be closer to the average velocity of its neighbors: $${\mathbf{v}}_i^{\text{new}}={\mathbf{v}}_i+c\sum _j({\mathbf{v}}_j-{\mathbf{v}}_i)W({\mathbf{p}}_i-{\mathbf{p}}_j,h)\text{\hspace{1em}(33.4.8)}$$ where $c$ is a positive viscosity coefficient, typically a small value like $0.01$.

Vorticity Confinement: The iterative nature of PBD can introduce numerical damping, causing the fluid to lose rotational energy and appear overly viscous. Vorticity confinement can be used to counteract this by calculating the vorticity (curl of the velocity field) ${\mathbf{\omega }}_i=\nabla \times {\mathbf{v}}_i$ at each particle and applying a corrective force ${\mathbf{f}}_{\text{vorticity}}=\epsilon (\mathbf{N}\times {\mathbf{\omega }}_i)$ to reintroduce small-scale rotational details. This force is then integrated to update particle velocities.