As the name suggests this approach is based on simulating particles, or in other words tracing the state of the fluid at a finite number of dynamic points ( also known as position based dynamics). Only external forces are integrated and other forces are solved through constraints between particles. A cloth-piece for example can be realised by apply a distance constraint between neighboring particles, keeping them at the same distance from each other. This approach also allows interaction between different type of particles by defining constraints between these particles. A collision with a sphere for example can be also realised by a distance constraint between the sphere-particle and all other particles. Similar constraints are used to create a particle system with fluid like behavior which is outlined below.
Source code and some technical details can be found here.
Particles are tracked throughout the simulation. Explicit force integration would require to model and add all forces in order to apply them to each particle. As some forces depend on time this approach tends to be unstable unless very small timesteps are used. An alternative approach is to separate forces into external forces which can be applied directly and internal forces which need be modeled through constraints between particles, which are solved after external forces have been applied. The algorithm proposed by Clavet05 follows this idom:
Step 2 applies a viscosity impulse to all particles correcting their velocities before it’s used to predict the next position. In Step 4 the internal forces are modeled through a spring constraint which can be thought of as a particle-density depended impulse to model the cohesive attraction as well as incomprehensibility of fluids. Both steps are discussed below in greater detail.
As fluids are a continuous medium (at least for the scale we are interested in) and particles are inherently discrete
in nature smoothed particle hydrodynamics(SPH) interpolation is used. Based on an interaction
radius h and a contribution function $\omega$ (kernel-function) the properties of each particle is distributed.
For example to evaluate the density at the position of a particle we have to consider all particle within the
interaction radius und sum their $\omega$ depended contributions.
Not only the density needs to be blurred but also the velocity to create a smooth velocity field. After all external forces, here only gravity, have been applied to each particle individually their contribution to other particles within their interaction radius needs to be accounted for. Thus we need to calculate a impulse for each particle pair which adds the forces contribution to the other particle. For this impulse we only consider the part of the velocity that is directed towards the other particle. In other words we subtract the neighbor velocity from the particles velocity and project the result on the normalized vector between the particles (see figure).
In order to cover a wide variety of viscous fluids two control parameters $ \alpha,\beta $ are introduced. These control how much of the impulse is added to the other particle.
struct Particle{
vec3 prev; // previous Position
vec3 p; // Position
vec3 v; // Velocity
}
applyVisco(float dt){
for (auto& par: allParticles){
for( auto neigh: getNeighbors(par)){
vec3 norm = normalize(par.p - neigh.p);
float q = length(par.p - neigh.p) / h; // interaction radius h
float u = dot((par.v - neigh.v), norm);
if(u > 0 ){ // if particles are moving towards each other
vec3 dv = dt * (1 - q) * (alpha*u + beta* u*u) * norm;
par.v -= dv*0.5f; neigh.v += dv*0.5f;
} } } }
The impulse dv is then scaled by time and the a linear kernel function making the contribution of
neighbors closer to the particle more relevant. Neighbors are defined to be all particle within
the interaction radius $ h $ thus $ q \in [0,1) $ is the domain of the kernel function $ (1 - q) $.
Based on the corrected velocity a prediction for the particle position can be made which is the basis to solve constraints between particles.
To model cohesive attraction and incomprehensibility we calculate the density $ \rho $ at every particles position and
a displacement term is generated from that.
As before the contribution of other particles within h need to be considered. If the calculated
density $\rho $
is lower than a user defined rest density $ \rho_0 $ neighboring particle will be pulled towards the particle.
If the density $\rho$ is higher the neighboring particles will be pushed away.
for(auto& par: allParticles){
<</* Compute density of particle and resulting displacements*/>>
<</* Apply displacement due to pressure */>>
}
However if only a single density is considered a small amount of particles may form a cluster at rest density and thus prohibiting these particles from interacting with the rest of the fluids particles. To avoid this clustering a second near density $\rho_{near}$ is calculated which exhibits a exclusively repulsive impulse to prevent particles from clustering hence the name double relaxation.
<</* Compute density of particle and resulting displacements*/>>
for(auto& neigh: getNeighbors(par) ){
float q = length(par.p - neigh.p) / h;
float q2 = (1 - q) * (1 - q);
rho = rho + q2;// density
rho_near = rho_near + q2 * (1 - q); // near density
}
// Pressure, basically a spring force with rest position at rho_0
float Pr = (rho - rho_0) * k;
// near Pressure, exclusively repulsive
float P_near = rho_near * k_near;
Note that the near density $ \rho_{near} $ is calculated using a cubic kernel $ (1-q)^3 $ where as the density $ \rho
$ is
calculated using a quadratic kernel $ (1-q)^2 $. As can be seen in the figure below the fall of in cubic kernel is
faster thus
particle with a large relative distance q barely contribute to the near_density value in term resulting a minor
displacement.
The calculated displacements need to be applied as a radial force to each neighbors depending on it’s distance. The near pressure $\rho_{near}$ is applied using a quadratic kernel whereas the the pressure $\rho$ is weighted by a linear kernel. Once again this ensures that particle further away are effected mainly by the pressure $\rho$.
<</* Apply displacement due to pressure */>>
vec3 dx;
for( auto& neigh: getNeighbors(par) ){
float q = length(par.p - neigh.p) / h;
vec3 norm = normalize(par.p - neigh.p);
float Pr_cur = Pr * (1 - q);
float Pr_near_cur = P_near * (1 - q) * (1 - q);
glm::vec3 D = dt * dt * (Pr_cur + Pr_near_cur) * norm;
neigh.p += 0.5* D;
dx -= D;
}
par.p += dx;
After the internal forces have been processed external constraints like collision constraint need to be solved. Since the last constraint solved may undo previous correction of other constraints collision should be processed last to ensure that the final position does not include any obstacle intersections.
Based on the previous and current position the initial velocities are corrected to ensure more accurate predictions in the next step.
As evident from the previous sections we frequently need to get a list of neighbors of a given particle. One could
simply go through all particles and test if the distance is less than the interaction radius h. However this would
result in a quadratic runtime for both the double relaxation as well as the velocity smoothing.
To reduce this complexity a spacial hash map is used. The hashfunction ensures that a certain region (box) is mapped to the same hashvalue thus they are put in the same bucket. The bucket spacing is at least as large as the interaction radius $h$ which ensures that we only need to check the particle that are in the 27 buckets around the particle of interest.
The hashmap is based on a std::multimap with a custom hashfunction.
inline int64_t hashFunction(vec3 pos) const{
u_int16_t x = pos.x * iGridSpacing;
u_int16_t y = pos.y * iGridSpacing;
u_int16_t z = pos.z * iGridSpacing;
if (pos.x < 0.0) { x--; }
if (pos.y < 0.0) { y--; }
if (pos.z < 0.0) { z--; }
int64_t out = 0;
out += x;
out = out << 16;
out += y;
out = out << 16;
out += z;
}
Given a function $$ f : \mathbb{R}^3 \to \mathbb{R} $$ the marching cube algorithm creates a surface in the domain at a given image value. In our example we use, once again, a continuous particle-density function by using the SPH interpolation. Then we define a $\alpha$ goal-density above zero where the fluids surface will be drawn. The region of space that includes the particles is divided into cubes, like a 3D grid. For every point in that grid the corresponding density is calculated using a quadratic kernel and stored in a 3D array.
float Fluid::getDensity(vec3 pos) const{
float sum=0;
for(auto& neigh: getNeighbors(pos)){
float d = length(neigh.p-pos);
sum += (1-(d/h))*(1-(d/h));
}
return sum;
}
Once all densities have been calculated they are compared to the user-defined property $$ \alpha $$. Every vertex with a density higher than $ \alpha$ is marked the others are left unmarked. Now we go through every cube and look for to two vertices a,b where one of them is marked and the other is not. When we found such vertices we know that our goal-density $$ \alpha $$ has to be on the edge between these two vertices and as such the surface that we wish to create also has to pass through that edge.
As a cube has 8 vertices we can represent all possible configuration as a 8bit number. The 8bit number can be used as an index into a look up table that precomputed all edges which the vertices lay on for each configuration. Linear interpolation based on the density value on the cubes-vertices gives us the exact position on the edges of interest. Now the vertices are at the correct position and can be stored for rendering.
As could be seen in the demo video the particles exhibit a clearly fluid like behavior. Further i was able to simulate enough particles to create a somewhat interesting scene. Also demonstrated was simple interactivity with a sphere as it is the easiest to implement. However by keeping track torque and orientation one can implement interaction with any rigid body.
Due to multithreading and the spacial hash map i was able to simulate about 9000 particles on an 4-core i7-4790K at 25-40fps. The surface creation with the marching cube algorithm (15x40x80 grid) takes about 12ms each frame to compute thus the number of renderable particles drops to about 5000.
As could be seen in the demo fairly large particle spacing (interaction radius) was used in order to reduce the number of particles needed to fill a certain volume. In combination with the suboptimal collision constraint for the bounding volume this reintroduces many problems that the approach originally set out to solve.
If a particle moves outside the bounding volume their position is clamped by the collision constraint.
This leads to a velocity singularity on the edge of the bounding volume. Assume the are looking at the edge along the
x-axis and the y-axis which ensures that all particle satisfy par.p.x > bMin.x and par.p.y > bMin.y.
This is done by setting the x and y values of each particle that is outside to par.p = vec3(bMin.x,bMin.y,par.z).
Based on this and the previous position the new velocity is calculated assuming that this is the second timestep
where the particle was corrected due to the collision constraint the resulting velocity will be zero in the
x and y dimension. Further as all particle along that edge have the same x,y position their displacement due to
the density calculation will only be along the z-axis. This prevents particles from flowing into lower density
positions. As the amount of trapped particles depends on how many particle get corrected it also introduces a
timestep dependency.
Non the less I was able to achieve somewhat fluid like behavior while learning about position based dynamics, basics SPH-fluid systems, marching cube and corresponding surface creation from a particle system, acceleration structures like the spacial hash map as well as multithreading with openMP.