SPH Fluid Simulation

Recalled from the previous post, there are 2 major approaches used to describe the fluid motion, Eulerian and Lagrangian. Eulerian is grid-based, measuring properties at fixed points in space. Physical values are stored on grids. Lagrangian is particle-based, measuring properties as particles drift through the flow. Physical values are stored on each particle. Both methods have pros and cons. The prevalent methods used nowadays are mostly the hybrids. (e.g. FLIP is the hybrid between particle-based and volume-based).

Let's start from the particle-based approach and take SPH (Smoothed Particle Hydrodynamics) for example. Via discrete particles located arbitrary in space , SPH utilizes weighted summation from adjacent particles to interpolate values and derivative of continuous field.

SPH uses smoothing kernel (W) to determine which adjacent particles have how much contribution of values on one particular point in space. The selected smoothing kernel should be normalized, positive, even and last but not least, no kernel interaction outside the range of the radius "h". The derivatives of the field could be approximated with analytical differentiation:

Since the physical value at point j doesn't directly depend on other components in space, thus, the gradient of the field only has effect on the smoothing kernel. Same for Laplacian.

With above SPH equations, we can start to approximate each term in Navier-Stokes's Equation. Since the main idea of Lagrangian method is the particle itself moves with the fluid, thus the advection term can be ignored, which simplify the equation a lot. For other terms (e.g. density, pressure, viscosity, external forces, etc), they all can be approximated with SPH approach. I won't paste the equation for each term here since they all can be found from lots of material online, and the main purpose of this article is just for getting myself familiar with the big picture and techniques behind the scene. However, there is a thing worth mentioning,  if just using SPH equation above to directly approximate something like pressure and viscosity force, we will confront an asymmertric problem. Think about a situation, if there are only 2 particles with contact in space, the pressure of one particle only uses the pressure force of another, and since the pressure for each partilce may be different, the result pressure will not be symmetrical. These can be fixed with few ways to rewrite the equation.

OK. With the above knowledge, let's talk about the advantages and disadvantages of SPH. There are several benefits of using particle-based approach over pure grid-based method. First of all, the equation is more intuitive and easier to implement. Second, since SPH stores mass on particles, it guarantees the conservation of mass w/o extra calculation. Third, SPH doesn't have to track fluid boundaries; it supports free surface effect. But there still has to be a way for mesh generation for rendering since SPH doesn't contain topology info.

The obvious drawback of using particle methods is the need of large numbers of particles to get more accurate result, which enhances computational cost. The other thing is the density summation would cause problems when neighboring particles have different rest densities (i.e. different masses). This would lead to incorrect pressure and have weird gaps (break into pieces?). Some solutions have been introduced which involves alternative density computation, formulas rewriting and new surface tension model.

What if we need the boundary to limit the movement of particles? One way could be adding boundary particles to prevent penetration and inverse the velocity. Another way is to use Ghost Particles. There is a technical paper about Ghost SPH in the upcoming Siggraph 2012: http://www.cs.ubc.ca/~rbridson/docs/schechter-siggraph2012-ghostsph.pdf

Navier-Stokes Equation

Here is a summarized note I took after reading some material relevant to fluid effect.

This is the compact version of Navier-Stokes Equation excerpted from Jos Stam's paper. For simulation in Computer Graphics, physical accuracy isn't the primary but the appearance and performance is. In order to get a faster simulation and realistic at the same time, here it assumes constant density, which means the fluid is incompressible. Basically, fluid effect in CG is mostly based on this assumption.

The first equation is the condition of incompressibility. With constant density, the net flow in and out of a region should sum to zero, which means the divergence of velocity should be zero, no sinks or sources.

The second equation is based on Newton's second law F=ma. Before going further, I spent some time trying to understand 2 different descriptions for fluid mechanics, Eulerian and Lagrangian respectively. I recommend the following video on YouTube, which does a great job on explaining these 2 ideas. It also provides a clear demonstration about how to deduce material derivative, which is quite important to help understand the second equation here.

http://www.youtube.com/watch?v=mdN8OOkx2ko&feature=related

It's obvious to see that the left hand side of the assignment, combining with the first term on the right hand side is the material derivative. It says that the change of a point moving within a velocity field is the sum of the change at the fixed position and the one due to the movement.  The first term on the right hand side is what we call "Advection" (imaging if we drip a pigment drop into a moving liquid, how the pigment will move within it).

OK. Back to F=ma. 

Here we represent the mass as density since unit size is used; then expand the material derivative to get the following:

Well... divide each side by density, we get equation 2 :)

It's apparent that the sum of (F)orces corresponds to the rest 3 terms in equation 2:  pressure, viscosity and other external forces.

Pressure is used for enforcing incompressibility. Since the gradient points to the greatest rate of increase and fluid should flow from high pressure to low, negative sign is used for reversing.
Viscosity is used for simulating viscous fluid. Laplace operator is used here for getting smoothing result. External forces are something like gravity, buoyancy, etc.

So lots of research are around solving this equation, leading to various algorithms, such as SPH (Smoothed Particle Hydrodynamics), PIC/FLIP (Particles in Cell / Fluid Implicit Particle). It's kind of complicated and better has separate articles for these.