# myphysicslab – collisions with contact

it may surprise you that getting objects to rest on the floor is a challenging problem! this simulation shows rectangular objects colliding in 2 dimensions. it is based on the rigid body collisions simulation, but it deals with steady contact forces where objects press against the floor, the wall or each other.
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 gray object's thrusters. the keys j,k,l,i -- and also the arrow keys -- control the green object's thrusters. you can also set gravity, elasticity (bounciness), and damping (friction). you can choose from one to six objects. if you don't see the simulation try instructions for enabling java. scroll down to see the math!

try rearranging the blocks so they lean against each other at various angles, or stacking the blocks. you'll see that they are very slippery and hard to get to stack! this is because there is no contact friction in this simulation.

the contact forces are shown as red lines that appear at the corners of the objects. contact forces come in pairs that are equal and opposite, but only one of each pair is shown here.

about friction and damping: in this simulation, the contacts between objects have no friction whatsoever, so they can slide against each other easily. the damping parameter that you can modify (you may have to click on the "show controls" button to see the damping parameter) has to do with the motion of the objects through the space -- damping makes it as though the objects are moving through a thick syrup, or that there is friction between the background surface and the objects.

why is this a challenging problem? contact forces are tricky to figure out. other forces are much simpler because they only depend on the position of one or two objects. for example, the force of a spring between two objects depends only on the position of the end points of the spring. but the contact forces between objects depends on all the other forces on the objects and how they rest against each other in a rather complex way.

## how to calculate contact forces -- overview

the software (a class called `contactsim`) described here for calculating contact forces works closely with software that handles collisions, see rigid body collisions. that software (a superclass called `thruster5`) takes care of collisions. as a result there should be no collisions occuring when the contact force code is active, because the collision code detects collisions, backs up in time to just before the collision and applies an appropriate impulse to cause the objects to bounce (it reverses their velocities). the superclass also calculates the effects of external forces acting on the objects such as gravity, thrust, rubber band, and damping. (interestingly, it is those external forces that cause the contact forces: without them there would be no forces pushing objects into each other.)

the idea is to calculate the exact amount of force needed at each contact to just barely prevent the objects from penetrating. these forces are calculated in the `evaluate` method of the `contactsim` class, which is called repeatedly by the differential equation solver to find the accelerations of the objects. the differential equation solver uses that information to advance the state of the world through time. on entering this `evaluate` method we are given the current positions and velocities of each object, which is required information for figuring out where contact points are and the geometry of contact forces between objects.

first we find all the resting contacts between objects. the criteria are: a corner must be very close to an edge and moving very slowly (in the direction perpendicular to the edge). using these criteria, we assemble a list of contact points. for each contact, we store information such as: which are the two objects that are in contact, what is the normal (perpendicular vector) to the edge, where the point of contact is relative to each object, etc.

2 bodies in contact
consider the two objects in contact in the figure. the contact point on body 1 is the corner labeled p1, the contact point on body 2 is on the edge and labeled p2. let d be the distance between the two points, so that d = |p1p2|. let n be the normal vector pointing out from body 2 which is perpendicular to the edge. in the diagram below we see the two contact forces that develop at the contact point. there is a contact force f n which acts on body 1 at the point of contact to prevent it from penetrating into body 2. by newton's 3rd law there is also an equal but opposite force f n which acts on body 2 at the point of contact to prevent it from penetrating into body 1.

2 bodies with contact forces
what we want to find are the (currently unknown) contact forces, such as f n, which will just barely keep the objects from accelerating into each other, counteracting all the other forces (such as gravity) that are tending to push the objects together.

we focus on the contact distance d, the gap between the points that are in contact. because the objects are in contact, we have that d = 0 (within numerical tolerance). a value of d < 0 would indicate that the objects are interpenetrating. a value of d > 0 would indicate the objects are not in contact but separated by the distance d.

the rate of change of the contact distance is the velocity d'. because the objects are in resting contact we have that d' = 0 (within numerical tolerance). this was one of the criteria for selecting the set of resting contact points. note that the objects can be sliding against each other, because we only care about motion in the direction of the normal (perpendicular to the edge).

what we are most interested in is d'', the acceleration of the gap. without contact forces we will find that d'' < 0 for most of the contact points, meaning that without any contact force, the objects will fall into each other and interpenetrate.

what makes computing the contact forces somewhat tricky is that they are all interrelated. this is because you can have a complex arrangement of objects resting on each other. the solution involves finding a set of equations that captures exactly how each contact force affects each gap. then we solve the equations to find forces that exactly set d'' = 0 at each contact point. actually, we also allow d'' > 0 which means the objects are separating and therefore there is no contact force at that particular contact.

once we have found the exact force at each contact that will preserve the resting contacts (keeping the objects from accelerating into each other) we then apply those forces to each object. this results in changes to the accelerations of the objects. we then return control to the differential equation solver which continues to evolve the simulation over time using those new accelerations. but now the accelerations are modified so that the objects will remain in resting contact as desired.

## calculating contact forces -- in detail

this software (the `contactsim` class) for calculating contact forces is a subclass of the rigid body collisions software (the `thruster5` class). the following description builds on the explanations at that web page, particularly for the external forces. the algorithms used here are based on two papers by david baraff:

[baraff-1] david baraff, an introduction to physically based modeling: rigid body simulation ii - nonpenetration constraints. these are lecture notes from physically based modeling: principles and practice (online siggraph '97 course notes).
note that there is a typographic error on page d62: in several places where you see "a force of j nj(t0)" it should read "a force of fj nj(t0)".

[baraff-2] david baraff. fast contact force computation for nonpenetrating rigid bodies. computer graphics proceedings, annual conference series: 23-34, 1994. 12 pages.

assume that we have a set of n contacts in some particular (but arbitrary) order. define the following variables:
• n = number of contacts.
• di'' = acceleration of contact distance of the i-th contact.
• a = length n vector of accelerations of contact distances, so that ai = di''
• f = length n vector of contact force magnitudes (to be determined)
• b = length n vector giving contact distance accelerations due to external forces (gravity, thrust, etc.)
• a = an n × n matrix that specifies how contact forces affect acceleration of contact distances
• ai j = the element of a in the i-th row and j-th column.
our goal is to set up and solve (subject to certain constraints) the matrix equation
 a = a f + b (1)
the matrix entry ai j tells how much di'' changes from a unit change in fj. if we do the matrix multiplication in equation (1) we see that the acceleration of the i-th contact distance is given by
 di'' = ai 1 f1 + ai 2 f2 +. . . + ai j fj + . . . + ai n fn + bi (2)
this says that the acceleration of the i-th contact distance is a function of the contact forces f1, f2,..., fn and the external forces in bi. we can figure out the ai j and bi (see below). the di'' need to come out zero or positive, because the objects can't accelerate into each other (they are rigid bodies that don't interpenetrate). so the only unknowns are the contact forces fi, and we can solve this system of simultaneous equations to find them.

5 bodies, 4 contact points, 8 contact forces

to illustrate how the a matrix and the b vector are calculated, consider the figure above. we have 5 bodies which are in resting contact at 4 contact points and the resulting contact forces are shown. (note that the normals ni are always perpendicular to an edge and by convention point out from the object; this determines which forces have minus signs). the contact points are numbered in accord with the contact forces, so that contact force f1 n1 is at contact 1, etc. focusing on contact 2 (between bodies 1 and 2), from equation (2) we can write
 d2'' = a2,1 f1 + a2,2 f2 + a2,3 f3 + a2,4 f4 + b2 (3)
consider what is affecting the acceleration of the contact distance d2''. anything that accelerates body 1 or body 2 can affect d2''. here is a list of these sources of acceleration of d2'' and where they wind up in equation (3).
• for all the non-contact forces (gravity, thrust, etc.) currently acting on body 1, we calculate what is the resulting acceleration of the corner of body 1 at contact 2; this is added in to b2.
• similarly, for all the non-contact forces (gravity, thrust, etc.) currently acting on body 2, we calculate what is the resulting acceleration of the point on body 2 at contact 2; this is added in to b2.
• regarding contact force f1 n1, for a unit change in f1 we calculate what would be the resulting acceleration of the point on body 2 at contact 2; this is put in a2,1.
• regarding contact force f2 n2, for a unit change in f2 we calculate what would be the resulting acceleration of the corner of body 1 at contact 2; this is put in a2,2.
• regarding contact force f2 n2, for a unit change in f2 we calculate what would be the resulting acceleration of the point on body 2 at contact 2; this is added to a2,2.
• regarding contact force f3 n3, for a unit change in f3 we calculate what would be the resulting acceleration of the corner of body 1 at contact 2; this is put in a2,3.
• because contact forces f4 n4 and f4 n4 do not directly affect bodies 1 or 2 we set a2,4 = 0.
i hope this gives you some idea of how the a matrix and the b vector are calculated. complete details on these calculations are given below in the sections calculating the a matrix and calculating the b vector.

once we have found the a matrix and b vector, we are ready to solve equation (1) for the unknown force magnitudes in the f vector, subject to the following 3 constraints:
 ai >= 0,   fi >= 0,   a · f = 0 (4)
the constraints in words are:
• ai >= 0 means we disallow interpenetration (which occurs when ai < 0). we require all the relative normal accelerations between bodies at contact points to be either zero (in contact) or positive (separating).
• fi >= 0 means we require contact forces to be zero or positive because they can only push, not pull.
• a · f = 0 means that, at every contact point, if there is a contact force then acceleration is zero or if there is acceleration (separation) then there is no contact force. the dot in the equation is the vector dot product.
the paper fast contact force computation for nonpenetrating rigid bodies. by david baraff goes into full detail about how to solve equation (1) subject to the constraints (4). i'll give a quick description here of how the algorithm described in that paper works.

the algorithm begins by setting all the contact forces in f to zero. the basic idea is to add in one contact force at a time, adding just enough force to maintain the constraints, and readjusting the other forces as necessary. we ignore the contact points we haven't yet considered, gradually increasing the set of contact points that obey the constraints. here is a very simplified sketch of how it works:
• begin by considering only contact 1, ignoring all other contact points (which are likely violating the constraint that ai >= 0). we calculate the amount of force f1 at contact point 1 which is just sufficient to make a1 = 0.
• next add in contact 2. solve for f2, adjusting f1 as necessary, so that the three constraints are valid for a1, a2 and f1, f2.
• continue adding in contacts and contact forces one at a time, while adjusting previous forces as necessary to maintain the three constraints (4).
at each point we are solving a subset of the matrix equation (1) involving only those contact points that we have considered so far. by the final step we have added in all of the contact points and therefore have a complete solution of equation (1) which obeys the constraints (4).

now that we know the contact forces, the last step is to apply them to the bodies, which gives the final set of body accelerations that the `evaluate` method then returns to the differential equation solver. the geometry of where each contact force is applied on each body comes into play here to determine how that contact force affects the linear and angular acceleration of a given body. this process is very similar to that described on the rigid body collisions web page for treatment of forces such as the thrust force.

the differential equation solver then continues to evolve the simulation over time using those new accelerations. these accelerations have however been precisely calculated so that the bodies that are in resting contact remain in resting contact without interpenetrating.

## calculating the a matrix

our goal here is to calculate the entries of the a matrix of equation (1). as described above, the i, j-th entry in the a matrix, ai j, specifies how the j-th contact force, fj, affects the acceleration of the i-th contact distance, di''.

2 bodies in contact
for a particular contact point i, let the bodies be numbered 1 and 2, with body 2 specifying the normal vector ni pointing out from body 2 towards body 1 (as in the diagram above). let p1 be the point on body 1 that is in contact with the point p2 on body 2. let di be the distance between p1 and p2. if we regard p1 and p2 as vectors from the origin to those points on bodies 1 and 2, then their difference (p1p2) is also a vector and we can write
 di = ni ⋅ (p1 − p2) (5)
because the dot product with the normal gives us the component of the vector (p1p2) in the direction of n. note that di is therefore a scalar quantity. keeping in mind that all of the variables in equation (5) are functions of time, we can take the derivative twice as follows: di' = ni' ⋅ (p1p2) + ni ⋅ (p1' − p2')
di'' = ni'' ⋅ (p1p2) + 2 ni' ⋅ (p1' − p2') + ni ⋅ (p1'' − p2'')
because this is a point of contact, we have that p1p2 = 0 so the above simplifies to
 di'' = ni ⋅ (p1'' − p2'') + 2 ni' ⋅ (p1' − p2') (6)
the second term, 2 ni' ⋅ (p1' − p2'), is velocity dependent (so you can immediately calculate it without knowing the forces involved), and is therefore part of bi. see calculating the b vector below. so the fj dependent part of equation (6) is just
 ni ⋅ (p1'' − p2''). (7)
a contact force fj nj only affects di'' if that contact force operates on body 1 or body 2. if that is not the case, we know that ai j = 0. assume now that fj nj affects body 1 or body 2. but keep in mind that fj nj may be a contact force that is not at the i-th contact point (see the earlier diagram of 5 bodies in contact for example). we can find ai j from the geometry of the situation as follows.

let nj be the vector normal at the j-th contact point. assume that the j-th contact involves body 1. then the contact force is fj nj if the contact involves a corner of body 1 or −fj nj if the contact involves an edge of body 1. from the geometry we can calculate the effect of the force on body 1, and therefore on the acceleration of p1''. if the contact involves body 2, then the effect is on p2''. then using equation (7) we will arrive at the contribution of fj nj to ai j.

3 bodies in contact
define the following additional variables:
• x1 = center of mass of body 1
• r1 = vector from center of mass of body 1 to p1 at the i-th contact point
• rj = vector from center of mass of body 1 to the j-th contact point
• ω1 = magnitude of the angular velocity of body 1
• ω1 = angular velocity of body 1, a vector that is perpendicular to the plane, so ω1 = {0, 0, ω1}.
• v1 = x1' = velocity vector of center of mass of body 1
• m1 = mass of body 1
then p1 = x1 + r1 and, recognizing that all these variables are functions of time, we can take derivatives:
 p1' = x1' + r1' = v1 + ω1 × r1 (7a)
here we used the knowledge that r1 is only changing by rotation at a rate of ω1. an elementary bit of calculus then gives the result that r1' = ω1 × r1 where the × indicates vector cross product. taking another derivative: p1'' = v1' + ω1' × r1 + ω1 × r1'
 p1'' = v1' + ω1' × r1 + ω1 × (ω1 × r1) (8)
the term ω1 × (ω1 × r1) is velocity dependent, so it goes into the b vector. the fj dependent part of p1'' is therefore
 v1' + ω1' × r1 (8a)
our task now is to find out how much a unit change in fj affects this.

because v1' is the linear acceleration of body 1, we have by newton's first law v1' = (total force on body 1) /m1 therefore the fj dependent part of v1' is fj nj /m1 because we are working in 2d, we can write angular acceleration as ω1' = τ1 / i1 where τ1 = torque on body 1, and i1 = rotational inertia of body 1 about the center of mass. (if you are working in 3d, angular acceleration has another term; see [baraff-1] for more information). suppose that the j-th contact occurs at the point pj, and the vector rj goes from center of mass of object 1 to pj. then the force fj nj contributes a torque of rj × fj nj. so the fj dependent part of ω1' is (rj × fj nj) / i1 and we can write equation (8a), the fj dependent part of p1'', as fj nj / m1 + (rj × fj nj) × r1 / i1
= fj (nj / m1 + (rj × nj) × r1 / i1)
and the fj dependent part of di'' in equation (6) is then fj ni ⋅ (nj / m1 + (rj × nj) × r1 / i1) therefore, looking back at equation (2) for context, we have the dependence of di'' on fj as
 ni ⋅ (nj / m1 + (rj × nj) × r1 / i1) (9)
which goes into ai j.

keep in mind that we assumed that fj nj affected body 1. if instead it is fj nj affecting body 1, then that changes the sign of expression (9). if fj nj affects body 2 instead of body 1, then the effects are on p2'' in equation (6), so there is a sign change and also we use m2, r2, and i2 in expression (9) and rj is computed relative to body 2. also, in the case where fj nj is at a contact point between bodies 1 and 2, then fj nj affects body 1 and fj nj affects body 2 and you wind up using expression (9) twice and putting the sum in ai j.

## calculating the b vector

here we calculate the b vector of equation (1). the b vector specifies what the contact distance acceleration di'' would be in the absence of contact forces. it depends on non-contact forces such as gravity, thrust, rubber band force, damping, etc., but also has a velocity dependent component.

note that it is these forces that cause the contact forces to exist; the contact forces develop in reaction to the other forces that are pushing the objects together, so contact forces are also known as "reaction forces".

in the context of solving equation (1), the b vector is regarded as a "constant" vector, because at the moment in time that we are solving equation (1) for the contact forces, the b vector is a known quantity.

2 bodies in contact
we continue with the same scenario and variables defined in the previous section calculating the a matrix. that is, we consider the i-th contact, and we number the bodies as 1 and 2, where body 2 determines the normal vector ni of the contact. we are seeking an expression for bi in equation (2). we start from equation (6) which was derived in the previous section: di'' = ni ⋅ (p1'' − p2'') + 2 ni' ⋅ (p1' − p2') we first consider the term 2 ni'⋅ (p1' - p2'). it is dependent only on current velocity, not acceleration, and so is independent of any contact forces being applied, and therefore belongs in the b vector. recall ni is the normal of body 2. therefore, if body 2 is stationary (a wall or floor for example) then its normal does not change, so ni' is the zero vector and this term is zero.

derivation of ni'.
recall that ni is the normal vector to body 2 in the scenario of the previous section. we know that body 2 is rotating with angular velocity ω2. therefore, by an elementary bit of calculus we can write ni' = ω2 × ni where we treat angular velocity as a vector ω2 perpendicular to the plane so that ω2 = {0, 0, ω2} here is the derivation of that result: let ni be defined as follows (because we are in 2d the z component is zero): ni = {nix, niy, 0} we know that the vector ni is rotating at a rate ω2. ignoring any acceleration, we could write the vector ni as a function of time like this: ni = |ni| {cos(ω2 t), sin(ω2 t), 0} where |ni| is the magnitude of the vector ni, and t = time. the first derivative is then ni' = |ni| {−ω2 sin(ω2 t), ω2 cos(ω2 t), 0} which is equivalent to: ni' = {−ω2 niy, ω2 nix, 0} which can also be expressed as a cross product of the two vectors: ni' = ω2 × ni
now we can use that expression for ni' and equation (7a) which gives an expression for p1' (and similarly for p2' as well) to write that
 2 ni' ⋅ (p1' − p2') = 2 (ω2 × ni) ⋅ (v1 + ω1 × r1 − (v2 + ω2 × r2)) (10)
the above expression is part of equation (6), and is added in to bi. but only in the case where body 2 is a moving object; if body 2 is stationary (a wall or floor for example), then ni' is zero.

next consider the term ni ⋅ (p1'' − p2'') of equation (6). we need to include in bi the effect of any non-contact forces on the accelerations of the contact points p1 and p2. we will start from equation (8) which was derived in the previous section: p1'' = v1' + ω1' × r1 + ω1 × (ω1 × r1) it turns out that the accelerations v1' and ω1' that are due to non-contact forces have already been calculated for us by the earlier processes (see `evaluate` method of `thruster5` class). those accelerations are passed in to our `evaluate` method and so we simply plug them into equation (8) to get p1'' (and p2'' has the similar expression). so we can write this out as
 ni ⋅ (p1'' − p2'') = ni ⋅ (v1' + ω1' × r1 + ω1 × (ω1 × r1) &minus (v2' + ω2' × r2 + ω2 × (ω2 × r2))) (11)
which is part of equation (6) to be added in to bi. keep in mind that in the above equation, v1' and ω1' are the accelerations of body 1 due only to non-contact forces (and similarly for v2' and ω2').

we now have in expressions (10) and (11), the contact-force-independent parts of equation (6) which is what goes into bi.