SPH

Particle Based Fluid Simulation

1. Introduction
This project is implemented particle-based approach using SPH technique to simulate fluids with free surfaces. We derive the viscosity and pressure force fields from Navier stokes equation and added a term to model surface tension effects. For interactivity purpose, we used smoothing kernels proposed in paper “Particle Based Fluid Simulation for interactive application” by Matthias Müller, David Charypar and Markus Gross. Then we used these particles to render the surface of the fluid using the marching cube based surface reconstruction method.
We are able to simulate nearly 4000 particles at interactive rates. However, fluid surface construction takes a very long time.
2. Algorithm
2.1. Smooth Particle Hydrodynamics
SPH is an interpolation method for particle systems where each field quantity is calculated at discrete particle location. These particles have a spatial distance known as the "smoothing radius", represented by h in equations, over which SPH distributes the quantities using kernel function. Hence, we can get the value of any physical by summing the relevant properties of all the particles which lie within the range of the kernel.
So, we use a function as mentioned below to find a value of a scalar quantityA:

where,
j iterates over all the particles
, is the mass of the particle,
, is position of particle,
, is density,
, is the field quantity at j

The function W (r, h) is called the smoothing kernel with core radius h.
2.2. Simulation
2.2.1. Equations
We use Navier-Stokes equation to describe the motion of fluid in our simulation

where,
g is an external force density field and µ the viscosity of the fluid.

2.2.2. Derivations
2.2.2.1. Pressure Term
Pressure is calculated using ideal gas state equation. We used the following modified equation:

where,
is the rest density to calculate the pressure.
We calculate density before calculating pressure. In doing so, we have kept mass as constant throughout our simulation. However, we calculate density at each time step since density is varying. The formula to calculate density is:

2.2.2.2. Viscosity
We calculate the viscosity following the same symmetrization behavior for stability by using equation:

However, in this case we symmetrize the viscosity using velocity differences.
2.2.3. External forces:
Our application supports external forces such as gravity and collision forces. These forces are directly applied to the particles at the time of time integration. To add support for various boundaries, we follow below steps:
• We check the intersection of particle with the boundaries of the container at the time of calculating velocity derivative.
• If the particle is colliding with the container surface, we calculate the normal of the surface.
• Then , we simply reflect the velocity component of the particle by multyplying it with the normal.
In this way, we have implemented 3 types of boundaries:
• Rectangular
• Circular
• Spherical-Bowl shape
These boundaries were implemented by using particle-cube intersection, particle-circle intersection, particle-sphere intersection.We additionally implemented arbitrary boundaries using particle-polygon intersection. But the method we described above is not handling the arbitrary boundaries efficiently. Hence, we have disabled that functionality from our application and plan to implement it as a part of future implementation.
2.2.4. Smoothing kernel:
As we know, stability, accuracy and speed of the SPH method highly depend on the choice of the smoothing kernels. We have used the following spiky kernel for pressure computations:

We have taken viscosity kernel dependent on velocity field. Hence, Laplacian kernel is used for calculating viscosity term,

We have also defined the smoothing kernels to calculate the surface tension. The smoothing kernels which we have used have the following formula:

m_gradient_Poly6Kern = 945 / (32 * PI * pow( h, 9) )

is first applied to the gradient of color field. Then, we applied

m_Lap_Poly6Kern = 945 / (8 * PI * pow(h, 9) )

for the calculation of surface tension.
2.2.5. Time Integration
We have used leap frog integration scheme for derivative calculation.
2.2.6. Fluid Surface Construction
2.2.6.1. Marching Cube
To visualize free surfaces, we rendered the fluid sim using iso surfaces. We use marching cubes algorithm to triangulate the iso surface. To implement this algorithm,
• We started by initializing a fixed grid based on the min and max boundary conditions for each scene.
• Then, we identified the cells which contained the surface using the distance between the particle position and grid minimum boundary.
• During simulation step, we search for all the cells that contain surface particles and for each cell identified to contain the surface, the triangles are generated.
• We update the grid at each simulation step to identify the cell containing the surface.

2.2.3. Implementation and Results
We have developed two versions of our application:
• Standalone C++ and opengl application.
• Maya Plugin for the fluid simulation.
The details of our plugin and opengl application are provided here.
Mostly, we generated our results with 3000-4000 particles with surface resolution nearly 0.3-0.5. However, it is very slow and when the particles are increased it gets much more slower. Hence, a better implementation technique for particle generation will be required for faster results. Here, are some of our results rendered in Maya(mental Ray and maya software)

References:
1. Particle-Based Fluid Simulation for Interactive Applications.   JMatthias Müller, David Charypar and Markus Gross.
2. Construction of Non-blobby Surface from Particles.   ,Takahiro et al
3. Realtime Particle Based Simulation.   , Stefan Auer