MyPhysicsLab – Rigid Body Collisions

This applet simulates rectangular shaped objects colliding in 2 dimensions. You can set gravity, elasticity (bounciness), and damping (friction) using the slider controls. You can choose from one to six objects. The mass of the objects is initially set to 1.0, but the mass of the green object is adjustable.

Click near an object to exert a rubber band force with your mouse. With the keyboard you can control four "thrusters". The keys S,D,F,E control the blue object. The keys J,K,L,I (and also the arrow keys) control the green object. (If the keys don't work, try clicking near an object first - this ensures that key strokes are passed to the simulation instead of the buttons or scrollbars).

If elasticity is less than 1, and gravity is greater than zero, then the objects eventually settle onto the floor. Because this simulation does not implement resting contact the simulation will get 'stuck'.

Energy Bar Graph

To check the correctness of the simulation, you can click the "show energy" checkbox and look at the display of the pre- and post-collision momentum and energy. If damping = 0 and elasticity = 1, then the energy should not change. The momentum should not change when objects collide, however it will change when an object collides with a wall (walls are assumed to be infinitely massive, so we don't model their movement).

The bar graph shows the gravitational, rotational and translational energy. Note that if gravity is zero, then there is no gravitational potential energy.

See the section on Energy and Momentum for how these quantities are calculated.

Physics Of Motion for Rigid Bodies

Ignoring collisions for the moment, we develop the equations of motion for the objects. The forces on the objects are

For each body, there are three variables to specify position:

In addition, each body has a velocity corresponding to each of these positions:

The equations of motion for the body involve the total force F = (Fx, Fy) and torque τ on the body as follows:

Fx = m x''
Fy = m y''
τ = I θ''

where m = mass and I = moment of inertia about the center of mass (moment of inertia is defined in the next section). Note that we indicate vectors, such as F, with bold and an overline.

Thrust Force

We now build up the equations of motion by adding in one force at a time, beginning with the thrust force. Let T = (Tx, Ty) be the thrust force vector which operates at the point P on the body. The thrust force will accelerate the body according to

m x'' = Tx
m y'' = Ty

In these equations it doesn't matter where on the body the thrust force is applied. The point P can be anywhere on the body, yet because the body is rigid the thrust accelerates the entire body. On the other hand, for rotational movement it matters a great deal where on the body the thrust is applied.

Moment of inertia is the equivalent of mass for rotational physics. It measures how difficult it is to rotate a body about a given point. Since our rectangle bodies rotate freely about their center of mass, we use center of mass as the point for calculating moment of inertia. From a physics textbook you can find the equation for the moment of inertia about the center of a thin rectangular plate. It is given by

I = m (width2 + height2)/12

Let R = (Rx, Ry) = be the distance vector from center of mass to P. The torque at the point P is given by the vector cross product R × T = Rx Ty - Ry Tx so we have

I θ'' = R × T

Actually torque is a vector, but since we are working in 2 dimensions we know that torque is always perpendicular to the plane and so we aren't using vector notation for it. The true vector cross product results in a vector. So the above cross product corresponds to the following 3 dimensional calculation:

(Rx, Ry, 0) × (Tx, Ty, 0) = (0, 0, Rx Ty - Ry Tx)

You can see that the result is always perpendicular to the plane, with zero x and y components. The general vector cross product of two vectors is defined as:

(x, y, z) × (u, v, w) = (y w - z v, -x w + z u, x v - y u)

Rubber Band Force

The rubber band force in the simulation operates along the vector from point P on the body to the mouse position, call this vector L. The rubber band force is B = (Bx, By) = s L where s is the stretchiness constant. The treatment of this force is identical to the thrust force, so we add it in to get the following equations.

m x'' = Tx + Bx
m y'' = Ty + By
I θ'' = R × T + R × B

Gravity Force

Let g = the gravitational constant. Gravity causes a force on the center of mass of -m g in the vertical direction leading to

m y'' = Ty + By - m g

Because gravity works on all points of the body, no torque is generated by gravity.

Damping Force

Damping (friction) causes a force opposite to the direction of motion. The faster you go, the more friction resists your motion. So the magnitude of the damping force is proportional to the velocity. Let k = the proportional damping constant. Adding this to our equations of motion gives

m x'' = Tx + Bx - k x'    (eqn 1)
m y'' = Ty + By - m g - k y'    (eqn 2)
I θ'' = R × T + R × B - k θ'    (eqn 3)

In a more realistic simulation, there may be different damping constants for rotational versus translational motion. But here, we use the same constant for both.

Numerical Simulation

Equations (1-3) are the equations of motion for one of our rectangular bodies. To use the Runge-Kutta method for solving a set of differential equations we need to convert these 3 second order equations into 6 first order equations. Define the velocity variables

Then we have the 6 first order differential equations

x' = vx
y' = vy
θ' = ω
vx' = (Tx + Bx - k vx)/m
vy' = (Ty + By - m g - k vy)/m
ω' = (R × T + R × B - k ω)/I

These equations are now in the form needed for solving numerically with the Runge-Kutta method. Each rectangle body will have its own set of 6 variables (3 positions x,y,θ and 3 velocities vx,vy,ω) and 6 first order differential equations.

How to Calculate Energy and Momentum

This section explains how the energy and momentum of the objects are calculated. See the description of the energy bar graph for how to observe these quantities in the simulation.

If there is no loss of energy to friction (damping = 0) or during collisions (elasticity = 1) then the energy of the system should not change.

A collision between objects should not change the angular momentum. However, a collision with a wall will not preserve angular momentum because the super massive walls are not included in the calculations of the momentum of the objects.

Gravitational energy is given by m g h where h = height of the object's center of mass above the floor.

Translational energy is 1/2 m v2 where v is the velocity vector for the object's center of mass.

Rotational energy is 1/2 I ω2 where I is the moment of inertia and ω is the angular velocity.

Linear momentum in the horizontal direction is given by m vx, in the vertical direction by m vy.

Angular momentum is measured with regard to a particular point in space, for example the origin. It is given by:

I ω k + m r × v

where r is the distance vector from the origin to the object's center of mass. The vector k is the unit z vector, which points out of the x-y plane. You can see that the angular momentum has two components: the spinning component I ω k and the rotation about the origin m r × v.

Collision Detection

At each step in the simulation, we check to see if there is a collision. The bodies can collide with each other or with a boundary wall. For the rectangular shapes we are using it is simple geometry to determine if a collision has occurred by checking if any vertex is within a wall or foreign body. If so, we use a binary search to back up to an earlier time just shortly before the collision occurred. We then make the approximation that the collision takes place at this exact time, and calculate the resulting velocities as described below. The Colliding Blocks simulation further describes these aspects of collision handling.

Physics of Collision for Rigid Bodies in 2 Dimensions

Handling collisions is the most challenging part of this simulation. The explanation here is fairly condensed, so you may want to read some other descriptions as well.

Suppose a vertex on body A is colliding into an edge of body B at the point P. Define the following variables

We now use a standard formula for the velocity of an arbitrary point on a rotating and translating rigid body to get the pre-collision velocities of the points of collision (which is the point P on each body).

vap1 = va1 + ωa1 × rap
vbp1 = vb1 + ωb1 × rbp

Similarly we have the final post-collision velocities vap2 and vbp2 as

vap2 = va2 + ωa2 × rap
vbp2 = vb2 + ωb2 × rbp

Here we are regarding the angular velocity as a 3 dimensional vector perpendicular to the plane, so that the cross product is calculated as

ω × r = (0, 0, ω) × (rx, ry, 0) = (-ω rx, ω ry, 0)

Now we can find an expression for the velocity with which the colliding points are approching each other. We call this the relative velocity. Let vab1 be the initial (pre-collision) relative velocity and vab2 be the final (post-collision) relative velocity. We define the relative velocities as follows

vab1 = vap1 - vbp1
vab2 = vap2 - vbp2

Using the formulas given above for velocity of a point on a rigid body we can expand these to

vab1 = va1 + ωa1 × rap - vb1 - ωb1 × rbp    (eqn 4)
vab2 = va2 + ωa2 × rap - vb2 - ωb2 × rbp    (eqn 5)

Let the vector n be normal (perpendicular) to the edge of body B that is being impacted, and pointing outward from body B. Also let the vector n be of length 1. Then we can find the relative velocity in the direction of the normal n by using the dot product:

relative normal velocity = vab1n = vab1x nx + vab1y ny

Note that for a collision to occur this relative normal velocity must be negative (that is, the objects must be approaching each other). Let e be the elasticity of the collision, having a value between 0 (inelastic) and 1 (perfectly elastic). We now make an important assumption in the form of the following relation

vab2n = -e vab1n    (eqn 6)

This says that the velocity at which the objects fly apart is proportional to the velocity with which they were coming together. The proportionality factor is the elasticity e.

Collision Impulse

To resolve the collision, we will use the concept of an impulse.   An impulse is the change in momentum of an object when a large force is applied over a very brief period of time. We imagine that during the collision there is a very large force acting for a very brief period of time. If you integrate (sum) that force over that brief time, you get the impulse.

Why do we need this strange concept of an impulse? Why not just use the familiar concept of force as in F = m a ?   The answer is that we do not know what the forces are during the collision. With a supercomputer and some very complex software we could model the forces that occur during a collision. We would need to know details about the materials of the bodies, their exact geometry, how they deform under stress, how the stress propagates through the body, etc.

This is far beyond what our simple simulation can do. Luckily, we can assume that the collision happens so quickly that the position and orientation of the bodies do not change during the collision. Instead, all that changes is the velocities of the bodies. Since a change in velocity is a change in momentum (remember momentum = velocity times mass) we have the concept of an impulse.

We are assuming no friction for our collision, so the only force during the collision is in the direction perpendicular to the edge, which is given by the vector n. (Friction would cause a force parallel to the edge as well). Let the net impulse of the collision be j n where j is a parameter to be determined. Body A experiences an impulse of j n while body B experiences the equal but opposite impulse of -j n. The impulse is a change in momentum. Momentum has units of velocity times mass, so if we divide the impulse by the mass we get the change in velocity. We can relate pre and post collision velocities as

va2 = va1 + j n / ma    (eqn 7)
vb2 = vb1 - j n / mb    (eqn 8)

The change in angular momentum of body A from the impulse j n is given by rap × j n. You can think of the impulse as being applied at the point P, so it generates an instantaneous torque there. We then divide by the moment of inertia to convert the change in angular momentum to a change in angular velocity, so the post-collision angular velocities are

ωa2 = ωa1 + (rap × j n) / Ia    (eqn 9)
ωb2 = ωb1 - (rbp × j n) / Ib    (eqn 10)

Solving for the Impulse Parameter

Now we can put all these various equations together and solve for the impulse parameter j. If you aren't interested in the details of solving the equations just skip down to the expression for j.

We start with equation (6), and then expand using our definition of relative velocity in equations (4-5).

vab2n = -e vab1n
(vap2 - vbp2) ∙ n = -e vab1n
[va2 + ωa2 × rap - vb2 - ωb2 × rbp] ∙ n = -e vab1n

Let's expand the left hand side, using the impulse relationships in equations (7-10).

[(va1 + j n / ma) + (ωa1 + (rap × j n) / Ia) × rap - (vb1 - j n / mb) - (ωb1 - (rbp × j n) / Ib) × rbp] ∙ n = -e vab1n

We recognize that the left hand side contains the quantity vab1n so we move that to the right side.

[j n / ma + (rap × j n) × rap / Ia + j n / mb + (rbp × j n) ×rbp / Ib] ∙ n = -(1 + e) vab1n

Note that n is assumed to be normalized so that nn = 1. Also, to simplify the various vector products, we can use the triple scalar product rule

(A × B) ∙ C = (B × C) ∙ A

to derive the following identity

(A × B) × AB = (A × B) ∙ (A × B) = (A × B)2

(Here squaring a vector means taking its dot product with itself.) We can then simplify to

j (1 / ma + 1 / mb + (rap × n)2 / Ia + (rbp × n)2 / Ib)= -(1 + e) vab1n

Dividing leads to our final expression for j.

              -(1 + e) vab1n
j = ————————————————————————————————————————————    (eqn 11)
    1/ma + 1/mb + (rap × n)2/Ia + (rbp × n)2/Ib

We can now calculate j at the time of collision using equation (11), and therefore calculate the post-collision velocities using the earlier equations (7-10).

We can use this same expression for calculating collisions with a wall by assuming that the mass of the wall is infinite. So let mb → ∞ and Ib → ∞ and equation (11) becomes

     -(1 + e) vap1n
j = ————————————————————    (eqn 12)
    1/ma + (rap × n)2/Ia

Multiple Impacts

The above analysis only handles the case of a single corner and edge impact. There are several other simultaneous multiple impact cases, such as:

If there are simultaneous impacts among unrelated objects, these can be easily handled separately. For example, suppose objects A and B collide at the same time that objects C and D collide, then we can deal with each collision separately because they don't affect each other.

For the case of two adjacent corners of object A impacting either a single wall or another object's edge, we change the impact point to be the midpoint between the two corners and then handle it as other collisions.

The more complicated cases with corners of object A impacting different walls or multiple objects is not handled correctly by this simulation. What the code does is handle each impact separately, which will probably be wrong because these collisions are not independent.

Resting Contact

You will notice that certain settings of gravity and elasticity cause the simulation to eventually get stuck. Specifically, this happens when gravity > 0 and elasticity < 1. The problem occurs when an object is settling down onto the floor. The velocity of the object is becoming smaller with each collision because the elasticity is less than one. Eventually the velocity becomes so small that it can't keep gravity from pulling the object down past the floor. This is what results in the simulation getting 'stuck'.

The right thing to do at this point (which this simulation does not do) is recognize that the object is in resting contact with the floor. This would change the equations of motion for the object to include a contact force from the floor, which prevents the penetration of the object. Some of the references cited earlier have sections on calculating resting contact forces.