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:

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.

The text Pfeiffer and Glocker, Multibody Dynamics with Unilateral Contacts, Wiley, 1996) contains a complementarity formulation of general rigid body dynamic systems and a section of impressive application results obtained with a straight-forward time-stepping method. You can also visit Pfeiffer's web page here. Brogliato, Non-smooth Mechanics: Models, Dynamics, and Control, 2nded., Springer, 1999, is also worth checking out.

---

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.

---

January 9, 2010