Calculating impulse due to rigid body collision with friction

Started by
15 comments, last by erwincoumans 16 years, 6 months ago
Sorry, I don't understand the question. Let me describe how you use my above formulas to deal with multiple simultaneous contacts. Maybe this helps:

for ( all contacts c )  Vec3 t1, t2;  PlaneSpace( c.Normal, t1, t2 );   // Find two mutally ortohgonal vectors    ApplyNormalImpulse( c );  ApplyFrictionImpulse( c, t1 );  ApplyFrictionImpulse( c, t2 );


One gotcha is the clamping of the normal and friction impulse. The straight forward way is to only to apply the normal impulse when then you have an approaching relative velocity at the contact point, but this naiive approach will result in jitter (see e.g. Erin Catto's Box2D slides for a more sorrow explanation). What actually is happening is that you are refining you global solution (the true contact impulse/force) by applying sequential delta impulses. So you clamp against the accumulated impulse and not the delta impulse. The same holds for the friction impulse, which is dependent on the normal impulse. So here again you accumulate the impulse and clamp this against the accumulated normal impulse.


Maybe you look at the mentioned papers by Erin Catto and let me know if you have further questions...


HTH,
-Dirk
Advertisement
DonDickieD - I think there may have been a little misunderstanding about the last question. No problem though, after doing a bit more reading and testing I think the answer is what I had assumed earlier. However, I do again have to say thank you for all the help, honestly any piece of information can be helpful.

As far as the friction equation goes I believe I have finally figured it out, below is the algorithm. The main different from what I first posted is that the friction impulse is handled separately then the contact impulse. If anyone can verify what I'm doing is correct I would appreciate it.

...First we calculate the frictionless impulse

if(|Vrt| > u*J) //Dynamic friction
Jt = -u*J;
else //Static friction
Jt = -|Vrt|;

Jt = Jt/{1/m1 + 1/m2 + n . [(r1 x t)/I1] x r1 + n . [(r2 x t)/I2] x r2}

Jt - The scalar friction impulse
Vrt - The tangent component of the relative velocity
m1 - The mass of rigid body 1
m2 - The mass of rigid body 2
t - The tangent 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

Then you apply the friction impulse similar to the frictionless impulse.

v1 = v1 + (Jtt)/m1
v2 = v2 + (-Jtn)/m2
w1 = w1 + (r1 x Jtt)/I1
w2 = w2 + (r2 x -Jtt)/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

This gives results that are really close to what I want, although there is one outstanding issue. The friction impulse will not apply to rotation around the relative contact point vector. This happens because of how the relative velocity is calculated. For example, if I have sphere spinning along the Y Axis on the ground then it will continue to spin forever.

Either this means that the above equation is wrong, or I just need to handle the friction for the rotation around the relative contact point separately (which is a simple task). As you can see below the cross product will basically cancel out any rotation along the axis created by the relative contact point.

Vr = (v1 + r1 x w1) - (v2 + r2 x w2)

Vr - The relative contact velocity

- Chase
Quote:Original post by DonDickieD
ApplyNormalImpulse( c );
ApplyFrictionImpulse( c, t1 );
ApplyFrictionImpulse( c, t2 );


I'm pretty sure I managed to implement accumulated/clamped friction impulses with just one friction calculation... though I'm not exactly sure now! If so, it's in the jiglib source.
Quote:Original post by KiroNeem
This gives results that are really close to what I want, although there is one outstanding issue. The friction impulse will not apply to rotation around the relative contact point vector. This happens because of how the relative velocity is calculated. For example, if I have sphere spinning along the Y Axis on the ground then it will continue to spin forever.

The Guendelman paper describes a technique to handle Spinning & Rolling friction. The idea is very similar to the static & kinetic friction impulse, but this time a torque impulse is used - will only affect angular velocity and not linear velocity.


Quote:Original post by MrRowl
Quote:Original post by DonDickieD
ApplyNormalImpulse( c );
ApplyFrictionImpulse( c, t1 );
ApplyFrictionImpulse( c, t2 );


I'm pretty sure I managed to implement accumulated/clamped friction impulses with just one friction calculation... though I'm not exactly sure now! If so, it's in the jiglib source.


Thats true, I've seen the jiglib source as well. The idea is pretty simple, just project the velocity vector onto the contact surface to get the tangent component. The friction impulse is applied along the reverse tangent velocity direction.
I'm myself little confused why Dirk's implementation splits the friction impulse into 2 orthogonal components!!?


Here is my implementation, actually 2 implementations one based on the Guendelman paper and the other similar to JigLib. (both lot easier to read as currently only single dynamic object is considered)
Both approaches seem to give very similar results, I wonder if I have implemented the Guendelman approach incorrectly ??


////////////////////////////////////////////////////////////////////////////MrRowls Approach//////////////////////////////////////////////////////////////////////////	gePt3Df Jn; float jn;	//Normal Impulse Vector & Magnitude	float vn = Rvel.Dot(N);	//Normal Speed	gePt3Df Vn = N * vn;	//Normal Velocity		if(vn<0)	//Approaching velocity	{		jn = vn*(-(1+elasticity));		jn /= massInv + ((cfg.iITwInv*(R.Cross(N))).Cross(R)).Dot(N);		Jn = N*jn;		dV += (Jn*massInv);		dL += R.Cross(Jn);		dw = cfg.iITwInv*dL; 		//Rvel needs to be recomputed with the normal impulse having taken effect.		Rvel = (cfg.iAngVel+dw).Cross(R);			Rvel += cfg.iLinVel + dV; 		vn = Rvel.Dot(N);		Vn = N * vn;				//Normal Velocity		gePt3Df Vt = Rvel - Vn;			//Tangent Velocity		float vt = Vt.Length();			//Tangent speed		if( vt > 1e-7f)		{			gePt3Df Jt; float jt;		//Tangent Impulse Vector & Magnitude			gePt3Df T = Vt/vt;		//Tangent			jt = vt;			jt /= massInv + ((cfg.iITwInv*(R.Cross(T))).Cross(R)).Dot(T);						if (jt > jn*friction)				//Normal impulse is not pressing down hard enough.				//So the tangent impulse should allow slide.				jt = jn*friction;			Jt = T * (-jt);			dV += Jt*massInv;			dL += R.Cross(Jt);			dw = cfg.iITwInv * dL; 		}	}



////////////////////////////////////////////////////////////////////////////Guendelman Approach//////////////////////////////////////////////////////////////////////////	float vn = Rvel.Dot(N);	if(vn<0)	{		gePt3Df Vn = N * vn;	//Normal Velocity		//Compute J such that Rvel's tangent component vanishes		gePt3Df J;		TMatrixf K(massInv,0,0,0, 0,massInv,0,0, 0,0,massInv,0, 0,0,0,1);		TMatrixf RSkew(0); RSkew.CreateSkew(R);		K += RSkew.Transpose() * cfg.iITwInv * RSkew;		J = -Vn*elasticity - Rvel;		J = K.Inverse() * J;		float jn = J.Dot(N);		//Impulse normal component		float jt = (J - N*jn).Length();	//Impulse tangent component				if (jt > friction*jn)		{			gePt3Df Vt = Rvel - Vn;	//Tangent Velocity			float vt = Vt.Length();	//Tangent Velocity magnitude			gePt3Df T = Vt/vt;	//Tangent			gePt3Df NUT = N - T*friction;			float j = -vn*(elasticity + 1);			j /= (K*NUT).Dot(N);			J = NUT*j;		}		dV += J*massInv;		dL += R.Cross(J);		dw = cfg.iITwInv*dL; 	}
Quote:I'm pretty sure I managed to implement accumulated/clamped friction impulses with just one friction calculation... though I'm not exactly sure now! If so, it's in the jiglib source.


You did :-) You basically split the relative linear velocity at the contact into a normal and tangential component as in the Guendelman paper. I follow more the PGS approach implemented in terms of impulses and had some concerns with the warmstarting of the solver when not projecting onto a basis in the tangent plane. Actually I use central friction in practise so this is not a real issue anyway.

For those who are interested:
For central friction you only apply normal impulses at the contact points and the friction impulse is applied at the center of the contact patch/geometry. Since this doesn't account for twist (relative rotation around the contact normal) between the faces in contact you need to apply an additional (angular) twist impulse. For a contact patch of 4 points this gets you from 12 constraints (or impulses) down to seven.
I remember now why I use two friction impulses. Basically I can't see how you can accumulate the friction impulse when the basis is changing during the iteration?
You see what I mean? In my opinion you need to project onto a staionary basis, while otherwise you are accumulating vectors (length) with changing direction which is not correct in my opinion.

Anyway it is good for the quality of the friction to align the first tangential axes with the relative velocity in contact plane (if it exists).
Quote:Original post by DonDickieD
I remember now why I use two friction impulses. Basically I can't see how you can accumulate the friction impulse when the basis is changing during the iteration?
You see what I mean? In my opinion you need to project onto a staionary basis, while otherwise you are accumulating vectors (length) with changing direction which is not correct in my opinion.

Indeed. During the iterations, the velocity changes, and if you just one one friction direction, there will be no friction correction orthogonal to the normal and friction direction. This will result in very minor drift. This is noticeable in the Bullet CcdPhysicsDemo with a stack of cylinders (you can disable/comment out the second friction direction to see this).
Quote:
Anyway it is good for the quality of the friction to align the first tangential axes with the relative velocity in contact plane (if it exists).

Indeed, so in Bullet we use the projected normal onto the separating plane as first friction direction, and the orthogonal direction (normal.cross(firstFriction)) as second, unless the projected normal length is zero. In that case you still better use 2 orthogonal friction axis (during the iterations, the velocity changes and makes the projected normal non-zero).

Note that the straighforward accumulated impulse clipping results in 'box' shaped friction model, which in certain cases can cause objects to slide in unnatural directions: this happens when the friction gets clipped onto one of the box axis but not the other. This effect is reduced when using the projected normal as first friction direction.

Also note that we are making a lot of approximations. We can use the current normal impulse for the magnitude for the friction impulse. However, we have the 'chicken-and-egg' problem here, because friction and normal impulse have a two-way effect, ideally we solve a bigger system. However, the relaxation/iterative process usually converges fine and hides all those approximations.

To make a long story shorter: indeed, impulse based friction correction the same as normal collision correction, just using a different axis, and using an impulse magnitude dependent on its related normal collision impulse. That is why the friction has to happen at the end of all iterations. An alternative is to use the approximation of previous frames' normal impulse magnitude. Then the challenge becomes how to come up with an approximation the very first frame (otherwise a single-frame collision would not result in any friction).

To finish with the snippet the Bullet currently uses:
	btVector3 frictionDir1 = vel - cp.m_normalWorldOnB * rel_vel;	btScalar lat_rel_vel = frictionDir1.length2();	if (lat_rel_vel > SIMD_EPSILON)//0.0f)	{		frictionDir1 /= btSqrt(lat_rel_vel);		addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);		btVector3 frictionDir2 = frictionDir1.cross(cp.m_normalWorldOnB);		frictionDir2.normalize();		addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);	} else	{		//re-calculate friction direction every frame, todo: check if this is really needed, because we use persistent contacts that keep the same normal		btVector3	frictionDir1,frictionDir2;		btPlaneSpace1(cp.m_normalWorldOnB,frictionDir1,frictionDir2);		addFrictionConstraint(frictionDir1,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);		addFrictionConstraint(frictionDir2,solverBodyIdA,solverBodyIdB,frictionIndex,cp,rel_pos1,rel_pos2,rb0,rb1, relaxation);	}

Hope this helps,
Erwin
http://bulletphysics.com

This topic is closed to new replies.

Advertisement