What good is Rigid Body Dynamics?
Many of the things we come in contact with on a daily basis are
rigid for all practical purposes. Wrenches, chairs, and sidewalks
all flex and vibrate, but these effects do not strongly impact us.
You ask, "So what?!" My answer is that I would like to endow robots
with the "knowledge" and skill to autonomously manipulate things so that we
can send them into dangerous environments (like the surface of the moon,
at the bottom of an ocean, in a nuclear power plant, or on a battle field)
to build and repair things. The first step in developing robots with
the required knowledge and skill is the development of human scientific
knowledge of how bodies interact with each other on a macro-scale.

The field of rigid body dynamic (more generally, multibody dynamics) is all
about designing mathematical models and algorithms to predict the motions of
bodies and the contact forces that arise between them when. The two most
exciting applications of rigid body dynamics are in robotics and computer games,
but there are also many other important engineering applications, like assembly
and machine design. In robotics, the goal is to build a robot with the
capability to plan and autonomously carry out dexterous manipulation
tasks - like doing the dishes. More and more, computer games contain physics
engines to improve realism - for example, dropping a stone into the gears of
a machine could cause jamming, thus stopping the knife blades from swinging
across your path, and allowing you to escape the collapsing building.

Applications of rigid body dynamics to robotics and computer games differ in
a critical way. In computer games, the physics engine must produce physically
plausible body motions in hard-real-time. This requires approximations of the
"correct" dynamic equations needed in robotics applications. In robotics,
the physics engine is used to answer the question, "Given the input (motor
torques, gravity, etc.), tell me how the robot will move and what contact
forces will arise in the process. This is the forward problem. For example,
if a robot moves its arm in a specified way that causes it to hit a box on a
table, the solution to the forward problem provides a prediction of where
the box will come to rest, whether it will slide or tumble, and the
magnitudes and directions of the forces that arise.

Once it becomes possible to accurately predict the consequences of robot
actions, it also becomes possible to plan the activities of a robot to achieve
a desired goal. Given the current state of the robot and its environment, and
a task specified as a goal state of the robot and its environment (the dishes
are put away and the robot is in the closet charging itself), planning is
equivalent to finding the set of robot actions that transform the robot and
environment to a goal state. This inverse problem is extremely challenging
and often counter-intuitive (two reasons I like studying rigid body dynamics).

Rigid Body Dynamics with Dry Friction:
Basic Mathematical Character
In the early 1980's Lodtstedt published the first paper I know of in which
the equations of motion of a system of rigid bodies in unilateral contact was
formulated as a complementarity problem. After writing the Newton-Euler
equations (a system of ordinary differential equations (ODEs)) for the bodies,
complementarity conditions were introduced to prevent interpenetration (these
constraints were written at the acceleration level, as opposed to the velocity
or position levels). The Newton-Euler equations were solved for the body
accelerations, which were then eliminated from the complementarity conditions,
yielding a complementarity problem whose unknowns are contact forces.
A complementarity condition is a system of inequalities and equations of the
following form:

0 =< w(z) _|_ z >= 0

where w and z are vectors of equal length, "_|_" implies
perpendicularity and inequalities are applied element by element. Roughly
speaking, in rigid body dynamics, the w and z vectors represent
relative accelerations (or velocities) between the bodies at the contacts
and contact forces (or impulses), both expressed in coordinate frames attached
to the contact points (assumed isolated). Letting w_i and z_i
denote the ith elements of the vectors w and z, respectively,
the complementarity constraint enforces the fact that a nonzero contact
force can only exist (z_i > 0) if the contact is being
maintained (w_i = 0), and conversely, a contact force may not
exist (z_i = 0) if the contact is breaking (w_i > 0).
Complementarity formulations of rigid body dynamics generally come in two
flavors: mixed and standard. Mixed formulations have equations and
complementarity conditions. Standard formulations have only
complementarity conditions. Relating back to Lotstedt, he started with
a mixed problem which included the Newton-Euler equations. He obtained
a standard problem by solving the Newton-Euler equations for the body
accelerations and substituted them into the complementarity conditions.

Simulation of dynamic rigid body systems using the complementarity formulation
requires the solution of a differential complementarity problem (DCP), which
can be seen as analogous to an ordinary differential equation (ODE). Time-stepping
methods for DCPs require the solution of one or more complementarity problems
for each time step, just as the numerical solution of ODEs require several
evaluations of the derivative function (a.k.a., the right-hand side). However,
it was noted by Lotstedt, that the acceleration-based complementarity formulation
does not always admit a finite-force solution (The same kind of problem was found
by Painleve in 1895, long before complementarity theory was born). This problem
was solved independently in the mid 1990s by Stewart and Trinkle, and by Anitescu
and Potra, by reformulating the complementarity problem in terms of velocities
and impulses. The main difference between these two methods in their original
form is in the formulations of the non-penetration constraints. In the
Stewart-Trinkle (ST) time-stepper, the constraints are written in terms of positions,
while in the Anitescu-Potra (AP) time-stepper, they are written in term of
velocities. This gave the ST time-stepper built-in constraint stabilization,
which the AP method had to handle outside the time-stepper.

There are many software packages for multibody dynamics simulation. For a list
short list with short descriptions of the basic underlying simulation approach
go to my
sim_packages page.

The animations below (click on the images) were computed using various
time-stepping methods that solve one linear complementarity problem per
time step. The complementarity problems were solved by the PATH
algorithm (Todd Munson, Steve Dirkse, and Michael Ferris).
Here's a video of several spheres interacting in a pyramidal cone with
several fixed sphere poking through. The Coulomb friction coefficient
varies from plane to plane. The red plane is the only frictionless
surface. This animation was produced with Matlab.
Here are several videos of identical "jacks" falling from the sky. The
two images below are the final frames of simulations using
friction coefficients 0.0 and 0.55. Toward the end of the
second video, there were are about 450 potential contacts, and the LCP
being solved for each time step was on the order of 3000 (requiring
about 20 seconds per time step on a 2GHz, Pentium IV laptop running Redhat
Linux 8.0). This simulation was computed with the Umbra software
package built at Sandia National Laboratories by Eric Gottlieb, Fred
Oppel, and Patrick Xavier.
Here's a chain of jacks swinging and hitting a frictional floor. There is
no friction in the joints. This simulation was computed with Umbra.
Here's a vehicle with wheels at the ends of articulations on uneven
terrain. In the first simulation, the articulations are locked and
the vehicle eventually tips over. In the second one, the
articulations are crudely controlled to avoid tipping.
Here is a pawl (part of a MEMS ratchet mechanism) being pushed by two
kinematically controlled fingers across the surface of a fixture plate
into a hole. When it enters the hole, it pushes a cantilever beam out
of the way. The beam is designed to seat the pawl against the three
"fixels" on the "top" side of the hole after the pushing fingers are
removed, but the current implementation failed before that point,
since the motion became kinematically infeasible when the beam
bottomed out. This simulation was computed with Umbra.
Here's a video of a peg, with off-center center of mass, passing through
an orienting device under the influence of gravity. The part must exit
the device with the center of mass down regardless of its entering orientation.
Matlab was used to simulate the dynamics and optimize the device geometry.
Here are some papers related to rigid body dynamics simulation and design problems:
theory,
velocity-base time-stepping,
theory and examples with torsional friction,
the quasistatic problem,
the design and manipulation planning with intermittent contacts.

There are other ways to simulate dynamic rigid body systems. The most
popular are known as penalty methods (used by nearly all software packages
developed for engineering applications (Working Model, Adams, SIMPACK, etc.).
Penalty methods are formulated as ordinary differential equations; the
Newton-Euler equations of the bodies. The contact forces and non-penetration
constraints are handled by virtual springs and dampers. As a pair of bodies
overlaps, the spring force develops to push them apart. While this
approach is fairly easy to implement and can take advantage of well
developed ODE solvers, the down side is that when the interacting bodies
are stiff, the ODE solver must take very small steps.

The other popular way are constraint-based methods. This methods are
essentially the same as complementarity methods in the following sense.
They use the Newton-Euler equations and the same non-penetration
constraints as in the complementarity approach. However, rather than
formulate a DCP, then formulate a DAE (differential equation with
algebraic constraints (see Haug 1982)). The algebraic constraints are
obtained by making assumptions about which contacts will be maintained,
which will begin to slide, which will separate, etc. In other words,
the choices about which complementarity constraints will be active
over the next time step(s) is assumed by logic tests in the solver
software (e.g., contact A will be maintained, because the normal
component of its contact force is large). The main problem with this
method is that it does not leverage complementarity theory or
complementarity problem solvers.