SPH Liquid Simulation by Peter Kutz
Hi, here's some more information about my 3D SPH liquid simulator, which features weakly-incompressible SPH, pressure, viscosity, surface tension, gravity, boundaries using ghost particles, anisotropic kernels for surface reconstruction, surface reconstruction using marching cubes, OBJ export, and much more.
Images
Wet Paint
Wet Paint: A frame from a simulation involving 44,000 particles. Rendered in my own ray tracer.
Water Dyed Blue
Water Dyed Blue: A frame from a simulation involving 44,000 particles. Rendered in my own ray tracer.
Liquid Mercury
Liquid Mercury: High surface-tension droplets containing 9000, 1500, and 250 particles. Rendered in my own ray tracer. The droplets would look a little smoother if I had used the anisotropic kernels for the surface reconstruction.
As of now, I've only rendered stills in my ray tracer, but I plan on rendering animations as well. My ultimate goal would be to make animations that are close to indistinguishable from real videos. Before taking the time for rendering, I'll make sure to set up some nice simulations. I might also want to use more optimizations so rendering doesn't take too long. These images were made using Monte Carlo path tracing, which is slow at this quality. I could lower the quality, or use photon mapping (which I've already implemented) to make rendering faster.
Some Technical Details

Basics:
I first implemented almost everything in the original SPH fluids paper, Particle-Based Fluid Simulation for Interactive Applications (2003). I implemented the density and pressure computations, pressure force, viscosity force, gravity, a uniform grid with cells the size of the smoothing length for quick lookup of neighboring particles, and all of the kernel functions mentioned (poly6, spiky, and viscosity) as well as the gradients and Laplacians of some of those. I implemented symmetric force computations from another paper to comply with Newton's laws and cut the number of required computations in half.

Weakly Incompressible SPH:
Very low compressibility is an important property of water and other liquids. After looking through some more modern papers, I decided to implement weakly-compressible SPH from Weakly compressible SPH for free surface flows (2007). Weakly-compressible SPH uses a different computation for pressure (a Tait equation of state) and is more incompressible and realistic than the original SPH, but it also can require smaller time-steps because it effectively uses very stiff springs for the pressure. The paper also introduced a new and improved method for surface tension, which I also implemented.

Boundaries:
Boundaries are still an open problem in SPH fluid simulation. There are many available methods but none are perfect. A few popular methods are simply reflecting particle velocities at boundaries, using boundary particles, and using ghost-particles. One major problem with SPH is the erroneous low density at solid and free boundaries due to using an isotropic kernel to compute the density. This is a pain. In my simulations, this caused particles to clump at the boundaries to achieve the proper density, and that in turn caused the density of the inner particles to get too high and a gap formed between the outer layer of particles and the inner layers. This looked bad and caused unrealistic behavior. (Surface tension helps flatten the surface a little, but does not fix the floating layer or unrealistic behavior problems, and surface tension would round off the sharp corners and edges of the simulation in a box without ghost particles.) I implemented ghost-particles for the boundaries of the box in my simulation. The ghost particles are temporary copies of only the particles that are near the boundaries. They are almost identical to the original particles, except with the normal component of velocity reflected (in the free-slip condition). I also have the necessary ghost particles at corners and edges, not just flat surfaces (figured this out, not based on any paper in particular). I also reflect the velocity (using a coefficient of restitution) if the particles actually leave the box. I also added a cylindrical boundary with ghost particles.

Surface Reconstruction:
Surface reconstruction is another important part of realistic-looking fluid simulation using SPH. I started by implementing metaballs in my own ray tracer using ray marching. The results looked good but rendering was pretty slow. So I decided to implement marching cubes to polygonize the iso-surface. But just naively using metaballs results in blobby, bumpy looking fluid—even flat surfaces look bumpy. So I looked into better solutions and found a new paper called Reconstructing Surfaces of Particle-Based Fluids Using Anisotropic Kernels (2010) (which received the award for second-best paper at SCA 2010). This method creates elliptical kernels that do a much better job of smoothly and accurately representing the surface than existing techniques. It's important to note that the smoothing and anisotropic kernels are not propagated into the simulation—they are just used for rendering. I've implemented pretty much all of that paper including the diffusion smoothing step (which smooths out the fluid, and fixes the floating layer problem), the anisotropic kernels, the spatial data structure for storing the elliptical particles, and surface reconstruction using marching cubes. I found a few problems with the method as described in the paper that I had to work around, although it's possible that I just misinterpreted something.

Code:
I've implemented as much of the project as I can from scratch. For the OpenGL rendering, I'm using code that I wrote for another project that features a 3D scene graph, camera controls, and simple lighting. I used NIST's JAMA and TNT libraries for a singular value decomposition necessary to compute the anisotropic kernels. I'm currently using some marching cubes code from online, but have modified it to serve my purposes, in particular creating a continuous mesh instead of just a triangle soup, which is essential for shading and exporting to OBJ for rendering. I wrote the OBJ exporting myself, which was easy after I stitched the marching cubes results into a continuous triangle mesh. I used stb_image_write for outputting images, so I can play back animations.

Miscellaneous:
I have been initializing the particles in a uniform grid but giving them a random velocity at first to stir them, because uniform grids of particles can cause noticeable funny behavioral artifacts.

Poster
Poster Thumbnail

Click here to download the poster as a PNG (~2 MB).

Click here to download the full-resolution PDF (~8 MB).
Videos
I put together this short video showing a few animations from the simulator. I plan on rendering more, higher-quality videos soon.

Click here to see the video.
Resources
Here are some papers and other resources that I have used:

Particle-Based Fluid Simulation for Interactive Applications paper from 2003: Link

Weakly compressible SPH for free surface flows paper from 2007: Link

Reconstructing Surfaces of Particle-Based Fluids Using Anisotropic Kernels paper from 2010: Link

Particle-based Viscoelastic Fluid Simulation paper from 2005: Link

Basic marching cubes code: Link

Template Numerical Toolkit (TNT) (for singular value decomposition): Link

I particularly like the look of the water at 1:03 in this Lagoa Multiphysics video: Link

I've considered implementing predictive-corrective incompressible SPH: Link 1 Link 2
Other Work
Click here to see some of my other work, including some more information about my renderer, which I used to render the images of my liquid simulation.

© 2011 Peter Kutz