Buoyancy force appears to be calculated incorrectly

[vihren-milev]

In my opinion the calculation of the buoyancy is incorrect (file pfgen.cpp lines 120 and 121).

I think the correct calculation should look like this:

    force.y = liquidDensity * volume *
        (depth - maxDepth - waterHeight) / (2 * maxDepth);

Current code is missing the parentheses around "2 * maxDepth"

[ButchDean]

Have you been seeing any anomalies? The way you are presenting it is a different calculation.

[vihren-milev]

Buoyancy here is being calculated by multiplying the maximum buoyancy (volume times density) by the ratio between submerged volume and total volume. This ratio calculation on page 91 of the book is:
d = (y0 - yW - s) / 2s
Note that 2s (which is 2 * maxDepth) is in the denominator. The acual code correctly divides by two but multiplies by maxDepth which appears to be erroneous to me.

[ButchDean]

I can see what you mean. Why not issue a pull request?

[vihren-milev]

I am just lazy I guess - wanted to report the issue without getting involved at this time :)

[ButchDean]

Pull requests normally get much more attention than simply reporting the issue, as well as making it easier for the author to test. But it's up to you. :)

[vihren-milev]

I might do that... Meanwhile I am finding more issues:
The ParticleContact::resolveInterpenetration applies the interpenetration corrections in the same direction to both particles involved in the collision. In order to separate interpenetrating particles the corrections should be applied in opposite directions.

[vihren-milev]

Ignore that last comment - it is only related to the book. The actual code is correct...

[xwc2021]

But I uset that

float d =1 - (x - (waterHeight - maxDepth)) / (2maxDepth);
force.y = liquidDensity * volume *d;

because I want
when x =waterHeight + maxDepth , then d =0
when x =waterHeight - maxDepth , thed d=1
Geogebra is a good tool to visualize.
https://ggbm.at/VJGGQnSk

[lch32111]

@xiangwei71
But, what if the x > (waterHeight - maxDepth) ? In this case, the buoyancy force should be 0. But, your code will make a extra force for that case.

[xwc2021]

@lch32111
Thanks for your remind. I should clamp it.

[bingbingzhenbang]

force.y = liquidDensity * volume *
        (maxDepth  + waterHeight - depth) / (2 * maxDepth);

[bingbingzhenbang]

[bingbingzhenbang]

depth maxDepth both from center

[lch32111]

I think The explanation and implementation of buoyancy force on this book is not stable.
So, I searched for the better theory and implementation. And I happened to see the Erin Catto Demo and describtion on the Game Programming Gems 6. And other reference is also based on the kinds of theory.

As I'd like to implement particle buoyancy force, I approximate the Volume of a particle as Y-axis height of part of submerged particle. In addition, There is one thing you have to insert, the drag. If you do not implement of drag of buoyancy force, You will see the particle jumping up over the water high. To prevent the effect, You should apply the drag to the buoyancy force.

Consequently, It is working well. Take my code and test. You will get plausible result.

The result GIF:
The box is a particle. I just render it as a box without any meaning.
The blue line means water height.
the white line passing through the box is just showing the position of the particle.

	/**
	 * A force generator that applies a buoyancy force for a plane of 
	 * liquid parallel to XZ plane.
	 */
	class ParticleBuoyancy : public ParticleForceGenerator
	{
		/**
		 * Consider a particle as a sphere, and it has a radius
		 * With this variable, the submerged part of a particle is
		 * determined.
		 */
		real particleRadius;


		/**
		 * The height of the water plane above y = 0. The plane will be
		 * parallel to the XZ plane.
		 */
		real waterHeight;

		/**
		 * The linear Drag of Water
		 */
		real linearDrag;

		/**
		 * The density of the liquid.Pure water has a density of
		 * 1000 kg per cubic meter.
		 */
		real liquidDensity;

	public:
		/** Creates a new buoyancy force with the given parameters */
		ParticleBuoyancy(real radius, real linearDrag, real waterHeight, real liquidDensity = 1000.0f);

		/** Applies the buoyancy force to the given particle.*/
		virtual void updateForce(GPEDParticle* particle, real duration);
	};

void GPED::ParticleBuoyancy::updateForce(GPEDParticle * particle, real duration)
{
	// Calculate the submersion depth.
	real depth = particle->getPosition().y;

	// Check if we're out of the water.
	if (depth - particleRadius >= waterHeight) return;
	glm::vec3 force(real(0.0));
	real partialMass = particle->getMass();

	// Fully Submerged
	if (depth + particleRadius <= waterHeight)
	{
		force.y = liquidDensity * (particleRadius * 2) * real(10.0); // gravity 10
	}
	else // Partially Submerged
	{
		real V = (waterHeight - (depth - particleRadius));
		force.y = liquidDensity * V * real(10.0); // gravity 10
		partialMass *= V / (particleRadius * 2);
	}
	
	// http://www.iforce2d.net/b2dtut/buoyancy 
	// Refer to the tutorial of buoyancy. the drag is approximation.
	// water velocity is assumed by (0, 0, 0)
	// https://box2d.org/downloads/
	// Refer to the buoyancy demo of Eric Catto
	// the water linear drag is normally 5.0, according to the demo
	float dragForce = (partialMass * linearDrag) * (-particle->getVelocity().y);
	force.y += dragForce;

	particle->addForce(force);
}

[moekaz]

Hi, i might be wrong on this but looking at the book's formula i am pretty sure that the force for partial submersion should be:
force.y = density * volume * (waterHeight + maxDepth - depth) / (2 * maxDepth)
which is the opposite of what is calculated for d as if i want to calculate the amount of the box submerged i would get the center of the object (depth in this case) and subtract from that the maxDepth which will give me the position of the very bottom of the object. subtracting that value from the height of the water plane will then give me the distance from the plane to the bottom of the object which is equivalent to the amount of the box submerged and then obviously we would get the ratio of that to the total size of the object which is 2 * maxDepth