BIB-VERSION:: CS-TR-v2.0
ID:: CORNELLCS//TR92-1275
ENTRY:: 1993-10-14
ORGANIZATION:: Cornell University, Computer Science Department
LANGUAGE:: English
TITLE:: Dynamic Simulation of Non-Penetrating Rigid Bodies
AUTHOR:: Baraff, David 
DATE:: March 1992
PAGES:: 205
COPYRIGHT:: David Baraff 1992 - All Rights Reserved
ABSTRACT::
This thesis examines the problems and difficulties in the forward dynamic 
simulation of rigid bodies subject to non-penetration constraints. By adopting 
a simple but well-defined model of rigid body dynamics, we are able to focus 
on and gain insight into some of the inherent difficulties of rigid body 
simulation. Additionally, computationally practical solutions to some of the 
problems encountered in this thesis are presented.

Enforcing non-penetration constraints is essentially a two step process. The 
first step of the simulation involves the detection of potential contacts 
between bodies. This thesis presents collision detection algorithms for the 
dynamic simulation of bodies that are composed of both polyhedra and convex 
closed curved surfaces. The collision detection algorithms exploit temporal 
coherence to achieve fast running times and are a practical solution to the 
problem of collision detection during simulation.

The second step of the simulation involves the computation of the contact 
forces between bodies that maintain the non-penetration constraint. This 
thesis considers first the problem of computing contact forces between a pair 
of bodies that contact at a point without friction. A mathematical formulation 
for the contact force between the bodies is presented, and then modified to 
yield a formulation that is computationally practical for use in a simulator. 
After considering the dynamics of single point contacts, systems with multiple 
 contacts are considered both in terms of computational complexity measures 
and practical solution methods. The methods used in this thesis to compute 
constraint forces are also theoretically and practically compared with a 
popular method for preventing inter-penetration called the ``penalty method''.

After considering frictionless systems, this thesis considers systems of 
bodies that behave according to the classical Coulomb model of friction (which 
includes both sliding and dry friction). This leads us to consider systems in 
which there are no solutions to the classical constraint force equations, as 
well as systems which admit multiple solutions for the constraint force 
equations and whose subsequent behavior is thus indeterminate. Both 
computational and practical complexity results for simulating such systems are 
discussed.
END:: CORNELLCS//TR92-1275
BODY::
Dynamic Simulation of
Non-Penetrating Rigid Bodies
David Baraft
Ph.D Thesis
92-1275
March 1992
Department of Computer Science
Cornell University
Ithaca, NY 14853-7501
DYNAMIC SIMULATION OF?
NON-PENETIQATING R?IGID BODIES
A Dissertation
Presented to the F'aculty of the Graduate School
of Cornell University
in Partial F'ulfillment of the Requirements for the Degree of
Doctor of Philosophy
by
Da?rid Baraff
May 1992
o David Baraff 1992
ALL RIGHTS RESER\7ED
DYNAMIC SIMULATION OF NON-PENETRATING RIGID BODIES
David Baraff, Pli.D.
Cornell University 1992
This thesis exan?ines the proble?is and difficulties iii tile for?vard dynan?ic sin?ulation
of rigid bodies subject to non-penetration constraints By adopting a sin?ple but
?ell--Hdefined n?odel of rigid body dynan?ics, ?ve are able to focus on and gaiii
insight into son?e of the inherent difficulties of rigid body sin?ilation. Additionalky,
coinputationally practical solutions to son?e of the prol?len?s encountered in this
thesis are presented.
Enforcing non-penetration constraints is essentially a t??o step process. The first
step of the sin?ulation iiivolves the detection of potential contacts bet??en bodies.
This thesis presents collision (letection algoritlin?s for tlie dyiiainic sin?ulation of
bodies that are con?posed of 1)0th polyliedra and convex closed curved surfaces. Tlie
collision detection algoritlin?s exploit ten?poral coherence to achieve fast running
times aiid are a practical soluth?n to the problen? of collision detection during
sin?ulation.
The second step of the sin?ilation involves the con?putation of the contact forces
bet?veen bodies that n?aintain tlie non-penetration constraint. This thesis coiisiders
first tlie problen? of computing contact forces between a pair of bodies that contact
at a point without friction. A n?athen?atical forn?ilation fi?r the contact force
betweeii the bodies is presented, and then n?odified to yield a fi?rn??ilation that is
con?p'itationally practical for use in a sii?ulator. After considering the dynan?ics
of single point contacts, systen?s with multiple contacts arc considered both in
terms of computational complexity measures and practical solution methods. The
methods used in this thesis to compute constraint forces are also theoretically and
practically compared with a popular method for preventing inter-penetration called
the "penalty method".
After considering frictionless systems, this thesis considers systems of bodies that
behave according to the classical Coulon?? model of friction (which includes both
sliding and dry fliction). This leads us to consider systems in which there are no
soliitioiis to tlie classical constraint force equations, as well as systems which admit
multiple solutions for the constraint force equations and whose subseq'ient behavior
is tli'is iiideterminate. Botii computational and practical com??lexity results fi?r
sin?ilating such systems are discussed.
Table of Contents
2
3
4
5
Introduction
1.1 Thesis concerns and restrictions
1.2 Thesis overvie?v
A Simple Dynamics Problem with Non-penetration Constraints
2.1			Unconstrained i?otion .
2.2			Equality constrained n?otion . . .			.			. .			. .			. .
2.3			Inequality constrained n?otion
2.4			The penalty method . . . . . .			.			.			. .			. .
2.5			Exact			solution methods
2.5.1			The equality constrained problem			. .			.			. .			.
2.5.2			The inequality constrained			problem
Simulation Overview
3.1 Updating the current state
3.2			Collision detection .			.
3.3 Contact point determination
3.4			ImpuLse computation			. . .
3.5 Constraint/friction force determination
Summary of Previous Work
4.1 Collision/contact determination
4.2 Impulse computation
4.3 Penalty methods for constraint
4.4 Exact methods for constraint forces
4.5			Friction . .			. .			. .
Collision/Contact Determination
5.1 Convex polyhedra .
5.2 Convex curved surfaces .
forces
vii
2
4
5
5
(3
7
8
11
11
15
21
22
24
2(3
27
28
31
31
33
35
36
37
39
41
45
5.3
5.2.1 Extrei?al points of curved surfaces
5.2 2 Polyhedron/curved surface determination
Bounding boxes
5.3.1 The one-dimensional case
5.3.2 The three-dimensional case
Local Motion Constraints
6.1 Configuration Space
6.2 Projections onto S . .
6.3 Calculating a2c/ap2
6.4 Computing C?(r?(t0)) directly
6.5 Predicting collision times
6.6 Rolling contact
6.7 Ill-conditioned Jacobians
The Penalty Method
7.1 Equality constrained motion
7.1.1 The physically correct motion
7.1.2 Niotion due to the penalty method
7.2 Inequality constrained motion
Frictionless Systems
8.1 Existence and uniqueness of sohitious
8.2 Implementation issues
The Coulomb Friction Model
9.1 Coordinate geometry at a contact point
9.2 Dynamic friction conditions
9.3 Static friction conditions . . .
9.4 Impulsive flictional forces
6
7
8
9
10 Dynamic Friction
10.1 Indeterminacy
10.2 Inconsistency
11 An NP-complete Class of Configurations
12
Physical Models
12.1 A model of inconsistency . . .
12.2 A model of indeterminacy
12.3 Indeterminacy and inconsistency con??ined
46
51
53
o+ .			53
56
59
60
62
65
70
75
76
77
85
86
87
89
94
97
100
101
105
106
107
o+ .			.			108
108
111
o+ . . . 112
118
121
131
132
134
137
viii
12.4 Preferred sointions
138
13
14
Computing Valid Contact Forces and Impulses
13.1			Defining valid contact impnlses . .			.
13.2			Lemke's algoritlim
13.3			Continuity of solutions . . .			. .			.
Static Friction
14.1			Complexity of static friction . .			.
14.2			Approximation methods
14.2.1			L5tstedt's approximation . .			.
14.2.2			Dynamic friction approximation
14.2.3 Quadratic programming approximation
15 Future Directions
A Extremal Point Conditions for Parametric Surfaces
B Termination Conditions of Lemke's Algorithm
C Iterative Solution Methods for Static Friction
D Simulation Examples
Bibliography
ix
143
143
145
148
149
149
157
157
160
162
165
171
175
181
187
201
List of Figures
2.1			An equality constraint between a particle and a surface, .			. .			.			7
2.2			An inequality constraint between a particle and a surface.			8
2.3			Constraint forces parallel to the surface norn?al			VC(.T(i)) for an
equality constrained problen?.			.			13
2.4			Constraint forces parallel to the surface norn-ial			VC(x(t)) for an
inequality constrained problen?. .			.			.			16
3.1			The coi?putational steps required to con?pute w<i't?(?')			. 23
5.1 The i?inin?un? distance points I3a and Pb ??et????en two convex surfaces. 47
5.2			lifter-penetration between surfaces. 			. 						. 48
5.3			The deflnitk?n of the extreu?al points						49
5.4			Geon?etric conditions for the extren?al points			. 						50
5.5			The sweep/sort algorithin			54
5.6			A coherenc&based n?ethod of detecting			55
6.1
6.2
6.3
6.4
6.5
6.6
6.7
6.8
10.1
10.2
10.3
overlaps			.
Three configurations in world space and their representation as
points in C-space.
The reaction force Freadion			. . .			. . .
Change of the extrenial points as a function of the configuration.
C expressed in tern?s of the extren?al points. . . . . .
Side view of the implicit surfaces F and 0 expressed explicitly by
flinctions f and g. .			.
Frictionless oscillation of an ellipsoid A on a plane I?
Frictionless oscillation of a n?ore eccentric ellipsoid A on a plane B.
Frictionless oscillation of a superquadric ellipsoid A on a plane B.
A one contact point configuration between a thin rod A and a fixed
base B. .			. .			. . .			. . .			. .
An indeteri?inate configuration.
A potentially inconsistent configuration.			. . .			. . .			. .
xi
61
64
66
67
71
`9
81
82
114
117
118
11.1
11.2
11.3
12.1
12.2
12.3
12.4
12.5
14.1
14.2
14.3
A possibly inconsistent configuration.
A configuration whose consistency depends on an external force.
A configuration that is consistent if and only if the friction forces on
B sum to jiS. .			. . .			. . .			. .			. . .			.
A model of indeterminacy.
A configuration with no initial inter-penetration. .
A configuration with an initial inter-penetration depth of (11
A configuration with a single ii-npulse-free solution.
A configuration with two different non-impulsive solutions at ti'?e
zero.
Force space
Line N intersects L and 1?L.
Line iV intersects RL, but not L.
15.1 Non-polygonal contact region.
D.1
D.2
D.3
Jack falling down stairs.
Falling superquadric dice.
Collapsing structure.
xii
124
125
126
133
135
136
138
141
152
153
154
167
189
193
197
Chapter 1
Introduction
The research problem motivating this thesis has a long and detailed history, having
been a driving force behind the development of both continuous mathematics and the
high speed digital computer. The problem in q'lestion is the problem of dynam?cs:
the study of the motions of bodies. Calculus, for instance was invented l?y Newton,
as he pondered the motion of bodies under the influence of gravity. The first
high speed (at the time) digital comp?lter, the Eniac, was built to solve ballistics
problems.
The dynamics problems originally posed by Newton were all of an unconstrciined
nature; for example, the motion of two bodies subject only to gravitational forces
between them. In contrast, the dynamics problems considered in this thesis are
all of a constrained nature. In particular, this thesis will consider the problem
of simulating the dynamics of rigid bodies subject to non-penetration constraints
that is, constraints that require bodies to behave as if they were solid, and avoid
inter-penetration.
In order for a simulation to be computable, it is necessary to make simplifying
assumptions about the environment being simulated. In order to make sim?ilation
2
practical, we usually try to find the simplest assumptions that adequately model a
particular environment. As the requirements on our simulation are increased (more
complex phenomenon, greater accuracy), we must also increase the complexity of
our sin?ilation model. It is naturally assumed, therefore, that realistic simulation
is hard mainly because it requires a very complex simulation model. However, as
this thesis will show, even very simple simulation models yield complex results and
questions.
1.1 Thesis concerns and restrictions
Iii this thesis, we will study the problems involved in simulating tlie niotions of bodies
in accordance with relatively simple physical models of bodies' dynaniics. Empirical
correctness is not a major concern of this thesis. ?7e do not consider the dynamics
models used in this thesis to be empirically correct. That is not to say that our
dynamics models are ad hoc, or deliberately wrong a more apt description would
be to say that our "simple" dynamics models are incomplete descriptions of more
"complicated" dyiiamics models (which are themselves considered simplifications of
even more coniplex models). Given these simple models, our interest lies in the
ramifications of using these models as a basis for simulation. The vie??poiiit taken
by this thesis then is not unlike the viewpoint of rational mechanics[25]. Rational
niechaiiics involves logically deducing the dynamics of systems based on a set of
axioms. Rational mechanics gives a "correct result" in as much as the axiom set is
valid (with respect to the real world) and complete.
Unlike rational mechanics, this thesis will also be concerned with the practical
and computational difficulties that arise in trying to perform simulations based on
our simple models. If in fact our simple models are "subsets" of more complicated
models, then any findings concerning the difficulties of simulations that use our
3
simple models should serve as a rough lower bound on the simulation difficulties
that will be encountered when more complex models are used. Given that the
complex models will yield additional difficulties, it is sensible to first discover the
problems stemming from the simple models before moving on to work with the more
complex models.
A secondary concern of this thesis is offering practical solutions to the problems
encountered in actually using our simplified dynamics model for sin?ilation. For
some applications (for example, computer animation, virtual reality, or robot motion
planning research) a sin?ilator based on the simple dynamics model considered in
this thesis is an end-goal in itself. In these cases, it is important to have practical
and efficient simulation algorithms at hand. At present, the vast majority of work in
constrained dynamics involving non-penetration constraints uses a simple approxi-
mation technique that does not completely enforce the non-penetration constraints
(even discounting numerical tolerances). In this thesis, techniques that completely
enforce non-penetration constraints (within mimerical tolerances) are used. These
teA?niques yield significantly better results and have a firmer theoretical basis than
the simple approximation methods, b'it are considerably more complicated.
The bodies considered in this thesis will be restricted to be perfectly rigid bodies.
This implies that bodies propagate forces instantly, without delay, and that complex
"micromechanical" processes (such as tlie elastic response during a collision) are
replaced with simple "macromechanical" results. A more iniportant restriction
concerns tiie modeling of contact between bodies. Contact between bodies will
always be subject to the restriction that it be modeled as occurring at some finite
number of points between bodies. Friction, when considered, will be restricted to tlie
Coulomb model of friction. As we shall see, these few simple laws yield simulation
problems of considerable depth and complexity.
4
1.2 Thesis overview
Chapter 2 presents the prototypical problem that characterizes the problem of
dynamic simulation with non-penetration constraints.
Chapter 3 gives an overview of the simulation process, in terms of the major
computational steps required.
Chapter 4 summarizes previous work, as classified by the major computational
steps of simulation defined in Chapter 3.
Chapter 5 presents collision detection algorithms that are well-suited for the type
of environment encountered in dynamics simulation.
Chapter 6 develops the local motion constraint that arises between contacting
bodies to prevent inter-penetration.
Chapter 7 analyzes the penalty method, which is currently the most popular
method for preventing inter-penetration during simulation.
Chapter 8 formulates equations to solve for constraint forces for frictionless
systems.
Chapters 9 and 10 introduce a friction model and discuss tlie coniplications that
arise fiom its use. Chapter 11 discusses the computational complexity issues that
arise fiom the introduction of fliction.
Chapter 12 proposes a physical interpretation of the compiitational difficulties
arising from friction. Chapters 13 and 14 discuss principles and computational
methods for simulations with friction.
Chapter 15 concludes this thesis by discussing future directions for siumlations
with non-penetration constraints.
Chapter 2
A Simple Dynamics Problem
with Non-penetration Constraints
In considering the steps needed to solve dynamics problems with non-penetration
constraints, it is helpful to have a dynamics problem with non-penetration con-
straints firmly in mind. The simple dynamics problem presented in this diapter
illustrates the important issues of dynamic sin?ilation involving non-penetration
constraints. The problem is to simulate the dynamics of a point mass particle in
three-space (R3). The position r of the particle at tilne t in R3 is written r(t).
The velocity v of the particle at time t is v(t) , and the particle is assumed to have
unit mass. An arbitrary external force F(t) acts on the particle at tiine t.
2.1 Unconstrained motion
When there are no constraints on the particle's motion (that is, no obstacles for
the particle to encounter), the motion of the particle is simple. Using Newton's
5
6
F = ma, the differential equation
= F(t)			(2--H1)
describes the particle's motion. In numerical simulation though, this second order
differential equation is converted to a coupled first order differential equation by
writing
or in vector notation,
= v(t)
= F(t)
(I			x(t)			__			v(t)			(2--H2)
Wt			v(t)			F(t)
The variables x and v are called state variables, and the collection of x and v is
called the state of the system. The variable x is a spatial variable while v is a
velocity variable.
2.2 Equality constrained motion
Now 5UPP05C that S is a surface in R3, and supposed the particle is constrained
to always lie on S. Let S be modeled implicitly by a scalar function c; that is, a
point p lies on s if and only if C(p) = 0. Using C, the constraint on the particle
can be written
C(x(t)) = 0.			(2--H3)
Since this constraint can be written as an equality, this is an example of an equality
eonstrained dynamics problem (Figure 2.1).
7
Figure 2.1: An equality constraint between a particle and a surfae. The particle,
constrained by C(x(t)) = 0, must always lie on the surface S
2.3 Inequality constrained motion
Now suppose that the constraint of the previous problem is relaxed, so that the
particle is constrained to lie either on or "above" S (Figure 2.2). If points p lying
above S satisfy C(p) > 0, then the constraint that the particle lie on or above S
can be written
> 0.			(2--H4)
Since the constraint requires an inequality, this is an example of an inequality
constrained dynamics problem.
The inequality constrained dynamics problem is a superset of the unconstrained
and equality constrained dynamics problems, and is thus the most difficult problem
of the three. Clearly, the inequality constrained dynamics problem subsumes the
unconstrained dynamics problem; if at time t the particle satisfies C(x(t)) > 0
then for some period of ti'?e following t,the motion of the particle can be regarded
as unconstrained. Similarly, the equality constrained problem is a special case of
8
Figure 2.2: An inequality constraint between a particle and a surface. The particle,
constrained by C(x(t)) > 0, can move either on or "above" the surface 5
the inequality constrained problem the equality constraint C(x(t)) 0 can be
replaced by the constrah?ts C(.r(t)) > 0 and --HC(x(t)) > 0.
2.4 The penalty method
Constrained dynamics problems can be very difficult to solve. A very common
method for solving a problem with constraints is the penalty method. The penalty
method for constrained dynamics problems is based on a mimerical sol?ition method
for constrained optimization problems. In both fields, the penalty method converts
a constrained problem to an unconstrained problem where deviation fiom the
constraint is penalized; that is, in the new problem, satisfiiction of the constraint is
encouraged, but not strictly enforced.
In optimization, a typical equality constrained problem such as
minimize f(z) such that			g(z) = 0
(2--H5)
9
can be rewritten as an unconstrained problem
minimize f(z) + kg(z)2 as k			00.			(2--H6)
The term kg(z)2 is called the penalty function. The idea is that as k grows larger,
potential solutions for z n?ist make g(z)2 smaller, to minimize Equation (2--H6). In
the limit, as k goes to infinity, the solution of Equation (2--H6) must satisfy g(z) = 0
while minimizing f(z). In practice, z is obtained by solving Equation (2--H6) for
a series of increasing values of k until the series of sollitions converges (within
numerical tolerance) to a limit. ?lthough the method has a theoretically firm basis,
in practice, it is not a very robust numerical method. The main problem is that
as k grows, Equation (2--H6) can become very poorly conditioned arid difficult to
solve[14]. The main attraction of the method is that it a very simple way of turning
a constrained problem into an unconstrained problem.
The equality constrained particle/surface problem can be coiiverted to an un-
constrained dynamics problem in a similar fashion. Tlie idea is that the particle is
not explicitly constrained to lie on the surface. Ho?vever, if the particle drifts off tlie
surface, a penalty force acts on the particle, in a direction normal to tlie surface, so
as to p?ill the particle back to the surface. Typically, the penalty force is modeled
as a linear spring force; that is, the penalty force pulls tlie particle back towards
the surface with a a strength equal to some constant k times the distance of tlie
particle from the surface. If we let p(x(t)) denote tlie closest point on the surface
to x(t), then the penalty force is1
--H k (x(t) --H
If the surface is a plane, the closest point 011 the plane to the particle is unique. Otherwise,
the closest point on the surface to the particle is unique as long as the particle remains in some
neighborhood containing the surface; the extent of this neighborhood depends on the curvature
of the surface. Since the penalty method is intended to keel) the particle close to the surface, the
closest point on the surface to the particle is always assumed to be well-defined.
10
Note that as x(t) drifts farther from the surface, the force --Hk (x(t) --H p(r(t))) grows
larger in magnitude. The introduction of the penalty force reduces the dynamics
problem to simply an evaluation of the ordinary differential equation
d			x(t)			__			v(t)			(2--H7)
dt			v(t)			--Hk (r(t) --H			t F(t)
If the penalty method for dynamics were to completely en?ilate the penalty
method for constrained optimization, the simulation would be repeated with in-
creasing values of k until the behavior of the particle approached a limit. However,
the penalty method, as used by dynamics, chooses a single value for k. If the value
chosen for k is too small, it may do an inadequate job of enforcing the constraint; in
this case, the simulation would be rerun with a higher value of k. Unfortunately, as
in the optimization problem, ill-conditioning occurs as k grows larger, in the guise
of "stiffness" for the differential equation of motion. Stiff differential equations are
expensive and difficult to solve.
The penalty method can also be used for the inequality constrained parti-
cle/surface problem. The modification is straightforward; whenever the particle lies
below the surface, a penalty force acts upwards on the particle, as in the equality
constrained case. However, when the particle lies above the s?irf?we, no penalty force
acts on the particle.
For neither the equality constrained prol?lem nor the inequality constrained
problem is it obvious that the penalty method must give a correct solution, given
a sufficiently large value of k. It is clear however that except for a limited class
of simulations, the penalty method is an unsuitable method, because of the ill-
conditioning of the differential equations of motion. Further analysis of the penalty
method will be deferred to Chapter 7.
11
2.5 Exact solution methods
The penalty method of the previous section involves a formulation of the problem in
which constraints are only approximately satisfied. A completely different viewpoint
is taken iii this section, in that constraints are required to be satisfied "exactly"
(within nui?erical tolerances).
2.5.1 The equality constrained problem
The equality constrained problem of the particle that must remain in contact with
the surface S can be solved in two ways. The first solution method involves
converting the constrained problem to an unconstrained problem essentially the
method of Lagrangian dynamics. This can be done if a parametric representation
for S can be found.
Otherwise, the problem can be solved by introducing a constraint force Fc(t)
into the problem. The role of the constraint force is to act on the particle so that
it remains in contact with the surface. Additionally, the constraint force must be
workiess; that is, the work done on the particle by the constraint force must always
be zero. For the particle/surface example, this requirement forces Fc(t) to act on the
particle in a direction normal to tlie surface 5 at x(t). Thus, the constraint force's
responsibility is to adjust the acceleration of the particle normal to the snrfa?e, so
that the particle will remain on the surface. The normal to 5 at the point x(t) is
12
the vector VC(x(t)) where2
VC(x(t)) =
ac
Ox
a?; (x(t))
OaCz (x(t))
Then Fc(i) can be written Fc(t) = A(t)VC(x(t)) where A(t) is an unknown scalar.
The scalar A(t) can be either positive or negative. In Figure 2.3a, the constraint
force Fc(t) acts to oppose a large external force F(t). The constraint force acts in
the same direction as VC(x(t)); in this case, A(t) > 0. In Figure 2.3b, the external
force F(t) now pushes the particle away from the surface in the VC(x([)) direction.
In this case, the constrah?t force must act opposite VC((t)); hence A(t) ? 0.
The equation of motion for the particle then becomes
= F(t) + Fc(t) = F(t) + A(t)VC(x(t)).
(2--H8)
How is A(t) computed?
In equality constrained dynamics problems, it is assumed that the constraint is
initially met; that is, at time zero
Additionally,
C(x(O)) = 0.			(2 9)
= ?tC(x(0)) = VC(x(0)) x(O)			(2--H10)
must also be zero as well. If it is not, then the particle has a velocity that is
not completely tangential to 5 at the point x(0); this situation will require an
impuThive constraint force, and will be dealt with shortly. If however both C(x(t))
2 ?? is assumed that the function C modeling S is chosen sud? that VC(p) does not vanish
on
13
(a)
?t)VC (x (t))
F(t)
?t)VC (x (t))
Figure 2.3: Constraint forces parallel to the surface norn?al VC(.T(t)) for an equality
constrained problem. (a) To maintain the particle/surface constraint, A(1) > 0 so
that the constraint force Fc(t) = A(t)VC(x(t)) pushes the particle in the positive
VC(x(t)) direction. (b) The constraint force Fc(t) n?ist be applied opposite to
""`C(x(t)); thus A(t) < 0.
14
and C(x(t)) are zero for t = 0, then both quantities will remain zero as long
as C(x(t)) = d$2 C(x(1)) is zero. Intuitively, C(r(t)) is a measure of the particle's
displacement normal to S, and C(x(t)) is a measure of the particle's velocity normal
to S. Thus, 6(x(t)) is a measure of the particle's acceleration normal to S. The
constraint force A(t)VC(x(t)) is d?osen so that C?(x(t)) is always zero.
Starting from C(x(t)) = ?C(x(t)) x?(t) and taking its denvative
= d
dt (VC(x(t))
--H			ftVC(x(t))			x(t) + ""`C(?(t))			ydt?(t)			(2--H11)
--H x(t) (?2C(x(t))5?(t)) + VC(x(t))
where the term V2C(x(i)) is the matrix of second partial denvatives of C (that
is, the curvature of the surface), evaluated at the point x(t). Using Equa-
tion (2 8),
=			(i(t)			(v2C(x(t))?(t))) + (vC(x(t))
--H			(x(t)			(v2c('?(t))i(t))) + (vC(x(t)) (F(t) +
--H			i'(t)			(v2C((t))J?([)) + VC(x(t)) F(t)			(2--H12)
+ A(t)VC(x(i))
Setting C?(?(t)) to zero and solving for A(t),
A(t) = i([)(v2C(x(t))x(t))+VC(x(i))F(t) (2--H13)
VC(x(t)) VC(x(t))
Having found A(t), Equation (2--H8) can be solved numerically, and the motion of
the particle will satisfy the equality constraint.
15
2.5.2 The inequality constrained problem
The solution to the inequality case is only slightly more complicated for this simple
example. Suppose that initially the particle lies above the surface 5, and does not
come into contact with the surface until time ti. During the time interval zero to
t1 the constraint is said to be inactive and the motion of the particle satisfies
The first problem then is to determine t1 , that is, to determine when the particle
contacts the surface. This is the problem of collision detection.
After the particle contacts the surface at time t1, the constraint C(w(t)) > 0
becomes active. To prevent the particle from moving below the surface, a constraint
force must act. In what follows, it is assumed that ?C(x(t)) points "?ipwards"
(towards the region of space that the particle is permitted to move in) as in
Figure 2.4. As in the equality constrained problem, the constraint force is assumed
to do no work on the particle; so again, the constraint force can be written as
Fc(t) = A(t)VC(x(t)). Unlike the equality constrained problem, the constraint
force must be directed away from the surface, in the VC(;r(t)) direction. This
means that A(t) n?ist satisfy A(t) > 0. In Figure 2.4a, where the external force
F(t) pushes the particle dowuwards, the constraint force A?C(i(t)) is exactly the
same as in Figure 2.3a. However, when the external force F(t) of Fig?ire 2.41) acts
on the particle, no constraint force acts (that is, A(t) is zero), and the particle
breaks contact with 5, accelerating upwards. If A(t) was allowed to be negative,
as in Figure 2.3b, the particle, once in contact with the surface, could be held to
the surface for all time by computing a value for A(t) as in the equality constrained
problem.
16
(a)
?t)VC (x (t))
F(t)
F(t)
Figure 2.4: Constraint forces parallel to the surface normal ""`C(x(t)) for an
inequality constrained problem. (a) A constraint force Fc(t) = A(t)VC(x(t)) with
A(t) > 0 acts to prevent the particle from accelerating downwards. (b) The external
force F(t) accelerates the particle upwards, and contact is broken. No constraint
force acts on the particle.
17
First derivative considerations
If at time t1 (when the particle first contacts the surface)
C(x(t1)) = VC(r(t1)) x'(ti) = VC(x(11)) v(ti) < 0, (2--H14)
then the particle is said to have collided with the surface, To prevent the particle
from moving below the surfac?e, the velocity ?)(t) must undergo a discontimiity
at time ti; that is the only means of completely enforcing the constraint. The
discontinuity is said to be the result of an impulsive constraint force that acts normal
to the surface at the point x(11). An impulsive force has the dimensions of force
times time, or equivalently mass times velocity, and produces a discontinuity in the
velocity. Impulsive forces are the limiting effect of large forces acting over small
time intervals, as the force strength goes to infinity and the time interval goes to
zero. Thus, the impulsive constraint force will also be normal to the surfae, and
can be expressed as jVC(x(ti)) where j is a positive scalar
Because the particle and the surface are regarded as rigid bodies, the strength
of the impulsive force is found by using the empirical relation
= --Hev7?			(2--H15)
In this relation, V?+ is the velocity of the particle normal to the surface after the
impulsive constraint force is applied, and V? is the normal velocity prior to the
collision. The scalar 0 < c < 1 is called the coefficient of restitution and depends
on the material properties of the particle and the surface.
Since the particle has unit mass, its velocity v+(t1) after application of the
impulse jVC(x(t1)) is
v+(t1) = v(ti) + j""'C(x(t1)).			(2--H16)
(This equation can be fon?ally derived by defining v+(ti) = limAt??o v(ti + At)
when a force AVC(w(t)) acts on the particle from th?e t1 to time t1 + At, with A
18
a scalar, and lirnAt?o AAt = j.) The normal velocities V? and V?+ are
= VC(r(11)) v(t1) and V?+ = VC(x(t1)) v+(ti).
Substituting into Equation (2--i 5),
or
VC(x(t1)) v+(tj) = --HEVC(r(t1)) v(tj)
VC(x(t1)) (v(t1) t jVC(x(t?))) = --H?VC(x(t1)) v(t1).
Solving for j,
(2--H17)
(2--H18)
(2--H19)
--HEVC(r(ti))v(t1)--H?C(x(t1))v(t1)
J = VC(x(t1))
(?+1)VC(x(ti))v(ti)
VC(r(t1)) VC(x(t1))			(2--H20)
(? + 1)C(x(tj))
VC(x(t1)) VC(x(11))
Note that since c(r(t1)) < 0 and VC(r(t1)) VC(x(t1)) > 0, then j > 0. This
indicates that the constraint impulse acts upwards on the particle, as expected.
If ? > 0, then the particle will have avelocity along VC(.T(ti)) after the collision;
the particle has bounced off the surface, and is now moving away from the surface.
The constraint changes back to being inactive, and the motion of the particle is
considered unconstrained until sud? time as it contacts the surface again. Note that
this entire process happens over a zero length time interval; thus, the motion of the
particle is piecewise continuous.
Otherwise, if E = 0, then the particle has no velocity normal to the surface after
the collision. In this case, the constraint stays active, and a regular, non-imp?ilsive
constraint force can be considered.
19
Second derivative considerations
The case when C(x(t)) < 0 has just been discussed. If C(x(t)) > 0, the particle has
a velocity directed away from the surface, and the constraint immediately becomes
inactive after time t. Thus, in considering the constraint forces that prevent the
particle from moving below the surface at time t, we may assume
C(x(t)) = C(x(t)) = 0.			(2--H21)
In the equality constrained problem, C(w(t)) = 0 was maintained by choosing
A(t) sud? that C?(x(t)) = 0. For the inequality constrained problem, C(x(t)) > 0 is
maintained by choosing A(t) so that C?(?(t)) > 0, subject to the condition A(t) > 0.
Substituting from Equation (2--H12) into C?(x(t)) > 0, we obtain
= x(t) (V2C(x(t)).?([)) + VC(x(t)) F(t)
+ A(t)VC(x(t)) VC(i(t)) > 0.			(2--H22)
Since VC(x(t)) VC(x(t)) is always positive, C(i(i)) > 0 is equh?lent to
(V2C(x(t))Jt(t))+VC(x(t))F(t)
A(t) > --H			VC(x(t)) VC(x(t))			.			(2--H23)
Equation (2--H23) places only a lower bound on A(t), and does not specif,T A(t)
uniquely. However, a reasonable restriction to impose on A([) is that it be zero
whenever C?(x(t)) > 0. This is sensible because if C(r(t)) > 0 then C(x(t)) is an
increasing flinction and the particle will move away from the surface immediately
after time zero. NVhen this happens, the constraint becomes inactive and the
constraint force is non-existent; or equivalently, A(t) is zero. This restriction can be
written very simply as
A(t)C(x(t)) = 0			(2--H24)
20
whid?, along with C?(x(t)) > 0 and A(t) > 0 asserts that either C?(x(1)) > 0 and
A(t) = 0, or C?(xQ)) = 0.
?Vith this added restriction, A(t) is computed very easily. If A(t)'s lower bound
(Equation (2--H23)) is positive, then setting A(t) equal to its lower bound results in
A(t) > 0 and C?((1)) = 0; otherwise, A(t)'s lower bound is negative. Since the
denominator of Equation (2--H23) is always positive, it must be that
(V2C(x(t)):it(t)) + VC(x(t)) o+ F(t) > 0			(2--H25)
Setting A(t) = 0 satisfies A(t) > 0 and A(!)C?(x(t)) = 0. Using Equation (2--H22),
C(x(t)) satisfies
o+ (v2C('?(t))([)) + VC(x([)) o+
o+ (v2c(r(t))t(t)) + VC(x(t)) o+ (F(t) + A(i)vc(x(t)))
o+ (V2C(x(t))it(t)) + VC(x(t)) o+ F(t) > 0.
Note that durh?g the period of time that C?(x(t)) = 0, the value of A(t) is the
same as if the inequality constraint was in fact an equality constraint; during this
time, the particle's behavior is the same as if the particle had been constrained to
remain on the surface. However, when in the equality constrained problem A(t)
would have decreased below zero, in the inequality constrained problem A(t) stops
decreasing at zero and the particle breaks contact with the surface.
Chapter 3
Simulation Overview
Having described a simple, yet complete, problem with non-penetration constraints,
the overall structure of a simulation algorithm for the problem is presented, along
with the extensions needed to sin?ilate an arbitrary system of rigid bodies.
Overall, the sin?ilation process can be described as nothing more than the
numerical sollition to a first order ordinary differential equation (ODE). A first
order ODE has the fon? wdtY(t) = f(t,Y(t)). The ODE for a system of arbitranly
many rigid bodies is given in Section 3.1. The vector Y(t) is called the state of the
system, and wdty(t) = f(t,Y(t)) is called the state dernvative. Numerical methods
for solving such equations are well known. The simulator discussed in this thesis
uses an explicit solution method (as opposed to an inq?licit solution method) for
solving the ODE. (Numerical solution techniques for ordinary differential equations
are discussed by Shampine and Gordon[47J.) An explicit solution method requires
a value for Y(O) , which is called the initial state of the system, and a time tend
for which the solution method is to compute Y(tend) In computing Y(tend),the
solution method will need to be given values for the state derivative at various times
from zero to tend. The solution method is completely responsible for choosing values
21
22
of I at which to compute the state derivative, and for computing Y(tend) based on
the computed state derivatives.
Thus, from a mathematical perspective, the work of a simulation lies "merely" in
the task of periodically computing the state derivative. For ead? value of I for which
the state derivative wdty(I) is to be computed, the sequence of steps in Figure 3.1
is performed. The meanings of the steps in Figure 3.1 are explained below. The
computation of the state derivative will sometimes be referred to as a lime step;
successive time steps refers to successive computations of the state derivative for
different values of I. The inIe?va1 between two time steps is the difference in the
value of I between the time steps.
3.1 Updating the current state
NVhen the numerical ODE solver requests a value for ddty(i) at some time I0, it
begins by divulging the cuvvent state of the system; that is, its estimation of the
value of Y(t0).
For the particle/surface problem, the current state is simply the vector containing
x(t0) and v(i0). For a collection of n rigid bodies, the current state is somewhat
more complicated, but conceptually the same. The ODE for a collection of ?? rigid
23
input Y(t0)			output d
ThtY(to)
Exception Messages
I			Update			Penetration			Discontinuity			Constraint/friction
current			detected			force
state			Estimate: te			determination
I			detection			determination			computation
Figure 3.1: The computational steps required to compute wdty(t) During the
computation, various "exceptions" can be communicated to the ODE solver.
24
bodies has the form
c1(t)			vi(t)
qi(t)			-21wi(t)qi(t)
F1(t)/n?1
wi(t)			IF1Tl(t)
ft?(t)=ft			(3-1)
cn(t)			vn(t)
--H21w?(t)q?(t)
vn(t)			Fn(t)/?nn
wn(t)			I??1Tn(t)
where cj(t) is the position of the ith body's center of mass, q?(t) is the uiiit
cluaternion representation of the ith body's orientation[21,48], vj(t) and wj(t) are
the linear and angular velocities of the ith body, Fj(t) and Tj(t) are the total force
and torque acting on the ith body, and fl?j and Ij(t) are the ith body's mass
and inertia tensor. The term w?(t)q?(t) denotes quaternion multiplication between
o + coxi(t)i + Wyj(t)j + wzi(t)k and q?(t)
The values of all the position and velocity variables of the bodies are taken from
Y(t0) during the computation the denvative state at time to.
3.2 Collision detection
In the particle/surface problem, the motion of the particle is treated as uncon-
strained whenever the particle lies above the surface. Suppose that at some time to
the particle is above, but near the surface, and has a velocity dowuwards (towards
the surface). Based on the history of the particle's velocity and acceleration up
to time to, the ODE solver computes a position for the particle at some ti'?e
to + At. However, because the ODE solver has no knowledge of the non-penetration
25
constraint between the particle and the surface, it does not know that the particle
cannot continue to move downwards indefinitely. The ODE computes a state
Y(t0 + At) for the particle such that the particle lies below the surface; clearly,
this is incorrect. The ODE soh'er then requests the value of f(t0 + At, Y(t0 + At)),
unaware of the fact that Y(t0 + At) is an illegal (and hence incorrect) value for the
state at time to + At.
It is the collision detection step's responsibility to examine the current state
(or rather the proposed current state). If the state is illegal, that is, if the non-
penetration constraint has been violated, then at some time tc lying between 1o and
to + At, a collision must have occurred. The collision detector makes some estimate
te of tc,informs the ODE solver that the proposed state Y(t() + At) is illegal, and
gives the estimated value te. The value of f(t0 + At, Y(t0 + At)) is undefined
since the proposed state Y(t0 + At) is illegal, and computation of ?ddtY(to + At) is
aborted. The ODE solver proceeds forward from time to and computes a proposed
state for Y(t0 + te).
If te < tc, the proposed state Y(t0 + te) is legal and it is now known that
te < tc < to + At. Based on information fi?m the state Y(te), it may be possible
to produce a better estimate for tc, whid? the solver can aim for. Failing that, the
solver can at least bisect the time interval fron? te to to + At by estimating tc to
be ?21(te + to + At)
Conversely, if tc < te, then the proposed state Y(te) is illegal, and the solver
knows that to < tc < te. Again, either a better prediction for tc can be made
based on Y(te), or the tiine interval can be bisected. Eventually, the solver finds
te to within a suitable numerical precision; that is, a state is reached for which the
non-penetration constraint is satisfied within numerical tolerances and the velocities
of bodies are such that a collision is occurring.
One potential difficulty lies in the fact that if two successive states Y(t0) and
26
Y(t0 + At) are determined to be legal, it is still possible that a collision might
have occurred somewhere between these two states. This could happen if one body
passed completely through another between time to and time to + At. Consider for
example simulating the motion of a bullet shot through a thin pane of glass. In this
case, the odds are against the solver requesting TdtY(t) at a vahie of t for which the
bullet inter-penetrates the glass. However, for reasonably complicated simulations,
the nui?erical considerations in solving the differential equations of motion limit the
interval between time steps and the problem arises infrequently. For this reason,
we choose not to treat the case where a body passes completely through another
body in one time step. Possible solutions to this problem are using some type of 4D
space-time volui?e algorithms to detect collisions, or simply limiting the maximum
interval between time steps.
3.3 Contact point determination
Once a legal state Y(t0) has been achieved, the constraint forces at time to must be
computed. The first step in that computation is determining where the constraint
forces can arise. For the particle/surface problem, this determination is simple; if
the particle lies on the surface (within some numerical tolerance), a contactpoint is
said to exist between the particle and the surface. The location of this contact point
is simply the position of the particle. If a constraint impulse or a constraint force
is needed, it will act between the particle and the surface at the contact point. For
systems of rigid bodies, there may be arbitrarily many points of contact between
the bodies, and a constraint force may act at each point. Thus, it is necessary
to first determine all the contact points between the bodies before giving further
consideration to the constraint forces.
In this thesis, contact between bodies has been restricted to the case when it
27
can be modeled as occurring at some finite number of points between bodies. This
restriction does not apply to contact between polyhedral bodies. For polyhedral
bodies, contact areas are always polygonal (considering single points and line
segments as degenerate polygons). By imposing a non-penetration constraint at
each vertex of the polygonal contact area, inter-penetration is prevented over the
entire polygonal contact area. Additionally, any distribution of normal forces over
the contact area is equivalent to some assignment of normal forces acting only at
the vertices of the contact region. Thus, the definition of a contact point for the
case of polyhedral contact is limited to the vertices of the contact region.
Although the collision detection and contact point determination steps are listed
separately in Figure 3.1, the computation of each overlaps. As bodies are checked to
see if non-penetration constraints have been violated, contact points are determined
at the same time. Both steps can share a significant amount of computation, and
performing them separately would be inefficient. NVhen a proposed state is found
to be illegal, the work done so far in determining contact points is wasted; but it
would be much more wasteflil to implement the two steps completely separately.
3.4 Impulse computation
Once the contact points have been determined, the relative velocity at each contact
point can be computed. If the velocity at a contact point is such that the contact
is breaking, no non-penetration constraint need be formulated at the contact point;
equivalently, the non-penetration constraint at the contact point is inactive.
If the velocity at a contact point is such that a collision is taking place, an
impulsive force will arise. If it is assumed that no two collisions occur at exactly
the same time, the calculation of the impulse, for frictionless surfaces, is exactly
as described in Section 2.5. Niethods for dealing with simultaneous impulses, and
28
impulses involving friction are summarized in Chapter 4.
Since impulsive forces cause discontinuities in the velocity variables vj(t) and
wj(t) of the rigid bodies, special care must be taken when impulses are applied.
Numerical ODE solvers assume that the flinction Y(t) is smooth. Since wdty(t) =
f(t, Y(t)) is only piecewise continuous, it must be treated as such. Accordingly,
when an impulse is applied at time to, new values for the affected vi(t) and wi(t)
n?ist be computed, and the ODE solver must be informed that a discontinuity has
been encountered. The ODE solver is then restarted with a new state Y+(t0), where
Y+(t0) contains the new values of the affected vj(t) and wj(t). Note that if the ODE
solver begins by requesting a vahie for wdtY+(to), the collision detection and contact
point determination information computed for Y(t0) can be used in computing
since Y(t0) and ??+(to) represent the same geometric positionings of
bodies.
3.5 Constraint/friction force determination
This step is the easiest to describe, but the most difficult to perform. At this
point in the computation of wdty(to), all contact points with active non-penetration
constraints have been identified. Since any constraints requiring impulsive constraint
forces have already been taken care of, all that remains is the formulation and
sohition of a system of equations for the unknown constraint forces. NVhile constrahit
forces always act normally between contacting bodies, contacts involving friction
give rise to forces that act tangentially between bodies, to oppose slipping. For
contacts with friction, a system of equations for both the unknown constraint forces
and tangential friction forces needs to be formulated and solved.
29
Once the constraint forces and the friction forces have been computed, the right-
hand-side of Equation (3--H1) is trivially computed, and the computation of wdty(to)
is complete. The ODE solver computes a new (proposed) state Y(10 t At), and
computation of wdty(to + At) begins.
Chapter 4
Summary of Previous Work
In this chapter, previons work relating to the various sin?ilation steps outlined in
Chapter 3 is summarized.
4.1 Collision/contact determination
The two problems of collision detection and contact point determination are very
similar, and can really be considered a single problem, which we shall call the
collision/contact determination problem. There is an m?mense amount of litera-
ture concerning this problem, but the majority of collision/contact determination
algorithms fall into one of two groups. First, the collision/contact determination
problem from time t0 to t1 can be viewed as a single, continuous function of time.
Given this viewpoint, the basic problem to be solved is "at what time, and where, do
bodies first come into contact?" Second, the problem can be considered discretely,
at a sequence of time vakies to < to t At1 < to + At2 < ... ? to + At? ? t1. In
this viewpoint, the basic problem is "given the positions of bodies at time to + Ati,
where do bodies inter-penetrate and contact each other?"
The first approach, which could be called the "continuum view", presupposes
31
32
a specified motion of bodies over some time interval. Examples include Canny[6],
who describes an algorithm for determining the first collision between rigid poly-
hedral objects with constant angular velocity. Using the constant angular velocity
assumption, Canny reduces the problem of determining the first instant of collision
to the problem of finding roots of (relatively low order) polynomials. Gilbert and
Hong[18] relax the problem of collision detection to include arbitrary trajectories of
rigid convex polyhedra. Since no closed-form solution exists for time at which the
the first intersection occurs, an iterative numerical method is used to determine this
time. Their algorithm is based in part on previous work by Gilbert, Johnson, and
I&eerthi[19]. Von IIerzen, Barr, and Zatz[53] describe an algoritlim that determines
the first collision between parametrically defined time-dependent curved s'irfaces. In
this work, the motion of the bodies is not necessarily rigid (that is, the actual shape
of the surfaces may vary over time). Unfort'inately, a?gorithms of the above nature
cannot be used in the simulation problem as o+`iescribed in Chapter 3. All of the
above algorithms presuppose a specified motion for the bodies over time; however
such a path is exactly what the simulator is trying to compute, which means we
cannot make use of these sorts of collision detection algorithms.
The second approach, which could be described as the "discrete ",involves
the geometric analysis of bodies at a fixed instant of time. The goal under
this approach is to produce algoritlims that solve a single, static instance of a
problem with optimal asymptotic time complexity. NNTork of this sort includes an
approximately O(n log n) algorithm by Cameron and Culley[5] for computing the
distance between two convex polyhedra with n vertices. Gilbert, Johnson, and
I?eerthi[19] describe an algorithm for computing the minimum distance between
convex polyhedra with n vertices; the algorithin appears to have a running time
somewhat larger than 0(n) but smaller than O(n log n). Related work by Gilbert
and Foo[17] extends the algorithin to handle smooth convex shapes and concludes
33
that smooth representations become more efficient than polyhedral representations
in the neighborhood of n = 100. For the problem of determining disjointness of two
convex polyhedra with n vertices (without regard to the distance between them), an
0(n) algorithm is readily available, based on Megiddo's work on linear programming
problems with constant dimension[37,43].
Unfortunately, the above algorithms, although they ad?ieve small asymptotic
time complexity, do so at the price of large run time constants. For example,
the linear time algoritlim for deciding disjointness between convex polyhedra with
n vertices is not practical unless n is quite large. More importantly, the `ise
of discrete view algorithms essentially ignores any similarity that may exist with
the current collision/contact determination problem and previous collision/contact
determination problems. Proper use of previous information can result in very
sii?ple, yet efficient collision/contact determinatiou algoritlims; such algorithins are
presented in Chapter 5. The collision detection algoritlim in Gilbert, Johnson, aii?l
I&7eerthi[19] is structured so that it can use information previously computed to
obtain a faster running time, but this is not a main consideration of the algoritliin.
Likewise, work by Cundall[lOj on dynamic simulation of polyhedral blocks uses a
collision/contact determination algorithm that specifically exploits information 01)-
tained from the previous determination However, Cundall's approach is somewhat
more complicated than needed for dynamic sin??lation and is restricted to polyhedra.
The philosophy of collision/contact determination presented in Chapter 5 was largely
inspired by Cundall's approach.
4.2 Impulse computation
The computation of impulses between two bodies that collide at a single contact
point without friction is straightforward, as it involves the solution of a single linear
34
equation with one unknown. If equality constraints involving one of the colliding
bodies exist (as happens when an articulated rigid body undergoes a collision),
additional impulsive constraint forces are needed to maintain the equality constraint.
In this case, the impulses must be computed simultaneously, so a larger set of linear
equations must be solved; this is simply a special case of solving equality constrained
dynamics problems, and is not particularly difficult.
If the contact involves friction, analysis of the problem is more complicated, but
the actual computation itself is simple, for two-dimensional systems. ?Vang and
Mason[34,54] discuss the case of collisions with friction at a single contact point in
two dimensions, and give a complete treatment of the problem. Their analysis of the
problem yields simple forn?ilas for the impulses arising from a collision. However,
?Vang and Mason's results are difficult to extend to three dimensions, which can be
seen by considering I?eller[24], who presents a differential description of the three-
dimensional motion of bodies undergoing a collision with friction at a single contact
point. The differential equations appear easy to solve numerically, but very difficult
to treat analytically. I&7eller's analysis involves regarding the collision as occurring
over an arbitrarily small time interval, and examining the limiting behavior of the
system as the time interval is brought to zero. ?Ve shall employ the same sort of
analysis in considering impulses with friction in chapters 12 and 13. Aii animation
system described by Moore and ?VillIelms[39] ?ises penalty methods to compute
constraint forces, but computes constraint impulses directly for collisions involving
a single contact point. Moore and ?Vilhelms use a simplified description of the
Coulomb friction law to compute frictional impulses at a single contact point in
three dimensions. We will use a similar simplification, applied to a number of
contact points simultaneously (see below).
Collisions can also be regarded as occurring simultaneously at a number of
contact points. For example, if a cube is dropped onto a level plane so that all
35
four vertices of the bottom face of the cube strike the plane together, the collisions
between the vertices and the plane can be modeled as occurring at discrete instants
of time, by the above methods, or as occurring simultaneously. There are large
computational savings to be gained by modeling the collisions as occurring simulta-
neously. Cremer[8] discusses the advantages and disadvantages of the simultaneous
model of collisions, and Featherstone[13] shows how to formulate the necessary
equations for the frictionless case. The impact model for simultaneous collisions
used in this thesis agrees with Featherstone's formulation, for the frictionless case.
L5tstedt[30J computes simultaneous frictional impulses in three dimensions by using
a modification of the Coulomb friction law that causes impacts to dissipate as much
energy as possible. In general, it is unclear how to deal correctly with sin?iltaneous
impacts with friction, in either two or three dimensions. The simulator described in
this thesis computes simultaneous fi?ictional impulses for three-dimensional bodies,
based on a modification of the model proposed by L5tstedt. As stated, for frictionless
systems, the model described in this thesis is the same as Featherstone's. However
for systems with friction, the model does not always properly conserve energy
for impacts[3,54]. The general problem of correctly computing frictional impulses
in three dimensions, and simultaneous frictional impulses in either two or three
dimensions is not addressed by this thesis.
4.3 Penalty methods for constraint forces
A vast number of simulations[10,22,2(3,3G,38,39,50,51,55] have employed the penalty
method as described in Section 2.4 to enforce non-penetration constraints; appli-
cations include the simulation of deformable bodies, cloth, and articulated rigid
bodies. The penalty method is a very attractive model in some respects; it is
extremely simple to implement, and very versatile. As a numerical method, it is
36
nowhere near as complex as the methods we will consider throughout most of this
thesis. In particular, exploiting parallelism and/or hardware implementations with
the penalty method is much easier to imagine than with the methods we will be
discussing throughout most of this thesis. Practically speaking though, the use
of the penalty method in essentially dynamic conditions leads to stiff differential
equations and is not a suitable method for the type of simulation we are interested
in. The penalty method is reasonable only for situations that are essentially static;
that is, where there is little relative motion between bodies and little to no change
in constraint forces between bodies. ?Te will consider the additional question of the
theoretical correctness of the penalty method as a simulation algorithm in Chapter 7.
4.4 Exact methods for constraint forces
In comparison with the penalty method, relatively few dynamics systems use exact
methods to calculate constraint forces. Descriptions of the system of equations
in terms of a quadratic program for calculating n?iltiple constraint forces without
friction (for polyhedra) appear in Nilmister and Reeve[25] and Ingleton[23]. Cottle[7]
clarifies and simplifies some of Ingleton's results. Nilmister and Reeve, Ingleton, and
Cottles' work will be our basic starting point when we consider systems with multiple
contact points in Chapter 8. ?Ve will make use of the same quadratic programming
formulation and extend the contact geometry allowed to include smooth curved
surfaces.
All of the above work took place before the discovery of the importance of P
and NP time-complexities. As a result, the computational complexity of computing
n?iltiple contact forces was not of concern in any of the above work. In work
done subsequent to this discovery, Erdmann[12] discusses the problem of computhig
constraint forces without friction, and notes that the problem is simply solved by
37
an exhaustive search method, requiring exponential time. The earliest description
of a simulator using exact methods to calculate constraint forces (by quadratic
programming) appears to be by L6tstedt[29]. Lo?tstedt notes that the quadratic
program can be solved efficiently because it is convex. More recently, Feather-
stone[13j has described a heuristic algorithm for solving the quadratic programs
generated for configurations without friction that is intended to improve upon the
obvious exhaustive search method. However, Featherstone's heuristic requires that
configurations of bodies be statically determinate; that is, the normal forces that
prevent inter-penetration for the configuration must be unique. This is a very
restrictive condition though, which we will not want to impose on our simulations.
In this thesis, the computational complexity of computing n?iltiple contact forces
will be of great interest to us.
4.5 ?riction
Adding a model of friction to simulations using the penalty method is straight-
forward. Computing both friction and constraint forces using exact methods is
considerably more difficult. For two-dimensional problems, an obvious exhaustive
search method exists, but is not suitable to simulations with many contact points
because it requires exponential time in the nun??er of contact points. Mason and
Wang[35] discuss the computation of fiiction and constraint forces in systems with
only a single contact point (where the time complexity of tlie computation is not
an issue). L6tstedt[28], Erdmann[12], and Mason and NVang all discuss interesting
properties of the Coulomb friction model in relation to rigid body simulation. Their
discussions and examples of the unexpected consequences of introducing Coulomb
friction provide the starting point for much of the discussion on friction in this
thesis. In particular, Erdmann describes the problem of computing contact forces
38
in terms of configuration space (Chapter 6). We will borrow from Erdmann by
using configuration space to develop the dynamics of systems with a single point
of contact, although we will dispense with configuration space when we move on
to systems with multiple contact points and friction, as the configuration space
approach offers no additional insights into the problem. Additionally, L6tstedt[30]
proposes a modification of the Coulomb friction law for rigid body simulation in
conjunction with a practical computational solution method for the modified friction
model. We will examine L5tstedt's modification and subscqucnt solution method
in Chapter 14. Again, the computational complexity in computing multiple contact
forces with friction will be of considerable interest.
Chapter 5
Collision Contact Determination
In summarizing previous work on collision/contact determination in Section 4.1, we
noted that of the two groupings of algorithms (the discrete view and the continuum
view), the continuum view, which is comn?only used in robotics, is not applicable to
the dynamics problem we are interested in. The other approach, the discrete view
could be used to solve the collision/contact determination problem, but at the price
of ignoring any knowledge due to the inherent spatial contin'iity of the problem as a
function of time. Instead, a third approach, which we will call the "coherence ??
will be adopted and used to solve the collision/contact determination problem. The
use of coherence need not be confined to collision/contact determination; in fact,
critical use will be made of coherence in computing constrahit and fiiction forces. In
general, methods that exploit coherence should be considered throughout the entire
simulation process.
As stated in Chapter 3, the sin?ilation process consists of the repeated compu-
tation of wdty(t) at various times t. The numerical ODE solver is responsible for
choosing the values of t at which the state derivative is to be computed. For any
reasonably complicated simulation, the values of t chosen are such that the state
39
40
does not change greatly between successive values of t. As a result, there is almost
always great geometric coherence between successive time steps. It is difficult to
formally quantize the degree of coherence, and we will not attempt to do so. In
practice, the degree of coherence present during simulations has been such that
the algorithms presented in this chapter are much faster than their discrete view
counterparts. Additionally, their simplicity has greatly decreased the programming
difficulties that arise with any complex simulation system.
The coherence view attempts to capitalize on geometric coherence to achieve
both simpler and more efficient collision/contact determination algorithms.1 The
series of collision/contact determination problems that need to be solved during the
simulation are considered, according to the coherence view, as separate problems;
however, the problems are considered to be related, and wherever possible, informa-
tion from the solution of the previous problem is used to aid iii the sollition of the
current problem. Using this view, algorithms that are both faster and simpler than
their discrete view algorithms counterparts are possible.
In this chapter, we will discuss three separate problems of collision/contact
determination. The first problem we consider is pairwise collision/contact deter-
mination between convex polyhedra. Next, pairwise determination between convex
curved surfaces or a convex curved surface and a couvex polyhedron is considered.
Last, we will consider a simple coherence-based algoritlim for detecting all pairwise-
intersections among a group of rectangular solids aligned with the coordinate axes.
This algorithm is used at the beginning of the collision/contact determination to
quickly cull out pairs of objects whose rectangular bounding boxes do not overlap;
1 For those interested in building a practical, robust sin?jlation system with a rich variety of
object geometries, it is sometimes hard to say whidl is more valuable increased simplicity or
decreased running time. Given the tendency complicated geometrical algoritlims have to get stuck
and abort on the special cases that always arise in practice, it sometimes seems that slow but
simple algorithms end up saving time over a complicated algon.thm that requires the sin?ilation
to be run and re-run and re-re-run and ...
41
such a step is necessary to avoid performing 0(n2) pairwise collision/contact
determinations among a group of n bodies. Using the algorithms developed to solve
these three problems, we will be able to perform collision/contact determination
on any body that can be described as the union of convex polyhedra and strictly
convex closed surfaces. These polyhedra and curved surfaces will be referred to
as primitives. For curved surfaces primitives, only closed surfaces with strictly
positive curvature at every point will be allowed for collision/contact determination.
Note that although polyhedral primitives must also be convex, this is not really
a restriction, since any polyhedron can be decomposed into a union of convex
polyhedra. The general case of collision/contact determination between non-convex
curved surfaces is beyond the scope of this thesis; however, this limitation will not
be imposed when considering the dynamics of bodies in subsequent chapters.
The discussion of the collision/contact determination step in Section 3.2 iuade
mention of the fact that the collision/contact determination step was responsible
for estimating the time of a collision. This will be covered in Chapter (3.
5.1 Convex polyhedra
Our primary med?anism for exploiting coherence will be through the ?ise of witnesses
to the decision problem of inter-penetration between primitives. A witness is soniC
piece of information that can be used to quickly answer a decision problem. NN?e will
utilize coherence by caching witnesses from one tinie step to tlie next; hopeflilly a
witness from the previous time step will be a witness d?iriiig the current time step.
Since primitives are convex, a pair of primitives do not inter-penetrate if and only
if a separating plane between them exists. A separating plane between two objects
is a plane such that eadi object lies in a different half-space of the plane. A given
plane can be verified to be a separating plane between two polyhedra in time 0(n)
42
where n is the total number of vertices in the two polyhedra. Thus, a separating
plane is a witness to the fact that two convex polyhedra do not inter-penetrate.
We would like to structure our algorithms so that their running time decreases as
coherence increases. The basic philosophy behind utilizing coherence through cached
witnesses is to choose a very simple, fast strategy that proposes a new witness based
on the cached witness, and then quickly attempts to verify the new witness. If the
proposed witness is invalid, a second, more time-consuming step is executed to find
a valid witness. The idea is that the first step is fast enough, and succeeds often
enough to more than make up the time lost on the occasions when it does fail.2
In particular, the strategy should be extremely simple, so that it has both a small
asymptotic time complexity and run time constant. The cost of initially finding a
witness (for the very first time step of the sin?ilation, or the first time two bodies
become close enough to require more than a bounding box test) is unavoidable;
either exhaustive search algorithms or the discrete view algorithms described in
Section 4.1 can be used to initially compute these witnesses.
As mentioned in Section 4.1, algorithms by Cundall[10] and Gilbert, Johnson
and I?eerthi[19J come close to this philosophy. In particular, the algorithm Gilbert,
Johnson and I?eerthi describe could be initialized with the solution of the previous
problem to decrease the running time; however, the initialized algoritlim is still n?ich
more complicated than is needed for situations with a high degree of coherence. Cun-
dall's approach to collision/contact detection is much simpler, and in fact inspired
the coherence philosophy just presented. Cundall called a separating plane between
two convex polyhedra a common plane. For each pair of convex polyhedra requiring
consideration, he used numerical techniques to initially determine a separating plane
that was approximately equidistant from each polyhedron of the pair; if such a plane
2 In general, such a philosophy can be implemented hierard?ically, as is done for memory caching
on computers, but we will not bother to do so.
43
could not be found, then clearly the two polyhedra had inter-penetrated. He then
noted several important facts about the separating plane. If the polyhedra were
disjoint, this was easily shown because each polyhedra could trivially be shown to
lie some distance from the separating plane. Furthermore, if the polyhedra were
contacting, then the contact points were easy to find because they would have to lie
on a common plane of contact (hence the term common plane); this greatly reduced
the complexity of determining the contact points. Most important was Cundall's
realization that it was not necessary to compute separating planes from scratch at
each time step; a simple numerical iteration could compute a new separating plane
by using the old separating plane as a starting point, because of the geometric
coherence between successive time steps.
Cundall's approach is good, but it can be made even better. If a pair of convex
polyhedra are disjoint or contacting (but not inter-penetrating), then a separating
plane exists with the following property: either the plane enil?eds a face of one of
the polyhedra, or the plane embeds an edge of one of the polyhedra and is parallel
to an edge of the other polyhedra. NVe will call the face or edges in question the
defining face or edges. Instead of caching a specific plane in R3 as the witness and
numerically updating it, the defining face or edges are cached as the witness. At the
next time step, a separating plane is formed from the defining face or edges, and then
verified. The verification is simple; the vertices of the t?? polyhedra are checked to
see that they lie on opposite sides of the separating plane. Since this verification
takes 0(n) time (where n is the number of vertices), a linear time algorithm results
whenever the cached face or two edges form a separating plane. Clearly, in addition
to having a linear time cost, this algoritlim will be considerably faster than other
0(n) algorithms, such as the linear programming algorithm discussed in Section 4.1.
However, an additional improvement can be made so that this verification can
usually be performed in sublinear time. It is not always necessary to test all the
44
vertices of a polyhedron to verify that the polyhedron lies completely on one side of
the separating plane. Instead, for each polyhedron, the vertex lying closest to the
separating plane is cached. If at the next time step that cached vertex is closer to the
separating plane than all of its neighbors, then by convexity all of the vertices of the
polyhedron must lie on the same side of the separating plane. This gives the result
that disjointness or contact between convex polyhedra can usually be determined in
sublinear time when coherence is high.
On those occasions when the cached face or two edges fails to form a valid
separating plane, faces or edges adjacent to the previously cached face or edges can
be used to form a possible separating plane; however, this happens infrequently
enough that it may be simpler to start from scratch and compute a new separating
plane without using any prior knowledge. It should be obvious that the linear time
cost of verifying a separating plane is of a n?ich different order than the linear time
cost of using a technique like linear programming[43J to find a separating plane;
although both algorithms have 0(n) time complexity, their run time constants are
considerably different.
Once the separating place has been fo'ind, the contact region between the two
polyhedra is determined, assuming the polyhedra are not disjoint. Contact points
between the two polyhedra can only occ?ir on the separating plane. Given the
separating plane, the contact points can be quickly and efficiently determined by
comparing only those faces, edges, and vertices of the polyhedra that are coincident
with the separating plane.
However, if no separating plane can be found, then the two polyhedra must
be inter-penetrating. ?7hen two polyhedra inter-penetrate, it is almost always the
case that either a vertex of one polyhedron is inside the other, or an edge of one
45
polyhedron has intersected a face of the other.3 In this case, the inter-penetrating
vertex, or intersecting edge and face are cached as a witness to the inter-penetration.
Since this indicates a collision at some earlier time, the simulator will back up
and attempt to compute wdty(t) at some earlier time. Until the collision time
is determined, the first action taken by the collision/contact determination step
will be to check the cached vertex or edge and face to see if they indicate inter-
penetration. Thus, until the collision time is found, states in which the inter-
penetration still exists are identified as such with a minimum of computational
overhead. In Chapter (3, we will show how the collision time may be estimated once
an inter-penetration is detected.
5.2 Convex curved surfaces
The surfaces considered in this thesis are assumed to be twice-differentiable and
without boundary; in this chapter, surfaces are also assumed to be strictly convex.
Before we can consider collision/contact determination between curved surfaces, we
must decide how the curved surfaces are to be represented.
There are basically two choices surfaces can either be described by implicit
functions or parametric fanctions. In the implicit case, a surface is described as a
time-varying function F(J), t) where p is a point in world space. At time 1, point ])
is on the surface represented by F if
F(p,t) = 0.			(5--H1)
In the parametric case, a surface is described by a time-varying fanction S(u, v,
where u and v are scalars and S(u, v, t) is a point on the surface. If F and S
3An exception is the following. Stack two cubes of equal size atop one another so that their
contacting faces exactly coincide. Lower the top one. This produces an inter-penetration such that
no vertex is inside either cube, and no edge penetrates through any face.
46
describe the same surface, then
F(S(u,v,t),t)
is identically zero.
NVe will not attempt any formal manipulations of F and S; that is, it is assumed
F and S are given to us as "black-box" functions. The only allowable operations
on F or 5 are evaluating F or 5, or their first or second derivatives at specific
values of p and t (for F) and n, v and I (for 5). For simplicity, derivations will be
developed in the text only for surfaces described implicitly. Niatching derivations
for parametrically defined surfaces are found in Appendix A.
?Ve will consider two convex c'irved suffaces A and 1? that are described
implicitly by time-varying functions F(p, I) and G(J), I). Any point i)o lying on
A at time t satisfies F(p0, I) = 0. Furthermore, a point Pi that is inside of A
satisfies F(p1, I) < 0, while a point P2 that is outside of A satisfies F(p2, I) > 0.
The other surface B is similarly defined in terms of &, and we will sometimes refer
to F and G as both functions and surfaces.
The surface normal of F at point p and time I is written as VF(p, I), and is a
col'imn vector. Formally,
__ = F'(p,t) = vF(p,I)T
Since F is negative on the interior of A and positive on the exterior of A, the
surface normal VF(p, I) is always outwards pointing.
5.2.1 Extremal points ofcurved surfaces
First, let us consider the problem of verifying that two convex curved surfaces A
and B do not inter-penetrate. Since A and B are strictly convex smooth surfaces, if
F1
My
47
F2
A
Pam
Pb
G
Figure 5.1: The minimum distance points Pa and Pb between two convex surfaces,
they do not inter-penetrate then there is a unique pair of points Pa and Pb lying on
A and I? that minimize the distance between A and B (Figure 5.1). If these points
can be found, then they clearly demonstrate that A and I? are not inter-penetrating.
Furthermore, if the bodies are in contact, then Pa = Pb is the contact point between
A and B. Thus, in both cases, the pair of points Pa and Pb minimizing the distance
between A and B seems a good choice for a witness.
?Vhat about the case when A and I? inter-penetrate? In this case, the minimum
distance between A and B is zero. In detecting contact between surfaces, if the
depth 6 of inter-penetration between A and B falls below some numerical threshold
?,then A and B must be considered to be contacting, rather than inter-penetrating.
For instance, in Figure 5.2, suppose the iliter-penetration depth between F1 and G
is 6i similarly, let 62 be the inter-penetration depth between F2 and G. From tlie
figure, 6? <62. If
6i <? < 62
then F1 and & would be considered to be contacting, while F2 and G would
be considered to be inter-penetrating. Thus, when surfaces are not disjoint, the
minimum distance, which is always zero, is not a useful piece of information.
Although the points of intersection are not directly useful, one thought might be to
examine the angle between surfaces at intersection points; however, as Figure 5.2
48
F1
maximum
interpenetration
Figure 5.2: Inter-penetration between surfaces. Surface F1 intersects G at a sharper
angle than surface F2 does, yet the i?aximnn? inter-penetration depth 62 between
F2 and G is larger than the n?axin?im inter-penetration depth 6j between F1 and
shows, the angle of intersection is independent of the depth of inter-penetration. For
inter-penetration, the points of interest are the points that indicate the ma?irnuni
distance between A and B in their intersection. Thus, a suitable witness for inter-
penetration are the points of maximum distance between A and B.
?Ve can formalize and unify our definition of a witness for both tl-ie inter-
penetrating situation and the non-inter-penetrating situation by defining the general
concepts of the ertrerne distance between A and B, and the extremal points between
A and B. The extreme distance between A arid B is defined as follows. If A and
B are disjoint, then the extreme distance between A and B is just the normal
minimum distance between A and B. If A and B are in contact, then the extreme
distance is zero. If A and B have inter-penetrated, then the extreme distance is the
maMmum distance between A and fi in their intersection. The extremal points of
A and B are the two points Pa and p?,lying on A and B respectively, that realize
the extremal distance (Figure 5.3).
The extremal points can be computed by considering the following four necessary
F1
My
G
49
F2
mm
Pb
F3
Pa
Figure 5.3: The definition of the extremal points. The extremal points Pa and Pb
between F1 and G are the points of minimum distance. The same is true for F2
and G. The extremal points Pa and Pb between F3 and G are the points that
maximize the distance between F3 and G in their intersection.
conditions for a pair of points Pa and Pb to be extremal at time t:
1 : ?F(Pa,t) t A2V&(ft,t) = 0
2 : F(Pa,t) = 0
3: G(Pb,t) = 0
4 : (Pb ?pa) + AlVG(Pb,t) = 0
The variables A1 and A2 are unconstrained scalar values (with no relation to the A
of Section 2.5) and 0 is a column vector of zeroes. Condition i guarantees that
the surface normals of F and G at Pa and Pb are colinear. Conditions 2 and 3
guarantee that Pa and Pb are points on A and B, and condition 4 guarantees
that Pb'5 displacement from Pa is colinear to the surface normals. Figure 5.4 shows
the geometric intuition behind these four conditions; they are derived based on the
Lagrange multiplier formulation for constrained minimization (hence the dioice of
A for the multiphers of VG).
These four conditions form a non-linear set of eight eq?iations in eight unknowns.
However, the conditions are necessary but not sufficient for Pa and Pb to be extremal
points. Because Equation (5--H2) is non-linear, it will have n?iltiple solutions. For
example, the points that globally maximize the distance between A and D always
50
VG(pb,t)
F2
y
Pb
VG(pb,t)
F3
Pb
Pa
F1
My
VG(pb,t)
G
Figure 5.4: Geometric conditions for the extremal points. In order for Pa and Pb to
be extremal, the surface normals VP and VC must be parallel, and point directly
from Pa to Pb
satisfy Equation (5--H2), but are not a sohition we are interested in. Non-linear
equations are usually solved numerically using iterative methods. The numerical
method proceeds from some initial estimate of the solution to an exact solution
(within numerical tolerances). The initial estimate determines which solution is
computed[1 1].
If we are initially computing the extremal points, we need an initial estimate
sufficiently close to the extremal points to avoid computing a spurious solution. For
implicitly defined convex surfaces, the following simple device can be used. A line
segment passing though the centroids of A and B is constructed. The intersection of
this line segment and A is roughly computed; that is, a point p on the line segment
satisfying F(p, t) = 0 is computed. This can be done very easily using ,
and since only a rough approxhuation is needed, a small mimber of bisection steps
will be sufficient (ten or less). Similarly, an intersection between the line segment
and B is computed. Using the two intersections as initial estimates for Pa and
Pb, Equation (5--H2) is iteratively solved to find Pa and Pb exactly. For parametric
surfaces, a coarse mesh of surface points can be generated and stored for each
surface. The parametric coordinates of the minimum distance pair of mesh pohits
51
can be used as an initial estimate. The initial solution of Equation (5--H2) requires
a fair amount of work, in terms of generating an initial estimate of Pa and Pb, and
then solving Equation (5-2).
However, once the extremal points Pa and Pb are determined, they are cached
for the next step. In subsequent time steps, the cached values of Pa and Pb are
used as initial estimates in the numerical solution of Equation (5--H2). Because of
coherence, the cached values Pa and Pb will greatly reduce the number of iterative
steps needed to numerically solve Equation (5--H2). Furthermore, the accuracy of
the initial estimate can be substantially improved if the derivatives of the extremal
points are cached at each time step. Suppose that at time t0,the extremal points
Pa and Pb, and their derivatives Pa and Pb are cached. If the next time step takes
place at time to + At, Pa and Pb at that time can be estimated as
Pa+AtPa			and			Pb+AtPb			(5-3)
respectively. Improving the accuracy of the initial estimate of Pa and Pb increases
the overall speed of the algorithm. Sections 6.3 and 6.4 will discuss the computation
of Pa and Pb
5.2.2 Polyhedron/curved surface determination
Essentially the same method for collision/contact determination between a pair of
curved surfaces can be used for collision/contact determination between a curved
surface and a polyhedron. The concept of extremal points is again used, however,
different formulations are needed depending on whether the extremal point on the
polygon lies on a face, an edge, or a vertex.
Let the curved surface be denoted A, and be described implicitly by a function
F as before. Let B denote the polyhedron. Let the plane equation of a face
of B be implicitly described by a function G(P, t). Necessary conditions for the
52
extremal points to occur between this face and A are given by Equation (5--H2),
with the restriction that Pb actually lie inside the face. A simple initial step to
compute Pa and Pb is to test each face of B against A, stopping if valid extremal
points are found. Clearly, not all faces of B need be considered; faces of B whose
outwards surface normal points away from A can obviously be ignored. However, if
no extremal points exist between a face of B and A, then the edges of B must be
tested.
An edge of B with endpoints P0 and P1 can be parameterized as the set of
points P0 + c r, where r = P1 --H P0 and 0 < c < 1. Necessary conditions for Pa and
Pb = Po + cr to be extremal are
vF(Pa,t)Tr = 0
F(pa,t) = 0			(5--H4)
(Pa? (P0 + cr)) + AiVF(Pa,t) = 0
along with the condition that Pb be contained in the edge; that is, 0 < c < 1. If
no extremal points exist between an edge of B and A, then an extremal point will
exist between a vertex P0 of B and A. The conditions for a vertex P0 of B and
Pa to be extremal are
F(Pa,t) = 0			(5--H5)
(Pa --H P0) + AlVF(Pa,t) = 0.
Once an extremal point is found between a face, edge, or vertex, the extremal
points are cached, along with the face, edge or vertex on which the extremal point
occurred. Subsequent solutions of the extremal points then take advantage of
knowing which face, edge, or vertex to test to find the extremal points, along with a
good starting estimate of the extremal points. In all cases, the velocity j)a is cached.
If the extremal point lies on a face, Pb is also cached, and if the extremal point lies
on an edge, c is also cached.
53
5.3 Bounding boxes
To reduce the number of pairwise collision/contact determinations necessary, a
bounding box hierardw can be imposed on the objects in the simulation environ-
ment. If two bounding boxes are found not to overlap, no further comparisons
involving the contents of the boxes are needed. Given a collection of n rectangular
bounding boxes, aligned with the coordinate axes, we would like to efficiently
determine all pairs of boxes that overlap. A naive pairwise comparison of all pairs
requires 0(n2) work, and can clearly be improved. Discrete view algoritlims for
this problem in R3 exist with a time complexity of O(nlogn + k) where k is the
number of pairwise overlaps; a general result is that the problem can be solved in
time o(n1ogd?2n + k) for d-dimensional bounding boxes[43]. Again, a coherence
view algorithm exists with substantially better performance.
5.3.1 The one-dimensional case
Consider the problem of detecting overlap between one-dimensional bounding boxes,
aligned with the coordinate system. Sud? a bounding box can be described simply
as an interval [b, e]. Let us consider a list of n such intervals, with the ith interval
being [bi, e?] c R1. The problem is then defined to be the determination of all pairs
i and j such that the intervals [bi, e?] and [bj, ej] intersect.
The problem can be solved initially by a sort and s'veep algorithm. A sorted list
of all the bj and e? values is created, from lowest to highest. The list is then swept,
and a list of active intervals, initially empty, is maintained. ?Vhenever some value bj
is encountered, all intervals on the active list are output as overlapping with interval
i, and interval i is then added to the list (Figure 5.5a). NVhenever some value e? is
encountered , interval i is removed from the active list (Figure 5.5b). The cost of
this process is O(nlogn) to create the sorted list, 0(n) to sweep through the list,
54
Ia
(a)
`4
15
Ii
_$
$
A			IIII
b3 b6 b1			e6 e1 b5 b2 e3
b4 e5			e4			e2
Figure 5.5: The sweep/sort algorithm. (a) When b1 is encountered, the active list
contains intervals 3 and 6; interval 1 is reported to overlap with these two intervals.
Interval 1 is added to the active list and the algorithm continues. (b) When e3 is
encountered, the active list contains intervals 2, 3 and 5. Interval 3 is removed from
the active list.
and 0(k) to output each overlap. This gives a total cost of O(n log n t k), and is
an optimal algorithm for initially solving the problem.
Subsequent comparisons can be improved ?s follows. First, there is no need
to use an O(n log n) algorithm to form the sorted list of 1?? and e? vallies. It is
considerably more efficient to start with the order found for bj and e? values from
the previous time step; if coherence is high, this ordering will be nearly correct. A
sorting method called an insertion sort[46] is used to permute the "nearly sorted"
list into a sorted list. The insertion sort algorithm works by moving items towauls
the beginning of the list, until a smaller item is encountered. Thus, the second item
is interchanged with the first if necessary, then the third item is moved towards the
beginning of the list until its proper place is found, and so on; each movement of an
item indicates a change in the ordering of two values. After the last item on the list
55
Figure 5.6: A coherence-based method of detecting overlaps. The order produced in
Figure 5.5 is nearly correct for this arrangement of intervals. Only b4 and e5 need
to be exd?anged. When the exchange occurs, the change in overlap status between
interval 4 and 5 is detected.
has been processed, the list is sorted. Such a sort takes time O(?1 + s) where s is the
mimber of "swaps" necessary to reorder the list. For example, the only difference
between Figures 5.6 and 5.5 is that interval 4 has moved to the right. Starting
from the ordered list of bj and Cj values of Figure 5.5, only a single exchange is
necessary to sort the list for Figure 5.6. The insertion sort is not recommended as
a sorting procedure in general, since it may require 0(n2) exchanges; however, it is
a good algoritlim for sorting a nearly sorted list, which is what occurs in our highly
coherent environment. To complete the algorithm, note that if two intervals i and
j oveilap at the previous time step, but not at the current time step, one or more
exchanges involving either a bj or e? value and a bj or Cj value must occur. The
converse is true as well when intervals i and j change from not overlapping at the
previous tii?e step to overlapping at the current time step.
Thus, if we maintain a table of overlapping intervals at each time step, the table
can be updated at each time step with a total cost of O(n+s). An optimal algorithm
16			15
Ii
e4 e2
I			I I			t
b3 b6 b1			e6 e1 b5 b2 e3
I			I
IA
e5 b4
13
1?
7,
14
2
56
would have time complexity O(n+k), so our algorithm, with its O(n+s) complexity
is efficient when s --H k is small. Assuming coherence, the number of exchanges s
necessary will be close to the actual number k of changes in overlap status, and
the extra O(s --H k) work will be negligible. Thus, for the one-dimensional bounding
box problem, the coherence view yields an efficient algorithm of extreme (if not
maximal) simplicity that approaches optimality as coherence increases.
5.3.2 The three-dimensional case
Efficient discrete view algorithms for solving the bounding box intersection problem
in R3 are much more complicated than the sweep and sort method for the one-
dimensional case. However, these algorithms all have in common a step that is
essentially a sort along a coordinate axis, as in the one-d'n?eusional case. Each
bounding box is described as three independent intervals [1)?X), c?r)], [b?j?? c?j??] , and
[b(iZ), e?z)] which represent the intervals spanned on the three coordinate axes by the
ith bounding box. Thus, our first thought towards improving the efficiency of a
discrete view algorithm for coherent situations would be to sort a list containing the
and e?x) values, and similarly for the y and z axes. Again, such a step will
involve O(n + s) work, where now s is the total number of swaps involved in sorting
all three lists. Ho??ver, if we observe that checking two J?ounding boxes fi?r overlap
is a constant time operation, it follows that if we simply check bounding boxes z
and j for overlap whenever an exchange is made between values indexed by i and
j (on any coordinate axis), we will detect all changes in overlap status in O(n + s)
time.
Again, we can maintain a table of overlapping bounding boxes, and update it at
each tin? step in O(n + s) tin?. The extra work involved is again O(s --H k). For the
three-dimensional case, extra work can occur if the extents of two bounding boxes
57
change on one coordinate axis without an actual change of their overlap status. In
practice, the extra work done has been found to be completely negligible, and the
algorithm runs essentially in time O(n t k).
Chapter 6
Local Motion Constraints
In this chapter, the basic motion constraint that prevents inter-penetration from
occurring near contact points is derived. Tlie motion constraint at one contact point
does not prevent inter-penetration at another contact point, and additional motion
constraints are required. For this reason, the motion constrahit at each contact point
is said to be a local constraint, as it is responsible for preventing inter-penetration
from occurring only in the immediate vicinity of the contact point, and only over
some (possibly small) range of motion. This locality property, and the fact that the
constraint is meant to apply only over some small range of motion is to be implicitly
?ssumed in the rest of our discussion.
NVe will start by expressing the non-peuctration constraint between bodies iii
terms of the spatial variables of the contacting bodies. By differentiating this
constraint twice with respect to time, a constraint on the accelerations of the bodies
is produced. Chapter 8 will show how workless constraint f??rces can be computed
that satisfy this acceleration constraint for each contact pohit. The derivation of
59
60
these acceleration constraints is well known for certain special cases of contact, but
as far as we know a formnlation for the general case has been lacking in the literature.
6.1 Configuration Space
The problem of preventing inter-penetration at any one contact point between two
bodies can be viewed as a more complicated version of the simple particle/s?irface
problem of Chapter 2, by considering the contact in terms of the configurationspace,
or C-space, of the two bodies[1,12,31,32]. The C-space of a rigid body is a six-
dimensional space, representing the degrees of freedom of the body. A single point
in this ?space is specified by six coordinates; each separate coordinate represents
one degree of freedom of the body. Thus, each point of this C-space represents a
particular configuration (a position and orientation) of the rigid body. The cLspace
of a pair of rigid bodies is 12 dimensional; a single point of this C-space represents
a configuration for each body. ?`e will describe the motion constraint between two
bodies at a single contact point using this 12 dimensional C-space.
Points in CLspace will be denoted by a horizontal bar, s?ich as p, while points in
world space (R3) will be written simply as p. Both points p in R3 and points p in
?space denote column vectors. Partial derivative operators a/Op and a/Op yield
column vectors when applied to fijuctions.
At time to, consider two bodies A and B that contact at Pc and whose
configuration is represented by the point 7r0. By considering the contact in C-space,
we can employ the standard technique of viewing the constraint that A and B
remain in contact near Pc as a manifold S, containing xo, in some neighborhood of
o. The set of all points in C-space near To that represent configurations such that
A and B remain in contact (near Pc) is a manifold, S, containing 0 Nioreover S
is a "hypersurface" in C-space; that is, 5 has an 1 1-dimensional tangent space and
61
Configuration Space			x?2
World Space
A			B
Pc
Configuration x1			Configuration xo
Configuration x2
Figure 6.1: Three configurations in world space and their representation as points
in ?space. The points r?i and 2 lie on opposite sides of the hypersurface S.
a single norn?al direction at each point.
The hypersurface 5 partitions C-space near To into two regions. Points near
To on one side of 5 represent configurations where A and B are disjoint near P?;
points near ? on the other side represent configurations where A and B have
inter-penetrated near Pc (Figure 6.1).
Since we are interested in the behavior of A and B over tin?e, let TT(t) be the
configuration of A and B at time t. Then ro = fflT(to). N\7e will describe 5 i'?plicitly
in terms of a scalar fi'nction C of ?space; that is,
5=fplc(p)=ot.			(6-1)
62
We will adopt the convention that if C(p) <0, then p represents a configuration
such that A and B inter-penetrate near Pc - Conversely, if C(p) > 0, then p
represents a configuration such that A and B are disjoint near Pc- Since A and B
contact at Pc at time t0, we know that
C(rr(t0)) = 0.			(6--H2)
In Figure 6.1, C is positive for configurations lying "above" S; thus, valid (non-
penetrating) configurations are those that lie on or "above" S.
Given the above definitions, the ?space representation of the non-penetration
constraint for an arbitrary contact has the same form as the non-penetration con-
straint for the point/surface problem of Section 2.3. To prevent inter-penetration,
rx(t) n?ist be prevented from lying below S; thus
>--H 0			(6--H3)
enforces the non-penetration constraint. Comparing this equation with Equa-
tion (2--H4), the analogy between the particle/surface non-penetration constraint and
the non-penetration constraint for an arbitrary contact is obvious.
6.2 Projections onto S
Equation (6--H3) expresses the non-penetration constraint in terms of spatial variables.
We convert the non-penetration constraint into a constrahit on acceleration by
differentiating. Consider
d
= ?dtC(T(t))
The flinction C(w(t0)) measures relative velocity in the same manner as
of Equation (2--H10). If C((t0)) < 0, then A and B are colliding (near Pc), while if
(33
C(?x(t0)) > 0 then A and B are breaking contact at P?; otherwise C(x(t0)) = 0 and
A and B are neither colliding nor separating at Pc. Since collisions are assumed to
have been dealt with prior to this stage, and the non-penetration constraint would
not be active if A and I? were separating at Pc, we can assume C(r(t0)) = 0.
Since c(rr(10)) = C(r?(t0)) = 0, inter-penetration is prevented by the condi-
tion
= d2
dt2C(T(to)) > 0.			(6--H4)
The quantity c(rx(t0)) is a measure of the acceleration between A and B at Pc, in
the normal direction. If C?(7T(to)) is positive, then A and B are breaking contact,
while if C?(--H(t0)) is zero, then A and B are remaining in contact.
Given any net forces and torques applied to A and B, these forces and torques
can be translated into a net C-space applied force acting on the configuration point
x0. If the contact between A and B is frictionless, then a ?space reaction force
normal to the hypersurface S at 0 arises (if necessary) to prevent fflT(t) from moving
downward.1 The ?space reaction force is computed by projecting the C-space
applied force onto the C-space normal, and reversing its direction (Figure 6.2).
Thus, it is assumed by those accustomed to working in C-space that computing
the reaction force, in either C-space or world space, is a simple matter. However,
the computation of the net C-space applied force can be difficult, because it must
take into account inertial forces generated by the motion of A and B. Consider
C(?x(t0)). Applying the chain nile,
= a0; (7r,(to))Tt(to).			(6--H5)
1 The normal direction is defined by the the natural inner product induced by the kinetic energy
of the system. The projection of Figure 6.2 must also be performed using the natural inner product,
but the intuition behind Figure 6,2 is valid. Erdmann[12J contains a more detalled discussion of
the subject.
64
C-space normal
Fapplled
il-'
Figure 6.2: The reaction force Freaction The reaction force is computed by
projecting the applied force Fapplied onto the C-space normal of S
(R?ecall that both ac/ap and t are column vectors.) Differentiating again,
d2			aac.			ac
dt2 C(wt(t0)) =			0p0p(r([0))T(tO)			(t0)+ `(7(to))Tt(to)
= t(t0)TO2C2 (r?(t0))#(t0) + ac
ap			ap (rr(to))Tt?,.(to)
with
a2c
(rr(t0))
ap2
the (symmetric) 12 x 12 matrix of second partial derivatives of C at i'(to). The
right term of the last sum in Equation (6--H6) is the component of acceleration of
x(t0) normal to the cLspace surface S. The left term of the last sum represents the
inertial forces acting on A and B. In order to ew?l'iate this term it is necessary
to calculate the second partial derivatives of C at To. This has been done only
for the case of polyhedral contact[12], but appears not to have been done for the
more general case of contact between curved surfaces. Thus, the computation of the
65
reaction force when either there is no initial motion (x?(t0) = 0) or the contacting
bodies are polyhedral is simply accomplished. In the next section, the second partial
derivatives of C are computed for the more general case of contacting surfaces in
motion. The derivation for the second partial derivatives of C yields a result that
is not particularly tractable for use in a simulation environment. As a result, in
Section 6.4, a more tractable expression for Equation (6 6) is derived, that does not
explicitly involve the computation of the second partial derivatives of C. However,
an analytical expression for the second partial derivatives of C may be useful in
some contexts, so a suitable expression is derived in the next section.
6.3 Calculating a2c/a?--H2
In order to create a constructive definition of C for curved surfaces, we will use the
concepts of the extremal points between two curved surfaces, from Section 5.2.1.
Since surfaces need not be convex in this chapter, the extremal points Pa and p? can
only be defined locally, in the vicinity of a contact point. Tlie extremal points near
Pc are defined as follows. If A and B are (lisjoint near Pc, then the extremal points
minimize the distance between A and B near Pc If A and B are in contact at Pc,
then the extremal points are Pa = PI) = Pc. If A and B have inter-penetrated near
Pc, then the extremal points maximize the distance distance between A and B (near
Pc). The points Pa and PI) at a contact point vary according to tlie configuration
p of A and B; thus, Pa and p? may be considered as f?inctions of configuration
(Figure 6.3).
Again, as in Section 5.2, we will model the sufface of A and B near Pc implicitly
by functions F and G; however, F and 0 will now be ftinctions of configuralion
as opposed to tin?. That is, if the configuration of A and B is p, then a point P
is on the surface of A if F(P,p) = 0, and similarly for I?. Again, A and I? can be
66
I			I
? Pa ,?
K
I			Pb
Figure 6.3: Change of the extremal points as a function of the configuration.
modeled equally well l?y parametric functions; Appendix A details the derivations
for the parametric case. Partial derivatives of F and C with respect to p are easily
computed, assuming partial derivatives of F and G with respect to P are known in
some rest frame.
As in Section 5.2, the functions F and G are such that the surface normal of
F
VF = OF
Op,
is outwards pointing, and similarly for G. Thus, the function C can be written
simply as
C(p) = VG(Pb,P)o+ (Pa --HPb) (6--H7)
where Pa and Pb are the extremal points for the configuration p (Figure 6.4). Note
that V&(Pb,P-) is always parallel to the vector Pa --H Pb. Differentiating,
OC			0			0Pa			0Pb
--HVC(Pb,P) (Pa --H Pb) + V&(Pb,T))			--H			--H .			(6--H8)
= Op			Op			Op
67
F
My
VG(p?,p)			VG(pb,p)
G			Pa
C(p)>O
C(p)=O			C(i?)<O
Figure 6.4: C expressed in tern?s of the extren?al points.
Then
a2c			=			aac
2(p)
Op			OpOp(?)
O2			a			Opa _ Op?
--H 0p2VG(p1),?)?(pa?ft)+0pVG(ft,p)			Op			Op
a
Op			tV&(?,p). Opa?Op?b(
+--HvG(p1),p)			Op			Op			Op2			Op
O2			-			a			-			Opa _ Op?
Op2 VG(p1),p) (Pa?P1))t20pVG4)1),P)			Op			Op
+V&(p1),p).
0Pa			01)1)
Op2 --H ap2
Clearly, to evaluate the second partial derivatives of C, we n?ist be able to
differentiate Pa and p? in tern?s of the configuration p. This requires a n?ore
constructive definition for Pa and p?.
In Section 5.2.1, a set of necessary conditions was given for the extrei?al points.
Because we are defining Pa and p? only locally in this chapter, the conditions of
Section 5.2.1 are both necessary and sufficient. That is, if we say that Pa and Pb
68
are the points near Pc that satisfy the conditions
E1: VF(pa,P) + A2VG(Pb,P) = 0
E2 :			F(Pa,P) = 0			(6--H10)
E3:			G(PI),T))=O
E4:			(PI) --H pa) + A1VG(pI),p?) = 0
then this uniquely defines Pa and PI) in some neighborhood of Pc for a			given
configuration p. Equation (6--H10) is modified exactly as in Section 5.2.2 for the
case of contact between a polyhedron and a curved surface. Constraint fi?inctions
for contact between two polyhedra are easily formulated[12,13] if the contact involves
either two non-parallel contacting edges, or a vertex in contact with the interior of
a polygonal face; both these cases admit a well defined normal direction However,
contact between two vertices, a vertex and an edge, or two parallel edges do not
yield a well defined normal and may be considered as degenerate contact situations.
Palmer[42] discusses the complexity of simulations involving degenerate contacts in
general, degenerate contacts yield NP-hard simulation problems. NVe will assume
that a normal direction is defined at all contact points.
NNre can use Eq?iation (6--H10) to differentiate Pa and p? with respect to p.
Suppose we are given values for Pa and p? for some configuration p. (Typically,
we wish to differentiate Pa and p? for the coiifig?iration 7T?; fi?r this configuration,
Pa = PI) = Pc.) Given Pa and p?, values for A1 and A2 are easily calculated. Using
Equation (6--H10), we can calculate opa/ap and &p?/?p by making use of the implicit
function theorem for simultaneous equations[49]. This theorem asserts that under
proper conditions, Pa and PI) may be regarded as functions of the configuration p?;
the theorem also gives a constructive expression for the derivatives of Pa and PI).
NVe will write the Jacobian determinant of a set of vector flinctions H1(x1Tn)
OH1			OH1
Ox1
OHn			OHn
Ox1			OX?
69
through Hn(xi,... ,xn) as
O(H1,...,Hn)
O(xi,. .. ,xn)
Since Pa and Pb are points in R3, Pa and Pb are each three (scalar) functions of
configuration: Pax(P), Pay(T)) and paz(p?) and similarly for Pb The implicit ft?nction
theorem, applied to Equation (6 10) yields
0Pa			A			((3--H11)
where
K = 0??,p??pE?1?f1??$?)$, A1, A2)
(6--H12)
J =			O(E1, E2, E3, E4)
O(Pax,Pay,Paz,Pbx,Pby,Pbz,A1,A2)
and p? denotes the ith component of p. Similar expressions give the derivatives for
Pay Paz, Pb?, Pby and Pbz' The case when J is near-zero or zero is considered in
Section 6.7.
Second partial derivatives for Pa and Pb can be calculated by differentiating
Equation (6 11). For example,
02Pax			0			A
0Pi0Pj			Oi?)j			J
(6--H13)
Combining this result with Equation (6--H9) yields an analytical expression for the
second partial derivatives of C. Examining the dimensions of the equations involved
however, both A and J are the determinants of 8 x 8 matrices. Differentiating A
and J with respect to p requires computing the derivative of an 8 x 8 matrix deter-
minant. This is extremely impractical. Although A and J have considerable block
70
structure, block structure cannot be exploited in computing matrix determinants.
In its present form, the Jacobian determinants of K and J contain more than 1,000
terms of seven factors each; the derivatives contain far more terms. Additionally,
each of the 72 different second partials of Pax, Pay, Paz, Pbx, Pby, and Pbz would
need to be computed in this manner. Thus, although we have achieved an analytical
expression for Equation (6--H6), it is not practical.
6.4 Computing C?(x(to)) directly
In this section, we shall formulate C in a slightly different manner. Let us assume
a rotated coordinate system in which VF(pa, rr(t0)) and V&(Pb, ?(to)) are colinear
with the z axis, with VG(Pb, rr(t0)) directed positively along z. (??7e will employ the
standard right-handed coordinate system used to depict functions z = h(x, y), with
the vertical axis.) Next, we explicitly model A and f3 near Pc as scalar functions
f and g of r and y and configuration p. NVhere F was a function F(.'r,y,z,p), f
is a function f(x, y, p) and similarly for g. The justification for the existence of f
and g is the implicit function theorem. The formal definitions of f and g are given
by
F(x,y,f(x,y,p),p) = 0			and			G(,y,g(x,y,p),T?) = 0.			(6--H14)
Thus, the point (x, y, f(x, y, r?(t))) is a point on F at time t, and similarly for C.
In this new coordinate system, the extremal points Pa and Pb share the same x and
y coordinates at time to. At times t near t0, A and B do not penetrate as long as
the z value of Pa is greater or equal to the z value of Pb (Figure 6.5). ?Ve will use
this to redefine C as
C(p) = f(Pa1,Pay,T)) --H g(Pbx,Pby,P)
where Pax, Pay? Pbx, and Pby are the extremal points for configuration p.
(6--H15)
71
z =f(x,y)
(Pax `Pa? `f(Pax `Pa?' t))
(Pbx?Pby `?(pbx'Pb? ,t))
z =g(x,y)
Figure 6.5: Side view of the implicit surfaces F and G expressed explicitly by
functions f and g.
Since the z coordinates of the extremal points are no longer needed, let us rewrite
Equation (6--H10). Differentiating Equation (6--H14) we obtain
OF			OF			OF
Of _			Ox			Of _			Og			Of _			Op
Ox			OF'			Oy --H			OF'			Op			OF
Oz			Oz			Oz
and			(6--H16)
and similarly for g. Second derivatives of f and g are obtained by repeated
differentiation of Equation (6-14). Using f and g, condition E4 of Equation (6--H10)
may be written componeutwise as
Pbx			Pax
Pby			--H			Pay
g(Pbx,Pby,P)			f(Pax?pay,P)
+ A1
OG
Ox
00(Pb,P)
Oy
O&(ft?p)
Oz
f(Pax,Pay,P)--Hg(Pbx?Pby,P)
A1 =			OG
Oz
0
=			0			, (6--H17)
0
from which we obtain
(6--H18)
72
Multiplying Equation (6--H17) by --H1 and using Equation (6--H16) allows us to rewrite
E4 as
Pax Pbx
+(f(Pax,Pay,P?)?g(Pbx,Pby,P))
Pay --H Pby
ag
?(Pbx?Pby?P?)
a? x,Pby,P)
= 0. (6--H19)
This new condition has one less eQiation than E4 because of the reduction of
variables from F to f. In a similar fashion, A2 of condition E1 is eliminated
and condition E4 is rewritten as
???(pax,pay,p)			?(Pbx?Pby?P?)			= 0.
af			ag
0?(Pax?Pa??P)			?y(Pbx?Pb??P)
(6--H20)
Conditions E2 and E3 arc no longer required since f and g give explicit definitions
of the surfaces. At this point, it will be more useflil to treat the extremal points as
functions of time. We can define Pax, Pay? Pbx and Pby as the extremal points at
time t by the system
D1:
af
Og
a? ax,Pay,T%(t))			a? (Pbx, P?y,T(t))
=0
af			???p? fty,T(t))
?y(Pax?Pay?TY))			ay
D2 : (f(Pax,Pay,r'(t)) --H g(Pbx,Pby,T(l)))
ag
--H(Pb
a?			x,Pby,?r(t))
Og
?y(Pbx?Pby?TT(t))
Pax?Pbx			_
+
Pay --H Pby
(6--H21)
73
Now Pax, Pay? Pbx and fty can be differentiated as functions of time. The
implicit function theorem for simultaneous equations gives the result
a(D1, D2)			a(D1, D2)
Pax = --H 0(1,Pay,Pbx,Pby)			and			1)ay = --H ?(Pax, t,Pbx,Pby)			(6--H22)
J			J
where
and
O(D1,D2)			6--H23
?(Pax,Pay,Pbx,Pby)
af			af			(6-24)
at(Pax?Pay?T(t)) = ap(Pax?Pay?rT([))T#([)
(and similarly for g). Similar expressions are obtained for 1)bx and 1)??. The deriva-
tives of Paz and Pbz are simply af/at and Og/Ot. Equation (6--H23) redefines J, but
we will not have any flirther use for the 8 x 8 Jacobian matrix in Equation (6 12).
?Ve will consider the case when J is near-zero or zero in Section 6.7.
Since D1 and D2 do not involve Paz or Pbz, let Pa and Pb henceforth denote the
column vectors (??X,??y)T and (Pbx,Pby)T at time t. Similafly, let 1)a = (1)ax,1)ay)T
and Pb = (1)bx,Pby)T be the coliin?n vector denvatives of Pa and Pb at time t. The
quantities ?f/?p, af/? and i?(t) are also column vectors. From Equation (6--H15),
the condition c(r?(t0)) > 0 is simply
d2			d2
?t2f(Pa?T(to)) --H ????Y???%(tO)) > 0.			(Q-25)
Differentiating f with respect to time,
d			Of(Pa,(to))T1)a + Of			(6--H26)
wtf(Pa?%(to)) = 0P			Op (pa,x(to))Tt([o).
74
Then
d2
?t?f(Pa?x(to)) =
Of			rx(to))pa + O2f			?(to)?x(to) T
0?2(pa,			?p?p(Pa?
Of
+ (pa,x(to))p'a
+ 0))(Pa, x(t0))Pa + Of2 (Pa, x(t0 ))t(t0) Tt,(to)
Op
(6--H27)
.T0f			.			TO2f
= Pa 0p2(Px?T(tO))Pa + 2x(t0) 0p0p(Pa?%(tO))Pa
+ Of			T0f
0p(Pa?wT(to))TPa +t(t0) Op2t([o)
+ Of
??(pa,rT(to))TY(to).
OP
Fron? Equation ((3--H16) and the fact that ?F(pax,pay,x(to)) is colinear with the z
axis
This yields
(12
?t2f(Pa??T(to)) --H
Of			0
(Pa,WT(to)) =
OP			0
o+T Of			T O2f
Pa 2 (??, T(to))pa + 27L(t0) (P Y(t0))?a
OP			0POP
(6--H28)
o+			TOf			Of..
2t(to)+--HWT(!o).			((3--H29)
+rv(t0) Op			Op
A sii?ilar expression is obtained for g. Thus, neither pa nor i;b is required to
calculate c?(r?(t0)). Equation (6--H29) depends linearly on tt(to), which in turn
depends linearly upon the forces acting on A and B; this equation can be used
to calculate the reaction force between A and B. Since Pa and Pb are simply
calculated as the ratio of two 4 X 4 matrix determinants, the final expression is
quite practical for use.
75
6.5 Predicting collision times
The primary use for the derivation of C(x(t0)) is to formulate the equations for
the constraint forces at contact points. However, C?(c(t0)) can also be used to
help predict the time of a collision. When the collision/contact determination step
detects inter-penetration between two bodies, an estimate must be made of when the
two bodies collided. The work done by the simulator can be significantly decreased
by increasing the accuracy of the prediction of the collision time.
Consider the case when inter-penetration is detected between two bodies at time
t1. Let C(r(t)) represent the non-penetration constraint that has been violated in
the vicinity of the inter-penetration.2 Since the constraint is not satisfied at time
t1, we have C(wr(t1)) < 0. However, if the previous time step (before the inter-
penetration) occurred at time to, then C(T?(t0)) > 0. Computing the collision time
is thus a root finding problem, on some known interval: we are looking for a time
tc, wfth to ? tc ? ti such that C(rr(tc)) = 0.
Since we can easily calculate derivatives of first and second derivatives of C, we
model C(7??(t)) quadratically near to as
= C(r?(t0)) + C(r?(t0))(t --H t0) + C(r?(t0))(t --H t0)2 + O((t --H t0)3). (6--H30)
Ignoring the O((t --H t0)3) term, we solve for C(rr(t?)) = 0 to obtain the estima-
tion			__________________________
tc = to +
--HC(wr(t0))I?C(r(t0))2--H2C?(7T(t0))C(x(t0))
C?(TT(t0))
(6--H31)
2 ?? is not always clear exactly whidi constraint this should be. For example, if a vertex from
polyhedron A is inside polyhedron B, whid? face of B did it pass through to get there? In cases like
this, the most likely non-penetration constraint should be chosen. For most situations, the correct
non-penetration constraint is easy to determine. If several non-penetration constraints appear
to have been violated, for example, if n?iltiple vertices of A lie inside B, the non-penetration
constraint that appears to have been violated the earliest should be d?osen. For situations involving
multiple constraint violations, an incorrect choice is not fatal at worst, it will fail to yield a good
estimate for the collision time.
7(3
The estimate Ic should be the smallest real root of Equation ((3--H31) greater than
t0; if no such root exists, or C(x(t0)) is near zero, tc can be estimated using a
linear approximation of C(x(t)). The method is made robust by incorporating a
bisection step when convergence is slow[15,44]; that is, predicting tc = (t0 tt1)/2.
For constant acceleration, Equation (6--H31) gives an exact result for tc, if C(?x(t0))
is a linear measure of distance.
6.6 Rolling contact
As we have noted, C?(r?(t0)) is a measure of relative acceleration normal to the
contact surfaces of A and B. In Chapter 9, when friction is discussed, we will need
a description of both the relative slip velocity between A and B tangent to the
contact surface, and the derivative of the slip velocity with respect to time. These
quantities are expressed more naturally in world space than configuration space.
The derivation of these quantities is straightforward for bodies of any geometry,
and is certainly well known[4,41], but we will derive it here for completeness. NVe
will also show why the constraint of pure rolling motion between bodies of arbitrary
geometry is so much simpler to express than the non-penetration constraint we have
just developed.
Let Va, Vb, Wa and Wb denote the linear and angular velocities of A and B at
time to. Let ca and cb denote the centers of mass of A and B. Then Ca = Va and
Cb = Vb. Assun? that A and B do not break contact near Pc. The slip velocity, s,
between A and B at time to is therefore
5 = (Va + Wa X (Pc --H Ca)) --H (vb + W? X (Pc --H
(6--H32)
77
The derivative of the slip velocity is
d			d			d
?dtS = TtVa + ?dtWa X (Pc ca) + WaX (Pc --H Va)
(?dtVb+ftwb x (Pc?Cb)+Wb X (pc?Vb)). (6--H33)
Assumingthat s = 0, that is, that A and B are not initially slipping, A and B
can be made to maintain s = 0, and thus roll without slipping, by enforcing
d
Tts = 0.			(6--H34)
Since p? isfairly easily calculated (by the previoussection), the constraint that two
bodies roll without slipping is easily derived as aconstraint on the forces acting on
A and B. Notethat thisconstraint automaticallyenforcesnon-penetrationbetween
A and B. In orderfor A and B tointer-penetrate, they must attainsome non-zero
relative velocity normal to the contact surfaces; the constraint wdts = 0 prevents
this fromoccurring. ?NJliyshould this rollingwithout slippingconstraint be son?ich
sii?pler to express than simply the non-penetration constraint?
Themotionofrollingwithoutslippingbetweentwobodiesisahighlyconstrained
motion. Equation (6--H34) enforces three constraiiits on the relative acceleration
while Equation (6--H4) enforces only one constraint. Thus, rolling without slipping
has fewer degrees offi?eedom than simply moving without inter-penetrating. As a
result, we have more knowledge abo?it the motions of A and B, and the constraint
is formulated n?id? more easily.
6.7 Ill-conditioned Jacobians
As stated in Chapter 3, impact is assumed to have been dealt with prior to
the computation of constraint forces. However, what has really been dealt with
are situations where A and D have a non-zero approach velocity at Pc normal
78
to the contact surfaces. This is not the only case in which impact can arise.
Featherstone[13] gives an example where discontinuities in the curvature of the
contact surfaces require an impulse to be applied, even though the approach velocity
normal to the contact surfaces is zero. Given the equations developed to compute
C for curved surfaces, we can add yet another type of impact that might occur
between bodies.
Suppose we wish to simulate the frictionless oscillation of an ellipsoid A upon
a plane surface B (Figure 6.6). Then the flinction F describing A is a simple
quadratic function with non-zero curvature everywhere. B is described by a linear
function G and has zero curvature everywhere. By solving Equation (6--H25) for the
reaction force between A and B, and considering Equation (6 29), it can be shown
that for given initial conditions, the reaction force magnitude 1? is of the form
R(t) = O(IIi)c(t)II + 1)
(6--H35)
over the simulation, treating Pc as a fiinction of time.
Recall from Equation (6--H22) that Pc is calculated as the ratio of two deter-
minants. Let us calculate the denominator determinant J at a time to when A
and B are in contact (and thus f(P??x,Pay,T(to)) --H g(p?x,P?y,7?(to)) = 0). Eval-
uating Equation (6 -23), and since C (and thus g) have second partial derivatives
everywhere,
fii			f12			0			0
= f12			f22			0			0 = f11f22 --H 42 =
1			0			--H1			0
0			--H1			0			--H1
(6--H36)
where
af			af
fii = ax2(Pax?Pay?T(to))? f22 = ay2(Par?Pay?T(to))?
(6--H37)
79
b			c			d
500
400
R(t) 300
200
b			d
a
0			0.2			0.4
Time
0.6
Figure G.(3: Frictionless oscillation of an elli1)soid A on a plane B.
80
and
O2f
f12 = OxOy (pax,pay,x(to)).			(6--H38)
This gives J = ic where K is the curvature of the ellipsoid at the contact point Pc
For the case of an ellipsoid, the curvature ? is always positive, so J is always non-
zero. However, we can make ? arbitrarily close to zero by scaling the ellipsoid as
shown in Figure (3.7. As the curvature at the contact point decreases, the reaction
force increases, in a continuous fashion. The sin?ilator n?ist take small steps in
order to accurately sample how the reaction force R changes over time. However,
this is merely a practical problem of simulation. Niathematically speaking, nothing
unusual occurs.
Suppose however that we enable A to have zero curvature at certain points. NVe
can do this by letting A be a s?iperquadric ellipsoid, defined in some rest flame as
the set of points p satisfying
PttP4y+P4z = 1.			(6--H39)
Again, we will let A undergo frictionless oscillation in contact with B. Suppose
that a point of zero curvature on A comes into contact with B at time to. As t
approaches to, the value of J approaches , in a contin'ious flishion. As a result,
R(t) increases in a continuous fashion without bound, as t approaches to In the
limit, the reaction force R(t) approaches infinity as the contact point approaches
the zero curvature spot of A (Figure 6.8). Should this be considered as an impulsive
reaction between A and B?
The answer to this is "alI?ost". Clearly, the reaction force between A and B
at to cannot be considered to be a force, since it is infinite. However, in order
for the reaction force to behave impulsively, it must cause a velocity discontinuity
at time to. Even though the reaction force is unbounded, it does not cause such
a discontinuity in velocity, although it does cause a discontinuity in acceleration.
81
b			c			d
R(t)
2OOO?
1000w
a
OK .
T
0
c
b			d
v
0.1			0.2			0.3			0.4
Time
Figure 6.7: Frictionless oscillation of a more eccentric ellipsoid A on a plane B. The
reaction force J?(1) is stronger and more sharply peaked the closer the curvature ?
of A approaches zero,
82
a			b
c			d
6000
R(t) 4000
2000
a
0
b			c
0			0.2			0.4
Time
Figure (3.8: Frictionless oscillation of a superqnadric ellipsoid A on a plane B. The
reaction force i?(t) diverges as the contact point approaches a zero curvature spot
on the superquadric ellipsoid A.
83
That is, the reaction force at time t0 is an infinite force acting for an instant of
time, but in such a way that the measured impulse is zero. To see this, consider the
integral of the reaction force R over time. In order for the reaction force to behave
as an impulse, it is necessary that
to+?t
lim
At?o t0--HAt R(t)dt>O.
Since 1?(t)			O(IIpc(t)Ii + 1), consider
t0+At
hm			O(IIpc(t)II + 1) (it.
Ai?o t0--HAt
(6--H40)
Since this is an integral over time involving i)c(1), by the fundamental theorem of
calculus there exists a continuous function H of Pc such that
it%tt
O(IIj)c(t)II + 1) (it = 1I(pc(t0 + At)) --H H(pc(t0 --H At)). (6--H41)
Since pc(t) is a spatial measure (the position of the contact point at time t), pc(t)
is clearly continuous. Therefore,
t0+At			t0+Ai
lim			R(t)(it =			lim			O(Ii)c(t)II + 1)dt
At?o t0--HAt			At?O t0--HAt
--H Al?WWo H(1)c(t + At)) --H H(pc(t --H At)) (6 42)
--H 1I(pc(t0)) --H ]I(pc(to))
--H 0.
Thus, the reaction force is not impulsive. Viewed in another way, the reaction force
behaves in the limit as an impulse of zero magnitude. For want of a better term,
we could call such behavior a "zeroeth order" impact.
In practical terms, the problem of a zeroeth order impact is how to evaluate the
integral of the force over a time period containing the impact. Current ordinary
84
differential equation solvers are not well suited for evaluating this integral, even
though it is finite. From a sampling viewpoint, the reaction force must be sampled
infinitely often near time to to accurately evaluate the integral of the reaction force.
One approach is to analytically evahiate the integral in a region of time bounding
1o, but it is unknown how to do this efficiently in a simulation environment. The
s'n?ilation process described in this thesis does not deal with zeroeth order impacts.
Chapter 7
The Penalty Method
The remainder of this thesis will be devoted to the problem of computing the forces
that arise at contact points for systems with n?iltiple contact points. In this chapter,
we will examine the penalty method for computing constraint forces between bodies
in the absence of fliction. Since we already know that the penalty method usually
performs poorly in practice, why bother to consider it at all? The main reason for
considering this method is that in Chapter 12, we shall consider a model of frictional
behavior that is based on the penalty method. Prior to this, we need to consider
whether the penalty method for frictionless systems yields physically correct motion.
Validating the penalty method for the frictionless case does not necessarily validate
the model of frictional behavior proposed in Chapter 12, but it would clearly be
unreasonable to extrapolate from the penalty method if it were known to be flawed
for the case of frictionless systems. A secondary reason for considering the penalty
method is to gain insights into its correctuess as a si'?ulation method; considering
the large amount of simulation work that uses the penalty method, this is a topic
worth pursuing.
85
86
7.1 Equality constrained motion
To simplify matters, we will first consider the penalty method as applied to equality
constrained dynamics problems. This is actually not a large restriction, because sys-
tems with inequality constraints can be viewed as being subject to the corresponding
equality constraints on any time interval during which no constraint changes from
active to inactive or vice versa. For example, in the particle/surface example of
Section 2.5, the motion of the particle, whenever it is in contact with the surface,
is exactly the same as if it were constrained to always lie on the surface. Thus, we
would first like to study the motion produced by the penalty method for equality
constrained problems, such as the problem described in Section 2.4.
Rubin and Ungar[45] give an excellent proof that the penalty method does indeed
yield (in the limit) physically correct behavior for flictionless systems with equality
constraints. Ho??ver, their proof assumes that the external forces acting on bodies
are due to the potential energy of the system, rather than being arbitrary flinctions
of time. Niodifying Rubin and Ungar's proof to handle arbitrary time-dependent
external forces does not seem overly difficult, but their proof is fairly lengthy and
somewhat complicated, and proves a n?ich stronger result than we are interested
in. ?te will show that the penalty method produces (in the limit) physically correct
motion by comparing differential descriptions of the physically correct motion and
the motion produced by the penalty method. This is much simpler and more
intuitive than Rubin and Ungar's approach; note however that Rubin and Ungar
obtain a much stronger result about the limiting motion of the penalty method than
can be achieved using the simple differential comparison.
Let us consider a system with equality constraints from time to to t1. In
the previous chapter, the non-penetration constraint at a single contact point was
represented by a hypersurface in C-space, and the configuration r?(t) was required to
87
lie on or "above" this hypersurface. In restricting ourselves to equality constraints,
we are requiring ?x(t) to lie on the hypersurface. For a system with k bodies and n
contact points, the n non-penetration constraints are represented as n hypersurfaces
in a 6k-dimensional ?space. The configuration point x(t) must lie on all n
hypersurfaces. Equivalently, ?(t) must lie in the intersection of all the constraint
hypersurfaces. Since the constraint hypersurfaces are manifolds, the intersection of
the constraint hypersurfaces is itself a manifold, which can be expressed in terms of
a scalar function C1(p) as C1(p) = 0. We shall call this manifold the intersection
manifold.
7.1.1 The physically correct motion
We will use Lagrangian dynamics to derive the equation of motion for a system
constrained to satisfy c1(rv(t)) = 0 from time to to ti. The remainder of this
chapter assu'?es a knowledge of Lagrangian dynamics and differentiable manifolds.
In order to describe the motion of the system (that is, the fLinction fflT(t) from to to
t1), we will have to describe the configuration state 1(t), and velocity t(t), in terms
of a coordinate system on the intersection manifold. Let the intersection manifold
have dii?ension d, and let S(q) be a parameterization of the manifold. The function
S(q) invertibly maps points from some region of Rd to the intersection manifold.
(If the manifold cannot be covered with a single parameterization, the motion of
1(t) can be considered in suitably small intervals.) The motion of 1(t), restricted
to lie on the intersection manifold is then described in terms of a continuous path
q(t) in Rd. Given initial values for 1(t0) and 1(t0), the initial conditions q(t0) and
q(to) must satisfy
??S(q(to)) = dSq(t0)q(to) = 1(t0)
S(q(t0)) = 1(t0)			and
(7--H1)
88
where
dSq(t0) = as			(7--H2)
The function dSq(t) maps vectors from Rd to vectors that are tangent to the
intersection manifold at S(q(t))
We also need to describe the kinetic energy of the system. Given a system with
configuration p in ?space and velocity v, the kinetic energy T of the system can
be written in the form
T = 21vTM(T))?v			(7-3)
where M(p) is a square matrix with size equal to the dimension of C-space. The
matrix M(p) is symmetric positive definite. The kinetic energy is used to induce an
inner product on C-space; given vectors ti and w at point p of ?space, the inner
product of tt and w is written ? tt, `?v >p and is defined as
?ti,tV>?p= ?Thtt4-=?JM(p)?w.
Using the inner product, T can be written as
T = -21<Thv4.
(7--H4)
(7--H5)
Expressing the kinetic energy T in terms of the parametric coordinates of WL'(t),
we have
=?2l<)?s(q(t)),			--Hdts(?(t))>s?qu?? (7--H6)
--H			`qt'qt
(Henceforth, all evaluations of dS will take place at q(t), so we shall simply write
dS in place of dSq(t). Likewise, all inner products are assumed to be evaluated at
89
S(q(t)) unless otherwise qualified.) Given any external CLspace force F(t) acting on
the system, let Q(t) be the equivalent generalized force on q(t). If we assume that
the potential energy of the system is always constant, since Q(t) can include the
effects of forces due to potentials, the equation of motion for q(t) is given by
Wt aq' --H aq			(7-7)
Using Equation (7--H6), this yields
ft a$ (21?dSq(t),dSq(t)>) --H ; (-21<dSq(t),dsq(t)>) = Q(t)
with initial conditions for q(to) and q(to) as specified by Equation (7--H1).
(7--H8)
7.1.2 Motion due to the penalty method
Next, we will derive an equation of motion for the penalty method, which will then
be compared with Equation (7--H8). As defined in Section 2.4, the penalty method
works by imposing a penalty force on the system when r?(t) deviates from the
intersection manifold. The penalty force can be viewed as the gradient of a potential
energy function P(p) that increases the farther p lies from the intersection manifold.
Thus, the penalty force is conservative, and does not change the overall energy of
the system. If the potential energy function P is made steep enough, then for a
conservative system, rr(t) can never deviate more than a certain distance from the
intersection manifold; that distance is determined by the total energy of the system.
Thus, for the case when the only external forces acting on the system are due to
the potential energy of the system, increasing the penalty force (by increasing the
steepness of P) binds ?x(t) closer to the intersection manifold.
However, if time-dependent, and thus, non-conserw?tive, external forces act on
the system, increasing the slope of P does not automatically limit the deviation
of ?r(t) from the intersection manifold. If the external force varies over time, it is
90
possible that a steeper P may cause a "resonance", and allow ?x(t) to deviate more
from the intersection manifold. For example, the amount of work done by a given
external driving force on a linear spring system is not monotonic with respect to the
spring constant. In general though, as the spring constant is increased, the driving
force does less work. Similarly, in analyzing the penalty method, we will assume
that the maximum deviation of x(t) from the intersection manifold can be made
arbitrarily small, by increasing the steepness of P in such a way that resonance is
avoided. We would like to show that as the maximum deviation decreases, the path
traced in the vicinity of the intersection manifold converges to the physically correct
path lying exactly on the intersection manifold.
In the previous section, the position of r?(t) on the intersection manifold was
described by the parametric coordinate q(t). For the penalty method, we need
to describe ?x(t) in terms of a coordinate q(t) on the intersection manifold, and a
displacement from the intersection manifold. To accomplish this, ?? will describe a
coordinate system for the orthogonal complement of the intersection manifold. Let
the orthogonal complement of the intersection manifold have dimension k. At each
point p on the manifold, we define vectors ]Vi(P) through ATk(p) such that
and
?Ni(P),Ni(P)>p?= 1			(1 ? i ?
<Nj(?), Nj(P)4 = 0			(1 < i ? j < k).
If u is tangent to the intersection manifold at p, then
<Nj(fl), U>p?=0			(1 < i < k);
(7--H9)
(7--H10)
(7--H11)
that is, each Nj(P) is normal to the intersection manifold at p, as defined by
the inner product induced by the kinetic energy. Each Nj(P) varies smoothly
91
over the intersection manifold; since we are considering the motion of x(t) locally,
orientability of the intersection manifold is not a concern.
Given this orthogonal coordinate system, we will express points ?x(t) lying near
the intersection manifold using the parametric coordinate q(t) and a k-tuple of real
numbers (h1(t), . . . , hn(t)) as
k
= S(q(t)) + ?
i=1
(7--H12)
This coordinate system is one-to-one in some region of ?space properly containing
the intersection manifold, as long as the curvature of the manifold is bounded.
Since we can force x(t) to lie as close to the intensection manifold as necessary,
Equation (7--H12) is an acceptable way of defining the position of ;t([)
Using the notation of Equation (7--H12), it is simple to write an equation of motion
for the penalty method. Since the values hi(t) describe the displacement of rr(t)
from the constraint manifold, we can write the potential energy of the penalty
function in terms of llj (t). The assumption is that for any E > 0, a potential energy
function P(h1(t),. . . ,hk(t)) can be chosen1 so that 111i(t)I < e for all 1 < i < k
and to < t < ti.
The kinetic energy ?? of the system, as in Equation (7 G), is given by
From Equation (7--H12),
= -21<#(t),#;(t)>.			(7-13)
k			k
t(t) = dSq(t) + ? h??(t)N?(S(q(t))) + z hj(t)ATj(S(q(t)))
i=i			i=i
where N?(S(q(t))) = ?d?Nj(S(q(t))). Then the kinetic energy can be written as
(7--H14)
1 This does not mean to imply that such a d?oice is simple it's not. The choice of a suitably
strong penalty forces usually requires great care. This is simply one of the penalty method's many
problems.
92
= -21<dSq'(t) +
hj(t)ATj(S(q(t)))
dSq'(t) + ?k1hj(t)Nj(s(q(t)))
?21 (?dsq(t), dsq(t)>
+
+
%ki hi(t)Ni(S(q(t)))
k			k
+ ?? hj(t)I\Tj(S(q(t))), ?
k			k
+ <? h?(t)&?(S(q(t))), ? hi(t)Ni(S(q(t)))>)
z=1
(7--H15)
k k
+ <dSq(1), ? h?(t)N?(S(q([)))> + <d5q(t), ?
i=1
k
+ <? hj(t)ATj(S(q(t))), ?
i=1			i=1
NVe can sin?p'ify this expression for T' by noting that the vectors Ni(S(q(t)))
through N?(S(q(t))) are orthonorn?al with respect to the inner product. Thus
k			k			k
<? h?(t)N?(S(q(t))), ? h?(i)N?(S(q(t)))> = ?
i=1			i=1			i=1
(7--H16)
Also, since the vector dSq(t) is tangent to the intersection n?anifohl at q(t) for any
q(t), the inner product ?dSq'(t), A)(S(q(t)))> is aLso zero for any 1 < i ? k and
any q(t). Thus,
k
?dSq(t), ? hj(t)ATj(S(q(t)))> = 0.
i=1
Using these facts, we n?ay rewrite Equation (7--H15) as
21(<dsq(t),dsq(t)> + kjli(t)2
i=1
k			k
+ <? h?(t)N?(S(q(t))), ?
i=1
(7--H17)
(7--H18)
93
k
+ <dsq(t), ? hj(t)N?(S(q(t)))>
k
+ <? hj(t)N?(S(q(t))), ?
i=1			i=1
The equation of motion for q(t) is
--Hdt )(Tt?P)			$(Tt?P)=Q?			(7-19)
with the same initial conditions of the physically correct motion specified by
Equation (7--H1). Since the penalty energy function F is independent of both q(t)
and q(t), we can express Equation (7--H19) as
kit $?? --H aT'			(7--H20)
aq
This equation is independent of the penalty force --HVP; as expected, the penalty
force imparts no "tangential" acceleration to w?(t), along the intersection manifold.
All that remains is to show that the motion of q(t) according to Equation (7--H19)
is the same as the motion according to Equation (7--H8) in the limit as hj(t) 0
for each hi(t). Since the physically correct motion and the penalty method motion
both have the same inftial conditions, all ?? need to do is show that as the h1(t)s
converges to zero, Equation (7--H19) and Eq?iation (7--H8) yield the same differential
equation
First, consider aT?/aq Since we are not differentiating with respect to hj, any
terms containing hj(t) drop out as hj H 0 Differentiating Equation (7--H18) with
respect to q?, all that remains is
aT'
aq
$. --H21<dSq(t), dSq(t)> + k
? hj(t)
--H;q. (--H21<dsq(t) dsq(i)>)
(7--H21)
94
Now, consider ?T'/?q. Again, any terms containing hj(t) drop out, and we
have
aT' --H a --H21?dSq(t), dSq(t)> t
aq			aq			? hj(t)2
= --H)i (21<dsq(t),dsqU)>)
Thus, as the hi(t)s converge to zero, the equation of motion
--Hdt %`			aT' -
becomes
(7--H22)
(7--H23)
Wdt a$j (-21<d54(i),dsq(t)>) --H ; (21<dsq(t),dsq(t)>) = Q(t). (7--H24)
Since this is precisely the same as Equation (7--H8), the motion due to the penalty
method converges to the physically correct motion as the maximum deviation of i't(t)
from the intersection manifold converges to zero. The proof by Rubin and Ungar
establishes that the normal velocity, which is described by iii (I) through hk(t) in
our parameterization, also converges to zero as the steepriess of the penalty function
P is increased.
7.2 Inequality constrained motion
The penalty method yields a correct result (in the limit) on time intervals during
which no constraints A?ange. However, what happens at the instant that a constraint
does change? Specifically, how does the penalty method know when to "let go" of
a contact point, and allow the contact to break? (The converse of this involves a
constraint activating, which results in an impulse, but ?? are assu'?ing that impulses
have already been dealt with.) For non-penetration constraints, the penalty force
95
for a given non-penetration constraint is set to zero if the constraint has not been
violated; this prevents penalty forces from becoming attractive. Unfortunately1 this
does not always give the correct effect, because in some situations the penalty force
must act attractively to correctly maintain contact. Thus, one cannot simply release
a constraint when the penalty force would have become attractive. It would be nice
to speculate that the normal accelerations, that is, h'1(t) through 11k(t), converge
to zero as the deviation from the intersection manifold goes to zero, If that were so,
then it is possible that, given a very small step size At, examination of I'i'i (1) through
h"k(t) would indicate whether a constraint should be made inactive; presumably, an
increasing valtie of h"j(t) would indicate that the ith contact point should be released.
Unfortunately, Rubin and Ungar give an example that shows that in general h"1 (t)
through h"k(t) do not converge to zero.2
Thus, the penalty method must be used in conjunction with an algorithin that
decides when constraints should be released. The algorithms in the next few chapters
furnish precisely this information; but since they also calculate the constraint forces
as well, they are clearly a better simulation method than the penalty method. The
frictional model extrapolated from the penalty method in Chapter 12 is based on
the behavior of the penalty method over an interval in which constraints do not
change state.
2 Actually, they doii't converge to anything.
Chapter 8
Frictionless Systems
As we saw in the previous chapter, the problem of simulating non-penetration
constraints on intervals when there is no change in constraint status is an equality
constrained problem. Constraint forces for equality constrained problems can be
computed by solving a system of linear equations. However, for the general problem
of computing the constraint forces at time t0 for a system with n contact points, we
must allow for the possibility that some set of the constraints may become inactive
after time to. The inability of the penalty method to determine which constraints
should become inactive at time to is a major shortcoming of the method.
Let us consider a system of n contact points without fliction. At the ith contact
point, the unit surface normal between bodies is denoted ??. If this contact occurs
between bodies A and B, and 1?j points outwards from B (towards A) ?s in
Figure 6.4, the constraint force acting on A can be written as fNi?i where fNi > 0
An equal and opposite force fNi1?i acts on B. The magnitude of the ith constraint
force is thus fNi. The n-vector collection of constraint force magnitudes is denoted
/?; the ith element of !N is fAri. (In general, matrices and vectors will be denoted
by boldface type; for example, "!N" Components of matrices and vectors are
97
98
written in regular type; for example, "fNi".)
If a constraint function Cj(?(t)) is formulated for the ith contact point, then
C?i(T(to)) is a measure of the relative normal acceleration at the contact point, as
described in Chapter 6. For simplicity, we will denote this normal acceleration
c??(r?(t0)) by the quantity ?Ni. The n-vector collection of normal accelerations is
denoted a?. The vector aN can be considered a function of the unknown constraint
force magnitudes !N. Since acceleration is a linear function of force, a? can be
expressed, using the derivations of Chapter 6, as
aN = A!N --H b			(8-1)
for some n x n matrix A and fl-vector b.
The matrix A reflects both the inverse inertias of the contact bodies, and the
contact geometry. The vector b reflects relative normal accelerations at contact
points due to known forces; this inchides both inertial effects and acceleration due
to known external forces. The matrix A is well known to be positive semidefinite
(PSD) and symmetric. The vector b is always in the range of A, that is,
b = Ac			(8 2)
for some n-vector c; situations iii which this is not so represent unsatisfiable
kinematic constraints on the system and are not of concern to us.
Since inter-penetration is prevented by requiring the relative normal accelerations
at each contact point to be non-negative, the non-penetration constraint for the
entire system can be written as
aN?>D			(1<?i<n).			(8--H3)
NVe denote this as a? > 0 or using Equation (8--H1), as
A!N --H b> 0.
(8--H4)
99
Likewise, the constraint that each fNi > 0 is written as
>- 0.
(8--H5)
Recall from Section 2.5.2 that a constraint force, in addition to being non-
negative, also has to be zero when the constraint is being broken. As in Section 2.5.2,
this can be expressed as the constraint
(fNi)(aNi) = 0			(1 < i ? n).			(8--H6)
Since fN > 0 and aN = A!N --H b > 0 imply !?Ta? > 0, Equations (8-4)
through (8--H6) can be written as
aN = A!N --H b> 0			fN > 0			and			!?Ta? = 0.			(84)
Equation (8--H7) defines what is known as a linear cornplcmcntarnty problem
(LCP); a vector !N satisfying Equation (8--H7) is a sohition to the LCP. Thus,
the constraint force magnitudes for a frictionless configuration can be computed
by solving an LCP. NVe will call a set of normal forces satisfying Equation (8--H7) a
valid set of constraint forces.
LCP's are closely related to quadratic programming. A quadratic program (QP)
has the general form
minimize xTMx + cTx subject to
x
Mx > d
x>O
Any LCP can be written as a QP. Using
!NTa? = fNT(AfN --H b) = !NTA!N --H !NTb,
Equation (8--H7), phrased as a QP, is
A!N > b
>- 0
(8--H8)
(8--H9)
(8--H10)
minimize !NTA!N --H bT!N subject to
100
The objective function !NTA!N --H !NTb achieves zero for every !N that satisfies
Equation (8--H7), and no others. Additionally, the QP formulation can directly
represent equality constraints. If the ith constraint of the system is an equality
constraint, the condition fNi > 0 in Equation (8--H10) can be removed, and the
condition ?Ni > 0 can be changed to aN i = 0. The QP in this case is still of size
n. The same equality constraint can be indirectly represented in Equation (8--H7) by
a reformulation that adds the constraint ?Ni > 0 and replaces the variable fNi
with fNi = fN'i --H fN11i where fN'i > 0 and fN11i > 0. The wsulting LCP is of size
n + 1 and is positive semidefinite, but not positive definite.
In the general case, determining if an LCP has a solution or a QP can achieve
a certain minimum value are NP-complete problems [40,52]. The corresponding
problems of finding a solution to an LCP or finding a minimizer of a QP are NP-
hard. However, when the coefficient matuces of the QP or LCP are PSD, then the
problems have polynomial time solutions[40]. LCP's and QP's with PSD coefficient
matrices are known as convex problems. In practice, convex QP's and LCP's are
solved by efficient algorithms whose worst case behavior is exponential but whose
expected running time is polynomial. NVe can use sucli an algoritlim to compute fN
and determine the constraint forces for configurations without friction.
8.1 Existence and uniqueness of solutions
There are two basic Q?estions about !N that should be asked. First, does an fN
satisfying Equation (8--H7) always exist? Furthermore, if an !N exists, is it unique?
(Equivalently: do valid constraint forces exist for every frictionless configuration?
Are they unique?) For the case when the LCP of Equation (8--H7) is convex, and
b lies in the column-space of A, an !N satisfying the LCP always exists. If A,
in addition to being PSD, is non-singular, then A is positive definite (PD), and
101
!N is unique. Otherwise, although !N is not necessarily unique, aN = AfN --H b
is unique; Cottle[7j gives an elegant proof of this fact. The non-uniqueness of fN
arises in cases where the distribution of constraint forces among contact points is
underdetermined. The canonical example is of a chair resting on the floor on four
or more legs. Although the acceleration of the chair is unique, the distribution of
force on the legs is not.
8.2 Implementation issues
For frictionless configurations, the problem of determining a valid set of constraint
forces is well in hand. ?7e will conclude this chapter by discussing two implementa-
tion issues in computing !N.
The first issue involves avoiding computational effort by reusing information fiom
previous time steps. It is well known in the mathematical programming literature
that when Equation (8--H7) has a solution, it has at least one solution called a basic
solution. A basic sollition may be quantified by a subset L C ?1, 2,. .. , as follows.
Let B denote the matrix obtained from A by deleting each row and column whose
index is not in L; this yields a square sub-matrix of A. Similarly, let d denote the
vector obtained from b by deleting elements of b whose index is not in L
The special properties of a basic solution are as follows. First, the matrix B
is invertible. Second, the vector B--H1d satisfies B--H1d > 0. Last, let us expand
the vector B--H1d by inserting zero elements. The zeroes are inserted into B--H1d by
inverting the process of deleting elements from b to form d; the expanded vector
has a zero for every index i not contained in L. The newly expanded vector is then
a solution to Equation (8--H7).
The concept of a basic solution also applies to the QP of Equation (8--H10); a
solution to Equation (8--H10) can be obtained by inverting a matrix formed from the
102
coefficient matrices defining the QP. Most algorithms for solving LCP's and QP's
find a basic solution. In a simulation enviwument, the concept of a basic solution
can greatly reduce the amount of numerical work done at each time step. The major
effort of computing the constraint forces lies in solving for !N. However, given a
basic solution for !N at one time step, the simulator can usually exploit coherence
in L at the next time step. Although the coefficients given by A and b vary from
one time step to the next, the index set L denoting a basic solution usually remains
valid for some large number of time steps. NVhenever L is known, !N is easily found
by computing B--H1d. As a result, !N can be computed by solving only a linear
equation for most time steps.
Using L from the previous time step, we may occasionally produce a singular
B for the current time step or a non-positive vector B--H'd, indicating that L is
invalid for the current time step. This usually happens when conditions change
abruptly due to an impact. It also happens whenever an ?Ni becomes positive,
indicating separation at a contact point. Since this may indicate a discontinuity
in acceleration [30', the ordinary differential equation integrator is restarted at this
point. Whenever L becomes invalid, a solution for !N must be computed from
scratd?. The extra cost of forming B and d and computing B--H'd even when L
is invalid is still much less than tlie cost of always coniputing fN froni scratch.
Gilmore and Cipra[20] discuss methods for predictiiig the change in L from one
time step to the next.
The second implementation issue has to do with the forniulation of !N in terms
of a QP. In analysis, it is well known that a necessary, but not sufficient condition for
a function y(x) to achieve a minimum at xo is that the gradient of y vanish at ro
This is known as a first order optimality condition. For a QP, analogous first order
conditions called the Karush-Kuhn T?ickcr (NNT) conditions for optimality exist.
In the case of a convex QP, the NKT conditions are both necessary and sufficient;
103
the global minimum is achieved if and only if the NKT conditions hold[40]. For a
PSD symmetric matrix A, as is the case for frictionless configurations, the QP
minimize 2l!NTA!N --H !NTb subject to !N > 0			(8-11)
has KNT conditions
A!N--Hb=A, !NTA>0 and A>O.
(8--H12)
(The elements of the fl-vector A are referred to as the Lagrange niultipliers of
the minimization problem.) A solution of Equation (8--H11) is therefore a solution of
Equation (8--H10). The vector !N may be computed by solving either Equation (8--H10)
or Equation (8--H11). The form of the two different QP formulations will become
more important when we deal with static friction iii Chapter 14. In actual practice,
Equation (8--H10) is computationally preferable in dealing with frictionless systems.
In dealing with the less-than-perfect world of numerical solution and optimization
including the explicit constraint A!N --H b > 0 in the problem formulation ii-icreases
the likelihood of finding an !N that satisfies the most important constraint of non-
penetration.
Chapter 9
The Coulomb Eriction Model
At this point, we will extend our treatment of contact between bodies to the case
when frictional forces act at contact points. ?Ve will call the force that arises between
two bodies at a contact point a contact force. A contact force is the sum of a
workiess constraint force, and a friction force. The constraint force does no work on
the system and acts normal to the contact surface; tlie constraint force is sometimes
called the normal force. The friction force acts tangential to the contact surface
and may perform work by dissipating kinetic energy between contacting bodies. In
this chapter, the model of friction used thro'igliout the rest of this thesis is defined.
Following this chapter, we will consider the theoretical and practical complexity
of computing both normal and friction forces. NN'e will say that the contact forces
acting on a configuration with friction are valid if the friction forces satisfy the
conditions defined in this section, and the constraint forces satisfy the conditions of
the previous chapter.
1()5
106
9.1 Coordinate geometry at a contact point
At each contact point, let us define a local coordinate system as follows. Let n,
l? and t? denote mutually perpendicular unit vectors. The vector ? is normal to
the contact surface at the contact point while t4 and t% span the plane tangent to
the contact surface. As in Chapter 8, we adopt the convention that ?? points from
B towards A; for example, n would be VG(pb,P)/IIVG(pb,p?)II in Figure 6.4. ?Ve
have already defined fN and aN, as the contact force and relative acceleration in
the n direction. (In situations where we are discussing only a single contact point,
we will drop subscripts and write fN and a? in place of fNi and aNi.) Similarly,
we decompose the friction force along the t% and ?y axes into magnitudes fr and
f?. Following the same convention as for normal forces, the friction force acting on
A is fx?x t fyt%, with an equal and opposite friction force acting on B
Since friction acts to oppose any slipping motion, let us consider the slip velocity
tangent to the contact surface. In Equation (6--H32) we derived s, the slip velocity
between A and B at the contact point. Let s be projected onto the tangent plane
and then decomposed along the tx and ty axes into maguitudes Vx and v?:
Vx			= tx S			and			Vy = ty S.			(9--H1)
Thus, the tangential slip velocity (henceforth called the tangential velocity) is
Vxtx + Vyty.
?7hen the slip velocity s is zero, let the tangential slip acceleration be described by
magnitudes a? and ay where
=			ft8			and			ay = ty			d			(9--H2)
?dt8
If the tangential velocity is initially zero, then as long as the tangential acceleration
107
a?t? + a?t? remains zero (and contact is not broken), -the tangential velocity also
remains zero.
The Coulomb model of friction is an accepted empirical relationship between the
normal force magnitude fN and the friction force magnitude j? + j?. At all times,
the condition
+			< (I'fN)2
holds, where ? is a coefficient of friction that depends on material properties, and
may be different at each contact point. If the tangential velocity is non-zero, we say
that the friction force is dynamic; otherwise, the friction force is called static.
Typically, the coefficient it of static friction (jtstatic) is larger than the coefficient
? of dynamic friction (Jtdynamic). Aside from this, it is roughly independent of the
tangential velocity. Since we are computing contact forces for a specific instant of
time to, we know which coefficient of friction to use, and we will not bother to make
a distinction between fl-static and fl-dynamic (?Ve also will not bother to index fi over
different contact points.)
9.2 Dynamic friction conditions
?N-7hen dynamic friction arises, the friction force magnitude achieves its upper bound
of flfN Also, the friction force acts directly opposite the slip velocity. Nuowing
that ff + fy2 = (flfN)2 and that the friction force is opposite the slip,
or
Vxtx+vyt%
fxfx+fyt%=--HflfN
Hvy2
(9--H4)
Vx
fx = !ifN??y2			and			f? = --HilLyVy			(9-5)
?Vy2
Thus, given ii and the tangential velocity, the dynamic friction force can be written
108
solely in terms of fN - Since
fNfl+f??x+f?,t% = fNn-
?!ifNVxtx--HIIfNVyt%
Thvy2
=			IIVxt?--HIiVY?Y),			(9--H6)
Thvy2
the direction of the net contact force is predetermined given the tangential velocity
and ji.
9.3 Static friction conditions
The conditions between fx, fy and fN are more complex when the tangential
velocity is zero, resulting in static friction. There are two possibilities, according to
whether or not the tangential acceleration is zero. First, f? and f? may satisfy
=			= 0			and			j?2 +			< (ItfN)2			(9--H7)
indicating that the friction remains static. Second, fx and fy may satisfy
ff +			= (jifN)2			awl			(fx1x + fyJy) (a?i.? + ay?.y) ? 0,			(9--H8)
indicating that slipping is occurring, and that the static friction force has become
dynamic. In this case, the friction force must oppose the initial tangential acceler-
ation and have magnitude !ifN- The friction force direction does not have to be
directly opposite the tangential acceleration direction.
9.4 Impulsive frictional forces
When two bodies collide at a contact point with friction, we can consider the collision
as taking place over some small but non-zero time interval. We then consider the
109
limiting value of the product of the force acting between the bodies and the time
interval, as we let the time interval tend to zero. This limiting value is the impulsive
force of the collision. During the time interval of the collision, the normal and friction
forces must satisfy the conditions of the previous two sections. It is possible that
during this time interval, the direction of the relative tangential velocity between
the bodies may change. For planar collisions, this happens if the relative tangential
velocity passes through zero and reverses itself. ?Vang and Nlason[34,54] describe
the computation of the impulsive forces for all possible planar collisions involving
one contact point. For three-dimensional collisions, the velocity direction can vary
in the tangent plane of the collision, and the computational problem becomes much
more difficult.
For simulation purposes, the simulator described in this thesis makes the sim-
plification that the direction of the relative tangential velocity does not change
during a collision. For planar collisions, this amounts to ignoring the possibility of
reversed sliding, as considered by NVang and Niason. For threc-dunensional collisions,
this is simply an out-and-out approximation that is reasonable for some situations,
and probably unreasonable for others. This model is no more than an ad hoc
approximation that allows the sim'ilator to produce a (hopefully) reasonable looking
collision response behavior. However, the model, though incorrect physically, is
computationally attractive because it produces precisely the same sort of constraint
equations as are encountered for static friction.
The idea is the following. Given a collision at a contact point, let fN denote
the normal impulse magnitude, fx and f? the tangential frictional impulse, and
vN and vk the relative approach speeds in the the normal direction prior to the
collision, and after the collision, respectively. Similarly, let Vx+ and Vy+ denote the
tangential velocity components after the collision. The constraints on the normal
quantities are fN > 0, and ??N+ = --HevN where ? is the coefficient of restitution of
110
the collision. Note that this kinematic condition, which is Newton's hypothesis, is
known to violate energy conservation under some circumstances if both E and ? are
non-zero[3]. The standard condition j?? + < (JtfN)2 is imposed at all times on
the frictional impulse.
If the bodies are sliding after application of the impulse, then the frictional
impulse must achieve its upper bound, and must oppose the sliding velocity. This
is written as
+			= (iifN)2			and			(fxl% + fyt%) (vz+4 + Vy+t?) ? 0			(9--H9)
Otherwise, the bodies are not sliding, and the frictional impulse need only satisfy
its upper bound condition. This yields
--H			--H 0			and			fx2 +			? ([tfN)2			(9--H10)
Comparing these last two equations to Equations (9--H8) and (9--H7), it is clear that
both sets of equations can be solved by the same techniques. Accordingly, we will
be able to use the solution methods described in Chapter 14 to compute frictional
impulses.
If we are dealing witIi a number of contact points simultaneously, we will need
to add the following kinematic constraint. If at some colliding contact ?k > --H??`N
then the normal impulse at that contact n?ist be zero. In other words, if bodies
bounce away at a contact point with a speed in excess of that predicted by
the coefficient of restitution, the normal impulse at that point must be zero[13].
Satisfying these conditions is still exactly analogous to satisfying the static friction
conditions for a collection of non-colliding contact points, so again, we can compute
frictional impulses using the static friction methods described in Chapter 14.
Chapter 10
Dynamic Friction
Having defined a model of friction we now consider the complexity of computing
valid contact forces in the presence of friction. Most of the complexity results in
this thesis apply to configurations with dynamic friction only. Complexity results
for static friction are presented in Chapter 14. Until then, only configurations of
bodies with dynamic friction are considered.
Given that no static friction forces act at any contact points, the contact
force direction at each contact point is given by Equation (9 (3), since ii and the
tangential velocity are known quantities. Thus, as in Chapter 8, we are interested
in computing only the vector !N of constraint force `nagnitudes. Solutions for
!N for configurations with just dynamic friction can still be expressed as the LCP
of Equation (8--H7). However, as a result of the contact force having a tangential
component, the matrix A in Equation (8--H7) is no longer necessarily symmetric
PSD and Equation (8--H7) is no longer a conwx LCP.
As a result, several nice properties are lost. Although solutions to convex LCP's
can be found in polynomial time, non-convex LCP's do not share this property;
finding solutions for non-convex LCP's is NP-hard. What then, is the complexity
111
112
of solving the non-convex LCP's corresponding to configurations of bodies with
dynamic friction? Also, the uniqueness and existence properties associated with
the convex LCP's of Chapter 8 no longer always hold. With dynamic friction, it is
possible that there may be no solution !N, and thus no valid contact forces. We call
a configuration exhibiting this behavior inconsistent.1 It is also possible that the
acceleration produced by solutions of the LCP is not unique; such a configuration
is called indeterminate.
Configurations with one contact point that exhibit inconsistency and indetermi-
nacy have been discussed by Erdmann[12], L5tstedt[28], and ?Iason and Wang[35].
L5tstedt[30], realizing that indeterminacy and inconsistency major present difficul-
ties for a simulation process, proposes a modification of the Coulomb friction law
that eliminates both inconsistency and indeterminacy, and fi?rther causes the LCP
to always be convex. We will discuss Lo?tstedt's modification in Chapter 14, when we
consider static friction. In this chapter, we discuss two well known one contact point
configurations that exhibit inconsistency and indeterminacy. Chapter 11 uses these
configurations as building blocks in a proof about the computational complexity of
solving dynamic friction problems.
10.1 Indeterminacy
In Figure 10.1, body A is a thin rod of length two with a symmetric mass distribution
that contacts body B at a single contact point. Body B, the "base", is fixed. The
1 The use of the word "inconsistent" throughout this thesis requires some qualification. The fact
that configurations exist which have no valid contact forces is counter-intuitive to most people at
first. These configurations are Inconsistent" with the initial expectation that a valid set of (finite)
contact forces must exist for any configuration that can be physically realized. However, we do not
use "inconsistent" to mean "unnatural", nor do we mean to imply that inconsistent configurations
are in any way extraordinary physical occurrences. Readers who are familiar with "self-locking"
behavior of assemblies due to friction will find "inconsistent" configurations to be nothing special
at all.
113
example is two-dimensional, so the tangent direction is written simply as t. The
particular values of I, 0 and ji given in Figure 10.1 are chosen to simplify later
computations. A gravity force of --Hmgn acts on A. Let the angular velocity w and
gravity constant g be such that A is rotating clockwise and
jw2sinO--Hg= 1.			(10--H1)
Given these values, let the linear velocity of A be chosen so that the point Pa of
the rod is sliding tangent to B as shown. Since A and B are neither colliding nor
separating at Pa, a contact force is presumed to act. Since the tangential velocity
is non-zero, a dynamic friction force acts on A in the t direction, From Chapter 9,
the net contact force acting on A is
+ I'fNJ fN(? + [`?)
Because 1? is fixed and A is thin, the normal acceleration aN is simply
(10 2)
aN =			Pa.			(10--H3)
To compute aN, let a and ? denote the linear and angular acceleration of A,
w,the angular velocity of A, and r, the displacement of Pa from the center of mass
of A. Altho?igh this configuration is two-dimensional, we will regard a, Pa Pcz, fl
and r as 3-space vectors lying in the xy plane, and w and ? as 3-space vectors
parallel to the z axis. From Figure 10.1, r = cos 0, --H sin 0, 0). The acceleration
Pa may be expressed as the sum of three terms: the linear acceleration a, the
tangential acceleration ? x ?, and the centripetal acceleration w x (w x r). The
linear acceleration, a, is
a =
fNn+[ifNf +(--Hrng)?
ni
fN?+IifNt --H gn%.
rn
(10--H4)
114
Variables
contact point			Pa			contact point velocity
unit surface normal			?			unit surface tangent
A's mass			I			A's moment ofinertia
A's angular velocity			?			coefficient of friction
gravitational acceleration
Pa
n
m
0)
g
Relations
I = ?			0 = 720			16(cos2 0--H ticosOsinO) = --H2			(? = 3/4)
16
n
Pa
-mgn
(gravity)
B			Pa
Figure 10.1: A one contact point configuration between a thin rod A and a fixed
base B.
115
The torque on A about its center of mass is
r x (fNit + IifNft
which yields an angular acceleration of
rx(fNit+ [tfNt?
Then
(10--H5)
Pa = ? + ? x r +w x (w x r)			(10-6)
fNit + IifNt --H gn + TX(fNit + JifN1) X r + w X (w x
m			I
Taking the dot product of Equation (10--H4) with it,
it			= fNitit+[IfNitt --H J(it			_ fN --H g.			(10-7)
T?1
Taking the dot product of the tangential acceleration ? x r with it and using the
geometry of Figure 10.1,
and
n w X (w x 7) = Iwi2sino
it (? x r) = it			rX(fN)+I'fN?)
x r
fN(cos2O--Hjtcos0sin0)
Then, from the relations in Fig?ire 10.1,
it Pa --H L?q+ f?(cos2O--HiicosOsinO) + Iwi2sinO
rn			I
--H fN
--H7n(1+16(cos2O--H;cosOsinO))+?w2sh?O--Hg
fTy(1?2)+IwI2shlO?g
fN
--H			7?+(KI2sin0?g).
(10--H8)
(10-9)
(10--H10)
116
Finally, from Equation (10--H1),
Pa =
fN +(w?2sinO--Hg)
m
fN + 1.
m
(10--H11)
From Equation (8--H7), the constraints on fN are the LCP
fN > 0			+ 1 > 0			and			fN(?L + 1) = 0.			(10-12)
fl?			fl?
This LCP yields two solutions for f?; one solution is fN = 0, and the other is
fN = m - For the fN = 0 solution, ?? Pa = 1. In this solution, the centripetal
acceleration of Pa is stronger than the force of gravity pulling A down; thus, A
merely continues its rotation and the point Pa moves off of B (Figure 10.2a).
In the second solution, Lv = 7?? and n Pa = 0. A normal force of fl??
and a friction force of [tn?t act on A. The torque generated by friction balances
the centripetal acceleration of pa; as a result, A and B do not break contact
(Figure 10.2b). Note that the only valid values of f]\T for this configuration are
fN = 0 or fN = ni - Since the solutions produce different accelerations for A, the
configuration is indeterminate.
The fact that there is more than one possible behavior for A should not be
of great concern. ?Ve have already seen that the rigid body assumption allows
for non-unique distributions of contact forces, as in the case of a chair resting on
four legs. When the rigid body ?ssumption is relaxed to allow deformation, tlie
indeterminacy in the force distribution is resolved. In the c?se of the sliding rod,
sufficient knowledge about the flexibility and compressibility of the rod would serve
to determine the actual behavior of the rod. Since we have chosen to treat bodies
as perfectly rigid, we are in essence stating that this knowledge is not available to
us.
117
(a)f?=O			(b)j;v=m
0)
M
A			A
B
Pa
0)
8			?mt
mn
Fignre 10,2: An indetern?inate confignration (a) The contact force bet?veen 4 and
B is zero The point Pa rotates to the left and np, breaking contact with I?. (b) The
normal and friction forces balance gravity and centripetal acceleration; Pa nioVes
horizontally and maintains contact with B.
118
(a) no solution for fN			(b)fN =0
v
A			A
--Hmgn
B			?kfNt
fvn
v
B
Pa
Figure 10.3: A potentially inconsistent configuration. (a) For any fN > 0 the
point Pa is accelerated dowuwards into B. This configuration is inconsistent. (b)
The configuration has a unique solution of fN = 0 when gravity is removed; A
skims along the surface of B. This configuration is not inconsistent.
10.2 Inconsistency
To demonstrate inconsistency, consider the same sliding rod configuration, but
with a different angular velocity. Suppose that A has an angillar velocity of zero
(1w = 0), and a linear velocity v opposite [, as shown iii Figure 10.3. Then the
condition ?N =			P?a > 0 is
o+ j;a =			fN --H g >0.			(10-13)
However, if g > 0 (Figure 10.3a), then fN > 0 implies that ?? Pa < 0 Thus no
(positive) value of fN is sufficient to prevent Pa from accelerating downwards and
into B. This configuration is inconsistent.
Note that the value of g here is crucial. If g = 0, so that no external force acts
119
on A, then fN = 0 becomes the (unique) valid contact force (Figure 10.3b). Any
positive value of fN for this configuration results in inter-penetration, Figure 10.3b
corresponds to to Pa "skimming" horizontally over B, with neither a normal force
nor a friction force exerted on A. If g becomes even slightly positive however, the
configuration is inconsistent.
How can we interpret this phenomenon? Although our mathematical model
admits no solution for this configuration, we know that there must be some physical
answer. Still, how can we get around the fact that no matter what value we assume
for fN, the point Pa accelerates dowuwards? The resolution of the issue lies in
how we view the subject of inconsistency. ?Vhen two rigid bodies collide, no force,
no matter how large, is sufficient to prevent some degree of inter-penetration. For
rigid body mechanics, the need to prevent inter-penetration requires the extension
of the rigid body model to include the concept of impulsive forces which can
discontinuously change velocity to prevent inter-penetration. Impulsive forces are
added to the rigid body model to avoid the inconsistency that non-imp'ilsive fi?rces
are not always sufficient to prevent inter-penetration.
For the sliding inconsistent rod of Figure 10.3a, non-impulsive forces cannot
prevent inter-penetration because the constraints of dynamic friction result in a
fixed direction for the contact force, no matter how strong. Ii-i contrast, a contact
impulse applied to A and B can prevent inter-penetrat ion because the direction
of the contact impulse is not fixed. After the impulse is applied, the inconsistency
is removed and a non-impulsive contact force can be found. Mason and NVang[35?
discuss the computation of this frictional impulse and show that it is consistent
with the way in which frictional impulses are determined for collisions between rigid
bodies.
Interestingly, what appears to disturb people most about the configuration of
Figure 10.3a is that an impulse is necessary even though the relative normal velocity
120
between A and B is zero. Because of this, the configuration is seen as pathological,
and the application of an impulse is seen as an action of last resort. In the next
chapter, we will consider the computational complexity of finding non-impulsive
solutions, so as to avoid applying impulsive solutions as much as possible.
Chapter 11
An NP-complete Class of
Configurations
?Ve define the frictional consistency problem as the problem of deciding if a given
configuration with dynamic friction is consistent. In this section, we prove that
the frictional consistency problem is NP-complcte. ?Ve first show that the frictional
consistency problem lies in NP and then establish the frictional consistency problem
as NP-hard. Although the configurations constructed in this section seem contrived
(and arguably are), the NP-hardness result has grave implications even when
inconsistency is not a concern during simulation.
DEFINffIoN. An instance of the frictional consistency problem is a configuration
C of bodies lliat contact at n distinct contact points. The physical properties of each
body (mass, moment of inertia, linear and angular velocity, position and orientation,
and external forces) are described by rational numbers. The specifics of a contact
point (position, coefficient of friction, surface normal) are also described by rational
numbers. The relative motion between bodies at contact points with friction is non-
zero in the direction tangent to the contact surface and zero in the direction normal to
121
122
the contact surface. The notation CI = k means that configuration C is describable
in k bits. Clearly n < k.
THEOREM 1. The frictional consistency problem lies in NP
PROOF. Given an instance of C, an LCP of size n with the following two properties
exists.
(1) If C is consistent, then an n-vector !N that is a solution to the
LCP exists. The set of contact forces such that the magnitude of
the normal force at the ith contact point is fNi is a valid set of
contact forces for the configuration C.
(2) If C is inconsistent, the LCP has no solution.
The LCP is the LCP of Equation (8?7). The numencal quantities in the LCP are
computed from the rational entries of C in a total of 0(n3) arithmetical operations.
The LCP can therefore l)e constructed in time polynomial to k. The problem of
deciding if an LCP has a solution lies in NP[4o]. It follows from this that deciding
frictional consistency is also in NP. 5
In order to show that deciding frictional consistency is NP-hard, we will reduce
the NP-complete problem "subset sum" to the flictional consistency problem.
DEFINITION. An instance of the subset sum problem is a pair (A, S) where A =
....... , ?? is a set of positive integers and S is a single positive integen A subset
sum instance (A, S) is satisfiable if there exists a subset ?? C A such that
a = S.			(11?1)
aEA'
Deciding if an instance of the subset sum problem is satisfiable is an NP-complete
problem[iG].
123
To show that deciding frictional consistency is NP-hard we will take an arbitrary
instance (A, S) of the subset sum problem and construct (in polynomial time) a
configuration of bodies C. The configuration C will have the property that C is
consistent if and only if (A, 5) is satisfiable.
THEOREM 2. Deciding frictional consistency is NP-hard.
PROOF. We begin by presenting a slightly modified version of the configuration of
Figure 10.3. In Figure 11.1, body B is free to move vertically, either up or down,
but it is still not allowed to rotate or move horizontally. Let the mass Al, of body
1?, be Al = lOOm so that B is still massive compared to A. Assume that there is
no external gravity force acting on A, but that B is subject to a vertical external
force fb??, though its center of mass (Figure 11.1). A positive value for fb indicates
the force acts upwards, pushing B towards A
Since B may now move vertically, Equation (10--H13) is modified. Since there is
no external force on A, g = 0. In place of --HLv/m --H g > 0,
_			+ L --H kb > ?			(11-2)
m			Al			Al --H
The term fN/Al is the linear acceleration of B away from A due to the contact
force. The teH? fb/Al is the upwauls linear acceleration of B due to the external
force. Since Al = lOOm,
--H fN +L?Lb _			99Lv _ fb
m			Al			Al			100 m			Al			(11--H3)
Equation (10--H13) is then rewntten as
99fN _ Lb > 0			(11--H4)
100 m			Al --H
As before, if fb > 0, an inconsistency occurs. If fb = 0, then fN = 0 is the unique
valid solution and the configuration is consistent.
124
v
j?n
B
fNfl			fNn
Figure 11.1: A possibly inconsistent configuration. Body B n?ay n?ove vertically,
but is not allowed to move horizontally or rotate. If the external vertical force
magnitude fb is positive, the configuration is inconsistent. The contact force fN?
on A acts equally and opposite on
125
E1
v
A			B			A
Figure 11.2: A configuration whose consistency depends on an external force.
Body B is constrained by the fixed wedges and can only move horizontally. The
configuration is consistent only if B is not subject to a net horizontal force.
Now, consider the configuration of Figure 11.2. Body B of Figure 11.2 is
initially at rest and is positioned by four fixed triangular wedges that contact B
without friction. Body B is therefore free to move horizontally, but can neither
rotate nor move vertically. On either side of body B are thin rods E1 and 2
that have no angular velocity, and a linear velocity as indicated. The rods i
and 2 contact B in the same manner ?s the configuration of Figure 11.1 except
that the frames of reference for i and 2 are rotated by 900 with respect to
Figure 11.1. In the configuration of Figure 11.1, inconsistency occurred if fb was
positive, pushing B towards A. The same holds true for Figure 11.2. If B has
an acceleration left??rds (towards i), then inconsistency occurs. Likewise, if B
has an acceleration rightwards (towards 2), then inconsistency also occurs. Thus,
the configuration of Figure 11.2 is consistent only if the net horizontal acceleration
of B is zero. In this case, the rods i and 2 skim along the surface of I? as in
Figure 10.3b.
Now consider Figure 11.3, where a collection of thin rods ?.,... , R? have been
added. In addition, an external horizontal force with magnitude jiS acts on I?,
126
0)
0)			0)
E1
??ifi			??if2
v
B
E2
Figure 11.3: A configuration that is consistent if and only if the friction forces on
B sum to ?tS.
trying to accelerate I? to the right. Each rod Rj has mass fl?j. The configuration
between each rod Rj and B is the same as the configuration of Figure 10.2; thus
each rod Rj has angular velocity w and is subject to an external gravity force. Let
fNi be the magnitude of the non?al force between Rj and B. As in Figure 10.2,
the only valid solutions for fNi are Lvi = 0 and fNi = m?. If fNi = 0, then no
friction force acts between Rj and B. Otherwise, fNi = Thj and a friction force of
magnitude jtfl?j acts between Rj and B. The friction force pushes l?? to the right
and B to the left, with magnitude jilni. The friction force on B therefore acts to
oppose the external force of magnitude 1iS.
In order for the configuration of Figure 11.3 to be consistent, R must have no
net horizontal acceleration. This means that the friction forces exerted on B from
the n rods must sum to 1jS, balancing the external force applied to B. Thus, the
configuration is consistent if and only if
n
? [ifNi = [IS.
i=1
(11--H5)
Since each fNi is either 0 or fl?j, the configuration is consistent if and only if some
subset of fmi, . .. , sums to 5.
127
We can now perform the reduction from subset sum to show NP?hardness. Given
any set A = ....... , ani and any target sum S, construct the configuration of
Figure 11.3. The values chosen in the previous section for jt, 0 and w may introduce
irrational quantities into the descriptions of the configurations. Since the frictional
consistency problem must be described in terms of rational quantities, we may
assume that ii, 0 and w are suitably perturbed so that all quantities are rational
without otherwise disturbing the sollition properties of the configurations. This can
be done since the set of rational numbers is dense with respect to the set of real
numbers. In the constructed configuration, assign m? = a? for 1 < i < n and
let an external horizontal force of jiS act on as shown in Figure 11.3. By the
above discussion, the configuration is consistent if and only if there exists a subset
of frni, . . . , ni?1 that sums to 5. But since A = fm?1, . . , , the configuration
is consistent if and only if (A, 5) is satisfiable. N'\'e conclude that the problem of
deciding frictional consistency is NP-hard. [Zi
THEOREM 3. Deciding frictional consistency is NP-coniplete.
PROOF. The result follows immediately from Theorem 1 and Theorem 2. E
COROLLARY 1. Computing contact forces (if U?ey exist) for a configuration is NP-
hard.
PROOF. Since deciding if a set of contact forces exists is an NP-complete problem,
computing the contact forces (if they exist) is an NP-hard problem. c
A reasonable response to this corollary is to question its practical applicability
to the problem of computing valid contact forces. Suppose the possibility of an
inconsistent configuration occurring in a simulation is so unlikely that it can be dis-
counted completely. Can we construct an efficient algorithm that computes contact
128
forces only for the set of consistent configurations? The following corollary asserts
that a polynomial time algorithm for computing contact forces, even restricted to
consistent configurations, does not exist unless P = NP.
COROLLARY 2. A polynomial time algorithm for computing valid contact forces for
consistent configurations exists if and only if P = NP.
PROOF, Suppose that P = NP. Since deciding if an LCP has a solution lies in NP,
P = NP implies a polynomial time algorithm for finding the solution to an LCP.
Since valid contact forces for a consistent configuration of bodies can be found by
solving an associated LCP, valid contact forces are computable in polynomial time
if p = NP
Conversely, suppose that contact forces for consistent configurations can be
computed in polynomial time. Then there exists a machine Al and a polynomial
p with the following behavior. ?Vhenever Al is given a consistent configuration c
as input, Al outputs a valid set of contact forces within time p(ICI). Al's behavior
when C is inconsistent is undefined. Given any configuration C, not necessarily
consistent, Al can be used to decide consistency in polynomial time as follows.
Let C be input to Al and run for p(ICI) time. If Al fails to output within
this time, then C is inconsistent. Otherwise, Al has produced some output. Since
deciding frictional consistency is in NP, the validity of Al's output can be decided
in an additional amount of time that is also a polynomial fiinction of C?. If Al's
output is a valid set of contact forces, then clearly C is coiisistent. If Al's output
is invalid, then C must be inconsistent (else Al would have output a valid answer).
In any event, the consistency of C has been decided in polynomial time.
Since deciding consistency is NP-complete, we conclude that the existence of a
polynomial time algorithm for computing contact forces on consistent configurations
129
would imply that P = NP. 5
From a simulation viewpoint, this seems to be a depressing result. However,
there has been an implicit, unchallenged assumption in the above discussion. The
assumption is that we must determine whether or not a configuration is consistent, so
that we will know whether impulsive or non-impulsive forces should be applied. This
involves a further assumption: it is invalid to apply impulsive forces to configurations
that are consistent (and thus have a non-impulsive force solution). In the next
chapter, the argument will be made that this second assumption is false. Based on
this argument, we will reformulate the problem of computing contact forces so that
it is neither necessary nor even desirable to determine whether or not a configuration
is consistent. Following that, it will be shown how avoiding the need to determine
consistency leads to a practical polynomial time algorithm for computing contact
forces for dynamic friction.
Chapter 12
Physical Models
In this section, a physical model for both inconsistency and indeterminacy is
presented. ?7e will use this model to develop a computationally practical method
of dealing with dynamic friction in Chapter 13. The model in this section was
developed in order to understand the ??ehavior of inconsistent and indeterminate
configurations. The insights gained from considering this model also suggest a
general principle regarding indeterminacy during simulation.
The motivation of a physical model stems from the need to answer the following
basic questions. First, what should be the result of a simulation when inconsistency
is encountered? As ?? have seQn, the only resohition of inconsistent configurations,
such as Figure 1O.3a, is the application of an impulsive contact force. NVhat is the
physical process which produces such a force? Second, is there a simple physical
explanation for the type of indeterminacy that tlie configuration of Figure 10.2
exhibits?
In the physical world, there is of course no such thing as a perfectly rigid body.
For near rigid bodies, contact forces arise as a result of small elastic deformations in
the neighborhood of the contact area. The penalty method, as described in Chap-
131
132
ter 7, is based on a mathematical method from numerical optimization. However,
the penalty method can equally well be viewed (and is, by many) as an attempt
to model elastic deformations in a very simple form. In the penalty method, small
amounts of inter-penetration model the amount of deformation between bodies, and
the penalty force accompanying the inter-penetration models the elastic restoring
force. Typically, the normal force is modeled as a linear spring force --HAd, where A
is the spring constant and d is the amount of inter-penetration. As Chapter 7 shows,
the penalty method converges to a correct result as I? is increased. It therefore
makes sense to consider applying the penalty model, conceptually, to situations
involving friction, and studying the limiting behavior. This thesis does not propose
the penalty model as embodying a specific physical process. Rather, the penalty
model is used as an intuition to help visualize why impulses are necessary in some
situations, and to help guide our development of a definition of exactly when and
how impulses can be applied to a system without colliding contact.
12.1 A model of inconsistency
Figure 12.1 shows the behavior of the inconsistent configuration of Figure 10.3a
when the penalty method is applied. At time to, consider the tip of the rod, Pa
to be resting exactly on B, with zero inter-penetration. Since there is no inter-
penetration, the normal force is zero. Even though Pa is sliding along B, the
friction force is zero since the normal force is zero. Since the only force acting on A
is the external gravity force --Hmgr?, Pa accelerates downwards. At time to t At, Pa
has inter-penetrated B by an amount d1, so a normal force A?d1i? acts on A. Since
Pa is still sliding, a friction force of 1tI?d1t also acts on A. The net result, from
Equation (10--H11), is that this causes Pa to accelerate dowuwards (even faster than
before). As the penalty force continues to increase, it causes more inter-penetration
133
Time t0:			Timet0+At:
zero inter-penetration			d1 inter-penetration
Pa
A
A
--Hmgn			--Hrngn
B
??Pa
B
U
Kd1(n+ptt)
??Pa
Fignre 12.1: A model of indeterminacy. At time t(), only gravity acts on A. At
time 1o + At, the inter-penetration distance is d1 and both a penalty and a gravity
force act on A, cansing ])(1'5 downwards acceleration to ????C??s?.
134
between A and B; a form of positive feedback. Accordingly, both the friction and
the normal force increase, and the cycle continues. Since we are trying to model A
and B as rigid bodies, the penalty constant K must be allowed to be arbitrarily
large. The larger K is, the faster inter-penetration increases and the faster the
normal and friction forces build.
Recall that the friction force opposes the sliding motion of A across B. By
making I? arbitrarily large, the friction force brings Pa to rest (horizontally) in an
arbitrarily short time. For example, suppose A is adjusted so that Pa comes to rest
within time At. Then the amount of inter-penetration is O(At 2), since the distance
traveled by Pa depends quadratically on the time for which it travels. In the limit
as A goes to infinity, the contact force on A acts as an impulse and instantaneously
brings Pa to rest horizontally, without inter-penetration occurring.
The impulse that brings Pa to rest horizontally produces a configuration in which
A and B are colliding (in the normal direction) at Pa This requires the application
of a second impulse, to counter the dowuwards velocity acquired by Pa. Once the
impulses have been applied and Pa is no longer sliding across 1?, dynamic friction is
replaced with static friction. As a result, the net contact force direction is no longer
fixed, and the inconsistency no longer exists.
12.2 A model of indeterminacy
Consider the indeterminate configuration of Figure 10.2, which has solutions fN = 0
and fN = rn. Using the penalty method, the indeterminacy can be removed by
assuming some amount of initial inter-penetration between A and B. If the initial
inter-penetration between A and B is zero (Figure 12.2) then no normal force exists
at time to and the net acceleration of Pa is upwards. Contact is immediately broken;
this behavior is shown in Figure 10.2a. However, if the initial inter-penetration
135
Time t0:
zero inter-penetration
Time t0 +At:
A			A
fl Pa
fl Pa
B
Pa
--Hrngn
`I
B
Fignre 12.2: A configuration with no initial inter-penetration. Only gravity acts on
A initially. The centripetal acceleration of A pulls Pa away fron? B and contact is
broken.
136
Time t0:			Timet0+At:
d1 inter-penetration			d2 = d1 inter-penetration
A
A
8
u
fl?Pa=O
8
Kd1(n+?)			Kd2(n+?)
-(7
fl?Pa=O
Figure 12.3: A configuration with an initial inter-penetration depth of d1 . Both
gravity and a penalty force act on A. Body A slides and falls without breaking
contact with B.
produces a normal force magnitude of ??i,then the normal and friction forces prevent
A from breaking contact with B. In Figure 12.3, let the initial inter-penetration d1
be Tn/I?. Then the normal force magnitude at time 1o is n?. Since Pa is sliding on
a friction force acts on A as shown. As A falls, maintaining contact with B,the
inter-penetration varies smoothly, produce a varying normal force. At time to + At,
body A still inter-penetrates B by an amount d1 d2, and the behavior of the
configuration is that of Figure 10.2b. Thus, the initial amount of inter-penetration
determines which behavior occurs.
The numerical methods we will use to compute forces and impulses do not
137
model inter-penetration. Instead of determining behavior by initial d?oice of inter-
penetration, we can consider the normal force that existed at a contact point at
a time t immediately before 1o as the initial conditions at time to. These initial
conditions are used to determine the behavior of bodies. This brings up the obvious
question of what initial conditions to use for the first time step of a simulation, or
after a discontinuous velocity change of bodies caused by a collision. The choice of
initial conditions depends upon the application; if there is no basis for preferring one
set of initial normal forces over another, the solution found by the numerical routines
solving the contact force equations arbitrarily determines the behavior simulated.
?Ve will consider, in a moment, the validity of arbitrarily choosing a path for the
simulation to follow. A much deeper question that we will not attempt to address
is whether or not following only this one path to the end of the simulation reveals
useful information.
12.3 Indeterminacy and inconsistency combined
Let us now consider a configuration that combines indeterminacy and inconsistency
so that it is indeterminate whether or not an impulsive sohition is required. Consider
Figure 12.4. Once again, the combination of gravity and the angular velocity of Ri
is the same as in Figure 11.2 and Figure 11.3. Similarly, a horizontal acceleration of
B results in inconsistency. If E1 is ignored for the moment, then both fN = 0 and
fN = m1 are valid solutions for fN. The only valid solution for the configuration
as a whole though, is fN = 0; a magnitude fN = rn1 pushes B to the left, causing
inconsistency.
If we subscribe to the physical model presented in this chapter, then it is
clear that whether or not inconsistency occurs depends solely on the initial inter-
penetration between R1 and I?. The implication of our physical model then
138
E1
M
v
(0
ff
??ifN
B
Figure 12.4: A configuration with a single impulse-free solution. The magnitude fN
must be either 0 or rni to be valid. However, fN = rn? causes inconsistency. The
only valid solution is fN = 0
is the following: even though the configuration is consistent, that is, it has an
impulse-free solution, there is no a priori reason to prefer this impulse-free behavior
to the non-impulse-free behavior for this configuration. Stated differently, if we
have no bias on the initial conditions of the configuration (that is, the initial
inter-penetration distances), then tlie inconsistency resulting from fN = ??i, and
subsequent application of impulsive contact forces must be considered as acceptable
a behavior as the application of non-impulsive contact forces resulting from fN = 0
Even though the configuration in Figure 12.4 has only one valid sollition of contact
forces, (fN = 0), it has two possible behaviors and is thus indeterminate.
12.4 Preferred solutions
By accepting the physical model based on the penalty method, we accept the fact
that the two behaviors of Figure 12.4 must be considered equally valid. Let us
contrast this model with a principle proposed by Nilmister and Reeve[25] called
the principle of constraints. The principle states that when computing forces for
139
a configuration of bodies, impulsive forces should be used only if non-impulsive
forces do not exist for the configuration. In essence, the principle of constraints
is an assertion that non-impulsive solutions are preferable to solutions involving
impulses. This principle seems (at first glance) sound; however, let us consider
argui?ents against the principle of constraints.
From the discussion of the previous section, we can view the NP-hardness result
of computing contact forces in a different light. By reasoning that solutions involving
impulsive contact forces are acceptable even when a non-impulsive solution exists,
we have extended the possible behaviors of configurations of bodies. (A formal
definition of validity with regard to contact impulses has not been given yet,
and will be presented in the next chapter.) If we subscribed to the principle of
constraints however, we would first attempt to solve the LCP of Equation (8--H7), in
order to find a solution of non-impulsive contact forces. This is a search problem.
Since some configurations (such as the configuration of Figure 11.3) may possess
exponentially many impulsive sollitions and only a single non-impulsive sollition, it is
not surprising that an attempt to search for this single sohition leads to NP-hardness.
Thus, the first argument we might advance against the principle of constraints is
that by accepting a bias against solutions with impulsive forces, we are required to
perform possibly an exponential amount of work. In contrast, disregarding this bias
and accepting both impulsive and non-impulsive solutions turns out to require far
less work (Chapter 13).
This is not really an argument, but more of a statement of the consequences
of preferring non-impulsive solutions to impulsive solutions. A better argument is
that the physical model proposed in this section, which is logically consistent with
frictionless behavior, and explains both indeterminacy and inconsistency, has no
such preference. In comparison, the principle of constraints is nothing more than
a statement of the preference of non-impulsive behavior, without any justifying
140
explanation. Still, one might reject the physical model presented in this section,
and argue from a rational mechanics points of view. That is, impulses have been
added to our model of rigid body mechanics to deal with situations in which non-
impulsive forces cannot prevent inter-penetration. The principle of constraints is
therefore a good principle of constraints is preferable because it always seeks the
"simplest" solution; that is, it only extends the axioms of rigid body mechanics to
admit impulses when necessary.
One can refute this rational mechanics viewpoint very strongly with the following
argument. Suppose that at time t = 0, we are faced with an indeterminate
configuration, and suppose that there are two behaviors of the system past time
zero. Let us call these behaviors of the systems trajectornes A and B (Figure 12.5).
The principle of constraints asserts that we n?ist always avoid impulsive beha?rior,
whenever possible. As we have seen, it is NP-hard to determine which trajectory (if
either) of A or B avoids impulsive behavior at time zero. But suppose that neither
trajectory A and I? involves impulsive behavior at time zero, and that we arbitrarily
choose A. If at time t > 0 along trajectory A we encounter mconsistency, what
should we do? According to the principle of constraints, it would logical to back up
our simulation to time t = 0 and consider trajectory B,for trajectory B might not
require impulsive solutions for all t > 0. If so, then we see that in fact we really
should have chosen trajectory B to start with, for by doing do, we avoid impulses
forever, past time zero.
Clearly, this is enormously expensive, computationally, because it forces us to
explore potential trajectories, which might themselves branch, to find the trajectory
that is most free of impulses. The possible expense is for all practical purposes
unbounded, since we can make the choice of the correct trajectory depend upon
the outcome of some physical process which is arbitrarily hard to compute. Again,
it can be argued that computational expense should not be a factor in rejecting
inconsistency
no
inconsistency
141
(time)
A
t=O B
Figure 12.5: A configuration with two different non-impulsive solutions at time zero.
If trajectory A is chosen, inconsistency is encountered for some later time t > 0,
but if B is chosen, no inconsistency is encountered for t > 0.
the principle of constraints. However, suppose that an external time-varying force
F(t) acts during the course of the simulation. The decision as to whether trajectory
A or B best avoids impulses past time zero would then require knowledge of the
behavior of F(t) for time I >0. If we d?ose A or B based on F(t) for I >0, then
we violate causality, since events at times I > 0 determine our course of action at
time zero. In fact, we can create a paradox, by letting the fc?rce F(t) be determined
according to which trajectory is chosen, which is itself determined by F(t). Thus,
the principle of constraints in seeking to always use the smallest axiom set for rigid
body med?anics is an "unnatural" principle: it requires a physical system to be
able to look forward in time to determine its action at the present instant.
The above argui?ent clearly generalizes to any physical system with indetermi-
nacy. Any principle that asserts that one sort of behavior is more preferable than
another over the course of time leads to a violation of causality, and is thus an
unnatural principle. That is not to say however that we cannot prefer one behavior
for a particular instant. For instance, when we encounter an indeterminacy in a
simulation, we might be told that one trajectory is more probable than another. If
we regard indeterminacy as simply inade?iate knowledge of a system's state to the
142
point that we cannot predict its actions, we might in fact be given a "tip" as to the
correct behavior of the system by an oracular observer who has more knowledge of
the system's state than we do. Although we should regard all trajectories as equally
possible, they need not all be e?ially likely.
Given that we are willing to accept both impulsive and non-impulsive solutions
as valid solutions for a configuration with dynamic friction, the following question
arises. We know that searching for a particular type of solution, namely solutions
without impulsive forces, is NP-hard. How difficult is it to find either a solution with
or without impulsive forces? That is, how much time is required of an algorithm
that computes either an impulsive or a non-impulsive solution for a configuration,
without preference? Such an algorithm could not be used to decide if a configuration
was consistent, since the algorithm could return an impulsive solution whether or
not the configuration was consistent. Tlie use of such an algoritlim to compute
contact forces (and thus choose an arbitrary trajectory for the simulation) would
avoid the NP-complete problem of deciding consistency. In the next chapter, we
will consider an algorithm with this property.
Chapter 13
Computing Valid Contact Forces
and Impulses
Before computational methods can be considered, a definition of validity for contact
impulses is needed. In order for a definition of validity to be useful, all inconsistent
configurations should have a set of contact impulses satisfying the definition. An
analysis of the computational method presented in Section 13.2 will show that the
definition of validity presented satisfies this condition.
13.1 Defining valid contact impulses
In the penalty method interpretation of Figure 12.1 an impulse occurred because
no matter how strong the normal force became, it was insufficient to prevent inter-
penetration. As a result, after the contact impulse was applied, the relative velocity
of the bodies at the contact point was directed inwards. Since contact impulses may
need to be applied to configurations involving more than one contact point, validity
must be defined for a set of contact impulses. For example, in Figure 12.4, if the
fN = m1 behavior is chosen, a contact impulse should occur between E1 and 1?.
143
144
However, there should be no contact impulse between R1 and B.
We will call a set of contact impulses valid under the following two conditions.
The first condition is that the contact impulses must convert at least one of the
contact points with dynamic friction to static friction. The reasoning here is that
the impulse has occurred because the initial conditions of normal forces have caused
inconsistency, as in Figure 12.4, and the inconsistency is removed by converting a
contact point from dynamic to static friction.
The second condition is that every contact point at which a contact impulse
occurs must end up with a non-positive relative normal velocity; that is, after the
contact impulses are applied, bodies should not be separating wherever contact
impulses occurred. The justification for this is that the contact impulses occur only
when the normal force grows without bound to oppose inter-penetration. Intuitively,
valid contact impulses are the limiting result of increasing normal forces without
bound under the penalty method. If bodies are separating at a contact point after
contact impulses are applied, then the normal force at the contact point should not
have grown without bound into a contact imp?ilse.1
1 For planar configurations, this definition is valid, because a contact with dynamic friction
cannot change its direction of sliding without having zero sliding velocity. For three-dimensional
configurations, the sliding velocity direction can d?ange without a contact changing from dynamic
to static friction. Thus, for three-dimensional configurations, ?? would want to find contact
in?pulses that satisfy the second condition. ?Ve would then need to differentially apply the contact
impulses, and account for the change of sliding velocity direction; essentially the same problem as
dealing with fnctional impulses for three-dimensional collisions. Thus, the results of this chapter
are valid for planar configurations, and are a starting point fi?r three-dimensional systems. If one
adopts the ad hoc assumption that the sliding velocity direction does not change during a collision
for three-dimensional systems, then the results of this d?apter can be applied without complication
to three-dimensional systems.
145
13.2 Lemke's algorithm
Having defined validity for contact impulses, we can consider computational methods
for computing contact forces and impulses. Given that computing contact forces
alone is hard, how can either contact forces or impulses be computed efficiently?
We begin by considering the LCP of Equation (8--H7). One of the first algorithms
for solving linear complementarity problems was introduced by Lemke[27] and is
known as Lemke's algorithm. Lemke's algorithm is a pivoting method, similar to
the simplex method of linear programming. The algorithm is exponential in the
worst case, but has an expected running time polynomial in n. Lemke's algorithm
progresses, like the simplex method, by trying various descent directions. If an LCP
is PSD and has no solution then Lemke's algorithm will at some point encounter
an unbounded ray; that is, a descent direction along whidi one can travel infinitely
far without making any progress. Otherwise, if a PSD LCP has a solution, then no
unbounded ray exists for that LCP, and Lemke's algoritlim terminates by finding a
solution to the LCP. The algorithm is viewed as a practical solution method to the
problem of solving PSD LCP's.
However, for non-PSD LCP's, Lemke's algorithm is not guaranteed to terminate
correctly (although it still takes only expected polynomial time to do so). For a non-
PSD LCP, if there is no sollition, Lemke's algorithm terminates by encountering an
unbounded ray. Unfortunately, if there is a solution, the algorithm is not guaranteed
to find it. For non-PSD LCP's with sohitions, Lemke's algorithm terminates either
by finding a solution or by encountering an unbounded ray.2 As a result, Lemke's
algorithm is not suitable for solving non-PSD LCP's.
However, whenever Lemke's algorithm terminates by encountering an unbounded
2 Encountering an unbounded ray when there is a solution is roughly analogous to getting stuck
at a non-global minimum in a non-convex minimization problem.
146
ray, it has found an n-vector z with the property
(1 < i < n)			(13--H1)
z > 0			and			Zj > 0 h?phes (Az)j < 0
where (Az)i is the ith component of the vector Az. Why is this property of interest?
Suppose we let z represent a set of contact impulses in the same way that
!N represents a set of contact forces, and apply those contact impulses to the
configuration. By the definition of impulse, the d?ange in relative normal velocity
at each contact point is (Az)i. Since the relative normal velocity at each contact
point is assumed to be zero before application of the impulse, the vector of relative
normal velocities after applying the impulse is simply Az.
Thus, the contact impulses indicated by z satisfy the second condition of validity:
every contact point subject to a non-zero contact impulse Zj > 0 ends up with a
non-positive relative normal velocity (Az) < 0 In order to completely satisfy the
definition of validity, z must be scaled upwards from zero until it causes a contact
point with dynamic friction to be converted to static friction. A subtle point to note
is that if z specified contact impulses only at the frictionless contact points, then z
could not satisfy the first condition condition for valid contact impulses (although
it might satisfy the second). In order for z to represent valid contact impulses, z
must have the property that for at least one 1 < i < n the ith contact point has
dynamic friction and Zj > 0. A proof that Lemke's algorithm either solves the LCP
or finds a vector z such that oz represents valid contact impulses for some ? > 0 is
given in Appendix B. After application of the contact impulse represented by cz,
some real collisions (that is, collisions with positive normal approach velocity) occur
and must be dealt with.
The behavior of Lemke's algorithm exactly matches our new view of the problem
of computing contact forces. If the configuration has no valid contact impulse
solutions, Lemke's algorithm cannot terminate with a ray z and must therefore
147
find a valid contact force solution. For inconsistent configurations, no valid contact
force solution exists, so Lemke's algorithm must terminate with a ray z, providing
a contact impulse solution. For configurations with both a valid force and impulse
solution, Lemke's algorithm will terminate by computing one or the other.
Note that the definition of valid contact impulses does not add any new solutions
to frictionless systems. To see this, consider an arbitrary frictionless system, and
suppose that a vector z satisfies the condition of Equation (13--H1). Then
zTAz = ? zi(Az)i ? 0.			(13--H2)
i=1
But since the system is frictionless, A is PSD. Since a PSD matrix satisfies
xTAx > 0 for any vector x, Equation (13--H2) implies that
zTAz = 0.
(13--H3)
Since A is PSD, this in turn implies that Az = 0. Thus, any vector z that
fits the first condition of valid contact impulses is a vector in the null-space of A.
Since frictionless systems are determinate, given t?? vectors x and y such that
Ax = Ay, application of the impulses (or forces) indicated by x has the same
effect as the application of the impulses (or forces) indicated by y. Since Az = 0,
this implies that the application of the impulses indicated by z has absolutely no
effect at all on the system. Thus, the definition of valid contact impulses is such
that it adds no new solutions to frictionless systems.
Although Lemke's algorithm runs, practically speaking, in polynoniial time, this
is not a proof that finding either valid contact forces or impulses is a polynomial time
problem. From a practical standpoint, though, Lemke's algorithm provides an effi-
cient algorithm for computing valid contact forces or impulses. The computational
complexity of either solving an LCP or finding a vector z satisfying Equation (13--H1)
is unknown; but it would be surprising if this could not be done in polynomial time.
148
13.3 Continuity of solutions
One implementation issue we need to be concerned with is the continuity of the
solutions chosen for !N. Consider a configuration with two possible solutions !N
that produce different accelerations of the bodies. Suppose that we (perversely)
alternate between the two solutions at each time step, causing the accelerations of
the bodies to be discontinuous. This will cause the ordinary differential equation
solver to break down, since it expects to be integrating a continuous function.
In Section 8.1 where configurations were frictionless, this problem did not occur
because all solutions of Equation (8?7) produced the same acceleration. ?7hen
friction is introduced this is no longer true, and care must be taken to make !N (at
least) piecewise continuous over time. Recall that Section 8.2 introduced the concept
of computing !N by solving a non-singular linear system of equations as indicated by
an index set L. This computation method for !N may be regarded as a time saving
luxury for frictionless configurations, but it is essential in simulating configurations
with friction. Since A and b depend on spatial and velocity variables, they are
continuous functions of time (in the absence of impulses). As a result, the system of
equations produced from L is continuous as long as L remains constant. Because
of this, the !N vectors computed with this method are piecewise contimious, with
discontinuities occurring whenever L changes. The integrator is of course informed
whenever L changes.
Chapter 14
Static Eriction
The results in this thesis concerning the complexity of static friction are mud?
less comprehensive than the results for dynamic friction. ?Ve will indulge in
some speculation concerning complexity results for static friction, and then discuss
some methods for approximating static friction forces for simulation purposes.
The approximation methods are suitable for some applications (such as computer
animation), but may be unsuitable with regards to other applications.
14.1 Complexity of static friction
Consider once more the ith contact point of a configuration. The normal force
and tangential friction force magnitudes at this point are fivi fri and fyi' The
corresponding acceleration magnitudes are ?Ni, a?i and a??. Since we already
know that the presence of dynamic friction leads to NP-haulness, let us study
configurations with only static friction and no dynamic friction.
In Chapter 9, the conditions for static friction were presented. These conditions
149
150
can be written as
fN >0 and
= a?? = 0 or
aN' >0 and
and either			j?2j + fy2j = (IifNi)2 and			(14-1)
fNiaNi = 0 and
fr2i + f?j < (IifNi)2 f??a?? + f??a?? < 0.
The first questions we might ask are: does Equation (14--H1) always have a solution?
Does Equation (14--H1) always have a unique solution? L5tstedt[28] gives an example
of a one contact pohit configuration with static friction that does not have a unique
solution. The example given is the sliding rod of Figure 10.2, but with the linear
velocity adjusted so that the tip of the rod is motionless with respect to the ground.
Once again, the rod is free to either rotate off the ground or remahi in contact.
Thus, Equation (14--H1) does not necessarily have a unique solution. NVhether
Equation (14--H1) always has even one solution is unknown. So far, no examples
of configurations with static friction resulting in inconsistency are known (excepting
the trivial case when the configuration has unsatisfiable kinematic constraints). This
is particularly striking in light of the fact that inconsistency for dynamic friction
can occur for one contact point configurations. For static friction, we can however
prove that every one contact point configuration has a solution. The proof technique
makes use of geometrical constructs based on work by ?Vang and N'Iason[34].
THEOREM 4. All two-dimensional one contact point confignrations with static
friction have valid soltttions.
PROOF. For two dimensional configurations, we can describe the frictional force
and the tangential accelerations in terms of just fxi and a?i. The analogue of
151
Equation (14--H1) in two dimensions is therefore
fNi > 0 and			ari = 0 or
?Ni > 0 and
fNiaNi = 0 and			and either			Ifxil = ?fNi and			(14-2)
IifNi < fxi < f'fNi			fr?azi <0.
For a one contact configuration, we can express the normal acceleration aN and the
frictional acceleration a? in terms of the normal force fN and the frictional force
fx as
aN = BfN + Cfx + F
= CfN + Efx + G
where B, C, D, and E satisfy[28]
B>O, E>O and BE>C2.
(14--H3)
(14--H4)
Additionally, we can assume C < 0 by reversing the sense of the tangent direction
along which fx and a? are measured, if necessary. (If C = 0 then Equation (14--H3)
decouples and it is trivial to show Equation (14--H2) always has a solution.) From
this, it is easily shown that
--HE			--HC			(14--H5)
C			B
Using these facts, we show that Equation (14--H2) always has a solution.
If F > 0 then Equation (14--H2) has the trivial solution Lv = fx = 0 Otherwise,
F < 0. In this case, we claim a solution with Lv > 0, aN = 0 and one of the
following three cases always exists:
= 0 and --H[`fN < fx < IlfN
(ii) a? > 0 and f1 = --HJtfN
(iii) a? <0 and fr = [fN
152
(a)
N
fN
fx
N
jT
fx
Figure 14.1: Force space. (a) Possible solutions for P lie on N between L and RL.
(b) Possible solutions for P lie on N to the right of N's intersection with RL.
To verify this claim, we use a geometrical method based on the impulse space
described by Wang and Mason[34].
Consider Figure 14.1 where fN is plotted vertically arid fx is plotted horizon-
tally. We are seard?ing for a point P = (fN, fx) in this space that satisfies both
Equation (14--H2), and fN > 0 and a1v = 0. Several lines are identified in this space.
First, the line BfNtCfx+F = aN = 0 is labeled N. Next, the lines fx = [tfN and
fx = --H!tfN are labeled L and fiL respectively. Since F ? 0, B > 0 and C < 0
AT is either horizontal or has positive slope and crosses the Lv axis above the fx
axis. Because of this, N always iiitersects I?L and may (Figure 14.1 a) or may not
(Figure 14.1b) intercept L. In order for a point P to satisfy fN > 0 aN = 0 and
--H[tfN < fx < fN, P must lie on N above the horizontal axis, and on or between
L and RL. Points satisfying these conditions are indicated by the bold portion of
line N in Figure 14.1.
Suppose now that N crosses both L and RL, and consider Figure 14.2. The line
153
(a)
RL
______ _______			P2:;
5,
s
fx			fx
1fN
(c)
fN			fx
Figure 14.2: Line N intersects L and RL. (a) S intersects N and point P1 satisfies
case (i). (b) s lies above N between L and RL. Points lying below 5 have a? > 0;
point P2 satisfies case (ii). (c) 5 lies below N between L and RL. Points lying
above 5 have a? < 0; point P3 satisfies case (iii).
CfN + Efx + U = a? = 0 is labeled 5. Any point P on 5 results in zero tangential
acceleration. From the definition of 5 and since C ? 0, points lying above 5 satisfy
< 0 while points lying below 5 satisfy a? > 0. Suppose 5 intersects N between
L and RL (Figure 14.2a). Then the intersection point Pi falls under case (i) and
is a solution. If however 5 does not intersect A? between L and RL, then either 5
lies above N as in Figure 14.2b, or below N as in Figure 14.2c. For Figure 14.2b,
the point P2 is a solution under case (ii) since it satisfies fx < 0 and a? > 0. In
Figure 14.2c, P3 satisfies fx > 0 and a? ? 0 and is a solution under case (iii).
The only situation left to consider is if AT intersects only 1?L and not L
(Figure 14.3). Suppose that 5 lies above N's intersection with RL (Figure 14.3a).
Since points below 5 satisfy a1 > 0, the point P4 is a solution under case (ii)
Otherwise, 5 lies below N's intersection with i?L (Figure 14.3b). Since 5 has
slope --HE/C and N has slope --HC/B, from E?iation (14--H5), 5 has a steeper slope.
Therefore, 5 must intersect AT somewhere between  and i?L at P5; this solution
154
(a)
Ar
RL
N
fx
s
fN			fx
Figure 14.3: Line N intersects RL, but not L. (a) S lies above the intersection of
N and RL; point P4 satisfies case (ii). (b) S lies below the intersection of N and
RE. Since S's slope is greater than N's, 8' must intersect N somewhere between
L and fiL. The intersection point, P5 satisfies case
falls under case
Having verified our claim that that for F ? 0 one of cases (i), (ii) or (iii) always
holds, we conclude that all two-dimensional one contact point configurations with
static friction have valid solutions. E]
This result is easily exteudible to three dimensions. The same geometric intuition
is used, except that N and 8' are planes and the lines L and 1?L are replaced with
the cone f?2 + = ([tfN)2.
`A7e can speculate that Equation (14--H1) does in fact always have a solution.
Having shown that all one contact point configurations with static friction have a
solution, we can contemplate an inductive proof that all n contact point configura-
tions with static friction have a solution. Theorem 4 would be the base case, and
the inductive step might be as follows. Assume all n contact point configurations
with friction have solutions, and consider an n + 1 contact point configuration.
155
Let one contact point be arbitrarily called the "new" contact point. We will insert
an external force at this new contact point that acts equally and opposite on the
contacting bodies at this point, Henceforth we ignore this new contact point, and
consider the configuration to be simply an n contact point configuration. For any
external force we insert at the new contact point, there is a valid solution of contact
forces at the other n contact points, since all n contact point configurations have
solutions.
The idea then is that we slowly increase the strength of the external force acting
at the new contact point. As we do this, the contact forces at the other n contact
points continually readjust to form a valid solution. Eventually, we should reach
a point where the external force at the new contact point is just sufficient to
prevent inter-penetration without being strong enough to cause separation at the
new contact point. At this point, if we treat the external force as a contact force
at the new contact point, we will have a valid solution for the n + 1 contact point
configuration. This will complete the indiictive step. However, the above is only
speculation; there are still some unresolved difficulties associated with this proof.
What can we actually say about the complexity of computing contact forces?
First, note that the problem of computing contact forces for two-dimensional config-
urations with fiiction lies iii NP. Just as any sohition to the LCP of Equation (8--H7)
could be written in terms of the solution of a linear system, the same is true for a
solution to Equation (14--H2). As a result, any sollition of Equation (14--H2) must be
rational and of polynomial length[52]. Thus, finding a solution to Equation (14--H2)
lies in NP since a solution to the equation can be verified in polynomial time.
Unfortunately, for thre&dimensional configurations, the Coulomb friction con-
straint fx2i + < (ItfNi)2 brings up the possibllity that a solution might in-
volve irrational quantities. We therefore cannot say that finding a solution to
Equation (14--H1) lies in NP. Practically speaking however, we would not expect
156
finding a solution to Equation (14--H1) to be any harder than finding a solution
to Equation (14--H2). Laying aside the fact that solving Equation (14--H1) might
not lie in NP, can we at least express solutions to Equation (14--H1) in terms of
a convex minimization problem, as was done for frictionless configurations? Convex
minimization problems, even when not lying in NP, usually admit practical solutions
by some sort of descent procedure. Even if we could not exactly classify the difficulty
of solving Equation (14--H1), phrasing it as a convex minimization problem would go
a long way towards describing the inherent complexity. The answer to this question
is, unfortunately, "no" - The configuration described by L5tstedt[28] that does not
possess a unique solution has in fact exactly two distinct solutions. Since the solution
set of a convex minimization problem is always a single connected component, the
solution set of Equation (14--H1) cannot be posed in terms of a convex minimization
problem.
Should we believe then that computing contact forces for configurations with
static friction is NP-hard? NVe cannot say, but we can again speculate. If we were
able to find a configuration that had no solution, then we could probably form a
reduction similar to the reduction of Chapter 11 to show that computing contact
forces is NP-hard. Suppose however that our speculation that configurations with
static friction always have a solution is correct. As previously mentioned, computing
contact forces for two-dimensional configurations lies in NP. If all configurations
with static friction have solutions, the decision problem of whether or not a
two-dimensional configuration with static friction has a solution is trivially in P
(polynomial time), since the answer is always "yes". Would this mean that contact
forces could also be computed in polynomial time?
So far, an example of an NP decision problem lying in P that requires more
than polynomial time to actually find a witness is unknown. Indeed the existence
of such a problem would constitute a proof that P # NP, the reasoning being that
157
if P = NP, all NP witness-finding problems can be solved in polynomial time. If
in fact all configurations with static friction had solutions, showing that contact
forces could not be computed in polynomial time would be equivalent to showing
that P $ NP. Thus, if all configurations had solutions, it would probably be more
fruitful to attempt to construct a polynomial time algorithm for computing contact
forces than to attempt to disprove the existence of such an algorithm.
14.2 Approximation methods
We will conclude our discussion of static friction by considering three different
methods for approximating static friction forces. The first method is due to
L5tstedt[30] and appears to be the only published approximation method for contact
forces with static friction that does not involve the penalty method. The second
two methods are new. For all three methods, we express tangential and normal
accelerations as linear functions of the contact forces by writing
where
and
a = Af --H b			(14 6)
a = (CIArl, ?`i, (z??, . . . , av?, a1?,
f = (fNi,fxi,fy1,. .. ?fNn?fTfl?fyn)T
The matrix A is a 3n x 3n matnx and b is a 3n-vector.
(14--H7)
(14--H8)
14.2.1 L5tstedt's approximation
Lo?tstedt's approximation is a very simple modification of the Coloumb law of
friction, which very cleverly achieves two difficult goals. First, it allows dynamic
158
friction problems to be solved by convex QP's or LCP's. Second, it converts the
static friction conditions into a form that is also solvable by quadratic programming;
moreover, the quadratic program produced is convex. This is a very elegant
transformation!
For dynamic friction, L5tstedt sets the dynamic friction force to have a magni-
tude of ? times the normal force magnitude of the previous time step. This makes
the dynamic friction force a known external force at the current time step. Since
only unknown normal forces and known external forces act on the bodies now, the
configuration has a convex LCP. Thus, L5tstedt's modification does away with both
indeterminacy and inconsistency, and produces an LCP that can be efficiently solved
(since the LCP is convex).
For static friction, L6tstedt sets the upper bound on the magnitude of the static
friction force to j times the normal force magnitude of the previous time step. This
allows static friction and normal forces to be solved by quadratic programming in
the following manner. For a given time step, let fNi denote the normal force at the
ith contact point from the previous time step. L6tstedt's modification produces a
system of equations of the form
fNi > 0 and
?Ni > 0 and
fNiaNi = 0 and
Ii'fNi < fxi < liLvi and (14--H9)
?[LfNi ? fyi < I'fNi and
??i(jifNj --H Ifxil) = 0 and f??a?? < 0 and
ayj(JifNi --H Ifyji) = 0 and f??a?? ? 0.
The approximation is appealing because it still enforces a very important property
of static friction: the friction force achieves its upper bound and is opposite the
tangential acceleration whenever it is unable to prevent sliding. Still, the above
159
conditions still appear fairly complex, and certainly not the conditions of a quadratic
program.
It turns out however that the QP
fNi >0
minimize ?2l!TA! --H bT! subject to			P'fNi < fxi < fNi			(14-10)
lilNi < fyj < fNi
has NKT conditions that matd? the conditions of Equation (14--H9), so a solution
of Equation (14--H10) satisfies the conditions of Equation (14--H9). Thus, the fairly
complex conditions of Equation (14--H9) have the highly desirable property of being
solvable by quadratic programming. (Note that Equation (14--H9) does not appear to
be solvable as an LCP though.) Unfortunately though, when A is not symmetric,
Equation (14--H10) has an entirely different set of I&I&T conditions that do not match
Equation (14--H9). For configurations with dynamic friction, even if A is still PSD
(as often happens when the frictional coefficient [I is small) it is rarely (if ever)
symmetric. This means that we cannot easily extend this approximation to handle
the unmodified dynamic friction law, because dynamic friction causes A to become
asymmetric.
Additionally, there is the question of the decoupling of the t4 and t% axes, and
tlie replacement of fNi with fNi. The effects of decoupling of the axes can possibly
be minimized by a wise choice of vectors tx and Jy for each contact point (so that
one of fxi or fyi is near zero, for each contact point), but the replacement of fNi
with fNi is more serious. It is unclear how one produces a suitable fNi for the initial
time step (or after a discontinuity in the configuration's state). One possibility is to
pick arbitrary initial values for the fNi (such as zero), solve the system, and update
fNi with the newly computed fNi, and iterate. However, this is method will not
necessarily converge; for the one-point configuration of Figure 10.1, the magnitude
fNi will in fact diverge to infinity. Likewise, consider Figure 10.1 with a reversed
160
tangential velocity, so that the friction force is in effect accelerating the contact
point Pa "upwards". Suppose fNi is set to zero, and fNi is solved for; on the next
iteration, we set fNi = fNi Depending on the value of g and the value of 0, this
new value of fNi could be sufficiently large that it causes contact between the rod
and the ground to break. If this happens, the new value of fNi will be zero, since
contact is breaking, and the following iteration will start with fNi set to zero as
well. The iterative process then contimies, with fNi and fNi oscillating between
two distinct values.
It is difficult to quantitatively describe the effects of the approximation. In
situations where neither indeterminacy nor inconsistency occurs, perhaps the ap-
proximation has little effect on the outcome of a simulation. As a trivial example, if
we set j = 0, then the approximation has no effect at all. ?`e can imagine that as
we slowly increase p?, the effects of the approximation might not be apparent until
!l takes on a suitably large vahie. But for situations with inconsistency, Lo'?tstedt's
approximation produces a non-impulsive solution where such a solution did not exist
previously. In this situation, the approximation must clearly have a large effect on
the sii?ulation.
14.2.2 Dynamic friction approximation
The second approximation method we consider has the following virtues. It is
extremely sii?ple to implement, very robust, and handles dynamic friction. The
approximation is similar in concept to the penalty method of computing normal
forces based on inter-penetration depth. It avoids some of the difficulties that make
the penalty method unsuitable as a simulation method, but it does share the major
flaw of all penalty methods: it is an ad hoc description of a physical process.
The second method is as follows. In order to determine whether static friction or
161
dynamic friction should occur at a contact point, the simulation program must have
some threshold value E. If the tangential velocity (Vx2i + vy2?)1/2 > E then dynamic
friction occurs. Otherwise static friction occurs. Since the dynamic friction force is
calculated by
fxi = !ifNiVxi			fyi = IifNiVy?
#-+vy2?			??2i+vy2?
we approximate static friction as
$T?it%?y2?)
E
and			(14--H11)
fxi =			?fNi			?Vxi			= flfNi?VX?			(14--H12)
7m+vy2?
and similarly for fyi Thus, we approximate static friction as a dynamic friction
force that varies in magnitude from zero to an upper limit of IlfNi as the tangential
speed varies from 0 to c. The n?thod of Section 13.2 can be used to compute the
contact forces, since this is now entirely a dynamic friction problem.
In this approximation however, static friction occurs only when the tangential
velocity is non-zero. Bodies must acquire some small amount of "crawl" in order to
maintain a static friction force. This is similar to the penalty method, where bodies
must acquire some degree of inter-penetration for a sufficient normal force to exist.
However, in the penalty method, it is necessary to increase the spring constant I?
(Chapter 12) without bound as the mass of bodies increases. Our approximation
method does not suffer from this problem. If ( is made small enough, the "crawling"
behavior of bodies is not visible, no matter what masses or forces or exist. If ?
is made excessively small, the differential equations of motion may become stiff
otherwise, the method has a reasonable performance. The graphical simulations
produced using this approximation are fairly realistic.
162
14.2.3 Quadratic programming approximation
The last method we present attempts to model static friction by quadratic program-
ming. The basic idea is that by guessing the signs of the variables fxi and fyj in a
solution, we can solve for ! by solving a QP. Given a guess of the values
sgn(f??)			and			sgn(fyj)			(14--H13)
where sgn(r) = 1 if x > 0 and --H1 otherwise, we will try to find a solution f that
agrees with these signs. Thus, we treat sgn(f??) and sgn(fyj) as known and fixed
quantities for now. We will describe later how we produce a guess of these values,
and then discuss some shortcomings of the method.
Because we canuot represent the constraint fx2 + f2 < (liTNi)2 in a QP, we
decouple the t4 and t% axes as L5tstedt does, replacing fx2i+fy2 < (j"fNi)2 with
--H ?fNi < fxi < J1fNi			and			--H I'fNi ? fyj < liLvi			(14--H14)
Note however that for two-dimensional configurations, no decoupling of axes is
needed and this method does not making any approximations of the laws of static
friction.
?`e will also replace the condition axifri + ayjfyj < 0 with the couditions
f??sgn(f??) > 0			and			a1?sgn(f??) < 0			(14--H15)
f?,jsgn(fyj) >0			and			ay?sgn(fyj) <0.
This condition ensures that a??f?? + ayjfyj < 0 (and is actually somewhat stronger).
Note that the (seeming) tautologies fxi sgn(f1j) > 0 and fyi sgn(fyj) > 0 are
actually constraints that fxi and fyi have the same signs as sgn(f??) and sgn(fy?)
(which we are regarding as known quantities). Next, we add the standard constraint
that
fNi>0			and			a??>O			(14--H16)
163
Since we are treating sgn(f??) and sgn(j?,?) as known quantities, and ?Ni, a?? and
ay? are linear functions of the contact forces, all of the above conditions can be
expressed as linear constraints involving fxi, fyi and fNi'
The condition that the static friction force attains its upper bound when slipping
begins is written as
(I"fNi --H fxi sgn(f1?))(--Ha?? sgn(f??)) = 0
(!tfNi --H fyi sgn(fy?))(--Hay? sgn(fy?)) = 0.
Finally, we add the standard constraint on the normal forces that
(14--H17)
(fNi)(aNi) = 0			(14--H18)
Since each of fNi, 0Ni, !tfNi --H f??sgn(f??), [?fNi --H fy?sgn(fy?), --Ha1?sgn(f??),
and a?? sgn(fy?) are constrained to be positive, we can express the fact that
Equations (14--i 7) and (14--H18) hokl for all i as
fNi0Ni + (!"fNi --H fxi 5gn(fxi))(?Qri sgu(f??))
= 0.			(14--H19)
+(jifNi --H fyi sgn(fy?))(--Haxi sgn(fy?))
Since this is a quadratic function of fxi, fyi and fNi, we can find a solution to
Equation (14--H19) satisfying the linear constraints of Equations (14--H14), (14--H15)
and (14--H1(3) by solving a QP.
NVe still must describe how to guess the signs of fxi and fyi Iterative methods
for quadratic programming exist that are very similar to the Gauss-Seidel or Jacobi
iterative methods used to solve linear systems[40]. These iterative routines are mud?
simpler than routines that solve QP's by direct methods; in particular, it is extremely
simple to modify these iterative routines to find a solution to Equation (14--H1). The
reason that we do not use these modified iterative routines is that they have poor
convergence properties. However, we can run one of these modified iterative routines
for some period of time and examine the (hopefully) partially converged variables
164
fxi and fyi to guess their signs. (Appendix C discusses how to formulate an iterative
routine to solve Equation (14--H1).) Using the signs of the variables, we then form
and solve a quadratic program as described above. Coherence can be exploited (and
must, to provide continuity over time steps) by using the same concept of basic
solutions presented in sections 8.2 and 13.3.
This method can break down at any number of places. First, we still do not
know that a solution for static friction always exists. Even if there is a solution
it is possible that the iterative algorithm may not converge to a solution, so our
guess of the signs of fxi and fyj may not fon? a QP that leads to a solution.
Furthermore, the QP may be non-convex, and therefore impractical to solve (even
if only static friction occurs). Last, if we attempt to incorporate both dynamic and
static friction, we still need to be concerned about the NP-hardness due to just the
dynamic friction. Unfortunately, the form of the linear constraints of tile QP do not
allow us to solve it as an LCP using Lemke's method. No quadratic programming
algorithm is known that has tlie property of Lemke's algorithm: returning a valid
non-impulsive solution or indicating a valid impulsive solution. It is possible that
whenever the iterative algorithm fails to converge, it does so in such a way that an
analysis of the divergence indicates a valid set of contact impulses; it is not clear
how to perform such an analysis.
However, this last method yields a very acceptable visual result, when it works.
For large numbers of contact points (n 40), the last method sometimes breaks
down; either the iteration does not converge, or it converges to an incorrect solution.
For configurations with approximately 40 contact points or less, the method is
reasonably successful. Using this method, simulations with n 40, and involving
thousands of time steps have been riin successfiilly. Images from several such
simulations are shown in Appendix D.
Chapter 15
Future Directions
In this thesis, we have used a model of rigid body dynamics that is greatly simplified
to examine the difficulties of simulation of complex physical systems. Even using
a simplified model of rigid body dynamics, this thesis has shown that the problem
of dynamic simulation with non-interpenetration constraints is a difficult problem.
Beginning with the dynamics of frictionless configurations with a single point of
contact, we have seen that producing a computationally tractable solution method
for the contact force at the single contact point is non-trivial. Next, we considered
properties and solution methods for configurations with multiple contact points and
friction. The most difficult problems we encountered have centered on the behavior
of configurations with friction. ?N7e have shown that the computational complexity
of computing a non-impulsive solution for configurations without colliding contacts
is NP-hard; more importantly, however, we have shown that it is not necessarily
important to specifically compute such a solution, and that an absolute preference
for such solutions over the course of a simulation violates causality and is thus
unnatural. A model that allows expected polynomial time solution methods has
been proposed, and iterative computational methods for computing contact forces
165
166
for configurations with and without collisions have been proposed.
A large number of graphical simulations have been produced using the methods
described in this thesis, and images from several simulations are shown in Ap-
pendix D. Although these simulations yield results that seem visually correct, they
are still simulations of a system subject to a number of simplifying assumptions, or
restrictions. How severe have these restrictions been, and what is required to relax
them?
The major restriction in the work described in this thesis has been the assump-
tion that contact between bodies occurs at a finite nun??er of contact points. For
simulations with frictionless polyhedra, this restriction introduces no approximations
into the simulation, since frictionless systems are determinate. However, this
restriction prohibits simulations involving contact areas with curved boundaries,
such as a cylinder standing upright on a plane. The problem with such a simulation is
how to formulate the non-penetration constraint: since the cylinder could potentially
be tipped over from any direction, every point on the circular boundary of the
contact region (between the plane and the cylinder) must be included in a non-
penetration constraint. Currently, it is unclear how one should write or attempt to
solve the constraint equations for this case.
Interestingly, the same problem arises when considering the formulation of
non-penetration constraints between flexible bodies that contact with a one- or
two-dimensional contact region. Even if the contact region is polyhedral, inter-
penetration is not necessarily prevented by constraints on the allowed motion of the
vertices. A deformation that causes a fiat surface to become curved could result
in inter-penetration at an interior point of a contact region even though no inter-
penetration had occurred on the boundary of that contact region[2]. Additionally,
dealing with one- or two-dimensional contact regions also requires research into
contact determination algorithn? that can both detect and describe, in some fashion,
167
Figure 15.1: Non-polygonal contact region. The contact area between the cube
and the plane is a polygon, while the contact area between the cylinder and the
plane is a circle. The non-penetration constraint fi?r the cube is posed in teri?s of
the finitely i?any vertices describing the contact area boundary. For the cylinder,
infinitely many points are necessary to describe the contact area boundary.
168
one- and two-dimensional contact regions. Finding collision/contact determination
algorithms that can also deal with non-convex curved surfaces is also an interesting
research problem.
For simulations with friction, the restriction to finitely many contact points is
more severe. In essence, the tangential friction forces that should occur over a line
or area of contact are restricted to act only at the vertices of the contact region.
Again, the problem is how to fbrmulate the Coulomb friction laws over an entire
region of contact; this is a very open research problem.
Another problem is computing contact forces and impulses in the presence of
friction. The assumption of constant slip velocity direction for three-dimensional
collisions simplifies comp'itations, but is unrealistic and must be dealt with. Also,
?? do not have a guaranteed sol?ition method for configurations involving both static
and dynamic friction. There is also tlie philosophical issue of what a simulation
should do when faced with indeterminacy: this problem has been largely (if not
completely) ignored by this thesis. There has been some speculation that both
indeterminacy and inconsistency due to friction could be eliminated by discarding
the restriction that bodies are perfectly rigid. Depending on the model of non-
rigidity adopted, inconsistency and indeterminacy can either be partly or completely
eliminated. For example, if a standard mass point and spring system is used to model
deforiuable bodies, both inconsistency and indeteriuinacy are completely eliminated.
In such a model, a given contact force affects the instantaneous acceleration of only
a single mass point; the change in position of the mass point indirectly affects neigh-
boring mass points via the spring forces. Since a mass point's instantaneous normal
acceleration is unaffected by any tangential friction forces, neither indeterminacy nor
inconsistency is possible. Note that this model of flexibility eliminates inconsistency
and indeterminacy by giving the body an essentially unlimited number of degrees of
freedom (3n degrees of freedom for a body with n mass points). There also exists
169
a middle ground, in which bodies are allowed to deform, but with only a limited
number of degrees of freedom. Witkin[56] and Baraff and Witkin[2] discuss flexible
bodies that are limited to linear or quadratic deformations of a rest shape. Under
this model, inconsistency and indeterminacy are only partly eliminated; the limited
flexibility of the body eliminates inconsistency and indeterminacy in configurations
with one contact point, but allows inconsistency and indeterminacy in configurations
with two or more contact points. In general, one could speculate that the more
degrees of freedom a model permits, the less likely indeterminacy and inconsistency
can occur.
Appendix A
Extremal Point Conditions for
Parametric Surfaces
The derivations in this thesis have all assumed implicit descriptions of curved
surfaces. In this appendix, derivations for surfaces described parametrically are
given. NVe will begin by defining the parametric analogue of Equation (5--H2).
Let the surfaces A and B be modeled by the time-varying parametric functions
S(u, v, t) and T(?, v, t) - Componeiftwise,
Sx(?i, t?, t)			Tr(u, v,
S(?t, v, =			Sy(?t, v, I)			and			T(u, v, =			Ty(u, v, 1)			(A--Hi)
Sz(u, v, t)			Tz(u, v,
A normal vector SN(v,v,t) is given by
5N(n, v, t) = aasu (v, v, t) x
v, t)
and similarly for T. If the extremal points Pa and Pb at time t are defined in terms
of parametric coordinates (Va, Va) and (ub, vb) by
Pa = S(ua,va,t) and Pb = T(vb,vb,t)
171
172
then the analogue of Equation (5--H2) is
E1 :			SN(Ua,Va,t) t A2TN(ub,vb,t) = 0			(A--H4)
E2 :			(T(ub,vb,t) --H S(Ua,%'a,t)) t AlTN(ub,vb,t) = 0.
The analogue of Equation (5-4) is
SN(ua,va,t)T? = 0			(A--H5)
(S(ua,va,t) --H (P0 t cr)) + AlSN(fla,Va,t) = 0,
while the analogue of Equation (5--H5) is
(S(ua,va,t) --H Po) + AlSN(Ua,Va,t) = 0.
These equations can be solved to find the extremal points the during collision/
contact determination step.
In deriving an expression for C for force determination, the functions 5 and T
are treated as configuration dependent functions S(u, v, p) and T(u, v, p?-). Equa-
tion (A--H4) then becomes
SN(fla, va,Wr(t)) + A2Tv(v?, v?,?(t)) = 0
(T(ub, v?,?(i)) --H S(va, v0,Tit(t))) + AlTN('ub, v?,?(!)) = 0.
For force determination, it was also necessary to define functions f(:r, y, 1;) and
g(x, y, p) that explicitly describe the surfaces of A and B in a local neighborhood
of the extremal points. For parametric functions, the analogue of Equation (6--H14)
is
f(Sx(u, v, x(t)), Sy(u, v, TT(t)), mv(t)) = Sz(v', v, TT(t))
and similarly for g and T. As in Section 6.4, the existence of f and g follows from
the implicit function theorem. Partial derivatives of f and g can be expressed in
terms of partial derivatives of 5 and T. If p has k coordinates, then differentiating
173
Equation (A-8) with respect to u, v, and p? through Pk yields
Of OSz			+			Of OSy			OSz
0 Ou			Oy On			Ou
Of OSx			+			Of 05y			OSz
Ov Ov			Oy Ov			Ov
Of OSx			+			Of OS?			+			Of			OSz
OxO??			OyOp?			Op1			Op1
ffifOSr			+			Of OS?			+ Of			OSz
0 0Pk			0Y 0Pk			Opk			=			OPk
Thus, the first partials of f may be expressed in terms of the first partials of 5
by solving a simpler linear system (note that Equation (A--H9) is mostly diagonal).
Second partials of f are obtained by differentiating Equation (A--H9) and solving
another linear system. These derivatives are needed for Equation ((3--H29), along with
the derivatives with respect to time of the extremal points: applying the implicit
function theorem to Equation (A--H4) yields
O(E1,E2)			O(E1,E2)
O(t,Va,Ub,Vb)			and			Vq = --H O(flat,Ul),%`b)			(A--H1o)
=			J			J
(and similarly for ??b and ?)b) where
O(E1,E2)			A-il
0(na, Va, nb, vb)
Given the derivatives of the parametric coordinates, we can calculate Pa and p.b
From Equation (A--Hi), and since Pa = 5(na, va,(t))
___			___			Os?.
Pax = 0?(ua?va?x(t))va + 0?(na?va?([))va + 0p(na?va??(t))?(t) (A--H12)
and similarly for Pay, Pbx and Pby
Appendix B
Termination Conditions of
Lemke's Algoritlim
In this appendix, we will prove that when Leulke's algorithn? fails to find a solution
to Equation (8--H7), it finds a vector z satisfying Equation (13--H1); furthern?ore, for
son?e 1 < i < n (where n is the size of the LCP of Equation (8--H7)), the value Zj is
positive, and the ith contact point of the systen? has dy'ian?ic friction. The proof
is based on a proof by Nlurty[40, section 2.3, pages 86 88] that Len?ke's algorithrn,
when applied to a copositive plus LCP, 1 either finds a solution, or tern?inates with
an unbounded ray.
THEOREM 5. When Lernke `s algorithrn fails to find a sol?ttion to ll?e LCP of
Equation (8--H7), it generates a vector z > 0 with ll?e following two properties.
(1) For all 1 ? i < n, if Zj >0 then (Az)j ? 0
(2) For at least one contact point of the system with dynamic friction,
say llte ith contact point, Zj > 0
1 An LCP is copositive plus if for all x> 0 tile matrix A in the LCP satisfies xTAx> 0 and
(AT + A)x = 0 whenever xTAx = 0.
175
176
PROOF. Since the ordering of contact points is arbitrary, let the contact points
of the system be ordered such that contact points 1 to k are all frictionless, and
contact points k t 1 to n have dynamic friction. Let the matrix A and constant
vector b corresponding to this reordered system be partitioned by
MB
and			b=			b1			(W1)
C D			b2
where M is a k x k matrix, B is a k x --H k) matrix, C is an --H x k matrix,
D is an			--H			x			--H k) matrix, b1 is a k-vector, and b2 is an			--H k)-vector.
Since the first k points are frictionless, the matrix M is PSD, and the vector b1 is
in the range of M. Thus, as established in Chapter 8, the LCP
a?=M!?--Hbi>?O, !N>0			and			!NTa?=o
has a solution.
The proof that a vector satisfying property (1) exists is essentially a restatement
of the termination conditions of Lemke's algorithm as detailed by ?4urty. For
convenience, in this appendix the vector e is taken to be an appropriately sized
vector of ones; that is, Cj = 1 for all i. NVhen Lemke's algorithm fails to find a
solution to the LCP of Equation (13--H1), it finds vectors zh zk, wh, wk and scalars
%h and z0k satisfying
wk --H Azk --H b + z0ke
wh = Azh + z0he
wk > 0			wh > 0			zk > 0			zh > 0 (but zh $ 0),
Thh > 0, z0k > 0
(13-3)
and
(wk + Awh)T(zk + Azh) = 0
for all A > 0.			(13--H4)
177
From Equations (B--H3) and (B--H4), it follows that
Th			(wkT			?Th
(wh) z =0,			) zk=O, (w) z =0			and
From this, it follows that if zhj > 0 then w? = 0. But if zhj > 0, then
0 = w24 = (Azh)i t (Thhe)i
so
(wh)Tzk = 0. (B--H5)
(Azh)i = (Thhe)i = zohei = zoh < 0.
Thus, the vector zh satisfies property (1) of the proof.
In order for zh to satisfy property (2), at least one non-zero element of zh
must correspond to a contact point with dynamic friction. Let ?h, ?k and wit be
partioned as follows:
--H			vh			k			k			h			xh
h --H			uh			? =			,			and			z =			h			.
uk			y
The vectors vh, vk and xh have diu?ension k, while 11h, Uk and yh have dimension
n --H k. To show that property (2) is true for zit, we use proof by contradiction.
Suppose that yh --H 0; that is, suppose that all elements of zh that correspond to a
contact point with dynamic friction are zero. Then zh would have the form
h			h
x --H
h
y			0
Now suppose that ?h > 0. Then since zhis non-zero, for some 1 < j ? k, we
have z3h. > 0. As before, if bh. > 0 then wjh = 0. But then
0 = w,4 = (Azh)j t (%heb = (Azh)j + zoh (B--H10)
(Azh)j = mh ? 0.			(B--H11)
so
178
Since zh satisfies property (1),
1tT
(z ) Azh = ? zt4(Azh)i ? 0;
t=1
(B-12)
but (Azh). --H Thh <0 and > 0 for some 1 <j < k implies that (zh)TAzh <0.
Now, since (zh)T = (xh,o)T,
T			h			hT T
o>(zh) Az =((x ) 0 )
MB			xh			= (xh)TMxh.			(B-13)
CD			0
But this contradicts the fact that M is PSD; tinis z0h $ 0 leads to a contradiction.
Conversely, suppose that z0h = 0. Since (xh)TMxh = (zh)TAzh < 0 and M
is PSD, it must be that (xh)TMxh = 0. However, the fact that M is PSD means
that M is also copositive plus; this in turn implies that (MT + M)xh --H 0. Since
M is symmetric as well, it follows that Mxh = 0.
Since
wk = Azk --H b + %ke,
from the partitioning of A, wk, and zk
uk --H Mxk --H b1 + zoke.
Then since (xh)TM --H (Mxlz)T --H 0T,
(xh)TUk = (xh)T(Mxk --H b1 + ?ke)
h T			k			h T			T k
--H (x ) Mx?--H(x ) b1+(xh) ?e
--H (xh)Tbl + (xh)Tzoke.
Since wk and zh satisfy (w?)(zhj) = 0, it follows that (uk)Txh = 0, so
h TUk = 0 =			h T			h T k
(x )			--H(x) b1 + (x ) zje
(B-14)
(B-15)
179
or equivalently
hTk
--H (xh)Tb1 = --H(x ) z0e.			(B--H16)
Since xh > 0 e> 0 and z0& > 0, and xhj > 0 for some
hT
--H (x ) b1 <0.			(B--H17)
Letting I denote the kxk identity matrix, and (I --HM) the kx2k matrix fon?ed
by catenating I and --HAI, we have
(xh)T(II?M)=(xh)T
--H (xh)TM --H (xh)T > 0 and
?? T
--H(x ) b1 <0.
By Farkas' lemma[40, appendix 1] this implies that the system
(B-18)
(I --HAl)			aN			= --Hb1			and			> 0			(B-19)
fN
is infeasible; that is, no !N and a? satisfying Equation (B--H19) exists. This means
that
(I --HAl)			aN			--H aN --H M!N = --Hb1			and			> 0			(B-20)
or equivalently
aN
a?=M!?--Hbi			and			>0			(B-21)
fN
is infeasible. However, this contradicts the fact that Equation (B--H2) always has a
solution, since it is the LCP of a frictionless configuration. Thus, z0h = 0 leads to a
contradiction.
Since both z0h = 0 and ?h ? 0 lead to contradictions, the original assumption
that yh = 0 must in fact be false. Thus, at least one element z21t of zh corresponding
to a contact point with dynamic friction satisfies z2?? > 0, and zh satisfies both
properties (1) and (2). 0
Appendix c
Iterative Solution Methods for
Static Eriction
In this appendix, a Jacobi-like iterative niethod for solving Equation (14--H1) is
derived. As stated in Section 14.2.3, the iterative n?.thod is applied long enough to
estin?ate the signs of the unknown variables. Gauss-Seidel, Jacobi, and successive
overrelaxation (SOR) iterations have been developed for solving linear con?ple-
n?entarity problen?s[9,33]. The Jacobi-like iteration for solving Equation (14--i)
developed in this appendix can easily be made into either a Gauss-Seidel- or SOR-
like iterative method.
Linear equations
Consider the linear equation
Az + q 0
where A is an n x n mat?x, q is an n-vector, and z is the iii-iknown n-vector being
solved for. Iterative methods generate a sequence of solutions fz(0),z(1),z(2),... 1
that (hopefully) converge to a vector z satisfying Equation (C--H1).
181
182
The Jacobi iteration for computing z(k+l) is the following.
e(k) Az(k)+q
for i = 1 to n
(k+1) _ (k)
done
Note that if z(k) exactly solves Equation (C--H1), then e(k+l) --H Az + q = 0,
and z(k+I) --H z(k). In practice, the iteration is stopped when e(k) becomes suitably
small. The Jacobi iteration is known to converge if A is positive definite (PD); in
this case, the diagonal elements Ajj are all positive.
Linear complementarity
Now suppose that instead of solving Equation (C--H1), we wish to find z such that
Az + q > 0, z> 0			and			zT(Mz + q) = 0.
The Jacobi iteration for solving Equation (C--H2) involves a p?'0jcct?on of z(k), to
maintain the constraint z > 0. The projection of a scalar x so that it is non-
negative is written x+ and is defined as
= max(0,.?).
The Jacobi iteration for solving Equation (C--H2) is simply the Jacobi iteration for
solving Equation (C--H1), with the addition of a projection; that is,
e(k) Az(k)+q
for i = 1 to n
zfk+1) =			--H
done
This is called the projected Jacobi iteration.
183
If z(k) satisfies Equation (C--H2), then > 0 implies
--H (Az(k) + q)i = 0.
Thus 4k+1) = (z?k))+ = z(k). Conversely, if z?k) = 0, then
--H (Az(k) + q)i > 0,
so
= 0
assuming Aji > 0. (As a minimum, iterative methods generally assume the
diagonal elements of A are positive.) Thus, if z(k) satisfies Equation (C 2),
then z(k+l) --H z(k). If A is PD then the projected Jacobi iteration will always
converge[33]
Static friction equations
NVe will use the notation of Equations (14--H6), (14--H7), and (14--H8) in defining the
vector of accelerations, a, and the vector of forces, f, along with the known matrix
A and vector b. (To be consistent, we will write q = --Hb, so that a = Af + q.)
These equations define fri, fyi, ari and ay? for every contact point. If the
configuration of bodies is such that the ith contact point does not have static friction,
then the variables fri, fyi, ?ri, and ay? are left out of Equations (14--H7) and (14--H8);
the corresponding steps for computing fri and fyi in the Jacobi-like iteration are
ignored as well. For simplicity, if fNi, fri and fyi are the pth, qth and rth elements
of f, then the diagonal elements of A corresponding to these indices are denoted
DNi = App, and Dri = Aqq and Dyj = Arr. If the ith contact point has dynamic
friction, it is possible that DNi may be zero or negative if the value of !`i (the
184
coefficient of friction at the contact point) is non-zero, If this happens, then it
is trivial to apply a valid contact impulse to the system by imposing an impulse
at just the ith contact point, converting the dynamic friction to static friction.
Otherwise, the ith contact has static friction and DNi is positive. Thus, without
loss of generality, we may assume DNj > 0 in all cases. Additionally, both Dxi and
Dyj are positive[28]. The Jacobi-like iteration to solve Equation (14--H1) is
a(k) --H			+ q
for i = 1 to n
fNi(k+l) = (f1\tj(k) --H aNi(k+l)IDNi)+
Vf = ((fxi(k))2 + (fyj(k))2)112
= f o+(k)o+(k)
x?			+ fxjaxi
if Vf <JtifNi(k) or p> 0 then
f?j(k+1) --H f j(k) --H axi(k)IDxi
end
fyj(k+l) --H _			--H ?yi (k)/Dyj
Vf? = ((fxi(k+1))2 + (fyj(k+l))2)l/2
if Vf? > !tifNi(k) then
fxi(k+1)			f i(k+l)(flifNi(k)/vit)
fyj(k+l) --H f (k+l)(Jtif?vi(k)/vft)
end
done
The assignment to Lvi(k+l) is the same as in the linear complementarity case.
After this, the iteration proceeds as follows. The value Vf is the magnitude of the
friction force. The vahie p is the dot product of the friction force and the tangential
acceleration. If the magnitude of the friction force, Vf, has not exceeded its upper
bound of JiifNi(k), or if the friction and the acceleration are not opposite pointing
185
(p > 0), then the friction forces are modified. The assignments
fxi(k+1) = fxi(k) --H axi(k)/Dxi
fyj(k+l) --H			_ --H			(k) /Dyj
attempt to modify the friction forces so that the tangential accelerations will be
zero. Following this, the magnitude v11 of the (possibly) modified friction forces
is computed. If the friction force magnitude exceeds its upper bound of IIifNi(k),
the friction forces are suitably scaled (by IiifNi(k)Ivfl) so that their magnitude
exactly attains this upper bound. As in the previous two Jacobi algorithms, when
f(k) satisfies Equation (14--H1), then f(k+1) = f(k) As stated, the purpose of this
algorithm is to examine the signs of the variables after some number of iterations,
as a prediction to solve the QP of Section 14.2.3.
Appendix D
Simulation Examples
This appendix contains in?agcs fron? simniations produced using the methods
described in this thesis. Figure D.1 sho?vs a sinnilation of a jack rolling down
stairs. Figure D.2 shows a simulation of superquadric dice falling through a lattice.
Figure D.3 shows a simulation of a collapsing structure with friction,
187
Figure D.1: Jack falling down stairs.
189
191
Figure D.1 (Continued)
191
Figure D.1 (Continued)
Figure D.2: Falling superquadric dice.
193
Figure D.2: Falling superquadric dice.
193
Figure D.2 (Continued)
195
Figure D.2 (Continued)
195
Figure D.3: Collapsing structure.
197
Figure D.3: Collapsing structure.
197
199
Figure D.3 (Continued),,,
199
Figure D.3 (Continued),,,
Bibliography
[1] V.1. Arnold. Mathematical MeThods of Classical Mechanics. Springer-Verlag,
1978.
[2] D. Baraff and A. NVitkin. Dynan?ic sin?ilation of non-penetrating flexible
bodies. Submitted to Computer Graphics (Proc. SIGGRAPH), 1992.
[3] M.R. Bracli. Rigid body collisions. Jonnial of Applied Mechanics, pages 164--H
170,1989.
[4]
[5]
C. Cai and B. Roth. On the spatial motion of a rigid body with point contact.
In International Conference on Robotics and Automation, pages 686--H695. IEEE,
1987.
S.A. Cameron and R.I&. Culley. Determining the minimum trauslational dis-
tance between two convex polyhedra. In International Conference on Robotics
and Automation, pages 591-596. IEEE, 1986.
[6] J. Canny. Collision detection for moving polyhedra. IEEE Transactions on
Pattern Analysis and Machine Intelligence, 8(2), 1986.
[7] R.?V. Cottle. On a problen? in linear inequalities. Journal of the London
Mall?ematical Society, 43:378--H384,1968.
[8] J.F. Cremer. An Architecture for General Purpose Physical System &imulation.
PliD thesis, Cornell University, April 1989.
[9] C.?V. Cryer. The solution of a q?iadratic programming problem using system-
atic overrelaxation. SlAM Journal on Control, 9(3):385--H392, 1971.
P.A. Cundall. Formulation of a three-dimensional distinct element model Part
I. A scheme to represent contacts in a system composed of many polyhedral
blocks. International Journal of Rock Mechanics, Mineral Science and Geome-
chanics, 25,1988.
[10]
201
202
[11] J.E. Dennis, Jr. and R.B. Sd?nabel. Numerical Methods for Unconstrained
Optimization and Nonlinear Equations. Prentice Hall, Inc., 1983.
[12] M.A. Erdmann. On motion planning with uncertainty. Master's thesis,
Massachusetts Institute of Technology, 1984.
[13] R. Featherstone. Robot Dynamics Algorithms. Nluwer, 1987.
[14] k. Fletcher. Practical Methods of Optimization, volume 2. John ?Viley & Sons,
1981.
[15] G.E. Forsythe, M.A. Malcolm, and C.B. Moler. Computer Mell?ods for Math-
ematical Computations. Prentice Hall, Inc., 1977.
[16] M.R. Garey and D.S. Johnson. Computers and i??fractability. Freeman, 1979.
[17] E.G. Gilbert and C.P. Foo. Computing the distance between smooth objects
in three dimensional space. In International Conference on Robotics and
Automation, pages 158--H163. IEEE, 1989.
[18]
[19]
[20]
E.G. Gilbert and S.M. Hong. A new algorithm for detecting the collision
of moving objects. In Piternational Conference on Robotics and Automation,
pages 8--H13. IEEE, 1989.
E.G. Gilbert, D.W. Johnson, and 5.5. Neerthi. A fast procedure for computing
the distance between complex objects in three space. In P?ternational Confer-
ence on Robotics and Automation, pages 1883 1889. IEEE, 1987.
B.J. Gilmore and R.J. Cipra. Sin?ulation of planar dynamic mechanical sys-
tems with d?anging topologies: Part 1--Hcharacterization and prediction of the
kinematic constraint changes; Part 2--Himplementation strategy and simulation
results for example dynamic systems. In Advances in Design Automation, pages
369--H388. ASME, 1987.
[21] H. Goldstein. Classical Mechanics. Addison-NVesley, Reading, 1983.
[22]
J.P. Gourret, N.M. Thalmaun, and D. Thalmaun. Simulation of object and
human skin deformations in a grasping task. In Computer Graphics (Proc.
SICGRAPH), volume 23, pages 21-30. ACM, July 1989.
[23] NV. Ingleton. A problem in linear inequalities. Proceedings of the London
Mathematical Society, 16(3) :519--H536,1966.
[24] J.I3. Neller. Impact with friction. Journal of Applied Mechanics, 53(1):1--H4,
1986.
203
[25] W. Kilmister and J.E. Reeve. Rational Mechanics. Lougman's, 1966.
[26]
B. Lafleur, N. Magnenat-Thalmann, and D. Thalmann. Cloth animation with
self-collision detection. In T.L. Nunii, editor, Modeling in Computer Graphics.
Springer-Verlag, 1991
[27] C.E. Lemke. Bimatrix equilibrium points and mathematical programming.
Management Science, 11:681--H689,1965.
[28] P. L5tstedt. Coulomb friction in two-dimensional rigid body systems.
Zeitschrift fur Angewandte Mathematik un Mechanik, 61:605--H615,1981.
[29] P. L5tstedt. Mechanical systen? of rigid bodies subject to unilateral con-
straints. SlAM Journal of Applied Mathematics, 42(2):281--H296, 1982.
[30]
[31]
P. L5tstedt. Numerical simulation of time-dependent contact friction problems
in rigid body mechanics. SlAM Journal of Scientific Statistical Computing,
5(2):370--H393, 1984.
T. Lozano-Pe'rez. Automatic planning of manipulator transfer movements.
IEEE Transactions on Systems, Man and Cybernetics, SMC-11(10):681--H698,
1981.
[32] T. Lozano-Pe'rez. Spatial planning: A configuration space approach. JEEE
Transactions on Computers, C-32(2):108--120, 1983.
[33] O.L. Niangasarian. Solution of symmetric linear complementarity problems
by iterative methods. Journal of Optimization Theory and Applications,
22(4):465--H485, 1977.
[34]
[35]
[36]
M.T. Mason and Y. Wang. Modeling impact dynamics for robotic operations.
In International Conference on Robotics and Automation, pages 678 685. IEEE,
1987.
M.T. Mason and Y. Wang. On the inconsistency of rigid-body frictional planar
mechanics. In International Conference on Robotics and Automation, pages
524--H528. IEEE, 1988.
M. MeKenna and D. Zeltzer. Dynamic simulation of autonomous legged
locomotion. In Computer Graphics (Proc. SIGGRAPH), volume 24, pages
29--H38. ACM, August 1990.
[37] N. Megiddo. Linear time algorithms for linear programming in 1?? and related
problems. SlAM Journal of Computing, 12(4):759--H776, 1983.
204
[38] P.M. Moore. A flexible object animation system. Master's thesis, University of
California, Santa Cruz, 1987.
[39] P.M. Moore and J. Wilhelms. Collision detection and reponse for computer
animation. In Computer Graphics (Proc. SIGGRAPH), volume 22, pages 289--H
298. ACM, August 1988.
[40] N.G. Murty. Linear Complementarity, Linear and Nonlinear Programming.
lleldermann Verlag, 1988.
[41] Ju.I. Ne?rnark and N.A. Fufaev. Dynamics of NonhoThnomic Systems. Ameri-
can Mathematical Society, 1972.
[42] R.S. Palmer. Computational CompThxity of Motion and Stability of Polygons.
PhD thesis, Cornell University, January 1987.
[43] F.P. Preparata and M.I. Shamos. Computational Geometry. Springer-Verlag,
New York, 1985.
[44] NV.H. Press, B.P. Flannery, S.A. Teukolsky, and N\?.T. Vetterling. Numerneal
Recipes. Cambridge University Press, 1986.
[45] H. Rubin and P. Ungar. Motion under a strong constraining force. Communz-
cations on Pure and Applied Malliematics, 10:65--H87,1957.
[46] R. Sedgewick. Algornthms. Addison-NVesley, 1983.
[47] L.F. Shampine and M.I?. Gordon. Computer Solution of Ordinary Differential
Equations: The Initial Value Problem. Freeman, 1975.
[48] I&. Shoemake. Animating rotation with quaternion curves. In Computer
Graphics (Proc. SIGGRAPH), volume 19, pages 245 254. ACM, July 1985.
[49] A.E. Taylor and R.M. Mann. Advanced Calculus. John NViley & Sons, Inc.,
1983.
[50] D. Terzopoulos and K. Fleisclier. Deformable models. Visual Computer, 4:306--H
331,1988.
[51]
D. Terzopoulos, J.C. Platt, and A.H. Barr. Elastically deformable models.
In Computer Graphics (Proc. SIGGRAPH), volume 21, pages 205--H214. ACM,
August 1987
[52] S.A. Vavasis. Quadratic programming is in NP. Technical Report TR 90-1099,
Department of Computer Science, Cornell University, 1990.
205
[53] B. Von Herzen, A. Barr, and II. Zatz. Geometric collisions for time-dependent
parametric surfaces. In Gornpi?ter Graphics (Proc. SIGGRAPH), volume 24,
pages 39--H48. ACM, August 1990.
[54] Y. Wang and M.T. Mason. Two dimensional rigid body collisions with friction.
Journal of Applied Mechanics, (to appear).
[55] J. Wilhelms. Toward automatic motion control. IEEE Computer Graphics and
Applications, 7(4):11--H22, 1987
[56] A. NVitkin and W. Welch. Fast animation and control of nonrigid structures.
In Computer Graphics (Proc. SIGGRAPH), volume 24, pages 243--H252. ACNI,
August 1990.
205
[53] B. Von Herzen, A. Barr, and II. Zatz. Geometric collisions for time-dependent
parametric surfaces. In Gornpi?ter Graphics (Proc. SIGGRAPH), volume 24,
pages 39--H48. ACM, August 1990.
[54] Y. Wang and M.T. Mason. Two dimensional rigid body collisions with friction.
Journal of Applied Mechanics, (to appear).
[55] J. Wilhelms. Toward automatic motion control. IEEE Computer Graphics and
Applications, 7(4):11--H22, 1987
[56] A. NVitkin and W. Welch. Fast animation and control of nonrigid structures.
In Computer Graphics (Proc. SIGGRAPH), volume 24, pages 243--H252. ACNI,
August 1990.
