Usually, for every contact point, you solve one non-penetration and two friction constraints. The friction constraints depend on the non-penetration constraint (e.g. Coloumb friction). This is all you need. Other friction models exist, but I am not sure they have been published.

Some of the friction ideas I've seen around are quite simple. At the end of the day you will want to reduce the number of constraints as possible, and go with something straightforward.

It won't be a mountain if you are comfortable formulating and solving constraints. The hard part is, in my opinion, coming up with a simple model and/or figuring out how to reduce inputs to the constraint solver.

I know this isn't very specific info... I haven't done a custom solution myself. I always just did the super cheesy Coloumb approach Dirk described.

The formulation is almost identical to the normal response constraint. Just change normal N to tangent T. The two tangent axis are based on the normal constraint axis. Use any function you can think of that will take a vector as input and generate an orthonormal basis. The two new axes can be treated as your friction tangents.

Erin posted a function a while ago: https://box2d.org/2014/02/computing-a-basis/

That's pretty much it. It works OK. You get sliding artifacts as one axis will converge before the other... You also get 2 extra constraints (one for each tangent axis) per normal contact -- this gets expensive and also deteriorates the solver's convergence. Two extra constraint per normal constraint is a lot. This is why Dirk was mentioning other "unpublished" strategies that are probably nicer because they use less constraints, and avoid the artifact I mentioned.

I have an example of the cheesy Coloumb stuff here: https://github.com/RandyGaul/qu3e/blob/master/src/dynamics/q3ContactSolver.cpp#L152-L175

You can still get problems though. Like a sphere spinning -- there's no way for this tangent basis approach to account for actual spinning. All the Coloumb stuff is doing, is driving a specific contact point's velocity to zero, within the contact plane, and the forces are fed with inputs from the previous tick's normal force. Nowhere in this model is the angular velocity of the *body* taken into account (except when it contributes to the *linear* velocity of a contact point).

The problem with friction is the coupling of the friction forces to the normal (non-penetration) force. E.g. we have

dC1/dt = (v2 + omega2 x r2 - v1 - omega1 x r1 ) * n   and  lambda1 >= 0

dC2/dt = (v2 + omega2 x r2 - v1 - omega1 x r1 ) * s   and  -mu * lambda1 <= lambda2 <= mu * lambda1

dC3/dt = (v2 + omega2 x r2 - v1 - omega1 x r1 ) * t   and  -mu * lambda1 <= lambda3 <= mu * lambda1

The problem we solve is not linear anymore. Practically it is not an LCP in the presence of friction. So all proofs about convergence don't hold anymore. There are a couple of ways to solve this thougj. The first is what Randy suggests and use the normal force from the last frame. This works pretty bad in practice. A better approach is to use the current accumulated normal force (which is changing every iteration). This is mathematically not well defined (at least I haven't seen any proofs here), but it works great in practice.

The next issue with above formulation is that you are clamping against box and not a circle. There are two ways to solve this. The first is to align the first tangent direction with the relative velocity in the friction plane. This is very simple.

v = vn + vt <=> v = (v * n) * n + vt => vt - = v - (v * n) * n

Be careful when the relative velocity is zero. You need to fall back onto a geometric approach as described in Erin's blog post. A better way to solve the friction problem is to use circular clamping. You solve the the friction constraints as a 2x2 block system and then you clamp the length of vector.

Finally what really makes the difference for friction for me is warmstarting. If your warmstarting doesn't work you will get poor friction and poor simulation quality.

HTH,

-Dirk