In the past i have successfully implemented a 2d fluid simulations in dozens of languages - except for c++ (All suffered bad parameter tuning).
Now i have ported / re-implemented it into c++ and all math are identical to the paper i used (without the elastic spring stuff):
- Particle-based Viscoelastic Fluid Simulation by Clavet
All math works but i cannot find a suitable combination of parameters to let it behave more like water.
So i want to ask some professionals here, how to find / calculate the right properties to get it more like water.
How to choose the Rest-Density?
Stiffness and Near-Stiffness?
Should the smoothing length be the same as the cell size?
Which scaling distance for initial configuration (Boundary condition) so that it starts to merge without exploding
Using a timestep of 1 / 60 is fine? (In the paper they use 1 / 30).
The parameters from the paper do not work for me in 2D (The simulation simply explodes):
- Rest-Density = 10
- Stiffness = 0.004
- Near Stiffness = 0.01
- Linear viscosity = 0.3
There must be a way to visually see how these parameters will affect the particles, so that you choose parameters wisely.
Does someone know this stuff?
Right now i decide the cell size first and then use this as a smoothing radius directly and tune all parameters based on that like this (It works but it does not look great):
const float kDeltaTime = 1.0f / 60.0f;
// Decide boundary first based on a fixed width
const float kSPHBoundaryAspect = 16.0f / 9.0f;
const float kSPHBoundaryWidth = 10.0f;
const float kSPHBoundaryHeight = kSPHBoundaryWidth / kSPHBoundaryAspect;
const float kSPHBoundaryHalfWidth = kSPHBoundaryWidth * 0.5f;
const float kSPHBoundaryHalfHeight = kSPHBoundaryHeight * 0.5f;
// Decide cells second based on a fixed number of cells in the x dimension
const static Vec2f kSPHGridOrigin = Vec2f(-kSPHBoundaryHalfWidth, -kSPHBoundaryHalfHeight);
const int kSPHGridCountX = 64;
const int kSPHGridCountY = (int)(kSPHGridCountX / kSPHBoundaryAspect);
const int kSPHGridTotalCount = kSPHGridCountX * kSPHGridCountY;
const float kSPHGridCellSize = kSPHBoundaryWidth / (float)kSPHGridCountX;
struct SPHParameters {
// Smoothing kernel radius
float h;
// Rest density
float p0;
// Stiffness
float k;
// Near stiffness
float kNear;
// Linear viscosity
float alpha;
// Quadratic viscosity
float beta;
SPHParameters() {
h = kSPHGridCellSize * 1.0f;
p0 = 1000.0f;
k = 0.2f;
kNear = 0.3f;
alpha = 5.0f;
beta = 0.5f;
}
};
force_inline void SPHComputeDensity(const SPHParameters ¶ms, const Vec2f &position, const Vec2f &neighborPosition, float outDensity[2]) {
Vec2f Rij = position - neighborPosition;
float rijSquared = Vec2Dot(Rij, Rij);
// TODO: Make it branch-free
if (rijSquared < params.h * params.h) {
float rij = sqrtf(rijSquared);
float q = rij / params.h;
float oneminusq = 1 - q;
outDensity[0] += (oneminusq * oneminusq);
outDensity[1] += (oneminusq * oneminusq * oneminusq);
}
}
force_inline void SPHComputePressure(const SPHParameters ¶ms, const float density[2], float outPressure[2]) {
outPressure[0] = params.k * (density[0] - params.p0);
outPressure[1] = params.kNear * density[1];
}
force_inline bool SPHComputeDelta(const SPHParameters ¶ms, const Vec2f &position, const Vec2f &neighborPosition, const float pressure[2], const float deltaTime, Vec2f *outDelta) {
// TODO: Make it branch-free
Vec2f Rij = position - neighborPosition;
float rijSquared = Vec2Dot(Rij, Rij);
if (rijSquared < params.h * params.h) {
float rij = sqrtf(rijSquared);
float q = rij / params.h;
Vec2f n = Vec2Normalize(Rij);
float oneminusq = 1 - q;
float term = 0.5f * deltaTime * deltaTime * (pressure[0] * oneminusq + pressure[1] * (oneminusq * oneminusq));
*outDelta = Vec2Hadamard(term, n);
return true;
}
return false;
}
force_inline bool SPHComputeViscosityVelocity(const SPHParameters ¶ms, const Vec2f &position, const Vec2f &neighborPosition, const Vec2f &velocity, const Vec2f &neighborVelocity, const float deltaTime, Vec2f *outVel) {
// TODO: Make it branch-free
Vec2f Rij = position - neighborPosition;
float rijSquared = Vec2Dot(Rij, Rij);
if (rijSquared < params.h * params.h) {
float rij = sqrtf(rijSquared);
Vec2f n = Vec2Normalize(Rij);
float u = Vec2Dot(velocity - neighborVelocity, n);
if (u > 0) {
float q = rij / params.h;
float term = 0.5f * deltaTime * (1 - q) * (params.alpha * u + params.beta * u * u);
*outVel = Vec2Hadamard(term, n);
return true;
}
}
return false;
}
I have looked into dozens of other sph implementations and all do it differently. Some use the water volume for calculate the smoothing lenght, some calculate it from the number of max particles, some use a fixed radius, some use a big radius and scale it down for rendering only.
Video will come soon.
