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.

`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

`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

2 bodies with contact forces

We focus on the contact distance

The rate of change of the contact distance is the velocity

What we are most interested in is

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

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.

`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

[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*= number of contacts.*d*= acceleration of contact distance of the_{i}''*i*-th contact.**a**= length*n*vector of accelerations of contact distances, so that*a*=_{i}*d*_{i}''**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*a*= the element of_{i j}**A**in the*i*-th row and*j*-th column.

a = A f + b
| (1) |

d = _{i}''a _{i 1}f + _{1}a _{i 2}f +. . . + _{2}a _{i j}f + . . . + _{j}a _{i n}f + _{n}b
_{i} | (2) |

5 bodies, 4 contact points, 8 contact forces

To illustrate how the

d = _{2}''a _{2,1}f + _{1}a _{2,2}f + _{2}a _{2,3}f + _{3}a _{2,4}f + _{4}b
_{2} | (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
*b*._{2} - 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
*b*._{2} - Regarding contact force
*f*_{1}**n**, for a unit change in_{1}*f*we calculate what would be the resulting acceleration of the point on body 2 at contact 2; this is put in_{1}*a*._{2,1} - Regarding contact force
*f*_{2}**n**, for a unit change in_{2}*f*we calculate what would be the resulting acceleration of the corner of body 1 at contact 2; this is put in_{2}*a*._{2,2} - Regarding contact force −
*f*_{2}**n**, for a unit change in_{2}*f*we calculate what would be the resulting acceleration of the point on body 2 at contact 2; this is_{2}*added*to*a*._{2,2} - Regarding contact force −
*f*_{3}**n**, for a unit change in_{3}*f*we calculate what would be the resulting acceleration of the corner of body 1 at contact 2; this is put in_{3}*a*._{2,3} - Because contact forces
*f*_{4}**n**and −_{4}*f*_{4}**n**do not directly affect bodies 1 or 2 we set_{4}*a*= 0._{2,4}

Once we have found the

a >= 0, _{i}f >= 0, _{i}a · f = 0
| (4) |

*a*>= 0 means we disallow interpenetration (which occurs when_{i}*a*< 0). We require all the relative normal accelerations between bodies at contact points to be either zero (in contact) or positive (separating)._{i}*f*>= 0 means we require contact forces to be zero or positive because they can only push, not pull._{i}**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 algorithm begins by setting all the contact forces in

- Begin by considering only contact 1, ignoring all other contact points (which are
likely violating the constraint that
*a*>= 0). We calculate the amount of force_{i}*f*at contact point 1 which is just sufficient to make_{1}*a*= 0._{1} - Next add in contact 2. Solve for
*f*, adjusting_{2}*f*as necessary, so that the three constraints are valid for_{1}*a*,_{1}*a*and_{2}*f*,_{1}*f*._{2} - Continue adding in contacts and contact forces one at a time, while adjusting previous forces as necessary to maintain the three 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.

2 bodies in contact

d = _{i}n ⋅ (_{i}p_{1} − p_{2})
| (5) |

d'' = _{i}n ⋅ (_{i}p_{1}'' − p_{2}'') + 2 n' ⋅ (_{i}p_{1}' − p_{2}')
| (6) |

n ⋅ (_{i}p_{1}'' − p_{2}'').
| (7) |

Let

3 bodies in contact

**x**_{1}= center of mass of body 1**r**_{1}= vector from center of mass of body 1 to**p**_{1}at the*i*-th contact point**r**= vector from center of mass of body 1 to the_{j}*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}}.**v**_{1}=**x**_{1}' = velocity vector of center of mass of body 1*m*_{1}= mass of body 1

p_{1}' = x_{1}' + r_{1}' = v_{1} + ω_{1} × r_{1}
| (7a) |

p_{1}'' = v_{1}' + ω_{1}' × r_{1} + ω_{1} × (ω_{1} × r_{1})
| (8) |

v_{1}' + ω_{1}' × r_{1}
| (8a) |

Because

=

n ⋅ (_{i}n / _{j}m_{1} + (r × _{j}n) × _{j}r_{1} / I_{1})
| (9) |

Keep in mind that we assumed that

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

2 bodies in contact

Recall that

Now we can use that expression for

2 n' ⋅ (_{i}p_{1}' − p_{2}') = 2 (ω_{2} × n) ⋅ (_{i}v_{1} + ω_{1} × r_{1} − (v_{2} + ω_{2} × r_{2}))
| (10) |

Next consider the term

`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
n ⋅ (_{i}p_{1}'' − p_{2}'') = n ⋅ (_{i}v_{1}' + ω_{1}' × r_{1} + ω_{1} × (ω_{1} × r_{1}) &minus
(v_{2}' + ω_{2}' × r_{2} + ω_{2} × (ω_{2} × r_{2})))
| (11) |

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

� Maria Kovalenko, 2007 | PhysLab | Top | Next |