Jump to content
  • Advertisement
Sign in to follow this  
Finalspace

SPH parameters for water like behavior

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

If you intended to correct an error in the post then please contact us.

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 {
	// 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 &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));
		*outDelta = Vec2Hadamard(term, n);
		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);
			*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.

Edited by Finalspace

Share this post


Link to post
Share on other sites
Advertisement

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;
const float kSPHQuadraticViscosity = 0.1f;

//
// 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 this post


Link to post
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.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!