Position Based Fluids
As discussed in Position Based Simulation, position-based dynamics provided the foundations of simulating lagrangian fluids. Remember PBD is all about how to put constraints on particles and how to resolve them. Similarly, we can model how fluid particles interact by posing some constraints on them.
Density Constraint
The density constraint1 aims to model the compressibility on fluids by constraining the density (measuring the amount of mass per unit volume) of a single fluid particle.
Since we only care about local density around a given particle, the density constraint takes all its neighboring particle positions \(\mathbf{x}_1, \cdots, \mathbf{x}_\mathbf{n_j}\), and sum up their masses to \(\rho_i\). Noted that each particles are weighted by a smoothing function \(W\) that falls off radially from the center of the kernel, which we will discuss more in the smoothing kernel section. The ratio to its rest density \(\rho_o\) is then evaluated, indicating whether the fluid is being compressed locally (e.g. \(\rho_i > \rho_o\)). In a perfect equilibrium, \(\rho_i\) should equal to \(\rho_o\) and the constraint \(C_i=0\) will be satisfied, meaning this fluid particle is happy :)
Scaling up the above to a \(N\)-particle simulation, every particles will be assigned a density constraint, to ensure its local density remain the same at all time. This is why the bilateral condition \(C(\mathbf{x})=0\) is used here. Once we run the solver for a decent number of iterations, every fluid particles should be happily resting at their stable positions.
What if the density is lower than the rest condition (e.g. \(\rho_i < \rho_o\))?
There will be cases where a lone particle wandered into the void, where a few or even no neighboring particles are present in its surroundings. This posed a tensile instability in our simulation, which can cause clustering or clumping due to negative pressures.
This can result in extreme correction vectors and can break the simulation. One solution is to clamp the density constraint result to zero to prevent negative pressure.
However, this also means the cohesion nature of fluid particles is gone, particles almost always repel with each other. So by clamping the constraint, we are also limiting ourselves from having this cohesive behavior.
Luckily, in the paper "SPH without a tensile instability"2, an artificial pressure correction factor is proposed. Where \(\Delta\vec{q}\) is vector to a point at some fixed distance inside the smoothing kernel radius \(h\) and \(k\) is a small positive constant. \(\|\Delta\vec{q}\|=\{0.1h, \cdots, 0.3h\}\), \(k = 0.001\) and \(n = 4\) work well in my case.
Smoothing Kernels
Smoothing kernels are essential for integrating local density of a particle in an SPH simulation. A good kernel needs to fulfill a few requirements:
- Radially symmetric (i.e. only depends on distance)
- Finite support (i.e. particles should only affect up to a certain distance)
- Differentiable and integrates to 1
| Kernel Type | Function \(W(\vec{r}, h)\) |
|---|---|
| Poly6 (2D) | \(W_{poly6}(\vec{r}, h) = \frac{4}{\pi h^8} \begin{cases} (h^2-\|r\|^2)^3 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
| Poly6 (3D) | \(W_{poly6}(\vec{r}, h) = \frac{315}{64 \pi h^9} \begin{cases} (h^2-\|r\|^2)^3 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
| Spiky (2D) | \(W_{spiky}(\vec{r}, h) = \frac{10}{\pi h^5} \begin{cases} (h-\|r\|)^3 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
| Spiky (3D) | \(W_{spiky}(\vec{r}, h) = \frac{15}{\pi h^6} \begin{cases} (h-\|r\|)^3 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
In order to attract or repel neighboring particles, we need to compute the gradient of the smoothing kernel \(\nabla W(\vec{r}, h)\).
| Kernel Type | Gradient \(\nabla W(\vec{r}, h)\) |
|---|---|
| Poly6 (2D) | \(\nabla W_{poly6}(\vec{r}, h) = -\vec{r}\frac{24}{\pi h^8} \begin{cases} (h^2-\|r\|^2)^2 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
| Poly6 (3D) | \(\nabla W_{poly6}(\vec{r}, h) = -\vec{r}\frac{945}{32 \pi h^9} \begin{cases} (h^2-\|r\|^2)^2 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
| Spiky (2D) | \(\nabla W_{spiky}(\vec{r}, h) = -\frac{\vec{r}}{\|\vec{r}\|}\frac{30}{\pi h^5} \begin{cases} (h-\|r\|)^2 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
| Spiky (3D) | \(\nabla W_{spiky}(\vec{r}, h) = -\frac{\vec{r}}{\|\vec{r}\|}\frac{45}{\pi h^5} \begin{cases} (h-\|r\|)^2 & \text{if } 0\leq \|r\|\leq h\\ 0 & \text{else} \end{cases}\) |
Bloopers
| 1 | 2 | 3 |
|---|---|---|
 |
 |
 |