Hello,
I'm working on a SPH simulation demo based on the paper, Particle-based Viscoelastic Fluid Simulation, 2005 Clavet, and after translating the double density relaxation and viscosity functions into code, the simulation appears to work, with swirling movements and correct pressures through a U pipe. However, I noticed that my particles are not clumping together like they should in other videos, especially visible when gravity is disabled. I have tried tweaking the rest density, stiffness, near-stiffness and radius constants but I can't seem to get the desired parameters to make my particles clump.
Desired effect:
Here's my videos.
If anyone can offer me a clue of what is going on here I really apperciate it.
Project: https://drive.google.com/file/d/0B2macuhZeO18N3dJLUl5VU5ZVk0/edit?usp=sharing
Here are some relevant parts of my implementation:
public float kParticleRadius = 15.0f;
public float kGravity = 400.0f;
public float kLinearViscocity = 0.001f;
public float kQuadraticViscocity = 0.01f;
public float kRestDensity = 1000.0f;
public float kStiffness = 10.0f;
public float kStiffnessNear = 10.0f;
Update loop:
ApplyExternalForces(timeStep);
ApplyViscosity(timeStep);
AdvanceParticles(timeStep);
grid.UpdateNeighbours(listParticles);
DoubleDensityRelaxation(timeStep);
ResolveCollisions(timeStep);
UpdateVelocity(timeStep);
ColorAndRemoveStray();
ApplyViscosity():
void ApplyViscosity(float dt)
{
foreach (Particle p in listParticles)
{
foreach (Particle n in p.Neighbours)
{
float q = (p.Position - n.Position).Length / p.Radius;
if (q < 1.0f)
{
Vector2 vn = (n.Position - p.Position).Normalized();
float u = Vector2.Dot(p.Velocity - n.Velocity, vn);
if (u > 0.0f)
{
Vector2 i = 0.5f * dt * (1 - q) * (kLinearViscocity * u +
kQuadraticViscocity * u * u) * vn;
p.Velocity -= i;
n.Velocity += i;
}
}
}
}
}
DoubleDensityRelaxation():
void DoubleDensityRelaxation(float dt)
{
foreach (Particle p in listParticles)
{
// Sample neighbouring particle's density
p.Density = p.DensityNear = 0.0f;
foreach (Particle n in p.Neighbours)
{
float q = (p.Position - n.Position).Length / p.Radius;
if (q < 1.0f)
{
p.Density += (1 - q) * (1 - q);
p.DensityNear += (1 - q) * (1 - q) * (1 - q);
}
}
test = p.Density;
// Push and pull particles
// If particle's density > rest density, it will be pushed away
p.Pressure = kStiffness * (p.Density - kRestDensity);
p.PressureNear = kStiffnessNear * p.DensityNear;
Vector2 dx = Vector2.Zero;
foreach (Particle n in p.Neighbours)
{
Vector2 v = p.Position - n.Position;
float q = v.Length / p.Radius;
if (q < 1.0f)
{
Vector2 displacement = 0.5f * dt * dt * (p.Pressure * (1 - q) +
p.PressureNear * (1 - q) * (1 - q)) * v.Normalized();
n.Position += displacement;
dx -= displacement;
}
}
p.Position += dx;
}
}
