c 2000 Society for Industrial and Applied Mathematics
SIAM REVIEW
Vol. 42, No. 1, pp. 3–39
Rigid-Body Dynamics with
Friction and Impact∗
David E. Stewart†
Abstract. 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.
Key words. rigid-body dynamics, Coulomb friction, contact mechanics, measure-differential inclusions, complementarity problems
AMS subject classifications. Primary, 70E55; Secondary, 70F40, 74M
PII. S0036144599360110
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
∗ Received by the editors July 26, 1999; accepted for publication (in revised form) August 5, 1999;
published electronically January 24, 2000. This work was supported by NSF grant DMS-9804316.
http://www.siam.org/journals/sirev/42-1/36011.html
† Department of Mathematics, University of Iowa, Iowa City, IA 52242 (dstewart@math.
uiowa.edu).
3
4
DAVID E. STEWART
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
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
5
F
v
mg
N
θ
Fig. 1.1 Brick on a frictional ramp.
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).
1.1. 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
(1.1)
m
dv
= mg sin θ − µmg cos θ sgn v,
dt
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,
6
DAVID E. STEWART
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. Filippov around 1960. A differential inclusion
has the form
dx
∈ F (t, x),
dt
where F is a set-valued function. There are some properties that F should have. Its
graph { (x, y) | y ∈ F (t, x) } should 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 xT z ≤ C(1 + x2 ) 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 [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 cn , the friction force cf is the one that maximizes the
rate of energy dissipation −cTf vrel , where vrel is the relative velocity at the contact,
out of all possible friction forces allowed by the given normal contact force cn . To
be more formal, there is a set F C0 which is the set of possible friction forces cf for
cn = 1. The set F C0 is assumed to be closed, convex, and balanced (F C0 = −F C0 ).
So the maximum dissipation principle says that
(1.2)
cf
maximizes
− cTf vrel
over cf ∈ cn F C0 .
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
7
In the case of two-dimensional isotropic friction acting on a particle, F C0 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 cn + cf . The set of
possible contact forces is the friction cone, which is given by
(1.3)
F C = { n cn + cf | cn ≥ 0 and cf ∈ cn F C0 }.
Of course, as the contact point changes, so does the plane of the possible friction
forces. So we must allow F C0 and F C to depend on the configuration of the system:
F C0 = F C0 (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 C0 .
Note that F C0 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}.
1.2. 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 nT v(t+ ) = −ε nT v(t− )) [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: Ndecompr = −ε Ncompr
[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
nT 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
8
DAVID E. STEWART
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 ε = −(nT v + )/(nT 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
9
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
ek
collision time [µs]
140
b)
collision time [µs]
ek
collision time [µs]
ek
800
10
9
8
7
6
5
4
3
2
1
600
400
200
15
13
11
9
7
5
3
1
19
17
15
13
11
9
7
5
3
1
collision time [µs]
2000
1750
1500
1250
1000
750
500
250
4000
collision time [µs]
ek
ek
200
No. of micro-collisions
300
5
4
3
2
1
No. of micro-collisions
400
No. of micro-collisions
0
90
1000
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
100
No. of micro-collisions
e)
1
110
500
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
d)
120
80
1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
c)
2
130
No. of micro-collisions
a)
Simulation data
Experimental data
10
20
30
40
50
φ [deg]
60
70
80
90
3000
2000
1000
0
Collision time
Number of
micro-collisions
10
20
30
40
50
60
70
80
90
φ [deg]
Fig. 1.2 Results of Stoianovici and Hurmuzlu showing the effective kinematic coefficient of restitution (left), collision time and number of “micro-collisions” (right) for bars of lengths
(a) 100 mm, (b) 200 mm, (c) 300 mm, (d) 400 mm, and (e) 600 mm. Reprinted from
D. Stoianovici and Y. Hurmuzlu, ASME Journal of Applied Mechanics, 63 (1996), with
permission from the American Society of Mechanical Engineers.
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.
1.3. 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.
10
DAVID E. STEWART
Mass
m
Moment of inertia
J
µ
Coefficient of friction
l/2
(x, y)
l/2
mg
θ
velocity
F= µ N
(xc , yc )
N
Fig. 1.3 Painlevé’s problem.
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/ml2 .
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 (xc , yc ). 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
(1.4)
mẍ = F,
mÿ = N − mg,
J θ̈ = (l/2)[+F sin θ − N cos θ].
Of particular importance is the equation for yc , since the rigid-body condition is that
yc ≥ 0. If yc = 0 and ẏc = 0, then to prevent penetration we need ÿc ≥ 0. This can be
computed in terms of F and N since yc = y − (l/2) cos θ. Differentiating twice gives
ÿc = ÿ + (l/2) sin θ θ̈ + (l/2) cos θ θ̇2 .
Substituting F = +µN into these equations gives
l2
1
−
cos θ(µ sin θ − cos θ) N + (l/2) cos θ θ̇2 − g.
ÿc =
(1.5)
m 4J
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
(1.6)
0 ≤ ÿc
⊥
N ≥ 0.
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
11
(Note that “a ⊥ b” indicates that ab = 0 if a and b are scalars and that aT b = 0 if they
are vectors; for a vector “a ≥ 0” means that the components ai of a are nonnegative.)
Now consider the situation with large coefficient of friction µ and small J/ml2 . 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/ml2 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) − (l2 /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–
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.
The reader wishing to look at the extensive discussion of Painlevé’s problem is
invited to peruse the references [11, 12, 27, 28, 36, 49, 68, 76, 84, 82, 89, 91, 109].
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.
1.4. 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 (yc ), for example.
Complementarity problems typically have the following form: Find z ∈ Rn such
that 0 ≤ z ⊥ f (z) ≥ 0, where f : Rn → Rn is a given continuous function. Note that
this is equivalent to requiring that zi ≥ 0 and fi (z) ≥ 0 for all i = 1, . . . , n, and that
zi = 0 or fi (z) = 0 for all i = 1, . . . , n. Linear complementarity problems (LCPs)
are complementarity problems where f is affine: f (z) = M z + 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
12
DAVID E. STEWART
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
(1.7)
min f (x)
x
subject to gi (x) ≥ 0,
i = 1, . . . , m,
then provided a suitable constraint qualification [41] holds, there are Lagrange multipliers λi such that
m
∇f (x) − i=1 λi ∇gi (x) = 0,
λi , gi (x) ≥ 0,
(1.8)
λi gi (x) = 0
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 ⊥ M z + q ≥ 0,
where M is a copositive matrix (defined below) and q T y ≥ 0 whenever y solves the
homogeneous LCP 0 ≤ y ⊥ M y ≥ 0 [25, Thms. 3.8.6 and 4.4.12]. A matrix M
is copositive if y ≥ 0 implies that y T M y ≥ 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 K ∗ = { y ∈ X ∗ | y, x ≥ 0 for all x ∈ K }.
For a function f : X → X ∗ the complementarity problem CP (f ) is the problem
of finding z ∈ K such that f (z) ∈ K ∗ and f (z), z = 0. 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.
1.5. 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
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
13
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,
(1.9)
dv
∈ F (t, x),
dt
dx
= g(t, x, v),
dt
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
N
−1
i=0
v(ti+1 ) − v(ti )
over all N
Tand choice of 0 ≤ t0 ≤ t1 ≤ · · · ≤ tN ≤ T is finite. This supremum is
denoted 0 v 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]. 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”:
K∞ = x | x = lim tk xk , xk ∈ K .
(1.10)
k→∞
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
|µ|(E) = sup
µ(Ei ),
{Ei }
i
14
DAVID E. STEWART
where {Ei } 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
(1.11)
(1.12)
a(t) ∈ F (t, x(t)), Dt a.e.,
dµs
(t) ∈ F (t, x(t))∞ , |µs | a.e.
d|µs |
An alternative definition that is useful for handling problems of convergence is the
following: for each continuous function φ ≥ 0, not everywhere zero,
φ(t) dv(t)
∈ co
(1.13)
F (τ, x(τ )).
φ(t) dt
τ :φ(τ )=0
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).
• (C2) The asymptotic cone F (t, x)∞ is always pointed. (A pointed cone K is
a cone where K ∩ (−K) = {0}.)
• (C3) min{ z | z ∈ F (t, x) } is a bounded function of (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
h
sequence of functions
procedure. If it can
h {v }h>0 generated by some numerical
be proved that v ≤ M for some constant M and v h (0) = v0 for all h > 0,
then by the Helly selection theorem [93], there isa pointwise convergent subsequence
v hk (.) which converges to a function v(.) with v ≤ M . In this subsequence, the
measures Dv hk ⇀ Dv weak*. If v h satisfies a measure differential inclusion dv h /dt ∈
Fh (t, x(t)), and 0 < h′ < h implies that Fh′ (t, x) ⊂ Fh (t, x), then using the weak
formulation it is easy to show that dv/dt ∈ Fh (t, x(t)) for all h > 0, and consequently
dv/dt ∈ h>0 Fh (t, x(t)).
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.
2. 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
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
15
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.
2.1. 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) = 12 v T M (q)v−V (q), where 21 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:
d
dt
∂L
∂v
−
∂L
=0
∂q
or M (q)
dv
= k(q, v) − ∇V (q),
dt
where kl (q, v) = − 12 r,s [∂mir /∂qs + ∂mis /∂qr − ∂mrs /∂qi ]. With external and contact forces we can write this as
(2.1)
M (q)
dv
= n(q) cn + D(q) β + k(q, v) − ∇V (q) + Fext (t).
dt
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) ⊥
cn ≥ 0. We represent the set of possible friction forces through F C0 (q) = { D(q)β |
β ∈ Rd , ψ(β) ≤ µ }. The function ψ used for defining F C0 (q) should be convex,
positively homogeneous (ψ(αβ) = |α| ψ(β) for all α ∈ R), and coercive (ψ(β) → ∞
as β → ∞). Positive homogeneity implies that D(q)β ∈ cn F C0 (q) is equivalent to
ψ(β) ≤ µcn .
The maximal dissipation principle states that we choose β so as to maximize
the dissipation rate v T D(x)β over all β ∈ cn F C0 (x). Because of discontinuities, we
need to be careful to make some distinctions between v + (t) = v(t+ ) = lims↓t v(s),
v − (t) = v(t− ) = lims↑t v(s), and v(t). For inelastic impacts, for example, we require
that nT 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 β ∈ cn F C0 (q), that is,
(2.2)
min(v + )T D(q)β
β
subject to ψ(β) ≤ µcn .
16
DAVID E. STEWART
This is a convex program (with linear objective function), although it is nonsmooth.
If cn > 0, then using a Slater condition, there exists a Lagrange multiplier λ ≥ 0
such that ∂β h(β, λ) = 0 where h(β, λ) = v T D(q)β − λ(µcn − ψ(β)) and ∂f (z) is the
generalized gradient of f at z [19, Thms. 6.1.1 and 6.3.1]. Furthermore, 0 ≤ λ ⊥
µcn − ψ(β) ≥ 0. That is,
(2.3)
0 ∈ µ D(q)T v + + λ ∂ψ(β),
0 ≤ λ ⊥ µcn − ψ(β) ≥ 0.
If cn = 0, then β = 0 by the condition µcn − ψ(β) ≥ 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.
2.1.1. Formulations and Function Spaces. The Signorini contact condition that
0 ≤ f (q(t)) ⊥ cn ≥ 0 involves complementarity between a continuous function
t → f (q(t)) and a measure cn , so the complementarity condition can be represented
by f (q(t)) cn (dt) = 0 as well as the usual inequality conditions: f (q(t)) ≥ 0 for
all t and cn (E) ≥ 0 for
all Borel sets E ⊂ [0, T ]. The integrated maximum dissipation principle minβ (v + )T D(q)β subject to ψ(β) ≤ µcn 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 ψ(β) ≤ µcn
is equivalent to ψ(dβ/d|β|) |β| ≤ µcn or even ψ(dβ/dcn ) ≤ µ a.e. Note that the
derivatives dβ/d|β| and dβ/dcn are all Radon–Nikodym derivatives. The function λ
is a bounded Borel function; we can bound λ∞ by maxt µD(q(t))T v + (t)∞ /Cψ,
where Cψ = min{ z∞ | z ∈ ∂ψ(β), β = 0 }.
2.1.2. 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)cn +D(q)β | ψ(β) ≤
cn }. 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 cn for the normal contact forces and β to describe the frictional forces, and a
bounded Borel-measureable function λ where
dv
= n(q) cn + D(q) β − ∇V (q) + k(q, v) + Fext (t),
(2.4)
M (q)
dt
dq
(2.5)
= v,
dt
(2.6)
0 ≤ cn ⊥ f (q) ≥ 0,
(2.7)
0 ∈ µ D(q)T v + + λ ∂ψ(β),
(2.8)
(2.9)
0 ≤ λ ⊥ µcn − ψ(β) ≥ 0,
0 = n(q(t))T v + (t) if f (q(t)) = 0.
For partly or fully elastic impacts with coefficient of restitution 0 ≤ ε ≤ 1, we can
replace (2.9) with
(2.10)
n(q(t))T (v + + εv − ) = 0
if f (q(t)) = 0.
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
17
F
µN
v
‘‘regularized’’
friction law
Fig. 2.1 Regularized Coulomb law to avoid discontinuity.
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.
2.2. 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.
2.2.1. 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) cn and β over each timestep interval [tl , tl+1 ] is used. Two different numerical formulations are presented here.
The first is based on linear complementarity problems and uses a polyhedral approxi
mation 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) + µ di (q) | i = 1, 2, . . . , m }, where µ di (q) is a collection of direction vectors
in F C0 (q). Write D(q)
= [d1 (q), d2 (q), . . . , dm (q)]. The friction forces D(q)β are
18
DAVID E. STEWART
n
d7
d8
d6
ct
d5
d1
d2
d3
d4
polyhedral approximation
to friction cone
friction cone
Fig. 2.2 Polyhedral approximation to the friction cone.
where βi ≥ 0 and βi ≤ µ cn . The relationship between
β,
approximated by D(q)
i
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 di (q) = −dj (q). This is related to the assumption that F C0 (q) is a balanced
set: F C0 (q) = −F C0 (q).
The discretization of (2.4)–(2.9) is the problem of finding q l+1 and v l+1 (and the
l+1 and Lagrange multiplier λl+1 ) given q l and v l for a time-step
force integrals cl+1
n , β
of size h > 0 that satisfy the following conditions:
(2.11)
l l+1
M (q l+1 )(v l+1 − v l ) = n(q l )cl+1
n + D(q )β
+h[−∇V (q l ) + k(q l , v l ) + Fext (tl )],
(2.12)
(2.13)
(2.14)
(2.15)
l+1
l
q
− q = h v l+1 ,
0 ≤ cl+1
⊥ n(q l )T (v l+1 + εv l ) ≥ 0,
n
l )T v l+1 ≥ 0,
0 ≤ βl+1 ⊥ λl+1 e + D(q
0 ≤ λl+1 ⊥ µ cl+1 − eT βl+1 ≥ 0,
n
= 0 and βl+1 = 0 and
where f (q l + h v l ) < 0; if f (q l + h v l ) ≥ 0, then we set cl+1
n
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) ⊥ cn ≥ 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
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
19
0 ≤ f (q l+1 ) ⊥ cl+1
≥ 0. This does not work: the resulting discretization behaves
n
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 [tl , tl+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
l+1 l+1
of cn , β , and λl+1 , then (2.11)–(2.15) with ε = 0 can be reduced to a pure LCP
which has the form (superscripts suppressed)
T −1
T
0
cn
cn
n b1
n M n nT M −1 D
T b2 ⊥ β ≥ 0.
T M −1 D
T M −1 n D
e β + D
0≤ D
(2.16)
T
0
λ
λ
µ
−e
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].
To avoid approximations of the friction cone, we can use the integrated maximal dissipation principle: β l+1 maximizes −(v l+1 )T D(q l )β l+1 over all β l+1 satisfying
ψ(β l+1 ) ≤ µcl+1
n . The Kuhn–Tucker conditions for this problem replace (2.14), (2.15)
above, giving
(2.17)
(2.18)
0 ∈ D(q l )T v l+1 + λl+1 ∂ψ(β l+1 ),
l+1
0 ≤ λl+1 ⊥ µ cl+1
) ≥ 0.
n − ψ(β
Of course, (2.11) should be replaced by
(2.19)
l
l+1
M (q l+1 )(v l+1 − v l ) = n(q l )cl+1
n + D(q )β
+h[−∇V (q l ) + k(q l , v l ) + Fext (tl )].
Solving these systems requires more sophisticated methods since we have to solve for
an inclusion 0 ∈ D(q l )T v l+1 + λl+1 ∂ψ(β l+1 ) as well as for nonlinear complementarity
l+1
conditions 0 ≤ λl+1 ⊥ µcl+1
) ≥ 0.
n − ψ(β
The existence of solutions to the one-time-step conditions above can be shown
under mild conditions via the results in [129] for the formulation using β and extended
via [101] to the formulation(s) in β, using the nonlinear condition µcn − ψ(β) ≥ 0.
2.3. Practicalities. If we use the discretization (2.11)–(2.15), we can use Lemke’s
algorithm in an iteration: given an estimate q l+1,k of q l+1 , we can obtain an estimate
q l+1,k+1 by solving (2.11)–(2.15) with q l+1 replaced by q l+1,k . This iteration will
usually converge (and converge quickly), giving a solution of the mixed nonlinear
complementarity problem (2.11)–(2.15) [132].
Replacing (2.14)–(2.15) with (2.17)–(2.18), using the nonlinear condition µcn ≥
ψ(β) instead of the conditions µcn ≥ eT β 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 µcn ≥ ψ(β) with (µcn )2 ≥ β22 = β T β. The difficulty in
using this condition is that if cn = 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, section 2.2,
pp. 190–203], for example): replace (2.17) with
l T l+1
0 = cl+1
+ 2λl+1 β l+1 .
n D(q ) v
20
three initially
stationary balls
thrown
ball
thrown
ball
three initially
stationary balls
Note change of direction of travel
with each bounce due to the ball’s spin.
Note ‘‘jump’’ due to spin of first ball and
frictional impulses.
DAVID E. STEWART
(a) Plan view
(b) Elevation view
Fig. 2.3 Elevation and plan views of “billiards” problem with partly elastic impacts.
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
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
21
0.060
0.080
0.040
0.100
0.020
0.120
0.140 0.000
0.160
0.180
0.200
0.210
0.220
0.230
0.240
0.250
0.260
0.270
0.280
0.290
0.300
0.310
0.320
0.330
0.340
0.350
0.360
0.370
0.380
0.390
0.400
0.410
0.420
0.430
0.440
0.450
0.460
0.470
0.480
0.490
0.500
0.510
0.520
0.530
0.540
0.550
Fig. 2.4 Painlevé’s problem: impact version.
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.
3. 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.
22
DAVID E. STEWART
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))T v(t− ) ≥ 0
implies that v(t+ ) = v(t− ) under Erdmann’s condition.
The question of convergence is answered affirmatively in [129] for one inelastic
contact (ε = 0) in the case of one-dimensional friction (that is, F C0 (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].
3.1. 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 (v l+1 )T M (v l+1 − v l ) = 12 (v l+1 )T M v l+1 −
1 l+1
1 l T
l
− v l )T M (v l+1 − v l ) and substituting for M (v l+1 − v l ) in terms
2 (v ) M v + 2 (v
l+1
l+1
of cn and β . 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 [t0 , t1 ], t1 > t0 . This means that the computed velocities v h (.)
are bounded on [t0 , t1 ]. 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
h
closed graph,
it can be concluded from the uniform boundedness of v (.) 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.
3.2. 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 21 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
T
stage needs a uniform bound on 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, ∞).
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
23
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 cn is contained in the set { t | f (q(t)) = 0 }. This is
equivalent to showing that cn ({ t | f (q(t)) > 0 }) = 0. This can be done by showing
that cn ({ t | f (q(t)) ≥ ǫ }) = 0 for any ǫ > 0. Thus we have the complementarity
condition that 0 ≤ f (q) ⊥ cn ≥ 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 ⊥ cl+1
≥ 0 for f (q l + hv l ) < 0.
n
l
l
(If f (q + hv ) ≥ 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
β T D(q)T v + = −µcn D(q)T v + ∞
as measures. Since the numerical velocity functions v h converge pointwise, and
the measures β h and chn only converge weak*, we cannot directly take limits, even
though this property holds for the discrete formulation. First, it is shown that
nT 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 )cl+1
+ D(q l )β l+1 ) is
n
l+1
l
expanded two ways: one by expanding (v
− v ) 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 chn n(q h )Tv h − µD(q h )T v h ∞ is bounded below by the measure
cn 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:
d
1 T
v M (q)v + V (q)
2
=
1
1 +
(v + v − )M (q) dv + v T
2
2
d
M (q) v dt + ∇V (q)T v dt.
dt
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
24
DAVID E. STEWART
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 cl+1
is “large”; this is shown to be
n
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.
3.3. Other Aspects.
3.3.1. 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
l2
1
−
cos θ(µ sin θ − cos θ) N + (l/2) cos θ θ̇2 − g ⊥ N ≥ 0.
0 ≤ ÿc =
m 4J
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]!)
3.3.2. 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
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
25
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 Rn
is given by
dH (A, B) = max(max{ d(a, B) | a ∈ A }, max{ d(b, A) | b ∈ B }),
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.
26
DAVID E. STEWART
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 (nT 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
1
2η
u′′λ + √ G(uλ − PK (uλ )) + (uλ − PK (uλ )) = f (t, uλ , u′λ ),
λ
λ
where PK (z) is the projection of z onto K, G(u, v) = (uT v)u/u2 if u = 0 and zero
if u = 0, and η is related to the coefficient of restitution ε. Note that the quantity
uλ − PK (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).
4.2. 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
∂H
dq
=+
(q, p),
dt
∂p
dp
∂H
=−
(q, p).
dt
∂q
Symplectic integrators preserve the differential forms dqi dpi . They do not exactly
conserve the Hamiltonian (which is the total energy), but they nearly conserve a
“numerical” Hamiltonian Hh (q, p), where h > 0 is the step-size used for the symplectic
method [116]. The numerical Hamiltonian Hh (q, p), is usually obtained from H via
Taylor series.
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
27
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):
xn+1 = xn + h f ((tn + tn+1 )/2, (xn + xn+1 )/2).
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
1 2
p + ψ(q) + gq, where ψ(q) = 0 for q ≥ 0 and ψ(q) = +∞ for
system is H(q, p) = 2m
q < 0. This gives the differential equations and inclusions
dq
1
= p,
dt
m
dp
∈ −∂ψ(q) − g.
dt
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
h
(pn + pn+1 ),
2m
h
∈ pn − ∂ψ((qn+1 + qn )/2) − hg.
2
qn+1 = qn +
pn+1
Writing Nn+1 for the element we choose from −(h/2)∂ψ((qn +qn+1 )/2) (that is, Nn+1
is the normal contact impulse), we get the mixed complementarity problem
pn+1 = pn + Nn+1 − hg,
h
(pn + pn+1 ),
qn+1 = qn +
2m
0 ≤ Nn+1 ⊥ (qn + qn+1 )/2 ≥ 0.
Since symplectic methods do not conserve H(q, p) and the numerical Hamiltonian
Hh (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 qn < 0 and pn < 0, then Nn+1 > 0 is needed to ensure that qn +qn+1 ≥
0; by complementarity, qn+1 = −qn > 0; therefore pn + pn+1 = 2(2m/h)qn+1 and
pn+1 = −pn + (4mqn+1 /h) > |pn |. The effective coefficient of restitution will depend
on qn /(hpn ). 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 ≤ Nn+1 ⊥ qn + qn+1 ≥ 0 with
one involving the momenta 0 ≤ Nn+1 ⊥ pn + pn+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.
28
DAVID E. STEWART
Numerical solution for midpoint rule with contact (h = 0.01)
30
25
20
q(t)
15
10
5
0
5
0
1
2
3
4
5
t
6
7
8
9
1.6
1.8
10
(a) Midpoint rule
Bouncing ball simulation using Newmark–type algorithm
1.4
β = 0, γ = 1
1.2
h = 103
1
q(t)
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
1
t
1.2
1.4
2
(b) Newmark scheme
Fig. 4.1 Numerical trajectories using the midpoint rule and a Newmark scheme.
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
29
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.
4.3. 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
dv
µdynamic if v = 0,
∈ mg sin θ − mg cos θ Sgn v
m
µstatic
if v = 0,
dt
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 µ = µ(v2 ),
where µ′ (s) < 0 and lims→∞ µ(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
30
DAVID E. STEWART
F
F
µs N
µd N
µN
v
(a)
v
(b)
F
µs N
| F | = µ (v ) N
µd N
v
(c)
Fig. 4.2 One-dimensional friction models: (a) standard, (b) two-coefficient, (c) continuously dependent coefficient of friction.
µ′ (v2 ) < 0 actually makes the system unstable, in spite of the energy dissipation
due to µ(v2 ) > 0. Consider again the “brick on a ramp” problem of Figure 1.1 with
a small external force Fext (t):
(4.1)
m
dv
∈ mg sin θ − mg cos θ µ(|v|) Sgn v + Fext (t).
dt
Suppose that at time t = 0, v = 0. Then if we make Fext (t) sufficiently large,
so that mg sin θ + Fext (t) > mg cos θ µ(0), the brick will begin sliding. But then
m(dv/dt) = mg sin θ − Fext (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 µ = µ(v2 ) 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
31
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
contact 3
ω
ball 1
ball 2
table
contact 1
contact 2
Fig. 4.3 Two balls colliding inelastically.
m1
k1
m2
k2
m3
frictionless table
Fig. 4.4 Three masses problem.
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
(j)
force n(j) (q)cn , 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 k1 and k2 . The rigid limit is obtained by taking k1 , k2 → ∞. However,
32
DAVID E. STEWART
rod before impact
rod after impact
B
C
D
A
B
D
table
B
D
A
A
table
table
torsion spring
(a)
(b)
(c)
Fig. 4.5 Rod and table example of Chatterjee: (a) rigid bodies; (b) torsion spring approximation;
(c) elastic beam.
there are actually many rigid limits depending on the ratio k1 /k2 . If mass m2 is in
contact with mass m3 , and m1 collides with m2 from the left, then several outcomes
are possible. If k1 /k2 ≪ 1, then since the time constant for the m2 -m3 interaction
is much shorter than the time constant for the m1 -m2 interaction, m2 and m3 will
appear to be a single unit, and will leave together. On the other hand, if k1 /k2 ≫ 1,
then the impacts will appear to be consecutive: momentum is transferred to mass
m2 , which is then transferred to mass m3 , and only m3 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
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
33
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 nTD vD
> 0 and ND > 0 in the
limit as kspring → ∞, contradicting the usual assumption about complementarity:
+
⊥ ND ≥ 0. Nevertheless, it should be noted that the usual complemen0 ≤ nTD vD
tarity 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 ktable /kspring 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.
4.5. 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.
34
DAVID E. STEWART
REFERENCES
[1] E. Allgower and K. Georg, Numerical Continuation Methods: An Introduction, SpringerVerlag, Berlin, Heidelberg, New York, 1990.
[2] J. Angeles and A. Kecskeméthy, eds., Kinematics and Dynamics of Multi-Body Systems,
Springer-Verlag, Vienna, 1995.
[3] M. Anitescu, private communication, Argonne National Laboratories, Argonne, IL, 1997.
[4] M. Anitescu, J. F. Cremer, and F. A. Potra, On the existence of solutions to complementarity formulations of contact problems with friction, in Complementarity and Variational
Problems: State of the Art, SIAM, Philadelphia, 1997, pp. 12–21.
[5] M. Anitescu and F. A. Potra, Formulating dynamic multi-rigid-body contact problems with
friction as solvable linear complementarity problems, ASME J. Nonlinear Dynamics, 14
(1997), pp. 231–247.
[6] M. Anitescu, F. A. Potra, and D. E. Stewart, Time-stepping for three-dimensional rigid
body dynamics, Comput. Meth. Appl. Mech. Engrg., 177 (1999), pp. 183–197.
[7] D. Bainov and P. Simeonov, Systems with Impulse Effect: Stability, Theory, and Applications, Math. Appl., Chichester, New York, 1989.
[8] D. Bainov and P. Simeonov, Impulsive Differential Equations: Periodic Solutions and Applications, Pitman Monogr. Surveys Pure Appl. Math. 66, Longman, Harlow, 1993.
[9] P. Ballard, A counter-example to uniqueness in quasi-static elastic contact problems with
small friction, Internat. J. Engrg. Sci., 37 (1999), pp. 163–178.
[10] A. Bamberger and M. Schatzman, New results on the vibrating string with a continuous
obstacle, SIAM J. Math. Anal., 14 (1983), pp. 560–595.
[11] D. Baraff, Fast contact force computation for non-penetrating rigid bodies, in Proceedings,
SIGGRAPH ’94, ACM, New York, 1994, pp. 23–34.
[12] H. Béghin, Sur certain problemes de frottement, Nouv. Ann. Math. Sér. 2, 5 (1923/1924),
pp. 305–312.
[13] M. Benati and A. Morro, Equations for multibody dynamics, Atti Accad. Peloritana Pericolanti Cl. Sci. Fis. Mat. Natur., 68 (1990), pp. 291–308.
[14] B. Brogliato, Nonsmooth Impact Mechanics: Models, Dynamics and Control, Lecture Notes
in Control and Inform. Sci. 220, Springer, New York, 1996.
[15] S. L. Campbell and B. Leimkuhler, Differentiation of constraints in differential-algebraic
equations, Mech. Structures Mach., 19 (1991), pp. 19–39.
[16] C. Castaing and M. D. P. Monteiro Marques, BV periodic solutions of an evolution problem associated with continuous moving convex sets, Set-Valued Anal., 3 (1995), pp. 381–
399.
[17] A. Chatterjee, On the realism of complementarity conditions in rigid-body collisions, Nonlinear Dynamics, 20 (1999), pp. 159–168.
[18] A. Chatterjee and A. L. Ruina, A new algebraic rigid-body collision law based on impulse
space considerations, ASME J. Appl. Mech., 65 (1998), pp. 939–950.
[19] F. H. Clarke, Optimization and Nonsmooth Analysis, Classics Appl. Math. 5, SIAM,
Philadelphia, 1990.
[20] P. W. Cleary, Modelling granular flows with complex boundary geometries, in Computational
Techniques and Applications: CTAC93, D. Stewart, D. Singleton, and H. Gardner, eds.,
World Scientific, Singapore, 1994, pp. 148–156.
[21] P. W. Cleary and N. G. Barton, Use of particle methods to design mineral processing
equipment, E. Kreuzer and O. Mahrenholtz, eds., Z. Angew. Math. Mech., 76 (1996),
pp. 253–256.
[22] M. Cocu, E. Pratt, and M. Raous, Existence d’une solution du problème quasi-statique
de contact unilatéral avec frottement non local, C. R. Acad. Sci. Paris Sér. I Math., 320
(1995), pp. 1413–1417.
[23] M. Cocu, E. Pratt, and M. Raous, Formulation and approximation of quasistatic frictional
contact, Internat. J. Engrg. Sci., 34 (1996), pp. 783–798.
[24] E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations, Tata–
McGraw Hill, New Delhi, 1972.
[25] R. W. Cottle, J.-S. Pang, and R. E. Stone, The Linear Complementarity Problem, Comput. Sci. Sci. Comput., Academic Press, Boston, San Diego, New York, 1992.
[26] C. A. de Coulomb, Théorie des machines simples, en ayant égard au frottement de leurs
parties, et la roideur des cordages. Piéce qui a reporté le Prix double de l’Académie
des Sciences pour l’année 1781, Mémoires des Savans Etrangers, X (1785), pp. 163–332.
Reprinted by Bachelier, Paris, 1809.
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
35
[27] E. Delassus, Considérations sur le frottement de glissement, Nouv. Ann. Math. Sér. 4, 20
(1920), pp. 485–496.
[28] E. Delassus, Sur les lois du frottement de glissement, Bull. Soc. Math. France, 51 (1923),
pp. 22–33.
[29] F. Demengel and R. Temam, Convex functions of a measure and applications, Indiana Univ.
Math. J., 33 (1984), pp. 673–709.
[30] L. Demkowicz and J. T. Oden, On some existence and uniqueness results in contact problems with nonlocal friction, Nonlinear Anal., 6 (1982), pp. 1075–1093.
[31] A. Dontchev and F. Lempio, Difference methods for differential inclusions: A survey, SIAM
Rev., 34 (1992), pp. 263–294.
[32] G. Duvaut and J.-L. Lions, Inequalities in Mechanics and Physics, Grundlehren Math. Wiss.
219, Springer-Verlag, Berlin, Heidelberg, New York, 1976.
[33] E. Eich, C. Führer, and J. Yen, On the error control for multistep methods applied to
ODEs with invariants and DAEs in multibody dynamics, Mech. Structures Mach., 23
(1995), pp. 159–179.
[34] E. Eich-Soellner and C. Führer, Numerical Methods in Multibody Dynamics, B. G. Teubner, Stuttgart, 1998.
[35] C. M. Elliott, On the convergence of a one-step method for the numerical solution of
ordinary differential inclusions, IMA J. Numer. Anal., 5 (1985), pp. 3–21.
[36] M. Erdmann, On a representation of friction in configuration space, Internat. J. Robotics
Research, 13 (1994), pp. 240–271.
[37] M. C. Ferris and J.-S. Pang, Engineering and economic applications of complementarity
problems, SIAM Rev., 39 (1997), pp. 669–713.
[38] A. F. Filippov, Differential equations with discontinuous right-hand side, Amer. Math. Soc.
Trans., 42 (1964), pp. 199–231; original published in Mat. Sb., 5 (1960), pp. 99–127 (in
Russian).
[39] A. F. Filippov, Classical solutions of differential equations with multivalued right-hand side,
SIAM J. Control Optim., 5 (1967), pp. 609–621.
[40] A. F. Filippov, Differential Equations with Discontinuous Right-Hand Side, Kluwer Academic, Norwell, MA, 1988.
[41] R. Fletcher, Practical Methods of Optimization, 2nd ed., John Wiley, Chichester, New York,
1987.
[42] C. G. Franchi, A general formulation in rigid multibody dynamics, Internat. J. Numer. Meth.
Engrg., 38 (1995), pp. 1985–2016.
[43] J. Franklin, Methods of Mathematical Economics, Springer-Verlag, New York, Heidelberg,
Berlin, 1980.
[44] H. Frankowska, The maximum principle for an optimal solution to a differential inclusion
with end points constraints, SIAM J. Control Optim., 25 (1987), pp. 145–157.
[45] C. Führer and B. Leimkuhler, Numerical methods for differential-algebraic equations describing constrained mechanical motion, in Computational Ordinary Differential Equations (London, 1989), Oxford University Press, New York, 1992, pp. 275–284.
[46] J. Garcı́a de Jalón and E. Bayo, Kinematic and dynamic simulation of multibody systems,
Springer-Verlag, New York, 1994.
[47] Z. Ge and J. E. Marsden, Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators,
Phys. Lett. A, 133 (1988), pp. 134–139.
[48] C. W. Gear, B. Leimkuhler, and G. K. Gupta, Automatic integration of Euler–Lagrange
equations with constraints, J. Comput. Appl. Math., 12/13 (1985) pp. 77–90.
[49] F. Génot and B. Brogliato, New results on Painlevé paradoxes, European J. Mech.
A Solids, 18 (1999), pp. 653–677. Also available as INRIA Research Report 3366 via
ftp://ftp.inria.fr/INRIA/tech-reports/RR/RR-3366.ps.gz.
[50] C. Glocker and F. Pfeiffer, Multiple impacts with friction in rigid multibody systems,
ASME J. Nonlinear Dynamics, 7 (1995), pp. 471–497.
[51] S.-P. Han, J.-S. Pang, and N. Rangaraj, Globally convergent Newton methods for nonsmooth equations, Math. Oper. Res., 17 (1992), pp. 586–607.
[52] E. J. Haug, D. Negrut, and M. Iancu, A state-space-based implicit integration algorithm
for differential-algebraic equations of multibody dynamics, Mech. Structures Mach., 25
(1997), pp. 311–334.
[53] J.-B. Hiriart-Urruty and C. Lemaréchal, Convex Analysis and Minimization Algorithms
I and II, Grundlehren Math. Wiss., Springer–Verlag, Berlin, Heidelberg, New York, 1993.
[54] R. Horst and P. Pardalos, eds., Handbook of Global Optimization, Nonconvex Optim.
Appl., Kluwer Academic Publishers, Dordrecht, the Netherlands, 1995.
36
DAVID E. STEWART
[55] G. Isac, Complementarity problems, Lecture Notes in Math. 1528, Springer-Verlag, Berlin,
1992.
[56] J. Jarušek, Contact problems with bounded friction: Coercive case, Czechoslovak Math. J.,
33 (1983), pp. 237–261.
[57] J. Jarušek and C. Eck, Dynamic contact problems with friction in linear viscoelasticity, C.
R. Acad. Sci. Paris Sér. I Math., 322 (1996), pp. 497–502.
[58] J. Jarušek and C. Eck, Dynamic contact problems with small Coulomb friction for viscoelastic bodies. Existence of solutions, Math. Models Methods Appl. Sci., 9 (1999), pp. 11–34.
[59] C. Kane and M. Ortiz, Finite element analysis of non-smooth contact, Comput. Meth. Appl.
Mech. Engrg., 180 (1999), pp. 1–26.
[60] A. E. Kastner-Maresch, Implicit Runge–Kutta methods for differential inclusions, Numer.
Funct. Anal. Optim., 11 (1990), pp. 937–958.
[61] A. E. Kastner-Maresch, The implicit midpoint rule applied to discontinuous differential
equations, Computing, 49 (1992), pp. 45–62.
[62] Y. Kato, Signorini’s problem with friction in linear elasticity, Japan J. Appl. Math., 4 (1987),
pp. 237–268.
[63] J. B. Keller, Impact with friction, ASME J. Appl. Mech., 53 (1986), pp. 1–4.
[64] J. L. Kelley and T. P. Srinivasan, Measure and Integral, Volume 1, Grad. Texts in Math.
116, Springer-Verlag, New York, 1988.
[65] N. Kikuchi and J. T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM Stud. Appl. Math. 8, SIAM, Philadelphia,
1988.
[66] A. Klarbring, Mathematical programming in contact problems, in Computational Methods
in Contact Mechanics, M. H. Aliabadi and C. A. Brebbia, eds., Internat. Ser. Comput.
Engrg., Computational Mechanics Publications, London, 1993, pp. 233–263.
[67] A. Klarbring, Steady sliding and linear complementarity, in Complementarity and Variational Problems, M. C. Ferris and J.-S. Pang, eds., SIAM, Philadelphia, 1997, pp. 132–147.
[68] F. Klein, Zu Painleves Kritik der Coulombschen Reibungsgesetze, Zeit. Math. Physik, 58
(1909), pp. 186–191.
[69] M. Kunze and M. D. P. Monteiro Marques, Yosida-Moreau regularization of sweeping
processes with unbounded variation, J. Differential Equations, 130 (1996), pp. 292–306.
[70] M. Kunze and M. D. P. Monteiro Marques, Existence of solutions for degenerate sweeping
processes, J. Convex Anal., 4 (1997), pp. 165–176.
[71] V. Lakshmikantham, D. Bainov, and P. Simeonov, Theory of Impulsive Differential Equations, World Scientific, Singapore, 1989.
[72] S. Lang, Real and Functional Analysis, 2nd ed., Grad. Texts in Math. 142, Springer-Verlag,
Berlin, Heidelberg, New York, 1993.
[73] B. Leimkuhler and S. Reich, Symplectic integration of constrained Hamiltonian systems,
Math. Comp., 63 (1994), pp. 589–605.
[74] F. Lempio, Modified Euler methods for differential inclusions, in Set-Valued Analysis and
Differential Inclusions, A. B. Kurzhansky and V. M. Veliov, eds., Birkhäuser, Boston,
Basel, Berlin, 1993.
[75] F. Lempio and V. Veliov, Discrete approximations of differential inclusions, Bayreuth.
Math. Schr., 54 (1998), pp. 149–232.
[76] P. Lötstedt, Coulomb friction in two-dimensional rigid-body systems, Z. Angew. Math.
Mech., 61 (1981), pp. 605–615.
[77] P. Lötstedt, Mechanical systems of rigid bodies subject to unilateral constraints, SIAM J.
Appl. Math., 42 (1982), pp. 281–296.
[78] P. Lötstedt, Time-dependent contact problems in rigid-body mechanics, Math. Prog. Study,
17 (1982), pp. 103–110.
[79] P. Lötstedt, Numerical simulation of time-dependent contact and friction problems in rigid
body mechanics, SIAM J. Sci. Statist. Comput., 5 (1984), pp. 370–393.
[80] M. Mabrouk, Liaisons unilatérales et chocs élastiques quelconques: Un résultat d’existence,
C. R. Acad. Sci. Paris Sér. I Math., 326 (1998), pp. 1353–1357.
[81] M. Mabrouk, A unified variational model for the dynamics of perfect unilateral constraints,
European J. Mech. A Solids, 17 (1998), pp. 819–842.
[82] M. T. Mason and Y. Wang, On the inconsistency of rigid-body frictional planar mechanics,
in Proceedings, IEEE International Conference on Robotics Automation, IEEE, Piscataway, NJ, 1988, pp. 524–528.
[83] M. D. P. Monteiro Marques, Chocs inélastiques standards: un résultat d’existence, Sém.
Anal. Convexe, 15 (1985), pp. 1–32.
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
37
[84] M. D. P. Monteiro Marques, Differential Inclusions in Nonsmooth Mechanical Problems:
Shocks and Dry Friction, Prog. Nonlinear Differential Equations Appl. 9, BirkhäuserVerlag, Basel, Boston, Berlin, 1993.
[85] J.-J. Moreau, Sur les mesures différentielles de fonctions vectorielles et certain problèmes
d’évolution, C. R. Acad. Sci. Paris Sér. A, 282 (1976), pp. 837–840.
[86] J.-J. Moreau, Evolution problem associated with a moving convex set in a Hilbert space, J.
Differential Equations, 26 (1977), pp. 347–374.
[87] J.-J. Moreau, Application of convex analysis to some problems of dry friction, in Trends
in Applications of Pure Mathematics to Mechanics, Vol. II (Second Sympos., Kozubnik,
1977), Monogr. Stud. Math. 5 , Pitman, Boston, 1979, pp. 263–280.
[88] J.-J. Moreau, Standard inelastic shocks and the dynamics of unilateral constraints, in Unilateral Problems in Structural Mechanics, G. del Piero and F. Maceri, eds., C.I.S.M.
Courses and Lectures 288, Springer-Verlag, Vienna, New York, 1985, pp. 173–221.
[89] J.-J. Moreau, Une formulation du contact a frottement sec; application au calcul numérique,
C. R. Acad. Sci. Paris Sér. II, 302 (1986), pp. 799–801.
[90] J.-J. Moreau, Bounded variation in time, in Topics in Nonsmooth Mechanics, Birkhäuser,
Basel, Boston, 1988, pp. 1–74.
[91] J.-J. Moreau, Unilateral contact and dry friction in finite freedom dynamics, in Nonsmooth
Mechanics and Applications, J.-J. Moreau and P. D. Panagiotopoulos, eds., International
Centre for Mechanical Sciences, Courses and Lectures 302, Springer-Verlag, Vienna, New
York, 1988, pp. 1–82.
[92] J.-J. Moreau, Numerical experiments in granular dynamics: Vibration-induced size segregation, in Contact Mechanics, M. Raous, M. Jean, and J.-J. Moreau, eds., Plenum Press,
Proceedings of the 2nd Contact Mechanics International Symposium, 1994, Carry-LeRouet, France, New York, 1995, pp. 347–158.
[93] I. P. Natanson, Theory of Functions of a Real Variable, rev. ed., Frederick Ungar, New York,
1961; original translated from the Russian by L. F. Boron and annotated by E. Hewitt.
[94] I. Newton, Philosophiae Naturalis Principia Mathematica, Reg. Soc. Phreases, London, 1686.
[95] H.-D. Niepage, Inverse stability and convergence of difference approximations for boundary
value problems for differential inclusions, Numer. Funct. Anal. Optim., 9 (1987), pp. 761–
778.
[96] H.-D. Niepage and W. Wendt, On the discrete convergence of multistep methods for differential inclusions, Numer. Funct. Anal. Optim., 9 (1987), pp. 591–617.
[97] P. Painlevé, Sur le lois du frottement de glissemment, C. R. Acad. Sci. Paris, 121 (1895),
pp. 112–115. Articles under the same title appeared in C. R. Acad. Sci. Paris, 141 (1905),
pp. 401–405 and 546–552.
[98] J.-S. Pang, A B-differentiable equation-based, globally and locally quadratically convergent
algorithm for nonlinear programs, complementarity and variational inequality problems,
Math. Programming, 51 (1991), pp. 101–131.
[99] J.-S. Pang, Complementarity problems, in Handbook of Global Optimization, Noncovex Optim. Appl. 2, Kluwer Academic Publishers, Dordrecht, the Netherlands, 1995, pp. 271–
338, in [54].
[100] J.-S. Pang and L. Q. Qi, Nonsmooth equations: Motivation and algorithms, SIAM J. Optim.,
3 (1993), pp. 443–465.
[101] J.-S. Pang and D. E. Stewart, A unified approach to frictional contact problems, Internat.
J. Engrg. Sci., 37 (1999), pp. 1747–1768.
[102] J.-S. Pang and J. C. Trinkle, Complementarity formulations and existence of solutions of
dynamic multi-rigid-body contact problems with Coulomb friction, Math. Programming,
73 (1996), pp. 199–226.
[103] J.-S. Pang, J. C. Trinkle, and G. Lo, A complementarity approach to a quasistatic multirigid–body contact problem, Comput. Optim. Appl., 5 (1996), pp. 139–154.
[104] J.-S. Pang and J. A. Tzitzouris, private communication, manuscript in preparation.
[105] L. Paoli and M. Schatzman, Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales: Cas avec perte d’énergie, Math. Model. Numer. Anal., 27 (1993),
pp. 673–717.
[106] L. Paoli and M. Schatzman, Vibrations avec contraintes unilatérales et perte d’énergie aux
impacts, en dimension finie, C. R. Acad. Sci. Paris Sér. I, 317 (1993), pp. 97–101.
[107] M. F. O. S. Pereira and J. A. C. Ambrósio, eds., Computational dynamics in multibody
systems, Papers from the NATO Advanced Study Institute on Computer Aided Analysis
of Rigid and Flexible Mechanical Systems, Portugal, 1994, Kluwer Academic Publishers,
Dordrecht, the Netherlands, 1995.
38
DAVID E. STEWART
[108] J. B. Perry, P. Henning, J. Cremer, and G. Vaněček, Jr., Contact analysis in a virtual
environment, in Proceedings of the First Workshop on Simulation and Interaction in
Virtual Environments (SIVE 95), University of Iowa, Iowa City, IA, 1995.
[109] F. Pfeiffer and C. Glocker, Multibody Dynamics with Unilateral Contacts, Wiley Ser.
Nonlinear Sci., John Wiley, New York, 1996.
[110] S. D. Poisson, Traité de Mécanique, 2nd ed., Bachelier, Paris, 1833.
[111] F. A. Potra, Runge-Kutta integrators for multibody dynamics, Mech. Structures Mach., 23
(1995), pp. 181–197.
[112] L. Q. Qi and J. Sun, A nonsmooth version of Newton’s method, Math. Programming, 58
(1993), pp. 353–367.
[113] S. Reich, Symplectic integrators for systems of rigid bodies, in Integration Algorithms and
Classical Mechanics, AMS, Providence, RI, 1996, pp. 181–191.
[114] E. J. Routh, Dynamics of a System of Rigid Bodies, MacMillan, London, 1860.
[115] A. L. Ruina, Constitutive relations for frictional slip, in Mechanics of Geomaterials: Rocks,
Concrete, Soils, Z. P. Bažant, ed., John Wiley, Chichester, New York, Brisbane, Toronto,
Singapore, 1985, pp. 169–188.
[116] J. M. Sanz-Serna, Symplectic integrators for Hamiltonian problems: An overview, in Acta
Numerica 1992, A. Iserles, ed., Cambridge University Press, Cambridge, UK, 1992,
pp. 273–279.
[117] M. Schatzman, Sue une classe de problémes hyperboliques non linéaires, C. R. Acad. Sci.
Paris Sér. A-B, 277 (1973), pp. A671–A674.
[118] M. Schatzman, Le système différentiel (d2 u/dt2 ) + ∂ϕ(u) ∋ f avec conditions initiales, C.
R. Acad. Sci. Paris Sér. A-B, 284 (1977), pp. A603–A606.
[119] M. Schatzman, A class of nonlinear differential equations of second order in time, Nonlinear
Anal., 2 (1978), pp. 355–373.
[120] M. Schatzman, A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle, J. Math. Anal. Appl., 73 (1980), pp. 138–191.
[121] M. Schatzman and M. Bercovier, Numerical approximation of a wave equation with unilateral constraints, Math. Comp., 53 (1989), pp. 55–79.
[122] W. Schiehlen, Multibody system dynamics: Roots and perspectives, Multibody Syst. Dyn.,
1, Portugal, (1997), pp. 149–188.
[123] M. F. O. Seabra Pereira and J. A. C. Ambrósio, eds., Computational methods for rigid
and flexible mechanical systems, Papers from the NATO Advanced Study Institute on
Computer Aided Analysis of Rigid and Flexible Mechanical Systems, Tróia, Portugal,
1993, Nonlinear Dynam. 9, Kluwer Academic Publishers, Dordrecht, the Netherlands,
1996.
[124] A. A. Shabana, Flexible multibody dynamics: Review of past and recent developments, Multibody Syst. Dyn., 1 (1997), pp. 189–222.
[125] P. Shi, The restitution coefficient of a linear elastic rod, Math. Comput. Modelling, 28 (1998),
pp. 427–435.
[126] D. E. Stewart, A high accuracy method for solving ODEs with discontinuous right-hand
side, Numer. Math., 58 (1990), pp. 299–328.
[127] D. E. Stewart, Numerical methods for friction problems with multiple contacts, J. Austral.
Math. Soc. Ser. B, 37 (1996), pp. 288–308.
[128] D. E. Stewart, Existence of solutions to rigid body dynamics and the paradoxes of Painlevé,
C. R. Acad. Sci. Paris Sér. I Math., 325 (1997), pp. 689–693.
[129] D. E. Stewart, Convergence of a time-stepping scheme for rigid body dynamics and resolution of Painlevé’s problems, Arch. Rational Mech. Anal., 145 (1998), pp. 215–260.
[130] D. E. Stewart, Reformulations of measure differential inclusions and their closed
graph property, J. Differential Equations, submitted; also available online from
http:/www.math.uiowa.edu/faculty/comp-math-1998.htm.
[131] D. E. Stewart and J. C. Trinkle, An implicit time-stepping scheme for rigid body dynamics
with inelastic collisions and Coulomb friction, Internat. J. Numer. Methods Engrg., 39
(1996), pp. 2673–2691.
[132] D. E. Stewart and J. C. Trinkle, An implicit time-stepping scheme for rigid body dynamics
with inelastic collisions and Coulomb friction, Internat. J. Numer. Methods Engrg., 39
(1996), pp. 2673–2691.
[133] D. Stoianovici and Y. Hurmuzlu, A critical study of the applicability of rigid-body collision
theory, ASME J. Appl. Mech., 63 (1996), pp. 307–316.
[134] K. R. Stromberg, An Introduction to Classical Real Analysis, Wadsworth, Belmont, CA,
1981.
RIGID-BODY DYNAMICS WITH FRICTION AND IMPACT
39
[135] W. J. Stronge, Rigid body collisions with friction, Proc. Roy. Soc. London Ser. A, 431 (1990),
pp. 168–181.
[136] K. Taubert, Differenzverfahren für gewöhnliche Anfangswertaufgaben mit unstetiger rechte
Seite, in Numerische Behandlung Nichtlineare Integrodifferential- und Differentialgleichungen, R. Ansorge and W. Törnig, eds., Lecture Notes in Mathematics 395, SpringerVerlag, Berlin, New York, 1974, pp. 137–148.
[137] K. Taubert, Differenzverfahren für Schwingungen mit trockener und zäher Reibung und für
Regelungsysteme, Numer. Math., 26 (1976), pp. 379–395.
[138] K. Taubert, Converging multistep methods for initial value problems involving multivalued
maps, Computing, 27 (1981), pp. 123–136.
[139] J. C. Trinkle, J.-S. Pang, S. Sudarsky, and G. Lo, On dynamic multi-rigid-body contact
problems with Coulomb friction, Z. Angew. Math. Mech., 77 (1997), pp. 267–279.
[140] L. T. Watson, S. C. Billups, and A. P. Morgan, Algorithm 652: HOMPACK: A suite
of codes for globally convergent homotopy algorithms, ACM Trans. Math. Software, 13
(1987), pp. 281–310.
[141] J. M. Wendlandt and J. E. Marsden, Mechanical integrators derived from a discrete variational principle, Phys. D, 106 (1997), pp. 223–246.
[142] J. Yen, Constrained equations of motion in multibody dynamics as ODEs on manifolds, SIAM
J. Numer. Anal., 30 (1993), pp. 553–568.
[143] J. Yen and L. R. Petzold, An efficient Newton-type iteration for the numerical solution
of highly oscillatory constrained multibody dynamic systems, SIAM J. Sci. Comput., 19
(1998), pp. 1513–1534.