Rigid-Body Dynamics with Friction and Impact ∗

Rigid-body dynamics with unilateral contact is a good approximation for a wide range of everyday phenomena, from the operation of car brakes to walking to rock slides. It is also of vital importance for simulating robots, virtual reality, and realistic animation. However, correctly modeling rigid-body dynamics with friction is difficult due to a number of discontinuities in the behavior of rigid bodies and the discontinuities inherent in the Coulomb friction law. This is particularly crucial for handling situations with large coefficients of friction, which can result in paradoxical results known at least since Painlevé [C. R. Acad. Sci. Paris, 121 (1895), pp. 112–115]. This single example has been a counterexample and cause of controversy ever since, and only recently have there been rigorous mathematical results that show the existence of solutions to his example. The new mathematical developments in rigid-body dynamics have come from several sources: “sweeping processes” and the measure differential inclusions of Moreau in the 1970s and 1980s, the variational inequality approaches of Duvaut and J.-L. Lions in the 1970s, and the use of complementarity problems to formulate frictional contact problems by Lötstedt in the early 1980s. However, it wasn’t until much more recently that these tools were finally able to produce rigorous results about rigid-body dynamics with Coulomb friction and impulses.

1. Rigid Bodies and Friction.Rigid bodies are bodies that cannot deform.They can translate and rotate, but they cannot change their shape.From the outset this must be understood as an approximation to reality, since no bodies are perfectly rigid.However, for a vast number of applications in robotics, manufacturing, biomechanics (such as studying how people walk), and granular materials, this is an excellent approximation.It is also convenient, since it does not require solving large, complex systems of partial differential equations, which is generally difficult to do both analytically and computationally.To see the difference, consider the problem of a bouncing ball.The rigid-body model will assume that the ball does not deform while in flight and that contacts with the ground are instantaneous, at least while the ball is not rolling.On the other hand, a full elastic model will model not only the contacts and the resulting deformation of the entire ball while in contact, but also the elastic oscillations of the ball while it is in flight.Apart from the computational complexity of all this, the analysis of even linearly elastic bodies in contact with a rigid surface subject to Coulomb friction using a Signorini contact condition is not completely developed even now [22,23,56,30,57,62,65].Even if all this can be done, most of the details of the motion for the fully elastic body are not significant on the time-or length-scales of interest in many of the applications described above.For more information about applications of rigid-body dynamics, see, for example, [20,21,92] regarding granular flow and [11,108] regarding virtual reality and computer animation.
There are some disadvantages with a rigid-body model of mechanical systems.The main one is that the velocities must be discontinuous.Consider again a bouncing ball.While the ball is in flight, there are no contact forces acting on it.But when the ball hits the ground, the negative vertical velocity must become a nonnegative vertical velocity instantaneously.The forces must be impulsive; they are no longer ordinary functions of time but rather distributions or measures.While there has been considerable work on differential equations with impulsive right-hand sides, these are usually concerned with situations where the impulsive part is known a priori and is not part of the unknown solution.(The work of Bainov et al., for example, has this character [8,7,71].In these works Bainov et al. can allow for some dependence of the time of the discontinuity on the solution, but the way the solution changes at the discontinuity is assumed to be known, and problems like bouncing balls, where the ball comes to rest in finite time after infinitely many bounces, are beyond the scope of their approach.) The rigid body model with Coulomb friction has been subject to a great deal of controversy, mostly due to a simple model problem of Painlevé which appears not to have solutions.The list of papers on this problem is quite extensive and includes [11,12,27,28,36,49,68,76,77,78,84,82,88,89,91,109,131,128].The modern resolution of Painlevé's problem involves impulsive forces and still generates controversy in some circles.
In this article, an approach is described that combines impulsive forces (measures) with convex analysis.It develops a line of work begun by Schatzman [117] and J. J. Moreau [87,88,90,91] and continued by Monteiro Marques, who produced the first rigorous results in this area [83,84].Related work has been done by Brogliato, which is directed at the control of mechanical systems with friction and impact, and is based on the approach of Moreau and Monteiro Marques; Brogliato's book [14] gives an accessible account of many of these ideas.The intellectual heritage used in this work is extensive: convex analysis, measure theory, complementarity problems, weak* compactness, and convergence are all used in the theory, along with energy dissipation principles and other more traditional tools of applied mathematics.
A number of aspects of rigid-body dynamics nonetheless remain controversial and unresolved.These include the proper formulation of impact laws and how to correctly handle multiple simultaneous contacts.These are discussed below in sections 1.2 and 4.4.Neither of these issues affects the internal consistency of rigid-body models; rather, they deal with how accurately they correspond to experimentally observed behavior.
The structure of this article is as follows.This introduction continues with subsections dealing with Coulomb friction and discontinuous ODEs (section 1.1); impact models (section 1.2); the famous problem of Painlevé (section 1.3); complementarity problems, which are useful tools for formulating problems with discontinuities (section 1.4); and measure differential inclusions (section 1.5).Section 2 discusses how to formulate rigid-body dynamics, first as a continuous problem (section 2.1) and then as a numerical problem (section 2.2), and concludes with a discussion of practicalities and numerical results (section 2.3).Section 3 is about convergence and existence theory for the solutions of rigid-body dynamics problems with impact and friction.Section 4 discusses variants on the ideas presented in the preceding sections.In particular, there are discussions of how to treat rigid-body dynamics as a singular perturbation problem (section 4.1), how to apply symplectic integration methods and difficulties in using them (section 4.2), and how to handle velocity-dependent friction coefficients (section 4.3), multiple contact problems (section 4.4), and the treatment of extended elastic bodies (section 4.5).

Coulomb Friction and Discontinuous Differential
Equations.The Coulomb law is the most common and practical model of friction available.It is, however, a discontinuous law.Coulomb's famous law was derived from a great deal of experimental work that was published in 1785 [26] in his Théorie des machines simples (Theory of simple machines).While Coulomb's law still arouses controversy, and there are many variants on his basic law, it is a suitable starting point and has been successfully used in practice.In its simplest form, Coulomb's law says that the friction force is bounded in magnitude by the normal contact force (N ) times the coefficient of friction (µ); if the contact is sliding, then the magnitude of the friction force is exactly µN in the opposite direction to the relative velocity at the contact.As an example, consider a brick sliding on a ramp, as illustrated in Figure 1.1.
If the brick is sliding down the ramp (v > 0), then since N = mg cos θ, the friction force is F = +µmg cos θ.If the brick is sliding up the ramp (v < 0), then F = −µmg cos θ.The differential equation for the velocity v is m dv dt = mg sin θ − µmg cos θ sgn v, (1.1) where sgn v is +1 if v > 0, −1 if v < 0, and 0 if v = 0.The right-hand side is clearly a discontinuous function of the state variable v.If it were only discontinuous in t, then we could apply Carathéodory's theorem [24, section 2.1] to establish existence of a solution.But Carathéodory's theorem is not applicable.What is worse is that the typical behavior of bricks in this situation is that they stop and stay stopped; that is, we will have v = 0 in finite time, and v will stay at zero for at least an interval of positive length.Understanding the discontinuity is essential for understanding the solution.In fact, solutions do not exist for this differential equation as it is stated, since if v = 0 and dv/dt = 0, then we get 0 = mg sin θ − µmg sgn 0, which can only be true if sin θ = 0. To solve a discontinuous differential equation like this, we need to extend the concept of differential equations to differential inclusions [38,39,40], which were first considered by A. F.  [31,35,60,61,74,95,96,126,127,136,137,138].Of these, only [60,61,126,127] give methods with order higher than one.Excellent overviews of numerical methods for differential inclusions can be found in Dontchev and Lempio [31] or Lempio and Veliov [75].
Since the rigidity of objects is only an approximation, it is reasonable to consider approximating the Coulomb law for the friction force by a continuous or smooth law.The discontinuity in the Coulomb friction law has an important physical consequence: a block on an inclined ramp will not move down the ramp as long as the applied tangential forces do not exceed µN .If the Coulomb law were replaced by a smooth law, then the block would creep down the ramp at a velocity probably proportional to the tangential force divided by µN .Experimentally, very little if any creep is observed in typical situations with dry friction, which demands a friction force function that is discontinuous or very close to being discontinuous.On the other hand, using a continuous approximation for numerical purposes leads to a stiff ODE.Applying implicit time-stepping procedures then results in solutions that are very close to the solution obtained by applying the implicit method to the corresponding differential inclusion.In summary, the physics points to real discontinuities, and there is little advantage numerically in smoothing the discontinuity.The discontinuity is here to stay.
So far we have considered only one-dimensional friction laws where the set of possible friction forces is one-dimensional.For ordinary three-dimensional objects in contact, the plane of relative motion is two-dimensional, and so the set of possible friction forces is two-dimensional.In this case, to allow for complications such as anisotropic friction, we need a better approach.A better basis for formulating physically correct friction models is the maximum dissipation principle.This says that given the normal contact force c n , the friction force c f is the one that maximizes the rate of energy dissipation −c T f v rel , where v rel is the relative velocity at the contact, out of all possible friction forces allowed by the given normal contact force c n .To be more formal, there is a set F C 0 which is the set of possible friction forces c f for c n = 1.The set F C 0 is assumed to be closed, convex, and balanced (F C 0 = −F C 0 ).So the maximum dissipation principle says that In the case of two-dimensional isotropic friction acting on a particle, F C 0 is a disk of radius µ (the coefficient of friction) in the plane of contact.If n is the normal direction to the contact surface, then the total contact force is n c n + c f .The set of possible contact forces is the friction cone, which is given by Of course, as the contact point changes, so does the plane of the possible friction forces.So we must allow F C 0 and F C to depend on the configuration of the system: F C 0 = F C 0 (q) and F C = F C(q).Some pathologies should be prevented, such as having the normal direction n lying in the vector subspace generated by F C 0 .Note that F C 0 and F C are again set-valued functions.However, they are generally continuous ones on the boundary of the admissible region.In the interior of the admissible region, the normal and frictional contact forces must both be zero, so in that case, F C(q) = {0}.

Impact Models.
The behavior of impacting bodies is a topic in rigid-body dynamics that does not arise in formulating other problems in mechanics.It would normally be considered to be the result of the model, rather than an ingredient in building the model.However, for rigid-body dynamics, an impact is regarded as an atomic (i.e., indivisible) event and must be modeled as such.The use of measure differential inclusions below does not save us from having to decide.And it is a fact of nature that some materials behave more elastically than others on impact.Some seem to be inelastic with little or no "bounce," and some are very elastic and seem to lose very little energy in an impact.As a general rule, modelers use a coefficient of restitution, usually denoted here by ε between zero and one to describe the impact behavior of a pair of bodies or materials.For ε = 0 we have purely inelastic impacts, and for ε = 1 the impacts are purely elastic.
There are two generally pervasive approaches to modeling impact behavior.One (the Newtonian approach) relates pre-and postimpact velocities' normal components (typically [94]; the other (the Poisson approach) divides the impact into compression and decompression phases and relates the impulse in the decompression phase to the impulse in the compression phase: N decompr = −ε N compr [110,114].The value of the contact impulse for the compression phase of the contact should be determined by the impulse needed for inelastic impact.Each approach has been found to produce an increase in the total mechanical energy in certain situations!For difficulties with the Newtonian approach in its naive form (even with only one contact in two dimensions), see Stronge's article [135] and the references therein and also Keller's short article [63].
Whichever approach is used, the case of inelastic impacts is an important reference case that both approaches must handle: ε = 0.For one contact, this means that n T v(t + ) = 0 for any time t when there is contact.For a fuller discussion of how the Newton and Poisson approaches can be used and modeled, the reader is again referred to the excellent book of Brogliato [14].Another discussion of the Newton and Poisson formulations is given in Chatterjee and Ruina [18], who also discuss alternative collision laws.Note that both the Newton and Poisson impact laws have serious defects, which are discussed in [18,135], for example.
One of the abiding difficulties in this area is the lack of understanding of the physical mechanisms behind impact processes.It has long been believed that three mechanisms are responsible for energy dissipation in impact: (1) localized plastic deformation; (2) viscous damping in the material; and (3) energy transfer to elastic vibrations.Until recently it has not been clear which is the most important.As a result, models have lagged in terms of their physical accuracy and sophistication.Recent work has thrown fresh light onto this issue.
Recent experimental and simulation work, particularly by Stoianovici and Hurmuzlu [133], has highlighted the importance of elastic vibrations, although localized plastic deformation appears to play a role.Viscous damping seems to play very little direct role in the impact behavior.This is mostly because the time-scale of the impact (typically between 10 µs and 10 ms for metal objects of sizes between 1 cm and 1 m) is too short for significant viscous damping, except for very high frequency modes.
What Stoianovici and Hurmuzlu did was to drop steel bars onto a hard, massive block, record the impacts using high-speed video recorders, and identify contacts by measuring the current flowing from the bars to the block underneath.One of the main results of their experimental investigation was the way the kinematic coefficient of restitution ε = −(n T v + )/(n T v − ) depended on the angle of the bar relative to the upper surface of the block (φ). Figure 1.2 is taken from Stoianovici and Hurmuzlu [133] with the kind permission of the authors and the ASME Journal of Applied Mechanics.As can be seen in Figure 1.2, if the bar is dropped vertically (φ = 90 • ), then ε ranges between 0.8 and 0.9.As the angle is decreased, ε decreases until ε is between 0.1 and 0.2, at an angle that appears to approach 90 • as the bar becomes more slender.As φ is further reduced, increases in a sometimes erratic way until it reaches a value of around 0.6 for φ = 0 • .The results of Stoianovici and Hurmuzlu [133] also include simulation results, which show excellent agreement with the experimental results.The simulation results do not include plastic deformation, but viscous damping is included at the contact point.The viscous damping parameter was calibrated to fit the experimental results.
If we assume that plastic deformation and viscous damping are insignificant during impact, then an effective coefficient of restitution can be computed for a given body using only linear elasticity and taking the material stiffness (i.e., Young's modulus) to infinity to recover a rigid limit.This gives a coefficient of restitution ε = ε(q).However, to make ε(q) well defined, we need to assume that prior to contact there are no excited modes of elastic vibration.If the body undergoes a rapid series of impacts, or if other damping mechanisms are too slow, then later impacts will have excited modes of vibration, which will invalidate the assumptions.Significant energy could be transferred from the elastic modes back into rigid-body modes.While common experience suggests that these effects are probably not large, they may still be important for accurate simulations.Appropriate modeling of vibrational effects in a rigid-body framework has not yet been attempted.An open question here is, "In the rigid limit, do the phases (rather than just the amplitudes) of elastic vibrations play a significant role in the impact behavior of the body?" If the answer is "yes," then practical prediction will be extremely difficult, because the elastic vibrations have frequencies that are typically in the acoustic range of 100 Hz to 1 kHz, and future impact times will have to be computed to within a fraction of the period of these vibrations (i.e., probably within 100 µs) in order to accurately simulate the behavior.In robotics, where speeds often range up to 1 ms −1 , computing impact times to this accuracy would require knowledge of distances to within 0.1 mm, which is a very high accuracy requirement.If only the amplitudes of the elastic vibrations are important, then a plausible approach to modeling is to consider the motion of the body as consisting of rigid motions plus small high-frequency elastic vibrations.The elastic vibrations would generally decay while the body is in free flight but could change rapidly (instantaneously in the rigid limit!) during impact.This approach may allow the development of new low-order models that can accurately capture the motion of common objects.
In the remainder of the paper we will consider only the Newton and Poisson types of impact models.

Painlev é's Problem.
In 1895 P. Painlevé published a paper [97] in which he presented a simple rigid-body dynamics problem which appears not to have a solution.The scenario is much like getting chalk to screech on a blackboard.It is illustrated in Figure 1.3.
The parameter values that we are most interested in are large µ and small J/ml 2 .We let (x, y) be the coordinates of the center of mass; θ is the angle of the rod with respect to the horizontal (counterclockwise); the coordinates of the point where the rod contacts the table are (x c , y c ).If ẋ < 0 and θ = 0, then ẋc < 0 and the relative sliding velocity (of the rod with respect to the table) is to the left.Therefore, the friction force F should have magnitude +µN to the right.The equations of motion for (x, y, θ) are Of particular importance is the equation for y c , since the rigid-body condition is that y c ≥ 0. If y c = 0 and ẏc = 0, then to prevent penetration we need ÿc ≥ 0. This can be computed in terms of F and N since y c = y − (l/2) cos θ.Differentiating twice gives ÿc = ÿ + (l/2) sin θ θ + (l/2) cos θ θ2 .
We require that ÿc ≥ 0 and N ≥ 0 (since there is no adhesion to the surface-that is, there is no glue).If ÿc > 0, then contact is broken immediately and so (provided there are no impulsive forces) we must take N and also F to be zero.Conversely, if N > 0, then contact must be maintained and so ÿc = 0.Both cases are covered by the simple complementarity condition (Note that "a ⊥ b" indicates that ab = 0 if a and b are scalars and that a T b = 0 if they are vectors; for a vector "a ≥ 0" means that the components a i of a are nonnegative.) Now consider the situation with large coefficient of friction µ and small J/ml 2 .For µ large, the contact force (the vector sum of F and N ) is applied to the contact point at a shallow angle.This means that the torque from the contact forces is counterclockwise.If J/ml 2 is sufficiently small, the effect of this torque will overpower the effect of the vertical contact force N in preventing penetration.When this happens, the torque will cause the rod to rotate into the table.Mathematically, the coefficient (1/m) − (l 2 /4J) cos θ(µ sin θ − cos θ) < 0 in this situation.If (l/2) cos θ θ2 − g < 0 as well, then no value of N ≥ 0 will give ÿc ≥ 0, and penetration appears impossible to prevent.
This line of analysis has been used to argue that rigid-body dynamics and the Coulomb law of friction are inconsistent.For a recent example of this line of reasoning, see the otherwise excellent textbook of Pfeiffer and Glocker [109, section 5.3, pp.[61][62][63].However, this line of reasoning is flawed.The flaw is the assumption that all forces are bounded functions of time.In particular, it rules out the possibility that the horizontal component of the velocity could be brought to zero instantaneously by impulsive contact forces.If ẋc was driven to zero during the impact, then we would no longer require that F = +µN , the angle of the contact forces does not need to be shallow, and solutions can exist.This example also shows that it is important to allow for impulsive forces in rigid-body dynamics, even when there are no collisions.It also shows that care should be taken in developing formulations of rigid-body dynamics that can directly incorporate impulsive forces and discontinuous solutions.
Note that when Pfeiffer and Glocker [109, section 8.2] discuss impact laws with friction, they implicitly allow for this resolution of Painlevé's problem.However, this would require explicitly invoking the impact law in a situation without impacts.
On a final note, the Painlevé problem also has a bearing on formulating impact laws since applying the Newton impact law would imply that ẏc immediately after impact should be zero, even for the "perfectly elastic" case, ε = 1.However, simulations and asymptotic analysis with the table represented by a spring of stiffness k show that as k → ∞ there is a definite nonzero limit of ẏc after impact.This corresponds to ε = +∞ using Newton's impact law, since ẏc just before impact is zero.A Poisson impact law would be more realistic, since it would give a positive value for ẏc just after the impact for any ε > 0, due to the impulsive contact forces.

Complementarity Problems.
Complementarity problems give a uniform way of representing a very large range of conditions that would require modeling using discontinuous functions.For example, the contact conditions (1.6) of the previous section avoid the need for discontinuous or infinite functions that would arise if we tried to represent N as N (y c ), for example.
Complementarity problems typically have the following form: Find z ∈ R n such that 0 ≤ z ⊥ f (z) ≥ 0, where f : R n → R n is a given continuous function.Note that this is equivalent to requiring that z i ≥ 0 and f i (z) ≥ 0 for all i = 1, . . ., n, and that z i = 0 or f i (z) = 0 for all i = 1, . . ., n. Linear complementarity problems (LCPs) are complementarity problems where f is affine: f (z) = Mz + q, where M is a given matrix and q is a given vector.Not all (not even "almost all") linear complementarity problems have solutions.It should be understood that LCPs are truly nonlinear problems, in spite of their linear-combinatorial structure.For a comprehensive overview of LCPs, see Cottle, Pang, and Stone [25].A more recent overview of applications of linear and nonlinear complementarity problems is the article by Ferris and Pang [37].
Complementarity problems arise in many contexts; one of the most important of these is constrained optimization.If we consider an optimization problem then provided a suitable constraint qualification [41] holds, there are Lagrange multipliers λ i such that hold at the minimizing x.These are the Kuhn-Tucker (or Karush-Kuhn-Tucker) conditions for optimality [41,43].Note that the last two lines of (1.8) are actually complementarity conditions: 0 ≤ λ ⊥ g(x) ≥ 0. The existence of solutions and algorithms to compute them is known for a wide class of LCPs.The best-known algorithms for LCPs are pivoting methods, which are closely related to the simplex method for linear programming.The best known of these is Lemke's method [25, section 4.4].The most useful result for our purpose shows that Lemke's method computes a solution to the LCP 0 ≤ z ⊥ Mz + q ≥ 0, where M is a copositive matrix (defined below) and q T y ≥ 0 whenever y solves the homogeneous LCP 0 ≤ y ⊥ My ≥ 0 [25,Thms. 3.8.6 and 4.4.12].A matrix M is copositive if y ≥ 0 implies that y T My ≥ 0. The application of complementarity problems to general unilateral contact problems seems to have begun with the work of Lötstedt [76,77,78,79].In this work, Lötstedt was able to show the existence of solutions of certain LCPs provided the coefficient of friction was sufficiently small.Since then, there has been much work applying complementarity problems to contact problems [4,5,66,67,103,102,131,139,101], and the question of the existence of solutions to the complementarity problems has attracted a great deal of attention.
Complementarity problems can be extended to a more general setting which can include infinite-dimensional problems [55,99]: if X is a Banach space and K is a closed convex cone in X, then the dual cone is One way in which we will use this idea is to take X to be the set of continuous functions on [0, T ] so that X * is the set of finite-valued Borel measures on [0, T ].The cone K ⊂ C[0, T ] is the cone of nonnegative continuous functions, so K * is the cone of finite-valued nonnegative Borel measures.We understand complementarity between a measure µ and a function φ on [0, T ] to mean that µ(E) ≥ 0 for all Borel E, φ(t) ≥ 0 for all t, and µ, φ = φ(t) µ(dt) = 0.This we denote by the shorthand 0 ≤ µ ⊥ φ ≥ 0.

Measure Differential Inclusions.
General rigid-body dynamics requires a new approach which can combine impulsive forces with differential inclusions.A framework for this approach can be built using measure differential inclusions [87,90,91].For Moreau, this work developed out of a study of sweeping processes [16,69,70,85,86]: in a sweeping process, there is a set-valued function C(t) which is convex for all t and a "particle" at x(t) which is "swept" along by the C(t) so that x(t) ∈ C(t) for all t.If x(t) is in the interior of C(t), then x (t) will be zero, and if x(t) is on the boundary of C(t), then x (t) will be in the normal cone of C(t) at x(t).If C(t) is allowed to be discontinuous as a set-valued function, then x(t) may also have to "jump": x(t + ) − x(t − ) must then belong to the normal cone of C(t) at x(t + ).Measure differential inclusions (although not always called that) can be found in other contexts in the work of Schatzman [117,118,119], for example.
In a measure differential inclusion, x, v), (1.9) F (t, x) does not need to be bounded, and v(.) is only required to have bounded variation.The other properties assumed about F for differential inclusions should hold: F (t, x) should be a closed convex set, and the graph of F should be closed.The difficulty is in interpreting "dv/dt ∈ F (t, x)" when v(.) is not absolutely continuous or is discontinuous.Then "dv/dt" is not an integrable function, or perhaps not even a function at all, but a distribution or measure.We do, however, suppose that v has bounded variation on a finite interval [0, T ].That is, the supremum of We write µ = Dv to denote the distributional derivative of v(.).The expression "dv/dt" can be thought of as a Radon-Nikodym derivative of Dv with respect to the ordinary Lebesgue measure Dt.However, this will work only if Dv is an absolutely continuous measure with respect to Dt, which amounts to requiring that v(.) is an absolutely continuous function.If v(.) is not an absolutely continuous function we need to extend the notion of the Radon-Nikodym derivative.
Fortunately it is possible to split measures into absolutely continuous and singular parts: the Lebesgue decomposition of Dv is Dv = µ s + a Dt, where a(.) is a Lebesgue integrable function and µ s is a Borel measure that is singular with respect to the Lebesgue measure [64, pp. 111-113].The singular part µ s contains all the forces and impulses "at infinity"; this singular part is supported on a set that has Lebesgue measure zero.So we should require that a(t) ∈ F (t, x(t)) for Lebesgue almost all t.For µ s we want to look at the part of F (t, x) "at infinity."For closed convex F (t, x) we can use the asymptotic or regression cone [53].This asymptotic cone of a closed convex set K is the set of directions in K "at infinity": While we can't take the Radon-Nikodym derivative of µ s with respect to the Lebesgue measure, we can "normalize" µ s by using its variation |µ s | which is a nonnegative Borel measure: any finite-valued Borel measure µ has variation |µ|, which is another finite-valued measure defined by where {E i } ranges over all countable Borel partitions of a Borel set E. Then clearly µ(E) ≤ |µ|(E) for all Borel sets E, and µ is an absolutely continuous measure with respect to |µ|.So we define "dv/dt ∈ F (t, x)" to mean that a(t) ∈ F (t, x(t)), Dt a.e., (1.11) An alternative definition that is useful for handling problems of convergence is the following: for each continuous function φ ≥ 0, not everywhere zero, The definition based on (1.11), (1.12)I call the strong definition, and the definition based on (1.13)I call the weak definition.Under conditions (C1)-(C3) below, the two definitions are equivalent [130].
• (C1) The graph of F is closed, and F (t, x) is a closed convex set for all (t, x).
The strong definition is useful for proving properties about the solutions, while the weak definition is useful for obtaining existence results.Consider, for example, a sequence of functions {v h } h>0 generated by some numerical procedure.If it can be proved that v h ≤ M for some constant M and v h (0) = v 0 for all h > 0, then by the Helly selection theorem [93], there is a pointwise convergent subsequence v h k (.) which converges to a function v(.) with v ≤ M .In this subsequence, the measures Dv h k D v weak*.If v h satisfies a measure differential inclusion dv h /dt ∈ F h (t, x(t)), and 0 < h < h implies that F h (t, x) ⊂ F h (t, x), then using the weak formulation it is easy to show that d v/dt ∈ F h (t, x(t)) for all h > 0, and consequently While measure differential inclusions cannot satisfactorily treat all aspects of rigid-body dynamics, they give structure to a large part of it.The additional conditions can be specified in terms of complementarity conditions, usually between measures and functions.

Formulation and Simulation.
The paradoxes uncovered by Painlevé show the need for care in formulating the equations of rigid-body dynamics with contacts.When it comes to formulating numerical methods for simulating rigid-body systems, we must be even more careful.Current methods for handling rigid-body dynamics fall into several categories: 1.For each possible configuration of contacts, solve for the forces that could be generated and feed the result into an ODE or differential algebraic equation (DAE) solver.This method is vulnerable to Painlevé's problem as well as being very cumbersome.2. Formulate the problem as in case 1 but use a complementarity problem to decide at each step which contacts are active.(This is basically the approach of Pfeiffer and Glocker [109].)This is still vulnerable to Painlevé's problem.3. Use a penalty formulation of the no-interpenetration condition.This corresponds to approximating the rigid bodies by very stiff bodies.This avoids the theoretical existence questions and avoids impulses but raises new questions about singular perturbations and the accuracy of the computational methods.
A penalty approach can also be used to approximate the Coulomb friction law.All of these modifications give very stiff differential equations.4. Use a time-stepping formulation based on integrals of the contact forces rather than their instantaneous values.Complementarity conditions or optimization conditions are used to resolve whether contact is maintained or broken.The approach I will present here is approach 4: use a time-stepping formulation.This approach directly incorporates impulsive forces.In fact, in a single simulation with fixed step-size, there is no way of determining if there actually is an impulsive force, since only integrals are represented or computed.

The Continuous Problem.
In many ways it is easier to write down a numerical method for rigid-body dynamics than it is to say exactly what the method is trying to compute.Nevertheless, with the tools of measure differential inclusions, functions of bounded variation, and complementarity problems, this can be fairly straightforward for many situations.To keep matters simple at this stage, we consider a one contact problem.
To obtain the equations of motion we first need the equations of motion for a system without contacts.To do this we use generalized coordinates q, which can contain rectilinear coordinates of centers of mass as well as angles and orientation parameters and other types of coordinates.Associated with this are the generalized velocities v = dq/dt and the Lagrangian L(q, v) = 1 2 v T M (q)v−V (q), where 1 2 v T M (q)v is the kinetic energy (M (q) is the mass matrix) and V (q) is the potential energy.If there are only internal forces, the equations of motion are given by Lagrange's equations of motion: where With external and contact forces we can write this as The admissible region of configuration space is given by { x | f (q) ≥ 0 } for a suitable function f .The normal direction vector at q is n(q) = ∇f (q).The contact conditions that we need are essentially the Signorini contact conditions: 0 ≤ f (q) ⊥ c n ≥ 0. We represent the set of possible friction forces through The maximal dissipation principle states that we choose β so as to maximize the dissipation rate v T D(x)β over all β ∈ c n F C 0 (x).Because of discontinuities, we need to be careful to make some distinctions between v , and v(t).For inelastic impacts, for example, we require that n T v + (t) = 0 for all t where there is contact: f (q(t)) = 0. Also, in the maximum dissipation principle we should minimize (v + ) T D(q)β over β ∈ c n F C 0 (q), that is, min This is a convex program (with linear objective function), although it is nonsmooth.If c n > 0, then using a Slater condition, there exists a Lagrange multiplier λ ≥ 0 such that ∂ β h(β, λ) = 0 where h(β, λ) = v T D(q)β − λ(µc n − ψ(β)) and ∂f (z) is the generalized gradient of f at z [19, Thms.6.1.1 and 6 If c n = 0, then β = 0 by the condition µc n − ψ(β) ≥ 0; since ∂ψ(0) contains a neighborhood of the origin, there is a value of λ ≥ 0 such that 0 ∈ D(q) T v + + λ ∂ψ(0), no matter what v + and D(q) are.
For a simple particle, we can take ψ(β) = β 2 = β T β.Then as long as the columns of D(q) span the friction plane, we have a representation of the isotropic Coulomb law of friction.

Formulations and Function Spaces.
The Signorini contact condition that 0 ≤ f (q(t)) ⊥ c n ≥ 0 involves complementarity between a continuous function t → f (q(t)) and a measure c n , so the complementarity condition can be represented by f (q(t)) c n (dt) = 0 as well as the usual inequality conditions: f (q(t)) ≥ 0 for all t and c n (E) ≥ 0 for all Borel sets E ⊂ [0, T ].The integrated maximum dissipation principle min β (v + ) T D(q)β subject to ψ(β) ≤ µc n can be formulated where β is a measure.The main difficulty is that ψ(β) must be interpreted as a measure.Convex functions of measures were studied as measures in [29].Since ψ(.) is a positively homogeneous function, ψ(β) = ψ(dβ/d|β|) |β|, and the condition ψ(β) ≤ µc n is equivalent to ψ(dβ/d|β|) |β| ≤ µc n or even ψ(dβ/dc n ) ≤ µ a.e.Note that the derivatives dβ/d|β| and dβ/dc n are all Radon-Nikodym derivatives.The function λ is a bounded Borel function; we can bound λ ∞ by max t µ D(q(t)

Complete Formulations.
The continuous problem with one contact and inelastic impacts can be formulated as follows.The data of the problem consists of the mass matrix M (q), the contact constraint f (q), its gradient n(q) = ∇f (q), the matrix of direction vectors defining plane of the friction cone D(q), the coefficient of friction µ, and the function ψ(β) so that the friction cone is F C(q) = {n(q)c n +D(q)β | ψ(β) ≤ c n }.Then we wish to find the trajectory q(.) which should be absolutely continuous, the velocity function v(.) which should be of bounded variation, along with the measures c n for the normal contact forces and β to describe the frictional forces, and a bounded Borel-measureable function λ where For partly or fully elastic impacts with coefficient of restitution 0 ≤ ε ≤ 1, we can replace (2.9) with n(q(t)) T  Related models and theory for partly elastic impacts, at least for the frictionless case, can be found in the work of Mabrouk [80,81].Note that this is a Newton-style model of partly elastic impact.Interestingly, this formulation of the Newton approach always dissipates energy, unlike the formulation of the Newton approach discussed by Stronge in [135].A Poisson approach for partly elastic impacts is developed in Anitescu and Potra [5].A new Poisson-type formulation of impact recently developed by Pang and Tzitzouris [104] is based on complementarity problems and is guaranteed to be dissipative.
It has been argued that complementarity conditions between the normal contact force and the normal velocity of the form of (2.6) are not always valid [17].The phenomena discussed in Chatterjee [17] are associated with multiple simultaneous contacts and are discussed in section 4.4 of this article.

Numerical Methods.
Numerical methods for rigid-body dynamics without friction or impact have been studied extensively in the engineering and also mathematics literature.See, for example, [34,143,124,122,52,2,123,111,33,107,46,42,142,13], in reverse chronological order.One of the reasons for this is the need to simulate mechanical systems, especially the complex mechanical systems that arise in manufacturing processes and robotics.Even simple grasping problems involve problems of contact, impact, and friction.So far, most of the numerical methods are based on ODEs or DAEs or both.To avoid the difficulties with the discontinuity in Coulomb's law, for instance, a regularized version is used, as shown in Figure 2.1.In this paper, a different strategy is recommended.

Time-Stepping Methods for Problems with Impulses.
In order to handle problems like those arising with the Painlevé problem, a time-stepping approach which uses the integrals of the force functions (or measures) c n and β over each timestep interval [t l , t l+1 ] is used.Two different numerical formulations are presented here.The first is based on linear complementarity problems and uses a polyhedral approximation F C(q) to the friction cone F C(q).The second is a nonlinear complementarity formulation which uses ψ(β) directly.
The polyhedral approximation to the friction cone is the cone generated by { n(q) + µ d i (q) | i = 1, 2, . . ., m }, where µ d i (q) is a collection of direction vectors in F C 0 (q).Write D(q) = [d 1 (q), d 2 (q), . . ., d m (q)].The friction forces D(q)β are approximated by D(q) β, where β i ≥ 0 and i β i ≤ µ c n .The relationship between F C(q) and F C(q) is illustrated in Figure 2.2.It is assumed that for each i there is a j, where d i (q) = −d j (q).This is related to the assumption that F C 0 (q) is a balanced set: The discretization of (2.4)-(2.9) is the problem of finding q l+1 and v l+1 (and the force integrals c l+1 n , β l+1 and Lagrange multiplier λ l+1 ) given q l and v l for a time-step of size h > 0 that satisfy the following conditions: where f (q l + h v l ) < 0; if f (q l + h v l ) ≥ 0, then we set c l+1 n = 0 and β l+1 = 0 and solve the first two equations.Note that e is a vector of ones of the appropriate size.
This discretization is a partly implicit Euler method.Therefore it can give only O(h) accuracy at best.However, unlike conventional discretizations, it can handle impulsive forces, in particular Painlevé's problem.Note that the complementarity condition 0 ≤ f (q) ⊥ c n ≥ 0 does not appear explicitly in (2.11)-(2.15);(2.13) is essentially the differentiated form of this condition.Using the differentiated constraint only can result in the true constraint "drifting" into the inadmissible region, which is an effect that has been noticed in relation to DAE formulations of rigid-body dynamics with bilateral (i.e., equality) constraints [15].It is tempting to replace (2.13) with 0 ≤ f (q l+1 ) ⊥ c l+1 n ≥ 0. This does not work: the resulting discretization behaves as if it had a "random" coefficient of restitution when impacts occur.(The effective coefficient of restitution depends on the time within the time-step [t l , t l+1 ] that contact occurs.)Since (2.13) uses differentiated constraints, it may be occasionally advisable to project q l+1 back to the feasible region.This can be done without disturbing the time-stepping, since the time-stepping method is a one-step method.
If we ignore the dependence of M on q and substitute for v l+1 and q l+1 in terms of c l+1 n , β l+1 , and λ l+1 , then (2.11)-(2.15)with ε = 0 can be reduced to a pure LCP which has the form (superscripts suppressed) The 3 × 3 block matrix in this LCP is a copositive matrix; a solution to this LCP exists and can be constructed by Lemke's algorithm [25,Cor. 4.4.12].
Replacing (2.14)-(2.15)with (2.17)-(2.18),using the nonlinear condition µc n ≥ ψ(β) instead of the conditions µc n ≥ e T β and β ≥ 0, leads to highly nonlinear complementarity problems.Since ψ is nonsmooth, it may be desirable to replace it with smooth functions.For example, for isotropic friction, where ψ(β) = β 2 , we can replace the condition µc n ≥ ψ(β) with (µc n ) 2 ≥ β 2 2 = β T β.The difficulty in using this condition is that if c n = 0, the constraint qualification fails, and Lagrange multipliers may not exist for the maximal dissipation principle in this form.However, it can be restored by using a version of the Fritz John conditions (see [43,   These nonlinear complementarity problems can be solved using nonsmooth Newton methods, which converge at least locally (see, for example, [51,98,100,112]).These methods can be made global by using continuation methods (see, for example, [1,140] for overviews of practical homotopy/continuation methods and [101] for a theoretical analysis of the basis for using continuation methods for frictional contact problems).Some numerical results are presented graphically in Figure 2.3.This shows one ball being thrown and colliding with the first of three balls on a table.The coefficient  of restitution is 0.9, which is nearly perfectly elastic.Note that most of the momentum is transferred to the third ball on the table.Since the impact is slightly oblique, the balls do not remain in a line but move off in different directions.Also note some other effects: the thrown ball rebounds with a strong spin, and when it bounces, the spin reverses direction, resulting in an alternating pattern in terms of direction and speed of its bounces.This effect can be clearly seen when "superballs" (small solid rubber balls designed to bounce well, sold in toy stores) are made to spin rapidly.More subtly, the second of the balls on the table bounces up slightly on impact.An analysis of the frictional and normal impulses will reveal that in a similar situation with n colinear balls struck by a rolling ball, all of the even-numbered originally stationary balls will receive an upward frictional impulse.
Another numerical example is the Painlevé example in an impact situation which is illustrated in Figure 2.4.Any method which assumes that the horizontal component of the contact velocity does not change during the moment of impact cannot compute the solution of this problem correctly.The numerical results shown are for a rod that is initially spinning counterclockwise with stationary center of mass.It then falls while spinning until it hits the table.The solution shown is for inelastic impacts.The numbers shown are the times for each configuration in seconds.

Convergence and Existence.
The fundamental question is, "Do the numerically computed trajectories satisfying (2.11)-(2.15)converge to exact solutions of (2.4)-(2.9)?"If this question can be answered affirmatively, then it also gives an existence theorem for solutions to rigid-body dynamics.
Certain situations are harder than others to deal with; we have already seen how Painlevé's famous problem causes difficulties for analysis.The difficulties that arise with Painlevé's problem can be related to the fact that M (q) −1 F C(q) points outside the set of admissible velocities.This can be avoided if we assume Erdmann's condition that for any q ∈ ∂C and z ∈ F C(q), n(q) T M (q) −1 z ≥ 0 [36].Note that Erdmann's condition holds automatically for frictionless problems, because z ∈ F C(q) = cone(n(q)) implies that z = αn(q) (α ≥ 0) and n(q) T M (q) −1 z = α n(q) T M (q) −1 n(q) ≥ 0 since M (q) is positive definite.Also, Erdmann's condition implies that there cannot be any impulses without collisions; that is, n(q(t)) The question of convergence is answered affirmatively in [129] for one inelastic contact (ε = 0) in the case of one-dimensional friction (that is, F C 0 (q) is a onedimensional set), or if Erdmann's condition holds [36].This result includes Painlevé's problem.It also extends the fundamental results of Monteiro Marques [84] for inelastic frictional dynamics of a particle.
In this section we will review how to prove convergence (and thus existence) of solutions to rigid-body dynamics problems.This is based on the results in [129], which should be consulted for more details.A more complete summary of the proof is in [128].

The Easy Part.
The main part of the theorem can be handled by standard techniques once a few basic facts are established.The first is the existence of solutions to the mixed complementarity conditions (2.11)-(2.15) at each time-step.This is done in [129] by a combination of an approximation argument and the Brouwer fixed-point theorem.The numerical solutions are approximately dissipative.This is based on the result that for constant M , n, D, and linear V (q), the numerical solutions are exactly dissipative, which is proved using ) and substituting for M (v l+1 − v l ) in terms of c l+1 n and β l+1 .Using a generalization of the discrete Gronwall lemma to nonlinear ODEs, the energy of the computed trajectory is shown to be bounded on some sufficiently small interval [t 0 , t 1 ], t 1 > t 0 .This means that the computed velocities v h (.) are bounded on [t 0 , t 1 ].Since dq/dt = v, the numerically computed functions q h (.) are uniformly Lipschitz and thus equicontinuous.By the Arzelá-Ascoli theorem, there is a uniformly convergent subsequence.We restrict attention to the subsequence.
Assuming that the friction cones F C(q) are pointed and that q → F C(q) has a closed graph, it can be concluded from the uniform boundedness of v h (.) that the variations v h are also uniformly bounded.By the Helly selection theorem, there is a convergent subsequence where v h (.) → v(.) pointwise, and the differential measures dv h dv weak*.By the weak* closedness results in [130], the limit v(.) satisfies the measure differential inclusion in the continuous formulation.

The Hard Part.
The hard part involves nonstandard arguments and complementarity conditions between functions that converge pointwise and measures that converge weak*.
First, the limits q(.) and v(.) can be proved to be exactly dissipative in the sense that the total energy of the solution 1  2 v(t) T M (q(t))v(t) + V (q(t)) is a nonincreasing function of t.This could not be proved before, because the argument used at this stage needs a uniform bound on T 0 v h .From the exact dissipativity result, we can go on to show that the solution grows at most exponentially on any interval where the limit exists.By "bootstrapping" this argument with the local boundedness result, it can be shown that limits q(.) and v(.) exist on [0, ∞).The next step is to show that f (q(t)) ≥ 0 for all t.This is needed because the timestepping formulation only ensures that the differentiated constraint n(q l ) T v l+1 ≥ 0 at each step where f (q l + hv l ) < 0. The other step that is needed is to show that the support of the limiting measure c n is contained in the set { t | f (q(t)) = 0 }.This is equivalent to showing that c n ({ t | f (q(t)) > 0 }) = 0.This can be done by showing that c n ({ t | f (q(t)) ≥ }) = 0 for any > 0. Thus we have the complementarity condition that 0 ≤ f (q) ⊥ c n ≥ 0 between a (continuous) function and a Borel measure.This is the correct limiting contact condition.
The next condition to consider is the inelastic impact rule.The proof that the limits q(.) and v(.) satisfy the condition f (q(t)) = 0 ⇒ n(q(t)) T v + (t) = 0 in [129] is restricted to the one-contact case.Simple examples show that this condition cannot be directly generalized to multiple contacts.The proof is based on the result that for the one-contact case, n(q l+1 ) T v l+1 ≤ max(0, (n(q l ) T v l )) + O(h), which is obtained from the complementarity condition 0 ≤ n(q l ) T v l+1 ⊥ c l+1 n ≥ 0 for f (q l + hv l ) < 0. (If f (q l + hv l ) ≥ 0, then there are no contact forces, and the result holds trivially.) Finally, the Coulomb law in its maximal dissipation version must be satisfied.That is, we need to show that as measures.Since the numerical velocity functions v h converge pointwise, and the measures β h and c h n only converge weak*, we cannot directly take limits, even though this property holds for the discrete formulation.First, it is shown that n T v − µ D(q) T v ∞ is a right lower semicontinuous function of time.This is obtained from the discrete formulation, where (v l+1 − v l ) T (n(q l )c l+1 n + D(q l )β l+1 ) is expanded two ways: one by expanding (v l+1 − v l ) using the discrete equations of motion, and the other by multiplying the expression out and using the complementarity and optimization conditions.The next step is to show that any weak* limit ν of any subsequence of c h n n(q h ) T v h − µ D(q h ) T v h ∞ is bounded below by the measure c n n(q) T v + − µ D(q) T v + ∞ .The next step is to use this inequality as part of a detailed accounting of the energy in the limit.Part of these computations involve forming the differential measure of the energy: The first term comes from a result of Moreau for differential measures of quadratic functions [85].Corresponding formulas for the numerical approximations are developed.One of the difficulties in this step is attempting to take limits of the measures ((v h ) + − (v h ) − ) T M (q) dv h .Unfortunately, because the numerical acceleration measures dv h only converge weak* and not weakly, the usual weak lower semicontinuity arguments (based on Mazur's lemma) do not hold, and other means are needed to prove the desired inequalities.At this point, the argument breaks up into three cases.The first is where the limit v(t) is continuous; note that the measure (v + − v − )M (q)dv has its support on the discontinuities in v.This enables the proof of the validity of the Coulomb law for the limit at every point of continuity of v.Of course, the real interest is in the discontinuities of v.The second case is where the condition of Erdmann [36] holds: the friction cone transformed by M (q) −1 is strictly inside the tangent cone of the feasible region.That is, n(q) T M (q) −1 z > 0 for any z ∈ F C(q) and z = 0.The violation of this condition seems to be essential for Painlevé's and related examples.The argument used will not be described here, except that it uses the inelastic impact result.The final case includes the Painlevé case: it is the case where the friction plane is one-dimensional (dim span D(q) = 1).A general argument shows that the only way the Coulomb law can fail for the limit is if the numerical friction impulses have oscillations that are not O(h) in magnitude within a small time interval.In the one-dimensional case, the only way the numerical friction impulses can oscillate in this way is if the sign of (D(q) T v) 1 oscillates while c l+1 n is "large"; this is shown to be impossible.
While these results do not cover all cases of interest (especially two-dimensional friction planes), they do provide a satisfying resolution of some well-known "paradoxes."Further, to extend these results to two-dimensional friction planes, all that is needed is to rule out some pathological behavior by the numerical methods.

Uniqueness.
It is well known that solutions to rigid-body dynamics problems can have multiple solutions.For example, consider Painlevé's problem with slightly different starting values: θ < − 2g/l cos θ ensures that ẋc < 0, but the complementarity problem obtained using the analysis of section 1.3 is Since the quantity in brackets is negative and (l/2) cos θ θ2 − g > 0, there are two solutions to this complementarity problem, one of which corresponds to continued sliding, the other of which corresponds to the contact breaking (ÿ c > 0).There is, in fact, a third impulsive solution, which cannot be obtained from this analysis, just as the impulsive resolution of Painlevé's original problem cannot be obtained through his analysis.It turns out that the "continued sliding" solution is extremely unstable.
To see this, a singular perturbation analysis of the rigid limit of stiff contact must be carried out; as the stiffness k of the contact is increased, the time for a perturbation to double in magnitude decreases as O(1/ √ k).In general, nonuniqueness of solutions can be seen as the result of extreme instability in stiff approximations.To predict which of the possible solutions actually occurs requires knowledge of the microscopic details of the contact, which is not available at this level of modeling.In practice, there is little point in trying to predict which solution occurs in reality, because other unmodeled aspects of a real system will outweigh the effects of imprecise modeling of the contact phenomena.
The lack of uniqueness also signals another difficulty: the solutions are not necessarily continuous in the data.The best that can be done in these circumstances is to repeat the simulations, but with random disturbances, so that a set of possible outcomes can be determined.Determining all possible outcomes computationally is not feasible at present.Even if all of the time-stepping problems (2.11)-(2.15)had unique solutions, this does not imply that the continuous problem has unique solutions.(An example related to this is Ballard's existence and nonuniqueness proof for quasi-static contact problems; the solutions are nonunique for arbitrarily small friction coefficients µ > 0 [9]!)

Control and Optimization.
Much of the work done on rigid-body dynamics has a view to controlling mechanical systems; this is a clear motivation in robotics, for example.See [14].Classical optimal control methods such as the Pontryagin principle and the calculus of variations require a certain amount of smoothness, even for the nonsmooth versions in Clarke [19], which do not hold in general for these problems.Even if the normal contact force N is known, at least as a function of the configuration, and the map from control functions to solution trajectories is well defined and Lipschitz, the variational approach to optimal control is beyond current tools.Partial results in this direction have been achieved by Frankowska [44] based on the differential inclusion formulation, for example; however, these results assume that the right-hand side F (x) in dx/dt ∈ F (x) is a Lipschitz set-valued function in the Hausdorff metric.Note that the Hausdorff metric for closed bounded sets in R n is given by where d(x, C) is the distance from x to C. By contrast, the right-hand side for the Coulomb friction law is not even continuous in the Hausdorff metric on sets.
Since this area is of practical importance, it will see a great deal of interest.Several avenues are open: One is to regularize the Coulomb law and the contact conditions, and apply Pontryagin to the regularized system.This leads to stiff equations and a singular perturbation approach.Another is to attempt to handle discontinuous differential equations directly; since the computation of the adjoint or dual variables in the Pontryagin approach amounts to a differentiation, the adjoint functions can be expected to be discontinuous for discontinuous ODEs, even though the trajectories are continuous.Another approach is to apply a pattern search to these problems in order to compute optimal trajectories.None of these approaches is perfect, and much needs to be done.

4.
Open Questions, Other Ideas.Many people have worked on rigid-body dynamics from different points of view, and there are many aspects of these problems that deserve attention.Some of this related work is briefly discussed in this section.This is not an exhaustive discussion of these issues, but rather an introduction to some of the interesting unresolved issues in this area.
From an applications point of view, rigid-body dynamics is just the limit of dynamics of elastic bodies as the coefficients of elasticity go to infinity.Some work has been done on justifying impact laws from this approach.The "holy grail" for this point of view would be complete justification of models of rigid-body dynamics by singular perturbation analysis.This still seems a long way off.
Integrators for mechanical systems with bilateral constraints have become more refined with the development of DAE solvers and symplectic integrators.Rigid bodies can be approximated by very stiff elastic bodies and symplectic integrators can be applied to these problems.By using adaptive step-sizes (which requires some care to keep symplecticness!)solutions can be generated, at least for perfectly elastic impacts.However, I believe that for impact and contact problems, a tighter coupling between the contact conditions and the integrators should be developed.
The standard Coulomb model is inadequate in itself to describe a number of experimentally observed phenomena, such as velocity-dependent coefficients of friction.This brings new difficulties to the theory.The simplest model is a two-coefficients model with a static and a dynamic coefficient of friction.The extra discontinuity makes the theoretical analysis particularly difficult.Another model with a better experimental basis is to have µ = µ( v ), which is also better behaved theoretically, although there are a number of open questions.
Multiple contact problems predominate in applications; handling them well requires attention to theory, too.
Finally, dynamic problems with elastic bodies in contact with Coulomb friction are still beyond the reach of our theoretical tools, in spite of the considerable progress and partial or approximate solutions that have been found.4.1.Singular Perturbations and the Rigid Limit.Rigid-body dynamics is really the study of the limiting case, where the elasticity constants, such as Young's modulus, go to infinity.Ultimately, the justification of rigid-body dynamics should be via a singular perturbation theory for stiff elastic bodies.Such a theory has not yet been developed, although there are some partial results in this direction.
Paoli and Schatzman [105,106] have proven some singular perturbation results for problems involving particles without friction, but with partly elastic impacts.In [105,106] the limiting problem is assumed to have a convex admissible set K with nonempty interior and a C 2 boundary.The equations of motion in the interior of K are given by u = f (t, u, u ), while on the boundary the impact law is u (t + ) = −εu n (t − )+ u t (t − ), where u n (t) = n (n T u (t)) is the normal component of the velocity u (t) and u t (t) = u (t) − u n (t) is the tangential component of the velocity.To approximate this, Paoli and Schatzman use a penalty law where P K (z) is the projection of z onto K, G(u, v) = (u T v)u/ u 2 if u = 0 and zero if u = 0, and η is related to the coefficient of restitution ε.Note that the quantity u λ − P K (u λ ) acts as a measure of the amount of interpenetration occurring between the particle and the boundary of the admissible region.Paoli and Schatzman show that as λ ↓ 0, u λ → u with u a solution of the continuous problem with convergence strong in W 1,p for 1 ≤ p < ∞ and weak in W 1,∞ .Clearly there is a great deal to be done in this area: friction needs to be included, the convexity assumption on the admissible region should be weakened, and the whole theory should be extended to infinite-dimensional problems involving extended elastic bodies.Some relevant issues are discussed below in sections 4.4 (on multiple contacts) and 4.5 (on elastic bodies).

Symplectic and Variational
Methods for Impact.Symplectic integrators have become very popular for integrating the equations of motion of conservative mechanical systems.(See, for example, [45,48,73,113,116].)Symplectic methods are methods for Hamiltonian systems: given a Hamiltonian H(q, p) defined in terms of the coordinates q and momentum variables p, the Hamiltonian system is Symplectic integrators preserve the differential forms dq i dp i .They do not exactly conserve the Hamiltonian (which is the total energy), but they nearly conserve a "numerical" Hamiltonian H h (q, p), where h > 0 is the step-size used for the symplectic method [116].The numerical Hamiltonian H h (q, p), is usually obtained from H via Taylor series.
Adaptive versions of these methods can deal with certain kinds of singularities, such as those arising from terms of the form f (q) −α for large α which models "soft" contacts, and bilaterally constrained mechanical systems.Because of the possibility of unbounded right-hand sides, implicit methods are necessary.Take, for example, the simplest implicit (nonpartitioned) symplectic method, the implicit midpoint rule for dx/dt = f (t, x): Consider the case of a single particle in one dimension with only contact forces and gravitation, and a simple inequality constraint: q(t) ≥ 0. The Hamiltonian of the system is H(q, p) = 1 2m p 2 + ψ(q) + gq, where ψ(q) = 0 for q ≥ 0 and ψ(q) = +∞ for q < 0. This gives the differential equations and inclusions Note that ∂ψ(q) = {0} if q > 0, ∂ψ(q) = R + if q = 0, and ∂ψ(q) = ∅ if q < 0.
Applying the implicit midpoint rule to this system gives the numerical scheme Writing N n+1 for the element we choose from −(h/2)∂ψ((q n +q n+1 )/2) (that is, N n+1 is the normal contact impulse), we get the mixed complementarity problem Since symplectic methods do not conserve H(q, p) and the numerical Hamiltonian H h (q, p) is obtained assuming a smooth H(q, p), it does not follow that this method is conservative even in the limit as h ↓ 0. In fact, numerical results indicate that it is not conservative: if q n < 0 and p n < 0, then N n+1 > 0 is needed to ensure that q n +q n+1 ≥ 0; by complementarity, q n+1 = −q n > 0; therefore p n + p n+1 = 2(2m/h)q n+1 and p n+1 = −p n + (4mq n+1 /h) > |p n |.The effective coefficient of restitution will depend on q n /(hp n ).A typical numerical trajectory is shown in Figure 4.1(a) for q(0) = 1, p(0) = 0 with h = 0.01, g = 9.8213, and m = 1.Clearly energy is not conserved, even approximately.This should not be surprising: the Hamiltonian cannot in general be conserved by any fixed step-size unpartitioned method [47], and conservation of momentum is of little use when we have a priori unknown impulsive contact forces!If we replace the complementarity condition 0 ≤ N n+1 ⊥ q n + q n+1 ≥ 0 with one involving the momenta 0 ≤ N n+1 ⊥ p n + p n+1 ≥ 0, the resulting scheme does give perfectly elastic impacts.However, this new scheme cannot be interpreted as a standard method for ODEs.This shows how much care must be taken with numerical simulation of impact phenomena.Other schemes, such as the Newmark schemes, can also be applied to contact problems [59].However, these schemes also suffer from indeterminism in the effective coefficient of restitution.In Figure 4.1(b), results are shown for a Newmark scheme [59] with β = 0, γ = 1, and step-size h = 10 −3 for the same bouncing ball problem with several initial conditions to illustrate the nondeterminism.Even variational methods based on discrete Lagrangian principles [141] do not seem to behave correctly.

Static vs. Dynamic Friction.
It is an experimentally observed phenomenon that the coefficient of friction often depends on the sliding velocity as well as the nature of the materials in contact.This has a number of practical implications.For example, if you need to brake suddenly while driving, the common advice is to "pump" the brakes.(This advice does not apply to vehicles equipped with ABS, which "pumps" automatically.)This is because when the car begins to slide the coefficient of friction appears to decrease, and the car does not decelerate as quickly.(Many people perceive this effect as the car actually accelerating; in fact it is still decelerating, but not as quickly.) The standard Coulomb model of dry friction does not take this velocity dependence into account.There are several ways of incorporating velocity dependence into the Coulomb model.The simplest is to have two coefficients of friction: the static coefficient µ static and the dynamic coefficient µ dynamic .From the differential inclusion point of view, it is clear that µ dynamic ≤ µ static , for if µ dynamic > µ static , there could not be any motion until the friction limit µ dynamic N is overcome.Thus there is no sense having µ dynamic > µ static .This discontinuous model of friction leads to a number of theoretical difficulties.For example, consider the "brick on a ramp" problem with this model of friction.Treating the problem naively, the equations of motion become and the right-hand side contains the right-hand side for the original friction model (1.1) with µ = µ dynamic .This means that if we interpret this simply as a differential inclusion, the brick can begin sliding whenever µ dynamic mg cos θ < mg sin θ < µ static mg cos θ.Clearly this is missing something crucial about the two-friction model: there is a hysteresis effect.Once v = 0, the coefficient of friction is "stuck" at µ static until such time as the external force is sufficient to overcome the static friction force limit µ static N .This model is vulnerable to some strong criticisms, such as those of Ruina [115, section 9.3.1],who discusses the model from a continuum point of view.
Even when the hysteresis effect is understood and is used to remove most of the nonuniqueness associated with this model, it is still quite possible to get multiple solutions when the velocity function v(t) "grazes" v = 0; that is, v (t) < 0 for t < t * , v(t * ) = 0, v (t * ) = 0, but v (t * ) > 0. Then the friction coefficient may jump from µ dynamic to µ static and "grab" the solution, forcing it to stay at zero for an interval.Or perhaps the coefficient of friction might stay at µ dynamic and allow the velocity to move away from zero again.How to treat this nonuniqueness is an unresolved theoretical issue with this model.
Another model which has a stronger experimental basis is to write µ = µ( v 2 ), where µ (s) < 0 and lim s→∞ µ(s) = µ ∞ > 0. An example of this model is shown in Figure 4.2(c).The behavior of systems with this friction model is very similar to the common understanding of the two-friction-coefficients model.If v = 0, the fact that µ ( v 2 ) < 0 actually makes the system unstable, in spite of the energy dissipation due to µ( v 2 ) > 0. Consider again the "brick on a ramp" problem of Figure 1.1 with a small external force F ext (t): Suppose that at time t = 0, v = 0. Then if we make F ext (t) sufficiently large, so that mg sin θ + F ext (t) > mg cos θ µ(0), the brick will begin sliding.But then m(dv/dt) = mg sin θ − F ext (t) − mg cos θ µ(|v|).As v increases, µ(|v|) decreases, so dv/dt will increase, raising the acceleration still higher.This is an example of frictioninduced instability.If µ (|v|) approaches µ ∞ quickly, then µ (|v|) will be large and negative, giving a rapid transition to the regime where µ(|v|) ≈ µ ∞ .If we identify µ static with µ(0) and µ dynamic with µ ∞ , the two models produce very similar behavior, although the theoretical difficulties of the µ = µ( v 2 ) model are much less.From the point of view of experiments, the frictional instability means that measurements of these effects require strong viscous damping in order to obtain a proper relationship between µ and v. Nevertheless, there are a number of open questions: For example, in the case of Painlevé problems with velocity-dependent friction, can the sliding velocity be brought to zero instantaneously even if µ( v − 2 ) is below the threshold normally needed for the Painlevé effect?4.4.Multiple Contacts and Multiple Impacts.Multiple impact and multiple contact problems can be formulated in the same style as the one-contact formulations here.These can be found in, for example, [5,6,132,129].The essence of these formulations is that the friction cones in generalized coordinates are summed to give  the complete friction cone for the entire system: F C(q) = F C (1) (q) + F C (2) (q) + • • • + F C (p) (q) for a system with p contacts.Each contact has its own constraint f (j) (q) ≥ 0, coefficient of friction µ (j) , coefficient of restitution ε (j) , normal contact force n (j) (q)c (j) n , friction force D (j) (q)β (j) , and λ (j) .The complementarity conditions for a single contact are replicated across all of the contacts.
The difficulty with multiple contacts is the correct way of dealing with the continuous formulations of the impact laws.For example, for inelastic impacts, n (j) (q) T v + = 0 if f (j) (q) = 0 for all j does not hold.Consider, for example, a pair of balls colliding inelastically on a table as illustrated in Figure 4.3.If the coefficients of friction are all strictly positive, with ball 1 initially rolling and ball 2 initially at rest on the table, then on impact ball 1 will jump up due to a vertical frictional impulse at contact 3 on ball 1.This means that while contact 2 and contact 3 both behave inelastically, contact 1 appears not to.Correct and complete impact laws for multiple contacts and impacts are not known even for inelastic impacts.Other difficulties arise with formulating multiple contacts.One is to obtain even numerical formulations which are dissipative where the coefficients of restitution ε (j) are different for the different contacts.(One formulation put forward by Pfeiffer and Glocker in [50] was later found not to be dissipative [3].) As well as the difficulties in formulating laws for multiple contacts, there are a number of physical effects that should be considered as well.One of the most important of the properties that are lost in current rigid-body models is the concept of "relative stiffness."Consider, for example, three masses that can affect each other through springs.
The three masses in Figure 4.4 interact through the two springs with spring constants k 1 and k 2 .The rigid limit is obtained by taking k 1 , k 2 → ∞.However, there are actually many rigid limits depending on the ratio k 1 /k 2 .If mass m 2 is in contact with mass m 3 , and m 1 collides with m 2 from the left, then several outcomes are possible.If k 1 /k 2 1, then since the time constant for the m 2 -m 3 interaction is much shorter than the time constant for the m 1 -m 2 interaction, m 2 and m 3 will appear to be a single unit, and will leave together.On the other hand, if k 1 /k 2 1, then the impacts will appear to be consecutive: momentum is transferred to mass m 2 , which is then transferred to mass m 3 , and only m 3 will leave with a significant velocity.Clearly, models of impact with multiple bodies and contacts should take into account these stiffness-ratio effects.
High-speed photographs of rows of initially touching clear plastic disks under polarized light on impact clearly show that throughout the impact there are usually three or four disks that are simultaneously in contact.This seems to indicate that collisions with multiple contacts cannot be replaced by sequences of two-body collisions.When "stiffness ratio" effects are included, problems with multiple contacts become quite complex.The development of extended rigid-body models that can accurately model these situations is beyond the scope of this paper.
A further issue that is related to multiple simultaneous contacts and vibrations of stiff elastic bodies is raised by Chatterjee [17].Consider a uniform slender rod overhanging the edge of a frictionless table and then striking the overhanging edge of the rod, as illustrated in Figure 4.5.Treating the bodies simply as rigid bodies would give the picture (a), where a reaction impulse occurs at point B, while at D, the rod moves up off the table so that there is no impulse there.However, the rod is not perfectly rigid, and this can be modeled using a torsion spring at a point C (say, midway between B and D).Then since AC is considered rigid, immediately after A is struck, the vertical velocity at C is positive.There is an upward impulse on CD at C. Since CD is rigid, some calculations show that there must be an upward impulse at D to prevent D penetrating the table.The exact value of this impulse at D depends on the location of C, but does not depend on the stiffness of the torsion spring.Since the forces at D can only be upwards, later motion could only increase the total impulse applied at D. Thus the limit as the stiffness of the torsion spring goes to infinity does not recover the rigid-body solution.
Of course, the torsion spring model is only an approximation; a more realistic approximation is to model the rod as an elastic beam.Similar effects occur there, although the analysis has not been carried out yet, in part because of the difficulty of dealing with infinite-dimensional problems.(See the following section for more on these issues.)Simple experiments with a rod (like a knitting needle) and carbon paper verify that this effect is real.Thus we have both n T D v + D > 0 and N D > 0 in the limit as k spring → ∞, contradicting the usual assumption about complementarity: 0 ≤ n T D v + D ⊥ N D ≥ 0. Nevertheless, it should be noted that the usual complementarity conditions hold for all times and spring stiffnesses in (b) and (c).In addition, complementarity between the vertical displacement and the normal contact force remains fundamental, as noted in [17].
Can the rigid-body solution (a) be obtained as some sort of limit of (b) or (c)?The answer is yes, and again is related to stiffness ratios: if k table /k spring goes to zero as both stiffnesses go to infinity, the rigid body solution (a) is recovered.The reason that this is not seen experimentally in [17] is that it is much easier to build a stiff table than it is to find a stiff yet slender rod.Again we come to a question of how to deal with stiffness ratios and elastic vibrations.The problem is that as the material stiffness goes to infinity, the amplitudes of the elastic vibrations go to zero, but the amplitudes of the surface velocities do not.This means that complementarity conditions for rigid motions may not be correct even though complementarity conditions for the true surface velocities may be.This issue is more fundamental than the issue of which impact law to use, since all current impact laws must handle the inelastic case ε = 0 at least as a model of the compression phase of the contact.In short, there are a number of issues regarding multiple contact and multiple body impacts that need resolution, but new low-order models may be able to incorporate the effects of multiple contacts and elastic vibrations.

Elastica in Impact.
The problems in handling the infinite-dimensional problems associated with elastic (and elastoviscous, elastoplastic) materials in frictional contact are substantial.Currently, the best results are either for quasi-static problems where the forces are assumed to always be in equilibrium, or for special dynamic problems with "small" coefficients of friction.Good starting points for understanding these problems are the book by Kikuchi and Oden [65] and the book by Duvaut and Lions [32].However, the early formulations in terms of variational inequalities had to assume frictionless contact, that the normal contact force is given (Tresca friction), or some other simplifying assumption.The feedback between the friction force and the normal contact force was usually ignored in order to obtain a well-posed variational inequality.If the feedback is used to construct a fixed-point problem, then results can be obtained for coefficients of friction "small enough" by means of a contraction mapping argument.For large coefficients of friction, this approach fails.An alternative that uses a Leray-Schauder argument would be more attractive in this case, although there are a number of technical difficulties with compactness, even for quasi-static problems.Often a mollifier is used to modify the friction law to satisfy the compactness requirements, but this also makes the friction law nonlocal, which does not seem physically appropriate.The following references survey the state of the art in this area [22,23,30,56,57,58,62,65].
Work on fully dynamic infinite-dimensional contact problems includes that of P. Shi [125], Bamberger and Schatzman [10], Schatzman [120], and Schatzman and Bercovier [121].The work of P. Shi mainly concerns a one-dimensional wave equation with contact only at one end; the work of Schatzman et al. is mainly about the wave equation with a rigid obstacle, where contact can be made over the whole spatial domain.The problems in this area are characterized by a high level of technical difficulty.A full discussion of the issues in this area is beyond the scope of this article.

Fig. 2 . 3
Fig. 2.3 Elevation and plan views of "billiards" problem with partly elastic impacts.

Fig. 4 . 1
Fig. 4.1 Numerical trajectories using the midpoint rule and a Newmark scheme.
[137]uld be a closed set.The values F (t, x) should all be closed, bounded, convex sets.And F (t, x) should satisfy a condition to prevent "blow-up" in finite time, such as x T z ≤ C(1 + x 2 ) for all z ∈ F (t, x).Numerical methods for discontinuous ODEs need to use this differential inclusion formulation if high accuracy is desired.If this is not done, then the methods typically have first order convergence due to rapid "chattering" of the numerical trajectories around the discontinuities for simple discontinuous ODEs.The first published results on numerical methods for discontinuous ODEs and differential inclusions were those by Taubert[137]in 1976.Further work on numerical methods for differential inclusions and discontinuous ODEs includes where F is a set-valued function.There are some properties that F should have.Its graph { (x, y) | y ∈ F (t, x) and is called the variation of v on [0, T ].Then we can use Riemann-Stieltjes integrals to define φ(t) dv(t) for continuous functions φ; see, for example, [134, pp.281-284] or [72, Chap.X].By the Riesz representation theorem, the linear functional φ −→ φ(t) dv(t) is a continuous linear functional and is therefore equivalent to integration against a measure: φ(t) dv(t) = φ(t) µ(dt), where µ is a Borel measure [72, Thm.2.7, Chap.IX, pp.264-265].
Rod and table example of Chatterjee: (a) rigid bodies; (b) torsion spring approximation;(c) elastic beam.