Write a Fluid Solver

This was my final project for the Advanced Animation Technology class. We were given the freedom to choose any difficult problem and try to solve it.

Learning Fluids

(Left): Houdini’s FLIP Fluid Solver. (Right): Houdini’s Solver with the pressure solve taken out. My job was to replace this part of the solver.

Before I undertook this project I had a hard time trying to understand new research in computer graphics. I simply didn’t understand the math and physics behind the important technologies that research is building on. As I tried to read papers from siggraph I found that a lot of things went over my head. In order to help change this I implemented the fluid pressure solve in Houdini’s FLIP fluid solver. I based it on a paper written at BYU for people like me who had never done something like this before.
(http://people.sc.fsu.edu/~jburkardt/pdf/fluid_flow_for_the_rest_of_us.pdf)

Research

Before I started writing the code, I needed to understand some important concepts. These are the three main areas of research that I focused on:

1. Learn Vector Calculus

·  Taking multi-variable calculus has been extremely helpful in understanding important concepts like the gradient and divergence operators.

·  One important thing to keep in mind is that we are using a “discretized grid.” So when we take a partial derivative we don’t actually compute the limit as our division size approaches zero--we just compute differences between voxels. We can do that in three different ways: forward differences ((n+1) - (n)), backward differences ((n) - (n-1)) and central differences ((n+1) - (n-1))

2. Learn the Navier-Stokes Equations

·  Ensure that a fluid is volume-preserving and that it obeys laws of physics such as the conservation of momentum.

·  Read “Fluid Flow for the Rest of Us” by David Cline and “Notes on Fluids” from Robert Bridson. (To see my notes on the papers see navier_stokes.png)

3. Learn how the FLIP solver in Houdini is organized

·  The HDK documentation sometimes isn’t very straightforward, but the information is there. The basic idea is to extend the GAS_SubSolver class and implement the pressure solve inside the inherited SolveGasSubclass() function. One of the parameters is a pointer to the FLIPFluidObject, and from there we can grab fields from the voxel grid and manipulate them.

·  In the pressure solve, we need to know which fluid cells are bordered by solids. This involves using collision field data, which isn’t as straightforward as accessing other fields on the voxel grid.

I spent most of my time on the first phase trying to get a basic understanding of these important ideas. I have a very valuable resource in my professor Seth Holladay, who is doing research in computational fluid dynamics.

Implementation

The basic fluid solve algorithm is fairly straightforward. However, some of the details can get complicated. Here I’ve outlined my understanding of each step as given in “Fluid Flow for the Rest of Us”.

·  Store data in the MAC (Marker and Cell) grid: pressure is stored at the center of the grid cell, while the components of velocity are stored on the faces.

·  Calculate a time step using the CFL condition. This basically just makes sure that a particle never moves more than one voxel in a single time step.

·  Update the data structure.

·  Advance the velocity field. Trace a particle backwards using the current velocity and grab the velocity at that cell. Apply external forces like gravity.

·  Calculate pressure: This is the heaviest part of the fluid solve. We want to calculate the pressure such that the divergence of velocity in the fluid is zero. To do so we must solve the equation Ap=b. A is a called the Laplacian matrix and when multiplied with our unknown pressure vector p, basically calculates a special second derivative. b is a vector that contains the divergence of the velocity and each component of b is found using this equation:

b_i = (density*cell_width) / (delta_t) * (divergence of velocity_i)

·  The HDK has a Linear Programming library that I can use to solve the linear equation Ap=b, but the hard part is setting up the sparse matrix A and handling all the boundary conditions correctly. We have to make sure we don’t try to grab a nonexistant value outside of our voxel grid.

·  Once we have our new pressure values, we can use them to correctly calculate velocities that satisfy the ‘divergence-free’ constraint of the Navier-Stokes equation:

u_new = u_old - (delta t)/(density *cell_height) * (gradient of pressure)

·  Now we move the particles through the grid.

·  And repeat forever . . .

After researching and understanding the ideas behind the FLIP solver, I feel much more confident in my ability to comprehend difficult mathematical and physical concepts. I plan to continue my research by adding to my fluid pressure solve and investigating new techniques in computational fluid dynamics.