Noname manuscript No.
(will be inserted by the editor)

The Maximum Dissipation Principle in Rigid-Body Dynamics with
Purely Inelastic Impacts
Tobias Preclik1 , Sebastian Eibl1 , Ulrich Rüde1,2
1

arXiv:1706.00221v1 [cs.CE] 1 Jun 2017

2

Lehrstuhl für Informatik 10 (Systemsimulation), Friedrich-Alexander Universität Erlangen-Nürnberg, Cauerstr. 11, 91058
Erlangen, Germany
CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, Cedex 01, France

The date of receipt and acceptance will be inserted by the editor

Abstract Formulating a consistent theory for rigidbody dynamics with impacts is an intricate problem.
Twenty years ago Stewart published the first consistent
theory with purely inelastic impacts and an impulsive
friction model analogous to Coulomb friction. In this
paper we demonstrate that the consistent impact model
can exhibit multiple solutions with a varying degree of
dissipation even in the single-contact case. Replacing the
impulsive friction model based on Coulomb friction by
a model based on the maximum dissipation principle resolves the non-uniqueness in the single-contact impact
problem. The paper constructs the alternative impact
model and presents integral equations describing rigidbody dynamics with a non-impulsive and non-compliant
contact model and an associated purely inelastic impact
model maximizing dissipation. An analytic solution is
derived for the single-contact impact problem. The models are then embedded into a time-stepping scheme. The
macroscopic behaviour is compared to Coulomb friction
in a large-scale granular flow problem.

Key words Impulse (physics) – Coulomb friction –
collisions (physics) – rigid body dynamics – contact dynamics – impact dynamics – measure differential inclusions – complementarity problems

Correspondence to: Tobias Preclik
tobias.preclik@fau.de

1 Introduction
Simulating mechanical systems on computers requires
a model for describing the dynamics of the mechanical
parts. Models that can describe the deformation of the
mechanical parts require a high number of parameters
to describe the deformation. If the core of the mechanical parts can be assumed to not deform under the considered loads, the parameters describing the state of a
mechanical part can be reduced to that of a rigid body:
An invariant shape with an associated mass and principal moments of inertia as well as the spatial orientation,
position, linear and angular velocity of the shape. The
interaction of multiple such mechanical parts must be
described by another model determining the dynamics
of the mechanical parts in contact, where the contact
model usually allows compliance in a localized contact
region. If mechanical parts collide, the contact dynamics typically occur on a time-scale that is significantly
smaller than that of the motions between successive collisions. Resolving each such collision micro-dynamics in
computer simulations can become computationally expensive. Alternatively, the relation between the pre- and
post-collision state variables can be described by an impact model [29]. By condensing the impact dynamics to
a single point in time, only the response of the relative
contact velocities has to be specified. In order to instantaneously turn a colliding state into a non-colliding state
a contact reaction impulse must be applied, so that the
necessary discontinuities in the velocities can be effected.
The impact model must then be combined with a noncompliant contact model in order to determine the contact reaction forces and contact reaction impulses. Noncompliant contact models alone cannot resolve collisions.
For non-compliant contact models with Coulomb friction even non-colliding contact situations exist (shocks)
where an impulse becomes necessary to resolve the contact [26]. These paradoxical situations were first published by Painlevé [18]. The system including the noncompliant contact and impact models are often mathe-

2

matically described in terms of measure differential inclusions (MDI) [16].
The rigid body simplification combined with the
impact simplification considerably reduce the computational burden in simulations with many mechanical
parts. Applications range from robotics [10, 17], virtual
reality [22], physics-based animation [5] to granular
matter simulations [30]. In particular, granular dynamics can require a very large number of particles and
are insufficiently understood to date with and without
an interstitial fluid phase [15]. Optimizing mechanical
devices like powder mixers [7] or grinding mills [8, 14]
is of economical importance. Powder mixing is important in detergent, cosmetic, food and pharmaceutical
manufacturing, to name just a few applications. Understanding granular matter is also of crucial importance
for safety reasons: Assessing the stability of slopes is
important to prevent rock slides, land slides and snow
avalanches and getting the particle distribution right in
pebble-bed nuclear reactors is important to guarantee
safe and performant operation [30].
The construction of an impact model that in combination with a non-compliant contact model leads to a
consistent theory for rigid body dynamics is non-trivial.
Ideally, the solution of an impact model for a collision
is a limit point of a sequence of solutions of the collision based on a compliant contact model with increasing stiffness. The sequence of solutions is uniquely determined and the increasing stiffness decreases the duration of the collision towards an instantaneous event.
The solution of the collision based on a compliant contact model corresponds to the integral of the contact
reaction forces over the collision duration. Using such
an approach Stronge constructs an energetically consistent restitution hypothesis in [28], Mirtich solves rigidbody dynamics with impacts for virtual reality applications, where permanent contacts are treated as sequences
of collisions [12, 13]. And lately, Jia and Wang showed
in [11] how contact reaction impulses can be computed
from the limit of a contact model with linear normal stiffness and Coulomb friction for general collisions (central
or eccentric, direct or oblique) in three dimensions. The
authors established a condition that, if met, guarantees
solution existence. Whether solutions exist unconditionally and whether the impact model in combination with
a non-compliant contact model leads to a consistent theory remain open problems.
However, other impact models exist that possess
solutions unconditionally. Stewart showed in [24] for a
time-stepping scheme based on [1], that it converges to
a solution of an MDI as the time-step size decreases.
Recently, Gavrea et al. extended the result to systems
including joints in [6]. The MDI describes rigid-body
dynamics with Coulomb friction and purely inelastic
impacts (collisions and shocks), where the frictional impulses are required to directly oppose the post-impulse
relative contact velocities in the tangential planes

Tobias Preclik et al.

thus imitating Coulomb’s friction law in the case of
impulses. A consequence of this is that Stewart proved
that solutions exist for the MDI and thus resolved
paradoxical configurations in rigid-body dynamics with
non-compliant contacts and Coulomb friction, where
apparently no solutions exist even though no collisions
are present. Stewart made no attempt to show uniqueness of solutions. In fact in section 4.1 we present an
example demonstrating the existence of multiple solutions of a numerically constructed single-contact impact
problem. The non-uniqueness is in this single-contact
case directly related to the choice of the frictional
impact model. In this paper we construct an alternative
frictional impact model having a unique solution in
the single-contact case. The friction model is based on
the maximum dissipation principle [26] and takes into
account the coupling between the normal component
and the tangential components of the contact reactions.
The non-uniqueness due to redundant constraints in
the multi-contact case remains unaffected [19] as well
as non-uniqueness in the non-compliant contact model.
In section 2 we present integral equations describing rigid-body dynamics with impact and friction, where
the Coulomb friction model on the impulsive reactions
is replaced by a friction model based on the maximum
dissipation principle. In section 3 the model is embedded into an impulse-velocity time-stepping scheme for
numerically integrating multi-contact problems and an
analytic solution of the single-contact problem is established. Subsequently, section 4 presents results for singlecontact problems and the macro-scale behaviour of the
friction model in the simulation of a large-scale granular flow problem. The paper summarizes the results and
concludes in section 5.

2 Continuous System
Each particle i is associated with a co-rotating body
frame. The origin of the body frame corresponds to the
center of mass of the particle. Let xi be the position
function of the body frame in the inertial frame. The
position function is non-smooth when impulses act. The
derivative with respect to time is the discontinuous linear velocity function vi with left- and right-limits vi−
and vi+ . The orientation of the body frame in the inertial
frame can be represented by a unit quaternion. Instead
of mixing vector and quaternion algebra a quaternion
qw + qx i + qy j + qz k ∈ H describing the orientation of
the body frame of particle i at time t is represented as
a vector ϕi (t) = (qw , qx , qy , qz )T ∈ R4 . The orientation
function ϕi is non-smooth and the left- and right limit
of the derivative at time t is then [4]


−qx −qy −qz
1  qw qz −qy 
−/+
 ω −/+ (t) = Qii (ϕi (t))ω −/+ (t),
ϕ̇i (t) = 
i
2 −qz qw qx  i
qy −qx qw

The Maximum Dissipation Principle in Rigid-Body Dynamics with Purely Inelastic Impacts
−/+

where ωi (t) ∈ R3 is the angular velocity of the particle
before⁄after applying impulses at time t. We also introduce
the quaternion matrix function Qii for abbreviating the
notation.
The mass of the particle is denoted by mi and is
invariant with respect to time. The inertia tensor of the
particle in the inertial frame changes with respect to
the orientation of the body frame. It can be expressed
in terms of the constant inertia tensor in the body
frame Ii,0 . For time t the inertia tensor in the inertial
frame is given by
T

Iii (ϕi (t)) = R(ϕi (t))Ii,0 R(ϕi (t)) ,
where R(ϕi (t)) is the rotation matrix corresponding to
the orientation ϕi (t). The rotation matrix changes the
basis from the particle’s body frame to the inertial frame.
Choosing the body frame such that the axes match the
principal axes of the particle, the body frame inertia
tensor Ii,0 can be enforced to be diagonal.
−
In a system with νb particles, let x(t) (ϕ(t), v /+ (t),
−/+
and ω (t)) be the vertical concatenation of all particles’ positions (orientations, linear velocities, and angular velocities) at time t:
x(t) = vertcat xi (t) ∈ R3νb .
i=1..νb

Let
Q(ϕ(t)) = diag Qii (ϕi (t)),
i=1..νb

I(ϕ(t)) = diag Iii (ϕi (t)),
i=1..νb

and let M(ϕ(t)) be the block-diagonal mass matrix containing diagi=1..νb mi E3 in the upper-left quadrant and
I(ϕ(t)) in the lower-right quadrant, where E3 is the
3 × 3 identity matrix. Then given initial conditions at
time t0 , the state of the system at time t is described by
the integral equations

 
 Z t

x(t)
x(t0 )
v − (tp )
=
+
dtp ,
ϕ(t)
ϕ(t0 )
Q(ϕ(tp ))ω − (tp )
t0
! 
 X


−
v /+ (t)
v − (t0 )
−1 ∆p(tq )
=
+
M(ϕ(t
))
−
q
ω − (t0 )
∆L(tq )
ω /+ (t)
tq ∈Tq
t0 ≤ tq </≤ t

3

is to be understood as the vertical concatenation of all
single-particle cross-products:
ω − (t) × I(ϕ(t))ω − (t) = vertcat ωi− (t) × Iii (ϕi (t))ωi− (t).
i=1..νb

The appearance of the term stems from the fact that the
torque function τi corresponds to the time-derivative of
the angular momentum function Li (for non-impulsive
points in time), which in turn is the product of the timevarying inertia tensor and the angular velocity. Hence,
τi (tp ) =

dLi (t)
dt

t=tp

=

d
dt Iii (ϕi (t))ωi (t)

= Iii (ϕi (tp ))ω̇i (tp ) +

d
dt Iii (ϕi (t))

t=tp

t=tp

ωi (tp )

= Iii (ϕi (tp ))ω̇i (tp ) + ωi (tp ) × Iii (ϕi (tp ))ωi (tp ).
The forces, torques, linear impulses, and angular impulses at time t include components from non-impulsive
contact reactions λ(t) ∈ R3νc and impulsive contact reactions Λ(t) ∈ R3νc , where νc is the number of contacts
in the particle system. Each contact j involves a pair of
particles (i1 (j), i2 (j)). By convention let contact reactions act positively on the first particle i1 (j) and negatively on the second particle i2 (j). Each contact j is
also associated with a contact frame. Let the first axis
of the contact frame correspond to the contact normal
nj (t) pointing from particle i2 (j) towards particle i1 (j)
by convention, and let orthonormal vectors tj (t), and
oj (t) complete the contact frame. Let x̂j (t) denote the
position of the contact frame in the inertial frame. Then,
subsuming all forces and torques on particle i, which are
not due to contact reactions, as external forces fi,ext (t)
and external torques τi,ext (t), the equations
X
X
fi (t) = fi,ext (t) +
λj (t) −
λj (t),

τi (t) = τi,ext (t)

+

j=1..νc
i1 (j)=i

j=1..νc
i2 (j)=i

X

(x̂j (t) − xi (t)) × λj (t)

j=1..νc
i1 (j)=i

−

X

(x̂j (t) − xi (t)) × λj (t),

j=1..νc
i2 (j)=i

where Tq is the set containing all points in time tq , where
impulses are present, that is (linear) impulse ∆p(tq ) ∈
R3νb or angular impulse ∆L(tq ) ∈ R3νb is non-zero:

Tq = t ∆p(t) 6= 0 ∨ ∆L(t) 6= 0 .

define a wrench matrix function W relating the wrenches
to the contact reactions. This relation extends to impulsive reactions and linear and angular impulses:

 

f(t)
fext (t)
=
+ W(t)λ(t),
τ(t)
τext (t)

 

∆p(t)
∆pext (t)
+ W(t)Λ(t).
=
∆L(t)
∆Lext (t)

The terms f(t) ∈ R3νb and τ(t) ∈ R3νb are the forces
and torques acting on the particles. Note that the inverse of the mass matrix always exists, since it is symmetric positive-definite (SPD) - a property which it inherits from its diagonal blocks. The cross-product term

The impulsive and non-impulsive contact reactions are
then given implicitly as solutions of contact constraints.
The contact constraints are usually non-linear and underdetermined depending on the specific contact model
employed.

Z

t

+

M(ϕ(tp ))
t0

−1




f(tp )
dtp ,
−
−
τ(tp ) − ω (tp ) × I(ϕ(tp ))ω (tp )

4

Tobias Preclik et al.

The formulation of the contact constraints requires
the rigorous introduction of the contact position function x̂, the contact normal function n, the signed contact distance function ξ and the relative contact velocity
function δv. The latter is straightforward and for a contact j given by
−/+

δvj

−/+

−/+

−/+

−/+

(t) = vi1 (j) (t) + ωi1 (j) (t) × (x̂j (t) − xi1 (j) (t))
− vi2 (j) (t) − ωi2 (j) (t) × (x̂j (t) − xi2 (j) (t)).

It can be shown that
δv

−/+

T

(t) = W(t)

!
−
v /+ (t)
.
−
ω /+ (t)

The definition of the other three functions are difficult to state in sufficient generality. We confine ourselves
here to definitions that are at least well-defined for spherical particles. Let Si (t) be the set of points in the inertial
frame defining the shape of particle i at time t, and let
fS : R3 → R be the signed distance function associated
with the shape S. The signed distance function shall be
negative in the interior of the shape. Then, let
x̂j (t) ∈

arg min
fSi

2 (j)

(t) (y)≤0

fSi1 (j) (t) (y)

be the contact point between the pair of particles (i1 (j),
i2 (j)). If the boundary of the shape is sufficiently smooth
and the overlap sufficiently small, the contact position is
uniquely determined and the gradient of the signed distance function exists. Then the contact normal is given
by
nj (t) = ∇fSi2 (j) (t) (x̂j (t)).
The signed contact distance function is then simply
ξj (t) = fSi1 (j) (t) (x̂j (t)).
These specific definitions of the contact functions limit
the number of contacts νc to the number of particle
pairs ν2b (νb −1). To simplify the description of the contact
constraints, subscript n denotes the projection of a vector to the contact normal (e.g. λj,n (t) := nj (t)T λj (t) ∈
R) and subscript to denotes the vector of projections of
a vector to the contact tangential and contact orthogonal (e.g. λj,to (t) := (tj (t)T λj (t), oj (t)T λj (t))T ∈ R2 ).
Then the contact constraints for an inelastic contact with Coulomb friction are listed in Fig. 1. The
constraints can be classified into non-penetration constraints and Coulomb friction constraints. Both classes
can be subdivided into impulsive and non-impulsive
constraints. Impulsive constraints determine impulsive
contact reactions Λ(t) and non-impulsive constraints
determine non-impulsive contact reactions λ(t). Some
constraints are understood to be enabled only if a
precondition holds. For instance the restitution hypothesis should only constrain the solution if the contact
is closed (ξj (t) = 0). This precondition is indicated

by an arrow. The arrow originates from a constraint,
which enables the constraint if the precondition becomes active. In the case of the restitution hypothesis
the contact needs to close first and the arrow thus
originates from the impulsive Signorini condition. The
Signorini condition ensures that contact reactions are
non-negative (non-adhesive) if the contact is closed
and are zero if the contact is open. This relation
is expressed by the complementarity condition ⊥
and the corresponding inequalities. A similar chain of
non-penetration constraints exists for the non-impulsive
contact reactions. However, in that chain the contact
reaction in the worst case can only be determined after
the constraint on the acceleration level became enabled.
The impulsive and non-impulsive Coulomb friction
constraints require the contact reaction to reside within
a friction cone. The coefficient of friction µj determines
the aperture 2 tan−1 µj of the cone. The cone is aligned
along the contact normal. The friction cone condition
limits the Euclidean norm of the frictional reaction by an
upper bound proportional to the contact reaction in normal direction. The direction of the frictional reaction is
required to oppose the relative tangential contact velocity in the case of a sliding (dynamic) contact and its Euclidean norm must be at its limit. This is expressed in the
velocity-level equation. In the case of a sticking (static)
contact, the velocity-level equation is universally valid.
However, the zero slip enables the acceleration-level constraint. Then, the direction of the frictional reaction is
required to oppose the relative tangential contact acceleration.
The work performed by the frictional contact reac 
tion force of contact j for a non-impulsive time span t, t
is
Z
t

+
λj,to (t)T δvj,to
(t) dt,

t

where the Coulomb friction force performs no work to
the extent possible as expressed in the acceleration-level
constraint. However, if sliding is inevitable, the Coulomb
friction force maximizes dissipation by minimizing the
integrand through the velocity-level constraint, which
requires the friction force to directly oppose the relative
contact velocity. The velocity-level Coulomb constraint
can be formulated equivalently using the maximum dissipation principle as pointed out by Stewart in his review paper on friction and impact in rigid-body dynamics [26] at least if the normal reaction is considered to
be given [25]:
λj,to (t) ∈

arg min
kyk2 ≤µj λj,n (t)

+
y T δvj,to
(t),

(1)

where the objective function corresponds to the
(negated) rate of energy dissipation. Since the relative contact velocity at time t is independent of the
contact reaction at time t, the objective function is a
linear function of the frictional contact reaction at a
non-impulsive point in time t 6∈ Tq .

The Maximum Dissipation Principle in Rigid-Body Dynamics with Purely Inelastic Impacts

non-compliant

contact model

impact

model

Non-penetration constraints

Fig. 1:

Coulomb friction constraints

ξj (t) ≥ 0 ⊥ λj,n (t) ≥ 0

ξj (t) = 0
ξ˙j+ (t) = 0

λj,to (t)

ξ˙j+ (t) ≥ 0 ⊥ λj,n (t) ≥ 0

+
(t)
δvj,to

ξ¨j+ (t)

˙ +
δv
j,to (t)

≥ 0 ⊥ λj,n (t) ≥ 0

ξj (t) ≥ 0 ⊥ Λj,n (t) ≥ 0

ξj (t) = 0

ξ˙j+ (t) ≥ 0 ⊥ Λj,n (t) ≥ 0

Signorini condition,

−/+

(t) = U (t) + T
1
= U (t) +
2

−/+

+
δvj,to
(t)

restitution hypothesis,

(t)
−

v /+ (t)
−
ω /+ (t)

!T

!
−
v /+ (t)
M(ϕ(t))
,
−
ω /+ (t)

where the potential energy is not affected by impulses.
Insertion leads to the expression for the post-impulse
system energy
1
E + (t) = E − (t) + Λ(t)T W(t)T M(ϕ(t))−1 W(t)Λ(t)
2
 − 

!
v (t)
T
T
−1 ∆pext (t)
+ Λ(t) W(t)
+ M(ϕ(t))
ω − (t)
∆Lext (t)

T  − 

!
1
∆pext (t)
v (t)
−1 ∆pext (t)
+
M(ϕ(t))
+
∆Lext (t)
ω − (t)
∆Lext (t)
2
=

2

1
Λ(t)T A(t)Λ(t) − Λ(t)T b(t) + c1 (t)
2

in terms of the pre-impulse system energy and the impulsive contact reactions. Let A(t) be the Delassus operator and let b(t) condense the terms depending linearly

2

≤ µj λj,n (t)

+
(t)
λj,to (t) = −µj λj,n (t)δvj,to

λj,to (t) =
2

Λj,to (t)

When formulating the friction constraint on the
impulsive contact reactions, the situation changes
subtly but drastically: The impulsive contact reactions
now influence the post-impulse relative contact velocity.
The drastic consequence of this is that the maximum
dissipation principle and the Coulomb friction model
for that matter as it is formulated for non-impulsive
contact reactions in Eq. (1) cannot be transferred to
impulsive contact reactions without in-depth modifications if the property of maximizing the energy
dissipation is to be preserved. This statement stands in
contrast to common practice [3, 9, 26]. In particular the
+
term Λj,to (t)T δvj,to
(t) is a quadratic function of the
impulsive contact reactions and it does not reflect the
energy dissipated. Hence, impulsive frictional reactions
directly opposing the relative contact velocity in the
tangential plane also do not dissipate as much energy
as allowed by the friction cone condition in general.
−
The system energy E /+ (t) is the sum of the potential
−
energy U (t) and the kinetic energy T /+ (t):
E

5

2

2

˙ +
−µj λj,n (t)δv
j,to (t)

+
δvj,to
(t)

2

=0

≤ µj Λj,n (t)

+
(t)
Λj,to (t) = −µj Λj,n (t)δvj,to

friction cone condition,

frictional reaction opposes slip

on the impulsive contact reactions and let c1 (t) condense
the constant terms. A contact reaction Λj (t) complying
with the maximum dissipation principle should minimize
E + (t). Restating E + (t) in terms of the j-th contact reaction and assuming all other contact reactions to be
constant results in
Ej+ (Λj (t)) :=

1
Λj (t)T Ajj (t)Λj (t)
2



− Λj (t)T bj (t) − Ajj (t)Λj (t) + c2 (t),
where Ajj (t) corresponds to the j-th 3×3 diagonal block
of the Delassus operator and where j selects all columns
(elements) except for column j (element j). The diagonal
block can be determined to be [12, 20]
−1
Ajj (t) = (m−1
a + mb )E3

− (x̂j (t) − xa (t))× Iaa (ϕa (t))−1 (x̂j (t) − xa (t))× (2)
− (x̂j (t) − xb (t))× Ibb (ϕb (t))−1 (x̂j (t) − xb (t))× ,

where a = i1 (j) and b = i2 (j).
Then the impulsive contact reaction complying with
the maximum dissipation principle is
Λj (t) ∈

ξj (t) = 0

arg min

Ej+ (Λ),

(3)

kΛto k2 ≤ µj Λn

(3a)

ξj (t) ≥ 0 ⊥ Λn ≥ 0

(3b)

ξ̇j+ (t) ≥ 0 ⊥ Λn ≥ 0

(3c)

nj (t)T Ajj (t)Λ ≥ 0

(3d)

where Eq. (3a) corresponds to the friction cone condition, Eq. (3b) corresponds to the Signorini condition,
Eq. (3c) corresponds to the purely inelastic restitution
hypothesis, and Eq. (3d) is an additional constraint requiring that the contact reaction is not increasing the
contact pressure. The last constraint guarantees uniqueness for a single contact. It excludes non-zero solutions
if the contact opens by itself. The objective function is
a quadratic function of the contact reactions and it is
strictly convex since Ajj (t) is SPD.

6

Tobias Preclik et al.

For open contacts (ξj (t) > 0) the Signorini constraint
and the friction cone constraint restrict the feasible set
to the reaction Λj (t) = 0. The restitution hypothesis is
disabled and the pressure constraint is fulfilled.
If the contact is closed (ξj (t) = 0) the Signorini
condition reduces to Λj,n (t) ≥ 0 and the restitution
hypothesis is enabled. In order to determine whether
the restitution hypothesis is active (ξ˙j+ (t) = 0) or inactive (ξ˙j+ (t) > 0), the dependence of ξ˙j+ (t) on the impulsive contact reactions Λ(t) must become explicit. At
least for spherical particles the time-derivative of the
post-impact signed distance function can be expressed
in terms of the relative contact velocities:
ξ˙j+ (t) = δv +
j,n (t)
= nj (t)T (Ajj (t)Λj (t) + Ajj (t)Λj (t) − bj (t)).

(4)

Contacts fulfilling the property nj (t)T (Ajj (t)Λj (t)−
bj (t)) < 0, that is contacts where a penetration is imminent if no impulsive contact reaction acts, are termed
colliding in the following. Pre-impulse velocities, external impulses and impulsive reactions from other simultaneously colliding contacts determine if the contact is
in a colliding state.
Contacts fulfilling the property nj (t)T (Ajj (t)Λj (t)−
bj (t)) > 0, that is contacts where separation is imminent if no impulsive contact reaction acts, are termed
separating. For separating closed contacts, the reaction
Λj (t) = 0 fulfills all constraints and thus the restitution
hypothesis is inactive. The pressure constraint ensures
that no non-zero solutions exist.
For colliding closed contacts, the restitution hypothesis must be active, restricting the feasible set to the plane
of maximum compression defined by Eq. (4). Combined
with the friction cone condition, the feasible set forms a
conic section. The normal of the plane of maximum compression is Ajj (t)nj (t). Since Ajj (t) is SPD and since
the contact is colliding, the conic section is guaranteed
to be non-empty. Since the conic sections are non-empty
convex sets and since the objective function is strictly
convex, the optimization problem has a unique global
minimum. The pressure condition is fulfilled since it is
fulfilled for any point on the plane of maximum compression in the colliding case.
The unconstrained global minimum Λ0 is given by
∇Ej+ (Λ0 ) = Ajj (t)Λ0 − (bj (t) − Ajj (t)Λj (t)) = 0.
If the contact is colliding and closed and if Λ0 fulfills
all constraints, then Λj (t) = Λ0 and the post-impulse
relative contact velocity in the tangential plane is zero
corresponding to a sticking (static) post-impulse contact state. If Λ0 does not fulfill all constraints the postimpulse relative contact velocity in the tangential plane
is non-zero and the post-impulse contact state thus sliding (dynamic). The solution Λj (t) of the contact problem then resides on the boundary of the conic section.

The friction model for contact impulses from Eq. (3)
thus maximizes dissipation by minimizing the (kinetic)
energy analogously to the velocity-level constraint of
the Coulomb model for non-impulsive contact reactions.
This is in contrast to the alleged extension of the
Coulomb model to impulsive reactions, where the dissipation is not (sufficiently) maximized. At the same time
the contact impulses fulfill the friction cone condition
as in the Coulomb model and they fulfill the inelastic
restitution hypothesis. The friction model is lazy in the
sense that separating solutions (no contact reaction)
are preferred over any other solutions, which is guaranteed through the pressure condition. Furthermore,
static solutions (no work performed) are preferred over
dynamic solutions analogously to the acceleration-level
constraint of the Coulomb model.
3 Numerical Methods
3.1 Multi-Contact Problems
Integrating the equations describing the rigid-body dynamics can be approached in at least two different ways:
The event-driven approach aims to predict the next impulsive point in time tq ∈ Tq given an initial state at
time t0 . Then the simulation is integrated until tq assumingthe impulsive
contact reactions Λ(t) to be zero

for t ∈ t0 , tq . At time tq an impact problem given by
 +
  −



v (tq )
v (tq )
−1 ∆p(tq )
=
+ M(ϕ(tq ))
ω + (tq )
ω − (tq )
∆L(tq )
and Eq. (3) has to be solved. Then, the integration can
be restarted having a state fulfilling all constraints at
hand. The difficulty of this approach lies in the problem
of predicting the next impulsive point in time, which
is a priori unknown. If rigid bodies follow ballistic trajectories impact times can be predicted accurately and
efficient simulation codes exist [2, 12]. For the more general case, where e.g. impulsive points in time can not
only stem from collisions but also stem from self-locking
sliding frictional contacts [23], we know of no efficient
method to accurately predict the next impulsive point
in time. Furthermore, simulation codes necessarily stall
in situations, where the collision frequency increases unboundedly like in cases where a bouncing ball comes to
rest on a plane.
The second category of approaches for integrating
rigid-body dynamics are the time-stepping methods.
These methods proceed in time steps δt > 0 independent of the impulsive points in time. The methods do
not distinguish between non-impulsive and impulsive
contact reactions but implicitly solve for integrals of
the contact reactions λ̃. The integrals of the contact
reactions are then constrained to fulfill contact conditions at selected points in time. In the following an
impulse-velocity time-stepping scheme is used, which is

The Maximum Dissipation Principle in Rigid-Body Dynamics with Purely Inelastic Impacts

described in detail in [20]. It is similar to the schemes
by Anitescu [1], Tasora [31] and Stewart [27]. The
equations of motion are integrated using a discretization
similar to a semi-implicit Euler method, where positions
are integrated implicitly and velocities are integrated
explicitly:
x0 (λ̃) = x + δtv 0 (λ̃),
ϕ0 (λ̃) = ϕ + δtQ(ϕ)ω0 (λ̃) / ϕ + δtQ(ϕ)ω0 (λ̃) ,
2

!  
!
0
v (λ̃)
v
f˜ext
−1 
,
=
+ M(ϕ)
Wλ̃ +
ω
τ̃ext − δtω × I(ϕ)ω
ω 0 (λ̃)




where primes identify state variables at time t + δt in
contrast to state variables at time t omitting the prime.
The relative contact velocities are
!
v 0 (λ̃)
0
T
.
δv (λ̃) = W
ω 0 (λ̃)
The discrete non-penetration constraint for a contact j
is
ξj
+ δv 0j,n (λ̃) ≥ 0 ⊥ λ̃j,n ≥ 0,
δt
which if the gap is closed exactly (ξj = 0) states
that the relative contact velocity at the end of the
time step and in the direction of the contact normal
must be non-approaching and complementary to the
non-adhesive contact reaction in normal direction. The
ξ
term δtj acts as an error reduction term if penetrations
are present (ξj < 0). In that case it can be scaled
down when needed using an error reduction parameter
ε ∈ [0, 1[ to avoid introducing an excessive amount of
ξ
ξ
energy (ε min(0, δtj ) + max(0, δtj )). If a positive gap is
present (ξj > 0), the term ensures that contact reaction
remains zero if the contact would not close within the
time step. The friction cone condition
λ̃j,to

2

≤ µj λ̃j,n

is adopted as it stands. Instead of requiring the impulsive contact reactions to fulfill the conventional Coulomb
friction constraints as in
0
δvj,to
(λ̃)

2

0
λ̃j,to = −µj λ̃j,n δvj,to
(λ̃),

the maximum dissipation principle from Eq. (3) is used,
since the discretized relative contact velocity in the tangential plane depends on the contact reactions as before.
Therefore, the energy term is discretized, leading to
1
E 0 (λ̃) = U + T 0 (λ̃) = λ̃T WT M(ϕ)−1 Wλ̃
2

!
 
˜
v
fext

+ λ̃T WT 
+ M(ϕ)−1
ω
τ̃ext − δtω × I(ϕ)ω
!
 
1
v
T
−1
+z
+ M(ϕ) z + E
ω
2
=

1 T
λ̃ Aλ̃ − λ̃T b + c1 ,
2

7


T
T
where z = f˜ext
, (τ̃ext − δtω × I(ϕ)ω)T . For a single
contact j the energy term reduces to
Ej0 (λ̃j ) :=

1 T
λ̃ Ajj λ̃j − λ̃T
j (bj − Ajj λ̃j ) + c2
2 j

given all other contact reactions λ̃j . The relative contact
velocity in terms of λ̃j and λ̃j is
δvj0 (λ̃j ) := Ajj λ̃j − (bj − Ajj λ̃j ).
Hence, the contact constraint complying with the maximum dissipation principle is
λ̃j ∈

arg min
kλto k2 ≤µj λn
ξj
δt

Ej0 (λ).

(5)

+δv 0j,n (λ)≥0 ⊥ λn ≥0
nT
j Ajj λ≥0

Then a solution λ̃ of the multi-contact time-step problem satisfies Eq. (5) for all contacts j = 1 . . νc . The problem can be solved using a non-linear block Gauss-Seidel
or variants thereof as demonstrated in [20, 21]. The solution algorithm proposed there reduces the multi-contact
problem to the problem of solving a sequence of singlecontact problems. Hence, in the next section an analytic solution for the single-contact impact problem from
Eq. (5) (and Eq. (3) alike) is derived.
3.2 Single-Contact Problems
The single-contact problem for impulsive contact reactions complying with the maximum dissipation principle
from Eq. (3) constrains solutions to zero if contacts are
open (ξj (t) > 0). If contacts are closed, Eq. (3) has the
same structure as the discrete single-contact time-step
problem from Eq. (5):
x∗ ∈

arg min
kxto k2 ≤µxn
aT
n x−bn ≥0 ⊥ xn ≥0
aT
n x≥0

1 T
x Ax − xT b,
2

where µ ∈ R≥0 , A is SPD and




Ann Ant Ano
Ann
A =  Ant Att Ato  ∈ R3×3 , an =  Ant  ∈ R3 ,
Ano Ato Aoo
Ano
 
 
 
xn
bn
xt
x =  xt  ∈ R3 , xto =
∈ R2 , b =  b t  ∈ R3 .
xo
xo
bo
If the contact is separating (bn < 0) the solution is constrained to x∗ = 0 since the pressure condition (aT
nx ≥
0) prevents any non-zero solutions. If the contact is nonseparating (bn ≥ 0) the pressure condition is redundant
and the constraint set is formed by the intersection of the
plane of maximum compression and the friction cone.
This conic section is non-empty but can take on the

8

Tobias Preclik et al.

shape of an ellipse, parabola or hyperbola. The cases
where bn = 0 are special in the sense that the plane
of maximum compression includes the origin and thus
the conic section degenerates to a point, a ray or a degenerate hyperbola. These cases will be discussed after
the cases where bn > 0. The cases where bn > 0 and
µ = 0 are also special since the conic section degenerT
ates to a single (non-zero) point x∗ = (A−1
nn bn , 0, 0) . If
bn > 0 and µ > 0 the conic section is non-degenerate.
The unconstrained minimum of the objective function
is x0 = A−1 b. If it fulfills the friction cone condition it
is contained in the constraint set and it thus solves the
contact problem x∗ = x0 .
If the unconstrained minimum of the objective function does not fulfill the friction cone condition, the solution must be located on the boundary of the conic
section minimizing the objective function. Eliminating
the normal component using the equation for the plane
of maximum compression and switching to polar coordinates leads to

T
x = A−1
nn (bn − Ant xt − Ano xo ), xt , xo

T
= A−1
,
nn (bn − Ant r cos α − Ano r sin α), r cos α, r sin α

T
(a) The ellipsoidal contour lines of 12 xT
to Âxto − xto b̂ are solid
gray except the zero contour line which is solid black. The
conic section is a hyperbola and drawn as a solid blue curve.
The unconstrained minimum is marked by a black cross. The
red cross is the maximally dissipative solution along the hyperbola.

where r ∈ R≥0 and α ∈ R. The friction cone condition
then becomes
kxto k2 ≤ µA−1
nn (bn − Ant xt − Ano xo )
↔
↔

r(Ann + µAnt cos α + µAno sin α) ≤ µbn
q
r (Ann + µ A2nt + A2no cos(α − atan2(Ano , Ant ))) ≤ µbn .
|
{z
}
=fr (α)

For angles α where fr (α) ≤ 0, the inequality poses no
additional restrictions on the non-negative coordinate r.
For angles α where fr (α) > 0, the inequality defines
an upper bound on the coordinates r. In both cases the
component r satisfies
0≤r≤

µbn
.
max(0, fr (α))

Let I + 2πN be the set of angles for which fr (α) > 0
holds. I can be determined to be
(
p
]−∆α, ∆α[ + α0 if Ann ≤ µ A2nt + A2no
I=
, (6)
[0, 2π[
else
nn
where ∆α = acos √−A
2

µ

Let

Ant +A2no

and α0 = atan2(Ano , Ant ).



cos α
γ(α) = r(α)
sin α
describe the curve along the boundary of the conic secn
. Then the contact problem retion, where r(α) = fµb
r (α)
duces to
=fobj (α)

z
}|
{
1
T
T
α∗ ∈ 2πN + arg min γ(α) Âγ(α) − γ(α) b̂
2
α∈I
r∗ = r(α∗ ),

(b) The objective function fobj is plotted in solid blue. The
non-feasible region is indicated by gray stripes. The poles are
marked by dotted lines. The global minimum is marked by a
red cross.

Fig. 2: An exemplary contact problem, where A, b and
µ are such that the conic section is a hyperbola and
the objective function fobj has two local minima, but
only one global minimum. Only the global minimum is
dissipative.
where
2
Att − A−1
Ato − A−1
nn Ant
nn Ant Ano
Â =
−1
2
Ato − Ann Ant Ano Aoo − A−1
nn Ano


bt − A−1
nn Ant bn
b̂ =
bo − A−1
nn Ano bn




and

result from eliminating the normal component. The objective function fobj is univariate and 2π-periodic but no
longer strictly convex nor quadratic. Fig. 2 illustrates the
optimization problem for an exemplary (dynamic) contact. Fig. 2a plots the quadratic objective function for
points xto on the plane of maximum compression and
overlays the conic section. Fig. 2b plots the corresponding non-linear objective function fobj in comparison. An

The Maximum Dissipation Principle in Rigid-Body Dynamics with Purely Inelastic Impacts

iterative approach for solving this constrained minimization problem that is guaranteed to converge linearly is
derived in [20]. In the following an analytic approach is
presented.
Let
t(α) =

d
dβ r(β) sin β β=α
d
dβ r(β) cos β β=α



!
=

dr(β)
dβ
β=α
 dr(β)
dβ
β=α

cos α − r(α) sin α
sin α + r(α) cos α




be the unit vector tangential to the curve, where
q
r(α)2
dr(β)
=
A2nt + A2no sin(α − atan2(Ano , Ant )).
dβ
β=α
bn
Angles minimizing fobj must satisfy
T

t(α) (Âγ(α) − b̂) = 0.

(7)

(8)

with constants c0 , c1 , c2 , ϕ1 , ϕ2 ∈ R as specified in
appendix A. The trigonometric equation can be
transformed into a quartic equation by substituting y = eiα = cos α + i sin α. After solving the
quartic equation for yi , the corresponding angles
αi = atan2(Im(yi ), Re(yi )) have to be checked for
validity. Angles αi 6∈ 2πN + I are invalid as well as
angles not satisfying Eq. (7) or Eq. (8) for that matter.
Among all other candidates, the one with minimum
objective function value amounts to the 2π-periodic
solution α∗ . Finally,
 −1

Ann (bn − Ant r(α∗ ) cos α∗ − Ano r(α∗ ) sin α∗ )
x∗ =
.
γ(α∗ )
Back to the case, where bn = 0 and µ > 0: If the
conic
p section corresponds to a degenerate ellipse (Ann >
µ A2nt + A2no ), then the solution is exactly this point.
The conic sections
p that correspond to a degenerate parabola (Ann = µ A2nt + A2no ) are equivalent to the ray
( 
)

Ant
≥0
r
r∈R
.
Ano
Degenerate hyperbolas correspond to the conical combination of two rays (along the asymptotes)
(


r1

the friction cone. If that is not the case, then the objective function needs to be minimized along one of the
rays:
(

α0 − ∆α if (Ano , −Ant )T x0,to ≥ 0
,
α0 + ∆α else
(9)

T 


T
cos α∗
cos α∗ −1 cos α∗
r∗ = max(0, (
Â
)
b̂).
sin α∗
sin α∗
sin α∗

α∗ ∈ 2πN +

The structogram in Fig. 3 summarizes the analytic solutions of the various cases involved in solving the singlecontact problem.
4 Numerical Results

Insertion, trigonometric identity transformations and
µ2 b2
multiplication by r(α)n3 leads to a trigonometric equation
in the form of
c0 + c1 cos(α + ϕ1 ) + c2 cos(2α + ϕ2 ) = 0,

9

)



cos(α0 − ∆α)
cos(α0 + ∆α)
≥0
+ r2
r1 , r2 ∈ R
,
sin(α0 − ∆α)
sin(α0 + ∆α)

where α0 and ∆α as in Eq. (6). Thus degenerate parabolas are a limiting case of degenerate hyperbolas with
∆α = π. Minimizing the energy over degenerate parabolas and hyperbolas can be implemented by checking,
whether the unconstrained minimum x0 is contained in

4.1 Academic Single-Contact Problem
As an artificially constructed example consider a rigid
body composed of six mass points connected by massless
rods. The mass points are located symmetrically on the
axes (±dx ex , ±dy ey , ±dz ez ) and both mass points on
each axis concentrate the same mass (mx , my , mz ). The
center of mass thus coincides with the origin and the
inertia tensor is given by
I = 2 diag(my d2y + mz d2z , mx d2x + mz d2z , mx d2x + my d2y ).
Let the mass point in the positive x direction contact a
plane with normal n. Let t = n
 × ex and o = n × t span
the tangential plane. Let Q = n t o denote the contact
frame. Then
−1 ×
A = QT ((2(mx + my + mz ))−1 E − d2x e×
ex )Q.
xI

A rendering of this problem using dx = 10,
dy = dz = 1, mx = 0.03, my = 65, mz = 50,
µ = 3.7, n̂ = (0.25, 0.36, −0.9)T , n = n̂/kn̂k2 and
b = QT (0.06, 0.7, 0.23)T is shown in Fig. 4. The
Coulomb solutions can be calculated for example using
the polynomial root-finding approach described in
[3]. The contact problem has three dynamic Coulomb
solutions xi (i ∈ {1, 2, 3}) directly opposing the
post-impulse relative contact velocity in the tangential
plane and a unique maximally dissipative solution x∗ .
Approximate numerical values are listed in Tab. 1 for
reference. Fig. 5 plots the contour lines of the objective
function for points in the plane of maximum compression and the boundary of the constraint set in solid
blue. The red cross marks the maximally dissipative
solution along the ellipse. The red circles mark the three
Coulomb solutions.
4.2 Paradox Single-Contact Problem
Let a contact exist at time t0 between
a slender rod

of length l ∈ R>0 and a half-space x ∈ R3 nT x ≤ 0
with normal n = ey . Let t = ex and o = ez complete the

10

Tobias Preclik et al.

true
bn < 0

x∗ = 0
true
µ=0

T
x∗ = (A−1
nn bn , 0, 0)

true
P1

x∗ = x0
true

true

x∗ = 0

P2
α∗ , r∗ from Eq. (9)

bn = 0

false
Let S be the sol. set of Eq. (8).
α∗ ∈ 2πN + arg min fobj (α)
α∈I∩S

false

false

false

false r∗ = r(α∗ )

P1 = (x0,n ≥ 0 ∧kx0,to k2 ≤ µx0,n )

x∗,n = A−1
nn (bn − Ant r∗ cos α∗
−A
 no r∗ sin α∗ )

x∗,n


x∗ = r∗ cos α∗ 
r∗ sin α∗

P2 = (Ann > µ

p

A2nt + A2no )

Fig. 3: The structogram describing the analytic solution of the single-contact problem.

Fig. 4: Illustration of the academic single-contact problem. An object constituted out of six point masses symbolized as spheres connected by massless rods collides
with an inclined plane. The rods are aligned with the
coordinate axis and the point masses are arranged symmetrically on each axis.
≈ xT
i


0.12 −0.44 −0.065


0.29 −0.89 −0.62


1.56 −0.99 −5.7

≈ (Axi − b)T


0 0.40 0.060


0 0.082 0.057


0 0.0061 0.035

T
≈ 12 xT
i Axi − xi b

≈ xT
∗


1.6 −1.1 −5.8

≈ (Ax∗ − b)T


0 −0.057 0.034

T
≈ 12 xT
∗ Ax∗ − x∗ b

Fig. 5: Illustration of the academic contact problem including contour lines of the objective function, constraint
set and solutions.

−0.252
−0.393
−0.631

−0.634

Table 1: Approximate Coulomb solutions and maximally dissipative solution for the numerically constructed single-contact problem including relative contact velocities in the contact frame and objective function values.

Fig. 6: Illustration of the paradox single-contact problem. A slender rod is grazing along a plane.

The Maximum Dissipation Principle in Rigid-Body Dynamics with Purely Inelastic Impacts



contact frame, where ex ey ez = E3 . Let 0 < θ(t0 ) =
θ0 < π2 denote the initial angle between the plane and
the rod. Let the slender rod correspond to rigid body a
and let the half-space correspond to rigid body b. The
half-space is stationary such that vb (t) = ωb (t) = 0,
mb = ∞ and Ib = ∞E3 . The rod shall be centered
at the origin and aligned along the x axis of the body
frame. Then, with uniformly distributed mass ma , its
2
T
body-frame inertia tensor is Ia,0 = ml
12 diag(∞, 1, 1) .
3
T
The rod is located in the x-y plane x ∈ R ez x = 0
of the inertial frame with its center of mass at xa (t0 ) =
l
T
2 (cos θ0 , sin θ0 , 0) , such that the contact is closed initially (ξ(t0 ) = 0). Let ωa− (t0 ) = 0 and va− (t0 ) = −ex v0
with v0 ∈ R>0 , such that initially the contact neither
separates nor collides (ξ˙− (t0 ) = δv −
n (t0 ) = 0). Let the
contact position function x̂ track the lower tip of the
slender rod. Let gravity g = −ey g0 with g0 ∈ R>0 act.
A visualization of this setup is displayed in Fig. 6. Under
these conditions the contact problem is essentially planar
and corresponds to the paradox configuration published
by Painlevé [18, 26].
If t0 is a non-impulsive point in time (t0 6∈ Tq ),
−
ξ˙− (t0 ) = ξ˙+ (t0 ) and since ξ(t0 ) = ξ˙ /+ (t0 ) = 0, the
acceleration-level non-penetration constraint is enabled
from Fig. 1, where
−

−
˙ /+ (t0 ) and
ξ¨ /+ (t0 ) = δv
n

˙
δv

−/+

−

−

(t0 ) = v̇a /+ (t0 ) + ω̇a /+ (t0 ) × (x̂(t0 ) − xa (t0 ))
−/+

−

+ ωa (t0 ) × (ωa /+ (t0 ) × (x̂(t0 ) − xa (t0 )))
!
−/+
v̇
(t
)
a
0
= W(t0 )T
−/+
ω̇a (t0 )
= W(t0 )T M(ϕ(t0 ))−1 W(t0 )λ(t0 ) + g.
Since the contact is sliding at time t0 , the frictional contact reaction force is known to be
 
−/+
δvto (t0 )
1
= µλn (t0 )
.
λto (t0 ) = −µλn (t0 )
−/+
0
δvto (t0 )
2

Consequently,

  
µ

 λn (t0 )
λ(t0 ) = n t o µλn (t0 ) =  1  λn (t0 ).
0
0
Since gravity acts, the acceleration-level non-penetration
constraint is compelled to be active. Thus,
 
µ
−/+
T
˙
 1  λn (t0 ) − g0 = 0, λn (t0 ) ≥ 0,
δv
n (t0 ) = n A(t0 )
0
where A(t0 ) is given by Eq. (2):



1 + 3 sin2 θ0 −3 sin θ0 cos θ0
0


2
0
A(t0 ) = m−1
.
a −3 sin θ0 cos θ0 1 + 3 cos θ0
0
0
1 + 3(cos2 θ0 − sin2 θ0 )2

11

The coefficient is then

 
µ
nT A(t0 )  1  = m−1
a (1 + 3 cos θ0 (cos θ0 − µ sin θ0 )).
0
The equation is not solvable for non-negative λn (t0 ) if
the coefficient is negative. The sign of the coefficient depends on µ and θ0 . Thus, for a given angle 0 < θ0 < π2
the contact problem at hand has a non-impulsive solution if
1
cos θ0 + 3 cos
θ0
µ < µ∗ (θ0 ) =
,
sin θ0
where the lowest bound on µ is 34 at the angle
1
−1
− 35 . If this condition is not met, the assump2 cos
tion that t0 is a non-impulsive point in time is wrong.
Instead an impact problem has to be solved beforehand.
The impact model with purely inelastic impacts and
Coulomb-like friction as presented in Fig. 1 exhibits multiple solutions. The zero solution Λ(t0 ) = 0 is perfectly
valid, since no collision is taking place (ξ(t0 ) = ξ˙− (t0 ) =
0). However, it clearly does not lead to a post-impact
state with a non-impulsive solution.
Any other solution must be located on the plane of
˙ 0 ) = 0) and in the friction
maximum compression (ξ(t
cone. The impulse necessary to obtain a post-impulse
sticking contact state is determined by the vector equation δv + (t0 ) = A(t0 )Λ0 + va− (t0 ) = 0:

5 + 3 cos(2θ0 ) 3 sin(2θ0 )
ma 
Λ0 =
 3 sin(2θ0 ) 5 − 3 cos(2θ0 )
8
0
0
=

 
v0
 
 0 =
16
0
5+3 cos(4θ0 )
0
0

ma v0
(5 + 3 cos(2θ0 ), 3 sin(2θ0 ), 0)T .
8

It is easily verified, that Λ0 resides within the friction
cone for any 0 < θ0 < π2 and any coefficient of friction
requiring an impulsive solution (µ ≥ µ∗ (θ0 )):
Λ0,to

2

ma v0
(5 + 3 cos(2θ0 ))
8
m a v0
= µ∗ (θ0 )
3 sin(2θ0 )
8
m a v0
3 sin(2θ0 ) = µΛ0,n .
≤µ
8

=

Thus Λ(t0 ) = Λ0 is a second solution of the impact
problem effecting a slip-stick transition. In fact, there
also can be an infinite number of sliding solutions. For
example in the edge case where µ = µ∗ (θ0 ), all convex
combinations of the 0 solution and the Λ0 solution are
also solutions, where the slip directly opposes the frictional impulse.
When replacing the impact model with the Coulomblike friction by the impact model complying with the
maximum dissipation principle, the structogram in
Fig. 3 easily identifies the unique impulsive solution: If
µ ≥ µ∗ (θ0 ), then Λ(t0 ) = Λ0 . It is the same impulsive
solution as the sticking solution in the impact model

12

Tobias Preclik et al.

avg. number of particles
solid volume fraction
avg. number of contacts
avg. number of contacts per particle
max. penetration (µm)

Fig. 7: Top-down view of the simulation area. The particle channel is confined by solid walls in y and z direction.
A constant particle inflow at the left side is artificially
created. Within the channel an obstacle is present to
disturb the particle flow.
with the Coulomb-like friction. Conversely, if a nonimpulsive solution exists (µ < µ∗ (θ0 )), then Λ0 is not
contained in the friction cone and predicate P2 in the
structogram then always selects Λ(t0 ) = 0:
Ann = 1 + 3 cos2 θ0 = µ∗ (θ0 )3 sin θ0 cos θ0
q
> µ3 sin θ0 cos θ0 = µ A2nt + A2no .
Thus the impact problem results in the zero solution if
a non-impulsive solution exists as expected. The impact
model with Coulomb-like friction behaves the same way
in this respect: If µ < µ∗ (θ0 ) then the plane of maximum
compression intersects the friction cone only at 0.
The impulsive contact reaction Λ(t0 ) = Λ0 leads
in both impact models to a sticking post-impulse contact state, where non-impulsive contact reactions exist
for the subsequent contact problem. However, alternate
solutions of the impact model with Coulomb-like friction
do not necessarily lead to subsequent contact problems
with non-impulsive solutions.

4.3 Macro-Scale Behaviour
The macro-scale behaviour of the proposed maximum
dissipation friction model is compared to the Coulomb
friction model by numerical simulations of fast granular
channel flows. All simulations are carried out with the pe
module of the waLBerla multi physics framework freely
available at walberla.net. The algorithms are based on
time stepping methods presented in Sec. 3. The parallel
implementation is described in [20, 21]. For the computations the compute resources of the Regionales Rechenzentrum Erlangen (RRZE) are used. The Emmy cluster
comprises 560 compute nodes, each equipped with two
Xeon 2660v2 proccessors (10 cores, 2-way SMT) clocked
at 2.2 GHz and 64 GiB of RAM.
The simulation domain is a rectangular channel confined by solid walls in the y and z direction. Its dimensions are 15 cm × 6.5 cm × 2.0 cm. The channel is filled
with monodisperse spherical particles with a diameter

Coulomb

Max. Dissipation

5

2.47 · 105
57.2 %
3.20 · 105
2.59
46.5

2.47 · 10
57.2 %
3.20 · 105
2.59
39.1

Table 2: Key figures collected for both simulations. Samples were collected every 10 ms during the 10 s total simulation time.

of 0.47 mm and a density of 2.65 g/cm3 . These are arranged on a regular rectangular grid with a spacing of
1 mm. All particles are given an initial velocity of 1 m/s
in positive x direction. A random perturbation velocity
in y and z direction is applied with each component varying between −0.5 m/s and 0.5 m/s. To produce a steady
inflow a moving plane with infinite mass is added at the
left end of the channel (x = 0 cm). This plane moves
with 1 m/s pushing the particles into the channel. After the plane has moved a distance equal to two times
the radius of the particles the position of the plane is
reset and a new layer of particles is generated. This process is repeated throughout the whole simulation generating an inflow rate of roughly 1.4 · 106 particles per
second. The other end of the channel remains open and
the particles leaving the channel are deleted. All the
particles are also influenced by a gravitational acceleration of 9.81 m/s2 in negative z direction. In the region
[3 cm, 5 cm] × [3 cm, 5 cm] × [0 cm, 2 cm] a stationary object obstructing the channel is introduced. The obstacle
is composed of a cylinder with radius 1 cm located at
(4 cm, 4 cm, 1 cm) with its axis aligned along the z axis.
The right half of the cylinder is replaced by a box. A
cross-section of this setup is sketched in Fig. 7.
For the computation the domain is decomposed into
16 × 10 × 4 evenly sized rectangular subdomains. Each
subdomain is handled by one process totaling to 640 processes which are distributed onto 32 nodes of the Emmy
cluster. The time-step length is chosen to be 10 µs. This
guarantees a stable simulation. The total simulation time
is 10 s resulting in 106 time steps. After these time steps
all initial perturbations are resolved. The implementation of the collision resolution solver described in [20] is
used. Both friction models use a coefficient of friction
equal to 0.4.
Some general information about the simulation is
collected every 1000 time steps and averaged over the
complete simulation. Within the range of accuracy both
friction models result in almost identical values which
are summarized in Tab. 2. The average number of particles in the simulation domain is 2.47 · 105 which leads
to a solid volume fraction of 57.2 %. During every time
step an average of 3.20 · 105 contacts which equals 2.59
contacts per particle are resolved. The maximum penetration depth between two particles was 39.1 µm (8.31 %

The Maximum Dissipation Principle in Rigid-Body Dynamics with Purely Inelastic Impacts

0.0

13

Coulomb Friction

Coulomb Friction

Maximum Dissipation

Maximum Dissipation

0.5

1.0

0.0

0.5

1.0

1.5

2.0

velocity (m/s)

angular velocity (rad/ms)

Fig. 8: Velocity profile of the particles throughout the
channel. Around the obstacle the velocity of the particles
increases due to the narrowing of the channel. After the
obstacle collisions as well as friction slow down the fast
particles again.

Fig. 9: Angular velocity profile of the particles throughout the channel. One can see that rotations do not play
a major role in this setup as friction damps any angular
velocity rapidly. A certain amount of rotating particles
can be found at the inflow and near the stationary object
in the middle of the channel.

of the particle radius) in the Coulomb friction case and
46.5 µm (9.89 % of the particle radius) in the maximum
dissipation case.
To validate the accordance of both friction models
the last frame of each simulation is used and the particle configuration is analyzed. The particles are sorted
into 150 × 65 equally sized cells in the x-y plane ignoring
the z axis and the velocity (see Fig. 8) as well as the
angular velocity (see Fig. 9) is averaged over all particles within one cell. Both simulations show a very similar
velocity profile throughout the whole channel. The particles right after the inflow have very high velocities but
are damped rapidly. After a short almost homogeneous
region the particles gain speed again as the channel narrows around the obstacle. The particle velocities peak
at 1.42 m/s (Coulomb Friction) and 1.39 m/s (Maximum
Dissipation) respectively. The angular velocity comparison shows a similar picture. After a short initial region
the angular velocity gets damped heavily. At the boundary of the stationary obstacle and the walls, the angular
velocity stays high due to friction. Both friction models
result visually in a similar behaviour.
As a last check the radial distribution function
of the particles near the exit is analyzed to obtain
insight into the spatial arrangement. The region
[12 cm, 14 cm] × [0 cm, 6.5 cm] × [0 cm, 2 cm] is chosen.

All the inter-particle distances are calculated and
summed up in a histogram. The histogram is then
normalized by the total number of particles and the
expected
amount of particles within each individual bin
Rx
- x01 4πρr2 dr, with ρ being the overall particle density
and x0 , x1 being the boundaries of the corresponding
histogram bin. For both friction models the histogram
can be found in Fig. 10. They match almost exactly.
One can see a strong peak at a distance equal to the
particle diameter which indicates a dense packing. The
rest of the histogram does not reveal additional favored
distances suggesting an amorphous packing.
The numerical experiment described above shows
no significant difference between the two friction
models. The newly developed maximum dissipation
friction model thus reproduces the expected macro scale
behaviour.

5 Summary
In this paper, an alternative purely inelastic frictional
impact model is presented. In that impact model
the contact reaction impulses consistently maximize
dissipation leading to unique contact reactions for

14

REFERENCES
Coulomb Friction

12

10

g(r)

8

6

4

2

0
1.0

1.5

2.0

2.5

3.0

r/d

Maximum Dissipation

12

10

g(r)

8

6

4

2

and a weighted non-linear block Jacobi, where each
block corresponds to a single-contact problem. The
paper presents an analytic solution of the single-contact
problem that can act as a subsystem solver in the
non-linear block relaxation method. The analytic solution involves transforming the single-contact problem
into a quartic equation with complex coefficients. Back
transformation and filtering of the solutions of the
quartic equation leads to the unique contact reaction
maximizing dissipation. Additionally, the analytic
solution can be used to resolve two-particle collisions in
event-driven integrations of rigid-body dynamics.
The presented impact model is missing a restitution
hypothesis for partly elastic impacts. Even though applying Poisson’s hypothesis is straightforward, the extension of the impact model by an energetically consistent restitution hypothesis is not. Also proving or falsifying that the constructed time-stepping scheme converges to a solution of the corresponding integral equations remains an open problem. Besides addressing the
non-uniqueness in the impact model, the approach might
also be transferable to remove non-uniqueness in the
non-compliant contact model when non-penetration and
friction constraints on the acceleration-level are enabled.

0
1.0

1.5

2.0

2.5

3.0

r/d

Fig. 10: Radial distribution function of the particles near
the outflow. The radial distance r is given in units of the
particle diameter d = 0.47 mm. The dotted line g(r) = 1
represents the value for a completely amorphous material.

single-contact impact problems. An academic singlecontact impact problem is analyzed, demonstrating
the potential non-uniqueness of the impulsive contact
reactions if Coulomb-like friction constraints act instead. Furthermore, a paradox single-contact impact
problem was analyzed. It was found that the impact
model maximizing dissipation results in the solution
producing a slip-stick transition, whereas the impact
model with Coulomb-like friction in addition allows
the zero solution and possibly even an infinite number
of dynamic solutions, which not necessarily result in
configurations with subsequent non-impulsive solutions.
The paper also shows how the impact model based on
the maximum dissipation principle can be embedded
in an impulse-velocity time-stepping scheme. A numerical experiment is conducted to demonstrate that
changing from an impact model with Coulomb-like
friction to the maximally dissipative impact model in
the time-stepping scheme has a negligible influence on
the macroscopic behaviour of a granular channel flow
past an obstacle. The multi-contact problem is solved
using a blend between a non-linear block Gauss-Seidel

Acknowledgment
The authors would like to acknowledge the support
through the Cluster of Excellence Engineering of
Advanced Materials (EAM).
A Appendix
c0 = 23 µAnn (Ano bt − Ant bo )
c1 = (µ4 ((A2no (A2tt + A2to ) − 2Ant Ano Ato (Att + Aoo ) +
A2nt (A2to + A2oo ))b2n + 2(Ant bo − Ano bt )(A2no Ato + Ant
Ano (Att − Aoo ) − A2nt Ato )bn + (A2nt + A2no )(Ant bo −
Ano bt )2 ) + 2µ2 Ann ((2Ant Ano Ato − A2no Att − A2nt Aoo )
b2n +Ann (Ano (Att bo −Ato bt )−Ant (Ato bo −Aoo bt ))bn +
Ann (Ant bo − Ano bt )2 ) + A2nn ((A2nt + A2no )b2n − 2Ann
1
(Ant bt + Ano bo )bn + A2nn (b2t + b2o ))) /2
c2 = µAnn (( 14 (A2tt + A2oo ) + A2to − 12 Att Aoo )b2n + (Ano ( 12
(Att − Aoo )bo − Ato bt ) − Ant (Ato bo + 12 (Att − Aoo )bt ))
1
bn + 41 ((A2nt + A2no )b2o + (A2nt + A2no )b2t )) /2
2
2
ϕ1 = atan2(µ ((Ano Ato − Ant Aoo )bn − Ano bt + Ant Ano
bo ) − A2nn bt + Ann Ant bn , µ2 ((Ant Ato − Ano Att )bn +
Ant Ano bt − A2nt bo ) − A2nn bo + Ann Ano bn )
ϕ2 = atan2((Att − Aoo )bn − Ant bt + Ano bo , 2Ato bn − Ano
bt − Ant bo )
References
[1]

M. Anitescu and F. Potra. “Formulating Dynamic
Multi-Rigid-Body Contact Problems with Friction
as Solvable Linear Complementarity Problems”.
In: Nonlinear Dynamics 14.3 (1997), pp. 231–247.

REFERENCES

[2]

[3]

[4]
[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

M. N. Bannerman, R. Sargant, and L. Lue. “DynamO: A free O(N) general event-driven molecular
dynamics simulator”. In: Journal of computational
chemistry 32.15 (2011), pp. 3329–3338.
O. Bonnefon and G. Daviet. Quartic Formulation
of Coulomb 3D Frictional Contact. Technical Report RT-0400. INRIA, 2011.
J. Diebel. Representing Attitude: Euler Angles,
Unit Quaternions, and Rotation Vectors. 2006.
K. Erleben. “Stable, Robust, and Versatile Multibody Dynamics Animation”. PhD thesis. University of Copenhagen, 2004.
B. I. Gavrea, M. Anitescu, and F. A. Potra.
“Convergence of a class of semi-implicit timestepping schemes for nonsmooth rigid multibody
dynamics”. In: SIAM Journal on Optimization
19.2 (2008), pp. 969–1001.
A. Hassanpour et al. “Analysis of Particle Motion
in a Paddle Mixer using Discrete Element Method
(DEM)”. In: Powder Technology 206.1 (2011),
pp. 189–194.
C. Jayasundara et al. “CFD–DEM Modelling of
Particle Flow in IsaMills – Comparison Between
Simulations and PEPT Measurements”. In: Minerals Engineering 24.3 (2011), pp. 181–187.
M. Jean. “The non-smooth contact dynamics
method”. In: Computer methods in applied mechanics and engineering 177.3 (1999), pp. 235–
257.
Y.-B. Jia. “Three-dimensional impact: energybased modeling of tangential compliance”. In: The
International Journal of Robotics Research 32.1
(2013), pp. 56–83.
Y.-B. Jia and F. Wang. “Analysis and Computation of Two Body Impact in Three Dimensions”.
In: Journal of Computational and Nonlinear Dynamics (2016).
B. Mirtich. “Impulse-based Dynamic Simulation
of Rigid Body Systems”. PhD thesis. University of
California, 1996.
B. Mirtich and J. Canny. “Impulse-based simulation of rigid bodies”. In: Proceedings of the 1995
symposium on Interactive 3D graphics. ACM.
1995, 181–ff.
B. Mishra and R. Rajamani. “The Discrete Element Method for the Simulation of Ball Mills”.
In: Applied Mathematical Modelling 16.11 (1992),
pp. 598–604.
N. Mitarai and H. Nakanishi. “Granular Flow: Dry
and Wet”. In: The European Physical Journal Special Topics 204.1 (2012), pp. 5–17.
J. J. Moreau. “Unilateral Contact and Dry Friction in Finite Freedom Dynamics”. In: Nonsmooth
Mechanics and Applications. Ed. by J. J. Moreau
and P. D. Panagiotopoulos. Springer Vienna, 1988,
pp. 1–82.

15

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]
[26]

[27]

[28]

[29]
[30]

[31]

A. A.-F. Nuseirat and G. Stavroulakis. “A complementarity problem formulation of the frictional
grasping problem”. In: Computational Methods in
Applied Mechanics and Engineering (2000).
P. Painlevé. “Sur les lois du frottement de glissement”. In: C. R. Acad. Sci. Paris 121 (1895),
pp. 112–115.
C. Popa, T. Preclik, and U. Rüde. “Regularized
solution of LCP problems with application to rigid
body dynamics”. In: Numerical Algorithms 69.1
(2015), pp. 145–156.
T. Preclik. “Models and Algorithms for Ultrascale
Simulations of Non-smooth Granular Dynamics”.
PhD
thesis.
Friedrich-Alexander-Universität
Erlangen-Nürnberg, 2014.
T. Preclik and U. Rüde. “Ultrascale simulations
of non-smooth granular dynamics”. In: Computational Particle Mechanics 2.2 (2015), pp. 173–196.
J. Sauer and E. Schömer. “A Constraint-Based
Approach to Rigid Body Dynamics for Virtual Reality Applications”. In: Proceedings of the ACM
Symposium on Virtual Reality Software and Technology. 1998, pp. 153–162.
Y. Shen and W. Stronge. “Painlevé paradox during
oblique impact with friction”. In: European Journal of Mechanics-A/Solids 30.4 (2011), pp. 457–
467.
D. E. Stewart. “Convergence of a Time-Stepping
Scheme for Rigid-Body Dynamics and Resolution
of Painlevé’s Problem”. In: Archive for Rational
Mechanics and Analysis 145.3 (1998), pp. 215–260.
D. E. Stewart. Dynamics with Inequalities: impacts
and hard constraints. SIAM, 2011.
D. Stewart. “Rigid-Body Dynamics with Friction
and Impact”. In: SIAM review 42.1 (2000), pp. 3–
39.
D. Stewart and J. Trinkle. “An Implicit TimeStepping Scheme for Rigid Body Dynamics with
Inelastic Collisions and Coulomb Friction”. In: International Journal of Numerical Methods in Engineering 39.15 (1996), pp. 2673–2691.
W. J. Stronge. “Rigid body collisions with friction”. In: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering
Sciences. Vol. 431. 1881. The Royal Society. 1990,
pp. 169–181.
W. J. Stronge. Impact mechanics. Cambridge university press, 2004.
A. Tasora and M. Anitescu. “A Convex Complementarity Approach for Simulating Large Granular Flows”. In: Journal of Computational and Nonlinear Dynamics 5.3 (2010), pp. 1–10.
A. Tasora and M. Anitescu. “A matrix-free cone
complementarity approach for solving large-scale,
nonsmooth, rigid body dynamics”. In: Computer
Methods in Applied Mechanics and Engineering
200.5 (2011), pp. 439–453.

