Fluid Simulation Tutorial

This tutorial walks you through the advection and projection steps of simulating incompressible fluids on a grid. The algorithm follows the approach described in "SIGGRAPH 2007 Course notes on fluids" by Robert Bridson and Matthias Muller-Fischer.

This tutorial 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.

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!