Fluid Flow TutorialKarl Sims


Flow field A flow field is represented by a grid of velocity vectors. The example on the right shows rotational vortex motions about two centers.
A flow field can be visualized by advecting tracer particles or textures through it. Particles can be "pushed" through a flow field by reading the flow velocities to move their positions, but images are "pulled" by warping them in the local direction of flow. For the examples below, animation cycles were made by repeatedly inserting lines of particles or stamping a checker pattern in the tracer image over time.


Particles are pushed: each particle's position is updated by the flow velocity at the particle's current position, interpolating between flow vectors as appropriate. 
Textures are pulled: a tracer image is repeatedly warped by reading color values from upstream pixels in the negative direction of flow. The image is typically at higher resolution than the flow field.
 
Fluid momentum
The visualizations above show tracers moving through a static flow field, but the flow itself can also dynamically change over time. To simulate fluid momentum, the flow field can simply flow itself. Flow vectors are updated by reading interpolated values from upstream grid cells in the same way the tracer image is advected above. New vector values are pulled from the negative direction of flow. Momentum on its own causes circular flows to expand away from their centers which would create lower pressure areas analogous to those occurring at the center of hurricanes or tornadoes. To properly simulate incompressible fluids, we next need to even out the pressure differences caused by diverging and converging flow. 
A flow field warps itself by reading from upstream values. 
Divergence
Divergence is a measure of the local expansion in a vector field, and is calculated at a given point by summing the outward components of neighboring vectors. Negative divergence corresponds to contraction. Translation, rotation, and shearing motions have zero divergence because their inward and outward components are balanced.
Positive divergence  Negative divergence  Zero divergence  
Removing divergence
To enforce incompressibility and remove changes in pressure, we attempt to zero out the divergence throughout the flow field. This is done by adding flow away from high pressure converging areas, and toward low pressure diverging areas. A gradient is the slope or rate of change of a value across a grid. The gradient of the divergence is a vector field that points from high pressure to low pressure areas. So to remove the divergence from a flow field, we repeatedly increment it by the gradient of its divergence.
We can combine the divergence and gradient calculations, and if we find divergence values at the corners between grid cells instead of at their centers, we can use a simple 3x3 convolution for this. The left image below represents four local divergences in different colors. Each divergence value is the sum of the outward diagonal vector components of the 4 adjacent cells. The horizontal (x) gradient at the center cell is found by the two right divergence values minus the two left (green + blue  red  orange). The vertical (y) gradient is found by the two top minus the two bottom (red + green  orange  blue).
Divergence at 4 corners  Divergence x gradient  Divergence y gradient  
The divergence x and y gradients are shown above as 3x3 grids of vectors, and the neighboring flow vectors in those directions (dot products) can be summed for that gradient component of the center cell. For stability, the gradient is scaled by 1/8 before adding it to the flow field. Combining these x and y gradients using vector operations, the kernel code to update each flowfield vector might look something like this:
f'[x,y] = f[x,y] + (dot(f[x1,y1] + f[x+1,y+1], {1,1}) + dot(f[x1,y+1] + f[x+1,y1], {1,1}) * {1,1} + (f[x1,y] + f[x+1,y]  f[x,y1]  f[x,y+1]) * {2,2} + f[x,y] * 4) * 1/8
where f[x,y]
is the vector at location x,y
of the flow field, f'
is the updated flow field,
and {i,j}
are constant vectors. This step is performed many
times to spread the incompressibility across the flow field, perhaps
50 iterations per frame or more depending on the resolution and fluid
behavior.
The images below show an initial vector field in red with nonzero divergence, and the results in blue after divergence removal. On the left, the outward momentum forces from a rotation are cancelled giving more stable swirls. On the right, the effect of an external pulse at the center is spread throughout the grid and starts to generate a pair of vortices like a canoe paddle in water.
Removing divergence from the red vector field prevents swirls from spreading outward.  Removing divergence from an external pulse propagates it into a vortex pair. 
Boundary conditions
To find divergence values at the cells on the edge of the grid, we need to generate flow vectors from just outside the grid. Different methods for this can give different edge behaviors:
Obstacles can also be added to the grid with similar boundary conditions by forcing the flow, or its normal component, to zero at their locations.
External forces
Realistic looking fluid behavior can be generated by alternating between the fluid momentum and divergenceremoval steps described above. However some nonzero flow velocity needs to be set somehow, either procedurally or interactively. For example, adding linear flow at specific locations can create a squirting ink effect (see below), or tracking mouse or camera motion can be used to interactively push the fluid.
Other forces can be added to simulate various physical phenomena:
Summary of procedure
For each fluid simulation time step:

Multigrid efficiency
The divergence removal process can be sped up significantly using multiple grid resolutions. We can average the flow field down to halfres and remove the divergence from that, then add the difference between the processed and original halfres grid back to the fullres grid, and finally perform a few divergenceremoval iterations at fullres. This can be done recursively so most of the work is done at 1/4 or even lower resolution, which easily allows for realtime results when implemented on a modern GPU.
Accuracy
This method is meant to create useful visual results rather than accurately model physics, but here are a few tips that might reduce unwanted artifacts. For fast moving fluid, use multiple fluid simulation time steps per frame. Use higher resolution for tracer images than for the flowfield grid. Use smooth bicubic interpolation if possible when reading from tracer images and flow fields. For repeated warping of tracer images that aren't being modified in other ways, temporarily warp their coordinates instead, and then resample after multiple iterations.
3D
These methods are extendable to 3D, using volume data for flow fields and tracer images, and 3x3x3 convolutions for divergence removal, but of course it may be more difficult to generate 3D results in realtime, depending on your hardware and volume dimensions.
References
Jos
Stam, "Stable Fluids", SIGGRAPH 1999 Conference Proceedings
NavierStokes equations
Visit Flow, an interactive exhibit at MIT
that employs these techniques.
Back to other work by Karl Sims
© 2018, Karl Sims, All rights reserved.