Fluid Simulation Tutorial

Do you want to simulate fluids, but don't know how to get started?!? Below, is the tutorial I wrote while co-instructor of our Physically based animation course at UPenn. It walks you through using a staggered MAC Grid as well as your first advection step and projection steps. I provide computed values which you can then compare against your own code.

This tutorial assumes that you've read the "SIGGRAPH 2007 Course notes on fluids" by Robert Bridson and Matthias Muller-Fischer. We will use their notation as well as refer to concepts we they already describe with excellent detail and clarity in their notes.

This tutorial is meant to be used with the fluid framework basecode that I wrote for the corresponding class assignment on fluids. This code is the basis for my own simulations of smoke and fire as well as those of many of our awesome animation students who were able to extend it in various ways to make their own unique and excellent demos. The basecode includes an orbit camera, MAC Grid implementation, conjugate gradient solve for the pressure step, rendering of velocities, rendering of smoke density using sheets (similar to the one described in "Visual Simulation of Smoke" by Fedkiw, Stam, and Wann Jensen), recording options, and configurable grid size. If you like this writeup and use it make something cool, please tell me about it. It will make me happy. If you don't like it, well, it is free afterall. Seriously, though, if you have comments or something is unclear, please tell me about that too.

And one last thing before we get started:

The most important milestone for fluid simulation is creating a **stable, divergence free velocity field**. Once
the velocity field is stable, it is simple to advect arbitrary quantities within the grid, such as smoke
density, heat, color, or whatever. By divergence free, we mean

- the volume of each cell remains constant each time step
- the velocity entering each cell should equal the velocity leaving each cell
- if we sum the rates of change in velocity in every direction, they should sum to zero, e.g. \[ \nabla \cdot u = \frac{du}{dt} + \frac{dv}{dt} + \frac{dw}{dt} = 0 \]

Our first step is to understand the MAC Grid structure. The MAC Grid is a data structure for discretizing the space of our fluid. Quantities in the world, such as velocities and smoke density, will be stored at cell faces and centers. We use a staggered MAC Grid in particular because it allows us to estimate the derivatives (e.g. the rates of change in each direction) with more numerical stability (read the Bridson/Muller-Fischer notes to understand why!).

To start with something simple, create a 2x2x1 MAC grid as follows.

X faces are used to store the u-component of the velocity. | Y faces are used to store the v-component of the velocity. | Z faces are used to store the w-component of the velocity. |

In our small MAC Grid above, we have 4 cells, 8 Z faces, 6 Y faces, and 6 X faces. We will store velocities at the faces
and other quantities, such as smoke density, at the cell centers.
Specifically, if we define the velocities with (u,v,w), we store the u values at X-faces, v values are Y-faces, and
w values at Z-faces. The only tricky thing to remember is that *we will need to compute the full velocity (u,v,w) at each face location, but then only store one of the components there*. Implementation-wise, we have separate arrays, mU, mV, and mW for storing each velocity component. (i,j,k) indexing is used for each array of faces, for example,

- The front Z face indices are (1,1,0)(0,1,0)(1,0,0)(0,0,0). The back Z face indices are (1,1,1)(0,1,1)(1,0,1)(0,0,1).
- The Y face indices are (0,0,0),(0,1,0),(0,2,0) and (1,0,0)(1,1,0)(1,2,0).
- The X face indices are (0,0,0),(0,1,0) and (1,0,0),(1,1,0) and (2,0,0)(2,1,0).

An empty MAC Grid, where all values are zero, corresponds to a grid with no velocity in it. Our first step is to add some movement to shake things up. To simplest way to start is to add a non-zero velocity to one of the internal faces of the grid. Let's add 1 unit of velocity in the Y direction, e.g. mV(0,1,0) = 1.

- What is the velocity at the center of each cell? Either (0,0,0) or (0,0.5,0)

FOR_EACH_FACE currentpt = position at center of current face currentvel = velocity at center of current face oldpt = currentpt - dt * currentvel newvel = velocity at old location store one component of newvel depending on face typeComputing the velocities at arbitrary places in the grid is done by interpolating between the values stored at each face.

- What is the position of the Y-face? (0.5,1.0,0.5)
- What is oldpt? (0.5,1.0,0.5) - (0.1)*(0,1,0) = (0.5, 0.9, 0.5)
- What is the velocity at position (0.5,0.9,0.5)? (0,0.9,0)
- What is the new value of mV(0,1,0)? 0.9

- Because of traceback and interpolation, we never get a velocity that is larger than the values stored on the grid. This results in dissapation and loss of detail, but ensures that the system stays stable. Detail and dissapation can be counteracted by adding additional forces such as vorticity confinement.
- Parameter tweaking tip: make sure your timestep and velocities are not too large. In particular, you never want your traceback step to be larger than your grid size.
- What should you do when your traceback step takes you outside of the grid? If the position corresponds to a fluid source, just return the source value. You can also extrapolate the velocity from the nearest boundary, or set the value to zero (although this results in greater dissappation near boundaries)

We will compute pressure values at the center of each cell which ensure that the resulting velocity field is divergence free. |

- The coefficient for \(p_{i,j,k}\) is the number of fluid neighbors
- The coefficient for fluid neighbors is -1
- The coefficient for non neighbors is 0

- Check that the velocity field after this step is divergence free by computing the divergence at each cell.
- If you use a uniform grid size, \(\Delta x = \Delta y = \Delta z\)
- What is the velocity normal to a solid border? \(u_{solid}\)
- The fluid has free movement tangent to the boundaries.
- The resulting velocity field should be look smooth, no sudden changes!

Once we have a divergence free velocity field, we can advect other properties in the field. For
smoke, the fluid we simulate is the air, within which we would advect the smoke density
(**carefull!** the smoke density is *different* from
the air fluid density which we use in the projection step!) or temperature for buoyancy forces.
Putting it all together, a typical simulation step might look like

void SmokeSim::step() { double dt = 0.01; // Step0: Gather user forces mGrid.updateSources(); // Step1: Calculate new velocities mGrid.advectVelocity(dt); mGrid.addExternalForces(dt); mGrid.project(dt); // Step2: Calculate new temperature mGrid.advectTemperature(dt); // Step3: Calculate new density mGrid.advectDensity(dt); }To simulate fire, we can simulate two fluids simultaneously: one for the fuel sources and one for the air. The tricky parts are keeping track of the flame front (e.g. the boundary between the two fluids) and ensuring that mass is preserved with the two fluids coupled together (see "Physically Based Modeling and Animation of Fire"). To simulate explosions, we can modify the projection step so that instead of ensuring the divergence is zero, we solve for a non-zero divergence which estimates the explosion's blast wave (without explicitly simulating it, phew, see "Animating Suspended Particle Explosions"). To simulate water, we need to keep track of the water surface, say with a level set, and handle the interactions between the fluid and air (see the SIGGRAPH notes for a good introduction). And finally, best of luck and happy coding!