This document provides one way of thinking about divergence-free Eulerian fluids. It describes a different but similar technique to the 1999 stable fluids method and its more-tutorial presentation in 2003.
Assume we discretize 2D space using a square grid. We’ll store different properties of the fluid on different parts of the grid:
The resulting grid has three separate grids in the computer, of different sizes: an w×h grid of cell contents, an (w-1)×h grid of x fluxes, and an w×(h-1) grid of y fluxes. We’ll usually store these grids as flattened arrays to facilitate rapid access and linear algebra.
A 3×2 grid we’ll use in several examples:
This setup is called a staggered grid
or sometimes a MAC
grid
after the marker in cell
method that first popularized
it in 1965.
Divergence is a property measured in cells and is positive if more fluid is flowing out of the cell than into it, negative if more fluid is flowing in than out.
Incompressible fluids have divergence-free velocity fields: that is, all divergences are 0. Simple fluid simulation works in two steps: move the fluid ignoring the divergence-free constraint and then project the resulting velocities onto the nearest divergence-free flow.
The divergence of a 3×2 grid:
To remove divergence, we use the gradient of a pseudo-pressure field.
Real pressure has the property that the force it applies is the gradient of a scalar field and, like all gradients, is a purely divergence field. If we can subtract a purely-divergent flow from an existing flow we can make it divergence-free without changing any of the circulation that we want to keep.
True pressure involves terms like the speed of sound and is more complicated than we want to handle. Instead, we’ll make a pressure-like force that just removes divergence.
Given a pressure field
its gradient is a flux field created by subtracting adjacent cells of pressure
and the divergence of the pressure-induced flow is a sum of fluxes
which always equals n P_i - \sum_j P_j where j ranges over cell i’s 2–4 neighbors and n is the number of neighbors. We want to find the pressures that make that divergence equal to the negation of the divergences of the flow we are trying to make divergence-free \begin{align*} 2p_1-p_2-p_4 &= -d_1\\ 3p_2-p_1-p_3-p_5 &= -d_2\\ 2p_3-p_3-p_6 &= -d_3\\ 2p_4-p_1-p_5 &= -d_4\\ 3p_5-p_2-p_4-p_6 &= -d_5\\ 2p_6-p_3-p_5 &= -d_6\\ \end{align*} which we can pose as a linear system of equations \begin{bmatrix} 2&-1& & -1&&\\ -1&3&-1 & &-1&\\ &-1&2 & &&-1\\ -1&& & 2&-1&\\ &-1& & -1&3&-1\\ &&-1 & &-1&2\\ \end{bmatrix} \begin{bmatrix}p_1\\p_3\\p_3\\p_4\\p_5\\p_6\end{bmatrix} = \begin{bmatrix}d_1\\d_3\\d_3\\d_4\\d_5\\d_6\end{bmatrix}
The resulting system is known as the discrete Poisson equation and is particularly attractive because it is sparse (no row has more than 5 non-zero entries), symmetric (cell a_{ij} is the same as cell a_{ji}), weakly diagonally dominant (|a_{ii}| = \sum_{j≠i} |a_{ij}|), and positive definite (all eigenvalues are positive). While exact solutions can be computed in linear time, good-enough approximate solutions can be computed much more quickly with a few iterations of the conjugate gradient method
Real fluids preserve local kinetic energy, and hence vorticity or
circulation, except as diffused by viscosity. The Poisson-equation-based
removal of divergence does not do this: it makes velocities smaller to
remove diffusion, removing energy along the way. When energy is removed
because of how we discretize and solve a system, that energy loss is
called numerical damping.
Many approaches have been proposed for removing numerical damping, generally either by taking some kind of measure of energy or vorticity prior to performing advection and divergence removal and then re-inserting a similar level afterwards or by adding in some kind of particle-based vorticity tracker. There are also techniques that re-pose the entire fluid problem, for example in terms of circulation or turbulence.
All fluid properties (velocity, temperature, color, etc) move with
the fluid. This is called advection
; the movement of momentum
along the velocity field is called self-advection
.
All forms of advection need two concepts:
front buffer) and the new one we write to (sometimes called the
back buffer). Between time steps we swap the two.
We interpolate between grid locations and arbitrary points. Linear interpolation is generally used. Because we have a staggered grid, we need to interpolate on the correct one.
In a 2D grid, interpolating the velocity at point (2.3, 6.8) needs to consider both the x-velocity grid and the y-velocity grid.
In the x-velocity grid the surrounding four grid cells are (1.5,6), (2.5,6), (1.5,7), and (2.5,7) with weights 0.04, 0.16, 0.16, and 0.64 respectively.
In the y-velocity grid the surrounding four grid cells are (2,6.5), (2,7.5), (3,6.5), and (3,7.5) with weights 0.49, 0.21, 0.21, and 0.09 respectively.
Interpolation can be used to read a value at an arbitrary point, using the weights to average the values at several grid cells. It can also be used to add values at an arbitrary point, distributing weighted portions of the value to several grid cells.
Two forms of advection are common
Forward advection works as follows:
Forward advection conserves a value: the sum of the value across the grid is the same every time step. However, it can violate entropy by concentrating the value in a few cells, potentially causing instability if applied to velocity’s self-advection.
Backward advection works as follows:
Backward advection is unconditionally stable: the maximum magnitude of the value across the grid never increases as result of advection. However, it can violate the conservation of mass/momentum, potentially causing matter loss and numerical damping.
Since it was introduced in 1999, backward advection has been more popular than forward advection in computer graphics: we generally prefer stability over accuracy. However, both are used.