Problem with friction in a physics modeling computer program

[shadow12345]

I am writing a physics modeling program which simulates rigid bodies in a 3D world. I've come across the problem of properly implementing friction.

In general, the magnitude of the friction force is the friction coefficient times the magnitude of the normal force, and the friction force points in the direction opposing the tangential velocity of the contact point (the point of contact between the colliding bodies. In the computer program, only one contact point is taken into account when in reality when objects collide they do so over an area).

Friction = u * Normal

Now, I know the magnitude of the normal force (the impulse force, in this case). The problem I am running into is when objects come to a complete stop. You see, friction opposes the motion, but it cannot ever *reverse* the direction of motion. Handling friction for spheres that do not rotate, but just slide across a 0 degree incline, is a relatively easy problem. However, I have run into tons of problems for solving for the friction force for spheres that have both linear and angular motion. For right now, all spheres just collide with the earth, which has no motion and 'infinite' mass. I've chosen specific coefficients of friction which correspond to rubber on concrete. Because this is a computer program, static friction doesn't exist (rather, dynamic friction brings objects to rest and keeps them there). Here is the general algorithm for handling friction during an impulse:

-Compute the magnitude of the normal impulse (this works, I've implemented this for spheres that collide with other spheres, and is based on the relative normal velocity and the coefficient of restitution of the collision)

-Apply the tangential impulse force

-Check to see if the tangential velocity has been reversed. If it has, then that means I need to solve for a new friction force whose magnitude is less than the original friction force, but great enough to bring the relative velocity to zero.


edit:
What I've tried to do is resort to using energy methods, as that is the only way to relate linear motion and angular motion. Kinetic energy is in the same units, Joules, regardless of the type of motion it arises from, where as angular momentum is not in the same units as linear momentum. Subsequently, there is no way I can sum the linear and angular momentum in the tangential direction, and then compute a friction change in momentum (an impulse force, because a force is the time derivative of momentum, but because this is a computer simulation we take the limit of this such that the change in time approaches zero, and all forces are seen strictly as changes in momentum).

What I've tried to do is instead compute the amount of kinetic energy which stems from the tangential motion of the object. From that, I can compute the amount of negative work that needs to be done by the friction force to bring the relative tangential velocity to zero. It is negative work because friction works against the motion of the object.

 

[PerennialII]

The problem I am running into is when objects come to a complete stop. You see, friction opposes the motion, but it cannot ever *reverse* the direction of motion. Handling friction for spheres that do not rotate, but just slide across a 0 degree incline, is a relatively easy problem. However, I have run into tons of problems for solving for the friction force for spheres that have both linear and angular motion. For right now, all spheres just collide with the earth, which has no motion and 'infinite' mass.

I can really appreciate your problem ... I've once done a finite element based contact routine implementation with friction in it, and even though the algorithms etc. are relatively simple and straightforward, to get it working between two deformable subspaces in 3D in a numerical analysis is really something. Since what I did was a finite element scheme, I could essentially work with the mesh topology making the contact detection tolerably easy and could directly link it to the analysis incrementation routine ... then it was simply setting up a suitable traction - separation relationship to get the thing going. With respect to the stick formulation I adapted something I've seen quite common done in FEA software, by specifying a maximum allowable elastic slip, which is essentially used a local minimum to compute the stick/slip variables (so it's essentially a nonlocal Coulombian friction model, you take a small enough area from the possibly contacting regions to stabilize the numerical behavior). If we're speaking the same language so to speak the nonlocal approach to the problem was what simplified it in my case and made the whole thing converge way better numerically, however, I'm not completely sure that you should go that way if you're working with strict rigid bodies with well defined geometries, I'd think you ought to be able to handle that by controlling the iterations/convergence/accuracy of your numerical process ?

 

[shadow12345]

I know what you are talking about, as I've had to read the PhD thesis of several people in order to locate a numerical solution to the problem, and all of them refer to computing the 'maximum allowable elastic slip'. To be honest, I didn't understand this method quite well enough to give it a first shot.

And, I am using perfectly rigid bodies with well defined geometries and collision hulls. I am going the easier route with things because this project is just something I'm doing on my own time, and I've never written a physics modeling program before.

What I'm resorting to right now is trying to implement several ad-hoc methods for bringing the relative tangential velocity to zero.

The most obvious one is applying a tangent friction force, and if it reverses the tangential velocity then I reset the object's velocity, minimize the friction force, and try again.

The scenario here is a ball rolling across a flat incline. Ironically enough, this is about the most damn complicated thing to get working properly. First I've got to quantify the normal impulse force, which induces no torque on a sphere during perfectly idealized conditions, then I examine the tangential velocity to determine if friction is applicable. If then try to apply the friction force, which also induces torque, and I apply the adhoc method outlined above.

It nearly works, but somewhere in there I need to add something because the ball won't lose all of its kinetic energy, although it appears to be converging on doing just that.

On top of what I've outlined, I've started coding some 'static friction tracking' routines, which takes assumptions about what you know 'should' be going on during static friction, and tries to apply those to approximate the static friction force. One of these assumptions, for example, is a ball that rolls down an incline. If it doesn't slip, that means that at any instant the linear velocity down the incline plus the velocity from spin at the contact point exactly cancel each other, yielding no relative tangential velocity at the point of contact!

So, that's where I stand right now with this, and I appreciate your response.

 

[PerennialII]

And, I am using perfectly rigid bodies with well defined geometries and collision hulls. I am going the easier route with things because this project is just something I'm doing on my own time, and I've never written a physics modeling program before.

Yeah, I think the well defined (although in many collision cases very complex to track) geometry will in all likelihood give you the tools to get it working, getting the iterative process to work will just take its effort in order to run smoothly (contact analyses quite often seem to need some fine tuning).

The most obvious one is applying a tangent friction force, and if it reverses the tangential velocity then I reset the object's velocity, minimize the friction force, and try again.

The scenario here is a ball rolling across a flat incline. Ironically enough, this is about the most damn complicated thing to get working properly. First I've got to quantify the normal impulse force, which induces no torque on a sphere during perfectly idealized conditions, then I examine the tangential velocity to determine if friction is applicable. If then try to apply the friction force, which also induces torque, and I apply the adhoc method outlined above.

If I'm reading this right that sure ought to work, the rigid body gives you quite a bit of leeway in defining your ad hoc method (energy ones are usually the most stable ones if everything else fails), tracking accuracy over the algorithm & iterations that you ain't missing anything is sure the key thing in getting it to work properly, that everything interacts and connects properly, I can't even count the times when some minus somewhere or neglecting some coupling has caused much headache.

It's a really interesting problem in particularly from the algorithmic perspective ... let me know if I can give some more detailed help/info etc., your description is good but it's easy (especially from this end) to miss some of the "trap holes" of the implementation.

 

[shadow12345]

It's a really interesting problem in particularly from the algorithmic perspective ... let me know if I can give some more detailed help/info etc

I'd like to know more about energy methods. If you have any sources on that please give them to me (online pdfs preferred, but I'm willing to purchase items that have a listing on amazon as long as they run under 100 us dollars)

your description is good but it's easy (especially from this end) to miss some of the "trap holes" of the implementation.
...
that everything interacts and connects properly, I can't even count the times when some minus somewhere or neglecting some coupling has caused much headache.

I'm quite certain that my implementation at this point is correct. Don't get me wrong, I've already located 'trap holes' but I am pretty sure I've fixed them already (it all fits into a single C function that is effectively 50 lines of code). At this point, there's something small that is missing from the algorithm itself to completely bring the sphere to rest.

Once I'm done with this I'm going to work on it using energy methods, and then for different collision hulls.

Thanks again for your time.
edit:
how up north are you? I go to the university of maine in the united states

 

[PerennialII]

I'm quite certain that my implementation at this point is correct. Don't get me wrong, I've already located 'trap holes' but I am pretty sure I've fixed them already (it all fits into a single C function that is effectively 50 lines of code). At this point, there's something small that is missing from the algorithm itself to completely bring the sphere to rest.

Yep, it appears you've done a thorough job in debugging it. A condensed routine.

Once I'm done with this I'm going to work on it using energy methods, and then for different collision hulls.

Thanks again for your time.

Sounds great, I'm actually sometime next year supposed to do an update on one FEA research code with respect to contact. It currently houses simple "hard" contact (penetration is enforced rigorously within a strict iteration tolerance that is usually in micrometer scale) that in deforming contact sometimes produces convergence and accuracy issues, I'll try to overcome it by adding a softer progressive form of contact, a form of damping. Luckily I only need to modify the specific routine ... probably explains some of my interest .

edit:
how up north are you? I go to the university of maine in the united states

Add additional couple of hundred miles northwards and switch continent to the old one ... University of Helsinki.

I'd like to know more about energy methods. If you have any sources on that please give them to me (online pdfs preferred, but I'm willing to purchase items that have a listing on amazon as long as they run under 100 us dollars)

This is something that is quite well represented in the web ... are you thinking about the friction problem, the contact problem or changing the overall kinematical approach to an energy based one ?

 

[rayjohn01]

This is just a thought -- but what actually happens in practice -- it would appear to me that no surface is perfectly smooth and at some point the frictional force may be negated by an elastic collision , if a surface projection is then caught between two such barriers it will the indefinitely oscillate between them with energy being randomly distributed in the lattice . From a simulation viewpoint you may need to introduce a limit somewhere which just says STOP.
Ray.

 

[PerennialII]

This is just a thought -- but what actually happens in practice -- it would appear to me that no surface is perfectly smooth and at some point the frictional force may be negated by an elastic collision , if a surface projection is then caught between two such barriers it will the indefinitely oscillate between them with energy being randomly distributed in the lattice . From a simulation viewpoint you may need to introduce a limit somewhere which just says STOP.
Ray.

You're quite right about it ... many (or more like most) numerical simulation methods require an introduction of "discrete" measure and criteria for convergence/accuracy, which are inherent to the applied numerical process itself. A numerical analysis can be ruined by both requiring too much or too little accuracy (ok, these issues don't usually surface in simple linear problems, but they are a "piece of cake" anyways). So in nowadays (and many more to come) analyses we're forced to for example use non-local implementations in one form or another simply for the results to make any sense on the scale we're working on. Contact analyses where surface topology is introduced are btw a really cool thing in my mind !

 

[rayjohn01]

It reminds me of the problem of a mass sliding down a shaped surface ( say spherical ) numerically you have to decide when the object has reached contact under gravity. You cannot ask if dy=0 , only if dy<0 in which case you then have to reverse the momentum. The implication is that a sliding object is actually bouncing down the surface and is rarely in contact with it . Also the bounces get larger and it appears impossible to predict where it last leaves the surface. The implication for sliding friction I hesistate to imagine .
Ray.

But even that is for a perfect surface , rather than an atomic one with crystal structure -- imperfections and so on. It is hardly surprising that we use lubricants in a fluid or semi fluid state to isolate atomic level roughness and even out the atomic forces --- I quite believe that this area of study could quite consume a life -- i.e get quite 'stuck' in it .
PS how far north is north ?

 

[PerennialII]

It reminds me of the problem of a mass sliding down a shaped surface ( say spherical ) numerically you have to decide when the object has reached contact under gravity. You cannot ask if dy=0 , only if dy<0 in which case you then have to reverse the momentum. The implication is that a sliding object is actually bouncing down the surface and is rarely in contact with it . Also the bounces get larger and it appears impossible to predict where it last leaves the surface. The implication for sliding friction I hesistate to imagine .
Ray.

Yeah, in numerics all kinds of return algorithms are required for all sorts of things ... contact being one of the trickier ones. It's not at all uncommon that a solution deteriorates due to the the very elements you described, in order to attain a convergent solution quite a bit of model fine tuning is required on many parts of it. There are some really good routines out there ... at my institute we do a lot of work on dynamic deforming contact/impact (airplanes have for pretty obvious reasons suddenly arisen in importance in this field ...) and the models typically require some pretty harsh discretizations & simplifications, but in the end perform admirably well for example in multiple surface overlapping contact.

But even that is for a perfect surface , rather than an atomic one with crystal structure -- imperfections and so on. It is hardly surprising that we use lubricants in a fluid or semi fluid state to isolate atomic level roughness and even out the atomic forces --- I quite believe that this area of study could quite consume a life -- i.e get quite 'stuck' in it .
PS how far north is north ?

I've seen some simulation results trying to depict macroscopic, Coulombian, friction from atomic level QM etc. evaluations ... and believe that field is going to be stuck with it in many respects for a long time to come. Taking it even to the crystal level is an interesting prospect, let alone deeper than that.

btw ... currently north could be described as one of the outermost corners of `western´ Europe ... Helsinki.

 

[shadow12345]

EDIT: I've emailed the author and he explained to me what was going on with this :) Just wanted to save you some time in case you were going to respond to this

Hi again, I've perfected my technique for the mentioned scenario I gave above. I downloaded a very advanced physics demo (pretty much exactly what I have been trying to accomplish) and he does something 'funky' for determining the maximum friction that can be applied without reversing the tangential velocity.

I'm assuming that you understand C code. The rest of the function is omitted, as it's just rudimentary normal impulse handling, etc etc. The function's purpose is to process collisions between two objects. Here is a snippet of how the friction force is handled.

		numerator = tangent_speed;	//Not sure what this does 

    Scalar denominator;	//This looks similar to applying a normal impulse
    if (body1)
    {
      denominator = body0->get_inv_mass() + body1->get_inv_mass() + 
        dot(T, cross(body0->get_world_inv_inertia() * (cross(collision.R0, T)), collision.R0)) + 
        dot(T, cross(body1->get_world_inv_inertia() * (cross(collision.R1, T)), collision.R1));
    }
    else
    {
      denominator = body0->get_inv_mass() + 
        dot(T, cross(body0->get_world_inv_inertia() * (cross(collision.R0, T)), collision.R0));
    }

		if (denominator > 0.0f)
		{
			Scalar impulse_to_reverse = numerator / denominator;	//Numerator is the change in normal velocity? 
			Scalar static_friction = collision.static_friction;
      Scalar dynamic_friction = collision.dynamic_friction;

      Scalar impulse_from_normal_impulse = static_friction * normal_impulse;
      Scalar friction_impulse;

      if (impulse_to_reverse < impulse_from_normal_impulse)
        friction_impulse = impulse_to_reverse;
      else
        friction_impulse = dynamic_friction * normal_impulse;

			body0->apply_world_impulse(friction_impulse * T, collision.position);
      if (body1)
  			body1->apply_world_impulse((-friction_impulse) * T, collision.position);

Here is a link to the site with the demo if you want to play with the binary and source (source is C++ to be used with visual studio 6)

physics demo