# SPH parameters for water like behavior

This topic is 594 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

## Recommended Posts

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 {
float h;
// Rest density
float p0;
// Stiffness
float k;
// Near stiffness
float kNear;
// Linear viscosity
float alpha;
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 &params, 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 &params, 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 &params, 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));
return true;
}
return false;
}

force_inline bool SPHComputeViscosityVelocity(const SPHParameters &params, 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);
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.

Edited by Finalspace

##### Share on other sites

Yeah i know its a hard topic, but if you dig enough you get yourself successful and get pretty results ;-)

I am still trying to find good parameters, but i got more close after setting up a fixed constants independent from the boundary und use this as a reference:

//
// Boundary
//
const float kSPHBoundaryAspect = 16.0f / 9.0f;
const float kSPHBoundaryWidth = 50.0f;
const float kSPHBoundaryHeight = kSPHBoundaryWidth / kSPHBoundaryAspect;
const float kSPHBoundaryHalfWidth = kSPHBoundaryWidth * 0.5f;
const float kSPHBoundaryHalfHeight = kSPHBoundaryHeight * 0.5f;
const static Vec2f kSPHGridOrigin = Vec2f(-kSPHBoundaryHalfWidth, -kSPHBoundaryHalfHeight);

//
// Default constants
//
// @NOTE: kH must be chosen well, everything else depends on this!
// @NOTE: Particle spacing must be slightly smaller than kH, otherwise there will be no interaction!
//
const float kSPHDeltaTime = 1.0f / 60.0f;
const float kSPHKernelHeight = 0.5f;
const float kSPHParticleSpacing = 0.475f;
const float kSPHParticleCollisionRadius = kSPHKernelHeight * 0.5f;

// @TODO: Water like coeffecients
const float kSPHRestDensity = 30.0f;
const float kSPHStiffness = 50.0f;
const float kSPHNearStiffness = 70.0f;
const float kSPHLinearViscosity = 50.0f;

//
// Render constants
//
const float kSPHParticleRenderRadius = kSPHKernelHeight * 0.5f;
const float kSPHVisualPlaneLength = kSPHBoundaryHalfWidth;

//
// Grid
//
// @NOTE: Cell size must be greater than kH
//
const float kSPHGridCellSize = kSPHKernelHeight * 4.0f;
const int kSPHGridCountX = (int)(kSPHBoundaryWidth / kSPHGridCellSize);
const int kSPHGridCountY = (int)(kSPHBoundaryHeight / kSPHGridCellSize);
const int kSPHGridTotalCount = kSPHGridCountX * kSPHGridCountY;
const float kSPHGridWidth = kSPHGridCountX * kSPHGridCellSize;
const float kSPHGridHeight = kSPHGridCountY * kSPHGridCellSize;


##### Share on other sites

Well it would helped if i had implemented the correct formulas from the paper...

Rij is not the same Rji !

Especially when the vectors are used as directions... -_-

Now the parameter tuning can begin.

• ### What is your GameDev Story?

In 2019 we are celebrating 20 years of GameDev.net! Share your GameDev Story with us.

• 13
• 12
• 15
• 11
• 12
• ### Forum Statistics

• Total Topics
634152
• Total Posts
3015842
×

## Important Information

GameDev.net is your game development community. Create an account for your GameDev Portfolio and participate in the largest developer community in the games industry.

Sign me up!