Sign in to follow this  

conservation of momentum

This topic is 811 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

Hey, guys.

I'm reading a book https://www.crcpress.com/Game-Physics-Engine-Development-How-to-Build-a-Robust-Commercial-Grade/Millington/9780123819765 and I have some problems understanding this part

real ParticleContact::calculateSeparatingVelocity() const
{
Vector3 relativeVelocity = particle[0]->getVelocity();
if (particle[1]) relativeVelocity -= particle[1]->getVelocity();
return relativeVelocity * contactNormal;
}

void ParticleContact::resolveVelocity(real duration)
{
// Find the velocity in the direction of the contact.
real separatingVelocity = calculateSeparatingVelocity();

// Check whether it needs to be resolved.
if (separatingVelocity > 0)
{
// The contact is either separating or stationary - there’s
// no impulse required.
return;
}

// Calculate the new separating velocity.
real newSepVelocity = -separatingVelocity * restitution;
real deltaVelocity = newSepVelocity - separatingVelocity;

// We apply the change in velocity to each object in proportion to
// its inverse mass (i.e., those with lower inverse mass [higher
// actual mass] get less change in velocity).
real totalInverseMass = particle[0]->getInverseMass();

if (particle[1]) totalInverseMass += particle[1]->getInverseMass();

// If all particles have infinite mass, then impulses have no effect.
if (totalInverseMass <= 0) return;

// Calculate the impulse to apply.
real impulse = deltaVelocity / totalInverseMass;

// Find the amount of impulse per unit of inverse mass.
Vector3 impulsePerIMass = contactNormal * impulse;

// Apply impulses: they are applied in the direction of the contact,
// and are proportional to the inverse mass.
particle[0]->setVelocity(particle[0]->getVelocity() + impulsePerIMass * particle[0]->getInverseMass());

if (particle[1])
{
// Particle 1 goes in the opposite direction.
particle[1]->setVelocity(particle[1]->getVelocity() + impulsePerIMass * -particle[1]->getInverseMass());
}


I just can't understand how this part correlates with conservation of momentum. Ok, we have momentum of system of two objects. Why if we divide this value on mass of one of the objects we get its change of velocity? Isn't overall momentum is a summ of momentums? Can you please explain what is going on in this piece of code?

Share this post


Link to post
Share on other sites

Well, I'm not going to try to explain every line of the code but will simply concentrate on the physics.  The momentum of a system of two bodies, p,  is simply the sum of their individual momenta, i.e.

 

ptot = m1*v1 + m2*v2

 

Now, if those two bodies interact with each other and exert a force on each other, then according to Newton's 3rd law, the forces are equal and opposite leading to and equal and opposite change in momentum of two objects (and therefore conservation of total momentum), i.e. 

 

dptot = dp1 + dp2 = m1*dv1 + m2*dv2 = 0  (since they cancel each other out)

=> dp = m1*dv1 = -m2*dv2

 

So, the change in momentum, dp, is called the impulse which leads to a change in velocity of

 

dv1 = dp/m1  and  dv2 = -dp/m2

 

i.e. the change in velocity of an object is inversely proportional to its mass.  So, the above code is calculating some impulse that the two objects impart on each other, and then normalises the velocity inversely proportional to its mass to ensure momentum is correctly conserved.  

 

For example, the line 

particle[0]->setVelocity(particle[0]->getVelocity() + impulsePerIMass * particle[0]->getInverseMass());

says that the new velocity is equal to the old velocity plus the velocity change due to the impulse from the collision, i.e. 

 

vnew = vold + dp/m

 

and then almost the same, but with a minus sign for the other object in order to conserve momentum.  Looking at the code, the first lines are doing some kind of geometric or harmonic mean (always confuse those two up) to try and stop objects that have 'infinite mass' (i.e. unmovable objects) from actually moving, but otherwise it's all consistent.   Hope that helps to clear things up.

Share this post


Link to post
Share on other sites

Looks like someone already posted an answer while I was dealing with reactivation issues!  I'll go ahead and post this answer anyways, as I approach the answer from a slightly different angle.

 

You can tell this code respects conservation of momentum because the two objects receive the same impulse (in opposite directions, of course) and impulse is the change in momentum.

This code does two things.  First, it uses the parameters of the collision to determine how to change the speeds of the objects involved.  The output of that is deltaVelocity, the desired change in relative velocity between the two objects.

Second, it calculates how much the velocity of each object must change both to conserve momentum (or, equivalently, to have an equal impulse applied to both objects) and produce the correct change in relative velocity.  The changes of velocity deltaVelocity1 and deltaVelocity2 of the two objects are constrained as below (I treat the velocities as positive numbers in opposite directions):

deltaVelocity1 + deltaVelocity2 = deltaVelocity (correct change in relative velocity)
mass1 * deltaVelocity1 = mass2 * deltaVelocity2 (equal impulses)

The inverse mass of an object with mass m is 1/m, and the code calculates the impulse as deltaVelocity times the sum of the inverse masses (so, 1/mass1 + 1/mass2).  You can use that to calculate deltaVelocity1 and deltaVelocity2, and then plug those in to the above equations to verify that the code is correct.

Share this post


Link to post
Share on other sites

Thanks guys for your answers, but still I can't figure out one thing.

 

Look, the code says:

real impulse = deltaVelocity / totalInverseMass;

I see it like this:

 

We got two bodies with masses m1 and m2 so

M = m1 + m2

Next, we got separating velocity of whole system before impact

V = v1 - v2

And after

V' = v1'-v2'

 

If I understand everything right, in case of this code, total momentum should be

 

P = M(V' - V) = (m1 + m2)((v1'-v2')-(v1-v2))

 

WHOA, what a beast!

 

So, I just can't any connection between this monster and m1(v1'-v1)=-m2(v2'-v2) cuz

 

(m1 + m2)((v1'-v2')-(v1-v2)) just can't be equal to m1(v1'-v1) + m2(v2'-v2) !

 

What am I missing? Where I go wrong?

Share this post


Link to post
Share on other sites

You are applying equal and opposite impules. So when you make the momentum balance for the *whole* system it will not change. The momentum of individual particles can indeed change. But the sum will not.

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

Okay, I think I see where you've gone wrong.  First of all, m1*(v1' - v1) is the change in momentum of 1 and m2*(v2' - v2) is the change for 2.  So the sum, m1*(v1' - v1) + m2*(v2' - v2) is the total change in momentum which is equal to zero!!  This is different from the actual momentum, which can be any value but then doesn't change, so you're comparing two different things there in a way.

 

Second, I think there is some confusion because of their mathematics.  Basically, they are applying a trick to ensure that the total impulse is zero if both object have infinite mass (i.e. are immovable) OR apply all of the impulse to the non-infinite mass object.  So that line of code, i.e. 

real impulse = deltaVelocity / totalInverseMass; 

is just them normalising the impulses to achieve this.  It's a trick and a bit mis-leading.  So, P = M*(V' - V) is not correct because you have to sum up the individual momenta, i.e P = P1 + P2 = m1*v1 + m2*v2, and NOT  P = (m1 + m2)*....   

Share this post


Link to post
Share on other sites

I think I'm starting to get what is going on.

So, is this part of code

real impulse = deltaVelocity / totalInverseMass;

not equal to

real impulse = deltaVelocity * totalActualMass;

?

Geez, I was so confused about this line.

So, does it actually means "how much change of velocity do we get per unit of inverse mass"?

Can I rewrite it to

real impulse = deltaVelocity / totalActualMass;

And then get change of velocities like this

particle[0]->setVelocity(particle[0]->getVelocity() + impulsePerMass * particle[0]->getActualMass());
particle[1]->setVelocity(particle[1]->getVelocity() + impulsePerMass * -particle[1]->getActualMass());

?

I know why we need inverse mass, but is it physically the same?

Also, I'm confused about one more thing: inverse masses addition. In the source code getInverseMass getter has no special code, it just returns value of 1/mass. And there is no special code for inverse mass addition. Can we really add inverse masses as is?

Cuz, lets suppose we have two objects with masses m1=2 and m2=3. Their inverse masses should be mi1=1/2 and mi2=1/3. If we sum them as is we get 5/6. But it feels like it should be 1/(2+3)=1/5. Seems not to be equal to the previous calculation at all. Is it a bug or am I missing something?

Edited by Fennec

Share this post


Link to post
Share on other sites

I'm afraid your intuition is wrong about how to add inverse masses.  You can indeed add them as is; inverse mass is just a number, like speed or regular mass.  If I add 7 and 9, I get 16; likewise, if I add mi1 = 1/2 and mi2 = 1/3 I should get 1/2 + 1/3 = 5/6.  Therefore totalInverseMass = inverseMass1 + inverseMass2 = 1/m1 + 1/m2, which is not the same as 1/totalMass = 1/(m1 + m2).  For that reason, you're right that

real impulse = deltaVelocity / totalInverseMass;

is different from

real impulse = deltaVelocity * totalActualMass;

But you cannot rewrite the line to

real impulse = deltaVelocity / totalActualMass;

because totalActualMass = m1 + m2 is not the same as totalInverseMass = 1/m1 + 1/m2.
 
You can look at totalInverseMass as just a convenient trick to solve the system of equations that I described in my first post:

deltaVelocity1 + deltaVelocity2 = deltaVelocity (correct change in relative velocity)
mass1 * deltaVelocity1 = mass2 * deltaVelocity2 (equal impulses)

 

There is a little bit of physical intuition to be had about totalInverseMass as a quantity by itself, but (in my estimation) it's not too important and there are a few confusing subtleties.  You might be better off considering it as an algebraic step in the solution of the quoted equations.

Share this post


Link to post
Share on other sites

Okay, I've done some calculations, and got some results. So, code says

real impulse = deltaVelocity / totalActualMass;

this is equal to

 

$$impulse=\frac{(\vec{v_{1}} - \vec{v_{2}}) - (\vec{v_{1}'} - \vec{v_{2}'})}{\frac{1}{m_{1}} + \frac{1}{m_{2}}}$$

 

lets try to find WHY it equals to impulse

 

$$\frac{\vec{v_{1}} - \vec{v_{2}}-\vec{v_{1}'}+\vec{v_{2}'}}{\frac{m_{1}+m_{2}}{m_{1}m_{2}}} = (-(\vec{v_{1}'}-\vec{v_{1}}) (\vec{v_{2}'}-\vec{v_{2}}))(\frac{m_{1}m_{2}}{m_{1}+m_{2}}) = (\Delta\vec{v_{2}}-\Delta\vec{v_{1}})(\frac{m_{1}m_{2}}{m_{1}+m_{2}}) = \frac{m_{1}m_{2} (\Delta\vec{v_{2}} - \Delta\vec{v_{1}}) }{m_{1}+m_{2}} = \frac{m_{1}(m_{2}\Delta\vec{v_{2}})-m_{2}(m_{1}\Delta\vec{v_{1}})}{m_{1}+m_{2}}$$

In case we have equation $$m_{1}\Delta\vec{v_{1}}=-m_{2}\Delta\vec{v_{2}}$$ we can do this

$$\frac{m_{1}(m_{2}\Delta\vec{v_{2}})+m_{2}(m_{2}\Delta\vec{v_{2}})}{m_{1}+m_{2}}$$

 

So

 

$$\frac{(m_{1}+m_{2})(m_{2}\Delta\vec{v_{2}})}{m_{1}+m_{2}}$$ = $$m_{2}\Delta\vec{v_{2}}$$

 

OH. SWEET. LORD. It really is.

 

I still don't understand where this equation comes from

 

$$impulse=\frac{(\vec{v_{1}} - \vec{v_{2}}) - (\vec{v_{1}'} - \vec{v_{2}'})}{\frac{1}{m_{1}}+\frac{1}{m_{2}}}$$

 

But at least I figured out how exactly its connected to conservation of momentum.

 

It feels like the formula take its origin from both conservation of momentum and conservation of energy, but I can't see any clues that could point on this.

Edited by Fennec

Share this post


Link to post
Share on other sites

Sorry, what you are writing doesn't make any sense at all:

 

Impulse is defined as P = m * v. But you write it is P = v / m. Then you write P = v / ( 1 / m ). How do expect to make sense of such a mess? The code snippet shows a scalar impulse and your derivation a vector impulse. Finally, if I get you correctly you want to compute a contact impulse, but there is no contact normal. This is all too 'hand-wavy' without any context and obviously you have no idea what you are doing (I don't mean that with any offense!). I recommend to study some basic mechanics and math. Here is some good reference (see the simulation section 'collision response' near the bottom for you specific question). 

http://www.essentialmath.com/tutorial.htm

 

To derive these equations you also need d'Alemberts principle. Funnily there was a very good article here on GameDev.net recently:

http://www.gamedev.net/page/resources/_/technical/math-and-physics/particle-systems-using-constrained-dynamics-r4170

https://en.wikipedia.org/wiki/D%27Alembert%27s_principle

 

If you want to study those things I recommend some other books that might be more helpful to you: 

http://www.amazon.com/Essential-Mathematics-Interactive-Applications-Edition/dp/1482250926/ref=dp_ob_title_bk

http://www.amazon.com/Game-Physics-David-H-Eberly/dp/0123749034/ref=sr_1_2?s=books&ie=UTF8&qid=1443027836&sr=1-2&keywords=game+physics

 

Both the physics and math tutorial have plenty of free material on this topic. You can find the math material in the link above and the physics material here:

http://box2d.org/downloads/

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

To answer your question, one way to derive this equation would be something along the line like this. Given the masses m1, m2, the velocities before the collision v1, v2 and the collision normal n. We are looking for the post-velocities v1', v2' and the collision impulse j. 

 

Since momentum must be conserved (e.g. equal and opposite impulses) and we know that the impulse is acting in the direction of the collision normal n (d'Alembert) we can now relate incoming and outgoing momentum:

 

m1 * v1 + j * n = m1 * v1'  <=>  v1' = v1 + j / m1 * n

m2 * v2  - j * n = m2 * v2'  <=>  v2' = v2  - j / m2 * n

 

These are essentially Newton equations of motion in impulse form. We have two equations for three unknowns and we need one more equation that relates the incoming and outgoing velocity. We can use that both are related by the coefficient of restitution e:

 

e * (v2 - v1) * n =  (v2' - v1') * n

 

If you now substitute these equations you can solve for the impulse j. I leave this an exercise for you. (I marked vector quantities bold as a help and a little tip is that n * n = 1 since the collision normal is a unit vector). Your equations don't show any e so simply use e = 0 which simplifies the last equation to:

 

(v2' - v1') * = 0

 

Regarding your original question we are using two physical principles here. The first is conservation of momentum which means we are applying equal and opposite impulse (e.g. +j * n in equation 1 and -j * n in equation 2). It might be actually easier to think of Newton's 2nd law here (actio = reactio). The other one is d'Alemberts principle which states that constraint forces are not doing any work. Practically this means here that the impulse is acting in the direction of the contact normal e.g: P = j * n. The coefficient of restitution works essentially as a velocity constraint which relates the in and outgoing velocities.

 

Hth,

-Dirk

 

 

PS:

https://en.wikipedia.org/wiki/Coefficient_of_restitution

https://en.wikipedia.org/wiki/D%27Alembert%27s_principle

https://en.wikipedia.org/wiki/Momentum#Conservation

Edited by Dirk Gregorius

Share this post


Link to post
Share on other sites

I think you tried to say

 

(v2' - v1') * n = e * (v2 - v1) * n

 

instead of

 

(v2 - v1) * n = e * (v2' - v1') * n

 

Because velocity after impact is velocity before impact reduced by COF.

 

So, since we have those three equation and  = j * n

v1' = v1 + j / m1 * n

v2' = v2  - j / m2 * n

(v2' - v1') * = 0

 

((v1 + j / m1 * n) (v2 - j / m2 * n)n = 0

(v1 + j / m1 * n v2 + j / m2 * nn = 0

v1+ j / m1 * n * n - v2+ j / m2 *  = 0

 

since unit vector (which is contact normal) multiplied by itself equals to 1

 

v1+ j / m1 - v2+ j / m2  = 0

j / m1 + j / m2 = v2n - v1n

j (1 / m1 + 1 / m2) = v2n - v1n

 

but now I'm not so clear about

m1 * v1 + j * = m1 * v1'

m2 * v2  - j * = m2 * v2'

j  = (v2n - v1n) / (1 / m1 + 1 / m2)

or   = (v- v1) / (1 / m1 + 1 / m2) where v1 and v2 are projections on contact normal.

 

Hooray! It finally makes sense! Thank you so much for guiding me all the way through! Thank you for providing me all these books!

 

Maybe I'm not so clear about d'Alemberts principle but I'll go read about this first.

 

Thank you again! All that information was super helpful!

Share this post


Link to post
Share on other sites

Yes thanks, I fixed this in my post above for completeness. And I am glad you could figure it out!

 

The work principles are little more abstract to grasp, so just give yourself some time. Here is an example I always liked. Imagine you are cleaning your house with a basket of water. When you lift the basket up you do work W = F * h where h is the height you lift it up. If you now move around you are not doing any work anymore since the force is orthogonal to the direction you are moving. 

 

HTH,

-Dirk

Share this post


Link to post
Share on other sites

This topic is 811 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this