Calculating impulse due to rigid body collision with friction

Started by
15 comments, last by erwincoumans 16 years, 6 months ago
I am programming a simple physics simulation and have been following along with some books about the subject. Everything has gone well so far, although at the moment I'm stuck on how to calculate the impulse when two rigid bodies collide. I'll explain the algorithms that I'm using below. In a frictionless environment I would use this equation to find the impulse that would be applied to both rigid bodies in the collision. J = -Vr(1+e) / {1/m1 + 1/m2 + n . [(r1 x n)/I1] x r1 + n . [(r2 x n)/I2] x r2} J - The scalar impulse Vr - The relative closing velocity of the two bodies e - The coefficient of restitution m1 - The mass of rigid body 1 m2 - The mass of rigid body 2 n - The normal of the collision r1 - The vector from the center of mass of rigid body 1 to the point of collision r2 - The vector from the center of mass of rigid body 2 to the point of collision I1 - Inertia tensor for rigid body 1 I2 - Inertia tensor for rigid body 2 I apply the impulse to the objects as such. v1 = v1 + (Jn)/m1 v2 = v2 + (-Jn)/m2 w1 = w1 + (r1 x Jn)/I1 w2 = w2 + (r2 x -Jn)/I2 v1 - Velocity of rigid body 1 v2 - Velocity of rigid body 2 w1 - Angular velocity of rigid body 1 w2 - Angular velocity of rigid body 2 I have tested this and everything works correctly. The thing that I'm clueless about is when I then decide to throw dynamic/kinetic friction into the equation. So far from the books I'm reading this is what I have been able to figure out. This is the updated equations to update the rigid bodies. v1 = v1 + [Jn + (uJ)t]/m1 v2 = v2 + [-Jn + (uJ)t]/m2 w1 = w1 + {r1 x [Jn + (uJ)t]}/I1 w2 = w2 + {r2 x [-Jn + (uJ)t]}/I2 u - The coefficient of friction t - The tangent normal, this is a unit vector that is perpendicular to the contact normal. The tangent normal is calculated with this equation. t = [(n x Vr) x n] t = t/|t| I can understand the concept of what his happening in the above equations. I also know that I'm going to need to change the impulse calculation equation because of friction, although I'm not sure how this is done. The books I have only skim over this aspect and I'm a little unsure how to continue. I would imagine it is not difficult, although my searching has not turned up any useful information. So my question is, when two rigid bodies collide in a frictional envionrment how does one go about calculating the impulse? Thank for reading. :p
Advertisement
I've learned this looking at the source code of the MrRowl's physics simulators(thanks for that :) ) and that works very nicely. Here's what I do(this is from the rigid sphere dynamics simulator I did(it has no resting contact :( ):

//Solves a collision between 2 spheresvoid PHYSICALWORLD::SolveCollision(const COLLISIONINFO& c){	SPHERE *sA, *sB;        //get pointers to the spheres so that sA points to the sphere with mass        //and sB points to the sphere with 0 mass or if both have 0 mass return	if(sphere[c.indexA]->GetMass() > 0.f)	{		sA = sphere[c.indexA];		sB = sphere[c.indexB];	}	else if(sphere[c.indexB]->GetMass() > 0.f)	{		sA = sphere[c.indexB];		sB = sphere[c.indexA];	}	else		return;	bool massB = sB->GetMass() > 0.f;        //put them at the 'timeOfImpact position' and also linearly interpolate the velocities	sA->SetPosition(sA->GetPrevPosition() + (sA->GetPosition() - sA->GetPrevPosition())*c.timeOfImpact, false);	sA->SetOrientation(sA->GetPrevOrientation() + (sA->GetOrientation() - sA->GetPrevOrientation())*c.timeOfImpact,																									false, false);	sA->SetLinVelocity((sA->GetPosition() - sA->GetPrevPosition())/(sA->GetDeltaTime()*c.timeOfImpact));	sA->SetAngVelocity((sA->GetOrientation() - sA->GetPrevOrientation())/(sA->GetDeltaTime()*c.timeOfImpact));	if(massB)	{		sB->SetPosition(sB->GetPrevPosition() + (sB->GetPosition() - sB->GetPrevPosition())*c.timeOfImpact, false);		sB->SetOrientation(sB->GetPrevOrientation() + (sB->GetOrientation() - sB->GetPrevOrientation())*c.timeOfImpact,																									false, false);		sB->SetLinVelocity((sB->GetPosition() - sB->GetPrevPosition())/(sB->GetDeltaTime()*c.timeOfImpact));		sB->SetAngVelocity((sB->GetOrientation() - sB->GetPrevOrientation())/(sB->GetDeltaTime()*c.timeOfImpact));	}        //compute collision normal	VECTOR3 n = (sA->GetPosition() - sB->GetPosition())/(sA->GetRadius() + sB->GetRadius());	VECTOR3 rA = -n*sA->GetRadius();	VECTOR3 rB = n*sB->GetRadius();        //relative velocity at contact point	VECTOR3 vRel = sA->GetLinVelocity() + (sA->GetAngVelocity() ^ rA);	if(massB)		vRel -= sB->GetLinVelocity() + (sB->GetAngVelocity() ^ rB);        //normal relative velocity 	float nvRel = n*vRel;        //the bodies are going to penetrate	if(nvRel < 0.f)	{                //compute and apply normal impulse. there's no need to include                //the inertia tensor here because we're working with spheres only		float numerator = -(1.f+(sA->GetElasticity() + sB->GetElasticity())*0.5f)*nvRel;		float denominator = sA->GetInvMass() + sB->GetInvMass();		float j = numerator/denominator;		VECTOR3 jn = n*j;		sA->ApplyImpulse(jn, rA);		if(massB) sB->ApplyImpulse(-jn, rB);                //FRICTION IMPULSE                //compute the relative tangential velocity at the contact point		VECTOR3 tVel = vRel - n*nvRel;                //tangential speed		float tSpd = tVel.LengthSq();                //the tangential speed is greater than 0 then we must apply friction impulse		if(tSpd > 0.0000001f)		{			tSpd = sqrtf(tSpd);                        //friction normal to use in the impulse equation. Use the                        //minus sign so that it points in the opposite direction                        //of the tangential velocity vector			VECTOR3 fn = -tVel/tSpd;                        //compute the denominator of the impulse equation the                        //same way you do when computing for normal impulse                        //but you use the friction normal instead                        //the inertia is included here			denominator = sA->GetInvMass() + fn*((sA->GetInvInertia()*(rA^fn))^rA);			if(massB) denominator += sB->GetInvMass() + fn*((sB->GetInvInertia()*(rB^fn))^rB);                        //the friction impulse. the tangential speed is the numerator			float fj = tSpd/denominator;                        //as when working with forces, if the friction impulse                        //is greater than the normal impulse times the static                        //friction coefficient, which means a slide, the                        //friction impulse is the normal impulse times the                        //dynamic friction coefficient else there's no change                        //to do there			if(fj > j*(sA->GetStaticFriction() + sB->GetStaticFriction())*0.5f)				fj = j*(sA->GetDynamicFriction() + sB->GetDynamicFriction())*0.5f;                        //compute the friction impulse vector			VECTOR3 fjv = fn * fj;                        //apply the impulses that same way			sA->ApplyImpulse(fjv, rA);			if(massB) sB->ApplyImpulse(-fjv, rB);		}		sA->ApplyForce(gravity*sA->GetMass());		sA->ApplyDragForce(fluidDensity);		sA->Integrate(sA->GetDeltaTime()*(1.f - c.timeOfImpact));		sA->ZeroOutForces();		if(massB)		{			sB->ApplyForce(gravity*sB->GetMass());			sB->ApplyDragForce(fluidDensity);			sB->Integrate(sB->GetDeltaTime()*(1.f - c.timeOfImpact));			sB->ZeroOutForces();		}	}}


Thats all, Hope this helps.
.
IIRC a good derivation of a contact with friction can be found e.g. in the PhD thesis of B. Mirtich:

http://www.kuffner.org/james/software/dynamics/mirtich/index.html


The method in MrRowl's Jiglib comes IIRC from the Guendelman paper. His PhD and the original paper can be found here:

http://graphics.stanford.edu/~erang/papers/thesis/thesis.pdf
http://graphics.stanford.edu/papers/rigid_bodies-sig03/


Both good lectures...


HTH,
-Dirk
Quote:Original post by KiroNeem

J = -Vr(1+e) / {1/m1 + 1/m2 + n . [(r1 x n)/I1] x r1 + n . [(r2 x n)/I2] x r2}

J - The scalar impulse
Vr - The relative closing velocity of the two bodies
e - The coefficient of restitution
m1 - The mass of rigid body 1
m2 - The mass of rigid body 2
n - The normal of the collision
r1 - The vector from the center of mass of rigid body 1 to the point of collision
r2 - The vector from the center of mass of rigid body 2 to the point of collision
I1 - Inertia tensor for rigid body 1
I2 - Inertia tensor for rigid body 2

I apply the impulse to the objects as such.

v1 = v1 + (Jn)/m1
v2 = v2 + (-Jn)/m2
w1 = w1 + (r1 x Jn)/I1
w2 = w2 + (r2 x -Jn)/I2

v1 - Velocity of rigid body 1
v2 - Velocity of rigid body 2
w1 - Angular velocity of rigid body 1
w2 - Angular velocity of rigid body 2

I have tested this and everything works correctly. The thing that I'm clueless about is when I then decide to throw dynamic/kinetic friction into the equation. So far from the books I'm reading this is what I have been able to figure out.

This is the updated equations to update the rigid bodies.

v1 = v1 + [Jn + (uJ)t]/m1
v2 = v2 + [-Jn + (uJ)t]/m2
w1 = w1 + {r1 x [Jn + (uJ)t]}/I1
w2 = w2 + {r2 x [-Jn + (uJ)t]}/I2

u - The coefficient of friction
t - The tangent normal, this is a unit vector that is perpendicular to the contact normal.

The tangent normal is calculated with this equation.

t = [(n x Vr) x n]
t = t/|t|




I started out using the very same equations for friction, my source was the book "Physics for Game Developers - David M. Bourg". But I realized there is a problem with those equations; I was not able to simulate a tire(squashed sphere) rolling over a plane. Its easy to see why the equations would fail in this case:

J = -Vr.n(1+e) / {1/m1 + 1/m2 + n . [(r1 x n)/I1] x r1 + n . [(r2 x n)/I2] x r2}

It follows from the equation above that J is directly proportional to the normal velocity component. So if Vr is purely tangential to n - as in the case I described above, then J reduces to zero. And so does the change in linear and angular velocities.
Obviously this is not correct - if you gently place a rotating tire on the floor, you expect its linear velocity (forward) to increase and its angular velocity to decrease.

I've modified the impulse equation to handle this situation. Its best to treat the impulse as a vector rather than a scalar.

J = -Vr.n(1+e)n - Vr.t(u)t
ImpN = J; ImpN.Normalize();
J /= {1/m1 + 1/m2 + ImpN . [(r1 x ImpN)/I1] x r1 + ImpN . [(r2 x ImpN)/I2] x r2}

and ImpN = direction of impulse

Then,

dV = J/m
dL = rxJ

dL = change in angular momentum


These seem to work well enough..
"I was not able to simulate a tire(squashed sphere) rolling over a plane. Its easy to see why the equations would fail in this case..."

You can't expect equations derived for rigid bodies (i.e. a rigid sphere rolling on a rigid plane) to result in any rolling friction, since rolling friction comes from energy lost through compression/decompression of the objects.

The Guendelman paper describes a (presumably) more physically correct method for calculating rolling friction.

Quote:Original post by MrRowl
"I was not able to simulate a tire(squashed sphere) rolling over a plane. Its easy to see why the equations would fail in this case..."

You can't expect equations derived for rigid bodies (i.e. a rigid sphere rolling on a rigid plane) to result in any rolling friction, since rolling friction comes from energy lost through compression/decompression of the objects.

The Guendelman paper describes a (presumably) more physically correct method for calculating rolling friction.


Thanks MrRowl, can you please provide a link to necessary paper for simulating rolling friction?
I see your point, guess I need to spend more time reading papers than coding :)

KiroNeem, apologies if I put you off track a bit.
xissburg - Thanks for the code, it has given me a little better of an idea of what I should be doing. Although seeing as angular motion is not taken into account It will still not completely work for what I'm trying to do.

DonDickieD - I started to print that PDF you linked to without even looking at how many pages it was. I came back to my computer with about two hundred pages on my desk and my printer beeping at me that it needed more paper. :p

This actually looks like a very good reference and I've begun to read over it, thank you very much. Once I find what I'm looking for I will post back my findings.

Chand81 - Actually, I might have some insight for you on what you are trying to do. From what I understand there are several variations of friction, the main ones in question are.

Static - This is a force that opposes motion of a body that is resting on another. This is related to rolling friction, but it is static friction that will cause a ball to roll across the ground when it has angular velocity.

Dynamic/Sliding/Kinetic - This is a force that will dampen motion of a body that is moving across another.

Rolling - This is a force that will slow a ball or tire from rolling due to deformation of the objects. This does not actually cause the linear motion as static friction does.

Here are some resources I have been reading over
http://en.wikipedia.org/wiki/Friction
http://en.wikipedia.org/wiki/Rolling_resistance
http://webphysics.davidson.edu/faculty/dmb/PY430/Friction/rolling.html
Basically the friction impulse works like the collision impulse.

We seek a momentum conserving impulses...

v1' = v1 + M1^-1 * P1
w1' = w1 + I1^-1 * L1 = I1^-1 * ( r1 x P1 )

v2' = v2 + M2^-1 * P2
w2' = w2 + I2^-1 * L2 = I1^-1 * ( r2 x P2 )

...such that relative tangential velocity vanishes...

( v2 + w2 x r2 - (v1 + w1 x r1) ) * t = 0

...also knowing the direction of the impulses

P1 = -lambda * t
P2 = lambda * t = -P1

v := linear velocity
w := angualr velocity (read omega)
v' := linear velocity after impulse (read post-velocity)
w' := angular velocity after impulse
P := linear impulse applied at center of mass
L := angular impulse

r := vector from center of mass to contact point

w x r is the cross product of w and r ( read cross(w, r) )



If you have computed lambda this way you have found an impulse that eliminates the entire relative tantential velocity. Unfortunately there is a coupling to the normal impulse that means that you maybe cannot activate the full impulse, so you need to check this and clamp the impulse if required:

lambda < mu * lambda_n


The procedure for finding the impulses is usually the same. You blindly write down the impulse equation. Then you write the kinematic constraint that the post velocities must satisfy. Finally you usually also know the direction of impulse.


HTH,
-Dirk
Chand:
"Thanks MrRowl, can you please provide a link to necessary paper for simulating rolling friction?"

Danny is refering to he second papers in my first post...
"DonDickieD" - Thanks for your reply about how you went about handling friction. Based off of your implementation I was able to get something that is a lot closer to what I'm trying to do. Though, shouldn't time be taken into account when calculating how much friction acts on the tangent impulse? Or is that taken care of inherently by the fact that you are working with velocity instead of forces?

[Edited by - KiroNeem on September 29, 2007 3:22:42 PM]

This topic is closed to new replies.

Advertisement