Rigid Body Simulation Friction

Question

I'm currently working on a 3D Rigid Body simulation program. I have currently managed to get the rigid bodies colliding with the floor and bouncing correctly using the impulse. However, my problem is once they have bounced they accelerate constantly despite using a friction vector to try and slow them.

This is the code when you hit the ground

Rvector fDirection(m_Bodies[i].Vel.x,0.0,m_Bodies[i].Vel.z);
Rvector relativeVelocities = m_Bodies[i].Vel - floorVelocity;
fDirection.normalize();


Real impulse = -(1+e) * (Rvector::dotProduct(relativeVelocities,floorNormal)) 
/ (1/m_Bodies[i].mass + floorMass);
Rvector friction = fDirection*mu*gravity.length()*m_Bodies[i].mass;

Rvector collision_forces = Rvector(0,1,0)*impulse;
collision_forces += friction ;

m_Bodies[i].Vel += (collision_forces/m_Bodies[i].mass);

Thanks

Edit: Here is the integration code.

void RigidBodySimulation::eulerIntegration(float dTime)
{
    Rvector newVel;
    Rvector newPos;
    Rvector zero(0.0, 0.0, 0.0);
    Real one_over_mass;
    Rvector accel;
    for( unsigned int i = 0 ; i < m_Bodies.size(); i++)
    {   
        one_over_mass = 1/m_Bodies[i].mass;
        newVel = m_Bodies[i].Vel + m_Bodies[i].force*one_over_mass*dTime;
        newPos = m_Bodies[i].Pos + m_Bodies[i].Vel*dTime;
        accel = m_Bodies[i].force / m_Bodies[i].mass;
        m_Bodies[i].acceleration = accel;
        m_Bodies[i].newPos = newPos;  
        m_Bodies[i].Vel = newVel;
        m_Bodies[i].Pos = newPos;
    }
}

Answer 1

I have to say, this is a pretty terrible piece of code you have there, and I've been doing this for more than 10 years. You should get a basic textbook on dynamics (like Hibbeler).

Real impulse = -(1+e) * (Rvector::dotProduct(relativeVelocities,floorNormal))
               / (1/m_Bodies[i].mass + floorMass);

This equation sort of looks like you are trying to calculate the restitution impulse from the impact (although the calculation is wrong). First, you must understand that an impulse is not the same thing as a force. An impulse is the integral of the force over a certain interval of time. During an impact, you can assume that period of time to be really small and that's why you perform an instantaneous velocity change. And that's why you should have specified that the integration code has nothing to do with the collision calculation, because it is by-passed for that instant, or at least, it should be if you do an impulse-based calculation. This is what the actual calculation should look like:

Real momentum_before = Rvector::dotProduct(m_Bodies[i].Vel * m_Bodies[i].mass + floorVelocity * floorMass, floorNormal);
Real rel_vel_after = -e * Rvector::dotProduct(relativeVelocities,floorNormal);
// conservation of momentum in normal direction gives this:
Real body_vel_after = (momentum_before + floorMass * rel_vel_after) / (m_Bodies[i].mass + floorMass);
Real floor_vel_after = body_vel_after - rel_vel_after;

Which actually simplifies to one line as follows:

Real body_vel_after = ( (m_Bodies[i].mass - e * floorMass) * Rvector::dotProduct(m_Bodies[i].Vel, floorNormal)
                      + (1.0 + e) * floorMass * Rvector::dotProduct(floorVelocity, floorNormal) 
                      ) / (m_Bodies[i].mass + floorMass);

However, if you assume the floor to have infinite mass (or much larger than that of the body), then you would simply have:

Real body_rel_vel_after = -e * Rvector::dotProduct(relativeVelocities, floorNormal);
Real body_vel_after = Rvector::dotProduct(floorVelocity, floorNormal) + body_rel_vel_after;

It's that simple. But, under that assumption, you do not have conservation of momentum. But in any case, the restitution impulse from the impact can be calculated as:

Real impulse = m_Bodies[i].mass * (body_vel_after - Rvector::dotProduct(m_Bodies[i].Vel, floorNormal));

Now, because the restitution impulse is the integral of the normal force over the small time period of impact, the impulse from the friction during the impact can be calculated from that restitution impact. The friction force is equal to "mu" times the normal force, i.e., |Ff| = mu * |Fn|, this is also valid for the impulse, i.e., |If| = mu * |In|. So, you can compute it directly:

Real friction_impulse = mu * fabs(impulse);

But that's just the magnitude of the friction impulse. It's direction is opposite from the relative tangential velocity, which is:

Rvector tangent_rel_vel = relativeVelocities - Rvector::dotProduct(relativeVelocities, floorNormal) * floorNormal;

And it's direction is:

Rvector dir_rel_vel = tangent_rel_vel;
dir_rel_vel.normalize();

(Notice that I need to keep the tangential velocity intact, because it will be needed later)

At this point, you could compute the tangential velocity after impact as follows (again, under the assumption of an infinite-mass floor, otherwise, it is more complicated than that):

Rvector tangent_rel_vel_after = tangent_rel_vel - dir_rel_vel * friction_impulse / m_Bodies[i].mass;

However, what if the friction impulse causes the tangential relative velocity to get to zero? That's a problem, because, with the above formula, part of the friction impulse could wind up reversing the direction of the tangential relative velocity, which would mean that during the latter part of the impact, the friction force is actually acting in the direction of the velocity (not good). The most friction can do is stop the relative motion. So, you need to check for that condition:

Real tang_rel_vel_change = friction_impulse / mBodies[i].mass;
Rvector tangent_rel_vel_after = tangent_rel_vel - dir_rel_vel * tang_rel_vel_change;

if ( tang_rel_vel_change > tangent_rel_vel.length() ) 
    tangent_rel_vel_after = Rvector(0.0, 0.0, 0.0);   // stop relative motion.

At this point, all you need to do is combine the two final velocities:

m_Bodies[i].Vel = floorVelocity + tangent_rel_vel_after + body_rel_vel_after * floorNormal;

And that's it, at least, for this very simple problem (infinite mass of the floor). In reality, this impulse-based approach become increasingly difficult to deal with as you complicate things: two finite-mass objects, multiple objects, and actual rigid-body dynamics (because you are just doing particle dynamics here). The impulse-based approach is rarely seen anywhere beyond simple schoolyard examples of balls bouncing on the floor. Btw, you shouldn't really call this a "rigid-body" simulator since all you are really doing a particle dynamics (and 3D rigid-body dynamics is way more complicated than this). Also, your integration law is terrible, but that's a whole different story.

Answer 2

Friction is backwards, isn't it?

 Rvector fDirection(m_Bodies[i].Vel.x,0.0,m_Bodies[i].Vel.z);

This is in the direction of the velocity. You then multiply it by some constants then add it to the velocity (along with the impulse)

collision_forces += friction ;
 m_Bodies[i].Vel += (collision_forces/m_Bodies[i].mass);

So that will act to increase the velocity along the floor.

---

There are some other odd things:
your impulse term has: (1/m_Bodies[i].mass + floorMass) This is adding 1/mass to mass. Shouldn't it be (1/(m_Bodies[i].mass + floorMass))

And then in the integrator, you calculate acceleration from force, but then never integrate it, you also apply force to velocity directly. So what's the acceleration member for?